Hamiltonian Extrema of an Arbitrary Flux-Biased Josephson Circuit

Flux-biased loops including one or more Josephson junctions are ubiquitous elements in quantum information experiments based on superconducting circuits. These quantum circuits can be tuned to implement a variety of Hamiltonians, with applications ranging from decoherence-protected qubits to quantum limited converters and amplifiers. The extrema of the Hamiltonian of these circuits are of special interest because they govern their low-energy dynamics. However, the theory of superconducting quantum circuits so far lacks a systematic method to find these extrema and compute the series expansion of the Hamiltonian in their vicinity for an arbitrary nonlinear superconducting circuit. We present such a method, which can aid the synthesis of new functionalities in quantum devices.


I. INTRODUCTION
State-of-the-art superconducting quantum processors rely on nonlinear circuits operating at GHz frequencies. These devices typically consist of a combination of linear microwave circuitry and Josephson tunnel junctions (JJs) [1]. Specifically, a JJ behaves as a nonlinear inductance described by a potential energy function U (φ) = −E J cos φ, where E J is the Josephson energy and φ is the phase drop across the JJ, itself linearly related to the integral of the voltage across the element [2]. The JJ is also characterized by a capacitance C J in parallel with the Josephson element, modeling the influence of the tunnel barrier. This defines a charging energy E C = e 2 /(2C J ), where e is the electron charge. Such a model is rendered in Fig. 1(a), where a JJ is represented as a cross-in-box symbol. The cross symbol is associated with the nonlinear inductive component of a JJ, while the box represents its capacitive shunt. One or multiple JJs, together with superconducting wires, can be combined to form a superconducting loop which can be threaded by an external DC magnetic fluxΦ e . To ease the notation, in the rest of the manuscript we will use the normalized external DC magnetic flux, defined as where Φ 0 = h/2e is the magnetic flux quantum. This flux acts as a tuning parameter which controls the linear and nonlinear properties of the loop. By properly assembling such loops and choosing a particular value ofφ e , it is possible to implement a large variety of quantum systems, ranging from noise-protected qubits [3][4][5][6][7][8] to readout circuits as quantum-limited amplifiers [9][10][11][12][13][14][15] and control circuits such as parametric couplers [16,17]. Superconducting loops are usually modeled as N -body quantum mechanical systems, where N is the number of JJs in the loop. The k-th JJ is characterized by an inductive energy * sandro.miano@yale.edu † michel.devoret@yale.edu E J k and a charging energy E C k . In addition, each loop hosts a linear inductor representing its total geometrical inductance. As a consequence, any superconducting loop could host up to N resonant modes. Typically, a superconducting loop is coupled to external circuitry via two of its terminals [9,14]. For instance, a loop can be shunted by an external capacitor of charging energy E C , as in Fig. 1(b). This defines a reference (ground) node and a main active node (top node) for the loop, as well as two branches. Each branch corresponds to a path from the main active node to the reference node. In the rest of the manuscript, we will refer to the two-terminal device obtained by selecting two nodes of the loop as a dipole.
Note that a loop hosting N elements can generate up to N − 1 different dipoles. A general expression for the classical Hamiltonian of the capacitively shunted loop in Fig.  1(b) is given by a combination of the capacitive (kinetic) energy T cap and inductive (potential) energy U ind H(n, φ,φ e ) = T cap (n) + U ind (φ,φ e ) (2) where φ = [φ 0 , φ 1 , . . . , φ N ] and n = [n 0 , n 1 , . . . , n N ] are the vectors representing the phase and conjugate charge variables associated with the active nodes of the loop, respectively [1]. Note that the φ k are linearly related to the integral of the voltage between node k and the ground node. The phase φ 0 corresponds to the main active node of the dipole, while the others correspond to the internal nodes of the dipole. The Hamiltonian (2) can then be quantized by promoting the variables φ k and n k to the operatorsφ k andn k , respectively, satisfying the commutation relation [φ k ,n k ] = i. The resulting Hamiltonian operatorĤ is related to the classical Hamiltonian (2) via the relationĤ = H(n,φ,φ e ) whereφ = [φ 0 ,φ 1 , . . . ,φ N ] andn = [n 0 ,n 1 , . . . ,n N ] are the vectors of phase and charge operators, respectively. For a given external flux-biasφ e , the eigenstates ofĤ can be computed and further analyzed to yield the properties of the system. This comprehensive approach can be extended to an arbitrary superconducting circuit, and is currently at the core of many general-purpose quantum circuits simulators [18][19][20]. However, this description does not directly capture the number and positions of the local extrema of the classical potential energy U ind (φ,φ e ) in Eq. (2), the knowledge of which is very convenient when designing a device for a particular purpose. For instance, many topological protected qubits [4][5][6] are implemented with Hamiltonians whose ground state wavefunction, in the phase representation, can be imagined as spread across multiple minima of U ind . On the other hand, quantum-limited amplifiers and couplers typically require the presence of a single minimum for U ind to reliably implement the multi-photon parametric processes at their core. These devices are typically described by a series-expansion of their effective Hamiltonian around the minimum. A brute-force approach to determining the number and positions of the potential energy local extrema is to compute, for a fixedφ e , the equilibrium points of the circuit [21], defined as the phase vector φ = [φ 0 ,φ 1 , . . . ,φ N ] which satisfies the set of equations where ∇ φ is the gradient with respect to the components of φ. In the general case, this method requires solving a set of transcendental trigonometric equations with iterative root-finding algorithms. A search for all the extrema of a potential energy function is thus very challenging and time consuming.
Here, we propose a systematic method to describe the equilibrium properties of an arbitrary flux-biased superconducting dipole. The essence of our method is based on the relation where i k (φ) is the total current flowing into the inductive subnetwork through the k−th active node, as shown in Fig. 1(c). This last relation can be demonstrated by computing the instantaneous power entering the inductive subnetwork as where V k is the voltage between node k and the ground node. Being the total differential of U ind can be computed from Eq. (6) as from which expression (5) arises. Consequently, Eq. (4) imposes that, at each equilibrium point, i k (φ) = 0, thus only a persistent loop current [22]ī ℓ is allowed to circulate in the dipole, as represented in Fig. 1 consequence, the x-th element in the loop will have an equilibrium phase dropφ x , such that i In particular, the CPR of a JJ reads while that of a linear inductor reads where I C is the critical current of the JJ and L is the inductance of the linear inductor. We want to emphasize that all the electrical variables at equilibrium are time-independent: as a consequence, their properties are independent from the capacitive subnetwork of the dipole in Fig. 1(b).
In section II, we show how the equilibrium points of a single loop are univocally related to the CPR of a single, open branch constructed from a sequence of all the elements in the loop. Then, we show how a generic CPR can be expressed as a parametric relation. To obtain such a representation, the curvilinear parameter of choice is the phase drop across the JJ with the smallest critical current, which we name the free JJ for reasons that will become clear later. In section III, we study the multistability of a superconducting loop and derive a set of rules that quantify the properties of circuits with multiple equilibrium points. In section IV, we discuss the application of our method to describe the low energy properties of the quantum Hamiltonian (3) which, when E C ≪ E C k , can be assumed independent from the JJs' capacitances. Under such approximation, we introduce an effective potential energy function associated to U ind , whose series expansion coefficients can be computed by systematically combining those of the inductive elements forming the dipole. Finally, in section V we discuss the applications of our approach: optimization and synthesis of superconducting circuits, modeling the effect of fabrication uncertainties and tight-binding approximation for multi-body superconducting circuits.
We have implemented the results of this work in a Python package available on a GitHub repository [23], named 'Nonlinear Inductive Network Analyzer' (NINA). The repository includes examples for common superconducting dipoles, as well as a comparison between our method and the method used in previous literature to analyze a SNAIL [14].
Remarkably, the set of equations (11) does also describe the case where the two branches A and B are arranged in series, as depicted in Fig. 2(b). We call such configuration the equivalent open branch of the loop. From the point of view of the persistent currentī ℓ , the two branches A and B can be regarded to be in series. This observation has some important consequences. For instance, as represented in Fig. 2(c), a set of three different inductive elements can be arranged to form three different flux-biased dipoles. These dipoles are obtained by selecting different pairs of terminals from the same loop. Consequently, the equilibrium phases of each element will be the same for a common value ofφ e , regardless of the terminals arrangement. The three superconducting dipoles would however differ by the shape of their potential energy function, as explained in section IV.
To analyze the equivalent open branch in Fig. 2(d), we introduce the maximum DC current allowed through such branch, I F , corresponding to the critical current of the free JJ. While the phase of the free JJ is unbounded, as rendered in Fig. 3(a), the domain of variation of the phase across the other elements is constrained. The absolute value of the maximum phase drop across the linear inductor and the larger JJ will be given, respectively, by max{|φ L |} = β L and max{|φ J |} = arcsin β J , where Note that β J ≤ 1 by definition. In these last definitions, L is the inductance of the linear inductor, and I J the critical current of the big JJ. Such maximum values correspond, as in Fig. 3(b) and (c), to a flow of current through the elements equal to I F , i.e. to a phase drop across the small JJ,φ F = π/2. For a generic value of φ F , the phases of the constrained elements read where z ∈ Z is the integer multiple of π around which the CPR of the bigger JJ is evaluated. In the rest of the manuscript, we will only consider the case z = 0, as in Fig. 3(c). The properties of a superconducting array for z ̸ = 0 are outside the scope of this work, since they are associated to higher potential energy curves. These curves can still have equilibrium points with positive curvature, but of higher energy with respect to those of the potential energy associated with z = 0. As a consequence, the equilibrium points for z ̸ = 0 can be classified as metastable.
We can now establish the CPR of the branch, observing that the total phase dropφ is obtained by adding the phase drops across the elements, while the currentī can be expressed as the one through the free JJ. We obtain the following system of parametric equations which represents the implicit form of the CPR of the equivalent open branch in Fig. 2 with the curvilinear parameter beingφ F . Note that the total phase dropφ across an equivalent open branch corresponds to the external flux of the associated loop. In Fig. 4, theφ(φ F ), i(φ F ) andī(φ) are computed for β L = 1 and β J = 0.25. The description discussed so far can be extended to the case where two JJs are identical, i.e. β J = 1. In such case, only one of the two JJs should be elected to be the free one, while the other has to be treated as constrained.
To generalize the system (14) to the case of N larger JJs, it is necessary to include all their contributions to the total phaseφ. The simple implicit CPR (14) then becomes the general where β J k = I F /I J k . Since the phases across the constrained elements are analytical functions ofφ F as in (13), the set of equations (15) describes the equilibrium properties of a generic superconducting loop for an arbitrary value of external flux when imposingφ =φ e andī =ī ℓ . The only caveat is that the inverse functionφ F (φ), in general, does not admit an analytical expression. However, it can be approximated via numerical interpolation techniques, up to a desired numerical precision.
This method can also be extended to include nonlinear inductive elements with non-sinusoidal CPRs as nanowires [24] or JJs with higher Josephson harmonics [25].

III. MULTI-MINIMA SUPERCONDUCTING LOOPS
Some loop designs, for instance the fluxonium [3] and cos 2φ [4] qubits, are purposely designed to work in a multi-minima configuration, while others require the presence of a single operating point, as in the case of parametric couplers [16,17] and amplifiers [9,14,26]. To formulate the conditions under which an arbitrary flux-biased superconducting loop has one or multiple operating points, it is useful to acknowledge that its operating points can be determined by analyzing the CPR of its equivalent open branch. In particular, we note how, as displayed in Fig. 4,ī(φ) is a single-valued function if the component of ⃗ γ ′ (π) along theφ F axis is positive, or a multi-valued function if the same component of ⃗ γ ′ (π) is negative. We want to clarify that, even if Fig. 4 describes a particular associated branch, such observation applies to the general case: the generic CPR is a continuous curve, andφ(φ F = π) = π regardless of the branch complexity, as can be obtained from the first equation in (15). The component of ⃗ γ ′ (φ F ) along theφ F axis is given by where we have defined Consequently, the conditions for single-or multi-valued branches are These inequalities generalize those already known for common loops as rf-SQUIDs [27] and SNAILs [14] to an arbitrary superconducting loop. In the multi-minima scenario β tot > 1, for a given value ofφ e there will be a set of equilibrium free phases, each one corresponding to an equilibrum point of the loop. Formally, such solutions are those of the system With reference to Fig. 5, the solutions of this system are the yellow dots at intersection between theφ(φ F ) curve and the yellow lineφ =φ e . To count the number of operating points in an arbitrary loop for any value of external flux, it is useful to note how the function φ(φ F ) has the translational symmetryφ(φ F +2π) = 2π+ φ(φ F ), as can be verified from the first equation in (15). Consequently, the solutions of the system (20) can be mapped to those of an equivalent system which has the benefit of limiting the range ofφ F where to search for solutions. The transition between system (20) and (21) is also represented in Fig. 5, where all the solutions of (20) which fall outside the interval can be mapped to equivalent solutions in the interval φ F ∈ [−π, π) (enclosed in a solid oval). Moreover,φ(φ F ) being a bounded function forφ F ∈ [−π, π), only a finite number of n need to be considered to find all the solutions of the equivalent system (21). With reference to Fig. 5, the upper-and lower-bound on n for a given value ofφ e are respectively given by whereφ max is the maximum of the functionφ(φ F ) for φ F ∈ [−π, π) and ⌊·⌋ represents the floor function. As for each n ∈ {[N ↓ , N ↑ ], n ̸ = 0} there will be two solutions, the total number of operating points, including the one for n = 0, is given by Note that N tot from the last expression is the total number of operating points, both stable (local minima) and unstable (local maxima). The number of stable operating points is given byN In the next section, we will show how to apply the so far discussed technique to describe the low-energy properties of the Hamiltonian (2).

IV. SERIES EXPANSION OF THE EFFECTIVE HAMILTONIAN OF A CAPACITIVELY SHUNTED SUPERCONDUCTING DIPOLE
Common experimental devices based on capacitively shunted superconducting dipoles [3,4,7,14,16] are usually described by phenomenological, effective Hamiltonians which neglect the presence of the loop internal resonant modes. Under such approximation, the average values and the quantum fluctuations of the phases φ k describing the internal nodes of the dipole are considered purely dependent on the values and the quantum fluctuations of the phase φ 0 describing the main node of the dipole. This assumption is reasonable when the charging energy E C k of the k-th JJ in the loop is much smaller than the charging energy E C of the external capacitance shunting the dipole, and it becomes exact in the limit E C k → 0 [28]. Note that, when an arm of the dipole includes a JJ in series with a large inductance, neglecting the internal JJ capacitances can result in a branched Hamiltonian which wouldn't correctly describe the dynamics of the circuit. In such case, a Born-Oppenheimer approximation can be applied to (2), accurately accounting for the presence of vanishingly small JJs' capacitances [28]. In scenarios where the effect of the JJs' capacitances can be neglected, the Hamiltonian (2) can be replaced by the one-body effective Hamiltonian where E C is the charging energy of the shunt capacitance shunting the dipole and U eff is a phenomenological potential energy function which models the entire dipole as an effective, flux-tunable nonlinear inductor. A general expression for U eff can be obtained from the potential energy U ind in (2) by imposing that the phase of the k-th internal node of the dipole φ k is related to the phase of the main node of the dipole φ 0 via a function f k such that φ k = f k (φ 0 ), resulting in The function f k enforces the constraints imposed by Kirchhoff's laws to φ k in absence of the JJs' capacitances, and has an analytical form only for a narrow set of dipoles. If the functional form of U eff is known, the effective Hamiltonian (25) can be expanded, for a fixed value of external flux, around a flux-dependent equilibrium pointφ 0 . The charge and phase variables n 0 and φ 0 can then be promoted to the operatorsn 0 andφ 0 , respectively, yielding the series expansion of the effective Hamiltonian operator where are the Taylor expansion coefficients of U eff around the equilibrium pointφ 0 . Furthermore, it is possible to introduce the ladder operatorsâ andâ † satisfying the relationŝ where are the quantum ground state uncertainties in phase and charge number, respectively. In the basis of these ladder operators, the Hamiltonian (27) readŝ where is the natural frequency of the associated harmonic oscillator and are the rates of n-photon interactions. Expression (30) has been successfully applied to quantum circuits as the SNAIL [7,10,14,16,17]. However, in the case of an arbitrary superconducting dipole, the effective potential energy function U eff (φ 0 ,φ e ) might not have an analytical expression to begin with. To overcome such obstacle, many aforementioned flux-biased loops were arranged in configurations which exploited certain symmetries and approximations to grant an analytical expression for U eff .
In this section, we show how to express (30) without an a-priori knowledge of the function U eff (φ 0 ,φ e ). As described in section II, the equilibrium phase of each element in the loop is an analytical function ofφ F . Consequently, the n-th order potential energy expansion coefficient of the x-th element in the loop reads where U x (φ) is the potential energy function of the element x andφ x (φ F ) is the relation between the phase of the element x of the loop and the phase of the free JJ. The individual expansion coefficients can be combined to return those of the effective potential energy u n in the following way. As a first step, it is useful to describe how the two arbitrary nonlinear inductors A and B in Fig. 6(a) combine in series and in parallel, see details in appendix A. Given the expansion coefficients a n and b n , their parallel combination is described by the expansion coefficients u p n = a n + b n , while the series combination can be expressed as u s n = S n (⃗ a n , ⃗ b n ), where ⃗ a n = (a 0 , a 1 , . . . , a n ) and ⃗ b n = (a 0 , a 1 , . . . , a n ) are vectors whose components are the potential energy expansion coefficients of A and B up to order n, and S n is a rational function of such components. The function S n is derived in Appendix A and can be extended to an arbitrary number of elements, thus can be applied to both branches of the loop in Fig. 6 (b). Consequently, the expansion coefficients of the left and right branches of the dipole can be expressed as a function ofφ F (36) The expansion coefficients of the dipole read Notice how the relations (13) were obtained with the reference directions for phases and flux in Fig. 1(d), while to assembly the dipole potential energy expansion coefficients u n it is more practical to use the reference directions in Fig. 6(b). For a consistent computation of u n , the functionsφ x (φ e ) related to the left branch of the dipole in Fig. 1(d) acquire a negative sign which only affects the odd expansion coefficients. This consideration gives rise to the (−1) n factor in the first line of (36). The relation between the expansion coefficients (37) and external fluxφ e can be represented as a parametric relation inφ F , similarly to the arbitrary branch CPR (15). Such description can also be applied to the parameters of Hamiltonian (27), which are functions of the expansion coefficients u n according to expressions (31) and (32).

V. OUTLOOK
Under the one-body approximation for a generic fluxbiased superconducting dipole, the method discussed in this manuscript makes it possible to express the effective Hamiltonian expansion coefficients (31) and (32) as well as the external flux in the first line of (15) as analytical functions of the dipole design parameters and the free JJ equilibrium phase drop. This description will be at the core of a superconducting quantum circuits optimizer/synthesizer based on a gradient descend algorithm, which is currently under development. For instance, when designing a single-minima effective Hamiltonian to implement multi-photon parametric processes, one could specify constraints on the Hamiltonian parameters ω (31) and g n (32). For a given dipole, the optimizer will be able to find the set of optimal parameters for the inductive elements within the loop, the optimal value of shunt capacitance C opt and the optimal equilibrium phase for the free JJφ Fopt . From this values, the correspondent value of external fluxφ e (φ Fopt ) can be computed from the second equation in (14).
In scenarios where the one-body approximation fails to describe the physical effects of interest of an arbitrary dipole, the equilibrium points computed with our method can be used as an input to more sophisticated algorithms. In particular, a recent proposal [29] demonstrated how a variational tight-binding method can be applied to compute the properties of a large, multi-minima fluxbiased superconducting circuit, accounting also for the junction capacitances. Once the expansion coefficients around each minima of the circuit are known, this approach provides a reduction in complexity with respect to a brute-force diagonalization of the Hamiltonian (2), with a promising fidelity in the estimation of the parameters of interest. A Born-Oppenheimer approximation can also be used to obtain an effective one-body Hamiltonian which accounts for the multi-body nature of an arbitrary dipole [28].
Another natural application of our technique would be to model the repercussion of stray linear inductors and fabrication uncertainties on the Hamiltonian expan-sion coefficients of a superconducting dipole. Indeed, in the current approach the effective Hamiltonian (27) can hardly be computed in presence of stray inductors and asymmetry in arrays of JJs without recurring to the N+1 body Hamiltonian (2). Instead, our method provides a useful shortcut to correctly describe these scenarios.
We are currently characterizing experimental devices with engineered asymmetries in arrays of JJs to validate the predictions of our method, with promising results. These will be presented in a future work focused on the Hamiltonian engineering capabilities enabled by the theory developed in this manuscript.

VI. CONCLUSIONS
In this article, we have investigated the equilibrium properties of an arbitrary flux-biased superconducting loop. In particular, we demonstrated how the relation between the equilibrium points of a superconducting loop and its flux-bias can be expressed analytically as a parametric relation. We also derived a set of rules to count the number of local minima of the loop Hamiltonian. This approach can be generalized to circuits with more than one loop. As an immediate application of our technique, we showed how to compute the effective Hamiltonian Taylor expansion coefficients for an arbitrary flux-biased superconducting dipole shunted by an external capacitor. This enables quantitative analysis of yet-unexplored circuit topologies, overcoming the symmetry constraints and approximations which have limited the variety of devices investigated so far in literature. Our method is also suitable to implement a hardware-level Hamiltonian optimizer based on iterative algorithms: constraints on the quantities of interest, gradients and Hessians can all be expressed as analytical functions of the electrical parameters of the circuit. We believe that the method reported here could play an important role in the development of the next-generations of superconducting quantum devices, providing a shortcut towards the physical understanding, modeling and optimization of advanced fluxbiased Josephson circuits. The results from Section IV rely on the capability of computing the potential energy expansion coefficients of two arbitrary nonlinear inductors arranged in either parallel or series. Here, we show how to compute such coefficients as a function of those of the individual inductors. In general, for both configurations in Fig. 7, the node phase φ, as well as the inductors phases φ A and φ B are related to their voltages, respectively V , V A and V B via the relations whereφ,φ A andφ B are, in our context, the equilibrium phases imposed by the external DC flux-bias to the node and to each inductor. To simplify the notation in the following treatment, we introduce the perturbations around the equilibrium for each phase, namelỹ Consequently, the expanded potential energy functions for the parallel and series configuration read while those of the nonlinear inductors A and B arẽ Note that the expanded potential energy functions are marked with a ∼ symbol to distinguish them from the non-expanded ones. The latter are used in Appendix B to characterize an arbitrary dipole in presence of a timedependent external flux.

Parallel expansion coefficients
To compute the expansion coefficients u p n in (40), we notice how, with reference to Fig. 7(a), A and B being in parallel implies V A = V B = V . Consequently, from Eq. (38) and (39) results that all the phase perturbations are the same in the parallel configuratioñ Thus, the parallel potential energy function reads Computing the n-th order derivative of this last identity, and evaluating it forφ = 0, results in the expression (34) from the definitions in (40) and (42). The generalization to a parallel configuration of an arbitrary set of nonlinear inductors {A, B, C, . . . } trivially reads

Series expansion coefficients
The series configuration of A and B in Fig. 7(b) is described by the set of constraints where the first identity states the voltage conservation V = V A + V B , while the second identity imposes current conservation arising from the series arrangement. The phase perturbations of each inductor can be expanded as a function of the node phase as are the nonlinear participation ratios [30] of A and B. We now define the algebraic rules for the expansion coefficients dp An dφ where we implicitly mean that the derivative is taken before evaluating the corresponding functions forφ = 0. This abuse of notation will ease the elaboration of the S n function. By computing the n-th order derivatives with respect toφ of the first line of (46), and combining them with the expressions (47), we obtain the constraints on the participation ratios The currents through each nonlinear inductor are related to their potential energy functions via the relations Combining these last expressions with the constraints in the second line of (46), and computing the derivative with respect toφ, we obtain which, together with the first of (49), provides the closed expression for the linear participation ratios Higher order participation ratios can be retrieved by applying the algebraic rules (48) to the identity (51). Their expressions up to third order are We can now retrieve the series potential energy expansion coefficients u s n as a function of a n and b n . The series potential energy function can be expressed as from which the expression of u s 1 can be retrieved as Applying the algebraic rules (48) recursively to this last identity, and keeping in mind the constraints (49), u s n can be easily computed. Here we show such coefficients up to fifth order Notice how these expressions are invariant under the swap of A and B, as the properties of their series are independent from the order of arrangement. Consequently, they can be computed, for instance, just for A and then trivially extended to include B as well. The rational functions in (56) define the series combination S n in (35), and can be extended to an arbitrary set of series nonlinear inductors {A, B, C, . . . }. Said x n the n-th order expansion coefficient of the x-th nonlinear inductor in the series, its linear participation ratio p x1 can be expressed as where l, m ∈ {a, b, c, . . . }. Higher order participation ratios for the x-th element can be computed from (57) by applying the algebraic rules (48). As an example, the second-order participation ratio of the x-th element in an arbitrary array reads For an arbitrary array, the expansion coefficients (56) up to fifth order generalize as which, together with p xn , define the rational function S n for an arbitrary array in (36). When the expressions (60) are applied to the arm of a loop which doesn't contain the free JJ, they don't present any singularity as a function ofφ F . This is easily demonstrated by noticing that x 2 (φ F ) > 0 for any constrained element, and the denominators of u s n are multinomial functions of x 2 with positive coefficients. On the other hand, when a branch contains a free JJ, u s n can have singularities if the branch is multi-valued.

APPENDIX B: EXPANSION COEFFICIENTS WITH TIME-DEPENDENT EXTERNAL FLUX
In this section, we analyze the properties of a capacitively shunted superconducting dipole in presence of a time-dependent external magnetic flux. Thus, we introduce the time-dependent component of the external magnetic fluxφ e (t), such that the total flux φ e (t) threaded to the loop forming the dipole reads φ e (t) =φ e +φ e (t).
We will consider the case where max {φ e (t)} ≪ 2π, thus the time-dependent component of the external magnetic flux can be treated perturbatively. Introducing the maximum frequency ofφ e as ω max φe , the analysis can be performed with two different approaches. In the case where ω max φe ≪ ω, where ω is the natural frequency of the harmonic oscillator associated to the capacitively shunted dipole as in Eq. (31), it is reasonable to assume that the effective potential energy (26) acquires a slow time dependence, corresponding to an adiabatic modulation of its minima. This is a realistic assumption, for instance, to model the effects of low-frequency flux noise on the device. In the case where ω max φe ≈ ω, the external magnetic flux can strongly modify the dynamics of the system, for instance, by driving parametric processes. As a consequence, the effective potential energy (26) acquires an additional functional dependence onφ e (t). We will refer to this scenario as "flux pumping". In the following treatment, we will refer to the circuit in Fig. 7(a).

Sensitivity to low-frequency flux noise
In realistic experimental scenarios involving DC fluxbiased superconducting loops is important to account for the presence of a non-deterministic component of the total magnetic flux. Such component can be modeled as a stochastic signal with a low-frequency power spectral density, typically 1/f noise [31], which can arise from both magnetic impurities in the vicinity of the loops and room temperature electronics providing the DC flux-bias to the device. A major detrimental effect of such noise is a slow, random modulation of the effective Hamiltonian parameters in Eq. (30), which can cause dephasing events in flux-biased superconducting quantum devices. The first-order sensitivity of the Hamiltonian parameters to low-frequency flux noise can be evaluated by computing their derivatives with respect to the external DC flux φ e . Here, we show how to compute such derivatives and relate them to the potential energy expansion coefficients of the nonlinear inductors A and B in Fig. 7(a). The sensitivity to slow-flux modulations of the natural frequency in Eq. (31) can be expressed as while that of the n-photons interaction rates in Eq. (32) reads where Inserting this last expression in Eq. (63) results in The derivative of the potential energy expansion coefficients with respect toφ e can be expressed as By noticing that where U A (φ A ) and U B (φ B ) are the potential energy functions of the left and right branches of the dipole, respectively, we derive the algebraic rules for computing the derivatives of the expansion coefficients with respect to the external DC flux-bias Inserting these last expressions in Eq. (66) results in which can be further processed to explicit the dependence on the free JJ equilibrium phaseφ F as Note how this last expression was obtained with the reference directions for phases as in Fig. 7(a). As a consequence, the functionφ A (φ F ) has opposite sign with respect to the one that can be obtained following the method for computing the equilibrium phase drops in Section II, which relies on the reference directions for phases as in Fig.2(a). The first-order sensitivity to slowflux modulations of the dipole' expansion coefficients in Eq. (70) can be inserted in the expressions (62) and (65) to retrieve analytical expressions for dω dφe and dgn dφe . Note how the latter are analytical functions of the free JJ equilibrium phase drop and the electrical parameters of the circuit, thus their relation to the external DC flux-bias φ e still has the form of a parametric curve of curvilinear parameterφ F . The n-th order sensitivity to slow-flux modulations can be obtained by computing the n-th order derivative with respect toφ e of expressions (62), (65) and (70) and applying the algebraic rules in (68).

Flux pumping
In this section, we analyze the case of a fast timedependent flux driveφ e (t) applied to the nonlinear dipole in Fig. 7(a). The potential energy of such driven dipole can be expanded in Taylor are the expansion coefficients of order n + l of the driven dipole potential energy function U (φ,φ e ,φ e ). Following the results in [32][33][34], the flux-driveφ e has to be allocated among the nonlinear inductors A and B according to the spatial distribution of the magnetic vector potential on the device. In particular, the phase drops φ A and φ B across the nonlinear inductors in Eq. (38) acquire a functional dependence onφ e of the form where α is the allocation factor that can be computed following the procedure detailed in [33]. Note how, according to these last expressions, the difference between the phases across B and A reads resulting in the expression (61) for the total magnetic flux threaded to the loop. The last identity can be demonstrated by expressing the second equation in (11) with the reference directions for the phase drops in Fig. 7(a). With the expressions in (73), the driven dipole potential energy function can be expressed as U (φ,φ e ,φ e ) = U A (φ − αφ e +φ A ) thus the dipole expansion coefficients in (72) read where are the driven expansion coefficients of the nonlinear inductors A and B, respectively. From the functional dependence of U A and U B on the right side of expression (75), the driven expansion coefficients in (77) resulting in a compact final expression for the driven dipole expansion coefficients (72), which reads u nl = (−α) l a n+l + (1 − α) l b n+l .
This last expression can be inserted in the driven dipole potential energy Taylor expansion (71) to obtain an expression as a function of the bare expansion coefficients of the nonlinear inductors A and B. The sensitivity to lowfrequency flux noise of the driven dipole potential energy expansion coefficients in (79) can be computed following the same treatment explained in the previous section as