Simulating dirty bosons on a quantum computer

The physics of dirty bosons highlights the intriguing interplay of disorder and interactions in quantum systems, playing a central role in describing, for instance, ultracold gases in a random potential, doped quantum magnets, and amorphous superconductors. Here, we demonstrate how quantum computers can be used to elucidate the physics of dirty bosons in one and two dimensions. Specifically, we explore the disorder-induced delocalized-to-localized transition using adiabatic state preparation. In one dimension, the quantum circuits can be compressed to small enough depths for execution on currently available quantum computers. In two dimensions, the compression scheme is no longer applicable, thereby requiring the use of large-scale classical state vector simulations to emulate quantum computer performance. In addition, simulating interacting bosons via emulation of a noisy quantum computer allowed us to study the effect of quantum hardware noise on the physical properties of the simulated system. Our results suggest that scaling laws control how noise modifies observables versus its strength, the circuit depth, and the number of qubits. Moreover, we observe that noise impacts the delocalized and localized phases differently. A better understanding of how noise alters the genuine properties of the simulated system is essential for leveraging noisy intermediate-scale quantum devices for simulation of dirty bosons, and indeed for condensed matter systems in general.


I. INTRODUCTION
Scientific investigation of dirty bosons [1,2] seeks to characterize the phase diagram (phases of matter, transitions) of strongly interacting bosons in a disordered environment.The combination of disorder and interactions present in these systems is a double-edged sword: one the one hand, it induces interesting quantum physics in the form of an exotic new disorder-induced phase of matter at low temperatures, namely, the Bose glass; on the other hand, it makes the study of these systems a challenging problem in condensed matter physics.Interest in dirty bosons traces back to experiments of superfluid 4 He in porous media [3,4], along with the theoretical work that followed, which identified the Bose glass phase [5][6][7][8][9].While bosons typically display superfluidity and Bose-Einstein condensation in two dimensions and higher (often associated with a global quantum coherence), the Bose glass phase is characterized as gapless, compressible, short-range correlated, and localized.Therefore, upon sufficient disorder, a bosonic system can experience a delocalized-to-localized phase transition accompanied with extreme changes in its physical properties.It is analogous to the Anderson metal-to-insulator transition for noninteracting electrons [10][11][12][13].
The problem of dirty bosons has also been considered theoretically and numerically with classical simulation methods (e.g., quantum Monte Carlo and the density matrix renormalization group) [5][6][7][8][9].However, despite extensive investigation into dirty bosons, both experimentally and with simulation on classical computers, several fundamental questions remain open: Is there an upper critical dimension for the superfluid-to-insulator transition [9]?Can there be an intermediate scenario between the delocalized and localized phases with a nonergodic delocalized phase similar to what is found for the Anderson localization problem on certain geometries [71][72][73]?Is there a connection between this problem in infinite dimension and the many-body localization phenomenon [66,74,75]?Moreover, some bosonic phases of matter are hard to simulate classically due to the existence of a sign problem in quantum Monte Carlo simulations.Quantum computers, with their ability to efficiently simulate quantum systems [76][77][78][79] and finely tune system parameters (e.g, disorder strength), offer a promising new platform for investigating such phases, e.g., Bose metals [80,81] or superglasses [82,83].
Already, a profusion of recent work has successfully demonstrated simulations of varied properties and behaviors of quan-arXiv:2210.08386v1[cond-mat.dis-nn]15 Oct 2022 tum systems [84][85][86][87][88][89][90][91] on currently available quantum hardware.Here, we demonstrate how quantum computers can be used to simulate dirty bosons.We focus on the toy model of hard-core bosons on a lattice in a random potential of tunable strength in one and two dimensions.We demonstrate how adiabatic state preparation [92] can be used to prepare the ground state of the system on the quantum computer in order to study its properties.In order to ensure adiabaticity, long evolution times are generally required.Simulating this evolution generally requires deep quantum circuits, as quantum circuit depths tend to grow linearly with evolution time [93,94].This poses a problem for current quantum computers.This is because currently available quantum hardware is noisy, leading to short qubit decoherence times and high gate error rates, which sets an effective limit on the depth of quantum circuits that are feasible to execute [95].
In one dimension, we show how recently developed circuit compression techniques for free-fermionic circuits [96][97][98][99] can be used to compress simulation of arbitrarily long time evolution of the dirty bosons into highly compact quantum circuits, thereby ensuring the adiabaticity of the evolution as well as the feasibility of successful circuit execution on current quantum hardware.We were therefore able to obtain quantitative experimental results from IBM's quantum processors on the nature of the superfluid-to-insulator transition by extracting its critical properties such as the correlation length exponent.
In two dimensions, the compression trick is no longer applicable and we therefore restrict our results to large-scale classical state vector simulations, which emulate quantum computation on classical computers, leveraging NVIDIA cuQuantum library on graphics processing units [100].Classical emulation of an ideal (i.e., noise-free) quantum computer allowed us to validate our quantum computational approach to obtaining the ground state physics of dirty bosons in two dimensions.
Finally, we perform simulations of dirty bosons in one and two dimensions via classical emulation of a noisy quantum computer, based on a depolarizing channel of Pauli operators.This allowed us to explore how the noise inherent in current quantum hardware affects the physical properties of the final quantum state of the simulated system.In particular, our results suggest that scaling laws control how hardware noise modifies observables versus its strength, the circuit depth, and the number of qubits.Moreover, we observe that noise impacts the delocalized and localized phases differently.We believe our results on the noise-induced physics will prove useful for simulating bosons on near-future quantum computers, which are expected to comprise more and better quality qubits.
The rest of the paper is organized as follows.In Sec.II, we define the models, physical observables of interest and the strategies used to simulate them.In Sec.III, we present and discuss our simulation results for dirty bosons in one and two dimensions, including an analysis of the observed noiseinduced physics.We conclude by summarizing our findings and giving perspectives in Sec.IV.

A. Model: Interacting bosons on a lattice
We consider a toy model of hard-core bosons at half-filling on a lattice, described by the following Hamiltonian, where b †  and b are the bosonic creation and annihilation operators on lattice site , respectively.They obey the usual bosonic commutation relation b is the local density operator on site  with the constraint n ≤ 1 due to the hard-core nature of the particles.The constraint is equivalent to an infinite repulsive interaction, preventing more than one boson being on a given lattice site.Hence, the local Hilbert space of the lattice site  is given by |0  and |1  , corresponding to an empty and occupied lattice site , respectively.The first term of Eq. ( 1) describes the hopping on nearest-neighbor lattice sites.Its negative sign favors delocalization on the lattice.The second term describes a sitedependent random chemical potential drawn from a uniform distribution   ∈ [−, +] with  characterizing the strength of the disorder in the lattice.Disorder may favor localization of the particles, leading to a competition with the kinetic term.We consider two lattices: (i) a one-dimensional linear chain with open boundary conditions, shown schematically in Figure 1a and (ii) a two-dimensional square lattice with periodic boundary conditions, shown schematically in Figure 1b.We discuss the physics of each model and existing literature in their respective result section.
Up to an irrelevant constant, the Hamiltonian (1) can be expressed in terms of Pauli operators by a Matsubara-Matsuda mapping [101], equivalent to interacting spin-1/2 in the XY plane subject to a random spin-dependent magnetic field along the Z direction.The form of Eq. ( 2) directly translates into a quantum computing approach.We consider two physical quantities to investigate the ground state properties of the dirty bosons.First, the local density, Then the off-diagonal two-point correlation, Due to the random nature of the models considered, we denote disorder-averaged observables over  samples realizations by (• • • ).With no disorder ( = 0), the bosons are delocalized, but with the onset of any finite value of disorder, the bosons will localize.(b) In two dimensions, we simulate a square lattice with periodic boundary conditions, with bosons initialized in a checkerboard pattern.Here, we show the 16-site lattice that was simulated with the noisy quantum circuit emulator.Up to a disorder strength of  ≈ 9.6, the bosons are delocalized, but at higher disorder strength, the bosons localize.

B. Ground state preparation
We are interested in the ground state properties of the Hamiltonian (1) in one and two dimensions.We obtain the ground state by performing an adiabatic state preparation [92] .We start by preparing the system in a product state |Ψ ini corresponding to a checkerboard pattern of occupied and empty lattice sites, depicted in Figure 1.A parent Hamiltonian to this state is, with  () = 0 if site  is initially empty and  () = 1 otherwise.We define the Hamiltonian interpolating between Eqs. ( 5) and (2), with  ∈ [0, ]. is fixed and corresponds to the total time for the unitary evolution interpolating between the two Hamiltonians (we set = 1), where |Ψ fin is exactly the ground state of ĤXY in the adiabatic limit  → +∞.T indicates a time-ordered exponential.The evolution (7) conserves the number of particles in the system, which is a symmetry of the final Hamiltonian [102].
In practice, the total evolution time  should be chosen such that the quantum state remains in the instantaneous ground state of the interpolating Hamiltonian Ĥ ,  throughout the time evolution.The Landau-Zener formula states that this condition is fulfilled for  Δ −2 min with Δ min the minimum spectral gap between the ground state and the first excited state along the interpolating path [103].Here, because the time evolution conserves the number of particles, Δ min corresponds to the gap within this particle-conserving sector of the interpolating Hamiltonian [104].The choice of  is discussed later following the expected physics of each model investigated.

Quantum circuit
On a digital quantum computer, one needs to discretize the unitary evolution (7) into  steps with a time step  such that  =  steps , Thanks to the local nature of the interpolating Hamiltonian, Eq. ( 8) can be broken down into one-and two-site unitaries by a Suzuki-Trotter decomposition [105].At first order, one gets, with  ( 2 ) corrections not taken into account.The single-site diagonal unitary is defined as, and the two-site off-diagonal unitary reads, The single-site unitary is equivalent to a parametric one-qubit rotation gate R Z () = exp(− Ẑ/2) and the two-site unitary to a parametric two-qubit gate R XY () = exp[( X X + Ŷ Ŷ )/4] with the angles  according to Eqs. ( 10) and (11).
From the state |000 . . .0 , the initial state |Ψ ini with a checkerboard pattern of empty/occupied sites is easily obtained by applying individual X gates on the desired occupied lattice sites.For each step, the one-qubit gates of Eq. ( 9) can be applied in parallel.The two-qubit gates of Eq. ( 9) are arranged to minimize circuit depth and applied in parallel accordingly: one layer for even and odd bonds, respectively, in one dimension.On the two-dimensional square lattice, four layers are needed to accommodate the four-fold connectivity of each lattice site.The local density of Eq. ( 3) is easily computed from bitstrings obtained in the computational basis.To compute the off-diagonal two-point correlation of Eq. ( 4), we make a basis rotation before measurement such that the operator X is diagonal by applying Hadamard gates on the individual qubits.Expectation values are computed by performing the average over a finite set of bitstrings  shots .

The special case of one dimension
In one dimension, the interpolating Hamiltonian of Eq. ( 6) maps to a spinless free-fermionic system through a Jordan-Wigner transformation [106].Due to the free-fermionic nature of the system, the resulting quantum circuit for the state preparation can be compressed to a final depth that scales linearly with the number of qubits  [96][97][98][99], as opposed to a depth that scales linearly with evolution, which is generally the case [93,94].Hence, the quantum circuit depths are independent of simulation time, enabling simulation out to arbitrarily long times, ensuring adiabaticity of the evolution.Albeit smaller, the circuit keeps a similar structure after the compression, with consecutive layers of R Z () and R XY () gates and rescaled angles  following the compression scheme.We generate compressed circuits using the open-source software package fast free fermion compiler (F3C) [107,108], which in turn is built based on the QCLAB toolbox [109,110] for creating and representing quantum circuits.

Emulations on a classical computer
Currently available quantum computers are noisy, mainly due to relatively short qubit decoherence times and relatively high gate-error rates.This puts a limit on the size of quantum circuits that can be feasible run (i.e., return high-fidelity results) on real quantum computers, which in turn limits the system sizes that are feasible to simulate.To gauge how simulations of larger systems may perform on improved quantum computers of the near-future, such simulations may be emulated on a classical computer.So-called quantum circuit emulation uses a classical computer to simulate execution of quantum circuits by storing the many-body wavefunction of an -qubit system (a complex vector of size 2  ) and updating it according to the gates in the emulated quantum circuit.This will provide effective simulation results from an ideal (i.e., noise-free) quantum computer.It is also possible to emulate a noisy quantum computer computer by adding noise channels to the emulation process, which can be tuned with a parameter that effectively sets the probability of error.This allows us to study the effect of noise on simulation results.Here, we use qsim [111], accelerate with the NVIDIA cuQuantum SDK [100], to emulate both ideal and noisy quantum computers on a classical GPU backend.All our classical emulations were performed on the NERSC Perlmutter supercomputer.

Simulations on a quantum computer
For smaller systems sizes (around 10 or fewer qubits), it is possible to obtain reasonable results from current quantum computers.Here, we use IBM's 'ibmq_brooklyn' quantum processing unit (QPU), comprising qubits implemented with superconducting circuits.While quantum circuits emulated on a classical computer may utilize arbitrary unitary gates, circuits executed on QPUs may only use so-called native gates, a small, universal set of gates that can be natively implemented on the quantum computer.The native gate set on IBM's QPU includes the √︁ X gate, which applies the square root of the Pauli-X gate; the R Z () gate, which rotates the qubit about the -axis by an angle ; and the two-qubit, entangling CNOT gate.This is a universal gate set because any unitary gates can be decomposed into a sequence of these three native gates.Note that circuits compiled down to the native gate set can be significantly longer than their counterparts using arbitrary unitary gates.For instance, the two-qubit R XY () gate of Eq. ( 11) decomposes as, where R X () = exp(− X/2) rotates the qubit about the axis by an angle  and which can be further decomposed using √︁ X =   /4 R X (/2) together with R Z (), with , (14) the Hadamard gate.Global phases are irrelevant.

Review of the existing literature and expectations
In one dimension, the bosonic model of Eq. ( 1) can be mapped to a spinless fermionic model with local terms by a Jordan-Wigner transformation [106], where ĉ †  and ĉ are the fermionic creation and annihilation operators on lattice site , respectively [112].They obey the usual fermionic anticommutation relation ĉ is the local density operator on site .Unlike the hard-core bosons, effectively interacting through an infinite repulsive interaction, the fermions are non-interacting in Eq. ( 15).The Hamiltonian is quadratic, and its eigenstates are easily constructed by filling one-particle orbitals.In the absence of a disordered potential ( = 0), the system displays superfluidity and quasi-long-range order with critical properties captured by Tomonaga-Luttinger liquid physics [44].For any nonzero disorder strength  > 0, all the single-particle orbitals will experience Anderson localization, and therefore, all the eigenstates of Eq. ( 15) will also be localized [6,7], including its ground state.The localization is characterized by a localization length scale  diverging as  ∼  −2 as  → 0 [113], such that the critical behavior is recovered at  = 0.In the other limit, at strong disorder, the localization length decreases as  ∼ (ln ) −1 [114].
Hence, a localized phase exists for  > 0 and delocalization only at  = 0, which is also the transition point.The phase diagram is depicted in Figure 1a.In order to stabilize a delocalized phase for hard-core bosons at finite , one can add a nearest-neighbor interaction term of the form 2  n n in the Hamiltonian, with strength .If the interaction is attractive and in the range  ∈ (−1, −1/2), such that disorder is irrelevant in a renormalization group sense, delocalized Tomonaga-Luttinger liquid physics can be robust up to  0.74 [6,7,63].Note that in that case, the localization length diverges differently as the delocalized phase is approached, with  ∼  2/(2 −3) where  = /2 arccos(−) the so-called Tomonaga-Luttinger liquid parameter [115]-for the case  = 0 considered in this work we have  = 1.

Gap scaling and state preparation protocol
For the state preparation to be accurate, it needs to be adiabatic, which requires prior knowledge of the minimum spectral gap Δ min between the instantaneous ground state and the first excited state during the evolution protocol.One can then choose a total evolution time  Δ −2 min according to the Landau-Zener formula.
Instead of Δ min which is cumbersome to evaluate, we consider the gap Δ of the final state.We perform exact classical simulations diagonalizing the final Hamiltonian.We show in Fig. 2(a) the average gap computed over more than  samples = 10 5 random disorder realizations versus the system size for various disorder strengths .In absence of disorder, the gap follows Δ ∼  −1 with  the system size, which is compatible with the critical properties of the Tomonaga-Luttinger phase.We find the same scaling in presence of disorder, see Fig. 2(a).However, the gap is widely distributed over several orders of magnitude.We illustrate this for  = 8 with the cumulative distribution function of the gap CDF(Δ) in Fig. 2(b).We observe that as the disorder strength increases, there exist random realizations with smaller and smaller gaps.
In the following, we decide to take a total evolution time  = 10 4 (along with a time step  = 0.05), which will lead to an adiabatic state preparation for gaps Δ min 10 −2 .For  = 8 and the disorder strengths  ≤ 10, Fig. 2(b) states that the probability to find such a disorder configuration with a gap smaller than ≈ 10 −2 is less than ≈ 0.1%.Hence, at most, a disordered sample for every one thousand will not be prepared adiabatically.

Local density of particles
We first consider the local density of particles of Eq. (3).For each disordered sample of length , we collect  expectation values and compute the probability distribution ( n ) with n ∈ [0, 1].The quantity displays no size effect (data not shown) and we consider  = 6.Exact diagonalization data is plotted in Fig. 3(a).In absence of disorder, the system being at half-filling is reflected by a delta-peaked distribution at n = 1/2.As the disorder strength increases, the distribution spreads out until at larger disorder the distribution takes a U-shape with the maximum probability for n → 0, 1.At strong disorder, the bosons will be predominantly present on lattice sites with a negative chemical potential   < 0 to minimize the energy of the system, in line with a localization 0 0.25 0.5 0. picture.In fact, it has been found to lead in the thermodynamic limit to a chain breaking mechanism as  → +∞ [74], holding in presence of nearest-neighbor interaction and for excited states in the context of the many-body localization phenomenon [116].
Using the same procedure, we plot the  = 6 qubits experimental probability distribution of the local density of particles in Fig. 3(b).We considered a hundred random realizations for each disorder strength.For each disordered sample, the expectation value of the local density is computed from  shots = 2 13  measurements according to Eq. (3).The general features of the exact data are reproduced experimentally: From small to strong disorder, the maximum shifts from half-filling to the extrema values n → 0, 1 to form a U-shaped distributionthough the U shape doesn't extend all the way to n → 0, 1.As we will see in Sec.III C, the discrepancy is due to inherent hardware noise which favors delocalization.

Delocalization-localization transition
We now turn our attention to the two-point off-diagonal correlator of Eq. ( 4).For a system of length , we take the reference lattice site to be /2 and note the distance with the other lattice site by .We plot in Figs.4(a), 4(b), and 4(c) the correlator versus the distance  for various disorder strengths  for three cases.The first case is from an ideal benchmark based on exact diagonalization data on  = 64.The second case is from a perfect simulation on  = 6 based on exact diagonalization.The third case is the experimental data on  = 6 from the quantum computer.
In absence of disorder, the system is critical and displays a power-law decay of the form  ( = 0, ) ∼ 1/ √  for  → +∞ and  → +∞.The expectation is that, close to the transition, the disorder will modify this scaling by introducing a length scale , To isolate the scaling function F (), we start by plotting in Figs.4(d), 4(e), and 4(f) the correlator divided by its disorderfree value.For a second order transition, the length scale will diverge algebraically as the critical point is approached, i.e.,  ∼  − with  the correlation length critical exponent.Hence, plotting the data as a function of the variable   should lead to a collapse onto a single function F ().The goal is now to find the value of  leading to the best data collapse.We employ the procedure of Ref. 117 which Taylor expands the unknown scaling function, with an extra exponential term accounting for the rapid decay observed in Fig. 4(d) that would otherwise require a high expansion order  to be captured in practice.In Eq. ( 17),  ≡  (, )/ ( = 0, ) and  () ≡   .For a given value of the exponent , we perform a least-square fitting of the data with parameters ã,  0 ,  1 , etc.The quality of the fit for  pairs of data points {  ,   } is measured by the chi-square, with Δ  the error on the data point   .We repeat the procedure for a grid of the exponent  and find the value leading to the best data collapse as the one minimizing  2 in Eq. (18).We use  = 3 in practice and plot the results in Figs.4(j), 4(k), and 4(l).The value of  minimizing  2 is then used in Figs.4(g), 4(h), and 4(i) to rescale accordingly the  axis.We theoretically expect  = 2 and find that the value is perfectly recovered in the ideal benchmark case [6,7,113].
Although less satisfying, we can obtain a data collapse for the perfect simulation on  = 6 with an exponent  ≈ 1.96.
The procedure that was carried on the ideal and perfect data is repeated on the experimental data, see Figs. 4(c), 4(f), 4(i), and 4(l).Each data point corresponds to an average over a hundred disordered realizations, where the correlator of Eq. ( 4) has been computed from  shots = 2 13 measurements.As expected, we observe a decay of the correlator with the distance and the hierarchy of the different curves with respect to the disorder strength is broadly respected when compared to the perfect simulation results.Proceeding with the data collapse, we find that the chi-square is minimized for  ≈ 1.86.The resulting scaling curve is qualitatively similar to the perfect simulation one.The difference is attributed to the inherent hardware noise.We study its effect on the physics of dirty bosons in Sec.III C. The correlation length exponent  ≈ 1.86 is close to the expected theoretical value of  = 2. μ = 0, 0.5, 1,…,4

Review of the existing literature and expectations
In two dimensions, the physics of the model ( 1) is different than its one-dimensional counterpart.In the thermodynamic limit, the larger dimensionality allows for a spontaneous breaking of the continuous  (1) symmetry of the model.Up to a critical disorder strength  c 9.6 [60,61], the symmetry is spontaneously broken and the bosons form a Bose-Einstein condensate with true long-range order, unlike the one-dimensional case which was limited to quasi-longrange order.We note that on a finite system,  (1) cannot be spontaneously broken-though one can engineer quantum states breaking it explicitly [118].At larger disorder strength, the system enters a Bose glass phase where condensation and superfluidity are destroyed and localization takes place.The phase diagram is depicted in Figure 1b.The transition between the two phases is second order, with critical exponents that have been estimated by classical quantum Monte Carlo simulations [35,37,45,57,61].A peculiarity of the superfluid-toinsulator transition is the conjecture that the dynamical exponent  equals the dimensionality  of the model [9], which is in line with quantum Monte Carlo estimates in  = 1, 2, and 3.

Delocalized versus localized bosons
a. Scaling of the gap-To choose the correct total time evolution for the protocol to be adiabatic, we consider the spectral gap Δ between the ground and first excited states in the half-filling sector of a  = 4 × 4 system [119].We perform exact diagonalization calculations on the final Hamiltonian for various disorder strengths  and plot the cumulative distribution function CDF(Δ) in Fig. 5(a) using over  samples = 10 5 random realizations.As the disorder strength increases, there exist random realizations with smaller and smaller gaps.For  19) for  = 4 × 4 as a function of the disorder strength .Each data point is averaged over a hundred random realizations whose ground state is obtained by an adiabatic state preparation circuit.From quantum Monte Carlo simulations, we expect the transition between the superfluid and insulating phases to take place at  c ≈ 9.6 [60,61].
= 15 (one of the largest disorder strengths considered), we find that the minimum gap is ≈ 0.1, which requires a total time evolution  = 100 for adiabaticity.We use this value along with a time step  = 0.005 in the following.

b. Local density of bosons-
We now perform simulations using a state preparation protocol based on a quantum circuit for  = 4 × 4. At the end of the circuit, we compute the local density of bosons of Eq. ( 3).We show their distribution in Fig. 5(b) for various disorder strengths .In absence of disorder ( = 0), all individual lattice sites are at half-filling.As the disorder strength increases, the distribution broadens toward a U shape, similar to the one-dimensional case studied in Fig. 3.

c. Coherent density of bosons-
We consider the momentum distribution function based on the two-point offdiagonal correlation function of Eq. ( 4), with  the momentum and   the distance between lattice sites  and  [120].In translationally invariant systems, Bose-Einstein condensation (BEC) typically occurs when a macroscopic fraction of bosons occupies the single momentum component  = 0, and Λ =0 plays the role of the order parameter for the symmetry-breaking BEC phase.In systems that are not translation invariant (e.g., because of disorder) some bosons may condense in other momentum components  ≠ 0, and it is more rigorous to consider the one-body density matrix together with the Penrose-Onsager criterion [121][122][123] to assess condensation.Nevertheless, despite subtle differences (see, e.g., Ref. 124 and references therein), Λ =0 , known as the coherent density, has been found to provide accurate information for the problem of interest in two dimensions [60,61].
For  = 4×4 and  = 6×4, we compute the coherent density as a function of the disorder strength.We plot the results in Fig. 5(c), where one observes that at fixed system size, it decays as the disorder strength increases.As  → ∞, one expects Λ =0 ( <  c ) → constant and Λ =0 ( ≥  c ) → 0 in the delocalized and localized phases, respectively.With only two system sizes, we cannot perform a finite-size scaling analysis.We note that at criticality ( =  c ), a system of linear length  will display a behavior of the form Λ =0 ( =  c ) ∼  −− , with  and  the dynamical ab anomalous exponents [60,61], respectively.

Noise model
We now investigate the effect of the noise on the physics of dirty bosons.To that end, we consider a stochastic Pauli error model with quantum channels E  and E  applied after each one-and two-qubit gates, respectively [125], with ρ the density matrix of the system, σ0  ≡ Î , σ1  ≡ X , σ2 ≡ Ŷ , and σ3  ≡ Ẑ are Pauli matrices. ∈ [0, 1] is a parameter controlling the error rate, which for simplicity is taken the same for all qubits and gates.In practice, the parameter  in Eq. ( 20) can be made dependent of the qubits (, ) and indices (, ) and fit a specific quantum device through a noise reconstruction protocol [126,127].This noise model has been been considered successfully in other condensed matter studies [117,[128][129][130].In particular, the noise-depth tradeoff in adiabatic quantum circuits was studied in Ref. 130-here, the depth will be kept fixed.
We simulate the noise model of Eq. (20) using "quantum trajectories" [131,132], which allow to work with pure states instead of the density matrix ρ.Many instances of the same circuit are simulated and each one of them is stochastically subject to the addition of random Pauli gates with a probability ∝  following Eq.(20).The probability that no Pauli gate is added after a one-or two-qubit gate is (1 − ).Here, each shot (output bitstring) is a different stochastic random Pauli realization.
In one dimension, we simulate systems up to  = 20 using a total time evolution  = 5 × 10 4 along with a time step  = 0.05, i.e., a total of 10 6 time steps compressed into a circuit depth  ().For each disorder strength, the data is averaged over  samples = 2 10 random realizations, and observables are computed from  shots = 2 10 measurements, each corresponding to a different quantum trajectory.In two dimensions, we simulate a system size of  = 4 × 4 with periodic boundary conditions.We perform a total time evolution  = 10 2 using a time step  = 0.025, corresponding to a total of 4, 000 steps.The data is averaged over  samples = 10 2 random realizations and observables are computed from  shots = 10 2 measurements.

Noise-induced modifications on the physical properties
It was argued in Ref. 117 that a noise model such as Eq. ( 20) on an adiabatic state preparation protocol was to introduce a length scale  noise ∼ 1/ in the system where  is the circuit depth and  the noise strength.For instance, a two-point correlation function probing critical properties with two local operators separated by a distance  > 0 in real space would be modified by the noise as  (, ) ∼  (,  = 0) × exp(−/ noise ) [133].Another interpretation of the result is that the average number of random Pauli operators in the circuit volume  is .Each of them will reduce the correlation by a positive multiplicative factor  < 1, amounting for a total  (, ) ∼  (,  = 0) ×    in the corresponding circuit volume.This is equivalent to the above result with  noise = −(  ln ) −1 .Along the same lines, for a local observable  defined such that it takes a finite value in a perfect simulation and zero for a random state (broadly defined as what an excessive amount of noise would lead to), one expects that  ( ) ∼  (  = 0) × exp(− ε ) with ε ≡ − ln  > 0. We now seek to characterize the effect of the noise on the physical properties and see how it fits with previous findings.

The local density of particles
First, we consider a quantity based on the local density of Eq. ( 3), which measures how the lattice site  is close to perfect (non)occupation.One gets N  → 1 for perfect (non)occupation n → 0 or 1, while N  = 0 for n = 1/2.In absence of disorder ( = 0), one expects N  = 0 ∀ since we work at half-filling and the Hamiltonian (1) has a particle-hole symmetry, resulting in n = 1/2 ∀.Because N  is already zero in absence of noise, we do not expect noise to affects its value.On the other hand, for a finite disorder strength ( > 0), the quantity N  will be finite, and we can study its dependence as a function of the noise strength .
Results in one dimension for  = 5 are shown in Fig. 6(a) Here,  has to be interpreted as the depth  of the circuit, which scales linearly with  after the free-fermionic compression (b) Two-dimensional data for  = 5 based on a hundred independent disordered samples and  shots = 10 2 random noise measurements for each.For  → 0, we observe an exponential decay ∼ exp(− ε ) with ε = 30(2) × 10 3 a fitting parameter.and show a behavior of the form, with  () a scaling function.For  ≡  1, we find that  () exp(− ε ) with ε = 15.7(2) a fitting parameter, in line with the expectations for a local observable.Here,  has to be interpreted as the circuit depth  following the free-fermionic compression.
In two dimensions, we plot data for  = 4 × 4 in Fig. 6(b).Similar to one dimension, the noise suppresses exponentially the finite value of the observable of Eq. ( 3): As a function of the noise strength , we find exp(− ε ) for  1.Here, the circuit depth is included in the fitting parameter ε = 30(2) × 10 3 , which explains why it is orders of magnitude larger than in one dimension.This highlights the power of the compression scheme in one-dimension to limit the effect of noise by reducing the circuit depth to  ().From general grounds, we expect Eq. ( 22) to remain valid in two dimensions with the modification  → , where  is the circuit depth.

The coherent density of bosons
We now turn our attention to the momentum distribution function of Eq. ( 19) at  = 0.In one dimension, we consider data at criticality ( = 0) and deep in the disordered phase ( = 5).A scaling law involving the noise strength  and the system size  can be obtained, albeit with different scaling variables in the delocalized and localized phases, see Figs.  19) at  = 0. Data normalized by the noiseless case with  = 0. Left column: One-dimensional where each data point is averaged over  samples = 2 10 random disordered samples, and each of those obtained from  shots = 2 10 measurements from noisy simulations.Error bars are smaller than the symbols and not shown.Right column: Two-dimensional case case where each data point is averaged over a hundred random disordered samples, and each of those obtained from  shots = 10 2 measurements from noisy simulations.(a) Data at fixed disorder strength  = 0.For  2  1, we observe a behavior of the form ∼ exp(− ε  2 ) with ε = 3.28(5) a fitting parameter.Here, one of the factors  has to be interpreted as the depth  of the circuit, which scales linearly with  after the free-fermionic compression (b) Data at fixed disorder strength  = 5.For  1 we observe a behavior of the form ∼ exp(− ε ) with ε = 24.5(3) a fitting parameter.Here,  has to be interpreted as the depth  of the circuit, which scales linearly with  after the freefermionic compression.(c) Data at fixed disorder strength  = 5 in the delocalized phase (d) Data at fixed disorder strength  = 15 in the localized phase.In both cases, for  → 0, we observe an exponential decay ∼ exp(− ε ) with ε = 60(3) × 10 3 ( = 5) and ε = 75(2) × 10 3 ( = 15), a fitting parameter.and 7(b).In particular, we find, with G() a scaling function.For  1, we find that  () exp(− ε ) with ε = 3.28(5) ( = 0) and ε = 24.5(3)( = 5), a fitting parameter.In the localized scaling,  has to be interpreted as the circuit depth and in the delocalized scaling, one of the factor  has to be interpreted as such.In one dimension, the circuit depth is linear with  thanks to the free-fermionic compression of the circuit.In absence of noise, for a critical behavior, Ref. 117 found a similar to scaling to Eq. ( 23) for a two-point correlator  (, ).Since Λ =0 is a sum of two-point correlators over all distances , our result fits previous findings.In the localized case, individual correlators can be thought of as local observables extending in real space over a scale corresponding to the localization length .Hence, for system sizes such that  , the coherent density Λ =0 should behave as a local observable with a scaling only involving the circuit depth, fitting the data of Fig. 7(b) and Eq. ( 23).We believe the two behaviors of Eq. ( 23) can be reconciled by a single scaling function G( ) where  is substituted by the length of the system for  .Hence, for   and  → 0 where we expect  ∼  −2 , one would get a scaling function G(   −2 ).However, this regime is out of reach of our simulations for testing.
For completeness, we also perform simulations in two dimensions for  = 4 × 4 and  = 5, 15, see Figs. 7(c) and 7(d).While the lack of multiple system sizes does not permit to draw conclusions on the -dependence of the noise-induced physics, we find that, similar to one dimension, the coherent density is exponentially suppressed as a function of the noise strength.Although we cannot verify it, we expect Eq. ( 23) to remain valid in two dimensions, as it was derived from general arguments (one of the parameters  would need to be replaced explicitly by the circuit depth ).

A. Summary
In this work, we considered the dirty boson problem in one and two dimensions using quantum computing techniques and a programmable quantum computer.In all cases, adiabatic state preparation was used to prepare the ground state of hardcore bosons hopping on a lattice, subject to an on-site random chemical potential of tunable strength.
In one dimension, the model maps to a free-fermionic system, and we leveraged this property to compress the depth of the circuits for the adiabatic state preparation such that it is of the order of the lattice size [96][97][98][99] (as opposed to of the order of evolution time).Upon compression, the circuits were suitable for the current generation of noisy quantum devices.We obtained experimental results on a superconducting quantum computer for  = 6 qubits.We considered the probability distribution of the local density of bosons displaying a distinctive U-shape at strong disorder.We also studied an off-diagonal two-point correlation function and performed a scaling analysis together with a data collapse to extract critical properties of the transition.We estimated the correlation length exponent  ≈ 1.86, in good agreement with the expected theoretical value of  = 2 [6,7,113].
In two dimensions, the model does not map to free fermions and the depth of the circuits cannot be compressed anymore.As the depth is directly proportional to the duration of the adiabatic time evolution protocol, the circuits are too large to be reliably executed on the current generation of quantum hardware.Instead, we performed large-scale classical state vector simulations on graphics processing units [100], which emulate the performance of a quantum computer.Our results from emulation of an ideal (i.e., noise-free) quantum computer val-idated our algorithmic approach to obtaining the ground state physics of two-dimensional disordered systems, demonstrating promising prospects for the next generation of quantum computers with more qubits of better quality.
Finally, we simulated dirty bosons via emulation of a noisy quantum computer, based on a depolarizing channel of Pauli operators with a parameter  controlling the strength of the noise.We sought to understand how noise affects the physical properties of the final quantum state.In the noiseless limit ( → 0), we found that an observable  is generically modified by the noise as  ( ) ∼  (  = 0) ×  −  with  a parameter.In one dimension, we performed a systematic analysis with the system size  and the circuit depth, scaling as .We found that the noise-induced physics is described by scaling relation  ( ) =  (  = 0) × G() with  ()  −  for  1.For local observables, such as the local density, we have  ≡ , while for global quantities, such as the coherent density of bosons, we have  ≡  2 and  ≡  in the localized and delocalized phases, respectively.The different behaviors are understood through the existence of the localization length  making non-local physical observables (e.g., a correlation function) effectively local over this length scale as long as  .We also advanced a scenario reconciling the two regimes as one tunes the delocalization-to-localization transition.We believe that understanding how noise modifies the genuine properties of different phases is fundamental for leveraging noisy intermediate-scale quantum devices for condensed matter, extending the work of Ref. 117.

B. Perspectives
For an adiabatic state preparation protocol to be successful, its total duration  should scale as  Δ −2 min with Δ min the minimum spectral gap between the ground state and the first excited state of the interpolating Hamiltonian [103].In the thermodynamic limit, the Bose glass phase is gapless, meaning that the gap at the end of the interpolation goes as Δ → 0. However, on a finite system, the gap will be finite, but can be exponentially small for some disorder configurations (we studied the distribution of the gap in the main text); this is a general feature of glassy systems [134].The existence of arbitrarily small gaps is also an issue in classical methods like quantum Monte Carlo: it works at a finite temperature Θ-the algorithm's complexity scales linearly with Θ in the absence of a sign problem [135]-which needs to be much smaller than Δ in order to probe ground state properties.Nevertheless quantum Monte Carlo is still the method of choice to simulate thousands of interacting particles [60,61,63], meaning that not perfectly capturing the ground state of configurations with small spectral gaps is not so bad when considering a statistical ensemble over random disordered realizations.Whether this is the case for an adiabatic state preparation remains to be seen.Inspired by the quantum approximate optimization algorithm [136][137][138], one could leave instead the parameters of the evolution gates in Eqs.(10) and (11) as variational parameters to find an adiabatic shortcut in the state preparation.
In addition to the Bose glass phase and fundamental questions that remain pending on its exact nature (see the discussion in the introduction), one could use a quantum computer to simulate other bosonic phases of matter which are hard to simulate classically due to the existence of a sign problem in quantum Monte Carlo simulations.For instance, one could study Bose metals [80,81] or a superglass phase [82,83]-though frustration might not be a necessary ingredient in the latter [139].In particular, we believe our results on the noise-induced physics will prove useful for such future studies, once more and better qubits become available.

FIG. 1 .
FIG. 1. Illustrations of the simulated systems along with their phase diagrams with respect to disorder strength  in one (a) and two (b) dimensions.Diamonds represent lattice sites, whose random coloring represents their random on-site potentials.The blue spheres represent the hard-core bosons, which are allowed to hop between lattice sites connected with an edge.The depicted positions of bosons represents their initial configuration in all simulations.(a) In one dimension, we simulate a chain of lattice sites with open boundary conditions, with bosons initialized at every other lattice site.Here, we show the chain with 6 lattice sites that was simulated on IBM's quantum processor.With no disorder ( = 0), the bosons are delocalized, but with the onset of any finite value of disorder, the bosons will localize.(b) In two dimensions, we simulate a square lattice with periodic boundary conditions, with bosons initialized in a checkerboard pattern.Here, we show the 16-site lattice that was simulated with the noisy quantum circuit emulator.Up to a disorder strength of  ≈ 9.6, the bosons are delocalized, but at higher disorder strength, the bosons localize.

10 FIG. 2 .
FIG.2.Exact diagonalization results for the one-dimensional case of Eq.(1).The data include  samples = 2 17 random disorder configurations.(a) Average spectral gap Δ between the ground state and the first excited state in the half-filling sector as a function of the system size  for different disorder strengths .We find that Δ ∼  −1 .(b) Cumulative distribution function of the spectral gap Δ for various disorder strengths  at fixed system size  = 8.

9 FIG. 3 .
FIG. 3. Probability distribution for  = 6 of the local density of particles of Eq. (3) for various disorder strengths.(a) Exact diagonalization data computed over  samples = 2 13 random disorder realizations.(b) Experimental data from the quantum computer over a hundred random disorder realizations for each disorder strength with  shots = 2 13 measurements for each realization.

FIG. 4 .
FIG.4.First row: Ideal benchmark based on exact diagonalization for  = 64 averaged over  samples = 2 13 random realizations.Data points close to the boundaries are dismissed.Error bars smaller than the symbols and not shown.Second row: Perfect simulation based on exact diagonalization for  = 6 averaged over  samples = 2 13 random realizations.Third row: Quantum computer result for  = 6 averaged over a hundred random configurations with  shots = 213 measurements each.The legend is the same as in the second row.First column: Average correlator versus the distance for different disorder strengths .Second column: Average correlator rescaled by the disorder-free data versus the distance for different disorder strengths .Third column: Average rescaled correlator versus the rescaled distance using the exponent  minimizing the chi-square in the fourth column.Fourth column: Chi-square  2 characterizing the quality of the data collapse as a function of the correlation length exponent  (smaller is better).The vertical dashed line is positioned at the minimum  2 value for which  leads to the best data collapse ( = 2 is the expected theoretical value).

6 FIG. 5 .
FIG. 5. (a) Cumulative distribution function of the spectral gap Δ between the ground state and first excited state in the half-filling sector.We consider a fixed system size  = 4 × 4 at half-filling for various disorder strengths .Data obtained by exact diagonalization by averaging over  samples = 2 17 random disordered samples.(b) Probability distribution for  = 4×4 of the local density of particles of Eq. (3) for various disorder strengths.Data computed from quantum circuit for adiabatic state preparation over a hundred random disorder realizations.(c) Average coherent density of bosons Λ =0 of Eq. (19) for  = 4 × 4 as a function of the disorder strength .Each data point is averaged over a hundred random realizations whose ground state is obtained by an adiabatic state preparation circuit.From quantum Monte Carlo simulations, we expect the transition between the superfluid and insulating phases to take place at  c ≈ 9.6[60,61].

5 FIG. 6 .
FIG.6.Quantity of Eq. (21) based on the local density of particles at fixed disorder strength  for different noise strengths  and system sizes .(a) One-dimensional data for  = 5 based on  samples = 2 10 independent disordered samples and  shots = 2 10 random noise measurements for each.Error bars are smaller than the symbols and not shown.For  1, we observe a behavior of the form ∼ exp(− ε ) with ε = 15.7(2) a fitting parameter.Here,  has to be interpreted as the depth  of the circuit, which scales linearly with  after the free-fermionic compression (b) Two-dimensional data for  = 5 based on a hundred independent disordered samples and  shots = 10 2 random noise measurements for each.For  → 0, we observe an exponential decay ∼ exp(− ε ) with ε = 30(2) × 10 3 a fitting parameter.

FIG. 7 .
FIG.7.Coherent density based on the momentum distribution function of Eq. (19) at  = 0. Data normalized by the noiseless case with  = 0. Left column: One-dimensional where each data point is averaged over  samples = 2 10 random disordered samples, and each of those obtained from  shots = 2 10 measurements from noisy simulations.Error bars are smaller than the symbols and not shown.Right column: Two-dimensional case case where each data point is averaged over a hundred random disordered samples, and each of those obtained from  shots = 10 2 measurements from noisy simulations.(a) Data at fixed disorder strength  = 0.For 2  1, we observe a behavior of the form ∼ exp(− ε  2 ) with ε = 3.28(5) a fitting parameter.Here, one of the factors  has to be interpreted as the depth  of the circuit, which scales linearly with  after the free-fermionic compression (b) Data at fixed disorder strength  = 5.For  1 we observe a behavior of the form ∼ exp(− ε ) with ε = 24.5(3) a fitting parameter.Here,  has to be interpreted as the depth  of the circuit, which scales linearly with  after the freefermionic compression.(c) Data at fixed disorder strength  = 5 in the delocalized phase (d) Data at fixed disorder strength  = 15 in the localized phase.In both cases, for  → 0, we observe an exponential decay ∼ exp(− ε ) with ε = 60(3) × 10 3 ( = 5) and ε = 75(2) × 10 3 ( = 15), a fitting parameter.