Demonstration of long-range correlations via susceptibility measurements in a one-dimensional superconducting Josephson spin chain

Spin chains have long been considered an effective medium for long-range interactions, entanglement generation, and quantum state transfer. In this work, we explore the properties of a spin chain implemented with superconducting flux circuits, designed to act as a connectivity medium between two superconducting qubits. The susceptibility of the chain is probed and shown to support long-range, cross chain correlations. In addition, interactions between the two end qubits, mediated by the coupler chain, are demonstrated. This work has direct applicability in near term quantum annealing processors as a means of generating long-range, coherent coupling between qubits.


I. INTRODUCTION
Superconducting quantum information platforms have reached a level of maturity where tens of individual qubits, comprising a computational device, can provide proof of principle demonstrations of quantum simulations, quantum algorithms, and basic error correction functionality [1]. As these devices, and the tasks they seek to address, scale in size and complexity, so does the need for realizing qubit networks with increased dimensionality and expanded connectivity. These two desired features of future quantum processors prompt the development of long-range, qubit coherence preserving interactions [2,3]. Quantum spin chains have been proposed as an effective medium for qubit interactions with these desired properties [4][5][6][7][8][9]. In this article, we explore the possibility of long-range interactions supported by quantum spin chains for superconducting qubits [10][11][12]. This architecture has direct application in recently proposed quantum annealing platforms based on superconducting capacitively shunted flux qubits [13,14], rf-SQUIDs [15], fluxmon qubits [16], and fluxonium qubits [17].
Quantum annealing is emerging as a promising paradigm for near term quantum computing [18][19][20][21]. An initial Hamiltonian, whose ground state is straightforward to prepare, is transformed continuously to the problem Hamiltonian. The prepared state of the problem Hamiltonian is located in the vicinity of the true ground state and represents a useful solution of the optimization problem. In the limit of weak coupling to the environment, adiabatic quantum computing has been shown to be immune to dephasing in the energy basis, making it a particularly attractive candidate for near term, noisy quantum computing platforms [22]. Commercial quantum annealers, based on superconducting Josephson flux qubits [23][24][25], have recently become available to the larger community and are beginning to make their mark as a valuable research tool, see e.g., Refs. 26-28. There are strong motivations for improving upon the performance of quantum annealing processors [29], in particular with respect to how their constituent qubits interact with one another. Increasing both the graph dimensionality of qubit networks [2,30], and improving connectivity [31], the degree to which the qubits are coupled to one another, would greatly reduce physical hardware overhead by increasing the types and sizes of optimization problems that can be natively embedded. Existing quantum annealing processors based on superconducting qubits possess either nearest neighbor [32] or a combination of inter-and intra-unit cell interactions [33] between qubits. Commercial annealers, possessing this combination of inter-and intra-unit cell interactions, currently rely on minor embedding [34,35], a procedure of extending logical qubits over multiple physical qubits to implement problems that require higher dimensionality or connectivity than the processor's hardware natively allows.
As each connection made to a qubit introduces additional noise and decoherence channels, expanding qubit connectivity in quantum annealing processors must be balanced against the need to maintain the qubits' coherence properties. Developing quantum annealing processors that support improved qubit coherence times would allow greater functionality in computation. Higher precision flux control, afforded by improved coherence, is required by many computational problems of interest [36][37][38]. In general, just how much of a computational advantage greater qubit coherence provides in quantum annealing processes is itself an open scientific question [39,40]. Furthermore, more coherent quantum annealers will enable diabatic annealing protocols that require a greater degree of qubit coherence throughout the annealing process [41][42][43].
These two, often competing, improvements -creat-ing qubit networks with higher dimensionality and expanded connectivity and maintaining qubit coherencecall for further development of long-range qubit interactions. One proposed scheme that accomplishes this dual need is utilizing spin chains as the qubit interaction medium [4,6,7]. Gapped spin chains, in the context of semiconducting quantum dots [5,8,9], have been proposed to support long-range, Ruderman-Kittel-Kasuya-Yoshida (RKKY) type qubit interactions [44]. Recent progress in this direction includes a demonstration of adiabatic quantum state transfer along a linear array of four electron spin qubits [45]. In addition to the possibility of supporting coherent coupling between two distant qubits, the spin chain architecture lends itself to higher connectivity schemes. Multiple qubits can be simultaneously interacting with a single 1-D chain [5]. Additionally, paramagnetic trees, formed by spin chains forking into multiple paths, offer another possible scheme for higher qubit connectivity [10][11][12]. This work demonstrates the viability of this coupling scheme in the context of superconducting Josephson qubit hardware. In the following we discuss long-range coupling mediated by quantum spin chains in a hardware independent fashion. This is a more natural language to describe long-range coupling as a consequence of the system's underlying quantum phase transition [46][47][48][49]. Following this discussion, we demonstrate a realization of the quantum spin model with superconducting circuits. To accomplish this, we design a system of two qubits, coupled together through a chain of seven spin units. The spin chain, shown in Fig. 1, is realized by a one-dimensional array of seven tunable rf-SQUIDs [13,[50][51][52][53][54] inductively coupled to their nearest neighbor through the SQUIDs' main loops. Each end coupler is inductively coupled to a tunable, capacitively shunted, superconducting flux qubit [55][56][57][58]. Finally, to illustrate the viability of mediating long-range, coherent qubit interactions with our device, we characterize the non-local susceptibility of the coupler chain, demonstrate long-range qubit-qubit interactions, and identify the parameter region where both long-range correlations exist and the detrimental effects of low frequency flux noise are negligible.
The Hamiltonian for the quantum spin chain is the one-dimensional Ising model. Incorporating the two end qubits, it can be written as with and In the previous equations, ∆ qi /2 (∆ ci /2) and qi /2 ( ci /2) are the transverse and longitudinal components of the qubits' (couplers') spin while J cici+1 and J qicj represent the coupling strength between adjacent coupler units and between qubits and their nearest coupler unit. For the remainder of the article, we will assume the coupler units are operated homogeneously, that is ci = c , ∆ ci = ∆ c and J cici+1 = J cc . Virtual excitations of the coupler chain can be integrated over to derive an expression for the couplerchain-mediated effective qubit-qubit interaction strength, J eff q1q2 . By considering the qubit-adjacent coupler unit interaction, J qicj , to be a weak perturbation to the coupler Hamiltonian, the interaction energy can be calculated to second order as the shift of the ground state energy of the coupler Hamiltonian. As the operating temperature of the device will be much less than the coupler chain excitation energy, it is reasonable to assume that the coupler chain remains in its ground state and the cross chain interactions are supported by virtual excitations [6,9,59]. This is reminiscent of the RKKY interaction whose longrange interaction between magnetic impurities is mediated by virtual excitations of conduction electrons above the Fermi surface [44].
By taking these above stated approximations into account it is possible to derive an expression for the chain mediated effective coupling strength between the end qubits [59]. The effective Hamiltonian is where Ω c is the energy gap between the coupler chain ground and first excited state and |0 c represents the collective ground state of the unperturbed seven unit coupler chain. Note that an exact expression for the effective coupling between qubits contains the integrals of frequency dependent connected coupler correlation functions. The ground state connected correlation function in Eq. (5) is an approximation assuming a large excitation gap, Ω c , and that the coupler chain excitation frequencies are sufficiently degenerate [7,8]. This approximation is strictly valid in the case where the transverse field on each coupler unit is much larger than the exchange interaction between coupler units. Writing this expression in terms of the zero temperature, bulk susceptibility of the response function in the Lehmann representation, χ c1c7 , [7] the effective interaction can be expressed as The main objective of this work is to measure the quantity χ cicj as a function of J cc /(∆ c /2), the ratio of the inter-coupler longitudinal coupling strength, proportional to σ z ci σ z ci+1 , to the individual coupler unit transverse field strength, oriented along σ x ci , for the homogeneously tuned chain. Long-range coherent coupling a. b.

FIG. 2.
Single Coupler Behavior. a) The ratio of the Josephson inductance to the geometrical inductance, βc, dictates the shape of the potential energy of the tunable rf-SQUID coupler circuit. When the geometrical inductance dominates, βc 1, the potential energy landscape is approximately harmonic. When the Josephson inductance dominates, βc 1, the energy landscape becomes double-welled with each minimum representing oppositely circulating current states. Due to the large energy barrier between the states, moderate changes in fz do not change the current state of the circuit. The coupling is optimized when the geometrical and Josephson inductances are approximately equal. This results in a wide, shallow energy minimum where even slight changes in fz can induce strong fluctuations between the oppositely circulating current states of the coupler circuit. b) Both the single coupler transverse field, ∆c/2, and the intercoupler interaction energy, , are displayed as a function of the coupler fx when fz = Φ0/2. These parameters are calculated in single coupler simulations and then transcribed into spin model parameters. The equality of these two terms appearing in the transverse field Ising model for fx 0.14 Φ0 signals the location of the quantum critical point, in the vicinity of which we expect long-range correlations to emerge. In addition, the dependence of βc is displayed as a function of the coupler's fx for fz = Φ0/2. By design, the optimum coupling point, βc ≈ 1, coincides with the coupler fx value where we expect critical behavior in the coupler chain.
becomes possible when the spin chain is tuned to the vicinity of its quantum critical point [46,47]. In the case presented here, this occurs when the strength of the transverse fields of the coupler spins and inter-unit longitudinal coupling energies between the nearest-neighbor coupler spins become comparable. The coupler chain susceptibility is determined by measuring the response of the longitudinal fields of the coupler units along the chain when a small longitudinal field, δ , is applied to the end coupler unit. It is shown that the response truly becomes long-range, that is entirely cross chain, for J cc /(∆ c /2) 1, where the system approaches and enters its ordered phase.
The one-dimensional transverse field Ising spin model can be realized by multilevel superconducting Josephson circuits. With the assumption of negligible state occupation of higher levels, the two lowest energy levels of the circuit define the qubit subspace where the transverse and longitudinal components of the unit's spin can be determined as a function of the applied magnetic flux. The individual qubit and coupler circuits utilize inductive couplings for implementing the inter-unit interactions J q1c1 , J c7q2 , and J cici+1 . Coupling of this type for single unit coupler circuits has been demonstrated in flux qubits [13,[50][51][52], phase qubits [53,54], and fluxmons [16]. The design choice of independent coupler circuits, as opposed to direct coupling between qubits, is particularly appealing for use in annealing processors where it is necessary to independently control the qubit properties and coupling strengths. This single unit method of identifying both the qubit and coupler's spin components is not as universally applicable as more general means such as the Schrieffer-Wolff transformation [60], particularly in the strong coupling regime. However, as we will restrict our analysis to the weak coupling limit, the results of the two methods coincide [16,[61][62][63]. With these assumptions, the behavior of the physical device can be mapped to the one-dimensional, transverse field Ising spin model.
Each coupler circuit can be approximately characterized by its susceptibility, χ, which is the change in current induced by a biasing flux. Assuming the coupler remains in its ground state, this is equivalent to the curvature of the ground state energy with respect to the flux in the coupler's main loop.
In Eq. (7), I z c represents the ground state expectation value of the current in the coupler's main loop and E 0 c is the ground state energy of the coupler unit. The character of these two quantities, I z c and E 0 c , is determined by f x , the magnetic flux in the coupler's small loop. The coupler circuit's susceptibility, χ, is optimized as the unit's β c ≡ L c /L eff = 2πL c I (c) c /Φ 0 ≈ 1, where the coupler's local potential minimum is highly sensitive to biasing flux. As shown in Fig. 2, this occurs in the same f x region where J cc /(∆ c /2) ≈ 1, a design choice made to optimize the generation of long-range correlations across the device. In addition, the device can be operated in a regime such that the coupler's minimum excitation energy, larger than 5 GHz, is much greater than the temperature of the system, approximately 400 MHz, the strength of the qubit-coupler interaction, which is below 1 GHz, and the typical qubit excitation frequency, approximately 2 GHz. This ensures the ground state properties of the coupler dictate its behavior, entanglement between qubits mediated by the coupler is supported, and fast (coupler) and slow (qubit) modes can be separated to preserve the qubit subspace.
The same notions of susceptibility and design constraints can also be applied to a long but finite chain of couplers. The current induced in coupler j when a flux is applied to coupler i is expressed as χ cicj , the inter-chain susceptibility. The design constraints are slightly more involved for the chain when compared to a single coupler. For example, the length of the chain has a closing effect on the size of the gap as the fundamental mode The displayed signal is the difference in effective fz = Φ0/2 point when the source unit is placed at fz = Φ0/2 ± 20 mΦ0. As the strong coupling between the tunable resonator and the qubit (coupler) makes it challenging to model the resonator response, we therefore resort to image processing techniques to determine the effective symmetry point of the qubit (coupler) unit. This sets the uncertainty of the extracted flux signal to be on the order of the pixel size of the scan, 2.4 mΦ0. In the case of both the coupler only and qubit to qubit measurements, the cross chain signal becomes non-zero in the vicinity of the coupler chain's predicted critical region.
frequency of the chain decreases with length. This introduces a trade-off between the physical range of interaction and the need to preserve the excitation energy gap of the chain [64].
Taking the same concept of single coupler susceptibility to hold for the coupler chain, as well as adhering to the extra constraints introduced by the many-body chain system, it is possible to construct the form of the effective qubit-qubit interaction, Eq. (6), in terms of circuit parameters. Assuming the coupler gap is much greater than the qubit working frequencies allows one to separate the device spectrum into 'slow' qubit-like states and 'fast' coupler-like states. Invoking the Born-Oppenheimer approximation restricts the coupler spectrum to its unperturbed ground state [65]. Further restricting our analysis to the weak coupling limit, M qc /L c 1, allows us to write the qubit-chain interaction, in general a complicated nonlinear quantity, as an inductive interaction between qubit currents [16,66], given by The final line of Eq. . The symbol I p ci refers to the persistent current of the i th coupler, which, when operated at f z = Φ 0 /2, is simply the current dipole moment 0|Î z ci |1 [61]. Now that the long-range, effective interaction between the two qubits is expressed in terms of circuit parameters, it is possible to measure the response function, χ c1c7 , of the coupler chain in a quantitative manner. This task is accomplished by performing two similar measurements. Firstly, only the behavior of the coupler chain units are considered. This measurement is performed with both flux qubits placed at a magnetic flux bias operating point where the circuit has its maximum transition frequency which is much greater than the operating frequencies of the remaining chain units. This decouples the qubits from the coupler chain dynamics and allows the coupler susceptibility to be characterized. Secondly, with knowledge of the chain susceptibility, the two flux qubits are brought into an interacting flux operating point and cross chain qubit-qubit interactions are demonstrated. Finally, to quantitatively measure the effective qubit-qubit interaction strength, J eff q1q2 , we revisit the coupler chain only measurements in more detail to extract the cross chain susceptibility, χ c1c7 . We find that our measurements of J eff q1q2 through susceptibility measurements agree well with full device simulations of the qubit-qubit spectral line splitting.

II. RESULTS
The coupler chain device consists of two capacitively shunted, tunable flux qubits and seven tunable rf-SQUIDs, all equipped with individual readout resonators. Device characterization, circuit parameter extraction, and wiring details are described in the Supplementary Materials [59]. The device is fabricated using the architecture described in Ref. [67], and consists of two separate chips -called the qubit layer and the inter- for Coupler 1, the effective coupling between the two qubits is calculated via Eqs. (9) and (10) and displayed. This is compared to J eff q 1 q 2 calculated from the splitting of the otherwise degenerate qubit transitions as predicted by full device simulations incorporating reasonable levels of flux noise. The strong agreement between the two different measures of J eff q 1 q 2 indicate that the optimal operating region of the chain is coupler 0.15 fx ≤ 0.18 Φ0, where there is significant coupling strength and, as shown in Fig. 5, both the detrimental effects of flux noise and the qubit-coupler state mixing is negligible. poser layer. The interposer layer, seated on the device package's printed circuit board cavity and wirebonded to the exterior control lines, holds the flux bias lines. The qubit layer, hosting the qubits, couplers, resonators, and co-planar waveguide, is indium bump bonded atop. The indium bumps provide structural stability, common ground paths between layers, and a conduit for microwave signals originating on the interposer layer, running through the bumps, and continuing on the qubit layer. This device environment allows greater flexibility than planar devices for distributing flux bias lines, represents a step towards full 3-D integration, and supports an electromagnetic environment suitable for quantum annealing controls [59].
Each unit, qubit or coupler, possesses a meandering resonator terminated in an rf-SQUID for purposes of readout and calibration. These resonators, when their terminating rf-SQUID is biased to a flux sensitive region, act as magnetic flux detectors, capable of discerning the qubit or coupler unit's persistent current state. When the terminating rf-SQUID is biased to its flux insensitive operating point, the resonator is exclusively sensitive to the unit's energy level occupation through the resonatorunit dispersive interaction [68]. Being able to operate in these two modes alleviates the need for multiple readout structures, further freeing up space on chip.
Gaining full flux control of a device of this size is a difficult task for a number of reasons. Complete individual control of each unit requires 27 flux bias control lines corresponding to the 27 Josephson flux loops. Current in one control line provides magnetic flux for its target Josephson loop but also couples to nearby loops. Hence, it is necessary to determine the full 27×27 element mutual inductance matrix before one can expect adequate control of this device. In addition, these inductive elements need to be determined while in the presence of spurious interactions between units. Strong inter-unit interactions can easily mask the linear line-loop inductive interaction. In order to address these points, scalable, device independent, automated methods have been developed and implemented to characterize the bias line to circuit flux inductive matrix to within acceptable errors for device control [69].
To explore the behavior of long-range qubit interactions mediated by the coupler chain, it is necessary to characterize the inter-coupler susceptibility, χ cicj . To isolate the coupler chain dynamics, we first flux bias the two end qubits to their high frequency, uncoupled state. Every coupler unit is operated such that its main loop is flux biased at one-half magnetic flux quanta and its small loops are uniformly biased with f c x . In this configuration, Coupler 7's f z is swept across its half quanta point for a range of uniformly biased coupler f c x values. Instead of directly measuring the current response in the target unit, we observe the shift of the target unit's effective main loop half-quanta point [59]. As the unit's main loop half-quanta flux operating point corresponds to its minimum transition frequency, the dispersive in-teraction with the unit's resonator provides an accurate determination of the unit's effective half-quanta point in the presence of strong inter-unit and unit-resonator interactions.
The magnitude of the induced fluxes for each coupler unit are displayed in Fig. 3. Recall that f c x simultaneously controls the unit's transverse field, ∆ c , in the Ising spin model picture as well as the magnitude of the unit's persistent current, thus the longitudinal coupling strength, J cc = M cc I z ci I z ci+1 [59].
For larger values of f c x , the flux propagation signal attenuates over short length scales, up to a few coupler units. As shown in Fig. 2b, this is the regime where the transverse field dominates and the chain system is in its paramagnetic state. For f c x ≤ 0.18 Φ 0 , long-range correlations are supported across the entire chain. This is where the critical region of the underlying Ising spin model is expected to be located.
The critical region of the Ising spin model is determined by single coupler properties. To further validate our results, full device simulations of the experimental protocol were performed. To perform simulations of a device with such a large Hilbert space, a hierarchical scheme is employed [71]. First, the low energy spectrum, eigenstates, and other operator eigenvalues are computed for individual units. These eigenstates then form the initial basis for computing the energy spectrum of two and three unit systems comprising subsections of the full device. Finally, the subsections are appropriately coupled together and the full device energy spectrum, eigenstates, and relevant operator eigenvalues are calculated. At each step of this procedure, necessary mode occupation numbers are included in the calculation such that the lowlying energy and other operator spectra are well converged. Simulating the identical procedure as the measurement protocol allows us to track the target unit's effective main loop half-quanta point yielding results that match well with the experimental outcome [59].
With the derived information on the chain susceptibility, we now place the two end qubits at flux operating points where it is possible to couple to the adjacent chain unit. Qubit 1, the target qubit, is operated such that its main loop is flux biased at one-half magnetic flux quantum and its smaller x-loop is flux biased such that its transverse field has strength ∆ q1 = 2.3 GHz, approximately where the qubit's potential becomes doublewelled. Qubit 2, the source qubit, is placed at its minimum ∆ q2 10 MHz, deep in its double well regime, and its flux bias f z is swept across its one-half magnetic flux quantum point. This measurement protocol is repeated for the different coupler chain operating points described in the coupler chain susceptibility experiment.
The results of the qubit-qubit interaction experiment are shown in Fig. 3. This figure displays the magnitude of the flux signal propagating along the coupler chain and ultimately into the opposite qubit. These results agree well with full device simulations [59] of the equivalent protocol. As shown in Fig. 3c & d, long- FIG. 5. Effects of Noise. a) Full device simulation illustrating the effects of low frequency flux noise on the qubit-like energy levels. Displayed are the means and standard deviations of the full device energy levels compiling ten separate calculations with the described random flux operating point offsets. The two lowest energy levels are identified as the two initially degenerate qubit levels, |+ − and |− + . In the limit of no coupling, or large coupler fx, they reproduce single qubit behaviour. As the coupler fx is lowered and cross chain coupling becomes significant, the previously degenerate qubit levels, shown in bold, split in the presence of cross chain coupling. Also shown, in lightly faded color, are the few next higher energy levels of the full device consisting of a mixture of coupler and higher qubit levels. Note that the coupler chain and the qubit levels have become comparable in frequency by coupler fx 0.15 Φ0. In this region, what were initially qubit-like energy levels, are now dressed by the coupler levels and the splitting of these two energy states no longer represents the effective qubit-qubit interaction strength [70]. b) Shown is the linewidth of the lower qubit-like level, calculated as the standard deviation of the transition frequency over multiple simulation runs in the presence of realistic flux noise. There is a region, fx = 0.15−0.18 Φ0, where significant cross chain coupling is present but where the effects of flux noise do not significantly broaden the lower qubit-like level's linewidth. For fx < 0.15 Φ0, the calculated linewidth suggests the coherence times of the qubit have deteriorated to the nanosecond scale. chain interactions become supported at approximately f c x ∼ 0.15 − 0.18 Φ 0 in both the coupler chain susceptibility and long-range qubit interaction experiments.
Furthermore, the full results of the coupler-only susceptibility measurements can be used to predict the strength of the effective qubit coupling, J eff q1q2 , mediated by the chain. Equations (9) and (10)  , determines the effective qubit interaction strength.
Shown in Fig. 4 are the effective one-half magnetic flux quantum points of Coupler 1's main loop as a function of Coupler 7's f z for various homogeneous coupler unit settings. As expected, for the flux regime where the transverse fields dominate, Coupler 1's effective half quanta point is unaffected by the f z of Coupler 7. As the transverse field is lowered and the coupler-coupler longitudinal coupling strength increases, the effect of Coupler 7's f z on Coupler 1 becomes more pronounced. The slope at the center of these curves, , are determined from single unit simulations [59]. A proper comparison of the qubit interaction strength, J eff q1q2 , derived from both full device simulations and the measured coupler susceptibility, is necessary to confirm that the coupler susceptibility is a valid measure of qubit coupling strength. As illustrated in Fig. 5, The effective qubit-qubit interaction strength is derivable from full device simulations by noting the energy level splitting between the two end qubits' previously degenerate energy spectrum. This is done by first setting both qubits to have the same finite transverse field and zero longitudinal field. As the cross chain coupling is increased by lowering the f c x of all couplers, the initially degenerate qubit states, |+ − and |− + , develop a splitting that, for weak coupling, is twice the coupling strength 2J eff q1q2 . A comparison can then be made between the qubit interaction strength, J eff q1q2 , as measured in the flux signal propagation experiment, and the qubit level splitting exhibited in full device simulations. As shown in Fig. 4c, the simulated results find strong agreement with the results of our coupler susceptibility measurements in the weak interaction limit. The quantitative divergence of these two results is expected in the strong coupling regime, for the coupler flux biases f c x < 0.15 Φ 0 . In this case, the weak interaction limit assumed in Eq. (8) breaks down and the effective qubit-qubit interaction can no longer be described as a linear inductive coupling. Additionally, in the strong coupling regime shown in Fig. 5, the qubit energy levels become dressed by the coupler levels. Hence, the two lowest energy levels can no longer be identified as purely qubit-like and their spectral distance no longer represents a simple qubit-qubit interaction.
Low frequency flux noise is expected to play a detrimental role in the coherence preserving properties of this long-range interaction. As demonstrated, flux signals can be transported and even amplified in certain coupler flux operating regimes. The question is then, is it possible to identify a flux operating regime for the coupler units where strong, long-range coupling is present, yet the detrimental effects of flux noise are not amplified across the device?
To answer this, full device simulations were performed with realistic values of low frequency flux noise. Flux noise has been measured, across many different platforms and frequencies [58,72,73], to be approximately 1/f α in nature, α ∼ 0.9, with magnitude 1 − 5 µΦ o / √ Hz. For moderate frequency measurements, the effect of this low frequency noise is to effectively add a small random flux offset to the flux operating point of the measurement. This small random flux offset is sampled from a Gaussian distribution whose standard deviation is determined by integrating the noise spectrum over the appropriate frequency range, from measurement repetition rate to pertinent experimental frequencies, as well as accounting for the circuit geometry. This amounts to a typical random flux offset in the tens of µΦ o .
Simulations of this type were performed repeatedly to determine the behavior of the device energy level structure in the presence of low frequency flux noise. As shown in Fig. 5, the energy spectrum of the device is highly susceptible to flux noise in its deeply coupled state. Significant line broadening occurs for the uniformly tuned coupler f c x between 0 -0.15 Φ 0 . This allows us to identify a region of flux operation, coupler f c x from 0.15 to 0.18 Φ 0 where significant long range interactions are present yet the detrimental effects of flux noise are still minimal.

III. DISCUSSION
Quantum annealing processors stand to benefit from higher dimensional qubit networks, expanded connectivity, and improved qubit coherence. Accomplishing this will require long-range qubit interactions that do not degrade qubit behavior. The use of spin chains as a quantum bus is a promising venue for this. Presented here is a preliminary step in this direction in the context of a superconducting Josephson system. To build on this idea, there have been proposals to generalize one dimensional spin chains, capable of entangling end qubits, to both paramagetic trees [10,11] and two dimensional spin networks capable of providing entanglement amongst a perimeter of qubits [11,12,74]. Such architectures could provide a scalable, coherent quantum annealing device with high graph dimensionality.
Susceptibility measurements in quantum systems, such as those performed in this study, have been considered as a possible measure of the system's entanglement [75].
The susceptibility experiment's close agreement with full device simulations, which also demonstrate qubit energy level splitting in the presence of expected noise levels, suggest the qubits can be prepared in an entangled doublet state. In this view, the susceptibility measurements presented here are a consequence and valid measure of a coherent long-range qubit interaction. This view, however, needs further experimental validation. Future experimental work will address measurement of the coherent coupling enabled by this method using spectroscopic characterization as well as adiabatic transfer protocols. Additionally, the detection of entanglement can be augmented by measuring other observables and witnesses for interacting quantum spin systems [76].
In closing, we have demonstrated long-range interactions in a superconducting Josephson spin bus by probing the device's response function. Simulations of the device, which agree well with measured quantities, predict significant long-range interaction simultaneous with satisfactory qubit coherence. This device has immediate application in near-term quantum annealing devices where both long-range and coherent qubit couplings are necessary for quantum computation speedup.

V. DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

VI. COMPETING INTERESTS
The authors declare no competing interests.

VIII. ACKNOWLEDGMENTS
We wish to thank the members of the Quantum Enhanced Optimization/Quantum Annealing Feasibility Study collaboration for various contributions that impacted this research. In particular, we thank Lorenzo Campos Venuti and Cyrus F. Hirjibehedin for helpful discussions on spin chains and device design and Ken Zick and David Ferguson for various discussions on experiments. We gratefully acknowledge the MIT Lincoln Laboratory design, fabrication, packaging, and testing personnel for valuable technical assistance.
The research is based upon work (partially) supported by the Office of the Director of National Intelligence  Cross Chain Susceptibility. a) The effective one-half magnetic flux quantum points for Coupler 1's main loop as a function of Coupler 7's f z for different homogeneous coupler f x settings. The response becomes non-zero in the vicinity of the couplers' f x = 0.15 Φ 0 . The offset from f z = Φ 0 /2 is due to slight mistunings of the couplers' longitudinal fields. The uncertainty in the extracted values is identical to those in Fig. 3. b) The curves in Panel (a), and the corresponding figures for the other units, were then fit to a sigmoid function. Displayed here is the midpoint slope extracted from those fits for all units for different homogeneous coupler f x settings. The error bars were generated by repeating this procedure many times with small offsets picked from a normal distribution with standard deviation of 1.2 mΦ 0 applied to the nominal effective symmetry point values. The error bars are the resultant standard deviation of the large array of slopes extracted from the fit sigmoid function. c) Using the slopes from Panel (b) for Coupler 1, the effective coupling between the two qubits is calculated via Eqs. (9) and (10) and displayed. This is compared to J eff q1q2 calculated from the splitting of the otherwise degenerate qubit transitions as predicted by full device simulations incorporating reasonable levels of flux noise. The strong agreement between the two different measures of J eff q1q2 indicate that the optimal operating region of the chain is coupler 0.15 f x ≤ 0.18 Φ 0 , where there is significant coupling strength and, as shown in Fig. 5 Effects of Noise. a) Full device simulation illustrating the effects of low frequency flux noise on the qubit-like energy levels. Displayed are the means and standard deviations of the full device energy levels compiling ten separate calculations with the described random flux operating point offsets. The two lowest energy levels are identified as the two initially degenerate qubit levels, |+ − and |− + . In the limit of no coupling, or large coupler f x , they reproduce single qubit behaviour. As the coupler f x is lowered and cross chain coupling becomes significant, the previously degenerate qubit levels, shown in bold, split in the presence of cross chain coupling. Also shown, in lightly faded color, are the few next higher energy levels of the full device consisting of a mixture of coupler and higher qubit levels.
Note that the coupler chain and the qubit levels have become comparable in frequency by coupler f x 0.15 Φ 0 . In this region, what were initially qubit-like energy levels, are now dressed by the coupler levels and the splitting of these two energy states no longer represents the effective qubit-qubit interaction strength [70]. b) Shown is the linewidth of the lower qubit-like level, calculated as the standard deviation of the transition frequency over multiple simulation runs in the presence of realistic flux noise. There is a region, f x = 0.15 − 0.18 Φ 0 , where significant cross chain coupling is present but where the effects of flux noise do not significantly broaden the lower qubit-like level's linewidth. For For the derivation of the effective qubit-qubit interaction strength, we assume that the coupler chain is operated in its paramagnetic regime such that its ground state is non-degenerate and the energy gap to its first excited state is sufficiently larger than the temperature of the system, the qubit-coupler interaction strength, and the qubit operating frequencies. As the chain approaches its critical region, where we expect long-range correlations to emerge, the gap does become smaller. But, as shown in Fig. 5 of the main text, the gap still exceeds these other energy scales.
Treating the qubit-chain interaction as a small perturbation on the coupler chain system allows us to write First order perturbation expansion shifts the longitudinal field of the qubits, where J qic adj = J q1c1 or J q2c7 . The effective qubit-qubit interaction, mediated by virtual excitations of the chain out of its ground state, can be calculated to second order in perturbation theory as where Ω ≡ ω 01 , the energy gap between the ground and first excited state. In Eq. (4) it is recognized that there is a band of N = 7 approximately degenerate excited states with gap Ω and finite matrix elements from the ground state. The higher coupler chain levels past this have, at best, exponentially (in energy difference) small σ z ci matrix elements from the ground state and can be neglected. Now the outer product of states, |n c n c | can be written as 1 − |0 c 0 c |. The effective coupling strength then reduces to

II. CIRCUIT MODEL
The Coupler Chain device consists of 7 tunable rf-SQUIDS coupled through their z-loop mutual inductance with the two end couplers coupled mutually to two capacitively shunted flux qubits. Shown in Fig. 1 are schematics of the qubit and coupler design. The qubit and coupler are divided into four and three floating islands respectively. Design work utilized Ansys Maxwell finite element electromagnetic simulations supplemented with simulations run in Sonnet to correctly account for higher frequency inductive effects. Effort was made to reproduce effective circuit parameters found in [1,2]. The results of these simulations are displayed in Table 1.

III. SIMULATION FRAMEWORK
Upon attaining all lumped element circuit parameters from the classical electromagnetic Ansys Maxwell and Sonnet simulations, quantum simulations of the lumped element single qubit and coupler units and full device were performed using the MIT-Lincoln Laboratory developed JJSim quantum circuit simulation package [3]. With this software package, energy spectra and operator matrix elements can be extracted. Single unit simulations are performed by identifying internal modes and evaluating the circuit Hamiltonian in the harmonic level basis for high enough levels of internal modes such that the operator expectation values in the low energy qubit subspace converge.
For full device simulations, individual units were grouped and coupled to one another in a hierarchical structure. The end qubits remain their own two subgroups. The coupler chain was partitioned into groups of couplers (1,2), (3,4,5), and (6,7). These individual groups then underwent convergence tests exploring how many group energy levels need to be included in simulations such that the low energy operator expectation values converge. Each of the 5 groups were then coupled together for the full device simulations. Again, convergence tests were performed at this final level. Finally, to validate our simulated results with our agreed upon energy level structure, energy levels at lower levels of the hierarchy were varied to double check the validity of the simulated results.
The JJsim lumped element simulations were then incorporated into the full simulation workflow, from classical electromagnetic simulations of the physical device to energy spectrum simulations of the derived lumped element circuits, to optimize the device design. Once the device design was finalized, an additional simulation platform calling JJSim was then developed. This enables an efficient means of adding variation in circuit parameter and flux value assignment in order to account for low frequency flux noise and fabrication uncertainties, as well as visual result presentation.

IV. DEVICE DETAILS
This device consists of two highly coherent superconducting tunable capacitively shunted flux qubits and seven similarly designed tunable rf-SQUIDs all equipped with individual readout. These are aluminum devices with Al/AlO x /Al junctions.
This device is hosted in a two-tier environment, illustrated in Fig. 3, similar to the one described in [4]. The interposer tier, sitting in the device package and wirebonded to a PCB which routes signals to the exterior control lines, holds the flux bias lines. The qubit tier, hosting the qubits, couplers, resonators, and co-planar waveguide, is Indium bump bonded atop. Both tiers are built out of a thick silicon substrate layer topped by the high quality Aluminum film. This device environment allows greater flexibility in distributing flux bias lines, represents a step towards full 3-D integration, and supports an electromagnetic environment suitable for quantum annealing controls. The device package is gold-plated copper with aluminum wirebonds and installed within a dilution refrigerator with base temperature of ∼ 20 mK.
Single unit coherence and annealing measurements have been performed in the two-stack environment. Coherence times measured match well with measurements taken on similar qubits in a single planar tier [6]. Single qubit annealing experiments, as shown in Fig. 5, were also performed yielding comparable transition widths as those performed on planar devices [7].
In addition to these encouraging results, the two-stack qubit environment also provides a 'clean' electromagnetic environment. Initial annealing experiments on single-tier qubit environments were hindered by long electromagnetic ring-down times following annealing flux control pulses. As shown in Fig. 4, this could be measured by tracking the flux-dependent frequency of the readout resonator. Settling times for the planar geometry devices could reach 100s of µs. Repeating these measurements in the two-stack environment yields settling times of approximately a microsecond. Our initial investigations into this effect broadly agree with the results of [8], i.e., that the long ring down times can be attributed to charge redistribution of the normal metal ground plane but that the superconducting ground planes on both tiers of the two-stack effectively shield the Josephson elements from this effect.
Each qubit/coupler unit is equipped with its readout resonator. These resonators are quarter wavelength meandering resonators terminated with an rf-SQUID. To model this component we consider a quarter wave transmission line resonator [9] terminated in a variable inductor whose value is dictated by the rf-SQUID properties [10]. Bare resonator frequencies came out slightly higher than design values by 0.5-2%. To investigate this discrepancy, three parameter fits, involving resonator length, junction critical current, and geometrical inductance in one case and a common phase velocity, junction critical current, and geometrical inductance in another were applied to the resonator spectra. In the first case, the resonator lengths came out consistently 2% lower than design values, while in the second case, the phase velocity turned out to be 2.4% larger than design. Both cases showed a similar spread in values of critical currents and geometric inductances. Qubit 2's circuit parameter values were extracted by matching spectroscopy measurements to single unit circuit simulations described above. Qubit 1 suffered an unresponsive capacitive drive line during the experimental run. Displayed in Fig. 6, the data is well fit to the simulation and extracted values agree well with design.
From coupler 'topographical' scans, i.e., probing the coupler resonator at a fixed frequency while varying the coupler's two flux bias currents as shown in Fig. 7, we are able to extract couplers' x-loop junction asymmetry, d = Jx 1 −Jx 2 Jx 1 +Jx 2 . Even small asymmetries of a few percent, due to fabrication imperfections, can appreciably shift the z-symmetry point of the qubit/coupler unit.

V. EXPERIMENTAL METHODS
Both the coupler-only susceptibility measurement and full device qubit interaction demonstration experiments were conducted in similar fashion. Either Coupler 7, for the coupler-only, or Qubit 2, in the full device, acted as the source unit. For uniformly tuned coupler units, all couplers at their z-symmetry point and homogeneously tuned x-flux settings, the source unit's z-flux is swept across its z-symmetry point, causing the source unit to transition from its 'left' to 'right' circulating current state. This acts to shift the effective z-symmetry point of the remaining units in a manner dependent on the coupler units x-flux.
The response of the target units is measured by observing the dispersive shift in the target units' resonator. The target units' resonators are maintained at zero flux or their high frequency point. The dispersive interaction between the unit and resonator can be used to ascertain the unit's z-symmetry point. So, for each source unit's z-flux setting, the target unit's resonator is probed over a range of frequencies for a range of target unit's z-flux settings. In this way, the target unit's effective z-symmetry point can be tracked as a function of the source unit's z-flux setting.
This procedure was applied iteratively starting with the unit adjacent to the source unit, and then continued down the chain. Applying this process iteratively also allowed us to fine-tune the target units' z-flux. As the coupler units' x-flux is lowered and the inter-coupler interactions become stronger, slight mis-tunings of the unit's z-flux away from their symmetry point can shift the effective symmetry point of nearby units. Applying the measurement/fine-tuning procedure iteratively allows the units to be tuned correctly even in the presence of slight imperfections further down the chain. Simulations using the platform described above were used to validate the experimental flux signal propagation results. The f z symmetry point of the target unit is defined as the unit's f z where the ground state expectation value of the z-loop current equals zero. To track the magnitude of the flux signal, we recorded the difference in the target unit's z-symmetry point when the source unit was 20 mΦ 0 on either side of the z-symmetry point. This protocol was followed in both simulation and experiment. Shown in Fig. 8, there is strong agreement between the two.
To measure the chain susceptibility, the entire data set was used, not just the source unit's ± 20 mΦ 0 z-flux points. The response of the target unit's f z symmetry point with respect to the source unit's f z was then fit to a sigmoid function. The slope at the midpoint, To complete the cross chain effective strength calculation, single unit simulations were performed calculating the coupler's ground state expectation value of the z-loop current as well as the qubit's 0| I z |1 matrix element at its f z symmetry point. This measurement of J eff was compared to a full device simulation of the qubit spectral line splitting and is displayed in Fig. 4c of the main text.  First, a small offset in fz is applied. Then the change in fx is applied. b) Shown are the dynamics of the Ising coefficients during the fx journey. Initially the qubit is dominated by the transverse field. As the fx journey continues the transverse field turns off while the longitudinal field increases. There is a minimum gap at 0.68 Φo where these two energy scales are equal and the nature of the qubit ground state transitions between theσx andσz basis. c) Single shot data showing the histogrammed output voltage of 10,000 runs as a function of initial qubit fz. The bi-modal distribution of output voltages correspond to the measuring frequency being on-resonance with the qubit in its 'left' circulating current state and then off-resonance with the qubit in its 'right' circulating current state. d) The output voltages are then thresholded and the ground state probability is displayed as a function of qubit fz. e) The center fz and transition widths are plotted as a function of fx journey time. The shift in center is a reflection of the qubit x-loop junction asymmetry. For a range of longer annealing times, the system thermally exits its ground state at different points along the anneal, shifting the effective fz symmetry point. These annealing pulses were conducting through low frequency lines, acting as low pass filters of approximately 20 MHz. The optimal annealing time of ∼ 1 µs results in a width of ∼ 3 mΦ0, similar to previously reported time and widths on these flux bias lines [5].
FIG. S6. Qubit 2 Parameter Extraction. Two-tone spectroscopy measurements were taken for Qubit 2 for a range of fx values. Using a rejection sampling method, multi-variable fits were performed on various circuit parameter. As shown in the table, all parameters match well with design values.
FIG. S7. Coupler X-Loop Junction Asymmetry. a) A 'topographical' scan of Coupler 3. By performing inversion symmetry analysis on these data sets [11], it is possible to extract the effective Φo/2 points. b) These points are well fit by the x-loop junction asymmetry model. c) Shown are the asymmetry parameters for the coupler units measured approximately 4 months apart during different dilution refrigerator cooldowns.