Efficient modeling of superconducting quantum circuits with tensor networks

We introduce an efficient tensor network toolbox to compute the low-energy excitations of large-scale superconducting quantum circuits up to a desired accuracy. We benchmark this algorithm on the fluxonium qubit, a superconducting quantum circuit based on a Josephson junction array with over a hundred junctions. As an example of the possibilities offered by this numerical tool, we compute the pure-dephasing coherence time of the fluxonium qubit due to charge noise and coherent quantum phase slips, taking into account the array degrees of freedom corresponding to a Hilbert space as large as$~15^{180}$. Our algorithm is applicable to the wide variety of circuit-QED systems and may be a useful tool for scaling up superconducting-qubit technologies.


I. INTRODUCTION
Superconducting qubits are a leading platform for quantum information processing [1,2]. These qubits are built from superconducting quantum circuits integrating linear elements, such as capacitors and inductors, together with the only known nonlinear and nondissipative circuit component: the Josephson junction. These circuits operate at milliKelvin temperatures where macroscopic electromagnetic degrees of freedom associated to currents and voltages in the circuit are described quantum mechanically [3,4]. In this regime, nodes (or branches) of the circuit are represented by bosonic fields with, in principle, infinite Hilbert-space dimension. The circuit topology defines linear and nonlinear interactions between these bosonic modes. Determining the low-lying excitations of the circuit in the presence of such interactions requires the diagonalization of the full circuit Hamiltonian. However, for circuits with more than a few nodes, this rapidly becomes intractable by exact diagonalization. With current devices integrating 10s [5] to 100s [6], 1,000s [7] and even 10,000s [8] Josephson junctions, finding new methods to efficiently model these devices is one of the challenges that the field is facing.
Most superconducting quantum devices operate in regimes where effective models with a reduced number of degrees of freedom are accurate enough to describe the physics of interest. However, these effective models are based on approximations that allow extracting only limited information about the system. Moreover, it is often not possible to trace back the original circuit parameters from the effective model and, when it is possible, these parameters have to be inferred indirectly from complex multivariate fits to the experimental data. This loss of information can be detrimental to circuit design.
In this work, we adapt to many-body superconducting quantum circuits a numerical tensor network method that we have introduced in Ref. [9]. We use this numerical toolbox to compute the relevant low-energy excitations of a large-scale superconducting circuit taking into consideration all of the degrees of freedom of a lumpedelement model of the device. We show how this gives access to information about the system that can be used, for instance, to estimate the device coherence times from first principles.
As an example of application of this method, we consider the fluxonium qubit [5]. This superconducting quantum circuit is made of a small Josephson junction shunted by an array of ∼ 100 Josephson junctions. Because of the large number of elements in the fluxonium circuit, this qubit is an ideal testbed for our numerical approach. Moreover, solving the complete fluxonium circuit Hamiltonian is challenging due to the short-and long-range linear and nonlinear interactions of the model, which is formulated under periodic boundary conditions. To benchmark our tensor network implementation, we develop an effective model for the fluxonium qubit that captures the essential circuit details and which can easily be solved by exact diagonalization. To assert the validity of the tensor network method, we first compare results obtained with this technique to those obtained with the approximate effective model in regimes where the latter approach is expected to faithfully describe the device. We then push the tensor network method to regimes where deriving an accurate effective theory is difficult. The effective model and the tensor network toolbox are used to investigate the charge dispersion of the fluxonium qubit in a broad range of parameters, confirming an existing theory [10] and clarifying its regime of validity. Finally, we use the tensor network method to estimate the puredephasing time of a realistic fluxonium device. We provide direct numerical evidence of the potentially harmful effects of charge noise in this system for certain circuit parameters.
This paper is organized as follows. In Sect. II, we summarize the tensor network method introduced in Ref. [9]. In Sect. III, we provide a tensor network implementation of the complete fluxonium-qubit Hamiltonian, describe an effective model for this qubit and compare results obtained with both approaches. Sect. IV discusses the interplay between charge noise and coherent quantum phase slips in the fluxonium qubit. The main result of this section is the direct numerical evidence of the charge dispersion in fluxonium devices, supporting a previously developed theory [10]. Sect. V is dedicated to the conclusions and to an outlook of the results of this work.

II. THE MULTI-TARGETED DMRG ALGORITHM
A useful strategy to determine the low-energy excitations of a quantum system is based on decomposing the many-body wavefunction into a series of tensors, each representing a single site (or mode). The form of the resulting wavefunction is called matrix product state (MPS) and has been known for some time [11]. For a review, see for instance Refs. [12][13][14][15][16]. The tensor decomposition applies to the full many-body wavefunction where σ i indexes orbitals (or levels) that belong to a finite-dimensional basis of states for the ith site. For a site representing a bosonic mode, a finite-dimensional basis for this site may be defined by truncating the site's Hilbert space. The probability amplitude c σ1σ2...σ N J in Eq. (1) is interpreted as a tensor with N J indices, N J being the number of sites. In order to obtain a MPS representation of |ψ , a series of tensor decompositions can be performed using the singular value decomposition (SVD). The SVD decomposes a tensor into two isometries, U and V , and a diagonal matrix D such that the original tensor may be reconstructed as U DV † . By performing successive SVDs on the full original tensor, one obtains a site-by-site representation of the wavefunction of the form [13] where A σi ai−1ai is the tensor of the MPS associated to the ith site. Here, an extra index a i appears corresponding to a link index that connects to an adjacent site. The dimension of this additional index is known as the bond dimension and is controlled by truncating the number of nonzero singular values that are kept in the diagonal matrix D of the SVDs. Effectively, this truncation leads to a compressed representation of the many-body state, leaving out small entries of the density matrix which are unimportant to understand the physical phenomenon of interest. Physical systems that can be modeled efficiently by a MPS with a much smaller bond dimension than the full wavefunction often involve short-range interactions and low dimensions [17]. Other cases can also be captured by a MPS at the price of using a larger bond dimension [12,13].
Equation (2) is represented in the left-normalized basis where the tensor A is determined from the U tensor of the SVD. The MPS can also be written with rightnormalized tensors (creating tensors from V † ). The most common gauge to choose is the mixed-canonical representation [13]. There, left-and right-normalized tensors are separated by one site where the D matrix has been contracted on to the site. This site is known as the orthogonality center, and represents the information passed between the left and right parts of the system.
In practice, the MPS is obtained by first constructing the Hamiltonian as a tensor network, known as a matrix product operator (MPO). Once the MPO is specified, an algorithm can be designed to converge from a starting initial state to the correct ground state. A well-known tensor network method to achieve this is the density matrix renormalization group (DMRG) algorithm [18,19]. This approach is found to be efficient for solving systems that are well captured by the MPS and can converge to the ground state in only a few iterations of the algorithm [17,20,21]. More importantly, the complexity of this algorithm scales linearly with the number of sites, making it possible to treat systems of sizes well beyond what is possible with exact diagonalization.
While DMRG is most commonly used to study groundstates, the analysis of superconducting quantum circuits requires us to determine several low-energy excitations. For example, in the case of a single superconducting qubit built using some large superconducting circuit, the ground state and the two first lowest energy excitations are needed to estimate the qubit frequency ω 01 and anharmonicity ω 12 − ω 01 , where ω i is the energy of the ith eigenstate of the circuit and ω ij = ω j − ω i . If n q such qubits are integrated on a chip, the number of excitations required to characterize the device typically scales as n 2 q . The conventional approach to compute excitations with DMRG is to add to the system Hamiltonian an energy penalty of the form i∈ex. Λ|ψ i ψ i |, with Λ > 0, where ex. denotes a set of previously determined excitations {|ψ i }. This energy penalty forces the previously determined low-energy excitations above the next excited state, which becomes the ground state of the modified Hamiltonian and for which standard DMRG can be run [13]. However, it can be noticed that this technique can miss excited states and suffers from convergence issues.
To remedy this problem, we have derived an extension of the DMRG algorithm that includes the excitations computed directly in the Lanczos step of the algorithm [9]. We extended the original MPS to a bundled MPS, where the orthogonality center has been given an additional index that identifies excitations in the system. By attaching this additional index to the state, we can derive an efficient tensor network update at each step of the DMRG algorithm that modifies the wavefunction of each excitation until the energy is variationally minimized to the correct eigenvalue. This procedure, that we name the 'multi-targeted' DMRG algorithm, is numerically stable and does not miss excitations or introduce numerical degeneracies in all tested situations.
Indeed, we have used this method to obtain tens or hundreds of excitations simultaneously, all in a single run Detailed circuit scheme including a "black-sheep" junction (center) shunted by a capacitance (top) and a junction-array superinductance with NJ junctions (bottom). Stray capacitances to ground are depicted in a lighter shade of blue. (b) Effective circuit in which the junction-array is modeled as a linear inductance. φi for i ∈ [0, NJ ] denotes the superconducting phase at every circuit node, while θi for i ∈ [1, NJ ] is the phase difference at every junction of the array. The superinductance (or fluxonium) mode is defined as the phase difference across the black-sheep junction: of the multi-targeted DMRG algorithm. This is where our newly developed technique differs significantly from the traditional DMRG approach for computing excitations, which needs to be run sequentially, once per required excitation. Furthermore, an important benefit of our multi-targeted DMRG algorithm is that the orthogonality of the computed excited states is guaranteed up to numerical precision. In contrast, in the traditional DMRG approach, the degree to which orthogonality conditions are satisfied within a set of computed eigenstates is determined by the accuracy of the associated eigenvalues. More information on the multi-targeted DMRG algorithm can be found in Ref. [9].

III. DMRG IMPLEMENTATION OF THE FLUXONIUM-QUBIT HAMILTONIAN
We choose the fluxonium qubit [5] as a testbed for the multi-targeted DMRG approach. Because of its relatively complex structure, with a Hamiltonian that includes periodic boundary conditions as well as short-and longrange linear and nonlinear interactions (see appendix A), this is an ideal test circuit for this numerical method. We note that non-multi-targeted DMRG has previously been used to study quantum phase transitions in Josephsonjunction rings [22,23] and the coherence properties of the current-mirror qubit [24].
The fluxonium qubit is a variation on the transmon qubit [25] in which a large shunt inductor is added to protect the device against low frequency charge noise [26]. Recent experiments have demonstrated long coherence times with this qubit [6,27,28]. The fluxonium circuit (see Fig. 1) consists of a small Josephson junction, referred to as the "black-sheep" junction, shunted by a superinductance, i.e. a circuit element with ef-fective impedance greater than the quantum of resistance R Q = h/(2e) 2 6.5 kΩ and self-resonance frequencies above 10 GHz [29][30][31][32][33]. Superinductances have been made using Josephson junction arrays [5,32], highkinetic-inductance superconductors [34,35] and granular aluminium [36,37]. Superinductances are also crucial to other qubit designs such as the noise-protected 0 − π qubit [30,38]. While a superinductance is in principle a multimode device, it can behave as a single-mode linear inductance under appropriate conditions [32,39,40]. The multimode structure of such a device has, however, important consequences [32,35], some of which are investigated below.
A. Setting-up the multi-targeted DMRG algorithm With the objective of determining the low-energy excitations of the full fluxonium device shown in Fig. 1 (a) using our multi-targeted DMRG algorithm, we first describe the associated circuit Hamiltonian. In this circuit, the black-sheep junction is described by its Josephson energy E J b and its capacitance C J b which may include a shunt capacitance. We take the superinductance to be realized by an array of Josephson junctions, with L Ji and C Ji being the ith junction inductance and capacitance, respectively. Moreover, a ground capacitance C 0i is associated to the ith circuit node. In the absence of circuit element disorder, these parameters take the constant values L J , C J and C 0 , respectively. We also define the junction plasma frequency ω p = 1/ √ L J C J and reduced impedance z = L J /C J /R Q . Following the standard circuit-quantization procedure [3], the Hamiltonian of the circuit of Fig. 1 takes the form (see appendix A) (3) In this expression, H 0i = 4E Ci (n i − n gi ) 2 − E Ji cos θ i is a noninteracting (or site) Hamiltonian for the ith array junction, where θ i is the phase difference across that junction and n i the conjugate charge. Moreover, n gi is an offset-charge parameter, E Ci is the effective charging energy of this junction and E Ji = ϕ 2 0 /L Ji is the Josephson energy with ϕ 0 = Φ 0 /2π where Φ 0 = h/2e is the flux quantum. In addition to the on-site energies, Eq. (3) includes a bilinear interaction ∝ n i n j arising from the ground, black-sheep and array-junction capacitances, that couples the sites with comparable strength and allto-all connectivity (see appendix A). Furthermore, the last term of Eq. (3) is a nonlocal interaction that depends on the external flux Φ ext = ϕ 0 ϕ ext and which results from the strongly nonlinear Josephson potential of the black-sheep junction. Because Eq. (3) includes a very large number of degrees of freedom and is therefore difficult to work with, this Hamiltonian is typically not directly employed in the literature to describe the fluxonium qubit. Instead, fluxonium devices are usually modeled by a phenomenological Hamiltonian that incorporates a single bosonic degree of freedom, φ = N J i=1 θ i , known as superinductance or fluxonium mode [5].
To obtain the low-energy excitations of Eq. (3) by means of a tensor network method, and in this way go beyond the usual effective model, the circuit Hamiltonian must first be converted to its matrix product operator form. Crucially, we noticed that the longrange cosine interaction is ideally suited to matrix product states and operators, preventing an increase of the bond dimension with the number of sites. This observation is one of the key findings of our work and extends to all circuit-QED Hamiltonians, from lumped-element models to black-box-quantization [41,42] and energyparticipation-ratio [43] formalisms. Indeed, we have successfully implemented a wide variety of such models and circuit Hamiltonians, results that will be reported elsewhere. On the other hand, the all-to-all capacitive interaction in Eq. (3) does not have an efficient MPO representation. However, this unfavorable interaction does not prevent an efficient implementation of the multitargeted DMRG algorithm, as the results that are presented below are obtained with a relatively small bond dimension using MPO compression techniques [44]. The efficient matrix-product-operator representation of the black-sheep Josephson potential in Eq. (3), and the possibility of handling an arbitrary capacitive coupling Hamiltonian by compression techniques, makes our DMRG implementation readily applicable to the wide variety of circuit-QED setups.

B. Effective single-mode theory
To assert the validity of our DMRG method, we derive in appendix B an effective single-mode theory from Eq. (3) that can be solved by exact diagonalization, and which goes beyond the standard treatment found in the literature. Under approximations controlled by the parameter regime of the device, we arrive at the Hamiltonian where the mode described by φ is closely related to the superinductance (or fluxonium) mode φ, and n the conjugate charge.
Here, E C , E L and E J are, respectively, effective capacitive, inductive and Josephson energies obtained from the classical normal-mode structure of the circuit. If the ground capacitances C 0i for i ∈ [1, N J ] can be neglected, then φ = φ and n = n = N −1 J N J i=1 n i , where n is the conjugate charge operator to φ. Otherwise, the φ mode includes corrections to φ that are linear in C 0 .
Although in the limit of large N J Eq. (4) reduces to the usual effective model for the fluxoniumqubit [see Fig. 1 (b)] [5], the parameters of Eq. (4) capture the full circuit's capacitance network and contain important corrections due to the nonlinearity of the array junctions. These corrections can lead to significant frequency shifts of the qubit transitions (see appendix B 3). Crucially, because of its single-mode nature, Eq. (4) can easily be diagonalized numerically by truncating the Hilbert space of the φ mode to finite dimension.

C. Comparison
Having derived the effective model of Eq. (4) which will be used as a benchmark, we are now in a position to demonstrate the results of our DMRG approach and to explore the capabilities of this method. To this end, we consider a device in the 'heavy fluxonium' regime [6,27,35] with a large shunt capacitance and a superinductance made of N J = 120 identical junctions where ω p /2π = 25 GHz and z = 0.03 [32]. See appendix B 2 for a qualitative description of the different regimes of the fluxonium qubit Hamiltonian. Each junction is modeled as a multilevel system using the 15 lowest energy eigenstates of the site Hamiltonian H 0i . We find that for low-impedance junctions, the site eigenbasis requires a smaller number of states to avoid truncation errors as compared to other local bases such as the charge basis. The DMRG implementation is thus defined in a product basis of local wavefunctions spanning a manybody Hilbert space as large as 15 120 and that has, a priori, no built-in information about collective modes of the system. Importantly, this choice of basis also makes our treatment readily extensible to other superconducting quantum circuits.  4), black dashed lines] as a function of the external flux Φ ext . We find excellent agreement between these two independent models. Importantly, this observation extends to all systems sizes and parameter sets that we have tested, from a few-sites fluxonium-like device to circuits with more than 200 junctions. These results provide supporting evidence of a successful DMRG implementation of the fluxonium qubit Hamiltonian. Moreover, this motivates applying the DMRG technique in regimes of parameters where deriving an effective model is not possible. Further numerical evidence is presented in appendix B.

D. Exploring the DMRG results
In addition to computing global properties of the circuit, such as its energy spectrum, the multi-targeted DMRG algorithm also gives access to local site properties and n-body correlators. These operators can give insights into the many-body structure of the fluxonium eigenstates. The purpose of this section is to motivate the use of our DMRG algorithm to explore some of these quantities.
As an example application, Fig. 2 (b) shows the mean photon-number population p i = ψ k |H 0i |ψ k / ω p of the ith site, for all sites (i ∈ [1, 120], vertical axis of each of the 6 density plots) as a function of Φ ext . These expectation values are computed for a given eigenstate |ψ k of the full fluxonium circuit, from the ground state (k = 0, bottom density plot) to the 5th excited state (k = 5, top density plot). Because of the absence of circuit-element disorder in these simulations, the results do not show any variations with site number. We observe that the photon-number population of the array junctions is relatively low for the ground state. The same is true for some excited states whose energies change rapidly with the external flux (fluxons). Note that energies are given with respect to the ground state energy, which is chosen to be always 0. In other words, the energy of the ith excited state as illustrated in Fig. 2 (a) corresponds to that of the transition |ψ 0 → |ψ i . Moreover, we note that the photon-number population of the array junctions is relatively high for excited states that have a weak frequency dispersion as a function of Φ ext (plasmons). We interpret these results with the help of Fig. 2 (c), which illustrates a portion of the local Josephson potential of an array junction and its single-site wavefunctions. From the point of view of this site (left panel), a fluxon state |ψ k involves a small displacement by α k /N J of the site's wavefunction (red) away from its noninteracting ground state position (light blue), where α k is a real number. With the current operator associated to the ith junction given by I i = I c sin θ i where I c is critical current, this displacement results in a circulating current for α k = 0. In addition to this mean-field displacement, plasmon states involve non-negligible population of the sites' excited states, as shown in Fig. 2 The above interpretation becomes clearer by considering the effective potential and wavefunctions obtained from the single-mode effective Hamiltonian Eq. (4), as shown in Fig. 2 The shape of the effective potential is determined by the cosine potential of the black-sheep junction and the inductive en- While fluxon states correspond in this picture to the lowest energy eigenstates associated to the local minima of the effective potential, plasmon states correspond to intrawell excitations (see also appendix B 2). The effective model potential connects to that of Fig. 2 for an excitation |ψ k localized in a single potential well. Thus, in this case, the displacement of the sites' wavefunctions adds to a collective value α k that approximately coincides with the position of a local minimum of the effective potential. This is examined further in Fig. 2 (e), which shows the expectation value of the phase drop φ 0 − φ i ≡ i j=1 θ j , obtained from DMRG and plotted as a function of the site number for the fluxon states |ψ 0 and |ψ 2 at Φ ext = Φ 0 /4 in Fig. 2 (d) [middle panel]. In this figure, the expectation value ψ k |(φ 0 − φ i )|ψ k is represented by the angle between the direction of a vector localized on the ith site with respect to the vertical direction. Thus, the total angle between the vectors belonging to the first and last sites can be identified with the positions of the local minima α 0 and α 2 of the effective potential of Eq. (4).
Overall, Fig. 2 shows that the multi-targeted DMRG algorithm correctly reproduces the results of the effective single-mode theory. It can also provide information that is not accessible from this theory. This comparison provides solid evidence of a correct DMRG implementation of the full circuit Hamiltonian of the fluxonium qubit. It also suggests that other circuit Hamiltonians can benefit from this numerical method. Moreover, the local physical quantities such as those illustrated in Fig. 2 (b), contain information about the energy-participation ratio of all circuit components for a given collective excitation. This information could be used to identify limiting dissipation channels and to understand the effect of circuit-element disorder. We return to these aspects in Sect. V.

IV. CHARGE DISPERSION AND COHERENCE TIME
We now proceed with a concrete application that shows how our DMRG implementation can be leveraged to pro-duce coherence-time estimates from first principles. In particular, we are interested in quantifying the coherence time of the fluxonium due to the combined effect of charge noise and coherent quantum phase slips [10,33].

A. Charge dispersion
In the fluxonium qubit, the black-sheep junction acts as a weak link that couples flux states of the superconducting loop. This mechanism makes quantum control of the flux degree of freedom possible but can also be a source of errors. In a semiclassical picture, the rate at which a quantum of flux can tunnel in and out of the loop through the black-sheep junction is proportional to the junction impedance, while the energy cost associated to the addition of a quantum of flux to the loop scales as 1/L. Since the tunneling of a flux quantum corresponds to a change of 2π in the phase of the superconducting order parameter, this phenomenon is known as coherent quantum phase slip (CQPS) [10,[45][46][47][48][49][50]. In experiments, fluxonium devices exploit a wide range of black-sheep junction impedances, ranging from relatively small in the heavy-fluxonium [6,27,35], to moderate in the fluxonium [5,10] and to large values for the lightfluxonium [51]. See appendix B 2 for a qualitative discussion of these parameter regimes. Ideally, the total amplitude for CQPS events is largely dominated by the contribution from the black-sheep junction. However, if the impedance of the array junctions is large enough, the added CQPS amplitude due to the superinductance can be non-negligible. In this limit, the junction array may be regarded as a "slippery" superinductance [33].
Reference [10] introduced an effective model describing the effect of CQPS events occurring in the superinductance of a fluxonium qubit. In this model, CQPS events due to the black-sheep junction are captured by a phenomenological single-mode fluxonium qubit Hamiltonian similar in spirit to Eq. (4). On the other hand, CQPS due to the superinductance enter in the effective Hamiltonian via the external flux. More precisely, the parameter Φ ext in Eq. (4) is replaced by Φ ext + m Φ 0 , where m is an integer-valued number operator that counts the number of CQPS in the superinductance. Since a CQPS event at any junction of the superinductance leads to a jump m → m ± 1, it can be interpreted as a 2π phase bias on φ.
To quantify the total CQPS amplitude resulting from the superinductance, we consider a realistic model of this composite circuit element with its N J islands and their independent offset charges [see Fig. 1 (a)]. As a consequence of the Aharonov-Casher effect, the fluxtunneling amplitude at a given array junction has a welldefined phase given by the offset-charge n gi associated to that junction [10,45,50,[52][53][54]. By adding coherently the contributions from the N J array junctions, the total CQPS amplitude (excluding the black-sheep junc- determines the charge dispersion of the ground state energy of the transmon Hamiltonian H 0i in terms of the reduced impedance z i and plasma frequency ω pi of the ith array junction [25,33,45,55]. Importantly, this result only holds in the low-impedance limit (z i 1). CQPS events in the superinductance can then be described by a phenomenological flux-tunneling Hamiltonian of the form removes (adds) a single flux quantum from the loop through any of the array junctions. In the limit of rare CQPS, |E S | E L , H CQPS can be regarded as a small perturbation to the fluxonium Hamiltonian. In this situation, first-order perturbation theory predicts a shift is a 2π-displacement operator whose expectation values are computed using the unperturbed eigenstates {|ψ i } with m = 0 [10]. For a homogeneous array , and the total charge dispersion of the qubit transition frequency is As the classical flux states of the loop are degenerate at Φ ext = Φ 0 /2, the effect of a nonzero E S is stronger close to this flux bias. Figure 3 shows the charge dispersion of the fluxon transition of a fluxonium device with parameter values chosen to be as close as possible to those of the experiment of Ref. [10]. The top panel shows the qubit transition frequency as a function of the external flux close to Φ ext = Φ 0 /2 for different values of the offset charge n gi ≡ n g , assumed to be the same on every junction of the array. Each sub-panel shows the DMRG results for a given value of the array-junction impedance. The lightest (darkest) transition in purple corresponds to n g = 0 (n g = 0.5). Since n g = 0.5 is a charge degeneracy point of the single-array-junction Hamiltonian H 0i , Cooper-pair transport between the circuit islands is relatively easier, leading to a stronger flux dispersion in comparison to the case of n g = 0 [56]. Dashed black lines show the qubit transition according to Eq. (4) which does not have an offset-charge parameter. Note that the offsetcharge dependence of the CQPS tunneling energy leads to constructive (|E S | > 0) and destructive (E S → 0) interference of CQPS events.
Qualitatively, charge dispersion increases rapidly with z due to the exponential scaling of Eq. (5). This is best illustrated by the bottom panel of Fig. 3, which shows the charge dispersion for Φ ext = Φ 0 /2 as a function of z. Light-blue circles (Full DMRG) correspond to a fully numerical estimation using DMRG for which the charge dispersion is computed by taking the difference between the energy of the fluxon transition for n g = 0  6), with both approaches showing a small but clearly visible deviations from the results obtain from fully numerical DMRG estimation (Full DMRG) at large z. Indeed, we observe a remarkable agreement between the estimation of the total charge dispersion from fully numerical DMRG and that predicted by Eq. (6), up to array-junction impedances as high as z 0.1. This provides evidence in support of the theoretical model introduced in Ref. [33]. Although not visible in Fig. 3, small deviations between the fully numerical DMRG estimation and those based on Eq. (6) are present for z 0.06. The largest truncation error for all simulations in Fig. 3 is of order 10 −11 , and the error tolerance on the eigenvalues are set to 10 −12 , guaranteeing the convergence of the fully numerical DMRG results to the same accuracy. DMRG being a variational method, we have verified that the convergence to the reported accuracy is also well be-haved. We noticed deviations of the same order of magnitude between the fully numerical DMRG estimation and the prediction of Eq. (6) for devices with different sets of circuit parameters.
On the other hand, the large relative difference between the full numerical multi-targeted DMRG estimation and those based on Eq. (6) in the range of z 0.1 is expected. Indeed, in this regime, Eq. (5) and the assumption that |E S | E L are both no longer valid [33]. Therefore, z 0.1 is a regime in which the DMRG method is at a clear advantage over effective theories.

B. Coherence-time estimations
Because of unavoidable charge noise, the value of δω ij fluctuates in time, leading to broadening of the qubit transition. For large charge dispersion, this effect can severely compromise qubit coherence. This observation is the basis of the experimental study of Ref. [10], where the reduction of the qubit coherence time around the flux sweet spot is taken as indirect evidence of CQPS events in the "slippery" superinductance. In support of the experimental observation and as a further example of the power of the multi-targeted DMRG algorithm, we show below that full DMRG simulation of a device with similar circuit parameters to those reported in Ref. [10] predicts the pure-dephasing coherence times to be dominated by the combined effect of charge noise and CQPS around Φ ext = Φ 0 /2. Moreover, the coherence-time values that we obtain with this method result very close to those measured experimentally.
In order to estimate the coherence times, we follow closely Ref. [10] assuming that the variables n gi are independent and randomly distributed. The probability density function of Re[E S ] can then be approximated by a Gaussian distribution with zero mean and standard deviation N J /2 0 [10]. Following this expression, the effective broadening of the qubit transition scales as √ N J , something which translates to the puredephasing rate 1/T ϕ,CQPS = |∆ω 01 |/4 √ N J [10,33]. To identify the domimant dephasing mechanism, we compare this timescale to that expected for 1/f flux noise by deriving in appendix C a multilevel pure-dephasing master equation of the form where Γ kl ϕ are time-dependent pure-dephasing rates proportional to the 1/f flux noise amplitude, σ kl = |ψ k ψ l |, and D[x, y] ρ = xρy † − {y † x, ρ}/2 is a generalized dissipator operator. By integrating Eq. (7), we define the flux-noise coherence time T ϕ,Flux by the implicit equation ρ 01 (T ϕ,Flux )/ρ 01 (0) = 1/e that we solve numerically. and C0 = 0, extracted from Ref. [10]. The 1/f flux-noise amplitude is taken to be AΦ = 1.2 µΦ0, which is a conservative value [25].
Figure 4 (a) shows the energy spectrum of the simulated device as a function of the external flux, results that should be compared to those of Ref. [10]. In contrast to the results in Fig. 2 (a), the difference between the DMRG and single-mode simulations for the parameters of Ref. [10] is sightly more noticeable due to the low plasma frequency of the array junctions ω p /2π = 12.5 GHz, around which ∼ 40 other additional circuit modes lie [35]. This makes any single-mode approximation invalid, except at low frequencies. Furthermore, Fig. 4 (b) shows the estimation of the device's coherence times using only the results from DMRG as a function of the external flux and close to the bias point Φ ext = Φ 0 /2. We find values which are very similar to the experimental observation (see Fig. 4 in Ref. [10]), thus providing further numerical evidence of the combined effects of charge noise and CQPS. This mechanism dominates over flux noise close to the device's flux sweet spot, resulting in sub-µs coherence times for the device parameters of Ref. [10], in agreement with the experimental observations.
Combined, the results of Fig. 3 and Fig. 4 illustrate the rich interplay between charge noise and CQPS in the fluxonium architecture. Added to the improved simulation capabilities provided by the multi-targeted DMRG algorithm, these findings motivate a systematic experimental study to understand these effects further.

V. CONCLUSIONS AND OUTLOOK
We have developed a multi-targeted DMRG algorithm to simulate large-scale superconducting quantum devices. As an example, we have applied this numerical technique to the fluxonium qubit. The fluxonium circuit integrates a large number of degrees of freedom with linear and nonlinear short-and long-range interactions that are subject to periodic boundary conditions. Combined, these features make this model a challenging target for our DMRG algorithm. To assert the validity of the DMRG simulations, we have developed a detailed single-mode theory for the fluxonium qubit. Finally, we have employed DMRG to investigate the combined effect of charge noise and coherent quantum phase slips in the fluxonium qubit, confirming the theoretical model introduced in Ref. [10] and reproducing some of the experimental findings of that work.
Having access to the expectation values of local and of n-body operators makes it possible to investigate the many-body properties of superconducting quantum circuits. This could help, for instance, in finding new approaches to encode quantum information nonlocally in protected subspaces by exploiting entanglement in these systems. Moreover, local information of large-scale superconducting quantum circuits may be used to evaluate the impact of dissipation channels and circuit-element disorder. This might also lead to a more detailed understanding of dissipation and decoherence mechanisms. Our numerical approach also has the potential to enable advancements in several areas of superconducting-qubit research. In particular, we envision future applications to the analysis of multi-qubit devices and the design of scalable superconducting-qubit architectures.

Hamiltonian without gate voltages
We derive the circuit Hamiltonian used in the DMRG calculations presented in the main text. We consider a fluxonium device where a black-sheep Josephson junction with capacitance C J b (including both shunt and junction capacitances) and Josephson energy E J b is shunted by a superinductance made of N J junctions, each of capacitance C Ji and energy E Ji with i ∈ [1, N J ]. We moreover assume that each circuit node of the superinductance is connected to ground by a stray capacitance C 0i . The N J + 1 node flux (phase) variables of the circuit are denoted by Φ i (φ i = Φ i /ϕ 0 ), where ϕ 0 = /2e is the reduced quantum of magnetic flux and i ∈ [0, N J ] [see also Fig. 1 (a)]. The circuit Lagrangian can then be written as [3] where Φ ext is the flux through the circuit loop. A more convenient basis is defined by the flux variables Θ The relation between the new modes and the original node fluxes can be expressed concisely by 0 · · · · · · 0 1 −1 0 0 · · · · · · · · · 0 1 −1 Under this change of basis, Eq. (A1) becomes where Note that the Σ mode does not enter in the potential energy. After a Legendre transformation, we arrive at the circuit Hamiltonian where q Θ ≡ ∂L(Θ,Θ)/∂Θ = C Θ ·Θ is a vector of conjugate charge operators, θ i = Θ i /ϕ 0 are phase operators and ϕ ext = Φ ext /ϕ 0 . In the presence of disorder in the circuit capacitances, the σ = Σ/ϕ 0 mode couples slightly to the θ i modes via the respective conjugate charge operators. Here, we neglect this capacitive coupling under the assumption of small circuit-element disorder and a large-frequency σ mode. The inverse capacitance matrix can thus be truncated to include only the θ i modes, i.e. C −1 Θ → C −1 Θ [0 : N J − 1, 0 : N J − 1], reducing Eq. (A4) to a Hamiltonian of N J interacting degrees of freedom. Note that the resulting pairwise θ i -θ j capacitive coupling has all-to-all connectivity and exhibits no particular structure in the θ i basis.

Accounting for charge dispersion
To model charge dispersion, we assume that each of the N J + 1 circuit islands is coupled to a local fictitious voltage source V i for i ∈ [0, N J ]. The associated terms in the Lagrangian can generically be written as where C gi is a gate capacitance for the ith circuit node. Equivalently, this can be expressed as where C g = diag(C g0 , C g1 , . . . , C g N J +1 ) and V = (V 0 , V 1 , . . . , V N J +1 ) T . In addition to Eq. (A5), the capacitance matrix of the circuit is modified to account for the gate capacitances as the conjugate charge operators are given by where ) is a constant of motion, and only N J of the N J + 1 offset charges are strictly independent. Using these expressions, the circuit Hamiltonian finally takes the form Omitting the σ mode and irrelevant constants, the above expression simplifies to are effective offset charges in the θ i basis and q i = [q Θ ] i . This Hamiltonian is equivalent to Eq. (3). Each of the bracketed terms in Eq. (A8) define a site Hamiltonian (H 0i for the ith array junction), while the remaining terms correspond to both linear and nonlinear all-to-all interactions between the sites. Note that, in the main text, we have used the Cooper-pair-number operators n i = q i /2e and the offset-charges n gi = q gi /2e, instead of q i and q gi , respectively.
We now derive an effective single-mode Hamiltonian for the fluxonium qubit that captures all circuit details. Because it is simple yet accurate, this model is used in the main text to assert the validity of the DMRG simulations in appropriate parameter ranges.
To obtain this effective model, we first consider a change of coordinates in which adiabatically eliminating the circuit modes other than the superinductance mode φ = N J i=1 θ i is simple. To find this appropriate change of coordinates, we reverse engineer the following Ansatz defining a new change of basis for k ∈ [1, N J − 1]. Note that Eq. (B1) acts as identity in the subspace of the σ mode and none of the σ-mode components of the capacitance matrix C Θ are included in Eq. (B2). The role of R (1) is to capacitively decouple a superinductance-like mode of the form from all other circuit modes, while leaving the σ mode invariant. Indeed, the new capacitance matrix with C (0) X = C Θ is block-diagonal in the absence of disorder. The first block has dimension 1 × 1 and corresponds to the φ (1) mode; the second block has dimension (N J − 1) × (N J − 1) and involves all circuit modes except φ (1) and σ; the last 1×1 block corresponds to the σ mode. By design, the first and second blocks of Eq. (B4) are exactly decoupled from each other, even in the presence of circuit-element disorder. In this case the first two blocks can be weakly coupled to the third block. Because the σ has a very high frequency for standard fluxonium circuit parameters, we neglect this coupling.
Similarly to R (1) , the matrix R (n) is composed by a n×n identity block for the modes labeled by k < n; a (N J −n+ 1)×(N J −n+1) block for modes labeled by k ∈ [n, N J −1]; and a 1 × 1 block for the σ mode. The coefficients {a which is a generalization of Eq. (B2).
The transformations R (n<N J −1) are designed to each decouple a single mode, while R (N J −1) decouples the last two modes n = N J − 1 and n = N J . Therefore, these N J − 1 successive transformations exactly diagonalize the upper N J × N J block of the capacitance matrix C Θ that does not include the σ mode. We can then invert these transformations arriving at the expression where the coefficient v ni quantifies how much the φ (n) mode couples to the ith Josephson junction of the array. Using Eq. (B7) and the definition φ = If C 0 = 0, it follows that V n = 0 for n ∈ [2, N J ], and φ (1) ≡ φ is the only mode that couples to the black-sheep junction. In other case, all modes are weakly coupled to the black-sheep junction, but this undesired coupling can be easily taken into account as we show in the following.
The relations Eq. (B7) and Eq. (B8) are now incorporated back to the potential energy of Eq. (A4). In order to trace out the unwanted degrees of freedom, we write the operator φ (n) for n > 1 in terms of the harmonic-oscillator ladder operators as φ (n) = √ πz n (a n + a † n ). Here, z n = L n /C n /R Q is the effective reduced impedance of the nth mode, given in terms of the effective inductance L n and capacitance C n . While C n can be readout directly from the block-diagonal capacitance matrix, the reduced inductance is determined by the product L −1 n = X T n · (M −1 ) T · L −1 · M −1 · X n , where X n is the mode vector associated to φ (n) and M = ( N J −1 n=1 R (n) ) T · R is a matrix that reverses the multiple changes of basis. The trace can then be performed straightforwardly by noticing that and thus tr n [e ixφ (n) ρ] = e −πx 2 zn/2 where we assume that the nth mode remains in its noninteracting vacuum state. Following to Eq. (B7) and Eq. (B8), we approximate n=2 e −πv 2 ni zn/2 , and cos(φ + ϕ ext ) tr n>1 [cos(φ + ϕ ext )] n=2 e −πV 2 n zn/2 . In Eqs. (B10) and (B11), tr n>1 indicates a trace operation over all circuit modes φ (n) , except for n = 1. Then, by renaming φ (1) → φ , we arrive at the effective single-mode Hamiltonian where E C is taken to be the charging energy E C = e 2 /2[C (1) X ] 00 of the φ mode and [φ , n ] = i. Note that Eq. (B12) is equivalent to Eq. (4) of the main text. Up to corrections of order N −3 J , Eq. (B12) reduces to are the effective inductive and Josephson-junction energies. Eq. (B13) corresponds to the original fluxoniumqubit model of Ref. [5]. Here, however, all energies entering Eq. (B13) are specified by a precise function of the circuit-element parameters.

Qualitative regimes of the fluxonium qubit
Despite the apparent simplicity of the effective fluxonium Hamiltonian Eq. (B13), its eigenstates can display a rich structure that depends on the parameter regime. For a systematic analysis, it is useful to redefine the parameters in Eq. (B13) in terms of the effective blacksheep junction plasma frequency ω b p = √ 8E J E C / and effective (reduced) impedance z b = π −1 2E C /E J . The potential energy of Eq. (B13) has a quadratic component given by the inductive term E L φ 2 /2 modulated by the cosine potential of the black-sheep junction and the external flux. Qualitatively, ω b p defines the characteristic energy of intra-well excitations within a given well defined by the Josephson potential, while z b is a measure of the tunneling amplitude between these wells. Figure 5 (a-c) shows the wavefunctions of the fluxonium qubit for different values of z b , taking ω b p /2π = 10 GHz and E L /h = 0.2 GHz constants and for Φ ext /Φ 0 = 0.35. Panel (a) corresponds to the case of a small effective impedance with z b = 0.1, in which tunneling between states localized in different wells is exponentially suppressed [27]. In this regime, the eigenstates of the fluxonium Hamiltonian are therefore localized within the deep potential wells of the potential-energy landscape. Excitations localized in a given potential well are approximately separated by the energy difference ω b p . For this reason, a transition between two of such states is called plasmon (or intra-well) transition. On the other hand, a transition between two states that belong to different potential wells is called fluxon (or inter-well) transition. Since the relative positions between potential wells shift significantly with Φ ext , fluxon transitions are highly sensitive to the external flux. In contrast, plasmon transitions are only weakly flux-sensitive. Since the lowimpedance limit requires the fluxonium mode φ to have a large effective capacitance (or "mass"), this regime is referred to as 'heavy-fluxonium' regime [6,27,35]. Figure 5 (b) shows an intermediate value of z b = 0.3, where the energy barrier (∝ E J ) between the potential wells due to the black-sheep junction has been reduced with respect to panel (a). Moreover, the effective capacitive energy E C has been increased, such that quantum tunneling between states localized in two neighboring potential wells is now non-negligible. This favors states that are delocalized across multiple potential wells and are the result of significant hybridization between plasmon and fluxon excitations. This intermediate regime for z b corresponds to the original fluxonium-qubit regime [5,33].
If the impedance of the black-sheep junction z b is increased further, the fluxonium wavefunctions can spread over many potential wells thanks to a lower E J and a larger E C . This situation is illustrated in Fig. 5 (c) where the distinction between plasmon and fluxon transitions is no longer useful and the spectrum is mostly determined by the harmonic part of Eq. (B13). The Josephson potential now leads to a weak flux sensitivity of the qubit transitions. Since the effective capacitance of the fluxonium mode needs to be lowered in order to make z b larger, this regime is known as the 'light-fluxonium' regime [51].
With the purpose of making the comparison above more precise, we now analyze qualitatively the energy spectrum of fluxonium devices from the heavy-to the light-fluxonium regimes. Fig. 5 (d) shows the result of the diagonalization of Eq. (B13) (light-blue lines) for the parameters of Fig. 5 (a). In this case, the low-frequency spectrum is highly sensitive to external flux, corresponding to a set of fluxon transitions. For a small z b , the lowfrequency spectrum around Φ ext /Φ 0 = 0.5 can be modeled by the weak coupling of two ground states {|m , |m+ 1 } with φ |m ∝ z Figure 5 (e) shows the frequency spectrum corresponding to the parameters in Fig. 5 (b). In this regime, the tunneling amplitude between different potential wells is stronger, leading to a larger hybridization gap between plasmon and fluxon transitions. However, for a moderate value of z b , the distinction between plasmon and fluxon transitions is still justified. As shown in Fig. 5 (f ), which corresponds to the spectrum associated to Fig. 5 (c), this distinction is no longer convenient to interpret the case of a large black-sheep junction impedance. Indeed, as z b is made significantly larger, plasmon and fluxon transitions undergo a very strong hybridization. In this limit, the fluxonium eigenstates become insensitive to external magnetic flux, leading to a reduced flux dispersion of the qubit transition [51,57].

Exploration of various parameter regimes
In this section, we provide further numerical evidence of the exceptional agreement between the DMRG simulations and the single-mode theory of appendix B 1. For this purpose, Fig. 6 shows an extension of the results in the main text, including the spectrum of a fluxonium device with N J = 180 array junctions and matrix elements of the phase and charge operators corresponding to the superinductance mode. As in the main body of the paper, the array junctions are modeled as multilevel systems including the first 15 eigenstates of the site Hamiltonian. To demonstrate that the agreement between these two approaches extends to all parameter sets for which the array junctions behave as weakly anharmonic oscillators, we compare Eq. (A4) and Eq. (B12) for various circuit design parameters. We also include the results obtained with an additional theory adapted from Ref. [35], where the nonlinearity of the array junctions is not taken into account. More precisely, we employ a single-mode approximation of the multimode Hamiltonian of Ref. [35], that we will refer to as 'linear theory'. The objective of this additional comparison is to highlight the effect of the nonlinearity of the array junctions which, as shown below, renormalizes the effective superinductance.
In particular, we test circuit Hamiltonians for various black-sheep junction capacitances (Fig. 7) and arrayjunction impedances (Fig. 8). The results of Fig. 7 demonstrate a very good agreement between the DMRG estimation (symbols) and the single-mode Hamiltonian Eq. (B12) [black dashed lines] from light-to heavyfluxonium parameter sets. Blue dotted lines correspond to the predictions of the linear theory. Overall, the latter estimations are in good agreement with the DMRG and the single-mode-theory results, although we find appreciable deviations for some of the flux-sensitive transitions. As we discuss below, these deviations are explained by the effect of the array-junction nonlinearity. Figure 8 shows a comparison between the results of DMRG, the single-mode theory and the linear theory for heavy-fluxonium parameter sets where the array-junction impedance is increased from z = 0.03 to z = 0.10. We observe that, in most cases, the single-mode theory of appendix B 1 provides an excellent estimation of the frequency of all fluxonium transitions determined by full DMRG. However, the predictive power of the effective single-mode theory weakens as z becomes larger (see Fig. 8 for z 0.08). We attribute this discrepancy to the unfavorable scaling of the multimode coupling in Eq. (A4) with z. This makes the approximation used to take the trace in Eq. (B10) and Eq. (B11) not completely justified. Although further refinement of the theory of appendix B 1 might be possible, the breakdown of the noninteracting approximation defines a parameter regime where the DMRG estimations are in principle out of reach of a simple theory.
Furthermore, Fig. 8 shows that the prediction of the linear theory of Ref. [35] [blue dotted lines] becomes increasingly inaccurate as z increases. Indeed, since the renormalization of the effective superinductance scales exponentially with the array-junction impedance [see Eq. (B10)], the frequency shifts of the qubit transitions due to the junction nonlinearity are more noticeable for larger z. These frequency shifts are more clearly appreciated for the flux-sensitive (or fluxon) transitions in Fig. 8, which are highly sensitive to the effective superinductance value. This also explains the relatively small deviations encountered in Fig. 7 between the linear theory and DMRG for z = 0.03.
Finally, we point out that we have also compared the result of the DMRG implementation to that of full exact diagonalization for fluxonium-like devices with a small number of junctions N J ∈ [2,6]. We find excellent agreement between the DMRG and the exact-diagonalization implementations for all circuit parameters, strengthening the validity of our DMRG algorithm. These numerical tests provide solid evidence of a successful DMRG In this section, we derive a master equation describing pure dephasing due to 1/f flux noise in the fluxonium qubit. Assuming weak system-bath coupling, the master equation is obtained from the standard integrodifferential equation where ρ(t) ⊗ ρ B is the system-bath density matrix, assumed to be separable at all times [58]. Assuming that the bath correlation functions are sharp around τ = 0, ρ(t − τ ) in Eq. (C1) can be approximated by ρ(t) with negligible error. This standard approximation conveniently leads to a Markovian master equation and allows us to extend the integral in Eq. (C1) to infinitely negative times. This last step is however not performed here in order to capture the Gaussian decay of the density matrix coherences in the presence of 1/f noise.
The system-bath interaction Hamiltonian can be obtained from the fluxonium circuit Hamiltonian assuming that Φ ext = Φ 0 ext + δΦ, where Φ 0 ext is the applied flux bias and δΦ represents fluctuations. To first order in δΦ, the interaction Hamiltonian can be written as [59] where H is the Hamiltonian of the fluxonium qubit and the derivative with respect to the external flux is evaluated at Φ ext = Φ 0 ext . Expanding Eq. (C1) in the eigenbasis {|ψ k } of the full circuit, we arrive at where we have introduced the matrix elements ∂ Φext H| kk Φ 0 ext = ψ k |∂ Φext H| Φ 0 ext |ψ k , and omitted the explicit time dependence of ρ(t) → ρ.
Tracing out the bath degrees of freedom leads to the socalled Bloch-Redfield equation [58]. This equation has, however, a number of disadvantages that can potentially lead to unphysical dissipation results. Thus, for prac-tical purposes, we use the rotating-wave approximation discarding terms for which ω ll + ω kk = 0. As shown below, this approximation reduces Eq. (C3) to a Lindbladform master equation. Assuming that the qubit has a set of nondegenerate energy transitions, this approximation is equivalent to the conditions l = k and l = k for ω kk = 0, and l = l for ω kk = 0. In this way, Eq. (C3) simplifies to We now assume that δΦ(t) can be modeled as a (real) stationary random process. This assumption is motivated by physical models of bistable two-level-system defects that are known to produce noise of type 1/f [60,61]. Furthermore, we make the usual assumption that the weight of the 1/f noise spectral density is negligible at the qubit transition frequencies such that it does not significantly contribute to the device's T 1 time. The pure-dephasing master equation is therefore derived from the third line and assume the general form where A Φ is the 1/f flux-noise amplitude, typically reported to be in the range 1 − 10 µΦ 0 [25]. It must be stressed that Eq. (C7) is an approximation to the spectral densities measured in the laboratory, which can scale as |ω| −µ with µ ∈ [0.6, 1.3] [25,62]. We proceed further by exploiting a simple mathematical fact. Using Eq. (C6) and Eq. (C7), we find that Here, ω ir is an infrared frequency cutoff in the order of 2π×1 Hz, introduced to regularize the cosine integral and motivated by physical reasons [63]. Since the time t over which we are interested in calculating the time evolution of the density matrix is small compared to the time scale set by ω −1 ir , we make use of the series expansion Ci(w) = γ + log(y) + ∞ k=1 (−y 2 ) k 2k(2k)! , where γ 0.58 is the Euler's constant, approximating Expanding the double commutators in Eq. (C5) and making use of Eq. (C10), we arrive at a pure-dephasing master equation of the form where Γ kl ϕ are time-dependent pure-dephasing rates given by σ kl = |ψ k ψ l |, and D[x, y] ρ = xρy † − {y † x, ρ}/2 is a generalized dissipator superoperator.
Equivalently, Eq. (C11) can be recast in the more familiar form where D[x] ρ = xρx † − {x † x, ρ}/2 is the standard dissipator superoperator. By projecting Eq. (C13), one has Thus, we verify that the decay of the coherences of the density matrix is proportional to the flux dispersion of the k ↔ l qubit transition, as expected for first-order dephasing processes. Since second-order corrections to the pure-dephasing rate at sweet spots are of order A 4 Φ and vanishing small, most devices are T 1 -limited at such operating points. Now, in order to produce an estimate of the pure-dephasing coherence time due to 1/f flux noise, we simply integrate Eq. (C14), arriving at the expression We note that expressions similar to Eq. (C16) have been derived previously in the literature [25,63]. However, these expressions do not include the correction ( 3 2 − γ) within brackets in Eq. (C16). Finally, we define the co-herence time T ϕ as the solution of the implicit equation ρ 01 (T ϕ )/ρ 01 (0) = 1/e. The solution of this equation has been used in Fig. 4 to produce an estimation of the pure-dephasing coherence times due to flux noise.