Efficient modeling of superconducting quantum circuits with tensor networks

Abstract

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 15¹⁸⁰. 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 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 task rapidly becomes intractable by exact diagonalization. With current devices integrating 10s⁵ to 100s⁶, 1000s⁷ and even 10,000s⁸ 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⁵, 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 long-range 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⁹. 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 flux- and 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 multi-targeted 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⁹. 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)¹⁰. For a review, see for instance refs. ^{11,12,13,14,15}. The full many-body wavefunction can be written as¹²

$$\left|\psi \right\rangle =\mathop{\sum}\limits _{\{{\sigma }_{i}\}}{c}_{{\sigma }_{1}{\sigma }_{2}\ldots {\sigma }_{{N}_{{\rm{J}}}}}\left|{\sigma }_{1}{\sigma }_{2}\ldots {\sigma }_{{N}_{{\rm{J}}}}\right\rangle ,$$

(1)

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}_{{\sigma }_{1}{\sigma }_{2}\ldots {\sigma }_{{N}_{{\rm{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 $\left|\psi \right\rangle$, 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¹²

$$\begin{array}{lll}\left|\psi \right\rangle &=&\mathop{\sum}\limits _{\{{\sigma }_{i}\},\{{a}_{i}\}}{A}_{{a}_{1}}^{{\sigma }_{1}}{A}_{{a}_{1}{a}_{2}}^{{\sigma }_{2}}\ldots {A}_{{a}_{{N}_{{\rm{J}}}-2}{a}_{{N}_{{\rm{J}}}-1}}^{{\sigma }_{{N}_{{\rm{J}}}-1}}{A}_{{a}_{{N}_{{\rm{J}}}-1}}^{{\sigma }_{{N}_{{\rm{J}}}}}\\ &&\times \left|{\sigma }_{1}{\sigma }_{2}\ldots {\sigma }_{{N}_{{\rm{J}}}}\right\rangle ,\end{array}$$

(2)

where ${A}_{{a}_{i-1}{a}_{i}}^{{\sigma }_{i}}$ 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¹⁶. 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 ω₀₁ and anharmonicity ω₁₂ − ω₀₁, where $\hbar$ω_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}_{{\rm{q}}}^{2}$.

The conventional approach to compute excitations with DMRG is to add to the system Hamiltonian an energy penalty of the form ${\sum }_{i\in {\rm{ex.}}}{{\Lambda }}\left|{\psi }_{i}\right\rangle \left\langle {\psi }_{i}\right|$, with Λ > 0, where ex. denotes a set of previously determined excitations $\{\left|{\psi }_{i}\right\rangle \}$. This energy penalty turns the next excited into an effective ground state for which standard DMRG can be run¹². 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 long-range linear and nonlinear interactions, the fluxonium qubit⁵ 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²³.

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)² ≃ 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²⁹, while flux-noise insensitivity can be reached for large values of the superinductance³⁰. Recent realizations of the fluxonium qubit have shown long coherence times³¹, making these devices very attractive for quantum information processing.

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: $\phi ={\phi }_{0}-{\phi }_{{N}_{{\rm{J}}}}=\mathop{\sum }\nolimits_{i = 1}^{{N}_{{\rm{J}}}}{\theta }_{i}$.

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, such as NbTiN, TiN, and granular aluminum. Experiments on nanowire superinductors have recently demonstrated inductance values comparable to some of the junction-array counterparts, with low levels of loss and residual nonlinearity^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²⁷. 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}_{{{\rm{J}}}_{{\rm{b}}}}$ and its capacitance ${C}_{{{\rm{J}}}_{{\rm{b}}}}$, which may include a shunt capacitance. We take the superinductance to be realized by an array of Josephson junctions, with ${L}_{{{\rm{J}}}_{i}}$ and ${C}_{{{\rm{J}}}_{i}}$ being the ith junction inductance and capacitance, respectively. Moreover, a ground capacitance ${C}_{{0}_{i}}$ 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₀, respectively. We also define the plasma frequency ${\omega }_{{\rm{p}}}=1/\sqrt{{L}_{{\rm{J}}}{C}_{{\rm{J}}}}$ and reduced impedance $z=\sqrt{{L}_{{\rm{J}}}/{C}_{{\rm{J}}}}/{R}_{{\rm{Q}}}$ of the array junctions. Following the standard procedure for quantizing a circuit³, the Hamiltonian of the circuit takes the form (see section “Fluxonium circuit hamiltonian”)

$$H=\mathop{\sum }\limits_{i=1}^{{N}_{{\rm{J}}}}{H}_{{0}_{i}}+\mathop{\sum }\limits_{j>i}^{{N}_{{\rm{J}}}}\hslash {g}_{ij}\ {\tilde{n}}_{i}{\tilde{n}}_{j}-{E}_{{{\rm{J}}}_{{\rm{b}}}}\cos \left(\mathop{\sum }\limits_{i=1}^{{N}_{{\rm{J}}}}{\theta }_{i}+{\varphi }_{{\rm{ext}}}\right).$$

(3)

In this expression, ${H}_{{0}_{i}}=4{E}_{{{\rm{C}}}_{i}}{\tilde{n}}_{i}^{2}-{E}_{{{\rm{J}}}_{i}}\cos {\theta }_{i}$, with ${\tilde{n}}_{i}={n}_{i}-{n}_{{{\rm{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}_{{{\rm{g}}}_{i}}$ is an offset-charge parameter, ${E}_{{{\rm{C}}}_{i}}$ is the effective charging energy of this junction and ${E}_{{{\rm{J}}}_{i}}={\varphi }_{0}^{2}/{L}_{{{\rm{J}}}_{i}}$ is the Josephson energy with φ₀ = Φ₀/2π, where Φ₀ = h/2e is the flux quantum. In addition to the on-site energies, Eq. (3) includes a bilinear interaction $\propto {\tilde{n}}_{i}{\tilde{n}}_{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 = φ₀φ_ext associated with this junction³. 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, $\phi =\mathop{\sum }\nolimits_{i = 1}^{{N}_{{\rm{J}}}}{\theta }_{i}$, known as superinductance or fluxonium mode⁵.

To go beyond the usual effective model and obtain the low-energy 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⁴³ 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⁴⁴. 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.

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

$$H^{\prime} =4{E}_{{\rm{C}}}{n}^{^{\prime} 2}-{N}_{{\rm{J}}}^{2}{E}_{{\rm{L}}}\cos (\phi ^{\prime} /{N}_{{\rm{J}}})-{E}_{{\rm{J}}}\cos (\phi ^{\prime} +{\varphi }_{{\rm{ext}}}),$$

(4)

where the mode described by $\phi ^{\prime}$ is closely related to the superinductance (or fluxonium) mode ϕ, and $n^{\prime}$ 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}_{{0}_{i}}$ for i ∈ [1, N_J] can be neglected, then $\phi ^{\prime} =\phi$ and $n^{\prime} =n={N}_{{\rm{J}}}^{-1}\mathop{\sum }\nolimits_{i = 1}^{{N}_{{\rm{J}}}}{n}_{i}$, where n is the conjugate charge operator to ϕ. Otherwise, the $\phi ^{\prime}$ mode includes corrections to ϕ that are linear in C₀. 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 low-frequency 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)⁵, 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 single-mode nature, Eq. (4) can easily be diagonalized numerically by truncating the Hilbert space of the $\phi ^{\prime}$ 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. ²⁷). 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₁. 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 parameter regimes of the fluxonium qubit is offered in the Supplementary Discussion.

Each junction is modeled as a multilevel system using the 15 lowest energy eigenstates of the site Hamiltonian ${H}_{{0}_{i}}$. Indeed, we find that for low-impedance junctions, truncation errors can be avoided using a smaller number of eigenstates of ${H}_{{0}_{i}}$ compared to other possible basis, such as the Cooper-pair-number basis⁴⁶. The DMRG implementation is thus defined in a product basis of local wavefunctions spanning a many-body Hilbert space as large as 15¹²⁰ 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¹⁸⁰, 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.

Fig. 2: A 120-junction superinductance heavy fluxonium as a function of Φ_ext.

a Energy spectrum of the Hamiltonians in Eq. (3) (DMRG) and Eq. (4) (single mode). b Mean photon-number population of the array Josephson junctions (sites) for every eigenstate $\left|{\psi }_{k}\right\rangle$ of the fluxonium circuit. c Single-junction picture of fluxon- and plasmon-like excitations. d Schematic of the effective potential energy and wavefuctions of the single-mode Hamiltonian for Φ_ext ∈ {0, Φ₀/4, Φ₀/2}. e Expectation value of the phase operator at every circuit node of the superinductance for the fluxonium eigenstates labeled by $\left|{\psi }_{0}\right\rangle$ and $\left|{\psi }_{2}\right\rangle$. Circuit parameters: ${C}_{{{\rm{J}}}_{{\rm{b}}}}=40\ {\rm{fF}}$, ${E}_{{{\rm{J}}}_{{\rm{b}}}}/h=7.5\ {\rm{GHz}}$, C_J ≃ 32.9 fF, and L_J ≃ 1.23 nH (from ω_p/2π = 25 GHz and z = 0.03 (ref. ²⁷)) and C₀ = 0. Single-mode model parameters: E_C/h ≃ 0.48 GHz, E_L/h ≃ 1.27 GHz (i.e., L ≃ 129.1 nH), and ${E}_{{\rm{J}}}={E}_{{{\rm{J}}}_{{\rm{b}}}}$.

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. 2b shows the mean photon-number population $\langle {p}_{i}\rangle =\langle {\psi }_{k}| {H}_{{0}_{i}}| {\psi }_{k}\rangle /\hslash {\omega }_{{\rm{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 $\left|{\psi }_{k}\right\rangle$ 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 $\left|{\psi }_{0}\right\rangle \to \left|{\psi }_{i}\right\rangle$. 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 $\left|{\psi }_{k}\right\rangle$ 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 {\theta }_{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 single-mode effective Hamiltonian Eq. (4), as shown in Fig. 2d for Φ_ext ∈ {0, Φ₀/4, Φ₀/2}. The shape of the effective potential is determined by the cosine potential of the black-sheep junction and the inductive energy $-{N}_{{\rm{J}}}^{2}{E}_{{\rm{L}}}\cos (\phi ^{\prime} /{N}_{{\rm{J}}})\simeq {E}_{{\rm{L}}}\phi {^{\prime} }^{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 that $\langle {\psi }_{k}| \phi | {\psi }_{k}\rangle \equiv \langle {\psi }_{k}| \mathop{\sum }\nolimits_{i = 1}^{{N}_{{\rm{J}}}}{\theta }_{i}| {\psi }_{k}\rangle ={\alpha }_{k}$ for an excitation $\left|{\psi }_{k}\right\rangle$ 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. 2e, which shows the expectation value of the phase drop ${\phi }_{0}-{\phi }_{i}\equiv \mathop{\sum }\nolimits_{j = 1}^{i}{\theta }_{j}$, obtained from DMRG, and plotted as a function of the site number for the fluxon states $\left|{\psi }_{0}\right\rangle$ and $\left|{\psi }_{2}\right\rangle$ at Φ_ext = Φ₀/4 in Fig. 2d (middle panel). In this figure, the expectation value 〈ψ_k∣(ϕ₀ − ϕ_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 α₀ and α₂ 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 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 “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,48,49,50,51,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 black-sheep 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³⁰ 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²⁸.

Manucharyan et al.⁹ 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 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 Φ₀, where m is an integer-valued 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 ${H}_{{\rm{CQPS}}}=({E}_{{\rm{S}}}\ {m}^{-}+{E}_{{\rm{S}}}^{* }\ {m}^{+})/2$, where the operator m⁻ [${m}^{+}={({m}^{-})}^{\dagger }$] 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}_{{\rm{S}}}=\mathop{\sum }\nolimits_{i = 1}^{{N}_{{\rm{J}}}}{\epsilon }_{{0}_{i}}{e}^{i2\pi {n}_{{{\rm{g}}}_{i}}}$ (refs. ^28,9), where

$${\epsilon }_{{0}_{i}}=8\sqrt{2}\ \hslash {\omega }_{{{\rm{p}}}_{i}}\exp (-4/\pi {z}_{i})/\sqrt{\pi {z}_{i}},$$

(5)

corresponds to the charge dispersion of the ground state energy of the transmon Hamiltonian ${H}_{{0}_{i}}$ in terms of the reduced impedance z_i and plasma frequency ${\omega }_{{{\rm{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}_{{{\rm{g}}}_{i}}$ associated to that junction^{9,47,52,54,55,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 $\delta {\omega }_{ij}={\rm{Re}}[{E}_{{\rm{S}}}]({\langle T\rangle }_{j}-{\langle T\rangle }_{i})/\hslash$ of the qubit’s i → j transition frequency, where $T=\exp (-i2\pi n)$ is a 2π displacement operator whose expectation values are computed, using the unperturbed eigenstates $\{\left|{\psi }_{i}\right\rangle \}$ with m = 0 (ref. ⁹). For a homogeneous array (${\epsilon }_{{0}_{i}}\equiv {\epsilon }_{0}$ for i ∈ [1, N_J]), one has $-{N}_{{\rm{J}}}{\epsilon }_{0}\le {\rm{Re}}[{E}_{{\rm{S}}}]\le {N}_{{\rm{J}}}{\epsilon }_{0}$, and the total charge dispersion of the qubit transition frequency is

$$| {{\Delta }}{\omega }_{01}| =2{N}_{{\rm{J}}}{\epsilon }_{0}| {\langle T\rangle }_{1}-{\langle T\rangle }_{0}| /\hslash .$$

(6)

As the classical flux states of the loop are degenerate at Φ_ext = Φ₀/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. ⁹. Figure 3a shows the qubit transition frequency as a function of the external flux close to Φ_ext = Φ₀/2 for different values of the offset-charge ${n}_{{{\rm{g}}}_{i}}\equiv {n}_{{\rm{g}}}$, assumed to be the same on every junction of the array. Each subpanel 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}_{{0}_{i}}$, 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. ⁵⁷). 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.

Fig. 3: Charge dispersion of a 40-junction superinductance fluxonium qubit as a function of the reduced impedance of the array junctions.

a Broadening of the fluxon transition around Φ_ext = Φ₀/2 for n_g ∈ [0, 0.5]. Color lines are obtained using the multi-targeted DMRG algorithm, while dashed black lines correspond to estimations using the single-mode Hamiltonian Eq. (4). b Total charge dispersion of the fluxon transition at Φ_ext = Φ₀/2 according to the DMRG calculation (circles) contrasted to the prediction of Eq. (6) with matrix elements evaluated by means of DMRG (triangles) or the single-mode model (dashed lines). Inset: data displayed in linear scale. Parameters: ${C}_{{{\rm{J}}}_{{\rm{b}}}}=7.5\ {\rm{fF}}$, ${E}_{{{\rm{J}}}_{{\rm{b}}}}/h=8.9\ {\rm{GHz}}$, ω_p/2π = 12.5, and C₀ = 0, according to ref. ⁹.

Qualitatively, charge dispersion increases rapidly with the array-junction impedance z due to the exponential scaling of Eq. (5). This is best illustrated by Fig. 3b, which shows the charge dispersion for Φ_ext = Φ₀/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. ²⁸. Although barely visible, 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⁻¹¹, and the error tolerance on the eigenvalues are set to 10⁻¹², 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. ⁹ 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. ⁹, 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. ⁹ assumes that the variables ${n}_{{{\rm{g}}}_{i}}$ are independent and randomly distributed. The probability distribution of ${\rm{Re}}[{E}_{{\rm{S}}}]$ can then be approximated by a Gaussian with zero mean and standard deviation $\sqrt{{N}_{{\rm{J}}}/2}\ {\epsilon }_{0}$ (ref. ⁹). The effective broadening of the qubit transition in presence of noise thus scales as $\sqrt{{N}_{{\rm{J}}}}$, translating to the pure-dephasing rate $1/{T}_{\varphi ,{\rm{CQPS}}}=| {{\Delta }}{\omega }_{01}| /4\sqrt{{N}_{{\rm{J}}}}$ (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. ⁹ predicts the pure-dephasing coherence times to be dominated by the combined effect of charge noise and CQPS around Φ_ext = Φ₀/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”)

$$\begin{array}{lll}{\partial }_{t}\rho \,=\,\mathop{\sum}\limits _{k}{{{\Gamma }}}_{\varphi }^{kk}\ {\mathcal{D}}[{\sigma }_{kk},{\sigma }_{kk}]\ \rho \\ \qquad+\,\mathop{\sum}\limits _{k > l}{{{\Gamma }}}_{\varphi }^{kl}\ \left({\mathcal{D}}[{\sigma }_{kk},{\sigma }_{ll}]+{\mathcal{D}}[{\sigma }_{ll},{\sigma }_{kk}]\right)\rho ,\end{array}$$

(7)

where ${{{\Gamma }}}_{\varphi }^{kl}$ are time-dependent pure-dephasing rates proportional to the 1/f flux-noise amplitude, ${\sigma }_{kl}=\left|{\psi }_{k}\right\rangle \left\langle {\psi }_{l}\right|$, and ${\mathcal{D}}[x,y]\ \rho =x\rho {y}^{\dagger }-\{{y}^{\dagger }x,\rho \}/2$ is a generalized dissipator operator. By integrating Eq. (7), we find the flux-noise coherence time T_φ,Flux by solving the implicit equation ρ₀₁(T_φ,Flux)/ρ₀₁(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. ⁹. In contrast to the results in Fig. 2a, the difference between the DMRG and single-mode simulations for the parameters of ref. ⁹ 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 = Φ₀/2. We find coherence time values that are very similar to the experimental observation (see Fig. 4 in ref. ⁹), 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. ⁹.

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}_{{{\rm{J}}}_{{\rm{b}}}}=7.5\ {\rm{fF}}$, ${E}_{{{\rm{J}}}_{{\rm{b}}}}/h=8.9\ {\rm{GHz}}$, z = 0.09, ω_p/2π = 12.5, and C₀ = 0, extracted from ref. ⁹. The 1/f flux-noise amplitude is taken to be A_Φ = 1.2 μΦ₀, which is a conservative value⁴⁶.

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. ⁹, 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.

Methods

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}_{{{\rm{J}}}_{{\rm{b}}}}$ (including both shunt and junction capacitances) and Josephson energy ${E}_{{{\rm{J}}}_{{\rm{b}}}}$ is shunted by a superinductance made of N_J junctions, each of capacitance ${C}_{{{\rm{J}}}_{i}}$ and energy ${E}_{{{\rm{J}}}_{i}}$ with i ∈ [1, N_J]. We moreover assume that each circuit node of the superinductance is connected to ground by a stray capacitance ${C}_{{0}_{i}}$. The N_J + 1 node flux (phase) variables of the circuit are denoted by Φ_i (ϕ_i = Φ_i/φ₀), where φ₀ = ℏ/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³

$$\begin{array}{lll}L({\boldsymbol{\Phi }},{\dot{\boldsymbol{{{\Phi }}}}})&=&\frac{{C}_{{{\rm{J}}}_{{\rm{b}}}}}{2}{({\dot{{{\Phi }}}}_{{N}_{{\rm{J}}}}-{\dot{{{\Phi }}}}_{0})}^{2}+\mathop{\sum }\limits_{i=1}^{{N}_{{\rm{J}}}}\frac{{{C}_{{\rm{J}}}}_{i}}{2}{({\dot{{{\Phi }}}}_{i}-{\dot{{{\Phi }}}}_{i-1})}^{2}\\ &+&\mathop{\sum }\limits_{i=0}^{{N}_{{\rm{J}}}}\frac{{C}_{{0}_{i}}}{2}{\dot{{{\Phi }}}}_{i}^{2}+\mathop{\sum }\limits_{i=1}^{{N}_{{\rm{J}}}}{E}_{{{\rm{J}}}_{i}}\cos \left[({{{\Phi }}}_{i}-{{{\Phi }}}_{i-1})/{\varphi }_{0}\right]\\ &+&{E}_{{{\rm{J}}}_{{\rm{b}}}}\cos \left[({{{\Phi }}}_{{N}_{{\rm{J}}}}-{{{\Phi }}}_{0}+{{{\Phi }}}_{{\rm{ext}}})/{\varphi }_{0}\right],\end{array}$$

(8)

where Φ_ext is the flux through the circuit loop. A more convenient basis is defined by the flux variables Θ_i = Φ_i−1 − Φ_i for i ∈ [1, N_J] and the cyclic mode ${{\Sigma }}=\mathop{\sum }\nolimits_{i = 0}^{{N}_{{\rm{J}}}}{{{\Phi }}}_{i}$. The relation between these modes and the original node fluxes can be expressed concisely by Θ = R ⋅ Φ, where ${\boldsymbol{\Theta }}={({{{\Theta }}}_{1},\ldots ,{{{\Theta }}}_{{N}_{{\rm{J}}}},{{\Sigma }})}^{T}$, ${\boldsymbol{\Phi }}={({{{\Phi }}}_{0},\ldots ,{{{\Phi }}}_{{N}_{{\rm{J}}}})}^{T}$, and R is the N_J + 1 × N_J + 1 matrix

$${\bf{R}}=\left[\begin{array}{lllllll}1&-1&0&\cdots \ &\cdots \ &\cdots \ &0\\ 0&1&-1&0&\cdots \ &\cdots \ &0\\ 0&0&1&-1&0&\cdots \ &0\\ \vdots &\vdots &\ddots &\ddots &\ddots &\ddots &\vdots \\ 0&\cdots \ &\cdots \ &0&1&-1&0\\ 0&\cdots \ &\cdots \ &\cdots \ &0&1&-1\\ 1&1&1&\cdots \ &\cdots \ &1&1\end{array}\right].$$

(9)

Under this change of basis, Eq. (8) becomes

$$\begin{array}{lll}L({\boldsymbol{\Theta }},{\dot{\boldsymbol{{{\Theta }}}}})&=&{\dot{{{\Theta }}}}^{T}\cdot \frac{{{\boldsymbol{C}}}_{{{\Theta }}}}{2}\cdot \dot{{{\Theta }}}+\mathop{\sum }\limits_{i=1}^{{N}_{{\rm{J}}}}{E}_{{{\rm{J}}}_{i}}\cos ({{{\Theta }}}_{i}/{\varphi }_{0})\\ &+&{E}_{{{\rm{J}}}_{{\rm{b}}}}\cos \left[\left(\mathop{\sum }\limits_{i=1}^{{N}_{{\rm{J}}}}{{{\Theta }}}_{i}+{{{\Phi }}}_{{\rm{ext}}}\right)/{\varphi }_{0}\right],\end{array}$$

(10)

where ${{\boldsymbol{C}}}_{{{\Theta }}}={({{\bf{R}}}^{-1})}^{T}\cdot {{\boldsymbol{C}}}_{{{\Phi }}}\cdot {{\bf{{{R}}}}}^{-1}$ is defined in terms of the capacitance matrix ${[{{\boldsymbol{C}}}_{{{\Phi }}}]}_{ij}={\partial }^{2}L({\boldsymbol{\Phi }},{\dot{\boldsymbol{{{\Phi }}}}})/\partial {\dot{{{\Phi }}}}_{i}\partial {\dot{{{\Phi }}}}_{j}$, for i, j ∈ [0, N_J + 1]. Note that the Σ mode does not enter in the potential energy.

After a Legendre transformation, we arrive at the circuit Hamiltonian

$$\begin{array}{lll}H\,=\,{{\boldsymbol{q}}}_{{{\Theta }}}^{T}\cdot \frac{{{\boldsymbol{C}}}_{{{\Theta }}}^{-1}}{2}\cdot {{\boldsymbol{q}}}_{{{\Theta }}}-\mathop{\sum }\limits_{i=1}^{{N}_{{\rm{J}}}}{E}_{{{\rm{J}}}_{i}}\cos {\theta }_{i}\\ \quad\,-\,{E}_{{{\rm{J}}}_{{\rm{b}}}}\cos \left(\mathop{\sum }\limits_{i=1}^{{N}_{{\rm{J}}}}{\theta }_{i}+{\varphi }_{{\rm{ext}}}\right),\end{array}$$

(11)

where ${{\boldsymbol{q}}}_{{{\Theta }}}\equiv \partial L({\boldsymbol{\Theta }},{\dot{\boldsymbol{{{\Theta }}}}})/\partial {\dot{\boldsymbol{{{\Theta }}}}}={{\boldsymbol{C}}}_{{{\Theta }}}\cdot {\dot{\boldsymbol{{{\Theta }}}}}$ is a vector of conjugate charge operators, θ_i = Θ_i/φ₀ are phase operators and φ_ext = Φ_ext/φ₀. In the presence of disorder in the circuit capacitances, the σ = Σ/φ₀ 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 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 $\mathop{\sum }\nolimits_{i = 0}^{{N}_{{\rm{J}}}}({C}_{{{\rm{g}}}_{i}}/2){({\dot{{{\Phi }}}}_{i}-{V}_{i})}^{2}$, where ${C}_{{{\rm{g}}}_{i}}$ is a gate capacitance for the ith circuit node. Equivalently, this can be expressed as

$${L}_{g}({\boldsymbol{\Phi }},\dot{{\boldsymbol{\Phi }}})=-{\dot{{\boldsymbol{\Phi }}}}^{T}\cdot {{\boldsymbol{C}}}_{{\rm{g}}}\cdot {\boldsymbol{V}},$$

(12)

where ${{\boldsymbol{C}}}_{{\rm{g}}}={\rm{diag}}({C}_{{{\rm{g}}}_{0}},{C}_{{{\rm{g}}}_{1}},\ldots ,{C}_{{{\rm{g}}}_{{N}_{{\rm{J}}}+1}})$ and ${\boldsymbol{V}}={({V}_{0},{V}_{1},\ldots ,{V}_{{N}_{{\rm{J}}}+1})}^{T}$. In addition to Eq. (12), the capacitance matrix of the circuit is modified to account for the gate capacitances as ${{\boldsymbol{C}}}_{{{\Phi }}}\to {\widetilde{{\boldsymbol{C}}}}_{{{\Phi }}}={{\boldsymbol{C}}}_{{{\Phi }}}+{{\boldsymbol{C}}}_{{\rm{g}}}$.

Defining d_Φ = C_g ⋅ V, the conjugate charge operators are given by

$${{\boldsymbol{q}}}_{{{\Theta }}}={\widetilde{{\boldsymbol{C}}}}_{{{\Theta }}}\cdot \dot{{\boldsymbol{\Theta }}}-{{\boldsymbol{d}}}_{{{\Theta }}},$$

(13)

where ${\widetilde{{\boldsymbol{C}}}}_{{{\Theta }}}={({{\bf{{R}}}}^{-1})}^{T}\cdot {\widetilde{{\boldsymbol{C}}}}_{{{\Phi }}}\cdot {{\bf{{R}}}}^{-1}$ and ${{\boldsymbol{d}}}_{{{\Theta }}}={({{\bf{{R}}}}^{-1})}^{T}\cdot {{\boldsymbol{d}}}_{{{\Phi }}}$. The circuit Hamiltonian then takes the form

$$\begin{array}{lll}H&=&{({{\boldsymbol{q}}}_{{{\Theta }}}+{{\boldsymbol{d}}}_{{{\Theta }}})}^{T}\cdot \frac{{\widetilde{{\boldsymbol{C}}}}_{{{\Theta }}}^{-1}}{2}\cdot ({{\boldsymbol{q}}}_{{{\Theta }}}+{{\boldsymbol{d}}}_{{{\Theta }}})\\ &-&\mathop{\sum }\limits_{i=1}^{{N}_{{\rm{J}}}}{E}_{{{\rm{J}}}_{i}}\cos {\theta }_{i}-{E}_{{{\rm{J}}}_{{\rm{b}}}}\cos \left(\mathop{\sum }\limits_{i=1}^{{N}_{{\rm{J}}}}{\theta }_{i}+\frac{{{{\Phi }}}_{{\rm{ext}}}}{{\varphi }_{0}}\right).\end{array}$$

(14)

Omitting the σ mode and irrelevant constants, the above expression simplifies to

$$\begin{array}{lll}H&=&\mathop{\sum }\limits_{i=1}^{{N}_{{\rm{J}}}}\left[\frac{{[{\widetilde{{\boldsymbol{C}}}}_{{{\Theta }}}^{-1}]}_{ii}}{2}{({q}_{i}-{q}_{{{\rm{g}}}_{i}})}^{2}-{E}_{{{\rm{J}}}_{i}}\cos {\theta }_{i}\right]\\ &+&\mathop{\sum }\limits_{j>i}^{{N}_{{\rm{J}}}}{[{\widetilde{{\boldsymbol{C}}}}_{{{\Theta }}}^{-1}]}_{ij}({q}_{i}-{q}_{{{\rm{g}}}_{i}})({q}_{j}-{q}_{{g}_{j}})\\ &-&{E}_{{{\rm{J}}}_{{\rm{b}}}}\cos \left(\mathop{\sum }\limits_{i=1}^{{N}_{{\rm{J}}}}{\theta }_{i}+{\varphi }_{{\rm{ext}}}\right),\end{array}$$

(15)

where ${q}_{{{\rm{g}}}_{i}}={[{\widetilde{{\boldsymbol{C}}}}_{{{\Theta }}}^{-1}\cdot {{\boldsymbol{d}}}_{{{\Theta }}}]}_{i}/2{[{\widetilde{{\boldsymbol{C}}}}_{{{\Theta }}}^{-1}]}_{ii}$ for i ∈ [1, N_J] are effective offset charges in the θ_i basis and ${q}_{i}={[{q}_{{{\Theta }}}]}_{i}$. This Hamiltonian is equivalent to Eq. (3). Each of the bracketed terms in Eq. (15) define a site Hamiltonian (${H}_{{0}_{i}}$ for the ith array junction), while the remaining terms lead to linear and nonlinear all-to-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}_{{{\rm{g}}}_{i}}={q}_{{{\rm{g}}}_{i}}/2e$, instead of q_i and ${q}_{{{\rm{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 $\phi =\mathop{\sum }\nolimits_{i = 1}^{{N}_{{\rm{J}}}}{\theta }_{i}$ is simple. To find this appropriate change of coordinates, we reverse engineer the following Ansatz defining an additional change of basis

$${{\bf{{R}}}}^{(1)}=\left[\begin{array}{lllll}1-\mathop{\sum }\nolimits_{k = 1}^{{N}_{{\rm{J}}}-1}{a}_{k}^{(1)}&1+{a}_{1}^{(1)}&\cdots \ &1+{a}_{{N}_{{\rm{J}}}-1}^{(1)}&0\\ -1&1&0&\cdots \ &0\\ \vdots &0&\ddots &\ddots &\vdots \\ -1&\vdots &\ddots &1&0\\ 0&0&0&0&1\end{array}\right],$$

(16)

where the constants $\{{a}_{k}^{(1)}\}$ are defined by

$${a}_{k}^{(1)}=\frac{\mathop{\sum }\nolimits_{i,j = 0}^{{N}_{{\rm{J}}}-1}({N}_{{\rm{J}}}{[{{\boldsymbol{C}}}_{{{\Theta }}}]}_{ik}{\delta }_{jk}-{[{{\boldsymbol{C}}}_{{{\Theta }}}]}_{ij})}{\mathop{\sum }\nolimits_{i,j = 0}^{{N}_{{\rm{J}}}-1}{[{{\boldsymbol{C}}}_{{{\Theta }}}]}_{ij}},$$

(17)

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⁽¹⁾ is to capacitively decouple a superinductance-like mode of the form

$${\phi }^{(1)}=\phi +\mathop{\sum }\limits_{k=1}^{{N}_{{\rm{J}}}-1}{a}_{k}^{(1)}({\theta }_{k}-{\theta }_{1}),$$

(18)

from all other circuit modes, while leaving the σ mode invariant. Indeed, the capacitance matrix now reads

$${{\boldsymbol{C}}}_{X}^{(1)}={[{({{\bf{{R}}}}^{(1)})}^{-1}]}^{T}\cdot {{\boldsymbol{C}}}_{X}^{(0)}\cdot {({{\bf{{R}}}}^{(1)})}^{-1},$$

(19)

where ${{\boldsymbol{C}}}_{X}^{(0)}={{\boldsymbol{C}}}_{{{\Theta }}}$ is block diagonal. Importantly, this transformation works in the presence of circuit-element disorder, in which case the block-diagonal form of ${{\boldsymbol{C}}}_{X}^{(1)}$ 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 ϕ⁽¹⁾ mode; the second block has dimension (N_J − 1) × (N_J − 1) and involves all circuit modes except ϕ⁽¹⁾ 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⁽ⁿ⁾}, for n ∈ [2, N_J − 1], with the general form

$${{\bf{{R}}}}^{(n)}=\left[\begin{array}{llllllllll}1&0&0&\cdots \ &\cdots \ &\cdots \ &\cdots \ &\cdots \ &\cdots \ &0\\ 0&1&0&\cdots \ &\cdots \ &\cdots \ &\cdots \ &\cdots \ &\cdots \ &0\\ \vdots &\ddots &\ddots &\ddots &\cdots \ &\cdots \ &\cdots \ &\cdots \ &\cdots \ &0\\ \vdots &\vdots &0&1&0&\cdots \ &\cdots \ &\cdots \ &\cdots \ &0\\ \vdots &\vdots &\vdots &0&1-\mathop{\sum }\nolimits_{k = n}^{{N}_{{\rm{J}}}-1}{a}_{k}^{(n)}&1+{a}_{n}^{(n)}&1+{a}_{n+1}^{(n)}&\cdots \ &1+{a}_{{N}_{{\rm{J}}}-1}^{(n)}&0\\ \vdots &\vdots &\vdots &\vdots &-1&1&0&\cdots \ &0&\vdots \\ \vdots &\vdots &\vdots &\vdots &-1&0&1&0&0&\vdots \\ \vdots &\vdots &\vdots &\vdots &\vdots &\vdots &\ddots &\ddots &\ddots &\vdots \\ \vdots &\vdots &\vdots &\vdots &-1&0&\cdots \ &0&1&0\\ 0&0&0&0&0&\cdots \ &\cdots \ &\cdots \ &0&1\end{array}\right].$$

(20)

Similarly to R⁽¹⁾, the matrix R⁽ⁿ⁾ 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}_{k}^{(n)}\}$ are defined as

$${a}_{k}^{(n)}=\frac{\mathop{\sum }\nolimits_{i,j = n}^{{N}_{{\rm{J}}}-1}\{({N}_{{\rm{J}}}-1+n){[{{\boldsymbol{C}}}_{X}^{(n-1)}]}_{ik}{\delta }_{jk}-{[{{\boldsymbol{C}}}_{X}^{(n-1)}]}_{ij}\}}{\mathop{\sum }\nolimits_{i,j = n}^{{N}_{{\rm{J}}}-1}{[{{\boldsymbol{C}}}_{X}^{(n-1)}]}_{ij}},$$

(21)

which is a generalization of Eq. (17).

The transformations ${{\bf{{R}}}}^{(n\,{<}\,{N}_{{\rm{J}}}-1)}$ are designed to each decouple a single mode, while ${{\bf{{R}}}}^{({N}_{{\rm{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 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 $\{{a}_{k}^{(n)}\}$ 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 so-called 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 single-mode 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

$${\theta }_{i}=\frac{{\phi }^{(1)}}{{N}_{{\rm{J}}}}+\mathop{\sum }\limits_{n=2}^{{N}_{{\rm{J}}}}{v}_{ni}{\phi }^{(n)},$$

(22)

where the coefficient v_ni quantifies how much the ϕ⁽ⁿ⁾ 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 {ϕ⁽ⁿ⁾} 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 $\phi =\mathop{\sum }\nolimits_{i = 1}^{{N}_{{\rm{J}}}}{\theta }_{i}$, we moreover have

$$\phi ={\phi }^{(1)}+\mathop{\sum }\limits_{n=2}^{{N}_{{\rm{J}}}}{{\mathcal{V}}}_{n}{\phi }^{(n)},$$

(23)

where ${{\mathcal{V}}}_{n}=\mathop{\sum }\nolimits_{i = 1}^{{N}_{{\rm{J}}}}{v}_{ni}$. If C₀ = 0, it follows that ${{\mathcal{V}}}_{n}=0$ for n ∈ [2, N_J], and ϕ⁽¹⁾ ≡ ϕ 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 ϕ⁽ⁿ⁾ for n > 1 in terms of the harmonic oscillator ladder operators as ${\phi }^{(n)}=\sqrt{\pi {z}_{n}}({a}_{n}+{a}_{n}^{\dagger })$. Here, ${z}_{n}=\sqrt{{L}_{n}/{C}_{n}}/{R}_{{\rm{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}_{n}^{-1}={{\boldsymbol{X}}}_{n}^{T}\cdot {({{\boldsymbol{M}}}^{-1})}^{T}\cdot {{\boldsymbol{L}}}^{-1}\cdot {{\boldsymbol{M}}}^{-1}\cdot {{\boldsymbol{X}}}_{n}$, where Xⁿ is the mode vector associated to ϕ⁽ⁿ⁾ and ${\boldsymbol{M}}={(\mathop{\prod }\nolimits_{n = 1}^{{N}_{{\rm{J}}}-1}{{\bf{R}}}^{(n)})}^{T}\cdot {\bf{R}}$ is a matrix that reverses the multiple changes of basis. The trace can then be performed straightforwardly by noticing that

$${e}^{ix{\phi }^{(n)}}={e}^{-\pi {x}^{2}{z}_{n}/2}{e}^{ix\sqrt{\pi {z}_{n}}{a}_{n}^{\dagger }}{e}^{ix\sqrt{\pi {z}_{n}}{a}_{n}},$$

(24)

and thus ${{\rm{tr}}}_{n}[{e}^{ix{\phi }^{(n)}}\rho ]={e}^{-\pi {x}^{2}{z}_{n}/2}$ where we assume that the nth mode remains in its noninteracting vacuum state. Following to Eqs. (22) and (23), we approximate

$$\begin{array}{lll}\cos {\theta }_{i}&\simeq &{{\rm{tr}}}_{n \,{>}\,1}[\cos {\theta }_{i}]\\ &\simeq &{x}_{i}\cos [{\phi }^{(1)}/{N}_{{\rm{J}}}],\end{array}$$

(25)

where ${x}_{i}=\mathop{\prod }\nolimits_{n = 2}^{{N}_{{\rm{J}}}}{e}^{-\pi {v}_{ni}^{2}{z}_{n}/2}$, and

$$\begin{array}{lll}\cos (\phi +{\varphi }_{{\rm{ext}}})&\simeq &{{\rm{tr}}}_{n \,{>}\,1}[\cos (\phi +{\varphi }_{{\rm{ext}}})]\\ &\simeq &{x}_{{\rm{b}}}\cos [{\phi }^{(1)}+{\varphi }_{{\rm{ext}}}],\end{array}$$

(26)

with ${x}_{{\rm{b}}}=\mathop{\prod }\nolimits_{n = 2}^{{N}_{{\rm{J}}}}{e}^{-\pi {{\mathcal{V}}}_{n}^{2}{z}_{n}/2}$. In Eqs. (25) and (26), ${{\rm{tr}}}_{n \,{>}\,1}$ indicates a trace operation over all circuit modes ϕ⁽ⁿ⁾, 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 ${\phi }^{(1)}\to \phi ^{\prime}$, we arrive at the effective single-mode Hamiltonian

$$\begin{array}{lll}H&=&4{E}_{{\rm{C}}}{n}^{^{\prime} 2}-\mathop{\sum }\limits_{i=1}^{{N}_{{\rm{J}}}}{x}_{i}{E}_{{{\rm{J}}}_{i}}\cos (\phi ^{\prime} /{N}_{{\rm{J}}})\\ &-&{x}_{{\rm{b}}}{E}_{{{\rm{J}}}_{{\rm{b}}}}\cos (\phi ^{\prime} +{\varphi }_{{\rm{ext}}}),\end{array}$$

(27)

where E_C is taken to be the charging energy ${E}_{{\rm{C}}}={e}^{2}/2{[{{\boldsymbol{C}}}_{X}^{(1)}]}_{00}$ of the $\phi ^{\prime}$ mode and $[\phi ^{\prime} ,n^{\prime} ]=i$. Note that Eq. (27) is equivalent to Eq. (4) of the main text. Up to corrections of order ${N}_{{\rm{J}}}^{-3}$, Eq. (27) reduces to

$$H=4{E}_{{\rm{C}}}{n}^{^{\prime} 2}+\frac{{E}_{{\rm{L}}}}{2}{\phi }^{^{\prime} 2}-{E}_{{\rm{J}}}\cos (\phi ^{\prime} +{\varphi }_{{\rm{ext}}}),$$

(28)

where ${E}_{{\rm{L}}}=\mathop{\sum }\nolimits_{i = 1}^{{N}_{{\rm{J}}}}{x}_{i}{E}_{{{\rm{J}}}_{i}}/{N}_{{\rm{J}}}^{2}$ and ${E}_{{\rm{J}}}={x}_{{\rm{b}}}{E}_{{{\rm{J}}}_{{\rm{b}}}}$ are the effective inductive and Josephson junction energies. Equation (28) corresponds to the original fluxonium-qubit model of ref. ⁵. 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

$${\partial }_{t}\rho (t)=-\frac{1}{{\hslash }^{2}}{\int\nolimits_{0}^{t}}d\tau \ {{\rm{tr}}}_{B}[{H}_{{\rm{int}}}(t),[{H}_{{\rm{int}}}(t-\tau ),\rho (t-\tau )\otimes {\rho }_{{\rm{B}}}]],$$

(29)

where ρ(t) ⊗ ρ_B is the system–bath density matrix, assumed to be separable at all times⁵⁸. 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 ${{{\Phi }}}_{{\rm{ext}}}={{{\Phi }}}_{{\rm{ext}}}^{0}+\delta {{\Phi }}$, where ${{{\Phi }}}_{{\rm{ext}}}^{0}$ is the applied flux bias and δΦ represents fluctuations. To first order in δΦ, the interaction Hamiltonian can be written as⁵⁹

$${H}_{{\rm{int}}}={\partial }_{{{{\Phi }}}_{{\rm{ext}}}}H{| }_{{{{\Phi }}}_{{\rm{ext}}}^{0}}\times \delta {{\Phi }},$$

(30)

where H is the Hamiltonian of the fluxonium qubit and the derivative with respect to the external flux is evaluated at ${{{\Phi }}}_{{\rm{ext}}}={{{\Phi }}}_{{\rm{ext}}}^{0}$. Expanding Eq. (29) in the eigenbasis $\{\left|{\psi }_{k}\right\rangle \}$ of the full circuit, we arrive at

$$\begin{array}{lll}{\partial }_{t}\rho &=&-\frac{1}{{\hslash }^{2}}\sum _{\begin{array}{c}k,k^{\prime} \\ l,l^{\prime} \end{array}}{\int\nolimits_{0}^{t}}d\tau \ {\partial }_{{{{\Phi }}}_{{\rm{ext}}}}H{| }_{{{{\Phi }}}_{{\rm{ext}}}^{0}}^{kk^{\prime} }{\partial }_{{{{\Phi }}}_{{\rm{ext}}}}H{| }_{{{{\Phi }}}_{{\rm{ext}}}^{0}}^{ll^{\prime} }{e}^{-i({\omega }_{ll^{\prime} }+{\omega }_{kk^{\prime} })t+i{\omega }_{kk^{\prime} }\tau }\\ &&\times {{\rm{tr}}}_{B}[\left|{\psi }_{l}\right\rangle \left\langle {\psi }_{l^{\prime} }\right|\delta {{\Phi }}(t),[\left|{\psi }_{k}\right\rangle \left\langle {\psi }_{k^{\prime} }\right|\delta {{\Phi }}(t-\tau ),\rho \otimes {\rho }_{{\rm{B}}}]],\end{array}$$

(31)

where we have introduced the matrix elements ${\partial }_{{{{\Phi }}}_{{\rm{ext}}}}H{| }_{{{{\Phi }}}_{{\rm{ext}}}^{0}}^{kk^{\prime} }=\langle {\psi }_{k}| {\partial }_{{{{\Phi }}}_{{\rm{ext}}}}H{| }_{{{{\Phi }}}_{{\rm{ext}}}^{0}}| {\psi }_{k^{\prime} }\rangle$, and omitted the explicit time dependence of ρ(t) → ρ.

Tracing out the bath degrees of freedom leads to the so-called Bloch–Redfield equation⁵⁸. 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 ${\omega }_{ll^{\prime} }+{\omega }_{kk^{\prime} }\ne 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^{\prime}$ and $l^{\prime} =k$ for ${\omega }_{kk^{\prime} }\ne 0$, and $l=l^{\prime}$ for ${\omega }_{kk^{\prime} }=0$. In this way, Eq. (31) simplifies to

$$\begin{array}{lll}{\partial }_{t}\rho &=&-\frac{1}{{\hslash}^{2}}\mathop{\sum}\limits _{k^{\prime} > k}{\int\nolimits_{0}^{\infty}}d\tau \ {\partial }_{{{{\Phi }}}_{{\rm{ext}}}}H{| }_{{{{\Phi}}}_{{\rm{ext}}}^{0}}^{kk^{\prime} }{\partial }_{{{{\Phi}}}_{{\rm{ext}}}}H{| }_{{{{\Phi }}}_{{\rm{ext}}}^{0}}^{k^{\prime}k}{e}^{i{\omega }_{kk^{\prime} }\tau }\,{{\rm{tr}}}_{B}[| {\psi}_{k^{\prime} }\rangle \langle {\psi }_{k}| \delta {{\Phi}}(t),[\left|{\psi }_{k}\right\rangle \left\langle {\psi}_{k^{\prime} }\right|\delta {{\Phi }}(t-\tau ),\rho \otimes {\rho}_{{\rm{B}}}]]\\ &&-\frac{1}{{\hslash }^{2}}\mathop{\sum}\limits _{k^{\prime}> k}{\int\nolimits_{0}^{\infty }}d\tau \ {\partial }_{{{{\Phi}}}_{{\rm{ext}}}}H{| }_{{{{\Phi }}}_{{\rm{ext}}}^{0}}^{k^{\prime}k}{\partial }_{{{{\Phi }}}_{{\rm{ext}}}}H{| }_{{{{\Phi}}}_{{\rm{ext}}}^{0}}^{kk^{\prime} }{e}^{-i{\omega }_{kk^{\prime}}\tau }\,{{\rm{tr}}}_{B}[\left|{\psi }_{k}\right\rangle \left\langle{\psi }_{k^{\prime} }\right|\delta {{\Phi }}(t),[\left|{\psi}_{k^{\prime} }\right\rangle \left\langle {\psi }_{k}\right|\delta{{\Phi }}(t-\tau ),\rho \otimes {\rho }_{{\rm{B}}}]]\\&&-\frac{1}{{\hslash }^{2}}\mathop{\sum}\limits_{k,l}{\int\nolimits_{0}^{\infty }}d\tau \ {\partial}_{{{{\Phi }}}_{{\rm{ext}}}}H{| }_{{{{\Phi}}}_{{\rm{ext}}}^{0}}^{kk}{\partial }_{{{{\Phi }}}_{{\rm{ext}}}}H{|}_{{{{\Phi }}}_{{\rm{ext}}}^{0}}^{ll}\,{{\rm{tr}}}_{B}[\left|{\psi}_{l}\right\rangle \left\langle {\psi }_{l}\right|\delta {{\Phi}}(t),[\left|{\psi }_{k}\right\rangle \left\langle {\psi}_{k}\right|\delta {{\Phi }}(t-\tau ),\rho \otimes {\rho}_{{\rm{B}}}]].\end{array}$$

(32)

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 (refs. ^60,61). The pure-dephasing master equation is then derived from the third line of Eq. (32), i.e.,

$$\begin{array}{lll}{\partial }_{t}\rho &=&-\frac{1}{{\hslash }^{2}}\mathop{\sum}\limits _{k,l}{\int\nolimits_{0}^{\infty }}d\tau \ {\partial }_{{{{\Phi }}}_{{\rm{ext}}}}H{| }_{{{{\Phi }}}_{{\rm{ext}}}^{0}}^{kk}{\partial }_{{{{\Phi }}}_{{\rm{ext}}}}H{| }_{{{{\Phi }}}_{{\rm{ext}}}^{0}}^{ll}\\ &\times &{{\rm{tr}}}_{B}[\left|{\psi }_{l}\right\rangle \left\langle {\psi }_{l}\right|\delta {{\Phi }}(t),[\left|{\psi }_{k}\right\rangle \left\langle {\psi }_{k}\right|\delta {{\Phi }}(t-\tau ),\rho \otimes {\rho }_{{\rm{B}}}]].\end{array}$$

(33)

To obtain a workable expression, we introduce the noise spectral density ${S}_{{{\Phi }}}^{1/f}[\omega ]$ for 1/f flux noise by the definition³

$${{\rm{tr}}}_{B}[{\rho }_{{\rm{B}}}\delta {{\Phi }}(t)\delta {{\Phi }}(t^{\prime} )]=\frac{1}{2\pi }\mathop{\int}\nolimits_{-\infty }^{\infty }d\omega \ {S}_{{{\Phi }}}^{1/f}[\omega ]{e}^{-i\omega (t-t^{\prime} )},$$

(34)

and assume the general form

$${S}_{{{\Phi }}}^{1/f}(\omega )=\frac{{A}_{{{\Phi }}}^{2}}{| \omega | /2\pi },$$

(35)

where A_Φ is the 1/f flux-noise amplitude, typically reported to be in the range 1–10 μΦ₀ (ref. ⁴⁶). 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

$$\mathop{\int}\nolimits_{0}^{t}d\tau \ {{\rm{tr}}}_{B}[{\rho }_{{\rm{B}}}\delta {{\Phi }}(t)\delta {{\Phi }}(t^{\prime} )]={\mathrm{lim}\,}_{{\omega }_{{\rm{ir}}}\to 0}-2{A}_{{{\Phi }}}^{2}\mathop{\int}\nolimits_{0}^{t}d\tau \ {\rm{Ci}}({\omega }_{{\rm{ir}}}\tau ),$$

(36)

where ${\rm{Ci}}(y)=-\mathop{\int}\nolimits_{y}^{\infty }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 the cosine integral and motivated by physical reasons⁶³. 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 ${\omega }_{{\rm{ir}}}^{-1}$, we make use of the series expansion

$${\rm{Ci}}(w)=\gamma +\mathrm{log}\,(y)+\mathop{\sum }\limits_{k=1}^{\infty }\frac{{(-{y}^{2})}^{k}}{2k(2k)!},$$

(37)

where γ ≃ 0.58 is the Euler’s constant to approximate

$$\mathop{\int}\nolimits_{0}^{t}d\tau \ {{\rm{tr}}}_{B}[{\rho }_{{\rm{B}}}\delta {{\Phi }}(t)\delta {{\Phi }}(t^{\prime} )]\simeq 2{A}_{{{\Phi }}}^{2}\ t\ [(1-\gamma )-\mathrm{log}\,({\omega }_{{\rm{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

$$\begin{array}{lll}{\partial }_{t}\rho &=&\mathop{\sum}\limits _{k}{{{\Gamma }}}_{\varphi }^{kk}\ {\mathcal{D}}[{\sigma }_{kk},{\sigma }_{kk}]\ \rho \\ &+&\mathop{\sum}\limits _{k>l}{{{\Gamma }}}_{\varphi }^{kl}\left({\mathcal{D}}[{\sigma }_{kk},{\sigma }_{ll}]+{\mathcal{D}}[{\sigma }_{ll},{\sigma }_{kk}]\right)\rho ,\end{array}$$

(39)

where ${{{\Gamma }}}_{\varphi }^{kl}$ are time-dependent pure-dephasing rates given by

$${{{\Gamma }}}_{\varphi }^{kl}={\partial }_{{{{\Phi }}}_{{\rm{ext}}}}H{| }_{{{{\Phi }}}_{{\rm{ext}}}^{0}}^{kk}{\partial }_{{{{\Phi }}}_{{\rm{ext}}}}H{| }_{{{{\Phi }}}_{{\rm{ext}}}^{0}}^{ll}\times 4{A}_{{{\Phi }}}^{2}\ t\ [(1-\gamma )-\mathrm{log}\,({\omega }_{{\rm{ir}}}t)]/{\hslash }^{2},$$

(40)

${\sigma }_{kl}=\left|{\psi }_{k}\right\rangle \left\langle {\psi }_{l}\right|$, and ${\mathcal{D}}[x,y]\ \rho =x\rho {y}^{\dagger }-\{{y}^{\dagger }x,\rho \}/2$ is a generalized dissipator superoperator. Equivalently, Eq. (39) can be recast in the more familiar form

$$\begin{array}{lll}{\partial }_{t}\rho &=&\mathop{\sum}\limits _{k}{{{\Gamma }}}_{\varphi }^{kk}\ {\mathcal{D}}[{\sigma }_{kk}]\ \rho \\ &+&\mathop{\sum}\limits _{k>l}{{{\Gamma }}}_{\varphi }^{kl}\ \left({\mathcal{D}}[{\sigma }_{kk}+{\sigma }_{ll}]-{\mathcal{D}}[{\sigma }_{kk}]-{\mathcal{D}}[{\sigma }_{ll}]\right)\rho ,\end{array}$$

(41)

where ${\mathcal{D}}[x]\ \rho =x\rho {x}^{\dagger }-\{{x}^{\dagger }x,\rho \}/2$ is the standard dissipator superoperator. By projecting Eq. (41), one has

$$\langle {\psi }_{k}| {\partial }_{t}\rho | {\psi }_{l}\rangle =-\frac{1}{2}\left[{{{\Gamma }}}_{\varphi }^{kk}+{{{\Gamma }}}_{\varphi }^{ll}-2{{{\Gamma }}}_{\varphi }^{kl}\right]\langle {\psi }_{k}| \rho | {\psi }_{l}\rangle ,$$

(42)

where

$$\left[{{{\Gamma }}}_{\varphi }^{kk}+{{{\Gamma }}}_{\varphi }^{ll}-2{{{\Gamma }}}_{\varphi }^{kl}\right]\propto {[{\partial }_{{{{\Phi }}}_{{\rm{ext}}}}(\hslash {\omega }_{kl}){| }_{{{{\Phi }}}_{{\rm{ext}}}^{0}}]}^{2}.$$

(43)

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}_{{{\Phi }}}^{4}$, most devices are simply T₁-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

$${\rho }_{kl}(t)={\rho }_{kl}(0)\exp \left\{-{A}_{{{\Phi }}}^{2}{({\partial }_{{{{\Phi }}}_{{\rm{ext}}}}{\omega }_{kl}{| }_{{{{\Phi }}}_{{\rm{ext}}}^{0}})}^{2}\ {t}^{2}\ \left[\left(\frac{3}{2}-\gamma \right)-\mathrm{log}\,({\omega }_{{\rm{ir}}}t)\right]\right\}.$$

(44)

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 $(\frac{3}{2}-\gamma )$ within brackets in Eq. (44). Finally, we define the coherence time T_φ as the solution of the implicit equation ρ₀₁(T_φ)/ρ₀₁(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.