Superfluid to Mott transition in a Bose-Hubbard ring: Persistent currents and defect formation

TL;DR

This paper investigates the non-equilibrium dynamics of superfluid currents in a Bose-Hubbard ring during the superfluid to Mott insulator transition, revealing persistent oscillations linked to defect formation and spectral gap properties.

Contribution

It introduces a combined variational and tensor network approach to analyze persistent currents and defect dynamics during the transition in a ring-shaped Bose-Hubbard model.

Findings

Persistent superfluid current oscillations inside the Mott phase.

Oscillation amplitude correlates with residual energy and defect creation.

Oscillation frequency matches the spectral gap in the Mott phase.

Abstract

We revisit here the Kibble-Zurek mechanism for superfluid bosons slowly driven across the transition towards the Mott-insulating phase. By means of a combination of the Time-Dependent Variational Principle and a Tree-Tensor Network, we characterize the current flowing during annealing in a ring-shaped one-dimensional Bose-Hubbard model with artificial classical gauge field on up to 32 lattice sites. We find that the superfluid current shows, after an initial decrease, persistent oscillations which survive even when the system is well inside the Mott insulating phase. We demonstrate that the amplitude of such oscillations is connected to the residual energy, characterizing the creation of defects while crossing the quantum critical point, while their frequency matches the spectral gap in the Mott insulating phase. Our predictions can be verified in future atomtronics experiments with neutral atoms in ring shaped traps. We believe that the proposed setup provides an interesting but simple platform to study the non-equilibrium quantum dynamics of persistent currents experimentally.

Equations

$$\hat{H}=-J\sum _{j=1}^L\left(e^{i\varphi /L}{b}_{j+1}^{†}{b}_j+H.c.\right)+\frac{U}{2}\sum _{j=1}^L{n}_j(n_j-1),$$

$$\hat{I}=-\frac{1}{\hslash }\frac{\partial \hat{H}}{\partial \varphi }=\frac{iJ}{\hslash L}\sum _{j=1}^L\left(e^{i\varphi /L}{b}_{j+1}^{†}{b}_j-H.c.\right).$$

$$i\hslash \frac{\mathrm{d}}{\mathrm{d}t}|\psi (t)⟩={\hat{P}}_{|\psi (t)⟩}\hat{H}|\psi (t)⟩,$$

$$|\psi ⟩=\sum _{k,m}{\left(Q^{(\lambda )}{C}^{(\lambda )}(U^{(\lambda )}{)}^{†}\right)}_{k,m}|{\Phi }_{L,k}^{(\lambda )}⟩|{\Phi }_{R,m}^{(\lambda )}⟩,$$

$$\hat{P}={\hat{P}}_M^{(2)}+\sum _{\lambda =1}^{M-1}\left({\hat{P}}_{\lambda }^{(2)}-{\hat{P}}_{\lambda }^{(1)}\right),$$

$${\hat{P}}_{L,\lambda }$$

$${\hat{P}}_{L+C,\lambda }$$

$${\hat{P}}_{R,\lambda }$$

$$i\hslash {\dot{C}}^{(\lambda )}$$

$$i\hslash {\dot{C}}_R^{(\lambda )}$$

$$I(t)=\sum _{\alpha ,{\alpha }^{\prime }}{c}_{{\alpha }^{\prime }}^{\star }{c}_{\alpha }{e}^{-i(E_{\alpha }-E_{{\alpha }^{\prime }})(t-t_0)/\hslash }⟨{\alpha }^{\prime }|\hat{I}|\alpha ⟩.$$

$${ϵ}_{\text{res}}\approx |c_2{|}^2\hslash \Delta .$$

$$\hslash \Delta (U_{\mathrm{f}})\simeq U_{\mathrm{f}}-5.97J+5.20\frac{J^2}{U_{\mathrm{f}}}.$$

$$E_0^{(1)}=J⟨\text{GS}|\hat{V}|\text{GS}⟩=0,$$

$$|s⟩=\frac{1}{\sqrt{L}}\sum _{q=0}^{L-1}(\hat{\mathrm{T}}{)}^q|0\underset{\text{length}s}{\underbrace{\mathrm{11..11}}}\mathrm{21..1}⟩$$

$$⟨s^{\prime }|\hat{V}|s⟩=-3({\delta }_{s,s^{\prime }+1}{e}^{-i\varphi /L}+{\delta }_{s^{\prime },s+1}{e}^{i\varphi /L}).$$

$$E_0^{(2)}=(J^2/U)\sum _{k\ne GS}\frac{|⟨k|\hat{V}|\text{GS}⟩{|}^2}{E_0^{(0)}-E_k^{(0)}},$$

$$E_0^{(2)}=(-4L)(J^2/U)$$

$$E_1^{(2)}=(J^2/U)\sum _{k\notin \{s\}}\frac{|⟨k|\hat{V}|1ex⟩{|}^2}{E_1^{(0)}-E_k^{(0)}},$$

$$E_1^{(2)}=-122.8(J^2/U).$$

$$\hslash \Delta =E_1-E_0\simeq U-5.97J+5.20(J^2/U).$$

$${\hat{I}}_k=\frac{iJ}{\hslash }\left(e^{i\varphi /L}{b}_{k+1}^{†}{b}_k-H.c.\right)$$

Full text

I Introduction

The Kibble-Zurek (KZ) mechanism, introduced to understand the defect formation at a symmetry-breaking phase transition, has been theoretically studied and experimentally verified in many different circumstances Kibble (1980); Laguna and Zurek (1997, 1998); Yates and Zurek (1998); Dziarmaga et al. (1999); Uhlmann et al. (2007); Dziarmaga et al. (2008); del Campo et al. (2010); Witkowska et al. (2011); Das et al. (2012); del Campo and Zurek (2014); Silvi et al. (2016); Weiss et al. (2018); Chuang et al. (1991); Bowick et al. (1994); Sadler et al. (2006); Monaco et al. (2009). The KZ dynamics, initially investigated in classical (finite temperature) critical phenomena has been further extended to quantum critical systems, where it has also notable connections with many-body state preparation and adiabatic quantum computation/quantum annealing Kadowaki and Nishimori (1998); Santoro et al. (2002); Polkovnikov (2005); Zurek et al. (2005); Dziarmaga (2005); Damski (2005); Lamacraft (2007); Pellegrini et al. (2008); del Campo et al. (2012); Anquez et al. (2016). Progress in the realization of quantum simulators has led to important experimental results in KZ physics, most notably with cold atoms Weiler et al. (2008); Lamporesi et al. (2013); Navon et al. (2015); Liu et al. (2018); Keesling et al. (2019) and trapped ions Pyka et al. (2013); Ulm et al. (2013). Despite the multitude of works, the initial proposal put forward more than 30 years ago Zurek (1985) — to observe the formation of defects through the changes of the persistent current in a superfluid ring — was never theoretically or experimentally addressed in quantum systems. The reason is twofold: Theoretically, the simulation of the dynamics of a quantum many-body system on a ring is a highly demanding task and has not been carried out so far. Experimentally, to date it has been difficult to realize a condensate with a ring-shape geometry.

In this paper we fill this gap theoretically, by investigating how persistent currents are modified on crossing the superfluid-insulator transition in a Bose-Hubbard ring. Besides its conceptual interest, we believe that our work is timely in view of the experimental possibilities offered by the newly born field of atomtronics (see, e.g., the focus issue Amico et al. (2017)). We consider the case where a toroidal trapping potential and a lattice modulation along the trapping ring is present, such that the system is described by the Bose-Hubbard model on a ring pierced by an external static gauge field (see Fig. 1). While the ground state of this system has been investigated extensively Cominotti et al. (2014); Arwas et al. (2017), previous studies of the non-equilibrium time-evolution focused on the mean-field regime (e.g. large particle densities), chain geometries without persistent current, or short time analysis Schützhold et al. (2006); Fischer et al. (2008); Cucchietti et al. (2007); Clark and Jaksch (2004); Bernier et al. (2012); Tuchman et al. (2006); Natu et al. (2011); Navez and Schützhold (2010). We present fully quantum, time-dependent results for the persistent current on a ring of up to $L=32$ sites, hence beyond the reach of exact diagonalization. We carry out our calculations using an approach which combines the Time-Dependent Variational Principle (TDVP) Haegeman et al. (2011, 2016) with Tree-Tensor Networks (TTN) Shi et al. (2006); Silvi et al. (2010); Gerster et al. (2014); Montangero (2018); Silvi et al. (2019). We show that the annealing from superfluid to a Mott insulator — which does not carry any current itself at equilibrium — leads to oscillations of the persistent current flowing in the ring, with the amplitude showing a strong dependence on the annealing rate and being related to the residual energy, a key quantity in the KZ mechanism. We also show that the frequency of the persistent current oscillations matches the spectral gap of the Mott insulating final state.

Induced by a classical gauge field, the magnitude of the persistent current of bosonic particles flowing on a ring strongly depends on the interaction between particles. At equilibrium, strong repulsion at integer filling typically leads to Mott localization, suppressing the persistent current, while weak interactions allow current flow in presence of an external field Cominotti et al. (2014). Persistent currents are of particular interest in the growing field of atomtronics Ramanathan et al. (2011); Wright et al. (2013); Labouvie et al. (2015); Haug et al. (2018), where atomic particle currents are used to mimic currents in electronic devices Seaman et al. (2007). For instance, atomtronic quantum interference devices (AQUIDs), being the analogue of superconductor based SQUIDs, allow to study persistent currents in highly controllable systems of ultracold atomic gases Greiner et al. (2002); Trotzky et al. (2012). In these systems of neutral atoms, currents are created by confining a Bose-Einstein condensate to an effectively one-dimensional system and driving the particles through gauge fields, implemented either artificially or by stirring, e.g. using a rotating potential barrier Dalibard et al. (2011); Aidelsburger et al. (2011). Thanks to the unrivaled tunability of interactions and potentials in combination with low decoherence, these systems provide experimental tools to study new collective phenomena, with possible applications in the development of high precision sensors, quantum simulation or quantum information processing Amico et al. (2017), where for instance the superposition of persistent current states could serve as an implementation of qubits Amico et al. (2014); Aghamalyan et al. (2015).

The paper is organized as follows: In Sec. II, we define the time-dependent Bose-Hubbard model considered throughout the paper and discuss some of its equilibrium features. Sec. III is dedicated to our Tree-Tensor Network algorithm. We briefly recapitulate the Time-Dependent Variational Principle and present the algorithm to carry out the time evolution of a TTN. Results are presented in Sec. IV, where we discuss the behavior of the persistent current, its connection to the residual energy, and the system size dependence of our findings. In Sec. V we provide details of our semi-analytical calculations for the oscillation frequency of the current using second order perturbation theory. We finally conclude the paper in Sec. VI with some additional remarks on our results and possible experimental realizations.

II Time-dependent Bose-Hubbard model on a ring

For a system consisting of $L$ sites the Bose-Hubbard model, pierced by a magnetic field, is described by the Hamiltonian

[TABLE]

where $J$ and $U$ are the hopping amplitude and on-site interaction, respectively, and the Peierls’ phase $\phi$ takes into account the flux $\Phi$ through the ring in units of the flux quantum $\Phi_{0}$ ($\phi=2\pi\,\Phi/\Phi_{0}$). From now on, we will simply refer to $\phi$ as the “magnetic field”, for simplicity. We work at fixed density of $\rho=1$ particle per site, where the model exhibits, for $\phi=0$, an equilibrium quantum phase transition between a Mott insulator and a superfluid at $(U/J)_{c}\approx 3.37$ Kühner et al. (2000). Still at equilibrium, but for $\phi\neq 0$, based on previous mean-field studies Oktel et al. (2007) and on a strong coupling analysis of the two-dimensional case Niemeyer et al. (1999), the critical value $(U/J)_{c}$ is expected to decrease as compared to the zero-field case, therefore extending the Mott insulating phase.

We consider an out-of-equilibrium dynamics, in which the initial superfluid in the presence of a magnetic field $\phi\neq 0$ is driven in time across the transition towards the Mott insulating phase by ramping-up the Hubbard interaction $U$ (or, equivalently, by ramping down the hopping matrix element $J$). This procedure is the reverse as compared to the typical Kibble-Zurek scenario discussed in Ref.Zurek (1985), where the transition from the disordered to the ordered phase is considered, corresponding to a Mott to superfluid transition for the Bose-Hubbard model. However, it has been shown in Ref.Gardas et al. (2017) that also the superfluid to Mott transition obeys Kibble-Zurek power law scaling within a limited range of annealing velocities. In contrast to the superfluid-Mott transition, excitations are created only after crossing the transition, meaning that there are almost no excitations within the superfluid phase. The quantity we investigate is the time-dependent current flowing on the ring, defined as the expectation value $I(t)=\langle\psi(t)|\hat{I}|\psi(t)\rangle$ of the current operator

[TABLE]

More in detail, we prepare the system in the ground state $|\psi_{0}\rangle$ of the Hamiltonian $\hat{H}$ within the superfluid phase, specifically for $U(t=0)=U_{\mathrm{i}}=2J$, with a given value of the external magnetic flux $\Phi<\Phi_{0}/2$, corresponding to $\phi=2\pi\Phi/\Phi_{0}$ in the interval $[0,\pi)$ — kept fixed during the dynamics —, and then anneal the value of the interaction $U(t)$ up towards a final value $U(t=t_{0})=U_{\mathrm{f}}$ deep inside the Mott phase, in a time $t_{0}$, at a constant rate $\gamma=(U_{\mathrm{f}}/U_{\mathrm{i}}-1)/t_{0}$. The ramp is followed by a final part of the evolution where the interaction is kept constant at $U_{\mathrm{f}}$. The time-evolution of the interaction is therefore given by (see also Fig. 5, top) $U(t)=U_{\mathrm{i}}(1+\gamma t)$ for $t\leq t_{0}$, while $U(t)=U_{\mathrm{f}}$ for $t>t_{0}$.

III The TTN-TDVP algorithm

In this section we provide details on the time-evolution algorithm, the Time-Dependent Variational Principle (TDVP) applied to a Tree-Tensor Network (TTN). The TTN Ansatz is well-suited for one- and two-dimensional systems with periodic boundary conditions, since the distance between the first and last sites scales with only $\mathcal{O}(\log(L))$ in the Tensor Network Tagliacozzo et al. (2009); Gerster et al. (2014, 2016, 2017). In particular we employ a two-tensor evolution scheme, allowing to adapt the bond dimension during the simulation dynamically. Starting from the TDVP as introduced by Dirac and Frenkel, it has been formulated for loopless Tensor Networks by Haegeman et al. Haegeman et al. (2011). Only recently this method has been further improved Lubich et al. (2015); Haegeman et al. (2016), overcoming problems with small singular values in the original formulation. The algorithm originally has been formulated for Matrix Product States (MPS) and has been extended for MPS with an optimized boson basis Schröder and Chin (2016), introducing additional tensors into the MPS.

In general, the idea of the TDVP algorithm for Tensor Networks is to project the change of the wave function onto the tangent space of the Tensor Network manifold $\mathcal{M}$ with given bond dimension $D$. This guarantees that each update of a tensor is optimal, meaning that the Euclidean distance between the exact evolution and its MPS approximation is minimal Haegeman et al. (2016). Formally the projection is introduced into Schrödinger’s equation by means of the projection operator $\hat{P}_{|\psi(t)\rangle}$, projecting on the tangent space of $\mathcal{M}$ at $|\psi(t)\rangle\in\mathcal{M}$:

[TABLE]

For the model considered in this paper we exploit the U(1) symmetry of the Hamiltonian explicitlySingh et al. (2010, 2011); Silvi et al. (2019) and we use a two-tensor integration scheme in order to allow for a dynamical adjustment of the bond dimension and the dimensions of the different symmetry sectors. For MPS this two-tensor algorithm was presented in Ref. Haegeman et al. (2016). It has been pointed out that TDVP cannot be formulated in the form of a differential equation as in Eq. (3) for the two-tensor scheme, because the projector does not keep the state within the Tensor Network manifold. However, the algorithm requires a discrete time step anyway and therefore we can perform a Singular Value Decomposition (SVD) after each update to bring the Tensor Network back to the variational manifold $\mathcal{M}$. We present this two-tensor integration scheme for a TTN, which can easily be generalized to arbitrary loopless Tensor Networks.

Let us introduce a decomposition of a Tree-Tensor Network consisting of $M+1$ tensors. Since the network is loopless, there are $M$ pairs of tensors, with the tensors of a pair being connected through a bond link. Then, for $\lambda=1,..,M$ labeling such a pair of tensors, as depicted in Fig. 2, the state represented by the network reads

[TABLE]

with $C^{(\lambda)}=C^{(\lambda)}_{L}C^{(\lambda)}_{R}$ being the matrix corresponding to the two-tensor center block: This matrix is obtained by multiplying the matrices corresponding to the tensors of pair $\lambda$, $C^{(\lambda)}_{L}$ and $C^{(\lambda)}_{R}$, where the tensors need be chosen such that in the network geometry $C^{(\lambda)}_{R}$ is closer to the pair $\lambda=M$ than $C^{(\lambda)}_{L}$. Note that for simplicity we use the same notation for tensors and the corresponding matrices formed by fusing tensor indices. The states $|\Phi^{(\lambda)}_{L,k}\rangle$ and $|\Phi^{(\lambda)}_{R,m}\rangle$ in Eq. (4) are product states in the local basis. They correspond to the two subsystems obtained by splitting the network between the matrices $C^{(\lambda)}_{L}$ and $C^{(\lambda)}_{R}$. In Fig. 2 we visualize this decomposition in an example. Note that even though the labels indicate a separation into a left and right part, the physical bipartition usually is of different shape.

Since the network is free of loops we can choose it to be isometrized with respect to the center block $C^{(\lambda)}$, such that $Q^{(\lambda)}$ and $U^{(\lambda)}$ have orthonormal columns $\left((V^{(\lambda)})^{\dagger}V^{(\lambda)}=\mathbb{1},V=Q,U\right)$. Using this decomposition it is possible to follow the derivations in Refs. [Schröder and Chin, 2016; Lubich et al., 2015, 2013], to find the projector onto the space of two-tensor variations:

[TABLE]

with $\hat{P}^{(2)}_{\lambda}=\hat{P}_{L,\lambda}\otimes\hat{P}_{R,\lambda}$, $\hat{P}^{(1)}_{\lambda}=\hat{P}_{L+C,\lambda}\otimes\hat{P}_{R,\lambda}$ and

[TABLE]

One of the key ideas of the TDVP in this form is to use a Lie-Trotter splitting Hairer et al. (2006) for Eq. (3), yielding one differential equation for each of the summands in the projector. These differential equations can be solved efficiently one after the other Lubich et al., 2015. We can therefore update the network according to every summand in the projector step by step. Since there are two types of projectors, $\hat{P}^{(2)}_{\lambda}$ and $\hat{P}^{(1)}_{\lambda}$, we get two structurally different equations to update the network. In particular the projector $\hat{P}^{(2)}_{\lambda}$ introduces the update of a two-tensor block $C^{(\lambda)}$, while the projector $\hat{P}^{(1)}_{\lambda}$ leads to an update of the tensor $C^{(\lambda)}_{R}$. Effectively these updates are governed by the equations

[TABLE]

where the application of the effective Hamiltonian on the right hand side corresponds to a tensor contraction of multiple indices. The effective Hamiltonians are constructed from the summands on the right hand side of Eq. (3) as depicted in Fig. 3 for the simple case of a two-layer tree (8 physical sites), where the center block is built from the left lower and upper tensors. Note that after evolving the two-tensor block an SVD has to be performed in order to bring the network back to the initial form. The full algorithm is summarized in Alg. 1, while some of the steps are shown graphically in Fig.4. In our implementation we chose the pair built from the uppermost tensors to be the one with $\lambda=M$, for which there is no backwards evolution. The tensors in the TTN are then evolved in time from left to right, bottom to top.

IV Results

Since the system starts in the ground state $|\psi_{0}\rangle$ at a finite value of the external flux $\phi\in[0,\pi)$ and inside the superfluid phase, it displays an initial persistent current $I(0)=\langle\psi_{0}|\hat{I}|\psi_{0}\rangle$ (see Fig. 1). Upon annealing the interaction upwards to enter in the Mott region of the equilibrium phase diagram, the current $I(t)$ drops towards zero (see Fig. 5), the expected value of the equilibrium current at zero temperature in the Mott phase. There is however a residual, time-dependent current which is approximately sinusoidal, $I_{\mathrm{AC}}(t)\sim I_{0}\cos(\omega_{0}t+\varphi)$, with a characteristic main frequency $\omega_{0}$. The amplitude $I_{0}$ depends on the annealing rate $\gamma$ and will be the key signal of the defects created by crossing the quantum phase transition: The slower is the annealing, the smaller is $I_{0}$, as suggested by the three curves in Fig. 5. As expected, this extrapolates well to the adiabatic limit ($\gamma\rightarrow 0$), where no oscillations are present, since the equilibrium current of the Mott insulator vanishes. On the other hand, the main frequency $\omega_{0}$ is essentially unaffected by the annealing rate, but increases when the final Mott interaction $U_{\mathrm{f}}$ is increased. As we will show, this behavior is related to the dynamical spectral gap of the final Mott insulating state, as the persistent alternating current is related to the dynamical oscillation between the Mott insulating state and higher excited states in which few holes and multiply occupied sites are created.

In the following, we analyze the current amplitude $I_{0}$ and its relation to the residual energy $\epsilon_{\text{res}}$, defined as the excitation energy above the final ground state energy, and the oscillation frequency $\omega_{0}$ in more detail. At the end of the section we will discuss the size dependence of the results.

We extract the oscillation amplitude by numerically calculating $I_{0}=(1/2)[\max\limits_{t\geq t_{1}}I(t)-\min\limits_{t\geq t_{1}}I(t)]$, where we take $t_{1}>t_{0}$ in order to neglect some possible transient behavior after the final value of $U_{\mathrm{f}}$ is reached. The results are presented in Fig. 6, quantifying the decay of the amplitude with decreasing annealing rate $\gamma$. As shown hereafter, this behavior is understood by the dependence of the amplitude on the occupation of higher excited states, and in particular the occupation $c_{2}$ of the first excited state of the final Hamiltonian. The probability $c_{2}$, on the other hand, behaves like in a two-level system undergoing a Landau-Zener dynamics, with decreasing excitation probability as $\gamma\rightarrow 0$. The amplitude $I_{0}$ is the main signature of the quasi-adiabatic driving through the transition. This is further corroborated by a relation between the current amplitude and the residual energy as detailed below. Denoting by $\{|\alpha\rangle\}$ the (many-body) eigenstates of the final Hamiltonian $\hat{H}_{\mathrm{f}}$, with associated energies $E_{\alpha}$, one can write for any time $t>t_{0}$:

[TABLE]

The constants $c_{\alpha}=\langle\alpha|\psi(t_{0})\rangle$ are the overlaps between the eigenstates of $\hat{H}_{\mathrm{f}}$ and the state $|\psi(t_{0})\rangle$ at the end of the annealing ramp. The diagonal matrix elements of the current operator can lead to a DC-component of the current, while off-diagonal elements can cause oscillations of the current. Let us now assume that the annealing ramp is slow, such that the system mainly remains in the ground state $|1\rangle$ and only one excited state $|2\rangle$, compatible with the translational symmetry, is slightly occupied during the ramp $(|c_{2}|\ll|c_{1}|\Rightarrow c_{1}=1+\mathcal{O}(c_{2}^{2}))$, while all the other excited states have negligible occupation. Then Eq. (7) simplifies and becomes $I(t)\sim|c_{2}\langle 1|\hat{I}|2\rangle|\cos(\Delta\cdot t+\varphi)+I_{\text{offset}}$, where $\varphi$ is an unimportant phase, while $\hbar\Delta=E_{2}-E_{1}$ is the energy gap between ground and first excited state. Here $I_{\text{offset}}$ is a DC-component, dominated by the ground state contribution, and vanishing deep in the Mott phase, as demonstrated in Fig. 1. On the other hand we can argue in the same way that the residual energy will be dominated by the contribution of the first excited state $|2\rangle$:

[TABLE]

In the slow annealing regime, this implies that both the current amplitude ($I_{0}\propto|c_{2}|$) and the residual energy ($\epsilon_{\text{res}}\propto|c_{2}|^{2}$) depend on the annealing rate $\gamma$ only through the occupation $c_{2}(\gamma)$ of the first excited state. In Fig. 6 this relation ($\epsilon_{\text{res}}(\gamma)\propto I_{0}^{2}(\gamma)$) is verified numerically for different $\gamma$, by comparing the residual energy to the square of the current (multiplied by a size dependent prefactor). Clearly, the agreement is quite good, in the slow annealing regime. Theoretically the propotionality prefactor is given by $\epsilon_{\text{res}}/I_{0}^{2}=\hbar\Delta/|\langle 1|\hat{I}|2\rangle|^{2}$, and has been obtained by fitting this ratio for slow annealing ramps to a constant.

As explained in Ref.Gardas et al. (2017) the residual energy shows an effective power-law behavior within a limited range of annealing velocities, which we find for $1/\gamma\lesssim 3$ (L=32) and $1/\gamma\lesssim 1.7$ ($L=16$). However, in this regime the simple picture of having only one excited state involved in the dynamics breaks down, as manifested by the disagreement of the current and the residual energy for $L=16$ and slightly indicated by the latest points for $L=32$. Very fast annealing leads to a rapid increase of entanglement, preventing us from studying this regime for large systems. Moreover, fast annealing leads to the appearance of additional frequencies in the current, with a more quantitative analysis, however, being beyond the scope of this work.

We now concentrate on the oscillation frequency, and in particular on the most relevant frequency, obtained by calculating the Fourier transform of $I(t)$ in an appropriate time window $[t_{1},t_{2}]$ $I(\omega)=\int_{t_{1}}^{t_{2}}\!\mathrm{d}t\,I(t)e^{-i\omega t}\;,$ where again $t_{1}>t_{0}$ as above in order to allow for a steady-state behavior to set in, and $t_{2}$ is the total simulation time. The Fourier transform of the current $I(\omega)$ displays a sharp carrier frequency $\omega_{0}$, where $|I(\omega_{0})|$ is maximum (see inset in Fig. 7). The oscillation frequency $\omega_{0}$ is basically insensitive to the external flux $\phi$ and to the annealing rate $\gamma$. However, it does depend on the final on-site interaction $U_{\mathrm{f}}$: it increases with $U_{\mathrm{f}}$, as shown in Fig. 7. We will now argue that, in the quasi-adiabatic limit, $\omega_{0}$ is essentially related to the gap $\Delta$ between the ground and the first excited state, hence depends only on the final Hamiltonian. A theoretical analysis using strong coupling perturbation theory for the ground and first excited state – presented with more details in Sec. V – yields:

[TABLE]

The corresponding theory curves including up to first and second order terms are plotted in Fig. 7, showing excellent agreement with the numerical results obtained from the Fourier transform of the current, hence confirming that $\hbar\omega_{0}=\Delta$. Furthermore we observe that the created defects are not single holon-doublon pairs only, as in this case zeroth order perturbation theory – taking into account only a single holon-doublon exicitation – would fit the numerical data. The necessity of going to second order shows that the defects consist of multiple holon-doublon pairs, sites with more than two particles, and several empty sites.

Let us now discuss the system size dependence of the results. First it is important to note that the equilibrium current of a perfect superfluid (interaction $U=0$) decreases as $1/L$ with the system size and therefore vanishes in the thermodynamic limit. For the oscillation amplitude of the current we can expect to find a similar behavior as long as the annealing is quasi-adiabatic, being the regime where we can relate the current amplitude to the residual energy. In Fig. 6 this behavior is indicated by the larger prefactor for $L=32$ as compared to $L=16$. On the other hand the energy gap at the critical point decreases with the system size (see Ref.Gardas et al. (2017)), leading to an increased population of excited states and correspondingly larger oscillations of the current. This effect is shown in Fig. 8, where the amplitudes for $N=16$ and $N=32$ are very similar for fixed annealing velocity, even though the initial values differ by a factor of $2$. However, keeping the annealing velocity fixed and increasing the system size will bring the dynamics out of the quasi-adiabatic regime. Therefore, if we stay in the quasi-adiabatic regime, the oscillations will vanish once we increase the system size.

V Strong coupling perturbation theory

In the previous section we compared numerical results for the oscillation frequency of the current with the spectral gap at the end of the evolution as obtained from perturbation theory. In the following we provide some details of our perturbative calculations for the gap $\Delta$ in the strong repulsion limit. To this end we expand the energies of both, the ground state and the first excited state — compatible with the symmetry — up to second order. From here on we focus on the case $L=32$, while the field is fixed to be $\phi=0.7\pi$ as before.

It is important to note that the Bose-Hubbard model is translationally invariant, implying that the initial state of the dynamics – a superfluid ground state – is an eigenstate of the translation operator $\hat{T}$, where we find the corresponding eigenvalue to be 1. Due to the translational invariance of the model this eigenvalue is conserved, i.e. the time-dependent state will be an eigenstate of $\hat{T}$ with unit eigenvalue at any time. Later we will use this to construct a suitable basis.

V.1 Zeroth and first order

In the infinite interaction limit $J/U=0$ the ground state of Eq. (1) is a Mott insulator with one particle per site $|\text{GS}\rangle=|11...11\rangle$, such that to zeroth order the energy for the ground state is $E_{0}^{(0)}=0$. The first excited energy level is degenerate with $E_{1}^{(0)}=U$ for all states with one doublon and one holon.

Since the ground state is non-degenerate, the first order contribution to the ground state energy is given by

[TABLE]

where $\hat{V}=-\sum_{j}\left(e^{i\phi/L}b_{j+1}^{\dagger}b_{j}+H.c.\right)$ is the perturbation. In order to calculate the first order contribution to the excited states we use degenerate perturbation theory. First, we construct a suitable basis, where in particular we restrict the Hilbert space to the subspace of eigenstates of the translation operator $\hat{\mathrm{T}}$ with unit eigenvalue. Let us define the basis states

[TABLE]

with $s=0,..,L-2$ the separation of the doubly occupied and the empty site. In this basis the perturbation matrix elements evaluate to

[TABLE]

The first order correction is given by the smallest eigenvalue of this $(L-1)\times(L-1)$ matrix, and is found numerically to be $E_{1}^{(1)}=-5.97J$ for $L=32$ and $\phi=0.7\pi$. Note that the exact eigenvalues of this tridiagonal Toeplitz matrix are known. However, we restrict ourselves to the numerical values for convenience in the next section.

V.2 Second order

The second order contribution for the ground state is given by

[TABLE]

where $|k\rangle$ are the energy eigenstates for $J/U=0$. In this sum only states $|k\rangle$ with one doubly occupied and one neighboring empty site can give a nonzero contribution. Since the operator $\hat{V}$ contains $2L$ operators contributing $-2$ each, we find for the second order correction of the ground state energy:

[TABLE]

Similarly we continue for the second order contribution of the excited state, which is

[TABLE]

where $|1ex\rangle$ (not to be confused with $|s=1\rangle$) is the first excited state we found by diagonalizing the perturbation matrix. Using the numerical result for this state we can evaluate the sum to obtain the energy correction ($L=32$)

[TABLE]

Adding up all the different contributions we find for the gap:

[TABLE]

VI Discussion and Conclusion

Our predictions are ready for experimental verification in future setups of cold atomic systems on ring lattices, where the Bose-Hubbard model including a gauge field can be naturally implemented with 87Rb atoms. Quasi-local currents can be measured following the read-out scheme presented in Ref. Trotzky et al., 2012. Due to the quasi-adiabatic annealing, the evolution times required for the observation of oscillations need to be longer than in commonly used sudden quench scenarios. However, the time periods in the order of $10\hbar/J$ needed here are still reachable in state-of-the-art experiments [Trotzky et al., 2012]. Moreover the non-decaying oscillations might eventually be used to demonstrate long-living coherence in next generation quantum simulators.

In conclusion, we have investigated the time-dependent behavior of the persistent current following a linear annealing procedure on a Bose-Hubbard ring of up to $L=32$ sites, where in particular we analyzed the crossover from a superfluid to a Mott insulator in the presence of a gauge field. We found that after an initial decay the current starts to oscillate around nearly zero current, being the ground state value at the end of the ramp. The current is nonzero due to the excitation of higher states and, in particular, results from the non-diagonal matrix elements of the current operator. In the slow annealing regime, where only one excited zero momentum state is occupied, the oscillation amplitude is proportional to the square root of the occupation probability and can therefore be related to the residual energy, characterizing the creation of defects. In a closed system the coherent oscillations are not expected to decay even for longer evolution times as long as the annealing is sufficiently slow. Instead, they persist according to Eq. (7), due to the occupation of few eigenstates and correspondingly a small number of frequencies. This can easily be understood from the special case where only the ground and first excited state are populated. After the annealing ($t>t_{0}$) the populations stay constant, while the relative phase between the two populated eigenstates changes in time, producing the oscillations of the current. On the other hand, fast annealing leads to the presence of many frequencies, which might result in the decay of oscillations due to averaging effects. Using perturbation theory up to second order, we have been able to compute the frequency of the oscillations – defined through the final Hamiltonian only – in the limit of strong final interactions and slow annealing, in very good agreement with the numerical findings. While in this work we focused on the case where translational invariance is not broken, it might be interesting for future research to include a localized barrier, breaking the translational symmetry and providing the possibility to realize current-based qubits Amico et al. (2014); Aghamalyan et al. (2015).

While finalizing this manuscript we became aware of the work by Bauernfeind and AichhornBauernfeind and Aichhorn (2019), presenting in detail the Time-dependent variational principle for general loopless Tensor Networks with an application to the Fork Tensor Product States tensor network used for Dynamical Mean-Field Theory calculations.

Acknowledgements

Research was partly supported by EU FP7 under ERC-ULTRADISS, Grant Agreement No. 834402, the Quantum Flagship project PASQUANS, the EU via ERC synergy grant UQUAM and via QuantEra QTFlag, the Austrian Research Promotion Agency (FFG) via QFTE project AutomatiQ, the Italian PRIN 2017, the DFG funded TWITTER project, and the US Air Force Office of Scientific Research (AFOSR) grant FA9550-18-1-0319. RF and GES acknowledge that their research has been conducted within the framework of the Trieste Institute for Theoretical Quantum Technologies (TQT).

Appendix A:

Numerical parameters/ convergence

In this section we discuss the numerical parameters we used for our simulation. For the results shown previously we use the bond dimension $D=60$, the local bosonic dimension $d=5$ — translating to the limitation of at most four particles per site —, and the fixed time step $\Delta t=2\times 10^{-3}\hbar/J$. In the following we focus on the convergence in the bond dimension for the largest system considered ($L=32$), since numerical errors due to the truncation of the local boson occupation and the finite time step were found to be small compared to the error due to the bond dimension. As shown in Fig. 9 the oscillation amplitudes obtained for $D=60$ compare well with $D=50$, while more significant differences are visible in comparison with $D=40$. This indicates that indeed $D=60$ is enough to obtain accurate results. Note that in contrast to equilibrium scenarios, where the ground state energy decreases monotonically with the bond dimension, dynamical quantities like the oscillation amplitude can show non-monotonic behavior as a function of the bond dimension.

In practice it turned out that alternatively we can check if the bond dimension is sufficiently large by comparing the local currents between sites $k$ and $k+1$, defined as the time-dependent expectation value of the operator (compare with Eq. (2))

[TABLE]

Considering that the model of interest is translationally invariant, we expect to find the same local current between any pair of neighboring sites. However, the TTN breaks the translational invariance, resulting in different local currents if the bond dimension is too small. As visualized in Fig. 10, the local currents do not agree very well for $D=40$, while the agreement is much better for $D=60$, confirming that $D=60$ yields reliable results.