Efficient modeling of superconducting quantum circuits with tensor networks

We use a tensor network method to compute the low-energy excitations of a large-scale fluxonium qubit up to a desired accuracy. We employ this numerical technique to estimate the pure-dephasing coherence time of the fluxonium qubit due to charge noise and coherent quantum phase slips from first principles, finding an agreement with previously obtained experimental results. By developing an accurate single-mode theory that captures the details of the fluxonium device, we benchmark the results obtained with the tensor network for circuits spanning a Hilbert space as large as 15180. Our algorithm is directly applicable to the wide variety of circuit-QED systems and may be a useful tool for scaling up superconducting quantum technologies.


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. Superconducting circuits operate at milliKelvin temperatures, where the 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 lowlying 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 task rapidly becomes intractable by exact diagonalization. With current devices integrating 10s 5 to 100s 6 , 1000s 7 and even 10,000s 8 Josephson junctions, finding methods to efficiently model these systems is timely.
In this work, we adapt to many-body superconducting quantum circuits a numerical tensor network method that we introduce in (Baker, T. E. et al. manuscript in preparation). As an example of this method, we use this numerical technique to compute the relevant low-energy excitations of a large-scale superconducting circuit known as fluxonium qubit 5 , taking into consideration all degrees of freedom of a circuit model of the device. Solving the complete fluxonium Hamiltonian is challenging due to the presence of longrange linear and nonlinear interactions between circuit modes. We show that the tensor network technique succeeds at modeling this circuit and gives access to information about the system that can be used, for instance, to estimate the qubit coherence times from first principles. Moreover, we develop an effective theory for the fluxonium qubit that we use to benchmark the results obtained with the tensor network algorithm in certain parameter regimes.
To highlight the benefits of our tensor network approach, we study the phenomenon of charge dispersion in the fluxonium qubit. To this end, we consider the device parameters reported in a experiment that analyzed the effect of charge dispersion on the coherence time of the qubit 9 . We perform numerical simulations in a range of parameters, including those of the experiment, confirming an existing effective theory for describing charge dispersion in this qubit, and clarifying the regime of validity of that theory. Moreover, we estimate the coherence time of the experimental device using our tensor network approach, finding agreement with the reported values. Finally, we push the tensor network method to parameter regimes where the accuracy of the effective theory is reduced, and quantify the deviations with respect to the tensor network result. Our findings are useful in the design of fluxonium devices with balanced fluxand charge-noise dephasing rates, and motivate further experimental work to probe charge dispersion in these circuits.
The manuscript is organized as follows. In section "The multitargeted DMRG algorithm", we summarize some important aspects of the tensor network technique employed in this work. In section "DMRG implementation of the fluxonium-qubit Hamiltonian", we provide a tensor network implementation of the complete fluxonium-qubit Hamiltonian, introduce an accurate single-mode description of this qubit, and compare the results obtained with these approaches as a way of benchmarking the tensor network results. Section "Charge dispersion and coherence time" discusses charge noise in the fluxonium qubit, and provides tensor network estimates of charge dispersion of the qubit transition supporting a previously developed theory 9 . We conclude in "Discussion" section and provide an outlook of our results.

RESULTS
The multi-targeted DMRG algorithm A useful strategy to determine the low-energy excitations of a quantum system is based on decomposing its many-body wavefunction into a series of tensors, each representing a single site (or mode). The resulting wavefunction is known as matrix product state (MPS) 10 . For a review, see for instance refs. [11][12][13][14][15] .
The full many-body wavefunction can be written as 12 where σ i indexes orbitals (or levels) that belong to a finitedimensional 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 ψ j i, a series of tensor decompositions can be performed using the singular value decomposition (SVD). By performing successive SVDs on the original tensor, one obtains a site-by-site representation of the wavefunction of the form 12 where A σi aiÀ1ai is a tensor associated to the ith site. Here, an extra index a i appears corresponding to a link that connects to an adjacent site. The dimension of this shared index is known as the bond dimension and is controlled by truncating the number of nonzero singular values that are kept following the SVDs. Effectively, this truncation leads to an approximate representation of the many-body state. Physical systems that can be modeled efficiently by a MPS often involve short-range interactions and low dimensions 16 . More general cases can also be captured by a MPS at the price of using a larger bond dimension 11,12 . 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 17,18 . This variational approach is efficient at solving for the ground state of systems that are well captured by the MPS and can converge in only a few iterations 16,19,20 . More importantly, the complexity of this algorithm scales linearly with the number of sites, making it possible to treat system sizes well beyond what is possible with exact diagonalization.
While DMRG is most commonly used to study ground states, 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 hω 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 P i2ex: Λ ψ i j i ψ i h j, with Λ > 0, where ex. denotes a set of previously determined excitations f ψ i j ig. This energy penalty turns the next excited into an effective ground state for which standard DMRG can be run 12 . However, this technique can in some cases miss excited states and suffers from convergence issues.
To remedy this problem, we have derived an extension of the DMRG algorithm that can find a number of excitations simultaneously, see (Baker, T. E. et al. manuscript in preparation). This procedure, that we name the "multi-targeted" DMRG algorithm, extends Eq. (2) to include an additional index that labels the excitation number. Energy is then variationally minimized to the correct eigenvalues and simultaneously for all excitations. This technique is numerically stable and, unlike standard approaches, does not miss excitations or introduce numerical degeneracies in all tested situations, including tens of excitations. Furthermore, an important benefit of our multi-targeted DMRG algorithm is that the orthogonality of the computed excited states is guaranteed up to numerical accuracy.
DMRG implementation of the fluxonium-qubit Hamiltonian Because of its relatively complex structure, with a Hamiltonian that includes periodic boundary conditions, as well as short-and longrange linear and nonlinear interactions, the fluxonium qubit 5 is a challenging system to explore with our numerical method. We note, however, that non-multi-targeted DMRG has previously been used to study quantum phase transitions in Josephson junction rings 21,22 and the coherence properties of the current-mirror qubit 23 .
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 effective impedance greater than the quantum of resistance R Q = h/(2e) 2 ≃ 6.5 kΩ and self-resonance frequencies >10 GHz (refs. [24][25][26][27][28] ). The inductive shunt makes the device insensitive to charge noise 29 , while fluxnoise insensitivity can be reached for large values of the superinductance 30 . Recent realizations of the fluxonium qubit have shown long coherence times 31 , making these devices very attractive for quantum information processing.
Superinductances are most commonly fabricated by connecting 10s to 100s of Josephson junctions in series, forming a linear array 5,27 . Other promising implementations of superinductances are based on nanowires built from disordered superconductors, Fig. 1 Circuit model of the fluxonium qubit. a Detailed circuit scheme including a "black-sheep" junction (center) shunted by a capacitance (top) and a junction-array superinductance with N J 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, N J ] denotes the superconducting phase at every circuit node, while θ i for i ∈ [1, N J ] 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:  [32][33][34][35] . Beyond fluxonium, superinductances are also exploited by other qubit designs, with the notable example of the noise-protected 0 − π qubit 25,36 . Regardless of the physical implementation, a superinductance is a multimode device that can in principle be approximated by a single-mode circuit element 27 . While the single-mode description of the superinductance is often enough to describe a qubit that incorporates this circuit element, the underlying multimode structure of the former device is inherited by the qubit Hamiltonian and can have important consequences 27,32,37,38 . Below, we investigate some of these aspects for the case of a fluxonium qubit based on a Josephson junction array. However, it is worth noticing that our theory and numerical methods are applicable to a wider variety of fluxonium devices, as nanowire superinductances can effectively be described as an array of Josephson junctions 33,39 .
Setting up the multi-targeted DMRG algorithm. With the goal of determining the low-energy excitations of the full fluxonium device shown in Fig. 1a by means of the multi-targeted DMRG algorithm, we first describe the circuit Hamiltonian. In this circuit, the black-sheep junction is specified by its Josephson energy E J b and its capacitance C Jb , 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 plasma frequency ω p ¼ 1= ffiffiffiffiffiffiffiffi =R Q of the array junctions. Following the standard procedure for quantizing a circuit 3 , the Hamiltonian of the circuit takes the form (see section "Fluxonium circuit hamiltonian") In this expression, H 0i ¼ 4E Ciñ 2 i À E Ji cos θ i , withñ i ¼ n i À n g 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 g i 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 /ñ iñj arising from the ground, black-sheep and array-junction capacitances, that couples all sites with comparable strength and full connectivity. The last term of Eq. (3) is a nonlocal interaction that results from the strongly nonlinear Josephson potential of the black-sheep junction, and that depends on the external flux Φ ext = φ 0 φ ext associated with this junction 3 . Because Eq. (3) includes a very large number of degrees of freedom and is therefore difficult to work with, this Hamiltonian is not directly employed in the literature to describe the fluxonium qubit. Instead, fluxonium devices are usually modeled by an effective Hamiltonian that incorporates a single bosonic degree of freedom, ϕ ¼ P NJ i¼1 θ i , known as superinductance or fluxonium mode 5 .
To go beyond the usual effective model and obtain the lowenergy excitations of Eq. (3) by means of a tensor network, the circuit Hamiltonian must first be converted to its MPO form. Crucially, we noticed that the long-range cosine interaction is ideally suited to MPSs 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, including the black-box-quantization [40][41][42] and energy-participation-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 multi-targeted DMRG algorithm, as the results reported in this manuscript are obtained with a relatively small bond dimension using MPO compression techniques 44 . The efficient matrixproduct-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.
Effective single-mode theory. To assert the validity of our DMRG method, we derive in section "Effective model for the fluxonium qubit" an effective single-mode theory from Eq. (3) that can be solved by exact diagonalization. This theory goes beyond the standard treatment found in the literature, accounting for the details of a circuit model of the device and some of the effects due to the multimode structure. Under approximations controlled by the parameter regime of the device, we arrive at the Hamiltonian where the mode described by ϕ 0 is closely related to the superinductance (or fluxonium) mode ϕ, and n 0 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 ϕ 0 ¼ ϕ and where n is the conjugate charge operator to ϕ. Otherwise, the ϕ 0 mode includes corrections to ϕ that are linear in C 0 . We note that Eq. (4) is a single-mode Hamiltonian obtained after tracing out the high-frequency modes of the fluxonium circuit that are assumed to be in vacuum. This constitutes an approximation useful to describe the lowfrequency properties of the fluxonium qubit.
Although in the limit of large N J Eq. (4) reduces to the usual effective model for the fluxonium qubit (see Fig. 1b) 5 , the parameters of Eq. (4) include important corrections that arise from the nonlinearity of the array junctions and the details of the circuit. We find, in fact, that these considerations can lead to significant frequency shifts of the qubit transitions (see Supplementary Fig. 3). Importantly, our theory is formulated for an arbitrary capacitance matrix and just like our numerical techniques, it is readily applicable to circuit models that generalize that of Fig. 1, for instance by including additional stray capacitances needed to describe a given device. Crucially, because of its singlemode nature, Eq. (4) can easily be diagonalized numerically by truncating the Hilbert space of the ϕ 0 mode to finite dimension.

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 explore the capabilities of this method. To this end, we consider a device in the "heavy fluxonium" regime 6,32,45 , where the shunt capacitance is large, and a superinductance made of N J = 120 identical junctions, where ω p /2π = 25 GHz and z = 0.03 (ref. 27 ). We choose to work in the heavy regime for a demonstration because of the recent interest in the heavy fluxonium due to its millisecond-long T 1 . It should be stressed, however, that there is nothing particular about this regime from a numerical point of view, and we show in section "Charge dispersion and coherence time" and the appendices of this work how the tensor network is equally successful at solving the fluxonium Hamiltonian in other parameter regimes. A comprehensive description of the various Each junction is modeled as a multilevel system using the 15 lowest energy eigenstates of the site Hamiltonian H 0i . Indeed, we find that for low-impedance junctions, truncation errors can be avoided using a smaller number of eigenstates of H 0i compared to other possible basis, such as the Cooper-pair-number basis 46 . The DMRG implementation is thus defined in a product basis of local wavefunctions spanning a many-body 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. Figure 2a shows the energy spectrum of the fluxonium device of Fig. 1 for both multi-targeted DMRG (Eq. (3), light-blue circles) and exact diagonalization of the effective single-mode theory (Eq. (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 many more junctions. Indeed, Supplementary Fig. 1 shows the result for a circuit with N J = 180 spanning a many-body Hilbert space of size 15 180 , while Supplementary Figs. 2 and 3 explore the result in distinct parameters regimes of the fluxonium Hamiltonian for moderate system sizes. These results provide supporting evidence of a successful DMRG implementation of the fluxonium-qubit Hamiltonian and motivate applying the DMRG technique in regimes of parameters, where deriving an effective model is challenging.
Exploring the DMRG results. In addition to computing global properties of the circuit, such as its energy spectrum, the multitargeted 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. 2b shows the mean photonnumber population hp i i ¼ hψ k jH 0i jψ k i=_ω p of the ith site, for all sites (i ∈ [1, 120], vertical axis of each of the six density plots) as a function of Φ ext . These expectation values are computed for a given eigenstate ψ k j i of the full fluxonium circuit, from the ground state (k = 0, bottom density plot) to the fifth 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. 2a corresponds to that of the transition ψ 0 j i ! ψ i j 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. 2c, 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 j i 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. 2c (right panel).
The above interpretation becomes clearer by considering the effective potential and wavefunctions obtained from the singlemode effective Hamiltonian Eq. (4), as shown in Fig. 2d for Φ ext ∈ {0, Φ 0 /4, Φ 0 /2}. The shape of the effective potential is determined by the cosine potential of the black-sheep junction and the inductive energy ÀN 2 J E L cosðϕ 0 =N J Þ ' E L ϕ 0 2 =2 of the array. 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 intra-well excitations (see Supplementary Discussion). The potential of the effective model connects to that of Fig. 2c by noticing  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. 2e, which shows the expectation value of the phase drop ϕ 0 À ϕ i P i j¼1 θ j , obtained from DMRG, and plotted as a function of the site number for the fluxon states ψ 0 j i and ψ 2 j i at Φ ext = Φ 0 /4 in Fig. 2d (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 and it can also provide information that is not accessible from this theory. The close comparison shown in the results from both methods demonstrates a correct DMRG implementation of the full circuit Hamiltonian of the fluxonium qubit. This fact also suggests that other circuit Hamiltonians can benefit from this numerical method. Moreover, the local physical quantities such as those illustrated in Fig. 2b, contain information about the energyparticipation 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 circuitelement disorder. We return to these aspects in "Discussion" section.

Charge dispersion and coherence time
We now proceed with a concrete application that shows how our DMRG implementation can be leveraged to estimate coherence time of the fluxonium qubit from first principles. In particular, we are interested in quantifying the impact of charge noise and coherent quantum phase slips (CQPSs) on the qubit coherence time 9,28 .
Charge dispersion. In the fluxonium circuit, the black-sheep junction acts as a weak link that couples flux states of the superconducting, loop allowing for the control of the qubit's degree of freedom. The tunneling of a flux quantum corresponds to a change of 2π in the phase of the superconducting order parameter and is therefore known as CQPS 9,47-52 . The rate at which a quantum of flux can tunnel in and out of the loop through the black-sheep junction defines the CQPS amplitude, which increases with the impedance of this circuit element. In experiments, fluxonium devices exploit a wide range of blacksheep junction impedances, ranging from relatively small values in the heavy-fluxonium regime 6,32,45 , to moderate and large values in the fluxonium 5,9 and light-fluxonium 30 regimes, respectively. Yet, in practice, the loop's superinductance is also built from Josephson junctions itself and CQPS can develop through the junctions of the array. Unlike the previous case, however, CQPS occurring in the superinductance can degrade qubit coherence and must be minimized by circuit design 28 .
Manucharyan et al. 9 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 blacksheep junction are captured by a single-mode Hamiltonian similar in spirit to Eq. (4). In addition, CQPS occurring through the superinductance enter in the same effective Hamiltonian via the external flux, such that Φ ext → Φ ext + m Φ 0 , where m is an integervalued operator representing the number of CQPS in the array. A CQPS event at any junction of the superinductance leads to a jump m → m ± 1, and is described by a flux-tunneling Hamiltonian of the form where the operator m − [m þ ¼ ðm À Þ y ] removes (adds) a single flux quantum from the loop. Here, E S is the total CQPS amplitude given by the superinductance and obtained, as an approximation, by coherently adding the contributions from each of the N J array junctions. In the limit of small array-junction impedance and large inductance, one has E S ¼ P NJ i¼1 ϵ 0i e i2πng i (refs. 28,9 ), where corresponds to 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 ω p i of the ith array junction 28,46,47,53 . We note that, as a consequence of the Aharonov-Casher effect, the flux-tunneling amplitude contributed by the ith array junction has a well-defined phase, proportional to the offset-charge n g i associated to that junction 9,47,52,54-56 .
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 δω ij ¼ Re½E S ðhTi j À hTi i Þ=_ of the qubit's i → j transition frequency, where T ¼ expðÀi2πnÞ is a 2π displacement operator whose expectation values are computed, using the unperturbed eigenstates f ψ i j ig with m = 0 (ref. 9 ). For a homogeneous array (ϵ 0i ϵ 0 for i ∈ [1, N J ]), one has ÀN J ϵ 0 Re½E S N J ϵ 0 , 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. Quite interestingly, in this regime, the qubit frequency can be a sensitive probe of many-body phenomena: CQPS interference. 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. 9 . Figure 3a 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 g i n g , assumed to be the same on every junction of the array. Each subpanel shows the DMRG results for a given value of the arrayjunction 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 that at n g = 0 (ref. 57 ). Dashed black lines show the qubit transition frequency according to Eq. (4). Note that the offset-charge 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 the arrayjunction impedance z due to the exponential scaling of Eq. (5). This is best illustrated by Fig. 3b, 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 and n g = 0.5. Black triangle symbols (Eq. (6) (DMRG)) are the result of Eq. (6) for which the matrix elements are evaluated, using the eigenstates obtained from DMRG for n g = 0. In contrast, the black dashed line (Eq. (6) (single mode)) is obtained by evaluating the matrix elements using the single-mode Hamiltonian Eq. (4). We find no significant difference between the DMRG (Eq. (6) (DMRG)) and the single-mode (Eq. (6) (single mode)) implementations of Eq. (6) in the entire range of z.
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 in Fig. 3, which corresponds to a junction impedance Z J ≃ 650 Ω or one tenth of the superconducting quantum of resistance. This provides evidence in support of the theoretical model introduced in ref. 28 . Although barely visible, small deviations between the fully numerical DMRG estimation and those based on A. Di. Paolo et al.
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-behaved. 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.
The difference between the full numerical multi-targeted DMRG estimation and those based on Eq. (6) increases in the range of z ≳ 0.1 in Fig. 3 (see also figure inset). In this regime, some of the assumptions on which the theory of ref. 9 is based are no longer valid, placing the DMRG method at a clear advantage over the effective theory. Furthermore, the fact that such deviations can be quantified with a tensor network method motivate further experimental and theoretical explorations to understand the physics of CQPS in fluxonium devices to a greater extent.
Coherence-time estimations. Because of unavoidable charge noise, the value of δω ij fluctuates in time, broadening the qubit transition, and deteriorating coherence. This observation is the basis of the experimental study of ref. 9 , where the decrease of the qubit coherence time around the flux sweet spot is taken as indirect evidence of CQPS events in the superinductance. To quantify this effect, ref. 9 assumes that the variables n g i are independent and randomly distributed. The probability distribution of Re½E S can then be approximated by a Gaussian with zero mean and standard deviation ffiffiffiffiffiffiffiffiffi ffi N J =2 p ϵ 0 (ref. 9 ). The effective broadening of the qubit transition in presence of noise thus scales as ffiffiffiffiffi N J p , translating to the pure-dephasing rate 1=T φ;CQPS ¼ jΔω 01 j=4 ffiffiffiffiffi N J p (refs. 9[,28 ). 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. 9 predicts the puredephasing coherence times to be dominated by the combined effect of charge noise and CQPS around Φ ext = Φ 0 /2. Indeed, we compare T φ,CQPS to the coherence time expected for 1/f flux noise by deriving a multilevel pure-dephasing master equation of the form (see section "Multilevel pure-dephasing master equation for flux noise") where Γ kl φ are time-dependent pure-dephasing rates proportional to the 1/f flux-noise amplitude, σ kl ¼ ψ k j i ψ l h j, and D½x; y ρ ¼ xρy y À fy y x; ρg=2 is a generalized dissipator operator. By integrating Eq. (7), we find the flux-noise coherence time T φ,Flux by solving the implicit equation ρ 01 (T φ,Flux )/ρ 01 (0) = 1/e. Our semi-analytical method leads to a slight correction to the pure-dephasing coherence time with respect to other approaches, see section "Multilevel pure-dephasing master equation for flux noise" for details. Figure 4a shows the energy spectrum of the simulated device as a function of the external flux, results that should be compared to those of ref. 9 . In contrast to the results in Fig. 2a, the difference between the DMRG and single-mode simulations for the parameters of ref. 9 is slightly 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. This makes any single-mode approximation invalid, except at low frequencies. Furthermore, Fig.  4b 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 coherence time values that are very similar to the experimental observation (see Fig. 4 in ref. 9 ), 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. 9 .
Combined, the results of Figs. 3 and 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.

DISCUSSION
We have reported the application of a DMRG algorithm to simulate large-scale superconducting quantum devices. We have used this numerical technique to study aspects of quantum coherence of the fluxonium qubit. To assert the validity of the DMRG simulations and interpret the numerical results, we have developed a detailed single-mode theory for the fluxonium qubit. We have employed this theory and the numerical method to investigate the combined effect of charge noise and CQPSs in the fluxonium qubit. Our results on the charge dispersion of the fluxonium qubit confirm a model introduced in ref. 9 , and reproduce some of the experimental findings of that work. Combined, these results are of significant value for the design of the next generation of fluxonium and other many-body superconducting quantum devices.
Having access to the expectation values of local and of n-body operators thanks to tensor network techniques makes it possible to investigate the many-body properties of superconducting quantum circuits. This could help, for instance, to nonlocally encode quantum information 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.

Fluxonium circuit Hamiltonian
In this section, 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 Jb (including both shunt and junction capacitances) and Josephson energy E Jb 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. 1a). 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 these modes and the original node fluxes can be expressed concisely by Θ = R ⋅ Φ, where Θ ¼ ðΘ 1 ; ; Θ NJ ; ΣÞ T , Φ ¼ ðΦ 0 ; ; Φ NJ Þ T , and R is the N J + 1 × N J + 1 matrix Under this change of basis, Eq. (8) 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, reducing Eq. (11) to a Hamiltonian of N J interacting degrees of Fig. 4 Coherence time of a 40-junction superinductance fluxonium qubit. a Energy spectrum according to DMRG and single-mode estimations as a function of Φ ext . The black dotted line corresponds to the plasma frequency of the array junctions. b Pure-dephasing coherence times for flux and charge (CQPS) noise as obtained from DMRG. Parameters: C Jb ¼ 7:5 fF, E Jb =h ¼ 8:9 GHz, z = 0.09, ω p /2π = 12.5, and C 0 = 0, extracted from ref. 9 . The 1/f flux-noise amplitude is taken to be A Φ = 1.2 μΦ 0 , which is a conservative value 46 .
A. Di. Paolo et al.
freedom. Note that the resulting pairwise θ i -θ j capacitive coupling has allto-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 capacitance for the ith circuit node. Equivalently, this can be expressed as where C g ¼ diagðC g 0 ; C g 1 ; ; C g N J þ1 Þ and V ¼ ðV 0 ; V 1 ; ; V NJþ1 Þ T . In addition to Eq. (12), the capacitance matrix of the circuit is modified to account for the gate capacitances as the conjugate charge operators are given by where e The circuit Hamiltonian then takes the form Omitting the σ mode and irrelevant constants, the above expression simplifies to where Θ ii for i ∈ [1, N J ] 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. (15) define a site Hamiltonian (H 0i for the ith array junction), while the remaining terms lead to linear and nonlinear allto-all interactions between the sites. In the main text, we have used the Cooper-pair-number operators n i = q i /2e and the offset-charge n g i ¼ q g i =2e, instead of q i and q g i , respectively.

Effective model for the fluxonium qubit
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 ϕ ¼ P NJ i¼1 θ i is simple. To find this appropriate change of coordinates, we reverse engineer the following Ansatz defining an additional change of basis where the constants fa ð1Þ k g are defined by for k ∈ [1, N J − 1]. Note that Eq. (16) acts as identity in the subspace of the σ mode and none of the σ-mode components of the capacitance matrix C Θ are included in Eq. (17). The role of R (1) is to capacitively decouple a superinductance-like mode of the form k ðθ k À θ 1 Þ; from all other circuit modes, while leaving the σ mode invariant. Indeed, the capacitance matrix now reads where C ð0Þ X ¼ C Θ is block diagonal. Importantly, this transformation works in the presence of circuit-element disorder, in which case the blockdiagonal form of C ð1Þ X is preserved. Only spurious couplings to the σ mode are not purposely eliminated, as these will be neglected later on. 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. (19) 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 σ mode has a very high frequency for standard fluxonium circuit parameters, we neglect this coupling.
While the transformation Eq. (16) isolates the most relevant mode of the circuit, we iterate recursively this procedure to decouple all remaining circuit modes in the capacitive interaction. Doing this will allow us to trace out such degrees of freedom later on. We proceed by defining an additional set of rotation matrices {R (n) }, for n ∈ [2, N J − 1], with the general form which is a generalization of Eq. (17). The transformations R ðn < NJÀ1Þ are designed to each decouple a single mode, while R ðNJÀ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 stress that these transformations work even in the presence of disorder in the capacitance matrix, eliminating all but the coupling to the σ mode. Indeed, it can be seen that these transformations recursively block diagonalize any capacitance matrix, where the coefficients fa ðnÞ k g should be determined for each case. This procedure is also a key difference between our strategy and previous approaches to finding multimode Hamiltonians, such as those in refs. 37,38 . In these works, the Hamiltonian of the fluxonium circuit is expanded in a predefined analytical basis of socalled difference modes. This basis is then used to analyze the effects of circuit-element disorder and multimode interactions by performing a series expansion of both the kinetic and potential energy parts of the Hamiltonian. In our case, we find a basis that, while accounting for disorder in the capacitance matrix of the circuit, diagonalizes the kinetic energy part of the Hamiltonian (with exception of the σ-mode components). The upside of our approach is that following the basis transformation, all multimode interactions are moved to the potential energy, making it relatively easier to eliminate high-frequency modes and to obtain a singlemode approximation. A downside is that our basis transformation is not analytical and depends on the parameters of the circuit. In contrast to former analytical approaches 37,38 , our theory is thus semi-analytical.
Following this change of basis, we invert these transformations to arrive at A. Di. Paolo et al.
where the coefficient v ni quantifies how much the ϕ (n) mode couples to the ith Josephson junction of the array. The differences with the theory in refs. 37,38 are more noticeable from this expression, where the modes {ϕ (n) } and the coefficients {v ni } in terms of which the single-junction operator θ i is written depend on the details of the circuit. Using Eq. (22) and the definition ϕ ¼ P NJ i¼1 θ i , we moreover have where V n ¼ P NJ i¼1 v ni . 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 potential energy of Eq. (11) is now rewritten using the relations Eqs. (22) and (23). 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Þ ¼ ffiffiffiffiffiffiffi πz n p ða n þ a y n Þ. Here, z n ¼ ffiffiffiffiffiffiffiffiffiffiffiffi L n =C n p =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 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 Eqs. (22) and (23), we approximate where with x b ¼ Q NJ n¼2 e ÀπV 2 n zn=2 . In Eqs. (25) and (26), tr n > 1 indicates a trace operation over all circuit modes ϕ (n) , except for n = 1. The trace operation introduces important corrections to the circuit Hamiltonian that vary exponentially with the effective impedance {z n } of the circuit modes and are a consequence of the full-cosine structure of the array junctions' potential. To the best of our knowledge, these corrections are not reported elsewhere. By renaming ϕ ð1Þ ! ϕ 0 , we arrive at the effective single-mode Hamiltonian where E C is taken to be the charging energy X 00 of the ϕ 0 mode and ½ϕ 0 ; n 0 ¼ i. Note that Eq. (27) where E L ¼ P NJ i¼1 x i E Ji =N 2 J and E J ¼ x b E Jb are the effective inductive and Josephson junction energies. Equation (28) corresponds to the original fluxonium-qubit model of ref. 5 . Here, however, all energies entering Eq. (28) are specified by a precise semi-analytical function of the circuit parameters.
Finally, we note that, while a single-mode theory is enough for the purposes of this work, the multimode properties of the fluxonium qubit could, in principle, be investigated using our theory by undoing the trace operation. However, these properties have been extensively studied before using other methods 37,38 .

Multilevel pure-dephasing master equation for flux noise
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 dτ tr B ½H int ðtÞ; ½H int ðt À τÞ; ρðt À τÞ ρ B ; (29) 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. (29) 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. (29) to infinitely negative times. This last step is however not performed here in order to capture the Gaussian decay of the coherences of the density matrix 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. (29) in the eigenbasis f ψ k j ig of the full circuit, we arrive at where we have introduced the matrix elements ∂ Φext Hj kk 0 Φ 0 ext ¼ hψ k j∂ Φext Hj Φ 0 ext jψ k 0 i, and omitted the explicit time dependence of ρ(t) → ρ.
Tracing out the bath degrees of freedom leads to the so-called Bloch-Redfield equation 58 . This equation has, however, a number of disadvantages that can potentially lead to unphysical dissipation results. Thus, for practical purposes, we use the rotating-wave approximation discarding terms for which ω ll 0 þ ω kk 0 ≠0. As shown below, this approximation reduces Eq. (31) to a Lindblad-form master equation. Assuming that the qubit has a set of nondegenerate energy transitions, this approximation is equivalent to the conditions l ¼ k 0 and l 0 ¼ k for ω kk 0 ≠0, and l ¼ l 0 for ω kk 0 ¼ 0. In this way, Eq. (31) simplifies to ∂tρ ¼ À 1 We now assume that δΦ(t) can be modeled as a (real) stationary random process. This assumption is motivated by physical models of bistable twolevel system defects that are known to produce noise of type 1/f (refs. 60,61 To obtain a workable expression, we introduce the noise spectral density S 1=f Φ ½ω for 1/f flux noise by the definition 3 tr B ½ρ B δΦðtÞδΦðt 0 Þ ¼ 1 2π and assume the general form where A Φ is the 1/f flux-noise amplitude, typically reported to be in the range 1-10 μΦ 0 (ref. 46 ). It must be stressed that Eq. (35) is an approximation to the spectral densities measured in the laboratory, which can scale as |ω| −μ with μ ∈ [0.6, 1.3] 46,62 . We proceed further by exploiting a simple mathematical fact. Using Eqs. (34) and (35), we find that Z t 0 dτ tr B ½ρ B δΦðtÞδΦðt 0 Þ ¼ lim ωir!0 À 2A 2 Φ Z t 0 dτ Ciðω ir τÞ; where CiðyÞ ¼ À R 1 y dx x À1 cos x is the cosine integral. Here, ω ir is an infrared frequency cutoff in the order of 2π × 1 Hz, introduced to regularize A. Di. Paolo et al.
the cosine integral and motivated by physical reasons 63 . Since the time t in 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 where γ ≃ 0.58 is the Euler's constant to approximate Z t 0 dτ tr B ½ρ B δΦðtÞδΦðt 0 Þ ' 2A 2 Φ t ½ð1 À γÞ À log ðω ir tÞ: (38) Expanding the double commutators in Eq. (33) and making use of Eq. (38), we arrive at a pure-dephasing master equation of the form where Γ kl φ are time-dependent pure-dephasing rates given by σ kl ¼ ψ k j i ψ l h j, and D½x; y ρ ¼ xρy y À fy y x; ρg=2 is a generalized dissipator superoperator. Equivalently, Eq. (39) can be recast in the more familiar form where D½x ρ ¼ xρx y À fx y x; ρg=2 is the standard dissipator superoperator. By projecting Eq. (41), one has hψ k j∂ t ρjψ l i ¼ À where Thus, we obtain 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 Φ , most devices are simply 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 integrate Eq. (42), arriving at the expression ρ kl ðtÞ ¼ ρ kl ð0Þ exp ÀA 2 Φ ð∂ Φext ω kl j Φ 0 ext Þ 2 t 2 3 2 À γ À log ðω ir tÞ : We note that expressions similar to Eq. (44) have been previously reported in the literature 46,63 . However, these expressions do not include the correction ð 3 2 À γÞ within brackets in Eq. (44). Finally, we define the coherence 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 time due to flux noise.

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