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.