Fusionize: Improving Serverless Application Performance through Feedback-Driven Function Fusion

Serverless computing increases developer productivity by removing operational concerns such as managing hardware or software runtimes. Developers, however, still need to partition their application into functions, which can be error-prone and adds complexity: Using a small function size where only the smallest logical unit of an application is inside a function maximizes flexibility and reusability. Yet, having small functions leads to invocation overheads, additional cold starts, and may increase cost due to busy waiting. In this paper we present Fusionize, a framework that removes these concerns from developers by automatically fusing the application code into a multi-function orchestration with varying function size. Developers only need to write the application code following a lightweight programming model and do not need to worry how the application is turned into functions. Our framework automatically fuses different parts of the application into functions and manages their interactions. Leveraging monitoring data, the framework optimizes the distribution of application parts to functions to optimize deployment goals such as end-to-end latency and cost. Using two example applications, we show that Fusionize can automatically and iteratively improve the deployment artifacts of the application.

I. INTRODUCTION
In recent years, serverless computing as the latest incarnation of cloud computing has become more and more popular [1], [2]. In serverless computing, developers rely on managed cloud platform services as application building blocks which are connected through application-specific glue code. Usually, this glue code is written as sequences of small, stateless functions [3], [4] -often invoked in an eventdriven model -running on top of Function-as-a-Service (FaaS) platforms such as AWS Lambda 1 or Google Cloud Functions 2 .
While some work has been devoted to sizing of functions, e.g., [5]- [7], the question of an optimal function size has not been solved yet: From a software design perspective, functions should be small and cover only an isolated piece of functionality to maximize flexibility and reusability. From a performance perspective, however, fine-grained functions lead to increased invocation overheads and cold starts [8], i.e., having only a single large function would be optimal. From a cost perspective, functions should prevent waiting for other functions to finish, since this leads to double billing as shown in Figure 1 [9]. Finally, from a platform perspective, functions need to comply with platform limits regarding resource usage or execution duration which are platform-dependent and may change over time.
In this paper, we propose to close this gap by differentiating between the developer-defined software design artifact (tasks) and the deployment artifact (functions) through a system called FUSIONIZE: During software design, developers follow a fine-granular programming model, composing tasks to build their application. Then, during the deployment process, the FUSIONIZE framework uses function fusion to transform the software design artifact into the actual deployment artifact. Starting from a basic deployment artifact, FUSIONIZE deploys the artifact and uses monitoring data to gradually transform the deployed artifact towards a more efficient version: an artifact with an optimal cost-performance profile, adhering to platform limits. The FUSIONIZE architecture can be used by both cloud users and providers. For this, we make the following contributions: applications at runtime using an optimization strategy (Section III). • We discuss a first heuristic for such an optimization strategy (Section IV). • We implement FUSIONIZE as an open-source prototype and present the results of extensive experimentation with two example applications (Section V). • We discuss the limitations of our approach and derive avenues for future work (Section VI).

II. BACKGROUND
In this section, we give a brief overview of FaaS as a paradigm, discuss influence factors on performance and cost, and finally function fusion as a concept.
FaaS is a popular paradigm for cloud software development and deployment in which developers write their source code in the form of small, stateless functions 3 and leave all operations aspects to the cloud provider. These functions are usually written in high-level programming languages such as Node.js [10]. At runtime, the cloud provider routes incoming requests or events to an instance of the corresponding function -either to a newly created one thus incurring the so-called cold start latency [8], [11] or reusing an existing one. In practice, typical FaaS applications are compositions of short-lived, transient, often highly concurrent function instances [2], [8].
With current FaaS providers, developers have to manually choose the amount of resources (i.e., memory or disk space) they wish to provision per function instance. Overall cost is then determined based on the resulting resource price multiplied with the execution duration of functions. Furthermore, there is usually additional cost per request. Depending on the composition (orchestration vs. choreography [8], [9]), FaaS applications often suffer from double billing when one function makes a blocking call to another one [9]. For optimizing costs, it is hence crucial to avoid double billing and the provisioning of unused resources.
Performance of FaaS functions again depends on the amount of provisioned resources, on cold starts, and the number of requests. Of these, the amount of resources has only a small effect whereas the performance costs for function calls are significant. For instance, in AWS function calls are routed through the API Gateway and in Apache OpenWhisk requests are usually queued in Apache Kafka before execution. On top of this, cold start latency depends on programming language and other parameters but is generally in the range of about one second. Furthermore, increasing load can cause cascading cold starts in a function composition [8], [11]. From a performance perspective, it is hence crucial to have as few invocations as possible which are not only costly but also each carry the risk of encountering a cold start. The mechanism to address this has been referred to as function fusion [6]: Two or more functions are rewritten (and redeployed) as one.
III. AUTOMATICALLY ADAPTING FUNCTIONS AT RUNTIME In this section, we give an overview of FUSIONIZE and describe how it can automatically adapt FaaS functions at runtime. We start by describing key terms and definitions (Section III-A) before discussing FUSIONIZE and its components (Section III-B) and the resulting programming model (Section III-C).

A. Terms and Definitions
As already introduced earlier, we use the term task to refer to the function written by the developer and the term function to refer to the executable deployment artifact. Each function contains a single fusion group, i.e., one or more tasks that are executed as part of that specific function. An application comprises one or more (FaaS) functions (and hence at least as many tasks and fusion groups). While an application is deployed, the assignment from tasks to fusion groups can change dynamically. We refer to a specific set of fusion groups, i.e., one possible deployment configuration of all tasks in an application, as fusion setup. In this paper, we use a short notation to describe fusion setups in the text and in figures: Tasks that are in the same fusion group are put in parentheses, separated with commas; the fusion setup separates fusion groups with hyphens. As order only matters for execution but not for deployment, we order all task names alphabetically. For example, in the fusion setup (A,B)-(C), the tasks A and B are in the same fusion group, and task C is in its own group.

B. The FUSIONIZE Approach
The optimal fusion setup depends not only on developer preferences (regarding cost and performance goals which may be in conflict) but also on runtime effects. A key design goal for FUSIONIZE was hence to implement a feedbackdriven, autonomous deployment framework that collects and uses FaaS monitoring data to iteratively fuse tasks to optimize the fusion setup.
As we show in a high-level overview of our approach in Figure 2, FUSIONIZE has two main components: the Fusion Handler and the Optimizer. In addition, there are also the Log Storage, responsible for storing monitoring data, and the API Endpoint which is the FaaS platform specific service accepting requests to FaaS functions, e.g., API Gateway on AWS.
The Fusion Handler is responsible for request dispatching from and to tasks. As there is one function per fusion group ( A in Figure 2), the Fusion Handler is implemented as a distributed component which is co-deployed with every function similar to the choreography middleware in [8]. From the FaaS provider perspective, the Fusion Handler is the endpoint for incoming calls to that function ( B ). Internally, the Fusion Handler calls the requested tasks either locally ( C ) or remotely ( D ).
We also collect monitoring data from each function. This data includes the execution durations of each of the tasks, memory usage of the function, and which tasks are called from each of the tasks in the fusion group. In our prototype, we rely  on the logging features of the FaaS platform to collect said monitoring data. In case of a FaaS platform not supporting monitoring or logging, the Fusion Handler could easily be extended to also measure the desired metrics and send them to an external, self-implemented monitoring application as the last activity before returning after a call.
The Optimizer retrieves monitoring data from the Log Storage ( E ), automatically derives the call graph of the application, and annotates it with execution information, e.g., latency values. In a second step, it uses an extensible optimization strategy module (we describe a first heuristic for this in Section IV) to derive an improved fusion setup and triggers a redeployment of all functions that have been changed ( F ). We propose to use an adapted version of the continuous sampling plan (CSP-1) [12], [13] to decide when to run the Optimizer: The original algorithm uses the quality of previous items in a factory production process to decide when to run the next quality inspection. In the adapted version, we propose to compare cost and performance metrics of the monitoring snapshots considered during the previous and the current Optimizer run. This way, the Optimizer will be run frequently for a newly deployed application but will still from time to time check performance and cost of "older" applications. Applications can hence adapt in the presence of cost or performance changes resulting from external factors such as load profiles or changes to platform limits.
Within the Optimizer, the "best" fusion setup can be determined in various ways. For instance, we could choose the fusion setup with the lowest median latency, which will be the overall fastest fusion group, or we could choose the fusion setup with the lowest 99 percentiles latency, which will be the fusion setup with the fewest cold starts. Similarly, we could consider various cost metrics. As part of the optimization strategy, application developers should here assign weights to different optimization goals.

C. The FUSIONIZE Programming Model
From a FaaS developer perspective, the programming model is similar to standard FaaS programming. This is the case since the Fusion Handler acts as a wrapper, simply exposing the interface of the respective task(s). The key difference is for outgoing requests: Tasks here do not make direct calls to the target FaaS function but rather request a call to the target task from the Fusion Handler. Overall, this means that existing FaaS applications can easily be migrated into the FUSIONIZE programming model.

IV. A HEURISTIC FOR OPTIMIZING FAAS FUNCTIONS
In this section, we introduce a heuristic which improves the fusion setup of FaaS applications using FUSIONIZE. Please note that the focus of this work is on the systems aspects of FUSIONIZE and not on the optimization strategy, which we leave to future work. Here, we propose a set of rules which will allow FUSIONIZE to create a good (but not necessarily the best) fusion setup.
Fusing tasks can have various effects on performance and cost. Consider the example of a setup in which one task calls another task and has to synchronously wait for the result before continuing. If the called task is in another fusion group this will lead to double billing. Thus, those two tasks might benefit from being in one fusion group. A task calling another task async sync sync

T1 T4
Invocation T2 T3 Response Fig. 3. Example call graph that could benefit from function fusion: Task T1 and T2 need to finish before a response can be sent, and could thus be fused to reduce costs. T3 is only called asynchronously, and should thus be moved to another fusion group.
asynchronously does not lead to double billing. However, the problem of cold starts and call overhead still remains. Another example is shown in Figure 3: T1 receives a request, makes a synchronous call to T2, which makes an asynchronous call to T3 (e.g., for logging purposes) which finally calls T4 synchronously. T1 and T2 would benefit from being in one fusion group to avoid double billing and cascading cold starts [11]. Within one fusion group, however, it does not matter whether a task is called synchronously or asynchronously (aside from concurrency aspects) since the function execution has to wait until all its tasks have completed execution. This means that T2 and T3 would benefit from being in two different fusion groups as only T1-T2 is on the critical path from the application perspective, i.e., also having T3 in that fusion group would unnecessarily increase the call latency for the client or function invoking T1, since the function needs to wait until every local task is finished before it can return.
Furthermore, not only tasks that call each other might benefit from being in the same fusion group. For example, if a task A asynchronously calls either task B or C (and never both) at a fixed rate that can be managed by a single function instance, the amount of cold starts might be reduced if both tasks B and C are handled by the same function instance. Sometimes, it may also be necessary to split tasks into different fusion groups so that the overall function complies with platform limitations such as maximum execution duration or available disk space. Depending on this, our tasks A, B, and C might end up either as one fusion group or as fusion groups (A)-(B,C).
Since there is no difference between calling a local task synchronously or asynchronously, overall request response latency of the original function can be reduced by keeping all long-running asynchronous tasks in separate remote functions. The total billed duration can be decreased by not waiting on remote function calls (which would lead to double billing) and by moving computationally intensive tasks to their own function so that multiple complex tasks do not block themselves. FUSIONIZE can therefore be especially effective if it can move computationally intensive tasks to their own fusion group so that they do not block the critical path. There is also the special case where a task uses up all available resources of a function, e.g., by completely filling up the available disk space. In this case, invocations can only complete if they are placed in their own fusion group.
Bringing all these aspects together, we propose the heuristic shown in Figure 4. Please note that the Optimizer relies on  Fig. 4. Our heuristic calculates the setup that currently performs best and tries to improve it in three ways: If there are asynchronous or no calls between tasks in the same fusion group, move one task into a different fusion group.
If there are synchronous calls between two tasks in different fusion groups, merge those fusion groups. If a new setup has not previously been deployed, it is considered for the next iteration. Alternatively, the best known setup is used. Note that some merges might lead to underprovisioning of resources, which would lead to failed calls. This is not yet handled by our simple heuristic. monitored invocation patterns and not on some kind of composition definition. This means that whenever the invocation patterns change, e.g., because an application has temporal dependencies or a new version has been released, the heuristic will iteratively adapt the fusion setup. Additionally, fusing two existing groups might lead to a function whose requirements exceed those provisioned by the platform. This can only be detected during runtime and will lead to the splitting of these fusion groups in the next run of the optimizer.

V. EVALUATION
To evaluate our approach, we show that FUSIONIZE can be implemented in practice using a proof-of-concept prototype for AWS Lambda (Section V-A). Using our prototype, we explore the effects of function fusion with FUSIONIZE through two example applications. The first application is a sample use case built to demonstrate all features of FUSIONIZE (Section V-B). It consists of multiple tasks that are either quick and called synchronously, or are computationally intensive and called asynchronously. The second case-study is a more realistic implementation of an IoT application, where input from different sensors is used to analyze traffic and stored in a serverless cloud database (Section V-C). We make our artifacts and prototype available as open-source 4 .

A. Prototype Implementation
We implement a prototype of FUSIONIZE for AWS Lambda. We focus on Node.js for this proof-of-concept as it is the most widely-used runtime on AWS Lambda [14] and an interpreted language, which simplifies dynamic loading of code. Please note, however, that our approach is extensible for other programming languages and FaaS platforms.
Internal task calls are implemented with an embedded handler that routes the call internally for fused tasks or externally over HTTPS for tasks that are deployed in a different fusion group. This requires only small changes to the existing programming interface. In our prototype, the Optimizer is implemented as two Lambda functions: The first periodically reads function logs from AWS CloudWatch and writes them to AWS S3. The other uses this information to update the fusion setup.
In our implementation, we chose not to redeploy the entire application whenever the fusion setup is changed. Instead, each function executable contains the source code for all tasks and updating the fusion setup only "rewires" the assignment of tasks to functions by updating an environment variable. This reduces the time required to deploy new fusion setups, but could be infeasible for some applications with large dependencies, e.g., large machine learning models. We have also implemented an extended prototype which generates distinct serverless function deployment packages and compared it to the first implementation (see Section VI-C) Due to our implementation of handling remote requests, the first remote request of a cold function often took ≥1s. All following requests were handled within ≤50ms. This limitation is implementation specific and could be sped up by saving this value in a location that is faster to access, e.g. a managed database service, object storage, or the function image. It could also be possible to use cold start reduction approaches to minimize the impact of these API calls [8].

B. Sample Tree Application
In the sample tree application, some computationally light tasks are called synchronously, while other computationally heavy tasks are called asynchronously. These tasks are called in a tree-like pattern, where one side of the tree consists of synchronous calls, and the other side consists of asynchronous calls. In the call graph shown in Figure 5, we mark all fusion setups that were tested by FUSIONIZE during our experiments.
In the initial fusion setup in all of our experiments, all tasks are in their own fusion group. This is the way that functions would be deployed in a normal serverless system to maximize flexibility and reusability. For all experiments presented here, we configured all functions to use 128MB of memory.
Request Response Latency: The focus of this evaluation lies on demonstrating that FUSIONIZE can be used to optimize the request response latency of invocations. To test this in a normal use case, we put our prototype under a load of one request per second (rps) for 1,000 seconds. This low request rate reflects that of common FaaS workloads [15] and leads to around three active function instances per function. After 1,000 requests, we run the Optimizer with the goal of reducing the median request response latency (rr) until it can find no further fusion setups to test. We then calculate median (rr med ) and average request response latencies (rr avg ) of invocations, sum up the billed duration of all functions that were triggered by an invocation, and plot their empirical cumulative distribution functions.
The results of the full experiment are shown in Figure 6. The Optimizer first changed the initial setup Tree 0 to Tree 1 by fusing tasks A and E together. In Tree 2 , task D is also put in this group. For the final fusion setup Tree 3 , task B is also added to this group leading to the final fusion setup (A,B,D,E)-(C)-(F)-(G). In Tree 3 , the whole synchronous call tree is in the same fusion group. Applying the heuristic reduced rr med from 3.1s to 3.0s and reduced the average billed . When initiating 300 cold starts for different fusion setups, the fastest setup significantly improves application performance: Almost all invocations are faster and cheaper than any invocation of the next-slowest fusion setup, as this setup avoid cascading cold starts in the application by merging tasks with synchronous calls. The fastest setup has an rravg of 3.7s, which is 2.4 times faster than the initial setup (9.2s) and 1.9 times faster than the next-best setup (7.1s). The average billed duration starts at 25.4s and goes down to 8.1s (313% decrease). duration per invocation from 10.5s to 4.1s. This experiment shows that by using function fusion, the median and average latency can be decreased while also reducing cost in a use case where some computationally complex asynchronous functions can be extracted into their own fusion groups.

Cold Start Request Response Latency:
We also run experiments to evaluate how FUSIONIZE can be used to minimize the request response latency of cold starts. In this case, we manually create cold starts by replacing the environment variables of the Lambdas and calculate the median request response latency of cold starts. We create 300 cold starts for every fusion setup, after which the Optimizer runs using the same method we described above. The results of the tree cold start experiments are shown in Figure 7. After testing Tree 0 (rr med =9,206ms, rr avg =8,988ms), the optimizer merges the tasks A and E (Tree 1 ). After also merging D into this fusion group (Tree 2 ), the Optimizer finished at Tree 3 (A,B,D,E)-(C)-(F)-(G) (rr med =3,721ms, rr avg =3,728ms), which is the same setup found in the full experiment. The order the fusion groups were tested in was the same as in the normal experiment, since the call graph is the same. The average billed duration of an invocation decreased from 25.4s to 8.1s.
Discussion: Overall, these results show that FUSIONIZE can successfully decrease the request response latency of normal invocations as well as cold start invocations in an artificial application. The fastest fusion setup we tested was the same for both use cases, which indicates that the heuristic we describe in Section IV can yield a noticeable performance increase over the initial setup.

C. Realistic IoT Application
To evaluate FUSIONIZE in a realistic use case, we present an IoT application, which is a common use case for serverless computing [16], [17]. In this application, sensors placed along roads generate temperature, volume, and air quality readings. These readings are then analyzed by different tasks and stored in an AWS DynamoDB table. An overview of the call graph of this application can be found in Figure 8. As is the case in the tree experiment, all asynchronous tasks calculate primes to emulate complex functions. Additionally, the tasks AS, CSA, DJ, and SE write information into DynamoDB, and the task CSL sends two queries to DynamoDB before writing data itself. The goal of this use case is to test whether FUSIONIZE can optimize more complex fusion setups. For example, one task is called by multiple other tasks (AS is called by CT and CA), and there is a chain of calls that are alternating between synchronous and asynchronous calls (I asynchronously calls CA, which synchronously calls DJ). Otherwise, the experiment setup is the same as in Section V-B. Request Response Latency: In the IoT experiment (cf. Figure 9), the Optimizer tried four fusion setups before arriving at the fastest tested setup IoT 5 ((AS)-(CA,DJ)-(CS,CSA,CSL)-(CT)-(CW,I,SE)). In IoT 1 , CA and DJ are fused together. Next, CS and CSA are fused together (IoT 2 ). In the next step, the optimizer adds CSL to this fusion group (IoT 3 ). Afterwards, CW is merged with SE (IoT 4 ) and finally with I (IoT 5 ). In this setup, all synchronous tasks are located in the same fusion group while asynchronous tasks are kept seperately. The median request response latency of the first four setups drops from 304ms in IoT 1 to 275ms in IoT 4 , while IoT 5 has an rr med of 197ms. Figure 8 shows a visual representation of IoT 1 and IoT 5 Cold Start Request Response Latency: The results for the cold start experiments are shown in Figure 10. The Optimizer tried the same fusion setups it did in the request response latency experiment and finished with the fusion setup IoT 5 , which has the lowest median request response latency of all tested setups (rr med =1,920ms, rr avg =1,922ms). All other setups had a similar distribution of cold start latencies, with a median between 3,423ms and 3,521ms. This can be explained by the order in which our optimization algorithm tried fusion setups: Only the fastest setup merges tasks CW, I, and SE into the same fusion group. In all other cases, task I has to synchronously wait for the completion of CW and SE, which leads to two additional cold starts that happen in the critical path. While the rr of the slower functions is virtually identical, the total billed duration per invocation has different distributions. Overall, invocations for the initial fusion group IoT 0 have the highest cost, since this setup creates the most cold starts. While the fastest fusion setup is also the cheapest one of the tested setups, this shows that request response latency and billed duration do not necessarily correlate.
Discussion: Overall, these results show that FUSIONIZE also works in a more complex application that models a realworld setup. As was the case in Section V-B, the fastest setup was the same for both the normal and the cold start test. This indicates to us that for many workloads, the heuristic of putting all synchronous invocations into the same fusion group and splitting off all asynchronous invocations is a good starting point from which further optimizations might be possible. In all cases tested here, the fusion groups tested were always faster than the previous fusion group. Our simple heuristic might find other local optimums if there are faster fusion setups that are only findable by changing two fusion groups at once. Since the average billed duration of an invocation in the fastest fusion setup is always lower than the initial setup, we argue that FUSIONIZE reduces the total cost of execution. The prototype also incurs other cost, e.g., the Optimizer function and API Gateway costs, which we argue is not relevant for our evaluation since they are a consequence of our implementation and independent of the architecture. Other FaaS providers might provide different methods to directly call their Functions that do not incur additional cost, while the number of requests and their runtime is likely to stay consistent across different FaaS providers.

D. Framework Overhead
FUSIONIZE adds a handler to every function that manages calling the different tasks. This adds an overhead for every   Fig. 9. In our test of FUSIONIZE on IoT application with 1,000 request per fusion setup, the first four improvements lead to a decreased distribution of request response latency. In the fastest setup, rr is reduced by a further 70%. While faster fusion setups also incur less cost, this difference is not as noticeable as in the tree application. function and task call. In an experiment where we called a single empty task once per second for 200 seconds, the handler on average ran for 1.3ms when the function instance was already warm (standard deviation σ=1.24), and on average ran for 36.6ms in cold starts (σ=23.4). While calling a task locally has almost no overhead, calling a task remotely requires additional time to send an HTTP request to another function. When FUSIONIZE first calls another function, the Base URL to call other functions needs to be determined. In our current implementation, this requires two API calls to API Gateway, which take around 1.1s in total. Consecutive HTTP calls take ≤ 50ms.
The Optimizer adds no additional overhead to function calls, since it runs inside its own function and only reads the Cloud-Watch logs written by every function call. Computing the next fusion setup for every 1000 invocations takes around one second, while extracting the invocation data from CloudWatch sometimes takes considerably longer depending on the number of cold starts due to limitations in the CloudWatch API.

VI. DISCUSSION & FUTURE WORK
FUSIONIZE abstracts even more operational concerns away from users than classical FaaS applications, since it also automates the task of sizing the functions. The comparatively simple heuristic we have implemented always found a faster fusion setup than the naive initial state. Nevertheless, some aspects of our architecture and implication warrant further research.

A. Platform Integration
Our prototype is implemented on top of AWS Lambda and runs inside the functions. This creates a performance overhead for every function invocation and limits the information available to the framework. The same architecture could be implemented as a part of the platform, which could lead to increased performance due to less overhead and additional information that can be used to, e.g., allocate tasks to functions or place related function instances in the same physical machine. For example, call graph analysis might be useful to preemptively start functions in anticipation of tasks that will need it. This could also solve platform limitations we encountered during our implementation. When extracting the call graph information from the logs, the amount of requests that can be sent to CloudWatch to query logs is limited. The creation of the call graph for cold-start heavy workloads took multiple minutes, which in some cases lead to timeouts of the handling Lambda function. If this task was managed by the platform, the call graph could be created directly. Additionally, the wait time when sending the first remote request could also be decreased by using provider knowledge. While integration into the FaaS platform could increase performance, our architecture also works when deployed by application developers and can be used until function fusion is supported by platforms. Because the prototype currently only works for Node.js, the Black Box Principle of the Serverless Trilemma [9] is violated, i.e., the application source code needs to be in JavaScript and must be available to FUSIONIZE.

B. Function Resources
In our prototype, all functions have access to the same amount of resources. As, e.g., Akhtar et al. [5] and Eismann et al. [18] have shown, giving different fusion groups access to different amounts of resources might further improve deployment goals, as the amount of available resources directly influences cost per invocation of a serverless function. A further avenue of research is the addition of hardware accelerators to serverless architectures. Some FaaS platforms offer optional support for hardware accelerators, such as GPUs that can be used by functions 5 . While some tasks can be massively sped up by using functions with access to hardware accelerators, they are also more costly to run. FUSIONIZE could be extended to handle the grouping of tasks to functions with or without hardware accelerators to further optimize applications.

C. Experiments
The use cases for our evaluation are not real-world serverless applications. We would have preferred to adapt a nontrivial open source serverless application in our experiments, but were not able to find any suitable candidates. Instead of image recognition, our prototype uses tasks that compute primes. This makes it easier to test different levels of CPU usage and does not rely on other services (e.g. object storage), which could influence latency.In our prototype, all fusion setups use the same deployment package. We also tested whether changing deployment packages to only contain the necessary tasks for their fusion group significantly changes our metrics and found that the average and median latencies as well as costs were within 3% of each other. We expect that the impact of changing the deployment unit might be bigger in cases where some tasks have a bigger deployment size than we used in our evaluation, which is around 10kB.

D. Programming Model
In previous work [19], we have presented an approach using transpiling to fuse tasks in fusion groups, rather than the fusion handler we use in our prototype. While not an issue in our experiments, this would decrease the deployment package to only the code that is necessary for the function group to run.
In both approaches, the code points that are suitable for function fusion, i.e., task entry points, need to be clearly marked by developers. Thus, they are mainly intended to be used when developing new applications and not to transform legacy applications into serverless applications. Spillner et al. [20] have presented a framework that transforms a Python application into (FaaS) functions, and we may integrate such an approach in FUSIONIZE to further abstract operational concerns.

E. Invocation-Dependent Fusion Groups
In our approach, fusion setups are determined only by information about previous invocations, leading to a performance profile of the application. These fusion setups are static in the sense that they only change after the Optimizer runs. Yet it may also be feasible to select a fusion setup based on the type of invocation. For example, the fusion handler could change its behavior if it detects a cold start, or if the input data adheres to certain properties. This dynamic behavior would add additional time to the duration of the fusion handler, so more complex computations to determine an optimal fusion setup might be voided by the additional duration this computation adds to the total function duration.

VII. RELATED WORK
In our previous work [19], we have explored the concept of function fusion. In this vision paper, code is dynamically transpiled for the different fusion groups. We found that loading the code dynamically during runtime offers more flexibility during prototyping and allows the function code to stay the same, so that switching of fusion groups does not require the code to be re-uploaded. Elgamal et al. [6] present an algorithm that minimizes the cost of functions while keeping the latency in a bound by using function fusion and placing the functions at specific edge locations using AWS Greengrass. Their approach uses AWS Step Functions, which they identify as major cost driver in their implementation. In comparison, this paper focuses on AWS Lambda without using edge services, e.g., AWS Greengrass. In this work, we assume that tasks are supposed to be deployed as serverless cloud functions, yet previous work has also considered the optimal placement of tasks over a varied set of compute services such as Container-as-a-Service platforms [7], virtual machines [21], or in a fog environment [22], [23].
While we target applications composed of multiple tasks in this paper, reducing latency and cost of executing a single serverless function, e.g., by reducing cold starts, has been discussed in multiple previous studies (e.g., [8], [24], [25]). Others have used statistical analysis [5] or machine learning [18] to predict optimal function parameters by only looking at some function configurations. This reduces the number of tests to find optimal function sizes. The CPU and memory footprint of functions can be reduced by sharing memory between functions running on the same virtual machine [26], by using application-level sandboxing [27], or by using unikernels [28]. These ideas are complementary to FUSIONIZE and may further improve performance of deployed fusion groups.
Ali et al. [29] present a method to optimize cost and minimize latency by batching multiple invocations into a single function execution. Such an approach especially useful if cold starts account for a significant part of the application duration, e.g., if the function needs to download big datasets during startup. This requires queuing requests, whereas applications optimized by FUSIONIZE can process events immediately.

VIII. CONCLUSION
In this paper we have presented FUSIONIZE, an approach to optimize FaaS deployments by combining tasks in fusion groups. Developers only need to write the application tasks following a familiar function programming model, while our framework automatically fuses different parts of the application into functions and manages their interactions. By abstracting the invocation of subsequent tasks, this framework can handle calls locally or remotely in another function. Leveraging monitoring data, FUSIONIZE optimizes the distribution of application parts to functions to incrementally optimize deployment goals such as end-to-end latency and cost.
Further, we have presented a heuristic for optimizing FaaS functions based on call behavior of the application. Using a proof-of-concept prototype of FUSIONIZE for the Node.js runtime on AWS Lambda, we have shown that our heuristic can improve request response latency and cost in both a sample application and an IoT use case. In future work, we hope to present further optimization algorithms using this prototype and to integrate this functionality at a platform level.


extfieldsofstudy:
Mathematics
Computer Science

text:
A micromorphic phase-field model for brittle and quasi-brittle fracture

The phase-field model for fracture, despite its popularity and ease of implementation comes with its set of computational challenges. They are the non-convex energy functional, variational inequality due to fracture irreversibility, the need for extremely fine meshes to resolve the fracture. In this manuscript, the focus is on the numerical treatment of variational inequality. In this context, the popular history-variable approach suffers from variationally inconsistency and non-quantifiable nature of the error introduced. A better alternative, the penalisation approach, has the potential to render the stiffness matrix ill-conditioned. In order to circumvent both aforementioned issues, a micromorphic approach towards phase-field fracture modelling is proposed in this manuscript. Within this approach, a micromorphic extension of the energy functional is carried out. This transforms the phase-field into a local variable, while introducing a micromorphic variable that regularises the fracture problem. This reduction the regularity requirements for the phase-field enables an easier implementation of the fracture irreversibility constraint through simple 'max' operation, with system level precision. Numerical experiments carried out on benchmark brittle and quasi-brittle problems demonstrate the applicability and efficacy of the proposed model for a wide range of fracture problems.

The efficacy of the novel micromorphic phase-field fracture model is demonstrated with benchmark problems in both brittle and quasi-brittle fracture. These include a single edge notched specimen under tension, the Winkler L-panel experiment [47], and a three-point bending experiment carried out in [48]. This manuscript is structured as follows: Section 2 introduces the reader to the phase-field model for fracture, its underlying energy functional and pertinent Euler-Lagrange equations. Subsequently, in Section 3, the micromorphic approach for phase-field fracture is introduced. The numerical benchmark problems are addressed in Section 4, followed by concluding remarks in Section 5.
2 Phase-field fracture model 2.1 The energy functional Figure 1a and 1b illustrates a discrete fracture and a phase-field regularised fracture in a 2D continuum respectively. The salient feature of the a discrete fracture model is the explicit representation of the fracture through displacement discontinuities. The phase-field fracture model belongs to the smeared representation of fracture, where the phase-field, ϕ ∈ [0, 1] interpolates between intact and fully broken (fractured) material states. The fractured continuum in both figures is assumed to occupy a domain Ω ∈ R dim (dim = 2 in this case). The external boundary Γ comprises of the Dirichlet and Neumann boundaries, represented by Γ u D and Γ u N respectively. Note that Γ = Γ u D ∪ Γ u N and Γ u D ∩ Γ u N = ∅.  . diffused (smeared) crack, with Dirichlet and Neumann boundaries indicated as Γ u D and Γ u N respectively. Figure adopted from [32] .
The energy functional pertaining to the phase-field fracture model [33] is given by, where, a degradation function g(ϕ) is attached to the fracture driving strain energy density Ψ + , with Ψ − and t u p being the residual strain energy density and the prescribed traction respectively. The last integral in the above equation represents the fracture energy, where G c and l are the Griffith fracture energy and the fracture length-scale respectively, while c w is a normalisation constant associated with the choice of the locally dissipated fracture energy function w(ϕ). Some commonly adopted phase-field fracture models are presented in the Table 1. The quasi-brittle phase-field fracture model requires additional parameters p, a 1 , a 2 and a 3 , as observed from the Table 1. These parameters are chosen in way, such that different traction-separation laws may be obtained. For a detailed derivation of these parameters, the reader is referred to [22]. Following [22], the constant a 1 is given by, Quasi-brittle [22] 2ϕ − ϕ 2 π (1 − ϕ) p (1 − ϕ) p + a 1 ϕ + a 1 a 2 ϕ 2 + a 1 a 2 a 3 ϕ 3 Table 1: Brittle and quasi-brittle phase-field fracture models.
where, the newly introduced parameters E 0 and f t represent the Young's Modulus and the material tensile strength respectively. The other parameters vary for different traction-separation laws, and given in Table 2.  Table 2: Quasi-brittle phase-field fracture model parameters for different traction-separation laws [22].
As demonstrated in Equation (1), the strain energy density is split into fracture driving and residual components, Ψ + and Ψ − . In this manuscript, the common consensus that a fracture occurs under tensile loading is adopted. Therefore, the fracture driving and residual strain energy densities are computed using the Amor split [49]. With this approach, the fracture is driven by the deviatoric part of the strain energy density together with a positive volumetric part of the strain energy density. The negative volumetric part corresponds to the residual energy. Mathematically, and Here, K and µ represents the bulk and the shear modulus of the material respectively, tr is a trace operator, • ± represents the positive/negative Macaulay brackets, and dev is the deviatoric strain. The corresponding Cauchy stresses are given by, and adopting the heaviside function H.
In order to simulate fracture initiation and propagation in continuum Ω, minimisation of the energy functional (1) is required, w.r.t. its solution variables u and ϕ. This results in the Euler-Lagrange equations. Along with appropriately defined test and trial spaces, they result in the following problem: using pertinent time-dependent Dirichlet boundary conditions u p on Γ u D and ϕ p on Γ ϕ D , and and Neumann boundary condition t u p on Γ u N . The trial and test spaces are defined as Note that the requirement ϕ ≥ n ϕ in (8c) ensures fracture irreversibility, with n referring to the previous time-step.
From Problem 1, it is clear that the fracture irreversibility constraint manifests in the form of a variationally inequality Euler-Lagrange equation (7b) with restrictive trial and test set (8c) for the phase-field. As mentioned in Section 1, several researchers have proposed different methods to treat the variational inequality problem. In the next section, a micromorphic phase-field fracture model is proposed, where the phase-field variable becomes a local quantity, while a micromorphic variable ensures regularisation of the problem. This enables a simpler 'point-wise' treatment of the fracture irreversibility constraint.

Energy functional and variants
In this sub-section, a micromorphic extension of the phase-field fracture energy functional (1) is carried out in the spirit of [40]. This results in, Here, d is the 'new ' micromorphic variable that regularises the fracture problem. From the above equation, it is observed that the micromorphic variable d is solely associated with the gradient term in the fracture energy integral, and the phase-field ϕ becomes a local quantity. Consequently, the regularity requirements on the phase-field w.r.t. the existence of derivatives is circumvented. Moreover, the micromorphic approach also introduces an additional energy term associated with the difference between the phase-field and the micromorphic variable. Theoretically, for the continuous problem, in the limit, the interaction parameter, α → ∞, the original energy functional (1) is recovered.

Euler-Lagrange equations
The set of Euler-Lagrange equations for the micromorphic phase-field fracture model is obtained upon minimising the energy functional (9) w.r.t. its solution variables u, ϕ and d. Along with appropriately defined test and trial spaces, it results in the following problem: using pertinent time-dependent Dirichlet and Neumann boundary conditions, u p on Γ u D , and t u p on Γ u N respectively. The trial and test spaces are given by, In Problem 2, the phase-field evolution equation (10b) is local. As such, ϕ can be computed 'point-wise' in the computational domain. In particular, the root(s) of the possibly nonlinear scalar equation, yield(s) the phase-field ϕ, if n ϕ < ϕ < 1. Using the locally computed phase-field, the complete problem for micromorphic phase-field fracture is assumes the form: using pertinent time-dependent Dirichlet boundary conditions u p on Γ u D , and Neumann boundary condition t u p on Γ u N . The trial and test spaces are defined as with the local phase-field ϕ computed using (12).
It is worth mentioning that the local phase-field evolution (12) is linear for brittle AT1 and AT2 fracture models. This is due to the quadratic nature of the degradation function g(ϕ) = (1 − ϕ) 2 coupled with linear/quadratic locally dissipated fracture energy function w(ϕ) = ϕ and ϕ 2 for the AT1 and AT2 models respectively. This yield explicit expressions for the local phase-field variable, for AT1, and ϕ = min max 2Ψ + + αd for AT2 model respectively. In the case of quasi-brittle phase-field fracture model, the rational nature of the degradation function (see Table 1) results in nonlinear scalar equation for the local phase-field variable, warranting the need for the Newton-Raphson method.

Discrete equations
The work in this manuscript is executed within the framework of the finite element method [51,52], using triangular (T3) elements for both, displacement and micromorphic fields. This allows assuming the displacement and the micromorphic fields at the nodes (ũ i ,d i ) as the primary unknowns, with the corresponding continuous fields (u,d) approximated as, In the above equation, N u i and N d i are the interpolation functions for the displacement and the micromorphic phase-field, associated with the i th node. The spatial derivatives of the interpolation functions N u i and N d i in a two-dimensional case are given by, Here, the subscripts , x and , y indicate spatial derivatives in x and y directions respectively. Using (18), the strain , and the gradient of the micromorphic variable ∇d are defined as, The discrete phase-field fracture problem is obtained upon inserting (17)(18)(19) in the Euler-Lagrange equations from Problem 3. Thereafter, (13a) and (13b) are assumed as the internal forces, and stiffness matrix derived from its derivative. This notation is consistent with [53], and allows the presentation of the phase-field fracture problem in the incremental iterative framework as: Discrete Problem 1. Compute the solution increment (∆ũ, ∆d) i+1 in the current iteration i + 1 using and update the solution fields, until a certain convergence measure is fulfilled. The local element stiffness matrices are computed as: and the local internal force vectors are computed as The external force vectors f int,u and f int,d are considered equal to zero. Computation of the material stiffness D in Voigt notation is presented in Appendix A.
The discrete problem 1 is non-convex. As such, conventional incremental solution techniques, like the Newton-Raphson method fails to achieve convergence in the softening regime, possibly due to an indefinite stiffness matrix. Efforts to circumvent this issue include the novel line search technique proposed in [34], the use of arc-length solvers [54][55][56], trust-region methods [57] and convexification via extrapolation of the phase-field for the momentum balance equation [36]. These aforementioned techniques are within the framework of monolithic solution techniques. In [39], the alternate minimisation solution techniques was proposed, since the phase-field energy functional is convex with respect to either of the coupled field, if the other one is held constant. In this manuscript, a monolithic solution technique is desired. As such, the convexification strategy proposed by [36] is adopted. With this choice, the phase-field for the momentum balance equationφ, is computed using, withd being the linearly extrapolated micromorphic variable from the last two converged (time) steps. It is important to note that this extrapolation-based computation of the phase-field applies only to the momemtum balance equation and not the micromorphic variable evolution equation. For the micromorphic variable evolution equation, the phase-field is computed using (12), restated below for clarity, Consequently, the discrete problem assumes the form: Discrete Problem 2. Compute the solution increment (∆ũ, ∆d) i+1 in the current iteration i + 1 using and update the solution fields, until a certain convergence measure is fulfilled. The local element stiffness matrices are computed as: and the local internal force vectors are computed as The local phase-fieldsφ and ϕ are computed using (21) and (22) respectively. The external force vectors f int,u and f int,d are considered equal to zero. Computation of the material stiffness D in Voigt notation is presented in Appendix A.

Numerical Study
In this section, numerical experiments are carried out on brittle and quasi-brittle phase-field fracture problems. The problems in brittle fracture include the single edge notched specimen under tension and shear, and the three-point bending test [39], whereas for quasi-brittle fracture, the Winkler concrete L-panel experiment [47] and the concrete three point bending experiment [48] are analysed. For each problem, the geometry, loading conditions as well as the additional model parameters are presented in the respective sub-sections. The phase-field fracture topology in the final step of the analysis, and the load-displacement curves are also shown therein. All problems are solved in a fully coupled (monolithic) sense, adopting the Newton-Raphson method. The iterative procedure is terminated when an error measure defined as ratio of the norm of the residual in the current iteration to that of the first iteration is less than a certain tolerance, tol. In this manuscript, tol = 1e − 3. The linear problem within each iteration is solved using the GMRES solver, preconditioned with an incomplete LU preconditioner. Moreover, for all numerical experiments, the interaction parameter α is parametrised as with β being a user-defined non-dimensional scalar.

Single Edge Notched specimen under Tension (SENT)
The single edge notched specimen [39] has been studied extensively under tensile and shear loading in the phase-field fracture literature. The geometry consists of a unit square (in mm) embedded with a horizontal notch, midway along height and equal to half of the edge length as shown in Figure 2. The notch is modelled explicitly in the finite element mesh. A quasi-static loading is applied at the top boundary in the form of prescribed displacement increment ∆u = 1e − 4[mm] for the first 55 steps, following which it is changed to 1e − 6[mm]. The bottom boundary remains fixed. The additional model parameters are presented in Table 3.

Single Edge Notched specimen under Shear (SENS)
In order to perform a shear test, the single edge notched specimen is loaded horizontally along the top edge as shown in Figure 3. The material properties remain same as presented in Table 3. A quasi-static loading is applied to the top boundary in the form of prescribed displacement increment ∆u = 1e − 4[mm] for the first 85 steps, following which it is changed to 5e − 6[mm]. Furthermore, the bottom boundary remains fixed, and roller support is implemented in left and right edges thereby restricting the vertical displacement. Figure 6a shows the load displacement curves obtained using β = 10, 100, 250. The curves obtained using β = 100, 250 is similar to that obtained in [59]. Moreoever, the phase-field topology at the final step of the analysis (see Figure 6b) is similar to that presented in Figure 5 of [59]. For β = 10, the non-local behaviour of the micromorphic model is not sufficient. This exaplains the loss of integrity of specimen at an earlier stage in the loading process (cf. Figures 6a and 4a).

Winkler L-panel
The concrete L-shaped panel studied by [47,60] is considered in this sub-section. Figure 7 shows the geometry as well as the loading conditions. The longer edges of the panel are 500 [mm] and the smaller edges are 250 [mm]. The loading is applied on the edge marked in blue, 30 [mm] in length, and is in the form of displacement increments of ∆u = 1e − 3 [mm]. The additional model parameters are presented in Table 4. Miehe [39] Ambati [58] 100 Kristensen [20] 10 250        Figure 9. Displacement-based load increments of ∆u = 1e − 3 [mm] is enforced throughout the simulation. The model parameters are presented in Table 5. Figure 10a presents the load-displacement curves using the micromorphic phase-field fracture model, with β = 100, 250. The range of the experimentally obtained load-displacement curves is exhibited by the shaded region. For both values of β, the load-displacement curves exhibit a good agreement with the experimental region. Moreover, the phase-field topology in the final step of the simulation in Figure 10b is similar to that from the literature [48].

Conclusions and Outlook
A novel phase-field fracture model is proposed in this manuscript, based on the micromorhic extension of the phase-field fracture energy functional. In this model, the phase-field variable is local, and a micromorphic variable is introduced for regularization. Combined with an irreversibility criterion on the phase-field variable, this results in a local variational inequality problem. The problem is then trivially solved pointwise (at integration points) throughout the computational domain. For brittle AT1 and AT2 models, an explicit expression for the local phase-field exists, however, for quasi-brittle fracture models, a nonlinear scalar equation needs to be solved. Furthermore, the local nature of the phase-field variable also provides the ease in implementing bounds, i.e., ϕ ∈ [0, 1] using the trivial 'min' and 'max ' operations. Unlike the penalization technique proposed in [35], fracture irreversibility is achieved with system precision without the need for a high penalty parameter. As regards to computational effort, since the phase-field variable is solved locally, the size of the finite element problem is identical to that of a standard phase-field implementation.  that ϕ ≈ d. This has been demonstrated in this manuscript with numerical experiments on the single edge notched specimen loaded in tension. When the interaction parameter is sufficiently high (≈ 250G c /l), the phase-field and the micromorphic variable profiles are similar for arbitrary sections in a computational domain. This essentially recovers the original phase-field fracture energy functional, while enabling the ease of local phase-field fracture irreversibility. Numerical experiments were also carried out in the context of quasi-brittle fracture, the Winkler L-shaped panel and the concrete three-point bending tests. It is observed that the micromorphic phase-field fracture model results in fracture topologies and load-displacement curves within the experimentally observed range. This demonstrates capability of the micromorphic approach to demonstrate both brittle and quasi-brittle phenomenon. Finally, the novel micromorphic phase-field fracture model opens a plethora of future research extensions, particularly, in multi-physics applications, composite laminates [61]. Other studies may include the implementation of a dissipation-based arc-length method [55,56] or quasi-Newton methods [62,63] for addressing the non-convexity of energy functional.

Software Implementation and Data Availability
The numerical study in Section 5 is carried out using an in-house finite element code, inspired by ofeFrac [64], and based on the Jem and Jive C++ libraries from the Dynaflow Research Group, the Netherlands. The simulation data would be made available in Github repository of the corresponding author (https: //github.com/ritukeshbharali).

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


s2fieldsofstudy:
Computer Science
Mathematics

text:
Mathematical Model of Strong Physically Unclonable Functions Based on Hybrid Boolean Networks

We introduce a mathematical framework for simulating Hybrid Boolean Network (HBN) Physically Unclonable Functions (PUFs, HBN-PUFs). We verify that the model is able to reproduce the experimentally observed PUF statistics for uniqueness $\mu_{inter}$ and reliability $\mu_{intra}$ obtained from experiments of HBN-PUFs on Cyclone V FPGAs. Our results suggest that the HBN-PUF is a true `strong' PUF in the sense that its security properties depend exponentially on both the manufacturing variation and the challenge-response space. Our Python simulation methods are open-source and available at https://github.com/Noeloikeau/networkm.

I. INTRODUCTION
Boolean networks (BNs) [1], [2] are connected systems of discrete nodes that exhibit logic-like behavior. They have been used to model gene regulation, design new digital circuits, and have even been studied as a toy model of the human brain. Here, we develop mathematical tools to study a class of BNs that are both challenging to measure experimentally, and useful as cryptography primitives known as Physically Unclonable Functions (PUFs). Specifically we model ultra-fast chaotic circuits known as Hybrid Boolean Networks (HBNs) implemented on Field-Programmable Gate Arrays (FPGAs).
In this context, the nodes of the network are FPGA logic elements, and the edges connecting them are routing resources that conduct the electrical signals representing the analog value of each node state. Unlike standard synchronous designs, nodes exchange signals and perform logical operations continuously -without the control of a global clock or input from the user. Each BN is therefore an autonomous, recurrent electrical circuit that self-oscillates, typically at the maximum frequency allowed by the hardware, which is in the GHz regime for a Cyclone V FPGA. In practice, a clock is used to set the initial state of the system, forming a BN that is a 'hybrid' of autonomous and synchronous components (HBN). This work was supported by the Department of the Army through award number W31P4Q-20-C-0003. DJG and AP have a financial interest in Verilock.id, which is commercializing the PUF technology described here.
These ultra-fast self-oscillations lead to a significant source of entropy generation in HBNs with nonlinear node functions, such as the exclusive-or (XOR) operation. The XOR function is maximally sensitive to individual bit flips, and hence a network of XOR gates can amplify the diverging effects of these bit flips exponentially over time. Thus, chaos has been observed in XOR-HBNs with as little as a single node having two self-loops [3], and more robustly in a system of three nodes [4]. Chaos is consistently observed above 8 nodes, and 16-node XOR-HBNs have been used as true random number generators with a 12.8 Gbit/s entropy rate [5].
Recently, XOR-HBNs have been used to develop PUFs on FPGAs [6]. PUFs act as 'digital fingerprints' by transforming stable environmental information into digital codes used for hardware identification [7]. The atomic-scale manufacturing variations present in the physical properties of FPGA logical elements and wires (such as impedance and length) are amplified and transformed into bits by the dynamics of the HBN. Given an N -bit initial condition or challenge specifying the binary state of each labeled node in the HBN class, the resulting final state or response can be used to uniquely identify the individual FPGA or instance that created it. A large number of challenge-response pairs (CRPs) is guaranteed, as there are 2 N possible N -bit responses for each of the 2 N challenges. As a result, the entropy of HBN-PUFs scales super-exponentially with network size, a characteristic of a 'strong' PUF.
Due to their super-exponential entropy scaling and subnanosecond dynamics, HBN-PUFs stress existing methods of measurement and analysis. Additionally, the FPGA prevents analog readout of the waveform for each node because their primary use is as a digital device. This characteristic significantly limits the amount of information available from experiments. Hence, a model bypassing these limitations is useful for analyzing and certifying the HBN-PUF as an entropy source.
Current HBN models only consider a subset of the realistic experimental conditions present on FPGAs -namely noise on the signal, time-delays on the wires, and differences in logic operation time [3], [4], [8]. Each of these prior works on HBNs have emphasized the importance of variations to these parameters on the resulting dynamics. Here, we develop a mathematical model connecting these parameter variations to the PUF statistics. Specifically, we extend the Glass model [9] to include time-delayed signals [10], noise, and inhomogeneous time-constants. Our unique contributions follow. • We introduce a mathematical model and Python library for simulating delay-differential network dynamics, specializing to HBN-PUFs on Cyclone V FPGAs. Drawing parameters randomly with means taken from the literature, we find that the model reproduces average experimental HBN-PUF uniqueness and reliability statistics. • The model predicts that PUF uniqueness depends exponentially on the inter-die variation of the FPGA parameters, suggesting that the HBN-PUF is a true 'strong' PUF. In Sec. II we introduce the HBN model and define the PUF statistics, then describe our simulation procedure. In Sec. III we simulate HBN-PUFs and compare our results to experimental data. We then use the model to estimate the functional form of the PUF statistics in terms of the amount of manufacturing variation and noise. Finally in Sec. IV, we interpret our results in the context of future work and conclude.

II. MODELING FRAMEWORK
We model HBN dynamics as a set of coupled first-order time-delay differential equations with noise, governing the node states x n ∈ [0, 1] for each of the N nodes in the network. This is summarized in Fig. 1. Each node n executes a logic function f n on the states of its neighbors y n . These neighbor states are time-delayed by ∆ nm , representing the finite propagation time required for the state of node m to reach node n. These time-delayed signals are also rounded toward 0 or 1, as in hardware logic gates, which convert the incoming analog voltages into binary values. Mathematically, where X m (t) = Θ(x m (t) − 0.5) is the Boolean state and Θ is the Heaviside step function, and pred(n) are the predecessors or fixed set of nodes with an edge directed toward node n.
The characteristic timescale τ n determines the time required for the node to transition from logic low to logic high or vice-versa, also known as its rise and fall times, respectively. Thermal/charge fluctuations are modeled as Gaussian noise n (t) ∼ N (0, ε 2 ) drawn randomly with zero mean and standard deviation ε at each time-step. We integrate Eq. (1) using a first-order Euler method mapping the state from time t to t + dt, given by where all times and delays are integer multiples of dt. The optimal PUF read-out time resides in the transient and is within an order-of-magnitude of the node timescale. Therefore, we integrate for only ∼ 10 1τ , and error accumulation due to the first-order approximation is negligible over this interval. Here, τ = (1/N ) n τ n is the mean node time constant. We express all future parameters in units ofτ = 0.25 ns, a value typical of Cyclone V FPGA logic elements [11].
In an actual FPGA, multiplexers are used to set the initial condition of each node in the network as shown in the top right panel of Fig. 1. These multiplexers are held fixed to the challenge for several thousand ns. During this time, the autonomous portion of the network (the XOR gates) stabilizes and synchronously evaluates the challenge x(0). Thus, the initial output of each multiplexer n is f n (y n (0)) (the node function evaluated over its input challenge bits). Experimentally, this is the first value observed in memory.
We observe that simulations explicitly modeling these multiplexers reproduce this behavior. Crucially, we observe that these more complex models produce nearly identical PUF statistics to the simpler models which omit multiplexers but account for this change in initial conditions. Therefore, we omit modeling multiplexers explicitly, and instead use f n (y n (0)) as the initial condition for each node in the simulation.
HBN-PUF challenge-response pairs (CRPs) are multidimensional arrays having indices summarized in Table I. Classes s are characterized by a particular network topology and set of parameters, representing the circuit netlist instantiated on FPGAs. Instances i of a class are characterized by a uniformly random perturbation of that class's parameters, representing the manufacturing variations present in each FPGA's physical structure. Challenges c are the binary vectors specifying the initial condition of an HBN-PUF instance. Each integration r of an initial condition is characterized by a different sequence of noise vectors randomly generated at each timestep. Each node n then responds with a Boolean state at time t. We label this data as X s i,c,r,n (t) ∈ {0, 1} Ns×Ni×Nc×Nr×N ×T . The PUF uniqueness and reliability statistics µ inter and µ intra are averaged pairwise differences between responses. Uniqueness measures the average difference between responses of separate instances to the same challenge. It arises due to manufacturing variations, and ideally µ inter = 0.5. Reliability measures the average difference between responses of the same instance to the same challenge. It arises due to noise, and µ intra = 0 in the ideal noise-free case (ε = 0). Let angle brackets denote averages X j j = (1/N j ) X j , and pairs of indices denote average Hamming distances For each class there exists an optimal time t s opt to select the response from the network timeseries. This is when the network has coupled enough to the chip manufacturing variation to be unique, while remaining unaffected enough by noise to be reliable. Mathematically t s opt is the point maximizing the difference between the statistics t s opt := arg max t≤T ∆µ s (t) .
Typically, the class index s is suppressed, as one is dealing with a single PUF class. However, here we will perform an average over classes, representing the typical statistics of a random network ensemble. In these cases an omission of s indicates an average over classes, i.e., µ = (1/N s ) µ s .
Simulated HBN-PUF statistics are procedurally generated in the following order: 1. class, 2. instance, 3. CRP. In the first step, class generation, an N -node network is drawn from the set of all random regular graphs of degree 3. The delay along each edge is drawn from the Uniform Random distribution over the interval ∆ s nm ∼ U(0, 10τ ). These edge delays represent the location assignment given by the FPGA compiler, which we find can vary from near-zero to several ns depending on the design and routing congestion.
The second step, instance generation, slightly adjusts the parameters of a given class. The node time-constants are τ s i,n ∼ N (τ , σ 2 ), where N is the Normal distribution with meanτ = 0.25 ns and standard deviation σ. The positivedefinite time-delays of instance i are given by ∆ s i,nm ∼ |N (∆ s nm , σ 2 )|. We interpret σ as a uniform manufacturing variation on all FPGA components, and set σ = 0.05 unless otherwise noted, a value typical of modern FPGAs [12]. The final step, CRP generation, integrates Eqs. 1-3 to predict the dynamics of the network and formats this data to match the experiment. The integration time T int = 42τ = 10.5 ns and timestep dt =τ /25 = 0.01 ns are fixed for all simulations. The noise amplitude is fixed to the timestep error ε = 0.01 unless otherwise noted. For each class, instance, and challenge, (3) is integrated N r times using the corresponding parameters. Each integration differs by the unique noise vector drawn as N (0, ε 2 ). The resulting timeseries is then decimated to match the measurement interval of the experimental data, which is 2τ = 0.5 ns. The first 2τ are discarded as they are not observable in experiment, resulting in 20 evenly spaced times over 10 ns from which Eqs. 4-7 are calculated. Figure 2 displays a comparison between actual HBN-PUF experiments on Cyclone V FPGAs (top row) and our model statistics (bottom row). We stress that no attempt is made to fit the parameters of the model to match the experimental data. Nevertheless, the null hypothesis -that the model does not predict the experiment -should be rejected, due to an RMS Z-score

III. RESULTS
√ std(sim) 2 +std(exp) 2 ∼ 1 ∀t, and std(.) are the standard deviations of the µ s (t). In other words, to within 72% probability (integrated PDF), the ensemble averages of the model and the experiment can be drawn from the same normal distribution. Additionally, the limiting behavior and curvature of both statistics match, demonstrating that the model can be used to make asymptotic predictions about HBN-PUFs. However, the individual simulated classes and instances display less diversity than their experimental counterparts. This is likely due to differences in the measurement times of nodes arising from skew on the FPGA clock fabric and registers, jitter on the external crystal oscillator, and correlated delays in the routing, none of which we account for in our random parameter draws. All of these factors contribute to the kinks in the experimental curves largely absent in the model. An understanding of these instance-level temporal correlations requires ps-scale intra-FPGA measurement tools such as a waveform capture device [11]. Measuring individual ∆ nm , τ n and calibrating the model is the subject of future work.
Next, we determine the extent to which the HBN-PUF model scales with manufacturing variation and noise, and how this information quantifies the PUF strength. To this end, we calculate the statistics µ inter and µ intra as a function of σ and ε by fixing the observed optimal response time t opt = 6 ns. These data are shown in Figure 3. In both cases, the data are well fit to a saturating exponential function. This provides mathematical evidence suggesting that the HBN-PUF is truly 'strong'. This is because, in addition to having an entropy that scales super-exponentially with N , the HBN-PUF uniqueness statistic µ inter scales exponentially with the inter-die manufacturing standard deviation σ. Furthermore, the exponential rate constant C fit to µ inter (σ) is roughly three times that of µ intra (ε). This demonstrates that the nonlinear dynamics amplify manufacturing variations to a greater degree than noise, ensuring effective HBN-PUF behavior. Additionally, exponential scaling of µ inter (σ) indicates that small uncertainties in the system parameters result in large deviations of the modeled response. Log-σ simulations (not shown) demonstrate that this uncertainty conservatively approaches the observed noise floor at σ min = 10 −3 , i.e. µ inter (σ min , = 0.01) ≈ µ intra (σ = 0.05, = 0.01). We therefore estimate that the dynamics of each node and edge in the HBN-PUF must be resolved to σ minτ ≈ 10 −13 s (sub-psscale) precision in order for a modeling attack to be feasible. FPGA clock fabrics routinely require ns-scale skew for large recurrent networks like the HBN. However, small, carefully placed networks may be measurable to ps-scale and below over short time periods, and this is the subject of future research.

IV. CONCLUSIONS AND FUTURE WORK
In summary, we accurately predict the class-averaged experimental statistics of HBN-PUFs using stochastic first-order time-delay differential equations and parameter values drawn randomly from the literature. These simulations predict that HBN-PUFs are exponentially sensitive to manufacturing variations to a greater degree than noise. In conjunction with the exponential entropy scaling, our results suggest that the HBN-PUF is a true 'strong' PUF. The introduced simulation library supports arbitrary network topologies and logical functions, and in the future we will combine it with 5-ps-resolution experimental measurements using a waveform capture device on the Cyclone V [11]. These measurements will be used to test the predicted sub-ps uncertainty bound required for a modeling attack. We will also present novel theoretical metrics analyzing strong PUF properties, specifically unclonability, entropy, and multidimensional chaos.