Interaction-driven dynamical quantum phase transitions in a strongly correlated bosonic system

We study dynamical quantum phase transitions (DQPTs) in the extended Bose-Hubbard model after a sudden quench of the nearest-neighbor interaction strength. Using the time-dependent density matrix renormalization group, we demonstrate that interaction-driven DQPTs can appear after quenches between two topologically trivial insulating phases -- a phenomenon that has so far only been studied between gapped and gapless phases. These DQPTs occur when the interaction strength crosses a certain threshold value that does not coincide with the equilibrium phase boundaries, which is in contrast to quenches that involve a change of topology. In order to elucidate the nonequilibrium excitations during the time evolution, we define a new set of string and parity order parameters. We find a close connection between DQPTs and these newly defined order parameters for both types of quenches. In the interaction-driven case, the order parameter exhibits a singularity at the time of the DQPT only when the quench parameter is close to the threshold value. Finally, the timescales of DQPTs are scrutinized and different kinds of power laws are revealed for the topological and interaction-driven cases.

We study dynamical quantum phase transitions (DQPTs) in the extended Bose-Hubbard model after a sudden quench of the nearest-neighbor interaction strength. Using the time-dependent density matrix renormalization group, we demonstrate that interaction-driven DQPTs can appear after quenches between two topologically trivial insulating phases-a phenomenon that has so far only been studied between gapped and gapless phases. These DQPTs occur when the interaction strength crosses a certain threshold value that does not coincide with the equilibrium phase boundaries, which is in contrast to quenches that involve a change of topology. In order to elucidate the nonequilibrium excitations during the time evolution, we define a new set of string and parity order parameters. We find a close connection between DQPTs and these newly defined order parameters for both types of quenches. In the interaction-driven case, the order parameter exhibits a singularity at the time of the DQPT only when the quench parameter is close to the threshold value. Finally, the timescales of DQPTs are scrutinized and different kinds of power laws are revealed for the topological and interaction-driven cases.
Introduction. The nonequilibrium dynamics of quantum many-body systems is a vibrant field of research, much less well understood than the phases of matter in thermal equilibrium. One way of bridging this gap is to drive a system away from an initial equilibrium state, for instance, by a sudden quench [1][2][3][4]. This allows one to define new concepts of nonequilibrium criticality, and to study how these are related to the equilibrium cases [5][6][7][8][9][10][11][12]. Experimentally, cold atom systems can routinely implement a quantum quench, where a closed quantum system prepared in its ground state shows a nontrivial unitary evolution after a sudden change of parameters. [13][14][15][16].
In this context, a particularly interesting generalization of equilibrium concepts is that of dynamical quantum phase transitions (DQPTs), where time and the Loschmidt echo play the roles of temperature and the partition function, respectively [10,17]. As the Loschmidt echo is often inaccessible experimentally, it is highly desirable, from a practical viewpoint, to identify observables whose dynamics are linked to those of the Loschmidt echo.
It has been shown that, in general, there is no exact correspondence between DQPTs and the equilibrium quantum phase transitions of a system [18][19][20][21]. On the other hand, various conditions have been identified under which a clear connection exists [17]. Rigorous examples are exactly solvable models that exhibit ground states with distinct topological invariants [22][23][24]. Here a quench across the phase boundaries always leads to DQPTs. Moreover, after starting in symmetry broken phases, the respective order parameter becomes zero at the times of the DQPTs [25][26][27]. Both phenomena have been confirmed in experiments [28][29][30][31][32].
In contrast, the link of DQPTs with order parameter dynamics and equilibrium phase transitions is much more ambiguous for nonintegrable models [33][34][35][36]. Only recently, such correspondences were found for quenches in a nonintegrable spin-1 XXZ chain, with respect to the symmetry protected topological Haldane phase and its nonlocal string order parameter [37,38]. Such a model has also been experimentally realized [39,40]. Further theoretical studies highlight the essential role played by symmetries in the dynamics of the string order parameter [41,42]. It is yet unclear if these correspondences still hold when the symmetry is lowered or the quench does not cross a topological phase boundary.
In order to investigate such a problem, in this work, we study a chain of bosons with on-site and nearestneighbor interactions [43][44][45]. Compared to the spin-1 XXZ chains, the model is less symmetric; for instance, particle-hole symmetry is missing. First, we consider quenches between topologically trivial phases. We find that strong nearest-neighbor interaction quenches above a threshold V dyn c can induce DQPTs even when no topology is involved. We call these "interaction-driven" DQPTs. Near the threshold, the order parameter develops nonanalytic signatures with temporal correspondence to the DQPTs. Second, we investigate quenches across topologically distinct phases starting from a Haldane insulator state. Here, quenches that reduce the interaction strength show no DQPTs on short timescales and no zeros of the order parameter, whereas quenches to stronger interactions do. The former, additionally, lead to a transient reduction of the entanglement entropy. Finally, we show that the time of the first DQPTs depends on the quench parameter roughly in a power law fashion for both kinds of initial states, while the detailed forms differ be-tween the interaction-driven and topological cases.
Model and methods. We consider the one-dimensional extended Bose-Hubbard model (EBHM) with on-site and nearest-neighbor interaction whereâ † i andâ i are bosonic creation and annihilation operators for site i andn i denotes the corresponding number operator. Throughout the paper, We fix the filling factor to N/L = 1 and the hopping strength to J = 1, setting energy and time units.
This model has attracted considerable attention over the past decades, resulting in a good understanding of its equilibrium phases [43,44,[46][47][48][49]. These include a superfluid phase for weak interactions, a Mott insulator (MI) for dominating U , a density wave (DW) phase for dominating V , and the topologically nontrivial Haldane insulator (HI), which occurs for intermediate parameters.
We focus on the case U = 5, where the MI-HI (HI-DW) transition lies at V eq c1 ≈ 2.95 (V eq c2 ≈ 3.525) [47]. We initialize the system in ground states corresponding to the MI (V i = 1.0) and HI (V i = 3.25) phases and subsequently drive it to nonequilibrium by a sudden change of V from V i to V f . The ground states are calculated by the density matrix renormalization group method, and the time-dependent variational principle is employed for the time evolution [50][51][52]. We have checked that the results are independent of the boundary conditions [53]. Also we note that the initial states are neither simple Fock states nor highly symmetric states such as at the Affleck-Kennedy-Lieb-Tasaki (AKLT) or Heisenberg points in the Haldane phase of spin-1 chains [37,41,42,54]. As shown below, the initial condition sensibly affects the DQPTs.
Quantities of interest. One way to generalize the notion of thermodynamic phase transitions to the nonequilibrium case is to consider the rate function as a dynamical analog of the free energy density, where |ψ(t) = e −iĤt |ψ(0) and |ψ(0) is the initial wavefunction of the system. DQPTs can be defined as nonanalytic points of this function in the thermodynamic limit [10]. Second, we study the time evolution of the order parameters to quantify the relation of the time-evolved states to the underlying ground state orders. For the MI and HI phases, these are, respectively, described by the nonlocal string and parity operators [43] O z string (i, j) = δn i i<k<j e iπδn k δn j , where δn i =n i − 1.Ô z string is inspired by the string operators corresponding to the z component of spin-1 XXZ chains, which can be considered as an effective model of the EBHM under the mapping δn i → S z i = −1, 0, 1. Order parameters are defined as the long-distance limits, i.e., O γ = lim |i−j|→∞ Ô γ (i, j) for γ ∈ {string, parity}. Below, we use i = L/4 and j = 3L/4.
We further propose another set of nonlocal operators that correspond to the x components of the spin-1 operator aŝ wherê and P ≤2 k is the projection onto site occupations n k ≤ 2. [55] In principle, the eigenvalues ofÔ k should all be of modulus 1, because otherwise the long-distance limit of expectation values of Eq. (4) will be either 0 or ∞ [56]. However, the above definition always yields a zero order parameter in the presence of occupation numbers n > 2. This problem cannot be mitigated by a straightforward generalization ofÔ k to n > 2, because there are no negative counterparts 2 − n < 0. More details are given in the supplemental material (SM) [57].
A simple solution consists in using the following projection string operator As shown in SM [57],P (i, j) exhibits the same asymptotic decay length ξ as the x component operators of Eq. (4), i.e.,P for γ ∈ {string, parity}. Hence, we obtain well-defined order parameters As the operators in Eq. (4) are off-diagonal in the Fock basis, these order parameters represent a new type of long-range phase coherence.  Results. We now turn to the numerical results. Fig. 1 shows the time evolution of the system initialized in an MI ground state (V i = 1.0) and quenched to V . As depicted in Fig. 1(a), the rate function develops pronounced kinks when V f is above certain threshold values. For the first maximum of the rate function we find V dyn c ≈ 4.2 (see Fig. 4), for the second one V dyn c ≈ 3.6. Both cases represent interaction-driven DQPTs, as the critical interaction strengths to observe DQPTs are not related to the topological MI-HI transition at V eq c1 ≈ 2.95. The former, V dyn c ≈ 4.2, is also far away from the HI-CDW phase boundary at V eq c2 ≈ 3.525. The corresponding dynamics of O z parity andÕ x parity are shown in Figs. 1(b)-(c), respectively. The two parity order parameters remain close to their initial values when V f is in the MI region (V f = 2.0). On the other hand, they decay by several orders of magnitude, if the phase boundary is crossed, i.e., V f ≥ V eq c1 ≈ 2.95. In particular, when V f is just above the threshold value of the DQPTs, V dyn c ≈ 4.2, the parity order parameters show minima at the time of the first kink of the rate function. ForÕ x parity the minimum is especially sharp. After the second kink, when λ drops to lower values again, the order parameters show revivals. For larger V f , there is no longer such a clear temporal relation (see SM [57] for a finer grid around V f = 4.5). As shown in SM [57], we confirmed similar behavior of the parity order parameters in a spin-1 chain. Contrarily, for the DQPTs for a quench starting from the topologically nontrivial HI phase (shown below), the temporal correspondence between the rate function and the order parameters persists for any quench, not just near V dyn Interestingly, a finite-size analysis reveals thatÕ x parity actually vanishes for all V f after a certain time as shown in the inset of Fig. 1(c), where the time is marked by a vertical blue line. This means that the correlation function Ô x parity (i, j) decays more rapidly with |i − j| than the projection string P (i, j) , and the quench destroys the corresponding long-range phase coherence. The time evolution shown in Fig. 1(c), therefore, indicates only shortrange correlations. By contrast, in the spin-1 model, O x parity does not vanish in the thermodynamic limit [57].
In order to confirm that λ develops nonanalytic features at certain critical times t * , we analyze the systemsize dependence of dλ/dt in Fig.2. For V f = 4.0 [ Fig. 2(a)], the first maximum around t ≈ 0.8 is smooth, while the second one converges slowly in L towards a sharp jump. For the stronger quench to V f = 5.0 [ Fig. 2(b)], dλ/dt changes signs at critical times t * /J ≈ 0.55 and 1.05. With increasing L, the regions around the jumps showing strong system-size dependencies diminish. For example, for V f = 5.0 at t * /J ≈ 0.55, there are large spikes in dλ/dt only for L = 54 and 180.
Next, let us turn to the case where the system is initialized in the HI phase (Fig. 3). We always find DQPTs when V is quenched to a larger value across the HI-CDW boundary, as shown in Fig. 3(a). Thus, we have c2 . On the other hand, for quenches across the MI-HI boundary to V f < V eq c1 , depicted in Fig. 3(b), there are no kinks on the given timescale. While it is possible that DQPTs occur at later times, we certainly see that the first few maxima of λ are smooth. This result is in contrast to Ref. [37], which always found kinks at the maxima of λ for a spin-1 XXZ chain initialized to the AKLT state within the Haldane phase. The discrepancy can be traced back to the lower symmetry of our initial state. As shown in SM [57], a spin-1 chain initialized in a less symmetrical state but still in the Haldane phase shows smooth first maxima of λ for corresponding quenches. The dynamics of the renormalized x string order pa-rameterÕ x string are depicted in Fig. 3(c). We find zeros ofÕ x string in close temporal relation to the kinks of the rate function. Similar results were obtained for the spin-1 chains [37]. In the presented time range, the zeros appear only when V is increased, but not when V is lowered. However, a longer time simulation up to t/J = 4.0 for a quench from HI to MI shows a zero ofÕ x string , indicating the existence of DQPTs [57].
As another interesting property of the time-evolving state, we consider the entanglement entropy S about the central bond, which we depict in Fig. 3(d). While it shows the typical approximately linear increase at later times, the short-time behavior depends on whether we lower or increase the interaction V . The entanglement entropy increases monotonically for increased V also at short times. A decreased V , on the other hand, leads to a decrease of entanglement in the system at short times. The weaker the quench, the stronger is the suppression of the entanglement. As shown in the inset of Fig. 3(d), the projections to n ≤ 2 increase when the entanglement entropy decreases, because of the lowered nearest-neighbor repulsion. The entanglement entropy of the ground state is lower in the MI phase along the parameters considered [46]. Thus, we infer that at short times the suppression of higher occupations reduces the entanglement entropy. Other possible observables, such as the doublon density, are likely to behave similarly. The following increase of the entanglement entropy is due to a slower build up of long-range correlations. Finally, we discuss time t * of the first maxima in λ and quench-induced heating ∆E/L with respect to the quenched interaction (see Fig. 4). We have checked that the latter is converged in the system size despite the open boundary condition. Here, ∆E = ψ(t)|H f |ψ(t) − E f 0 is the energy of the time-evolving state relative to the ground state of the quenched Hamiltonian H f . Since the system is isolated, the energy is constant for t > 0.
The induced energy densities (green) behave approximately linearly for V f V eq c2 for initial states in the MI [ Fig. 4 Fig. 4(b)] phases. The values and the slope are larger for the initial state in the MI phase. This larger slope is a consequence of a higher doublon density in the HI phase, which reduces the extra energy due to a given increase of V f . Next, we look at the times t * of the first maximum of λ, which is either smooth (empty circles) or nonanalytic (full circles). For both kinds of initial states, t * attains a maximum at the threshold value V dyn c , where the nonanalyticities start to appear. As V f increases beyond V dyn c , the time t * decays in a power law manner, as indicated by the fits in Fig. 4 (dashed red lines). Heuristically, such power laws can be understood in terms of the width of the spectra of |ψ i with respect to the eigenstates of H f (see SM [57]), which broaden as the quenches become stronger. Interestingly, starting from an HI state leads to an apparent divergence of the extrapolation of the power law at V i , t * ∼ (V f − V i ) −1.2 . Such an approximate power law with a divergence at the initial parameter is typical of DQPTs in noninteracting topological systems [10,11,58,59]. Thus, the present results further underpin the topological nature of quenches starting from the HI phase. On the other hand, for the interactiondriven DQPTs, which occur with initial states in the MI phase, t * is not well approximated by a simple power law in (V f − V i ).

(a)] and HI [
Conclusions. In this work we have analyzed the dynamics of the extended Bose-Hubbard model after sudden interaction quenches. We have contrasted initial states in the Mott and Haldane insulator phases. In case of the Mott insulator, DQPTs are induced by sufficiently strong quenches. These DQPTs correspondence to the equilibrium phase boundaries and bear no relation to topology. Near the threshold quench value, nonanalytic signatures of the parity order parameter accompany the DQPTs. Starting from the Haldane insulator, we find DQPTs and zeros of the string order parameter for quenches to larger nearest-neighbor interactions. While, for quenches to lower interactions, the short time behavior differs from a previous study [37], we expect that DQPTs occur for longer times. The discrepancy is attributed to different symmetries of initial states. Finally, we have shown that the timescales of DQPTs depend on the quench parameter in a power law manner, and its precise form differs between two types of DQPTs-either topological or interaction-driven. Experimental tests of these results by ultracold atoms could be feasible in the near future. On the theoretical side, further study of the timescales of the dynamics is an intriguing open problem.
Supplemental material for "interaction-driven dynamical quantum phase transitions in a one-dimensional boson system" We show that the renormalized order parameters, which are introduced in Eqs. (4)-(8) of the main text, provide a well defined way of determining the phase diagram. In Fig. 1, both the x-and z-components of the order parameters, O z string,parity andÕ x string,parity , are depicted as a function of V for increasing system sizes. Even though there is no SU(2) symmetry in the extended Bose-Hubbard model (EBHM) due to the lack of particle-hole symmetry, the order parameters of both components clearly distinguish the different phases. The Mott insulator (MI) is characterized by non-vanishing parity order parameters, while in the Haldane insulator (HI) phase, the string order parameters take finite values. Contrary to the renormalized x-components, the bare x-components of the order parameters are insufficient to distinguish different phases. In Fig. 2, we show the correlation functions Ô x string,parity (r) and the projection string operator P (r) over the distance r = |i − j| between two sites for selected points in the MI (full lines) and HI (dashed lines) phase diagram. All correlation functions decay exponentially, as in Eq. (7) of the main text. For the projection stringP (r) = i<k<i+r P n k ≤2 the decay length is given by i.e., for each site the projection to n ≤ 2 contributes the same factor. In the MI phase, the decay length of Ô x parity (r) is equal to ξ, whereas in the HI phase it is smaller than ξ. Therefore, in the MI phase, the decay of Ô x parity (r) is exclusively due to the projections in Eq. (5) of the main text, and we consider the system to have x-parity order in the renormalized sense. For Ô x string (r) , the decay length is given by the projection in the HI phase, while it decays faster due to its own fluctuations in the MI phase.
Naturally, it is desirable to have a definition of x-parity and x-string order that does not rely on projections.

arXiv:2106.10187v1 [cond-mat.quant-gas] 18 Jun 2021
However, this is not trivial. Note that the x operators are mainly products of acting on each site k, by analogy to e iπŜ x k [see Eq. (5) of the main text]. A straightforward generalization to n > 2, would require terms like |2 − n n| with 2 − n < 0. Another simple ansatz, adding the terms |n n| for n > 2, also fails. Here the reason is nonlocality. Because the filling factor is n = 1, every site with n > 2 requires that there are more sites with n = 0 than with n = 2. The application of k (Ô k + n k >2 |n k n k |) would then change the total particle number and therefore yield an expectation value of zero, if any site has n > 2.
In conclusion, on the one hand, such operators may exist that correspond to the x-string and x-parity operators of spin-1 chains without prescriptions of projections. On the other hand, it is unclear whether these operators can be expressed simply as products of local operators.

II. LONG-TIME DYNAMICS FOR THE HI-CDW QUENCH WITH
Here, for the case of an initial state in the MI phase, we provide a finer parameter grid around the threshold value V f V dyn c for longer times up to t/J = 4 with L = 120 (see Fig. 3). We clearly see the relation between the time evolution of the rate function [ Fig. 3(a)] and that of the order parameters [ Fig. 3

III. LONG-TIME DYNAMICS FOR THE HI-MI QUENCH
In order to show that DQPTs likely exist for HI to MI quenches, we show the results of a simulation up to t/J = 5 with V i = 3.25 and V f = 1.0. While the finite size effects are too strong to reveal a kink of λ, there is clearly a zero ofÕ x string close to t/J = 4. This indicates that DQPTs occur also for the HI-MI quench in long-time simulations. Furthermore, as in the inset of Fig. 2(c) of the main text, we see that the order parameter vanishes in the thermodynamic limit after a certain time.

IV. PHASE DIAGRAM AND DYNAMICS OF SPIN-1 XXZ CHAINS
A spin-1 XXZ chain with Hamiltonian can serve as an effective model for the extended Bose-Hubbard model, when the on-site interaction U is large. As a benchmark for how well the effective description works, we show a phase diagram and the corresponding order parameters in Fig. 5. String and parity order parameters for the spin-1 chain are given by where α = x, z. There is a rough correspondence D → 1/4 U and ∆ → 1/2 V , when J = 1 in both models, and if one replaces a † i → S − i . Thus, from U = 5 in the bosonic model, we obtain D = 1.25. Along this line, one moves from the large-D phase (corresponding to MI) for small ∆ through the Haldane phase into the Neel phase (corresponding to CDW) for large ∆. The phase transition points are at ∆ eq c1 = 1.3 ± 0.1 and ∆ eq c1 = 1.9 ± 0.1. Qualitatively, we have a very good agreement between the models, while numerical values obviously differ.
Finally let us briefly touch on the signs before J in the two models (Here, we always take J to be positive). We take the sign before J in the EBHM negative, while in Eq. (3) it is positive as in Ref. [2]. However, the difference is not essential. We have checked the case of an XXZ chain with the negative sign and found that-both for static and dynamical properties-the only difference is a staggered sign (−1) |i−j| that has to be included into the string and parity correlation functions. Therefore, the results are basically unaffected by the sign. A related discussion can be found in Ref. [3], distinguishing the Haldane phases with opposite signs as being protected by different symmetries. While for the positive sign the symmetry is a regular lattice inversion symmetry, for the negative sign this symmetry needs to be slightly modified. Let us also compare the dynamics of the EBHM to the spin-1 model. We fix J = 1 and D = 1.25. As an analog to the initial MI ground state we choose ∆ i = 0.5 in the large-D phase (see Fig. 6), and for the initial HI ground state we choose ∆ i = 1.6 (see Fig. 7). These initial values are based on Fig. 5. After starting in the large-D phase, we find the same kind of DQPTs as in the EBHM when the quench is strong. The threshold for DQPTs, 1.4 < ∆ th < 2.0, does not coincide with the phase transition points ∆ eq c1 and ∆ eq c2 [see Fig.6(a)]. More importantly, as shown in Fig. 6(c), a sharp minimum of O x parity appears at ∆ f = ∆ th at the same time as the first kink of the rate function λ. For ∆ f > ∆ th , e.g. ∆ f = 3.0, this minimum of O x parity occurs at later times than the first kink of λ. The finite size dependence of O x parity is depicted in the inset of Fig.6(c). In contrast to the projected x-parity operator of the EBHM in the main text, it converges in the system size.
Lastly, the dynamics starting from a state in the Haldane phase is shown in Fig. 7. We find qualitative agreement with the EBHM in the main text. In particular, as can be seen from Fig. 7(b) for quenches towards smaller ∆'s, there are no DQPTs on the time scale of the simulation, while the rate function does show several smooth maxima. This confirms that the present initial state, which has a lower symmetry than the ground state of the Heisenberg model (J = 1, D = 0, V = 1), has a significant influence on the timescales of DQPTs.

V. SPECTRUM
In order to show how the energy is brought into the system by a quench, we have a closer look at the energy spectrum of the initial state |ψ i with respect to the post-quench Hamiltonian H f , i.e., where H f |E f j = E f j |E f j . In the continuum limit, the overlap of the time evolved and the initial wavefunction can be written as The distribution p(ε) is given by As this is difficult to obtain in the continuum limit, we have calculated the overlap amplitudes c j by full diagonalization for L = 10, and approximated δ(ε − E f j ) by a Lorentzian of width η = 0.1. Figs. 8 and 9 show these spectra for V i = 1.0 and V i = 3.25, respectively, for various V f with the open boundary condition. As V f increases, the means gradually shift to higher values and the distributions become wider. For V f = 6, we see a clear signature of a gap formation in the spectrum. We checked that the same happens for the periodic boundary condition, where, due to the momentum conservation, the number of eigenstates with a finite overlap with |ψ i is roughly decreased by 1/L, and the peaks are more sparsely distributed.
A heuristic model that neglects the intricacies of the actual spectra is given by a box with mean E and width ∆E, The corresponding Fourier transform is the Loschmidt-amplitude, One can see that the times of its zeros scale as t * ∼ (∆E) −1 . Thus we have a most rudimentary example of how the width of p(ε) controls the times of DQPTs in a power law fashion. As the width of the more complicated spectrum of the EBHM increases for stronger quenches (see Figs. 8 and 9), and the critical times t * of the EBHM follow approximate power laws as a function of the quench parameter (see Fig. 4 of the main text), the spectrum of the EBHM likely controls the times of DQPTs in a similar way as p E,∆E (ε) in the above toy model.