Probing the many-body localization phase transition with superconducting circuits

Chains of superconducting circuit devices provide a natural platform for studies of synthetic bosonic quantum matter. Motivated by the recent experimental progress in realizing disordered and interacting chains of superconducting transmon devices, we study the bosonic many-body localization phase transition using the methods of exact diagonalization as well as matrix product state dynamics. We estimate the location of transition separating the ergodic and the many-body localized phases as a function of the disorder strength and the many-body on-site interaction strength. The main difference between the bosonic model realized by superconducting circuits and similar fermionic model is that the effect of the on-site interaction is stronger due to the possibility of multiple excitations occupying the same site. The phase transition is found to be robust upon including longer-range hopping and interaction terms present in the experiments. Furthermore, we calculate experimentally relevant local observables and show that their temporal fluctuations can be used to distinguish between the dynamics of Anderson insulator, many-body localization, and delocalized phases. While we consider unitary dynamics, neglecting the effects of dissipation, decoherence, and measurement back action, the timescales on which the dynamics is unitary are sufficient for observation of characteristic dynamics in the many-body localized phase. Moreover, the experimentally available disorder strength and interactions allow for tuning the many-body localization phase transition, thus making the arrays of superconducting circuit devices a promising platform for exploring localization physics and phase transition.

Figure 1

Schematic of a transmon chain realizing the Bose-Hubbard model with attractive interactions. A transmon is made of Josephson junctions (black crosses) and capacitor plates (blue rectangles). Transmons are anharmonic oscillators with natural frequency $\omega$ and anharmonicity $U$ and they interact with each other through capacitive interaction $J$. Nonidentical sites visualize fabrication scatter and disorder.

Figure 2

Normalized eigenenergies of the clean (a)–(c) and the disordered (d)–(f) Bose-Hubbard Hamiltonian at half-filling and $L=10$ for anharmonicities $U/J=0$ (a) and (d); $U/J=3.5$ (b) and (e); $U/J=20$ (c) and (f). The eigenenergies are scaled as $\epsilon =(E-E_{\min })/(E_{\max }-E_{\min })$. The disorder realization, drawn from a uniform distribution with the disorder strength $W/J=10$, is the same for each $U$. Black dots show the expectation value of the total anharmonicity operator ${\sum }_{\ell }{\widehat{n}}_{\ell }({\widehat{n}}_{\ell }-1)$ in the corresponding eigenstate.

Figure 3

Behavior of the mean energy $\overline{\omega }$ and corresponding standard deviation ${\sigma }_{\mathrm{transmon}}$ of the transmon potential as a function of magnetic field. The fact that the mean energy (blue curve) is close to the “middle” energy ${\omega }_{\text{mid}}=({\omega }_{\text{max}}+{\omega }_{\text{min}})/2$ (blue dashed curve) at high magnetic fields suggests that the uniform distribution provides a good approximation to the true distribution of ${\omega }_l$. The standard deviation (green curve) of the transmon disorder potential of Eq. (7) as a function of magnetic field also becomes close to the one of the uniform distribution ${\sigma }_{\mathrm{uniform}}=W/\sqrt{3}$ (dashed green curve) for strong magnetic field. In contrast, weak magnetic fields give the narrow distribution of transmon energies sharply peaked around the mean energy. Junction asymmetry is $d=0.1$, loop areas are assumed to be Gaussian distributed with standard deviation $1/5$ of the mean $\overline{A}$, and all the values are expressed in the units of unperturbed transmon energy $\sqrt{8E_{\mathrm{C}}{E}_{\mathrm{\Sigma }\mathrm{J}}}$.

Figure 4

The density of states $\rho \left(\epsilon \right)$ in the disordered Bose-Hubbard model of Eq. (8) as a function of normalized energy $\epsilon$ for three different disorder realizations with parameters $L=10$, $U/J=3.5$, and $W/J=10$. The location of the maximum varies due to interplay of anharmonicity and disorder. The realization 1 (blue) is the same as in Fig. 2.

Figure 5

The bipartite entanglement entropy $S$ (a) and the bipartite number uncertainty $F$ (b) as a function of the disorder strength $W$ of uniform disorder distribution for different system sizes $L=8\text{(blue),}10\text{(green),}12\text{(red),}\text{and}14\text{(yellow)}$ with $U/J=3.5$. The eigenstate, for which the properties are calculated, is the one closest to the maximum density of states. The results are averaged over 4000 disorder realizations, except in $L=14$ we have 264 realizations. Dots denote our data points, and the curves are polynomial fits for visibility. The standard error is denoted with the error bars, and for shorter chains they are smaller than the marker size. In the insets, we present the collapsed data from a finite-size scaling analysis with the ansatz $g\left[L^{1/\nu }(W-W_{\mathrm{c}})\right]$ without the shortest chain $L=8$. With this value of $U$, we obtain critical disorder strength $W_c/J=8.24$ and scaling exponent $\nu =1.126$ for the bipartite entanglement entropy and for the bipartite number uncertainty $W_c/J=9.36$ and $\nu =1.24$.

Figure 6

Average adjacent gap ratio $r$ as a function of the disorder strength $W$ for different sized systems with anharmonicity $U/J=3.5$. The reference values for Gaussian orthogonal ensemble ($r_{\mathrm{GOE}}$) and Poissonian distribution ($r_{\mathrm{Poisson}})$ are denoted by black horizontal lines. The results are averaged over 4000 realizations, except for $L=14$, which is averaged over 264 realizations. The eigenvalues are located at the maximum density of states, and the gap ratio is calculated over 16 adjacent eigenstates from each disorder realization. The error bars denote the standard error and are in general smaller than the marker size.

Figure 7

Phase diagram for the half-filled disordered Bose-Hubbard Hamiltonian (8) as a function of the on-site interaction $U$ and the disorder strength $W$. The transition point is estimated using finite-size scaling analysis for different sized systems $L=10,12,14$ for the bipartite entanglement entropy and the bipartite number uncertainty, shown in Fig. 5. The disorder is included only in the on-site energy ${\omega }_{\ell }$ and it is drawn from the uniform distribution $[-W,W]$. The data points are obtained as an average of the transition points given by the finite-size scaling analysis on the bipartite entanglement entropy and the bipartite number uncertainty and the error bars show the deviation between these two values. The white curve is a polynomial fit that improves visual clarity.

Figure 8

The bipartite entanglement entropy $S\left(t\right)$ [(a) and (d)], the temporal number fluctuation ${\mathcal{T}}_{\mathrm{even}}\left(t\right)$ of even sites [(b) and (e)], and the two-site number correlations ${\mathcal{C}}_r\left(t\right)$ [(c) and (f)] as a function of time $t$. The dynamical probes are calculated for interacting system with $U/J=3.5$ (solid lines) for system sizes $L=8\text{(blue),}10\text{(green),}12\text{(red),}14\text{(yellow),}\text{and}40\text{(purple)}$, and for the noninteraction situation $U/J=0$ with $L=12$ sites (dashed red line). The disorder distribution is uniform, and panels (a)–(c) correspond to strong disorder $W/J=15$ and panels (d)–(f) to weak disorder $W/J=1$. Results for system sizes $L\le 12$ are averaged over 2000, for $L=14$ over 1000, and for $L=40$ over 500 disorder realizations. For $L=40$, additional truncation of on-site Hilbert space was applied with maximal occupancy limited to $n_{\max }=4$ for $W/J=1$ and $n_{\max }=3$ for $W/J=15$, which is justified by the dilute initial state and rather short evolution times. Shaded regions depict the standard errors of the disorder averages. In panels (c) and (f), the yellow dots depict the points at which the corresponding correlations reach the predetermined reference value (dashed horizontal line). The insets in panels (d) and (e) describe the equilibrium values of the corresponding observables as a function of system size. Note also the different time scales in left and right columns.

Figure 9

Comparison of the density of states histograms for system with $L=12$ using $LDL$ matrix decomposition (solid blue) and Chebyshev expansion (dashed blue). Relative error percentage is shown as the red histogram and it tells how much the results obtained with the Chebyshev expansion differ from those obtained with the $LDL$ method. The Chebyshev method is the most accurate in intervals with large number of states. Around the maximum density of states the error is $\approx 1\%$. In the Chebyshev expansion, we use 50 terms and the stochastic trace is obtained with 30 random vectors.

Figure 10

Relative error percent of the bipartite entanglement entropy obtained from the exact diagonalization and Krylov method with $m=5$ as a function of evolution time calculated for a single, representative disorder realization. The system parameters of Eq. (8) are $L=12$, $W/J=15$, and $U/J=3.5$.

Figure 11

Evolution of the average bond dimension $\mathcal{D}$ in ergodic (green) and many-body localized (blue) phases.