Robust and Ultrafast State Preparation by Ramping Artificial Gauge Potentials

The implementation of static artificial magnetic fields in ultracold atomic systems has become a powerful tool, e.g. for simulating quantum-Hall physics with charge-neutral atoms. Taking an interacting bosonic flux ladder as a minimal model, we investigate protocols for adiabatic state preparation via magnetic flux ramps. Considering the fact that it is actually the artificial vector potential (in the form of Peierls phases) that can be experimentally engineered in optical lattices, rather than the magnetic field, we find that the time required for adiabatic state preparation dramatically depends on which pattern of Peierls phases is used. This can be understood intuitively by noting that different patterns of time-dependent Peierls phases that all give rise to the same magnetic field ramp, generally lead to different artificial electric fields during the ramp. Remarkably, we find that an optimal choice allows for preparing the ground state almost instantaneously. We relate this observation to shortcuts to adiabaticity via counterdiabatic driving. Our findings open new possibilities for robust state preparation in atomic quantum simulators.

Abstract. The implementation of static artificial magnetic fields in ultracold atomic systems has become a powerful tool, e.g. for simulating quantum-Hall physics with charge-neutral atoms. Taking an interacting bosonic flux ladder as a minimal model, we investigate protocols for adiabatic state preparation via magnetic flux ramps. Considering the fact that it is actually the artificial vector potential (in the form of Peierls phases) that can be experimentally engineered in optical lattices, rather than the magnetic field, we find that the time required for adiabatic state preparation dramatically depends on which pattern of Peierls phases is used. This can be understood intuitively by noting that different patterns of time-dependent Peierls phases that all give rise to the same magnetic field ramp, generally lead to different artificial electric fields during the ramp. As an intriguing result, we find that an optimal choice allows for preparing the

Introduction
The engineering of artificial magnetic fields for charge-neutral atoms in optical lattices has been a powerful tool to simulate lattice models with exotic phases including quantum Hall states and topological insulators [1][2][3][4][5][6][7]. More precisely, in these experiments a static artificial gauge potential (in the form of Peierls phases) is engineered in a particular choice of gauge (relative to the plain lattice without magnetic field). Typically, this choice is made based on experimental convenience. For a dynamic process, however, where these artificial gauge potentials are varied in time, this choice does not simply correspond to a gauge freedom anymore. This is because their temporal change generates an artificial electric field. After initial confirmation in a trapped quantum gas [8], such artificial electric forces were observed also in optical lattices [9,10] and predicted to lead to 'gauge-dependent' time-of-flight images of Bose Einstein condensates [11][12][13]. More recently, theoretical investigations showed that the engineering of time-dependent artificial gauge potentials can be employed for quantized charge pumping along tailored paths in two dimensional (fractional) Chern insulators [14,15] and for determining the dynamics of a wave packet in synthetic dimensions [16] and nonlinear systems [17]. With the recent advances in quantum gas microscope techniques [18][19][20][21][22][23][24][25], it becomes more and more important to explore the possibilities of controlling artificial gauge potentials in both space and time. In this paper, we show that this technique can be exploited for the optimization of adiabatic state preparation. Robust adiabatic state preparation is a prerequisite for the experimental investigation (quantum simulation) of interesting states of matter with atomic quantum gases.
As minimal lattice systems with artificial magnetic fields, flux ladders have recently drawn tremendous attention, including the experimental observation of chiral edge currents [25][26][27][28][29][30], the theoretical exploration of rich phase diagrams , the investigation of Laughlin-like states [55][56][57][58][59][60], the study of Hall effect [61][62][63][64][65] and other aspects [66][67][68][69][70][71][72]. In this work, we investigate the adiabatic preparation of the ground state in such ladder systems via continuously ramping up the corresponding Peierls phases. Comparing results for different patterns of Peierls phases, all giving rise to the same magnetic flux, we find that the degree of adiabaticity dramatically depends on this choice. As an intriguing result, the optimal choice of Peierls phases allows for an almost instantaneous preparation of the ground state. We show that for vanishing interactions, this effect can be related to counterdiabatic driving [73][74][75][76][77][78]. However, remarkably our approach works also for very strong interactions, where a simple explanation in terms of counterdiabatic driving is not possible.

Model
We consider interacting bosons in a two-leg ladder described by the Bose Hubbard with bosonic creation operatorâ † ℓ and number operatorn ℓ =â † ℓâ ℓ on site ℓ. The nearest-neighbor tunneling amplitude J ℓ ′ ℓ equals J along legs and J ⊥ along rungs, and it is accompanied by the Peierls phase θ ℓ ′ ℓ . U is the on-site repulsive interaction (a) Figure 1: Bose-Hubbard ladder, with interaction parameter U , tunneling amplitudes J (J ⊥ ) along the legs (rungs) as well as Peierls phases θ ℓℓ ′ either along rungs (a) or legs (b). θ ℓℓ ′ are symbolized by arrows and describe a uniform plaquette flux φ. (c) Phase diagram for non-interacting system. Upper inset shows the lowest Bloch band with single minimum in the Meissner phase, for J ⊥ = 2 and Peierls phases θ ℓ ′ ℓ (φ = π/2, η) with η = 0 (solid orange line) and η = π/2 (dashed green line). Lower inset shows double minima of the lowest band in the vortex phase with φ = 4π/5 and η = 0. The horizontal arrows indicate the paths for our state preparation via ramping artificial magnetic flux. energy. In the following, we use J, /J and lattice constant a as units for energy, time and lengths, respectively.
Due to the complex tunneling matrix elements, the accumulated net phase around one lattice plaquette is analogous to the Aharonov-Bohm phase experienced by a charged particle in a real magnetic field. Thus the Peierls phase θ ℓ ′ ℓ plays the role of a vector potential, and each set of time-independent Peierls phases {θ ℓ ′ ℓ } that gives the same plaquette flux reflects a gauge choice. A uniform flux φ can be realized, for instance, by using gauge potentials along rungs, θ ⊥ ℓ ′ ℓ (φ) [ Fig. 1(a)], or along legs, θ ℓ ′ ℓ (φ, η) [ Fig. 1(b)], with the phase η describing a continuous family of Peierls phases. However, when φ and η vary in time, θ ⊥ ℓ ′ ℓ (φ) and θ ℓ ′ ℓ (φ, η) no longer describe gauge choices, but different artificial electric fields.

Non-interacting case
Let us start with the non-interacting limit (U = 0), for which the phase diagram is shown in Fig. 1(c). For weak magnetic flux, the dispersion relation of the lowest band possesses a unique minimum and the ground state exhibits currents along the leg, resembling the screening currents of the Meissner phase (MP) of a superconductor. Increasing the flux beyond the phase boundary defined by J ⊥ = 2 sin(φ/2) tan(φ/2), the minimum of the dispersion relation splits into two minima and rung-currents appear in the ground state allowing the formation of vortices analoguous to the vortex phase of a type-II superconductor [30,36].
In order to study adiabatic state preparation, we take our initial state and target state as the ground states of the Hamiltonian with flux φ = 0 and φ = π/2, denoted as |ψ 0 and ψ π/2 , respectively. The tunneling amplitude along rungs is fixed at The darker color indicate higher densities and the size of the orange arrows along the bonds is proportional to the amplitude of probability currents. The dashed arrows F indicate directions of the average artificial electric forces. so that the target state lies in the MP, as is marked by the star in Fig. 1(c). By linearly ramping the Peierls phases from zero to final values given by either θ ⊥ ℓ ′ ℓ (φ) or θ ℓ ′ ℓ (φ, η), the flux is continuously increased from 0 to π/2 within the ramping time τ . The evolved state |ψ(τ ) is obtained by numerically solving the Schrödinger equation of the Hamiltonian for a finite system with M = 24 rungs under open boundary condition.
To quantify the degree of adiabaticity, we define the fidelity as the squared overlap between the evolved state and the target state, F = ψ π/2 |ψ(τ ) 2 . Fig. 2(a) shows the fidelities calculated by choosing artificial gauge potentials θ ⊥ ℓ ′ ℓ (φ) and θ ℓ ′ ℓ (φ, η) with η = {0, π/4, 3π/4} [cf. legend in Fig. 2(c)]. For gauge potentials on the rungs, we find fidelities close to 1 for ramping times on the order of τ = 300. For gauge potentials on the legs, this time scale strongly depends on η. Remarkably, it vanishes in the limit of η = 0, so that the ground state can be prepared by switching on the gauge potentials abruptly. This picture is confirmed also by looking at two other quantities characterizing the evolved state. One is the chiral current j c (τ ) scaled by its target value j target c [ Fig. 2(b)], which can be readily measured in experiment [28][29][30] and which plays a key role in charactering different phases in a ladder system [30,36,37,39,40,79]. The other is the excitation energy ∆E [ Fig. 2(c)], defined as ∆E = | ψ(τ )|Ĥ|ψ(τ ) | − E g , where E g is the ground state energy for the final Hamiltonian. Both measures reflect the degree of adiabaticity observed in the fidelity.
The ultrafast adiabatic state preparation can be explained by the fact that the ground state does not depend on the flux for the choice θ ℓ ′ ℓ (φ, η = 0). For the translationally invariant ladder, the single-particle Hamiltonian for quasimomentum where the vector of Pauli matrices σ acts on the sublattice degree of freedom given by the upper and lower leg. The Bloch states |ψ ± (k; η, φ) of both bands E ± (k) = h 0 (k) ± |h(k)| are described by k dependent vectors ±h(k)/|h(k)| on the Bloch sphere. In the MP the ground state lies to quantify the similarity between lowest-band eigenstates with and without magnetic flux φ. Remarkably, in the case of η = 0, the ground state wave function (k = 0) does not depend on the magnetic flux φ, . For a system of M rungs with periodic boundary condition, the quasimomentum k takes discrete values given by integer multiples of 2π/M . As the spectrum is shifted by η, the squared overlap O (π/2) = ψ π/2 |ψ 0 2 between the initial and the target states drops suddenly from 1 to 0 when the shift η becomes larger than π/M , as shown by the dashed line in Fig. 2(e). Since k is not a good quantum number anymore in the finite system with open boundary conditions, we observe a smooth decay of O(π/2) as a function of η, starting from a value close to 1 for η = 0 [O(π/2) = 0.995 for M = 24 rungs]. This behaviour explains that the minimal ramping time τ * required to reach F = 0.9 approaches zero when η drops to zero.

Comparison with counterdiabatic driving
The idea of choosing an optimal vector potential for adiabatic state preparation can be related to the concept of counterdiabatic driving [73][74][75][76][77][78]. To be general, let H p be a Hamiltonian depending on a parameter p and |ψ p the corresponding ground state. Starting from the ground state at p = 0, we wish to rapidly prepare the ground state of the target Hamiltonian H p=f . The idea of counterdiabatic driving is to consider a family of unitaries U p labelled by p, so that the evolved state exactly follows |ψ p(t) for a new Hamiltonian where the second term corresponds to the so-called counterdiabatic driving that could be realized via some external forces [73,74]. Our approach, in turn, corresponds to directly working in the rotated frame of reference with instantaneous eigenstate |ψ ′ p = U † p |ψ p governed by the Hamiltonian For an ideal choice of U p (e.g. the optimal choice of Peierls phases with η = 0 in this work), one can find a p-independent ground state |ψ ′ p = |ψ 0 , and thus it allows for the parameter ramp within arbitrarily short time. The advantage of our approach is that there will be no need for applying external terms to the system. Meanwhile, our protocol can be easily extended to the manybody system, as will be demonstrated in Sections 5 and 6.
The optimal choice (η = 0) of Peierls phases can also be understood intuitively by noting that the artificial electric fields generated during the ramp correspond to Faraday's law of induction, as portrayed in Fig. 2(i). In turn, for non-optimal choices with η = 0, additional electric fields are created as well during the ramp [causing the drift shown in Figs. 2(g,h)] that are not related to the time-dependence of the magnetic field via Faraday's law of induction. These non-Faraday electric fields could be compensated by a time-dependent scalar potential. This freedom of choosing η is a consequence of the fact that the experimentalist directly engineers the artificial gauge potential (via Peierls phases) rather than the artificial magnetic field. The counterdiabatic driving terms required for rapid state preparation for the non-optimal choices of Peierls phases would simply correspond to time-dependent scalar potentials subtracting the non-Faraday forces generated by η(t) = 0. Note, however, that the choice of η = 0 and the absence of non-Faraday forces is not always optimal, as will be seen in Section 6 when discussing parameter ramps leaving the MP.

Role of interactions
Now we simulate the interacting system at filling n = 1/2 per site by using the TeNPy library [80][81][82][83] and a matrix product operator based time evolution method (tMPO) [84,85]. The ground state overlap O(π/2) as a function of interaction strength U is plotted in Fig. 3(a). In the case of η = 0, the overlap O(π/2) exhibits non-monotonous behavior, reflecting a complex competition between manybody interactions and artificial magnetic flux. While the system features a Meissnerlike superfluid ground state for weak interactions [86], (in the thermodynamic limit) it undergoes a Berezinskii-Kosterlitz-Thouless (BKT) transition to a Mott-insulator state with single particles localised on the rungs as U is increased [33,39,87]. The critical parameter is found to be U c1 ≈ 4.2 for φ = π/2 and U c2 ≈ 10.4 for φ = 0 [86], which determines three regions (I: U < U c1 , II: U c1 < U < U c2 , and III: U c2 < U ) shown in Fig. 3(a), where we plot the overlap O(π/2) (blue dots connected by dashed line). In the weakly interacting region I, the overlap first decreases rapidly, before it slightly increases again. This behaviour is qualitatively reproduced by Bogoliubov theory (red dashed line) [86]. It can be related to the fact that the interactioninduced population of finite momentum modes initially happens much faster in the presence of magnetic flux (giving rise to an enlarged effective mass). However, for even stronger interactions the resulting momentum mismatch becomes smaller again [86]. For U c1 < U < U c2 , while the ground state with zero flux remains superfluid, the ground state with flux φ = π/2 already becomes a Mott insulator [86], and therefore the overlap decreases once more. After U > U c2 , the fact that both ground states present Mott-insulating phase gives rise to an increase again. Despite this nonmonotonous behavior, O(π/2) takes comparably large values for η = 0. This leads to rather short adiabatic preparation times also in the strongly interacting regime. In Figs. 3(d) and (e), we plot the fidelity F versus the ramping time for U = 5 and the hard-core limit U → ∞, respectively. Remarkably, for hard-core bosons (and η = 0), we find fidelities close to one already for very short ramping times on the order of 1 (in units of the tunnelling time). The fact that such rapid state preparation found for the strongly interacting system cannot be explained by the single-particle analysis presented in the previous sections. This short ramping time may be related to the fact that the larger overlap of quasimomentum distribution occurs for stronger interaction, as indicated in Figs 3(b) and (c).
In Figs. 3(f) and (g), we investigate the finite size effect for the optimal choice of η = 0. The fact that a many-body fidelity drops with the system size can be attributed to two effects. On the one hand, a finite-size gap (as present in the superfluid regime) separating the ground state from the first excited state decreases with ladder length M , leading to a reduction of adiabaticity. On the other hand, a decrease of the manybody fidelity with system size is expected already from the very fact that (at least for product states) the N -particle fidelity is given by the N th power of the singleparticle fidelity. In order to compensate for the latter effect, when comparing results for different system sizes, we use the single-particle fidelity F 1/N . In Fig. 3(f), we plot F 1/N versus the ramping time τ for η = 0 and U = 20 for different ladder length M . In Fig. 3(g), we extract the ramping time τ ′ above which a fidelity F 1/N ≥ 0.995 is achieved and plot it versus M for various interaction strengths U . Remarkably, we find that the ramping time increases very slowly both in the superfluid and the Mott-insulating regions (i.e. I and III). A noticeable increase is only visible in regime II, where the Mott transition occurs during the ramp.
For finite values of η, taking η = π/4, 3π/4 as examples shown in Fig. 3(a), O(π/2) takes small values until deep in the Mott regime, where the correlations between individual rungs are suppressed by interactions for both φ = 0 and φ = π/2. This can also be understood from the quasi-momentum distribution defined by Fig. 3(b,c) we can see that the distribution is centered around k = 0 for the initial state (φ = 0), and at k = −η for the target state (φ = π/2). Although the shift of quasi-momentum (for η = 0) causes difficulties in state preparations, the increase of interaction broadens the quasimomentum distributions, which results in gradually increasing overlap and a shorter adiabatic ramping time as indicated in Fig. 3(d,e).

Leaving the Meissner regime
So far, we considered parameter ramps within the MP. Increasing φ further gives rise to various phases [37, 39, 40, 44], including the biased ladder phase (BLP) in the weakly and intermediately interacting regime [37, 39-43], which is characterized by vanishing rung currents and the spontaneous Z 2 reflection symmetry breaking in the form of a density imbalance between both legs. In the following, we show that starting from the MP, the BLP can be efficiently prepared by choosing proper Peierls phase patterns (determined by η). Let us start with the non-interaction limit, where beyond a critical flux φ c , the system enters the vortex phase and the dispersion relation develops two degenerate minima. Since each minimum predominantly corresponds to the occupation of one of the legs, the degeneracy can be lifted by introducing a small bias potential (0.01J) between both legs, so that the ground state resembles that of the BLP. Despite the fact that the small bias softens the sharp transition at φ c ≈ 0.667π into a narrow crossover, we observe a sudden drop of the fidelity at φ c when linearly ramping up the Peierls phases with η = 0 [ Fig. 4(c)]. Here the dashed line represents the fidelity between the evolved state and the instantaneous eigenstate. As a remedy, one can vary η during the ramp in such a fashion that the overlap O(φ) remains maximal during the ramp. (For an infinitely large system without bias, this can be achieved by choosingη(t) = arccos J 2 ⊥ /4 cot 2 (φ/2) + cos 2 (φ/2) for φ > φ c , so that the right minimum of the dispersion relation always remains at k = 0 [ Fig. 4(b)].) In this case, the evolved state successfully follows the instantaneous eigenstate even after the critical point, as indicated by the horizontal blue line in Fig. 4(c). Thus, different from the previously discussed case, now the optimal choice of Peierls phases does not correspond to the situation where all the non-Faraday forces were absent during the ramp. Instead, the forces induced by η = 0 are actively employed for state preparation, as they induce shifts in quasimomentum that keep the system state at the minimum of the dispersion relation.
The scheme can also be applied to the interacting system. For instance, the transition to the BLP occurs at critical flux φ ′ c ≈ 0.8π for a 0.8-filling ladder at U = 2.0, J ⊥ = 3 [39]. As shown by the dashed lines in Fig. 4(d), using η = 0 leads to an essentially vanishing fidelity after the critical point, as the BLP has imbalanced distribution between positive and negative quasimomenta due to the broken reflection symmetry [40,41]. To compensate the quasimomentum differences between initial and final states during the ramp, the protocolη ′ (t) [shown in the inset of Fig. 4

(d)]
can be determined by maximizing the ground state overlap, and the corresponding F assumes rather large values as shown by solid lines in Fig. 4(d). Note that the finite value F = 0.78 found for τ = 100 indicates a near unity single-particle fidelity (0.78 ≈ 0.994 N ) for the system with number of particle N = 40 considered here. Higher fidelities can be achieved for longer ramping times.

Conclusion and Outlook
We have proposed to design the time-dependent artificial vector potentials in the form of Peierls phases for rapid adiabatic state preparation in optical lattice systems. Our approach is based on the fact that in such systems the experimentalist directly controls the vector potential rather than magnetic fields. We demonstrated that for a ladder with flux, this approach allows for an almost immediate state preparation for noninteracting bosons. Remarkably, we find very short ramping times also for strongly interacting bosons.
While the abrupt adiabatic preparation in the ladder is an extreme example, it highlights that tunning Peierls phases can be a very powerful tool for state preparation. Specifically, choosing optimal gauge potentials to maximize the overlap between the instantaneous eigenstate and the initial state helps to reduce adiabatic ramping time. It is an interesting open question for future research in how far this approach can be used for the preparation of strongly correlated states of matter, such as fractional Chern insulators.

Dynamics during the ramp
Each point in Fig. 2(a-c) in main text corresponds to the result at the end of a parameter ramp. To interpret the oscillation behavior, we plot the fidelity F = | ψ φ |ψ(t) | 2 and center of mass X during a single ramping process in Fig. 5. It shows that while the center of mass gets closer to the middle of the ladder, the fidelity always has a large value. Thus the oscillation of F is related to the Bloch oscillations of the atomic cloud, which are triggered by the non-Faraday synthetic electric fields that are generated during the ramp for non-zero η.

Meissner-like phases
The ground state chiral current can be used to characterize different phases in a ladder system, like the Meissner or vortex phase. Based on the continuity relation, the local current operators on legs and rungs are respectively defined as [36,39,72,79], probability currents exist only along the legs and behave like screening currents, thus the low-flux phase is identified as a Meissner phase, in analogy to that in a type-II superconductor. For large values of the flux, the system enters into a vortex phase, where finite rung currents emerge and form vortex structures. From Fig. 6 we can see that for J ⊥ /J = 2, φ = π/2, the system with finite size assumes a Meissner-like phase for various values of U .

Bogoliubov Theory
The Hamiltonian can be written aŝ Hereâ † 1,r (â 1,r ) andâ † 2,r (â 2,r ) are the creation (annihilation) operators on the rung r in the lower and upper leg respectively, J denotes the amplitude of nearest-neighbor tunneling along the legs, with θ 1,2 = −η ± φ/2 being the corresponding Peierls phases, so that the flux in each plaquette is φ and we consider φ = π/2 here.
For a two-leg ladder with M rungs, under periodic boundary conditions along the legs, the quasimomentum takes discrete value k = 2π Ma m with m = 0, ±1, ±2, · · · , ±M/2 and a being the lattice constant. By performing the Fourier transformationâ l,r = 1 √ M k e ikarâ l,k , l = 1, 2 the above Hamiltonians can be expressed in quasi-momentum representation aŝ with and periodic Kronecker symbolδ k,q vanishing unless k = q modulo reciprocal lattice constants 2π/a.

Diagonal basis
The single-particle Hamiltonian (8) can be diagonalized by choosing a different basis, i.e.

Bogoliubov approximation
For weak interactions and at low temperature, the number N 0 of particles occupying the single-particle ground state with quasi momentum k 0 remains of the order of total particle number N in a system of finite extent. Thus one can make the approximation which leads tob Keeping all the terms up to second order inb k =k0 , the Hamiltonian (18) becomeŝ with the coefficients where we have used v 2 k0 = 1/2 = u 2 k0 according to Eqs. (10), (11) and (17).
kb k and keeping the terms up to second order in b k , we arrive atĤ with Here we have introduced the total particle number per site n = N 2M , and the additional term − k>0 C −k comes from the commutation relationb † −kb −k =b −kb † −k − 1.

Diagonalization
To diagonalize the Hamiltonian (26), we perform the Bogoliubov transformation with quasiparticle annihilation (creation) operatorsρ k (ρ † k ). Requiring bosonic commutation relations for the quasiparticle operators, we have To get the expressions for µ, ν, we plug Eq. (30) into Eq. (26) and impose that Thus we have which leads to the solutions: 8.3.5. Bogoliubov ground state In the following, we follow Ref. [88] and construct the Bogoliubov ground state Ψ B 0 , which is defined as the state with no quasi-particle, i.e.ρ As the Bogoliubov transformation (30) connects the states with k and −k, the Bogoliubov ground state can be expressed as the states where n k particles are present in k states and n −k particles are in the −k states [88], i.e.
where |0 denotes the vacuum state. Substituting Eq. (38) into Eq. (37) and using the expression ofρ k = µb k − νb † −k according to Eq. (30), we have where we define C k n k ,−1 = 0. Since the basis {|n k , n −k } are orthogonal, we get with α k = −ν/µ for short. By setting n −k = 0 in the above equation Eq. (40), we have C k n k +1,0 = 0 (n k ≥ 0). The similar procedure forb 1,−k Ψ B 0 = 0 gives us C k 0,n −k +1 = 0 (n −k ≥ 0). Based on these observations, it turns out that all the 'off-diagonal' components vanish, i.e. C k n k +1,n −k = 0 (n −k = n k + 1). In the case of n −k = n k + 1, Eq. (40) gives us the following expression of the diagonal terms where C k 0,0 is determined from the normalization of the wave-function. Therefore, the Bogoliubov ground state is a state where pairs of particles with wave vector k and −k are excited.
We denote |n 1, n 2 , · · · as a state with n pairs of particles with non-zero quasimomentum k and −k, and |ψ 0 as the state with k = 0. In this case the Bogoliubov ground state takes the following form where Z = k>0 1 − α 2 k is the normalization factor.
The state |ψ 0 for k = 0 is a coherent stateb 0 |ψ 0 = ψ 0 |ψ 0 and reads where we have defined the vacuum state |vac for the real particles operatorsb k , i.e.
According to Eq. (42) we have the overlap of two ground states The overlap of coherent states O coh ≡ ψ ′ 0 |ψ 0 is obtained by using Eq. (43), 8.3.6. Occupation of finite momentum states In the Bogoliubov ground state Ψ B 0 , pairs of bosons are virtually excited to state with k and −k. The average number of virtually excited bosons with wave vector k is obtained from the Bogoliubov transformation (30) and the definition of Bogoliubov ground state (37), We denote N k =0 as the number of virtually excited particles, i.e. the number of particles in the state |k = 0 , 8.3.7. Results Now we apply the above expressions in our ladder system at 1/2 filling with J ⊥ = 2. We plot the analytic result for the overlap Eq. (44) for M -rung ladder with periodic boundary condition in Fig. 7(a), which shows qualitative agreement with the dip behavior in the weakly interacting regime from the DMRG simulations of finite system with open boundary conditions [ Fig. 7(c)]. Note that the DMRG results for the interacting regime have been divided into three regions. The beginning and the end of the grey shaded region are given by the BKT transition from a superfluid to a Mott insulator for φ = π/2 and φ = 0, respectively. By extracting from the finite-size scaling of peaks in quasimomentum distribution n max k M −3/4 [33,39,87], the crossing determines the BKT-transition points at U c1 ≈ 10.4 for φ = 0 and U c2 ≈ 4.2 for φ = π/2 [ Fig. 7(e,f)]. Overall, both the analytic and numerical results show that the overlaps decay exponentially with the system size for finite U , and approach 1 for the non-interacting case [ Fig. 7(b,d)].  To understand the dip in the weakly interacting regime, we plot the average number of particles with non-zero quasi momentum N k =0 according to Eq. (47), and the relative difference in the occupation of non-zero k-modes ∆n k = n k (π/2)−n k (0) n k (π/2)+n k (0) between φ = π/2 and φ = 0 in Fig. 8(a) and (b), respectively. We can observe that when switching on the interactions, the excited quasi momentum modes become occupied much faster in the presence of magnetic flux. This is related to the fact that the single-particle dispersion relation E − (k) acquires a larger effective mass with increasing flux [see Fig. 4(a) in the main text]. As a result, the momentum modes become occupied rather differently for both fluxes when U is switched on, as can be seen from Fig. 8(b). The slight increase of the overlap for even larger U can then be explained by the fact that the relative differences in the momentum distributions for both fluxes become smaller again.