Fractonic Quantum Quench in Dipole-constrained Bosons

We investigate the quench dynamics in the dipolar Bose-Hubbard model (DBHM) in one dimension. The boson hopping is constrained by dipole conservation and show fractonic dynamics. The ground states at large Hubbard interaction $U$ are Mott insulators at integer filling and a period-2 charge density wave (CDW) at half-integer filling. We focus on Mott-to-Mott and CDW-to-CDW quenches and find that dipole correlation spreading shows the light-cone behavior with the Lieb-Robinson (LR) velocity proportional to the dipole kinetic energy $J$ and the square of the density in the case of Mott quench at integer filling. Effective model for post-quench dynamics is constructed under the dilute-dipole approximation and fits the numerical results well. For CDW quench we observe a much reduced LR velocity of order $J^2/U$ and additional periodic features in the time direction. The emergence of CDW ground state and the reduced LR velocity at half-integer filling can both be understood by careful application of the second-order perturbation theory. The oscillatory behavior arises from quantum scars in the quadrupole sector of the spectrum and is captured by a PXP-like model that we derive by projecting the DBHM to the quadrupolar sector of the Hilbert space.

Introduction.-Dipole-conservingsystems are a simple example of the particle dynamics and the many-body phases being altered in a fundamental way by kinetic constraints .In addition to the immobility of single particles reminiscent of the fractonic dynamics, other novel phenomena such as the lack of thermalization, Hilbert space fragmentation, and quantum scars are all manifested in the dipole-constrained systems [3,4,11,12,18,21,25,26].More recently they have received a great deal of attention as ways to understand anomalous transport and relaxation phenomena in tilted optical lattices [19][20][21][22]26].
Over the years the optical lattice system has proven to be excellent platforms for probing non-equilibrium states of matter.A prototypical example of non-equilibrium probe is the quench dynamics where a sudden change of system parameters results in the ground state evolving according to the post-quench Hamiltonian.Some intriguing aspects of the post-quench dynamics have been examined in the past, ranging from light cone-like information spreading subject to the Lieb-Robinson bounds [27][28][29][30][31][32], dynamical quantum phase transition (DQPT) [33,34], and quantum scars [25,[35][36][37].The issues have been addressed in the framework of e.g.Bose-Hubbard model [27,28,30], transverse Ising model [34], and PXP model [25,36].
Motivated by recent experiments in tilted optical lattices, several interacting models embodying the dipole conservation in addition to the charge conservation have been proposed [15-17, 20, 24].An interesting ramification of one such model, called the dipolar Bose-Hubbard model (DBHM) [15][16][17], is the disappearance of conventional superfluid phase and the emergence of dipole condensate phase taking its place in the weak Hubbard interaction regime.The ground state phase diagrams and various low-energy correlations of this model have been worked out.Notably, single-particle correlations * hyunyong@korea.ac.kr are heavily suppressed in all phases of the model and twoparticle dipole-dipole correlations take over as a measure of (quasi-)ordering.Recent progress in experiments shows that DBHM and its fermionic cousin, the dipolar Fermi-Hubbard model, are among the most experimentally accessible models displaying fractonic quasiparticle behavior through enforcing the dipole symmetry [19][20][21][22]26].
Despite the growing importance of dipole-constrained models with roots in tilted optical lattice, the quench dynamics of DBHM has not been examined theoretically.Here we present the first thorough study of the quench dynamics over different phases of DBHM at integer and half-integer fillings.Due to the strict prohibition of single-particle dynamics, dipoles as low-energy excitations become the main channel of correlation spreading.The Lieb-Robinson (LR) bound for the Bose-Hubbard model, which scales linearly with the density, is replaced by a new bound scaling as the square of the density in DBHM.At half-integer filling where the ground state is a period-2 charge-density wave (CDW), dipole correlation spreading is bounded by a much smaller LR speed and a periodic (in time) revival, reminiscent of quantum scars.Effective models for the post-quench dynamics in both integer and half-filling fillings are derived in terms of a low density of dipole excitations and can explain the numerically observed LR bound quantitatively.Furthermore, a PXP-like model consistent with the scar-like features in the half-integer quench can be derived by taking into account quadrupole excitations, and explain the observed periodicity very well.
Model and methods.-Theone-dimensional DBHM is [15][16][17] where n x = b † x b x is the boson number at site x, and n = ∑ x n x /L (L=number of sites) is the average density.The key departure from BHM is the absence of one-boson hopping and the dipo-lar hopping (J) that takes its place.The model is invariant under both the global U(1) and the dipolar U(1) phase changes b x → e iθ b x , b x → e iθ x b x , and possesses two conserved quantities: the total charge Q = ∑ x b † x b x and the dipole moment x b x .We employ the density matrix renormalization group (DMRG) [38][39][40][41] and the time-dependent variational principle (TDVP) [42,43] calculations to explore the ground state and its quench dynamics.For DMRG simulations, we utilize the two-site and subspace expansion algorithms [44], focusing on a finite system with size L = 100 and limiting the local boson number at each site to 10.The maximum bond dimension for DMRG is set to χ DMRG = 500 ensuring an accurate representation of the ground state in the matrix product states representation.In the context of TDVP, we adopt both one-site and two-site algorithms, with the maximum bond dimension up to χ TDVP = 3000.This substantial increase in the maximum bond dimension allows for a more detailed exploration of the system's dynamics.We also incorporate the conservation of boson number Q and dipole moment D in both DMRG and TDVP simulations.It not only guarantees the conservation of associated U(1) symmetries but also greatly enhances the computational efficiency of the simulations [45].
Mott quench.-We obtain the ground state |ψ⟩ of the DBHM [Eq.( 1)] at U = U i and observe their evolution under the new Hamiltonian with U = U f as |ψ(t)⟩ = e −itH DBHM |ψ⟩.The final value of U f is chosen such that the equilibrium state corresponding to U = U f is also in a Mott phase.We refer to such quench as the Mott-to-Mott quench, or simply Mott quench.In the simulation J is set to unity.With |ψ(t)⟩ we examine the time evolution of relevant quantities such as correlators and fidelities.The single-boson correlation ⟨ψ(t)|b † x b x ′ |ψ(t)⟩ remains strictly zero except x = x ′ at all times due to the dipole constraint, indicating the fractonic nature of the single-boson particle in the DBHM.Instead, meaningful information is contained in the dipole correlator is the dipole operator.In the TDVP simulation we choose x 0 = L/2 at the center of the system 1 ≤ x ≤ L; results are unaffected by the choice of x 0 unless it is positioned too close to the boundary.
We begin by focusing on the integer-filling n x = n at large U/J where the ground state is a Mott state, faithfully represented as a product state |M⟩ ≡ ⊗ L x=1 |n⟩ x .The dipole correlation function in the Mott state is extremely short-ranged, but the quench triggers the spreading of the correlation with a well-defined propagation front in the shape of a light-cone as shown in Fig. 1 (a).A modern interpretation of this is in terms of the LR bound [46], whose existence has been rigorously proven for the conventional BHM [31,32,47,48] after many years of numerical observation to the effect [28][29][30][49][50][51].The light-cone spreading of dipole correlation in Fig. 1 (a) is highly suggestive of an LR bound in the DBHM as well, with the information carried in the dipole, not charge, sector.To make a quantitative statement on the LR bound of the DBHM, we extract the group velocity (v g ) and also the phase veloc- Theory TDVP ity (v p ) from the TDVP data by fitting the leading wave packets in the dipole correlation [52].See Supplementary Material (SM) for details on how to determine the velocities [53].The results are presented in Fig. 1 (b) as a function of U f at filling n = 1 and 2, strongly indicating that v g remains independent of U f , whereas v p exhibits a linear dependence on it.In terms of the filling factor dependence, v g appears to increase with the square of the filling factor, while v p remains independent of it.
The Mott quench dynamics can be comprehensively understood by developing an effective model deep inside the Mott phase U ≫ J.The low-lying excitations in the Mott phase are the two kinds of dipole excitations |l x ⟩ ∼ d x |M⟩ and |r x ⟩ ∼ d † x |M⟩ called l-dipoles and r-dipoles, respectively.Considering a Hilbert subspace consisting of the Mott state |M⟩ and the dipoles {|l x ⟩, |r y ⟩}, the effective Hamiltonian in this space can be derived [53] The parameters ρ k and λ k are given by [53] for general integer filling factor n. The operator γ † kσ creates Bogoliubov quasi-particles with pseudo-spin σ = l, r and momentum k.There is much resemblance of this effective model to the quasiparticle model for Mott quench in the BHM [29,30], with the key difference that dipoles rather than doublons and holons are the elementary excitations here.
The post-quench wavefunction |ψ(t)⟩ can be derived in exact form using the Bogoliubov Hamiltonian and allows a closed-form expression of the dipole correlation function at t > 0 [53], which shows excellent agreement with the TDVP simulations as shown in Fig. 1 (a).One can deduce the two propagation velocities analytically from the Bogoliubov model as follows [53]: When U ≫ Jn 2 , v g is maximized at k max = π/2.The two velocity expressions provide very good fits to the velocities extracted from the TDVP data, as shown in Fig. 1 (b).Furthermore, the quadratic dependence of the group velocity on the density v g ∼ n(n + 1) deduced from effective model captures the observed increase in v g by three times in going from n = 1 to n = 2.This contrasts with its linear dependence on n in the conventional BHM [30,32].The group velocity in the DBHM approaches zero as n → 0, whereas it remains finite in the conventional BHM [32,48].
CDW quench.-The correlation spreading at half-integer filling n+1/2 shows a number of features which distinguish it sharply from those in the Mott quench.Though the discussion is based on detailed numerics at n = 3/2, the results straightforwardly generalize to arbitrary half-integer filling.Firstly, the ground state at half-filling obtained by DMRG is a period-2 CDW state with an alternate occupation of one boson and two bosons per site.The LR velocity bounding the correlation spreading in the CDW quench scales as J 2 /U f and substantially smaller than the Mott-quench value which scales as J. Finally, the dipole correlation functions show a periodic revival in time that was absent in the Mott quench.Both these features are apparent in the plots shown in Fig. 2.
First we discuss the origin of the CDW ground state at halffilling.The Hubbard term at half-filling demonstrates extensive degeneracy, with any state with half the sites occupied by one boson and the other half with two bosons sharing the same Hubbard energy.The massive degeneracy is lifted at the second-order of J/U in degenerate perturbation theory, resulting in two-fold degenerate CDW ground states [53].Without loss of generality, we choose the CDW state on an open chain of length L to be |CDW⟩ ≡ ⊗ As in the Mott quench, low-lying excitations are those of land r-dipoles, created in equal numbers to preserve the total dipole moment.Due to the translation symmetry breaking of the CDW, however, the l-dipoles (r-dipoles) are created at odd (even) sites only, given by the change in the local occupation Degenerate perturbation theory leads to an effective Hamiltonian of the dipoles in the CDW state [53]: The superscript D is a reminder that only the dipole states comprise the low-energy Hilbert space, of order J 2 /U above the CDW ground state per dipole, used to construct the effective Hamiltonian.
In constructing the effective model we ruled out configurations where the (l, r) dipoles are adjacent, i.e.
They are in fact quadrupole and anti-quadrupole excitations, and cost an energy of order U more than two separately created dipoles.Ignoring the quadrupole events, the effective Hamiltonian (4) can be diagonalized with the dispersion ω k = (J 2 /U)(32.8− 24 cos 2k) ≥ 8.8J 2 /U.The factor 2 in cos 2k appears as a result of the unit cell doubling.The group velocity is deduced v g = max(2∂ k ω k ) = 96J 2 /U.This prediction, shown as black solid lines in Fig. 2, agrees very well with the TDVP results for the propagation boundary of the dipole correlation function.Being of order J 2 /U, the LR velocity is considerably smaller than the v LR ∼ J in the Mott quench and, moreover, depends inversely on U in marked contrast to the Mott quench at integer filling or the quench in the conventional Bose-Hubbard model where it is governed exclusively by kinetic energy J.
The other prominent features in the CDW quench, i.e. a periodic revival of the correlation, originates from quadrupole  [36,37,54].Motivated by the fidelity results, we consider a new subspace in which only the quadrupole states are kept along with the CDW ground state.The connection between dipole and quadrupole sectors are considered negligible, given that such connections are made by an intermediate state of energy 2U or more, and are suppressed by a factor of J/U in perturbation theory.On the other hand, quadrupole states are linked to the CDW by a single application of the dipolar kinetic term in DBHM as: where we introduce the quadrupole operator q x = d † x d x+1 .The algebra suggests that at every odd site x = 2a − 1, {|CDW⟩, |131⟩} forms a two-level system and at every even site x = 2a, {|CDW⟩, |202⟩} forms another two-level system.This structure can be effectively modeled by assigning pseudo-spin-1/2 operators (X, Z) acting on an effective qubit to every site in the lattice: The CDW state maps to |0⟩ = ⊗ L x=1 |0⟩ x .Projecting the DBHM to the Hilbert space of quadrupoles gives with CDW as the vacuum.The projector P x = (1 + Z x )/2 projects a local state to |0⟩ x .The P x−1 X x P x+1 in the first line means that if one tries to create a 131-quadrupole at the odd site x, one can only do so if both of the adjacent sites (x − 1, x + 1) are devoid of existing 202-quadrupoles.
Otherwise, creating a 131-quadrupole on top of an existing 202-quadrupole at the adjacent site annihilates the state altogether.In the second line, the X-operator tries to create a 202-quadrupole at the even site x, provided that its two adjacent sites are devoid of any existing 131-quadrupoles (hence P x−1 X x P x+1 ).Extra projectors P x−2 , P x+2 are used because creating two 202-quadrupoles at adjacent positions like x and x +2 leads to an occupation of 20402 and an additional energy cost of 4U compared to separate 202-quadrupoles.In the first line, the projectors at second-neighbor sites are omitted since generating two 131-quadrupoles at sites x and x + 2 yields an occupation of 13031 without incurring additional energy compared to separate 131-quadrupoles.We refer to the emergent projected Hamiltonian as the PXPQ model, with "Q" referring to the quadrupole excitations.The on-site energy cost U acts as an effective magnetic field polarizing the state toward the CDW [26,55].Translational symmetry is explicitly broken by two lattice spacings in the PXPQ model.We calculate the post-quench wave function |ψ(t)⟩ = e −it H PXPQ |0⟩ and its overlap with |1, 0⟩ = (2/L) 1/2 ∑ x∈odd |1⟩ x and |0, 1⟩ = (2/L) 1/2 ∑ x∈even |1⟩ x using the PXPQ model.For the same values of (J,U), we find very good agreement in the fidelity evolution F(t) as shown in Fig. 2 (b) with the correct period as found in the DBHM.Having identified a PXP-like effective model governing the dynamics of quadrupoles, we address the important question of the nature of the quantum scar state in the PXPQ model, by constructing a coherent state of quadrupoles or a quadrupole-condensate (QC) where Discussion -The quench dynamics in the dipolar Bose-Hubbard model reveals that correlation spreading is mediated by dipole excitations.Effective models for the dipole dynamics can explain the observed quench dynamics at both integer and half-integer fillings, though in detail they are sub-stantially different in that the Lieb-Robinson velocity is set by J for the integer quench and by J 2 /U in the half-integer case.The CDW ground state at half-integer filling bring dramatic changes in the low-energy dipole dynamics.Furthermore, quantum scar states exist in the form of quadrupole excitations in the CDW quench, manifesting themselves as periodic oscillations in the dipole correlation function.The scarlike features in the CDW quench are captured by an effective model resembling the PXP Hamiltonian.The existence of CDW ground state as well as the novel quench dynamics at half-integer filling can be probed in future tilted optical lattice setup.
Supplementary Material for "Fractonic Quantum Quench in Dipole-constrained Bosons "

I. MOTT PHASE EFFECTIVE THEORY
The Mott-to-Mott quench processes can be understood quite well using the effective model constructed deep within the Mott regime U ≫ J, where a dilute gas of left and right dipoles dominates the low-energy spectrum.The two kinds of one-dipole states introduced in the main text are given by: where with unmarked sites occupied by n bosons.The dipole pair then drifts apart by further action of dipole hopping.
In terms of the dipole operators, the dipolar hopping operators can be replaced by The Hubbard interaction in the dipole subspace becomes This assumes that the dipoles are far apart, and each dipole costs energy +U.The dipole creation/annihilation processes take place only when they are adjacent, as indicated by the pair-creation and annihilation terms in Eq. (S1).This, however, is a rare event in the case of dilute-diplon regime, and for the most part the Hubbard energy is simply given by Eq. (S2).In the same dilute-dipole regime, r and l operators can be treated as ordinary boson operators subject to the hard-core constraints (r † x ) 2 = (l † x ) 2 = 0.The constraints are, in turn, resolved by mapping the boson model to the fermion model through Jordan-Wigner transformation [29,30].We follow the same footsteps and arrive at the effective Hamiltonian.
In the momentum space, the effective Hamiltonian becomes where After the Bogoliubov transformation, where one obtains with ω k = (ρ 2 k + λ 2 k ) 1/2 describing quasiparticle dynamics deep in the Mott phase of DBHM.The ground state of the post-quench Hamiltonian is given by γ l,k |M ′ ⟩ = 0 and γ r,k |M ′ ⟩ = 0 in the quasiparticle picture, related to the pre-quench ground state (Mott state) |M⟩ by One can show that l k |M⟩ = r k |M⟩ = 0.The time evolution of the post-quench state follows as
Performing the Fourier transformation and the Bogoliubov transformation in sequence, one can get Here, θ k and µ k are parameters for Bogoliubov transformation given in Eq. (S4), and ω k is the spectrum of post-quenched Hamiltonian given in Eq. (S5).Summarizing all together, the dipole correlator can be expressed as where the parameters t 0 , ω, k, and σ represent the center, frequency, wavenumber, and width of the wavepacket, respectively, which can be determined through fitting at given x.We consider x values ranging from 15 to 20, which are sufficiently far from the center yet not close to the boundary.For instance, Figure S1 (b) and (c) display the leading wavepacket along with the fitting function at x = 15 and 16, respectively.From t 0 values obtained for various x, one can estimate the group velocity using the formula: v g = (x 2 − x 1 )/(t 0 (x 2 ) − t 0 (x 1 )).For the phase velocity, we utilize the formula v p = ω/k, where both ω and k are obtained from the fitting procedure.

IV. FILLING NUMBER DEPENDENCY OF GROUP VELOCITY AT LARGE U LIMIT
The group velocity from the effective theory is given in Eq. (3) of the main text.More specific form is given by In the limit of U ≫ Jn 2 , one can deduce k max ≈ π/2.At k = k max , the group velocity is given by V. GROUP AND PHASE VELOCITIES AT FIXED U Here, we depict the group velocity (v g ) and the phase velocity (v p ) as a function of J while keeping U fixed.In other words, we consider a quench process: (J 0 ,U) → (J,U).Figure S2 illustrates the velocities calculated from Eq. (3) of the main text for both n = 1 and n = 2 as a function of J/J 0 , with U/J 0 = 40 kept constant.Here, J 0 serves as the normalized energy scale.In both cases, one can check the group velocity increases linearly with J, whereas the phase velocity remains relatively unaffected by the variations in J. Here, we provide supplementary results from the TDVP simulations of quench dynamics within the Mott phase of the DBHM for two specific filling factors: n = 1 and n = 2. Additionally, we compare these results with the ones obtained from the effective theory.See Figs.S3 and S4.FIG.S3.Spreading of the dipole correlations at the filling n = 1.Here, the initial state is the ground state at U i /J = 100.The data is normalized such that the maximum value is adjusted to unity.
Additionally, we include TDVP data for the quench dynamics transitioning from a smaller initial on-site interaction strength U i to a larger quench strength U f within the Mott phase, as illustrated in Fig. S5.Our findings indicate a speed of the dipole correlation is consistent that from the reverse quench direction discussed in the main text, where v g is approximately 20.This observation corroborates the effective theory's prediction that the group velocity depends only on the hopping strength J.

VII. CDW PHASE PERTURBATION THEORY
The ground state space of H 0 at half-integer filling (ν = 3/2) is denoted W . Having identified the CDW state as reference vacuum state in W , other low-lying excited states in W can be identified with the creation of l and r dipoles above the CDW.A subset of states in W is derived by transferring two bosons from two next-nearest neighboring sites within a CDW state that FIG.S4.Spreading of the dipole correlations at the filling n = 2. Here, the initial state is the ground state at U i /J = 100.The data is normalized such that the maximum value is adjusted to unity.
In each scenario, a pair of 1, 2 (magenta) and 2, 1 (cyan) swaps to form 2, 1 and 1, 2 pairs, respectively.The transitioned pair to 2, 1 is designated as an l dipole, and the 1, 2 pair as an r dipole, each seamlessly embedded into the otherwise perfect CDW pattern.A notable distinction from the Mott phase with integer filling is that the l dipole is exclusively generated and located on odd-numbered sites, whereas the r dipole is restricted to even-numbered sites.Considering the dipole-moment conservation, the numbers of l and r dipoles equal in the basis of W .The ground state space of H 0 at half-integer filling (ν = 3/2) is denoted W , and a particular state in W is |ψ W α ⟩.For any two states within this subspace we have ⟨ψ W α |H 1 |ψ W β ⟩ = 0, since H 1 changes the occupation at a particulate site by two (e.g. 2 → 0 or 1 → 3) and lifts the state out of W .We thus need to employ H 1 to second order to lift the degeneracy.
The effective Hamiltonian for subspace W , as derived using second-order perturbation theory, is given by

FIG. 1 .
FIG. 1.(a) TDVP and effective model results of the post-quench dipole correlation function (real part) C d (x,t) at the filling n = 1 and 2 with the initial value U i /J = 100.The data is normalized such that the maximum value is adjusted to unity.The dashed line and the solid arrows represent the travel speed of the overall wave packet and the peak in the response, respectively.(b) Phase and group velocities for the dipole correlation at n = 1 and 2 as a function of the quench interaction U f /J deduced from the TDVP data such as shown in (a).

1 FIG. 2 .
FIG. 2. TDVP results of the dipole correlation C d (x,t) at the filling n = 3/2, with the initial U i /J = 100.The Data is normalized such that the maximum value is adjusted to unity.The black solid line indicates the wavefront expected from the theory (See text).
x+1 ⟩ in the occupation number basis.Undesignated sites have the occupation n x = n.The Mott state |M⟩ = | • • • n x • • • ⟩ serves as the vacuum.The dipole-hopping J-term in DBHM acting on |M⟩ creates a pair of l and r dipoles:

)
Figure S1,(a) illustrates the propagation of the dipole correlation, denoted as C d (x,t) = Re[⟨ψ(t)|d † L/2 d L/2+x |ψ(t)⟩], as a function of time at distances x = 5, 7, and 9.We fit the wavepacket appearing earliest to a Gaussian wavepacket, FIG. S1.(a) Dynamics of dipole correlation as a function of time at x = 5, 7 and 9.The leading wavepacket with its fitting function at x = 15 and 16 with fitting parameters t 0 , ω and k.
FIG.S2.Phase and group velocities for the dipole correlation at n = 1 and 2 as a function of J/J 0 with constant U/J 0 = 40.The green and red lines represent the calculated group (v g ) and phase (v p ) velocities derived from the effective model using Eq.(3) of the main text.

2 FIG
FIG.S5.Propagation of the correlations at the filling n = 2. the quench direction is reverse, i.e., from a smaller U i /J = 20 to larger (left) U f /J = 30 and (right) U f /J = 40.The data is normalized such that the maximum value is adjusted to unity.
is the normalization, and α (β ) is the fugacity parameter of the 131 (202) excitation.Here, a projector P rules out configurations that contain overlapping two 202 excitations or neighboring 131-202 excitations.The |QC⟩ state can be represented by a matrix product state with a small bond dimension χ = 3 and 4 in PXPQ and DBHM, respectively[53].Note that the amplitude |⟨QC|ψ(t)⟩| 2 grows substantially and reaches a maximum as the overlap |⟨CDW|ψ(t)⟩| 2 or |⟨0|ψ(t)⟩| 2 is most suppressed at periodic intervals [Fig.3(c) and (d)],showing there is a periodic transfer of weight from CDW to QC and back.The dependence of amplitude on fugacities strongly suggests that the quantum scar state approximates the QC state with small (α, β ) ≪ 1, corresponding to a dilute quadrupole density.