Automated design of superconducting circuits and its application to 4-local couplers

Superconducting circuits have emerged as a promising platform to build quantum processors. The challenge of designing a circuit is to compromise between realizing a set of performance metrics and reducing circuit complexity and noise sensitivity. At the same time, one needs to explore a large design space, and computational approaches often yield long simulation times. Here, we automate the circuit design task using SCILLA. The software SCILLA performs a parallelized, closed-loop optimization to design superconducting circuit diagrams that match predefined properties, such as spectral features and noise sensitivities. We employ it to design 4-local couplers for superconducting flux qubits and identify a circuit that outperforms an existing proposal with a similar circuit structure in terms of coupling strength and noise resilience for experimentally accessible parameters. This work demonstrates how automated design can facilitate the development of complex circuit architectures for quantum information processing.


INTRODUCTION
The promise of quantum computing to surpass the capabilities of classical computers relies on a robust and scalable underlying hardware architecture. In the context of quantum simulation and annealing, strong coupling, high connectivity, and many-body interactions between the qubits are necessary to accurately and efficiently represent the problem Hamiltonian 1-3 . Superconducting circuits have proven to be a particularly well-suited platform due to their design versatility 4 . Their quantum behavior arises from the interaction of modes that are set by effective inductances, capacitances, and nonlinear Josephson junction elements in the circuit 5 . In this way, it is possible to design a wide variety of qubits and qubit-qubit coupling schemes at the circuit diagram level and realize them in nanofabricated devices 4,6 .
A largely unexplored approach to meet circuit design challenges is to use computational automated discovery. In other fields of science and engineering, automated discovery and inverse design have emerged as a solution to a variety of design problems. All of these problems share the task to identify a set of parameters for which a system of interest yields desired target properties. In contrast to forward design methods, where system properties are estimated from system parameters through direct measurements or computation, inferring the system parameters from target properties is a much more laborious process. Typical automated discovery workflows consist of proposing a set of system parameters and determining the resulting system properties. The similarity between the obtained system properties and the targeted properties is then used as a quantitative measure to refine the system parameters iteratively. Parameter refinements are usually implemented via various optimization procedures or reinforcement learning agents. In the physical sciences context, automated discovery has been applied to nanophotonic on-chip devices 7,8 , complex quantum state generation in optical platforms 9,10 , entanglement creation and removal in superconducting circuits 11 , and further problems in many-body physics 12,13 and chemistry [14][15][16] .
In the field of superconducting circuits, a contemporary design challenge is the implementation of many-body interactions. Such interactions of more than two qubits commonly appear in effective spin models of quantum chemistry, quantum error mitigation schemes, and advanced driver Hamiltonians for quantum annealing 3,[17][18][19][20][21] . They are therefore valuable building blocks for future quantum simulators and quantum annealers. Interactions involving n qubits are called n-local and comprise Pauli operations in the system Hamiltonian that act on n qubits simultaneously. The interaction Hamiltonian has the form where 2M is the coupling strength and α i ∈ {x, y, z} denotes the type of Pauli matrix acting on each qubit. Several implementations for 3-and 4-local couplers in superconducting circuits have been proposed. Hamiltonian gadgets use ancilla qubits to encode the desired interaction in the effective low-energy spectrum of the system 22,23 . Alternatively, individual flux qubits can be coupled to a common coupler circuit that mediates the interaction via a tailored nonlinearity [24][25][26] . The modularity of the tailorednonlinearity approach, in particular, is reminiscent of existing 2local couplers and holds promise for practical implementation 27,28 . Further work is required, however, to find coupler circuits that operate in experimentally accessible parameter regimes, reduce flux noise sensitivity, and eliminate spurious couplings. The automated design task is to find a superconducting circuit diagram that fulfills a set of desired properties. To this end, we introduce a method for superconducting circuit closed-loop automated design (SCILLA) and implement it in software. As shown in Fig. 1 and further detailed in "Methods", it enables parallelized, closed-loop workflows that have access to design algorithms (generators) and can be connected to different circuit evaluators. The method is applicable to a wide range of search problems involving spectral properties and noise sensitivities of galvanically connected superconducting circuits.
In the following, we present the results of two circuit search instances. First, SCILLA is applied to the well-studied example of capacitively shunted flux qubits 29 , testing the key features of the closed-loop implementation. Given our knowledge about the optimal solution, we study the performance of the closed-loop algorithm in reaching the target. We then highlight the competitive advantage of our approach by executing a design workflow that identifies noise-resilient 4-local couplers. The best circuit identified by the software is analyzed to elucidate its underlying operational mechanism, derive general design principles, and show its promise for experimental implementation.

RESULTS
Benchmarking SCILLA by flux qubit design As a benchmark for the automated design software, we define the target to be a capacitively shunted (C-shunt) flux qubit. This is a design variant of the flux qubit that has been shown to yield improved reproducibility and coherence for quantum information processing applications [29][30][31] . It is particularly well-suited for strong qubit-qubit coupling in the quantum simulation and annealing context 28 . The circuit is given by a loop of two identical, large junctions and one smaller junction, with the small junction shunted by a large capacitance. It can be represented by a twonode circuit as shown in Fig. 2a, with the additional bottom node declared ground. The following Josephson junction energies and capacitances are typical for this type of qubit and are chosen as the target parameters: The intrinsic parallel capacitance of each junction is added to the listed capacitances.
The simulated transition energies between the ground state and the first and second excited states as a function of external flux through the circuit loop are shown in Fig. 2b. This energy spectrum is defined as the first target property for the circuit search problem. There is, however, a continuous family of two-node circuits that fulfills this spectrum. The degeneracy is lifted by adding the requirement that the circuit is symmetric, i.e. that the component parameter values E J11 , C 11 are close to E J22 , C 22 , respectively. Rather than enforcing a constraint a priori, circuit symmetry is added as a second target. In this way, the design task tests a case of multi-objective search. In addition, symmetry in the circuit diagram is a practically relevant property that can often be translated to the chip design, limiting noise-inducing effects such as currents in the ground plane. It can be advantageous to have symmetry as a target property, rather than a constraint, such that it is only enforced when the other properties can be fulfilled at the same time. A single scalar loss function L bench is constructed from a combination of the spectrum sub-loss L spec bench and the symmetry sub-loss L sym bench : Here, E 0i and E Ã 0i represent the candidate and target circuits' transition energy from the ground to ith excited state as a function of external flux. The RMS deviation Á j j j j ' 2 between these functions is evaluated over the period of one flux quantum Φ 0 and approximated numerically by computing the function values at discrete flux points. The parameters C max ¼ 100 fF and E max J ¼ h 200 GHz are the upper bounds for the respective parameter values as defined in "Methods". The weighing of the sublosses is chosen such that their values are in the same order for a typical random circuit, and the logarithmic scaling was observed to improve optimizer performance. Therefore, the loss function is constructed by starting from a scalar form that allows for individual sub-loss optimization and then empirically choosing weighting and scaling parameters. We refer to the Methods for more information on loss functions for automated superconducting circuit design.
The workflow specified in the circuit searcher starts with random sampling from the parameter space, followed by gradient-based optimization with the L-BFGS-B algorithm 32,33 . Ten such search jobs are executed independently for greater throughput, each with ten parallelized random samples and subsequent gradient-based refinement. The parameter space for the benchmark is six-dimensional and consists of the junction Fig. 1 Implementation of SCILLA, enabling automated circuit design. a Definition of the circuit design task, for which details about the general circuit architecture, parameter bounds, and design targets are provided. b Based on the general architecture, the design module evokes parameter generating algorithms to place components and choose component parameters in the circuit. c Calculation of circuit properties such as spectra or noise sensitivity. d Estimation of the agreement between computed properties and target properties. e Circuits that are close to the design target are identified. A database system facilitates asynchronous execution and parallelization of the workflow as well as refinement of design choices (closed-loop feedback) based on merit evaluations of previously proposed circuits.
energies and capacitances shown in Fig. 2a. No external flux needs to be specified, because the sole flux degree of freedom is varied to calculate the transition energy spectrum. Moreover, the Hamiltonian of general two-node circuits without linear inductances takes a simple form that does not require the general Hamiltonian construction procedure included in SCILLA. An optimized simulator for such two-node circuits is provided as a separate module and used here. The simplified Hamiltonian is provided in Supplementary Note 3 and the detailed implementation of the benchmark workflow is reported in Supplementary Note 2.
The result of 40 parallel circuit optimizations is shown in Fig. 2c, which represents a subset of four of the ten independent search jobs. A fraction of 10% of the runs (4 of 40) reaches close to the global minimum defined by the target circuit. High accuracy in both the spectral and the symmetry sub-targets is achieved in the yellow-shaded area in the plot, which is reached after 30-60 optimizer iterations for the successful runs. The final spectrum of one successful run (circuit A) is shown in Fig. 2b, matching the target spectrum accurately. Symmetry in circuit A is very high, the deviation between the parameters E J11 (C 11 ) and E J22 (C 22 ) being h × 5 × 10 −9 GHz (2 × 10 −8 fF).
The remaining refinement runs terminate at a higher loss value. Inspection of the final spectrum of one such run (circuit B) in Fig. 2b reveals that the spectrum is not close to the target. We identify the failure mechanism in that the loss function has local minima. In the surface projection of the loss function shown in Fig. 2d, some local minima are visible as bright, yellow clusters around the global minimum. Since the gradient-based optimizer used for refinement is local, it will naturally terminate after reaching one of the local minima.
In summary, nearly all module types and functionalities of SCILLA have been tested in the benchmark. In particular, it is verified that circuit simulation and merit evaluations are executed in parallel and asynchronously as intended. The average run time per search job is 6 h 41 min and evaluates 1.1 × 10 4 circuits in total. Therefore, when including the computational overhead of the closed-loop software, the average simulation time of a circuit is 2.3 s/circuit. This compares well to 1.7 s/circuit in a typical single, isolated evaluation of a circuit on the same hardware. As a more general point about automated circuit design, we observe here that the optimization is nonconvex even for a moderate search problem. A purely gradient-based optimizer is of limited use in this case, although the random starting point ensures that the optimum can be found in some trials. The success rate for finding good circuits is expected to decrease for larger circuits because of the increased parameter space dimensionality. The hardness of the search problem is also increased by more restrictive loss functions with more sublosses. For the 4-local coupler search, the gradient-based optimizer is therefore replaced with an evolutionary strategy, which is able to avoid local minima. In addition, the available computational resources and parallelization capability of SCILLA are used to the maximum extent in order to explore a large portion of the design space and refine as many trial points as possible.

Noise-insensitive 4-local coupler design
We turn to the design challenge of coupler circuits for 4-local interaction of flux qubits. First, we identify the key spectral property of such a coupler that yields the desired interaction when four flux qubits are coupled to its external flux degree of freedom via mutual inductance (see Fig. 3): a double-well profile of the coupler ground-state energy E c 0 versus external flux 24 . Note that this double-well flux spectrum of the coupler should not be confused with the inductive double-well potential of the flux qubit. The ground state 0 j i and excited state 1 j i of each flux qubit are given by persistent left-and right-circulating currents in the qubit loop, respectively. The circulating current change associated with the excitation (relaxation) of any of the flux qubits will shift the flux bias point Φ var of the coupler by a small flux offset + δΦ (−δΦ). It is assumed that all four flux qubits are identical, have no pairwise coupling, and have the same interaction strength with the coupler. If the double-well profile of the coupler spectrum has flux spacings 2δΦ between energy values differing by 2M as indicated in Fig. 3b, then each qubit transition will raise or lower the potential energy of the system by 2M. Applying external flux such that Φ var is biased at the center peak, the system energy separates into two energy manifolds associated with the parity of Fig. 2 Flux qubit benchmark. a Parametrized circuit layout to which the circuit designer is constrained. b Target spectrum generated from a C-shunt flux qubit design, showing the transition energy from the ground state to the first (solid) and second (dashed) excited state. The final spectra of a successful (orange) and a stuck (blue) refinement run are also shown. c Convergence of 40 refinement runs starting at randomly sampled initial circuit parameters. d Visualization of a slice of the loss surface corresponding to the benchmarking problem. The orthogonally projected path of a successful refinement run towards the global minimum is shown.

Fig. 3
Coupler spectral property enabling 4-local interaction. a Four superconducting flux qubits are coupled mutually inductively to a shared coupler circuit. The states 0 j i, 1 j i of each qubit are determined by the direction of the persistent circulating current in the qubit loop. By mutual induction, the magnetic field generated by the qubit's persistent current adds a small flux offset ± δΦ to the coupler circuit. b The coupler has a ground-state energy E c 0 Φ var ð Þ that depends on an external flux Φ var . If this energy spectrum follows a double-well shape with spacing 2δΦ as indicated in the figure, the 4-qubit states with even and odd excitation numbers separate into two energy manifolds and the coupler mediates an effective 4-local interaction term. The challenge is to find a robust coupler circuit with such a double-well spectrum. The behavior of the coupler energy over a small external flux range should not be confused with the inductive double-well potential of the flux qubit.
the qubit excitation number. This is equivalent to the coupler providing a 4-local interaction term The challenge of engineering a 4-local coupler, therefore, becomes that of finding a circuit with the described double-well spectrum in a narrow flux range under realistic parameter constraints. In order to preserve quantum coherence, it is important to limit sensitivity to sources of noise, particularly flux noise if there are additional flux degrees of freedom in the coupler circuit 4,34 . After a circuit fulfilling these properties has been found, simulations of the full system, including all four qubits and the coupler, need to be performed in order to confirm the validity of the 4-local coupling mechanism presented above.
In order to define the loss function, the maximum well-to-well spacing of the coupler spectrum needs to be determined. The bias points for the 4-qubit excitation manifolds, two of which are close to the minima of the double-well profile, are shown as red dots in Fig. 3a. The spacing 2δΦ between the wells and the peak needs to be small such that it can be bridged by the spin flip-and corresponding circulating current change-of a mutually inductively coupled flux qubit. Given the typical persistent currents and mutual inductances in strongly coupled flux qubit systems, we determine that the well-to-well spacing of the spectrum needs to be below 40 mΦ 0 28 . Here, we use that for a qubit-coupler mutual inductance M and flux qubit persistent current I qb p , the equation δΦ ¼ MI qb p can be used to determine the flux offset induced by the qubit in the coupler.
In addition to the spectral property, insensitivity to noise sources in the system needs to be enforced. An often-observed effect of a flux offset in an additional inductive loop of the coupler circuit is shown as the gray trace in Fig. 4a: The double-well spectrum tilts and becomes asymmetric, which would lead to unwanted, non-four-body terms in the interaction Hamiltonian. Therefore, the second objective is to mitigate such spectral shifts due to slow flux noise. Additional sources of noise, such as fast flux and charge noise, can also be quantified by calculating the respective transition dipole moments. As experimental precedent shows that slow flux noise in the circuit loops is of particular importance for the use case of quantum simulation and quantum annealing 28 and as state transitions are implicitly penalized by rewarding large energy gaps, the inclusion of fast noise processes is left for future work. It is also noted that the noise model is intentionally designed to capture noise in the coupler alone. The flux qubits, which are not included in the optimization, would also be subject to noise but that contribution can be treated independently. As the noise sensitivity of flux qubits is widely studied and understood 4 , it is not part of the loss function.
In order to write these targets in a unified loss function, several parameters of the double-well are determined. These parameters are visualized in Fig. 4a (dark red). In the computational routine, they are calculated in case a double-well feature is detected in a 40 mΦ 0 range around the bias point of the primary external flux of the circuit: first, the energy difference h peak between the wells and the center peak determines the 4-local coupling strength and needs to be maximized. Second, the minimal energy difference h split between the ground and first excited state of the coupler is to be maximized; if it is too close to the qubits' transition energy, the coupler can swap excitations with the qubits and break its operating principle of remaining in the ground state at all times. Third, asymmetry induced by flux offsets in the nonprimary loops of the coupler is calculated as the energy difference h sens between the minima of the left and right well. Flux noise is usually dominated by local two-level systems on the metal surface of the circuit, with a degree of flux noise correlation between loops that depends on the length of wire shared by them 4,34,35 . An upper bound on the effect of flux noise is determined by applying flux offsets individually to each nonprimary flux degree of freedom and summing the resulting asymmetries h sens,i . A scalar loss function combines the above parameters. It is given by the p-norm of three terms that quantify the peak height, excited-state splitting, and noise-insensitivity targets.
The hyperparameters h max peak , h max split specify the target values for peak height and excited-state splitting. They serve as cutoff values for the respective parameters of the detected double-well, as shown below. The noise sensitivity is determined as the ratio between the summed asymmetries and the peak height and can assume a maximum value of 1.
h sens ¼ min The hyperparameters of the loss function are chosen empirically Loss of the candidate circuits throughout the optimization, starting with 150,000 randomly sampled circuits and followed by refinement iterations using swarm optimization. Inset: circuit diagram of the best circuit (circuit C) after refinement. c Circuit schematic of the full 4-qubit system, including coupler circuit C. d First 16 energy eigenstates of the full system versus the commonly swept external flux in the qubit loops. Darker lines indicate a higher number of degenerate states. The states separate by qubit excitation parity at degeneracy, which corresponds to 4-local coupling between the qubits.
T. Menke et al.
by the trial of different hyperparameter sets. For the results presented in this work, they assume the values p = 4, h max peak ¼ 1:5 GHz , and h max split ¼ 10 GHz . The loss function is constructed such that joint optimization of two sublosses is favored over individual optimization. The loss function assumes a minimum value of À3 1 p if all sublosses are minimized. A loss of zero is assigned to the candidate circuit if no double-well is detected in the 40 mΦ 0 flux range around the bias point, the peak height is below 50 MHz, or the excited-state splitting is below 100 MHz. In addition, a zero loss is assigned if the simulation fails, times out, or has large Hilbert space truncation errors. Therefore, circuits with little promise for successful optimization are effectively removed from the optimization workflow.
The workflow implemented in the circuit searcher applies insights from the flux qubit benchmark to the more complex task of coupler design with flexible circuit diagrams. It starts with a random sampling of 15,000 circuits, which is about the limit beyond which database operations are observed to slow down. The best two of the 15,000 circuits are kept after filtering and are refined using the evolution-inspired swarm optimization module. Ten such jobs are executed independently for better utilization of the available computing resources. The random search space is chosen to span three-node circuits with connections between the nodes and to the ground on which capacitances, junctions, and inductances can be placed. Under the constraints on the placement of circuit components detailed in "Methods", components are randomly assigned among the available positions. After sampling and filtering, the network configuration is fixed and only the component parameter values are refined. The configuration of the circuit network is therefore flexible, and the number of inductive loops varies between sampled circuits. This constitutes a consequential extension of the fixed-configuration, two-node search space that has been used for the benchmark, allowing for more degrees of freedom in reaching the 4-local coupler target. To calculate the properties of such circuits, the simulation module that implements the general-purpose circuit Hamiltonian simulation is used. The workflow implementation is described in more detail in Supplementary Note 2.
The results of the random search and subsequent swarm optimization are shown in Fig. 4b. A double-well spectrum is detected in 7 of all 150,000 sampled circuits (0.005 %). The swarm optimization improves the merit of all filtered circuits to varying degrees. The final best circuit, named circuit C hereafter, has the following properties after (before) refinement: The double-well spectrum of circuit C emerges from the strong interaction of two rf SQUIDs at different flux bias points. The first SQUID is formed by the loop surrounding the external flux Φ var , which is biased at 0 Φ 0 in the absence of qubits. The second rf SQUID loop encompasses the fixed external flux Φ fix = 0.5 Φ 0 . For details on the derivation of the circuit C Hamiltonian, see Supplementary Note 5. We note that this computer-generated circuit has strong similarities with the 4-local coupler design proposed by Kerman 24 . The key differences are that the parameters are more easily accessible experimentally in our design and the two inductive loops are connected galvanicallyand thus more strongly-than for mutual inductive coupling, allowing for a large double-well peak despite smaller loop currents. Moreover, the asymmetry relative to peak height arising from flux offsets is reduced by a factor of 3.3 in circuit C. While the galvanic connection could in principle have been thought of by a circuit design expert, it also requires other changes to the circuit that we believe are unintuitive.
In total, three of the seven double-well circuits from the automated search results have a similar layout as circuit C, with two inductive loops coupled through an inductance or junction. The remaining four circuits are three-loop circuits and therefore have two nonprimary external flux degrees of freedom. In addition to circuit C, two more exemplary circuit optimization trajectories are highlighted in Fig. 4b: circuits D and E. Circuit D is a three-loop circuit with 1.50 GHz double-well peak and 3.12 GHz excited-state splitting. These desirable properties are, however, contrasted by an increased flux noise sensitivity due to the additional external flux in the third loop. Circuit E has two loops and a similar circuit network as circuit C but features a much larger excited-state splitting of 10.8 GHz and a much smaller double-well peak height of 0.24 GHz. The properties of circuits C, D, and E are listed in detail in Supplementary Note 4. These exemplary search results show that the automated design workflow finds a variety of circuits that fulfill the double-well sublosses to different degrees. It therefore supplies design options with different trade-offs, which can inform both theoretical understanding and practical implementation of the relevant class of circuits.
It remains to be demonstrated that circuit C in fact behaves as a 4-local coupler mediating an interaction H int ¼ M σ z 1 σ z 2 σ z 3 σ z 4 , and thus that our reduction of the design problem to a spectral property is valid. For this reason, a full system simulation of four qubits coupled mutually inductively to the coupler as shown in Fig. 4c is performed (see Supplementary Note 6). The energies of the four-qubit excitation manifolds versus the qubit flux around degeneracy are shown in Fig. 4d, illustrating how the states of different parities separate. A 4-local coupling strength of 2M = 573 GHz in the circulating current basis (σ z basis) of the qubits is extracted, which is lower than the double-well peak height used as a proxy for the coupling. Part of the reduction in coupling is caused by the flux points of the odd qubit excitation manifolds not lining up exactly with the minima of the double-well spectrum, as indicated conceptually by the red dots in Fig. 4a. Further reduction mechanisms could include inductive loading of the coupler by the mutual inductive coupling to the qubits and interactions with the coupler's excited state. Spurious terms of different locality, which are manifest as a splitting of states at the additional spectral crossing points in Fig. 4d, are in the MHzregime and thus negligible. In a physical implementation of the system, additional spurious 2-local couplings would arise from the closeness of the qubit circuits. These can be mitigated by physical separation of the interaction points between each qubit and the coupler and by routing the qubit wires such that mutual inductances are canceled. We conclude that despite additional effects reducing the coupling, SCILLA applied to the double-well loss function was able to design a coupler with several hundred MHz of 4-local coupling strength, even without performing costly full system simulations.
As predicted from the benchmark analysis, our success in finding a promising 4-coupler circuit rests on the careful definition of the loss function, choice of the design workflow, and exploitation of computational resources. Given the low success rate of sampling a circuit with a double-well spectrum at all, a critical step has been to explore a large portion of the design space before attempting to refine promising circuits. Averaged over the ten batches, sampling 15,000 circuits took 21. 51.2 ± 15.9 h (23.0 s/circuit). While the sampling time is relatively uniform, the swarm optimization runtime varies significantly with the circuit layout that is being optimized. In addition, it is observed that the simulation time per circuit is much longer in the swarm optimization, which is expected from the lower degree of parallelization, additional circuit simulations required by the merit evaluation, and a larger number of database operations in the iterative optimization. On top of the runtimes of the sampling and optimization steps, the closed-loop procedure requires an additional 3.30 ± 0.37 h that is mainly spent filtering the sampled circuits before swarm optimization. The filtering time is limited by the input/output operation speed of the used hardware and can be much lower on a different computing cluster.
Overall, SCILLA is able to accommodate the increased computational effort from adding a third node to the circuit and allowing for full flexibility for the circuit network. The added complexity is rewarded by the successful fulfillment of the target properties. This warrants further study of improving the efficiencies of circuit simulation and search algorithms.

DISCUSSION
We demonstrated a computer-driven and highly parallelized approach to the design of superconducting circuits. We evaluated the performance of the developed method SCILLA and discovered circuits for design challenges with multiple objectives. Our results demonstrate that SCILLA is successful in designing circuit architectures with superior performance than an existing proposal, namely a 4-local coupler with several hundred MHz of coupling strength, small unwanted coupling terms, small flux noise sensitivity, and experimentally accessible parameters. Given the modularity of the software implementation, it is straightforward to tackle new circuit design challenges that can be defined in terms of spectral properties and eigenstate expectation values. One such challenge is the design of qubit circuits with a specific sensitivity to control pulses and noise, which can be calculated from transition dipole moments. While the computation of static circuit properties is included in our method, the addition of dynamical time evolution simulations would present a promising future opportunity to include control pulse shaping.
The exponential scaling of the Hilbert space with the number of nodes makes naive scaling of the method by simply adding more nodes to the circuit impractical. However, simulation of larger systems is possible if they consist of small sub-circuits with an intercircuit coupling that is weaker than the intracircuit coupling 36 . In that case, the sub-circuit Hamiltonians may be solved individually with only the lowest eigenstates of the sub-circuits then coupled together, enabling the inverse design of other common architectures such as transmon-based multi-qubit processors. Large circuits can also be made numerically tractable by performing their simulation on a quantum computer 37 . In addition, inspired by Krenn et al. 9 , one can envision functionality to store discovered circuits as building blocks for larger architectures, thus enlarging the set of available circuit components. Future work also entails the translation of the best automatically designed circuit networks into chip designs in order to validate the targeted properties experimentally. The process of developing a chip geometry involves a multitude of layout choices that mitigate undesired effects, such as crosstalk between circuit elements and spurious coupling to the control circuitry. Different layout variants, however, still correspond to the same circuit network that was designed by SCILLA. While automation may also assist in some aspects of chip layout development-and indeed is already part of electronic design automation software for integrated circuits-it is the circuit network that defines most of the relevant circuit properties. We, therefore, believe that our work advances human-computer co-design in the larger cycle of quantum device engineering.

METHODS
The automated design workflow requires a circuit parameter generator, property evaluator, and merit estimator in order to fulfill a specified search task. In the following, it is shown how circuit search can be formulated as an optimization problem by taking the circuit parameters as an input list and calculating the resulting closeness to the target properties. We then describe the closed-loop implementation that handles parameter generation and optimization in a parallelized manner.
Throughout this work, the properties of superconducting quantum circuits are calculated by constructing and diagonalizing the circuit Hamiltonian. For the automated design problem and for spectral engineering problems more generally, the eigenenergies need to be calculated as a function of external degrees of freedom such as flux or charge offsets. Circuit quantization has been discussed extensively in previous work, treating the circuit as a network of lumped elements with canonically conjugate flux and charge variables 4,5,38,39 . We automatically determine the Hamiltonian of a broad class of circuits and solve it efficiently by choosing a mixed representation in the charge and harmonic oscillator bases 36 .

Circuit design as an optimization problem
Automated design of circuits requires a quantitative metric to assess the closeness of a candidate circuit to the target. While the target can take a wide variety of forms, in this work we discuss and illustrate the practicality of optimizing specifically toward spectral properties, symmetry in the circuit network, and flux noise sensitivity. It is straightforward to include more circuit properties as long as they can be obtained from Hamiltonian simulations in reasonable simulation time. In addition, multiple target properties can be combined in a unifying loss function for more balanced design procedures. We define the loss function as L : R d ! R, taking the circuit parameters as an input list x and returning a scalar value. It is defined such that a smaller function value corresponds to a circuit that better fulfills the target and therefore has higher merit. The specific functional form of L depends on the search task and needs to be carefully engineered for the optimization procedure to balance improvement of the sub-targets evenly. It is defined explicitly for each task in "Results". Since the target properties, in general, depend on the results of Hamiltonian circuit simulations, it is assumed that L has access to the simulation results for the input circuit.
The input x contains an ordered list of the capacitance, inductance, and junction energy for each circuit element between each pair of nodes. If there is no circuit element between two nodes, the respective parameter value is zero. Therefore, the ordered list of circuit element parameters fully defines the circuit network and is used to construct the circuit Hamiltonian. The input also contains a list of external flux values-one for each inductive loop in the circuit. Each flux may only take the values 0.0 Φ 0 or 0.5 Φ 0 , ensuring a symmetric spectrum as required for the applications presented in "Results". The circuit parameters are bounded, following typical experimentally accessible values 6,28 : up to 100 fF for capacitances and up to 300 pH for inductances. The lower bound for inductances is chosen as 75 pH in the 4-local coupler design. This ensures that each inductor is large enough for qubits to couple mutually inductively to it. For Josephson junction energies, we allow a range of h × 0-200 GHz for the flux qubit benchmark and h × 99-1982 GHz for the 4-local coupler design. The latter provides more design flexibility and corresponds to junctions with 5 μA/μm 2 critical current density, 0.2 μm junction width, and 0.15-4.00 μm junction length. The intrinsic parallel capacitance of each Josephson junction depends on its respective physical geometry and is implicitly assumed. Throughout this work, we denote the Josephson junction energy, capacitance, or inductance of a circuit element between two nodes i and j by E Jij , C ij , and L ij , respectively.
In addition to the constraints on parameter values, we identify rules for the placement of components in the circuit network. Capacitative elements may be placed in parallel to inductive elements, such as junctions and inductors, i.e. between the same two nodes. However, parallel placement of an inductance and a junction is not allowed to prevent the increased circuit complexity resulting from additional inductive loops formed in this way. The degree of connectivity, i.e. the number of connections that each node has to other nodes, is set by the nonzero parameters in the input x. A circuit with more nodes and higher connectivity allows for more complex spectral engineering in reaching the target. However, it is disadvantageous to increase the circuit size beyond the requirement of the objective. Larger circuits may demand more calibration and more control hardware in experiments. In addition, a larger T. Menke et al.
number of external flux degrees of freedom usually leads to a greater flux noise susceptibility. In the case of two-and three-node circuits-excluding the ground node-as discussed in this work, the circuit network is planar even for all-to-all connection between the nodes. When moving to larger circuits, however, one needs to constrain certain circuit components to zero in order to maintain network planarity and thus ensure realizability in a planar on-chip architecture. In principle, circuit complexity and experimental feasibility can also be added to the target with a suitable loss function.

Realizing closed-loop circuit design
SCILLA is a closed-loop implementation for accelerated computational circuit design available on GitHub (see "Code availability"). The procedure autonomously searches the circuit space for circuit architectures satisfying desired target properties. Similar approaches have already been successfully applied in the context of autonomous discovery and experimentation in chemistry and materials science 40,41 . SCILLA contains a general-purpose method to compute properties of superconducting circuits (not provided on GitHub). It is thus applicable to a wide range of circuit design applications.
The typical workflow in the closed-loop implementation of SCILLA is illustrated in Fig. 1. Each circuit search is started from a set of general instructions, which inform SCILLA about the size of the circuit and the number of components (see Fig. 1a). With these specifications, a design algorithm chooses the kinds of components to be placed at particular locations in the circuit graph, as well as component-specific parameters (see Fig. 1b). Then, properties of interest are computed for each designed circuit (see Fig. 1c). As the calculation of circuit properties is typically the most time-consuming step of this workflow, SCILLA computes only the properties which are requested in the general instructions. SCILLA supports the computation of circuit spectra as well as estimations of the flux and charge noise sensitivities. The computed circuit properties are then used to evaluate the loss L, determining how well a given circuit matches the desired target specifications (see Fig. 1d). Finally, the circuit composition which best satisfies the desired targets is determined (see Fig. 1e). The crucial element to "close the loop"-and thus enable autonomous circuit design-is to report the loss of a proposed circuit architecture to the design algorithm. Based on the feedback, the design algorithm can propose refined circuit architectures with improved target properties. Some design tasks require higher-level operations, such as computing a gradient during circuit design or probing the flux sensitivity of a considered circuit only if a predefined condition is met. Therefore, both the circuit design algorithms, as well as the circuit merit evaluators, are equipped with the capability to define new circuit evaluation tasks. Circuit evaluations requested during circuit design are prioritized to assure that individual circuit design tasks are completed quickly.
The workflow presented in Fig. 1 can be parallelized with little computational overhead by separating each step of the workflow into independent units (modules). Information about individual circuits, such as component parameters or circuit properties, are dynamically stored in a system of databases built on the SQLite database engine. The database implementation is the key component that allows the software to decouple each of the steps in the workflow and execute them individually. Details on the implementation are provided in Supplementary Information (see Supplementary Note 1). The implicit parallelization of individual modules maximizes the number of circuits evaluated in parallel and thus utilizes available computing resources to a higher capacity. Moreover, circuit properties can be computed asynchronously, which is crucial as the computational time required to evaluate a given circuit often depends on its component configuration and parameter values. With the implemented database system, all circuit parameters and properties which have been proposed and evaluated during the closed-loop sampling procedure are easily accessible. As data are stored in a standardized format, profound post-process analysis of the circuits is possible and has the potential to identify general principles for circuit design, which is demonstrated in "Results".
The design modules within SCILLA can efficiently search the space of small circuits with relatively few components that can be evaluated quickly, as well as larger circuits with more components for which property calculations are more time-consuming. As the hardness of the optimization depends on the search space size and loss function definition, a number of different design strategies are implemented in the closed-loop, ranging from random search strategies 42,43 for coarse, massively parallelizable sampling over gradient-based methods 32,33 to evolutionary strategies [44][45][46] .
A more detailed description of the supported methods is provided in Supplementary Note 1.
The closed-loop framework also enables the implementation of multistep workflows. The user declares a search procedure based on a chosen design algorithm, can define analyses on the circuits sampled during this first search, and then trigger another circuit design round, for instance, to refine the previous search. Arbitrary combinations of different design and analysis steps are possible and can be easily integrated into the desired multistep workflow. Examples of such workflows are provided in the circuit design applications presented in "Results", as well as in Supplementary Note 2.

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