Evolution of entanglement entropy in strongly correlated bosons in an optical lattice

We investigate the time evolution of the second-order R\'enyi entropy (RE) for bosons in a one-dimensional optical lattice following a sudden quench of the hopping amplitude $J$. Specifically, we examine systems that are quenched into the strongly correlated Mott-insulating (MI) regime with $J/U\ll 1$ ($U$ denotes the strength of the on-site repulsive interaction) from the MI limit with $J=0$. In this regime, the low-energy excited states can be effectively described by fermionic quasiparticles known as doublons and holons. They are excited in entangled pairs through the quench dynamics. By developing an effective theory, we derive a direct relation between the RE and correlation functions associated with doublons and holons. This relation allows us to analytically calculate the RE and obtain a physical picture for the RE, both in the ground state and during time evolution through the quench dynamics, in terms of doublon holon pairs. In particular, we show that the RE is proportional to the population of doublon-holon pairs that span the boundary of the subsystem. Our quasiparticle picture introduces some remarkable features that are absent in previous studies on the dynamics of entanglement entropy in free-fermion models. It provides with valuable insights into the dynamics of entanglement entropy in strongly-correlated systems.


I. INTRODUCTION
Entanglement is one of the most intriguing concepts of quantum mechanics.It describes non-local correlations incompatible with local realism [1], which is clearly demonstrated by the violation of the Bell inequality [2].Entanglement is also a key to understanding quantum many-body systems in diverse fields.It is considered, for example, to be the origin of thermalization in an isolated quantum many-body system [3][4][5][6][7][8] and the Hawking radiation from black holes [9,10].In particular, entanglement plays a central role in strongly correlated systems.The investigation of entanglement in strongly correlated systems is promised to give us deep insights into fundamental aspects of emergent phenomena, such as quantum phase transition and topological order [11][12][13][14][15].
Entanglement between quantum objects can be quantified by entanglement entropy.It has been a major subject of theoretical investigation in quantum field theory, as well as in strongly correlated systems.Dynamics of entanglement entropy in integrable systems have been intensively investigated since the pioneering work by Calabrese and Cardy [16].They proposed a clear physical picture for the dynamics of entanglement entropy in terms of quasiparticles.Specifically, the long-time dynamics of entanglement entropy can be understood as a result of excitation and propagation of entangled quasiparticle pairs.This quasiparticle picture has been confirmed numerically and analytically in a number of papers [17][18][19][20][21][22][23][24][25][26][27][28].
Despite recent developments of experimental techniques, measuring entanglement entropy remains challenging in condensed matter systems.A great advance has been recently made, however, in the system of ultracold bosonic atoms in an optical lattice.The secondorder Rényi entropy (RE), which is one of the measures of entanglement entropy, has been successfully probed by preparing two independent copies of the same state, letting them interfere, and counting the number parity of atoms in one of the copies by an atomic gas microscope [29,30].The time evolution of the RE after a sudden quench of atomic hopping has been observed in the superfluid (SF) regime by using this technique [31].In the strongly correlated Mott insulating (MI) regime, on the other hand, spreading of correlations after a quantum quench has been studied theoretically and experimentally [32][33][34].However, quench dynamics of entanglement entropy has not been well investigated in this regime.
In this paper, motivated by the current status of theory and experiment, we study the quench dynamics of entanglement entropy of bosons in a one-dimensional optical lattice.Our main focus is on the quench dynamics of the RE in the strongly correlated Mott insulating regime.The low-energy dynamics in this regime can be effectively described by fermionic quasiparticles known as a doublon and a holon, which correspond to an excess particle and a hole on top of the unit filling, respectively [32].We develop an effective theory to derive an analytical expression for the time evolution of the RE after a quench of atomic hopping.Furthermore, we derive a direct relation between the RE and correlation functions for doublons and holons, which enables us to obtain a physical picture for the dynamics of the RE in terms of doublon-holon pairs.We find that the obtained quasiparticle picture is consistent with the one proposed by Calabrese and Cardy in the space-time scaling limit.Moreover, it exhibits re-arXiv:2209.13340v5[cond-mat.quant-gas]14 Nov 2023 markable features in both the short and long time scales that are absent in their picture.
The organization of the paper is as follows: In Sec.II, we explain the model and setup for the quench dynamics of bosons in an optical lattice and introduce the RE.In Sec.III, we introduce the effective theory in the strongly correlated MI regime.In Sec.IV, we introduce the formalism to calculate the RE and study the RE for the ground state.In Sec.V, we study the time evolution of the RE after a quench.In Sec.VI, we discuss the physical picture for the time evolution of the RE.In Sec.VII, we extend the analysis to study the n-th order RE.In Sec.VIII, we examine the validity of the effective theory.We finally summarize the paper in Sec.IX.We set ℏ = 1 and the lattice constant unity throughout this paper.

II. MODEL AND SETUP
We consider bosons in a one-dimensional (1D) optical lattice at zero temperature.When the lattice potential is deep enough, the system is well described by the Bose-Hubbard model (BHM) [35][36][37] where bj ( b † j ) denotes the annihilation (creation) operator of a boson on the jth site and nj = b † j bj the number operator on the jth site.J denotes the hopping amplitude between nearest-neighbor sites and U > 0 the strength of the on-site repulsive interaction.We assume the periodic boundary condition.
The BHM (1) exhibits a quantum phase transition between the SF and MI phases [38][39][40][41][42][43][44]: When the total number of bosons N is commensurate with the number of total sites L, the ground state is a SF state for small U/J, while it is a MI state for large U/J.The SF-MI phase transition of the Kosterlitz-Thouless type occurs at U/J ≃ 3.28 for unit filling (N/L = 1) [41][42][43][44].The ground state is a SF state when N is incommensurate with L regardless of the value of U/J.
We suppose that the whole system consists of subsystems A and B. The RE for subsystem A is defined as [45] where ρA = tr B (ρ) is the reduced density matrix for subsystem A and ρ is the density matrix for the whole system.tr A(B) stands for trace over subsystem A (B). tr A (ρ 2 A ) quantifies the purity of the state ρA [46]: tr A (ρ 2 A ) = 1 if ρA is a pure state, while tr A (ρ 2 A ) < 1 if it is a mixed state.When subsystems A and B have no entanglement, ρA describes a pure state and we obtain S A = 0.When subsystems A and B are entangled, ρA describes a mixed state and we obtain S A > 0.
We follow the quench protocol of the experiments [30][31][32].Namely, atoms are initially localized one in each of the lattice sites.At the initial time (t = 0), tunneling of atoms is abruptly switched on by lowering the lattice depth and the state of the whole system |ψ(t)⟩ evolves following the Hamiltonian (1) as |ψ(t)⟩ = e −i Ĥt |ψ 0 ⟩, where |ψ 0 ⟩ is the initial state at t = 0.
The initial state can be written as where |ν⟩ j (ν = 0, 1, 2, . . . ) denotes the Fock state with ν atoms on the jth site.It corresponds to the ground state of the MI limit (J/U = 0).Since the initial state (3) is a product state, S A = 0 at t = 0. S A grows in time after the quench as tunneling of bosons creates entanglement between the subsystems.

III. EFFECTIVE THEORY IN THE STRONGLY-CORRELATED MOTT INSULATING REGIME
We assume that the lattice potential is slightly lowered and the value of J/U is set in the strongly correlated MI regime (J/U ≪ 1) at t > 0. The low-energy excited states in this regime can be described in terms of doublons and holons [32,47].Such a weak perturbation associated with the quench involves only the low-energy excited states.As a result, the time evolution of the system after the quench is considered to be well described by the effective theory based on the doublon-holon description.
Introducing the fermionic doublon and holon annihilation (creation) operators, dj ( d † j ) and ĥj ( ĥ † j ), respectively, the Hamiltonian (1) is approximately mapped to [32,47] Ĥ = P Ĥeff P , where Ĥeff is the Hamiltonian of the effective theory given by and P = j (1 − d † j dj ĥ † j ĥj ) is the projection operator, which eliminates double occupancy of a doublon and a holon on the same site.The derivation of Eqs. ( 4) and ( 5) is given in Appendix A.
The projection operator P can be safely neglected in weakly excited states of the strongly correlated MI regime since the system can be considered as a dilute gas of doublons and holons and the possibility of their occupation on the same site is quite low.
By the Fourier transform, Ĥeff can be written as where , and g k = 2 √ 2J sin(k).Doublon and holon have energy gap f d,k=0 = U/2 − 4J and f h,k=0 = U/2 − 2J, respectively.Note that the initial state |ψ 0 ⟩ in Eq. ( 3) corresponds to the vacuum state of dj and ĥj .
The quadratic Hamiltonian Ĥeff can be diagonalized by the Bogoliubov transformation where γd,k and γh,k denote the annihilation operators of quasiparticles, which we refer to as "bogolons" hereafter.
u k and v k are given by where we expand u k and v k in terms of J/U for later use.Substituting Eq. ( 7) into Eq.( 6), Ĥeff is diagonalized as where the dispersions of bogolons are given by They have energy gap The ground state |vac⟩ that satisfies γd,k |vac⟩ = γh,k |vac⟩ = 0 can be written as [48] It implies that doublon-holon pairs are condensed in the ground state from its similarity with the BCS wave function [49].
The time-evolving state after the quench is given as where we used |vac⟩, which can be obtained from dk |ψ 0 ⟩ = ĥk |ψ 0 ⟩ = 0.It shows that pairs of bogolons are excited by the quench.Equation (14) will be used to calculate the time evolution of the RE.
In terms of dk and ĥk , Eq. ( 14) can be written as Equation (15) indicates that doublons and holons are excited in pairs.The numbers of doublons and holons should be thus equal: Using this relation, the total number of original bosons is conserved as where N is the total number of original bosons at t = 0. Note that the number of doublons (holons) itself is not conserved because they are excited from the vacuum state |ψ 0 ⟩.

IV. R ÉNYI ENTROPY FOR THE GROUND STATE
Let us first calculate the RE for the ground state Eq. ( 13) before studying its time evolution.For a Gaussian state, which includes the ground state and a thermal state of a quadratic Hamiltonian, the RE can be conveniently evaluated using single-particle correlation functions [50].We first adopt the formalism to our system in Sec.IV A, and apply it to the ground state in Sec.IV B.

A. Rényi entropy for a Gaussian state
We consider a Gaussian state of doublons and holons |ϕ⟩ with which any correlation function of dj and ĥj factorizes according to the prescriptions of Wick's theorem.The reduced density matrix of |ϕ⟩ can be formally written as [50] The entanglement Hamiltonian ĤA has a quadratic form of dj and ĥj (j ∈ A), because Wick's theorem also holds for correlation functions concerning degree of freedom in A. Thus, it can be diagonalized as where ω A α and nA α are the spectrum and the number operator for the eigenmode α.L A is the size of subsystem A. Note that the number of the eigenmodes 2L A corresponds to the total number of degrees of freedom in subsystem A. Once the entanglement Hamiltonian is diagonalized in the form of Eqs.(18), the RE can be obtained as where is the occupation number of the eigenmode α.
We still need to determine ω α or f α .To obtain f α , we consider a matrix of single-particle correlation functions Here, C σ,σ ′ and F σ,σ ′ (σ, σ ′ = d, h) are L A × L A matrices of normal and anomalous correlation functions, respectively, that have matrix elements where we denote ĉi,d = di and ĉi,h = ĥi .f α can be obtained by diagonalizing M thanks to the relation [50] where U A is a unitary matrix.

B. Rényi entropy for the ground state
We calculate the RE for the ground state Eq. ( 13) employing the formalism in Sec.IV A. Evaluating the singleparticle correlation functions in Eqs.(21) and (22) with |vac⟩, we obtain From u k = O(1) and v k = O(J/U ) in Eqs. ( 8) and (9), we obtain C i,j = O[(J/U ) 2 ] and F i,j = O(J/U ), because the summations over k in Eqs. ( 24) and ( 25) do not change the order of J/U .
Using matrix C and F , the matrix of the single-particle correlation functions M in Eq. ( 20) can be simplified as Here, F is antisymmetric, i.e., F i,j = −F j,i , due to the anticommutation relation { di , ĥj } = 0. Note that this property holds for any state with which the average is taken.
We derive a formula that directly relates the RE and the single-particle correlation functions.From Eq. ( 26), we obtain (27) where ∥O∥ F = i,j |O i,j | 2 denotes the Frobenius norm.In deriving Eq. ( 27), we use On the other hand, using Eq. ( 23), we obtain Comparing Eqs. ( 27) and ( 28), we obtain a relation between f α and the correlation functions as The first term of the right-hand side in Eq. ( 29) is of the order of (J/U ) (29).Expanding Eq. ( 19) by f α − f 2 α and using Eq. ( 29), the RE can be obtained in a concise form: Remarkably, Eq. ( 30) allows us to gain a clear quasiparticle picture for the RE, which we will discuss in Sec.VI.In deriving the above formula, we use only C = O[(J/U ) 2 ] and F = O(J/U ) in addition to the relation F i,j = −F j,i , which holds for any state.As long as these conditions are satisfied, therefore, Eq. ( 30) holds for any Gaussian state, regardless of the explicit forms of C and F .We will take advantage of this fact to calculate the RE for the time-evolving state.
In the limit L → ∞, replacing the summations over k with integrals [ k → (L/2π) π −π dk] in Eqs.(24) and (25), we obtain Substituting Eqs. ( 31) and (32) into Eq.( 30), the RE is obtained as The above expression clearly shows that S A is independent of subsystem size L A and therefore follows the arealaw scaling, which is a characteristic feature of a gapped ground state of a short-range Hamiltonian [51].The ground state |vac⟩ indeed satisfies this condition.
Figure 1 shows a comparison of Eq. ( 33) with the numerical results obtained by the matrix-product-state technique.They agree well with each other.To obtain the numerical results, we calculate the ground state wave function of the BHM (1) by imaginary-time evolution using the infinite time-evolving block-decimation (iTEBD) algorithm [52] and evaluate the RE with it.In the numerical calculations using the iTEBD algorithm throughout this paper, we implement it keeping the Schmidt coefficients larger than 10 −7 , setting the dimension of the local Hilbert space being 5, and using the second-order Suzuki-Trotter decomposition with time step ∆t = 0.01/U .

V. TIME-EVOLUTION OF THE R ÉNYI ENTROPY
In this section, we study the time evolution of the RE.We calculate the RE for a single site (L A = 1) in Sec.V A by directly evaluating the reduced density matrix.We further calculate the RE for L A ≥ 1 in Sec.V B using the formalism in Sec.IV.

A. Rényi entropy for a single site
The reduced density matrix for subsystem A can be written in the basis of the Fock states of doublons and holons |n d , n h ⟩ j (n d , n h = 0, 1) as ρA = τ,τ ′ ∈{(0,0),(1,0),(0,1),(1,1)} where n d (n h ) denotes the number of doublon (holon) on the jth site.The matrix elements are given as where 1j denotes the identity operator for all the sites other than j.
The matrix elements can be calculated as where r(t) = ⟨ψ(t)| d † j dj |ψ(t)⟩ = ⟨ψ(t)| ĥ † j ĥj |ψ(t)⟩ represents the number of doublon (holon) per site The calculation of matrix elements r τ,τ ′ is straightforward once we express the operators |τ ⟩ j ⟨τ ′ | j in terms of the creation and annihilation operators of doublons and holons.The detail of the calculation is given in Appendix B. Substituting Eq. (36) into Eq.( 2) and noting that r(t) = O[(J/U ) 2 ], we find that the RE is proportional to the number of doublon (holon) per site in the leading order of J/U as To figure out why the RE is proportional to the number of doublon (holon) per site, we expand the wave function (14) to the first order of J/U as The second term illustrates that all the doublons and holons are excited in entangled pairs [32].Recalling that the RE quantifies entanglement between subsystems A and B, entangled doublon-holon pairs spanning the boundary between subsystems A and B should contribute to the RE.Since the number of doublons (holons) of subsystem A (site j) is equal to that of doublon-holon pairs spanning subsystems A and B, the RE is naturally proportional to the number of doublons (holons).
R s t L n W a 5 r e V W M X w 0 s f X + r 0 q n 2 c X B l + p P z y 5 K W P K 8 q u T d 8 p j W L Z S 2 v n Z 4 8 r q 1 v B l r z L B L 9 k z + L 9 g j u 6 M b G L U 3 5 S r D N 0 8 R o g 9 I / H z u T p C b j y c W 4 s l M M p q K + 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " r 9 < l a t e x i t s h a 1 _ b a s e 6 4 = " P P 5 i d 7 r s l L 7 k y < l a t e x i t s h a 1 _ b a s e 6 4 = " y P a + x g z w j 2. Time evolution of the RE for a single site (LA = 1).Equation ( 40) is plotted as the solid lines.The horizontal lines indicate the asymptotic values SA(t → ∞) = 32(J/U ) 2 .The dotted lines indicate the numerical results of the iTEBD calculations.The inset shows the Fourier transform SA(ω) in Eq. ( 43) for J/U = 0.01.
In the limit L → ∞, evaluating the summation in Eq. ( 37) by replacing it with an integral, we obtain Figure 2 shows the RE in Eq. ( 40) as a function of time.The RE rapidly oscillates soon after the quench and converges after a while.The constant value after the convergence increases as J/U increases.Equation (40) indicates that the frequency of the oscillation is equivalent to U , and The constant after the convergence is given as Figure 2 shows that Eqs. ( 40) and ( 41) agree well with the RE calculated by the iTEBD algorithm.The Fourier transform of Eq. ( 40) is given as where θ(x) is the step function.We assume that the frequency ω is positive to eliminate the contribution from the time-independent term in Eq. (40).We also assume for simplicity that S A (t) = S A (−t) for t < 0. The detail of the derivation of Eq. ( 43) is given in Appendix C.
SA (ω) has a peak at ω = U and its width is ∆ω = 12J, as shown in the inset of Fig. 2. The rapid oscillations with the frequency U are induced by excited bogolons in Eq. ( 14).Their excitation energy can be approximated as The peak position ω = U corresponds to the center of the energy band, while the peak width 12J corresponds to its bandwidth.
B. Rényi entropy for LA ≥ 1 The reduced density matrix for the time-evolving state |ψ(t)⟩ under the quadratic Hamiltonian Ĥeff can be formally written in the form of a thermal state as where ĤA (t) can be written as a quadratic form of dj and ĥj (j ∈ A), because the Bloch-De Dominicis theorem holds for the correlation functions evaluated by ρA (see Appendix D).It means that the time-evolving state is also a Gaussian state, and hence the RE can be calculated using the formalism in Sec.IV A.
Evaluating the single-particle correlation functions in Eqs.(21) and (22) with |ψ(t)⟩, we obtain Note that C i,i = r(t) and F i,i = 0. We find C i,j = O[(J/U ) 2 ] and F i,j = O(J/U ) from Eqs. ( 45) and (46).Thus, the RE for the time-evolving state |ψ(t)⟩ can be written in the same form as Eq.(30).
In the limit L → ∞, evaluating the summation in Eq. (46), we obtain C i,i = r(t) is given in Eq. ( 37).Using Eq.( 30), we obtain the RE for L A ≥ 1 as Setting L A = 1, the above equation indeed reduces to Eq. ( 40).Figures 3 (a) and (b) show S A (t) in Eq. ( 49) as functions of time in the short-time scale for a small subsystem (Jt, L A = O(1)) and long-time scale for a large subsystem (Jt, L A ≫ 1), respectively.In the former, the numerical results of the iTEBD calculations agree well with Eq. ( 49), as shown in Fig. 3 (a).Analogous to Fig. 2, S A (t) exhibits rapid oscillations with the frequency U and converges in the time scale of O(1/J).Meanwhile, in Fig. 3 (b), we find that S A (t) linearly increases for 6Jt < L A and is saturated to a constant proportional to the size of the subsystem L A for 6Jt > L A .The RE after a long time thus obeys the volume-law scaling.We can confirm these behaviors analytically in the asymptotic forms of Eq. ( 49) as The second asymptotic form indicates that the RE approaches to a constant with a correction of order O(L 4 A /t 3 ).Note that the oscillations with frequency U can be seen in the short-time scale even for a large subsystem, as shown in the inset of Fig. 3 (b).

VI. QUASIPARTICLE PICTURE
In this section, we discuss how the RE in the ground state and the time-evolving state can be understood in < l a t e x i t s h a 1 _ b a s e 6 4 = " r 9 < l a t e x i t s h a 1 _ b a s e 6 4 = " b C N / 6 q y W d P s 0 e N k v j h W r X I 6 i < l a t e x i t s h a 1 _ b a s e 6 4 = " G W r 0 Z z y 2 P T 9 J e O q W Z R p n I j M X A i g J/U = 0.01 < l a t e x i t s h a 1 _ b a s e 6 4 = " z m 8 B 9 M T f 7 q l p t H y R S / l d Y s 2 J/U = 0.02 < l a t e x i t s h a 1 _ b a s e 6 4 = " M s L q X D r q o P 1 k m < l a t e x i t s h a 1 _ b a s e 6 4 = " e g b p O K P D e A Y h z k 5 q s K O U l 8 g R J 3 3. Time evolution of the RE for LA ≥ 1.We set J/U = 0.01.Equation ( 49) is plotted as the solid lines.terms of doublon-holon pairs.First of all, recall that the RE for both |vac⟩ and |ψ(t)⟩ is expressed as where we neglect O[(J/U ) ].The above expression can be understood as follows: 2tr(C) = j∈A (⟨n jd ⟩ + ⟨n jh ⟩) represents the total number of doublons and holons in subsystem A, which is equal to the sum of the number of doublon-holon pairs spanning the boundary of subsystem A and twice the number of doublon-holon pairs within subsystem A [See Fig. 4 (a)].Meanwhile, ∥F ∥ 2 F can be written as Here, we neglect a correction of O[(J/U ) 4 ].2∥F ∥ a 6 2 O U 6 P v p U 2 z N N B y o x I 4 W c i + / 6 t q 0 G z j 4 E v 1 p 2 c b V W w 6 X h X y b j h M / x b y Q N 8 6 u u h l t z L R z g q 7 Z q / k / 4 p 1 2 Q P d Q G u 9 y T d p n r m E j z 4 g / v O 5 h 0 F + L R Z f j y X S i U g y 4 X 6 F F 4 t Y x i q 9 9 w a S 2 E E K O T q X 4 x R n O P c 8 C 3 4 h J M w P U g W P q w n h W w h L n / U E i f M = < / l a t e x i t > x < l a t e x i t s h a 1 _ b a s e 6 4 = " L p F U N I 4 w r 5   q 0 G z j 4 E v 1 p 2 c b V W w 6 X h X y b j h M / x b y Q N 8 6 u u h l t z L R z g q 7 Z q / k / 4 p 1 2 Q P d Q G u 9 y T d p n r m E j z 4 g / v O 5 h 0 F + L R Z f j y X S i U g y 4 X 6 F F 4 t Y x i q 9 9 w a S 2 E E K O T q X 4 x R n O P c 8 C 3 4 h J M w P U g W P q w n h W w h L n / U E i f M = < / l a t e x i t > x < l a t e x i t s h a 1 _ b a s e 6 4 = " L p F U N I 4 w r 5 5 a o 1 M M 6 i 4 p Y 5 h n j + y W t d k D u 2 N P 7 P 3 X W k 2 / R s d L g 2 a 1 q + V O J X I 8 n X v 7 V 2 X S L H D w q f r T s 8 A e V n 2 v O n l 3 f K Z z C 6 2 r r x + e t n N r 2 f n m A r t i z + T / k r X Y P d 3 A q r 9 q 1 x m e P U O Y P k D + / t w / w X Y y I S 8 n U p l U f D 0 V f M U g Z j C H R X r v F a x j E 2 n k 6 V w H J z j H R e h F m p J m p V g 3 V Q o F m k l 8 C W n p A 7 a q j e g = < / l a t e x i t > 2tr(C) < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 n E a F p z a z n W C G 1 m 6 d 7 c p d V G w J q 0 9 m q 6 v l q l U 3 T q D i m j G G c P 7 J K 9 s H t 2 x Z 7 Y + 6 + 1 a n 4 N z 0 u V Z q W p 5 X a p / 2 h 4 + e 1 f l U G z w N 6 n 6 k / P A j u Y 8 7 1 q 5 N 3 2 G e 8 W a l N f O T h + W Z 7 P j 9 c m 2 B l 7 J v + n 7 J H d 0 Q 3 M y q t 6 n u P 5 E 4 T p A x L f n / s n W E 3 G E 7 P x V C 4 V y 6 S C r w h h B G O Y p P d O I 4 N F L K F A 5 x 7 i G j e 4 l U J S X J q R 0 s 1 U q S X Q D O F L S J k P 9 U O V A w = = < / l a t e x i t >  51) indicates that they are responsible for entanglement between subsystems A and B, which is consistent with the fact that the doublonholon pairs are entangled.
The area-law scaling of the ground state can be indeed understood within this quasiparticle picture.Given that F i,j is nonzero only when |i − j| = ±1 in Eq. ( 32), doublon-holon pairs spread between nearest-neighbor sites in the ground state.It turns out that only the pairs adjacent to the boundary can straddle the boundary and the number of such pairs does not depend on the size of the subsystem L A .Hence, the RE obeys the area-law scaling.
As for the time evolution of the RE, it can be understood by studying the anomalous correlation function.
< l a t e x i t s h a 1 _ b a s e 6 4 = " G W r 0 Z z y 2 P T 9 J e O q W Z R p n I j M X A i g This propagation speed of doublon-holon pairs has been experimentally confirmed by measuring correlation functions [32,33].These pairs decay in time as shown in Fig. 5 (a) due to the factor 3Jt in the denominator of the second term of Eq. ( 47).
To understand the characteristic features of the longtime dynamics, namely, the linear growth for 6Jt < L A and the following saturation for 6Jt > L A , we focus on the matrix |F i,j (t)| 2 , which is visualized in Fig. 5 (b).One finds that, when 6Jt < L A , the contribution of pairs emitted at t = 0 decreases as the size of pairs grows.It leads to the linear growth of S A (t) for 6Jt < L A .When 6Jt > L A , pairs emitted at t = 0 spread beyond the subsystem size and their contribution to ∥F ∥ 2 F vanishes.This results in the saturation of S A .The smooth transition from the linear growth to the saturation of S A is due to the contribution of the sub-dominant peaks, i.e., pairs emitted at t > 0.
In addition to the propagating peaks, there is a single localized peak at n = 1 in Fig. 5 (a), which corresponds to the subdiagonal elements in Fig. 5 (b).It results in the second term in Eq. ( 49) and its height oscillates with the frequency U .This localized peak represents doublonholon pairs with unit separation spanning the boundary of subsystem A. They are excited by hopping of a boson to the nearest-neighbor sites.They repeat creation and annihilation with the frequency U and eventually decay.
We compare the above quasiparticle picture and the one proposed by Calabrese and Cardy for quench dynamics of a general free-fermion model [16] to highlight our original results.First of all, they proposed a quasiparticle picture for quench dynamics of a general freefermion model, which is only valid in the space-time scaling limit (t, L A → ∞ with L A /t fixed).In contrast, our picture is derived microscopically and not restricted within the space-time scaling limit.In the present work, the quasiparticle picture is derived not only in the spacetime scaling limit, but also in the short-time scale and/or small subsystems.We find, for example, that localized doublon-holon pairs with unit separation yield rapid oscillation of S A in the short-time dynamics.In addition, our quasiparticle picture for the ground state is indeed not included in their picture.Furthermore, our picture has some remarkable features even in the space-time scaling limit that are absent in their picture.First, while the dynamics of entanglement entropy are described in terms of quasiparticle pairs emitted only at t = 0 in their picture, we find that doublon-holon pairs emitted after the initial time also play crucial roles in the dynamics of the RE.In particular, the smooth transition from the linear growth to the saturation of S A can be explained by their presence.In addition, the second asymptotic form of the RE in Eq. ( 50) for Jt ≫ L A can be explained by these pairs emitted at t > 0. Second, we find that the doublon-holon pairs decay as they propagate.This is also absent in their picture.Meanwhile, we confirm their predictions in the space-time scaling limit in our system.Their quasiparticle picture predicts that the entanglement entropy grows linearly up to t ≃ L A /2 (in units where the speed of elementary excitations is unity) and then is saturated to a value proportional to L A [26].This is indeed consistent with Fig. 3 (b) and Eq.(50).
In closing this section, we remark on the possibility of experimental verification of our results on the dynamics of the RE.It may be difficult to experimentally confirm our predictions on the long-time dynamics of the RE in Fig. 3 (b) due to the limitation of the lifetime of an atomic gas and the finite size effect.The short-time dynamics of the RE in Fig. 3 (a) may be verified using the experimental setup in Refs.[30,31].
VII. nTH-ORDER R ÉNYI ENTROPY A /S (2) A L A = 400 L A = 500 L A = 600 FIG. 7. Time evolution of the ratio of the 1st-order RE to the 2nd-order one.We set J/U = 0.01.We note that all lines take almost the same value.
Figure 6 shows the 1st-order RE for the ground state |vac⟩ as a function of L A . S A is almost independent of L A , which implies that the 1st-order RE follows the arealaw scaling.Figure 7 shows the ratio of the 1st-order RE to the 2nd-order RE as a function of time.The ratio is almost constant in time and shows very little dependence on subsystem size L A .The RE in arbitrary order n ≥ 1 thus exhibits qualitatively the same behavior as the 2ndorder RE in the dynamics as well as in the ground state.Its behavior indeed reflects the quasiparticle picture described in Sec.VI.

VIII. VALIDITY OF THE EFFECTIVE THEORY
We have studied the RE within the effective theory.The effective theory ( 5) is integrable whereas the original BHM (1) and the effective Hamiltonian with the projection operator (4) are non-integrable.This difference may induce deviations in REs in the region away from the deep MI regime.Therefore, it is important to clarify how large J/U can be for the effective theory to be qualitatively valid.To this end, we compare the REs calculated by the effective theory and the iTEBD algorithm.
Figure 8 (a) shows a comparison of their time evolution.It clearly shows that the effective theory overestimates the RE.Furthermore, the overestimation increases in time during the linear growth and stops increasing after it.This overestimation may arise from the approximation in deriving the effective Hamiltonian (5), in which we ignore the projection operator and allow unphysical double occupancy of a doublon and a holon.This approximation may result in overestimation of the number of doublon-holon pairs and accordingly the RE, because doublon-holon pairs are responsible for the RE.The overestimation of the number of doublon-holon pairs is expected to increase during the linear growth of the RE, because doublon-holon pairs keep generating in this period as we observe in Sec.VI.This may result in the increase of the overestimation in the linear growth part in Fig. 8 (a).Note that the projection operator introduces non-integrable effects that correspond to interactions of bogolons.
Figure 8 (a) also shows that the overestimation increases as J/U increases.To understand it quantitatively, we calculate time-averaged overestimation where S eff A and S iTEBD A denote the RE obtained by the effective theory and iTEBD calculations, respectively, and M = t max /∆t the number of data points in the interval 0 ≤ t ≤ t max .We set U t max = 20.Figure 8 (b) shows a log-log plot of ∆S A versus J/U .Since the time-averaged overestimation scales as ∆S A ∝ (J/U ) 3.2-3.6, the effective theory is valid as long as (J/U ) 3 ≪ 1.Since our quasiparticle picture is valid within the effective theory, it is also valid under the same condition.This condition is consistent with a previous study based on the analysis of a density correlation function [47], in which it is claimed that the effective theory is valid as long as J/U ≤ 1/8.Note that we confirmed that the power of J/U is larger than 3 regardless of the values of ∆t and/or t max .

IX. SUMMARY
In summary, we have investigated the time evolution of the RE for bosons in a one-dimensional optical lattice when the system undergoes a quench from the MI limit to the strongly correlated MI regime.Developing the effective theory, we have derived a direct relation between the RE and correlation functions associated with doublon and holon excitations.Using this relation, we have calculated the RE analytically and obtained a physical picture, both in the ground state and during time evolution, in terms of entangled doublon-holon pairs.Our quasiparticle picture for the dynamics of the RE introduces t e x i t s h a 1 _ b a s e 6 4 = " H G N I R Z P 4 y k o 9 u d L k q 2 0 B 3 P 1 m F 1 0 r M X u 2 S 1 7 Y u + / 1 m p 4 N d p e 6 j T L H S 2 3 i u G j i c 2 3 f 1 U 6 z S 4 O P l V / e n Z R w p L n V S X v l s e 0 b 6 F 0 9 L X D k 9 b m c n a m M c s u 2 T P 5 v 2 C P 7 I 5 u Y N R e l a s M z 5 4 i R B 8 Q / / 7 c P 0 F t e x i t s h a 1 _ b a s e 6 4 = " H G N I R Z P 4 y k o 9 u d L k q 2 0 B 3 P 1 m F 1 0

L A = 6 <
r M X u 2 S 1 7 Y u + / 1 m p 4 N d p e 6 j T L H S 2 3 i u G j i c 2 3 f 1 U 6 z S 4 O P l V / e n Z R w p L n V S X v l s e 0 b 6 F 0 9 L X D k 9 b m c n a m M c s u 2 T P 5 v 2 C P 7 I 5 u Y N R e l a s M z 5 4 i R B 8 Q / / 7 c P 0 F + I R Z P x h K Z R D Q V 8 7 8 i i E l M Y 4 7 e e x E p r C O N H J 1 r 4 R h n O A + 8 C O P C l C B 2 U o W A r x n D l x D m P w B a v I 2 1 < / l a t e x i t > l a t e x i t s h a 1 _ b a s e 6 4 = " m y W T J x O + a z Z 8 8 j v T 9 d s M U 1 8 d / T 4 1 b o F I 2 6 T c o I Y u y R 3 b I 2 e 2 B 3 7 J l 9 / F q r 5 d X o e G n S L H e 1 3 C q H j m e 3 3 / 9 V 6 T S 7 O P h S / e n Z R Q W r n l e V v F s e 0 7 m F 0 t U 3 D k / b 2 2 u Z W C v O r t g L + b 9 k T + y e b m A 0 3 p T r N M + c I U g f I P 5 8 7 l 6 w u 5 Q U l 5 O p d C q 6 n v K / Y h h z W E C C 3 n s F 6 9 j E F r J 0 r o 4 T n O M i 8 C r M C H P C f D d V C P i a a X w L I f 4 J S P m M 0 Q = = < / l a t e x i t > (b) < l a t e x i t s h a 1 _ b a s e 6 4 = " P P 5 i d 7 r s l L 7 k y A x f 8 M 6 Z h b U g x 8 o 4 k u y e X b M X d s d u 2 C N 7 / 7 V X I + j h e 6 l T 1 l t a 7 u w O H o 9 u v v 2 r q l C W 2 P 9 U / e l Z o o S F w K s g 7 0 7 A + K c w W v r a 4 e n L 5 u J G s j H F L t g T + T 9 n D + y R U W 9 I 7 p k z 9 9 y c u a M 5 p u E J x p o h q a u 7 p 7 c v 3 B 8 Z G B w a H o m O j m 1 5 d s X V e V 6 3 T d s t a K r H T c P i e W E I k x c c l 6 t l z e T b 2 t F a e 3 + 7 y l3 P s K 1 N U X N 4 s a w e W M a + o a u C q M L G X H 6 V y W y + F E 1 Q 8 i P + E y g B S C C I t B 2 9 w S 7 2 Y E N H B W V w W B C E T a j w 6 N m B A g a H u C L q x L m E D H + f 4 w Q R 0 l a o i l O F S u w R f Q 9 o t R O w F q 3 b P T 1 f r d N f T H p d U s a R Z I /s l r X Y A 7 t j T + z 9 1 1 5 1 v 0 f b S 4 2 y 1 t F y p z R y O p 5 7 + 1 d V p i x w + K n 6 0 7 P A P p Z 9 r w Z 5 d 3 y m f Q q 9 o 6 8 e n 7 d y K 9 l k f Y Z d s W f y f 8 m a 7 J 5 O Y F V f 9 e s M z z R H T X o k e 7 p h T 5 + r V U L a j S 9 n P C s t r T C K c X O x r b e / 1 W Z P P s 4 / l L 9 6 d n H I R Y C r z p 7 d w K m e Q u t p a + e n j e 2 F r N T t W m 6 p l f 2 f 0 X P 9 M A 3 s K p v 2 s 2 m y F 4 i y p 1 U e c h I s 4 t S g S p d l d U y b d 7 m 5 2 p 0 3 Y 1 A / w B x y c C B E R / o S L P + D g J 4 g j i Y u E t 9 t N B M E 7 m Z l n n n m f d 5 6 Z 0 W x D u J K x h 5 D S 1 N z S 2 t b e E e 7 s 6 u 7 p j f T 1 r 7 p W x d F 5 T r c M y 1 FIG. 3. Time evolution of the RE for LA ≥ 1.We set J/U = 0.01.Equation (49) is plotted as the solid lines.(a) REs for LA = 2, 4, and 6 in the short-time scale.The numerical results of the iTEBD algorithm are indicated by the dotted lines.The dashed lines indicate 16(J/U ) 2 (LA + 1).(b) for = 400, 500, and 600 in the long-time scale.The dotted and dashed lines indicate the asymptotic forms of the in Eq. (50).The inset shows a magnification in the short-time scale for LA = 400.

2 F 2 F
is equal to twice the number of pairs within subsystem A [See Fig. 4 (b)].2tr(C) − 2∥F ∥ is thus equal to the < l a t e x i t s h a 1 _ b a s e 6 4 = " M 4 e 5 h 2 P 8 2 f 3 O W 4 Z R Y t s M D E y V 9 C 4 = " > A A A C Z H i c h V F N S w J B G H 7 c v s w s L Q m C I C Q x O s k Y U t F J 6 N L R j / w A E 9 n d R l t c d 5 f d V T L p D 9 S 1 6 N C p I C L 6 G V 3 6 A x 3 8 A 0 F 0 N O j S o d d 1 I U q q d 5 i Z Z 5 5 5 n 3 e e m Z E M V b F s x r o e Y W R 0 b H z C O + m b 8 k / P B I K z c 3 l L b 5 o y z 8 m 6 q p t F S b S 4 q m g 8 Z y u 2 y o u G y c W G p P r 3 7 I q r 8 Y J m G 7 Z b U h W P G 7 r F C 0 I X B i 8 5 L l d M 1 e B F 9 W C p u V + s c t f T b W t N 1 B y + b S p 7 l r 6 r a 4 o g a m O l 7 G + 5 Z m K x X o 4 n W Y o F k f g J 5 B A k E U b W j l 9 j C z u w o a E C E x w W B G E D C j x q m 5 D B 4 B C 3 D Z 8 4 l 5 A e 7 H P U E S N t h b I 4 Z S j E H t C 4 R 6 v N k L V o 3 a z p B W q N T j G o u 6 R M Y J I 9 s B v 2 w u 7 Z L X t k 7 7 / W 8 o M a T S 8 1 m t W W l j v l g a P R / N u / K p N m g f 1 P 1 Z + e B X Y x H 3 j V y b s T M M 1 b a C 1 9 9 f D k J b + w O u l P s U v 2 R P 4 v W I P d 0 Q 2 s 6 q t 2 l e O r p 4 j R B 8 j f n / s n W J 9 J y b O p d C 6 d z K T D r 4 h i D B O Y p v e e Q w b L y K J A 5 5 o 4 x h n O I 8 / S i D Q m j b d S p U i o G c a X k K Y + A L v c j Q o = < / l a t e x i t > L A < l a t e x i t s h a 1 _ b a s e 6 4 = " M 4 e 5 h 2 P 8 2 f 3 5 5 5 n 3 e e m Z E M V b F s x r o e Y W R 0 b H z C O + m b 8 k / P B I K z c 3 l L b 5 o y z 8 m 6 q p t F S b S 4 m g 8 Z y u 2 y o u G y c W G x H 3 j V y b s T M M 1 b a C 1 9 9 f D k J b + w O u l P s U v 2 R P 4 v W I P d 0 Q 2 s 6 q t 2 l e O r p 4 j R B 8 j f n / s n W J 9 J y b O p d C 6 d z K T D r 4 h i D B O Y p v e e Q w b L y K J A 5 5 o 4 x h n O I 8 / S i D Q m j b d S p U i o G c a X k K Y + A L v c j Q o = < / l a t e x i t > L A < l a t e x i t s h a 1 _ b a s e 6 4 = " M 4 e 5 h 2 P 8 2 f 3 O W 4 Z R Y t s M D E y V 9 C 4 = " > A A A C Z H i c h V F N S w J B G H 7 c v s w s L Q m C I C Q x O s k Y U t F J 6 N L R j / w A E 9 n d R l t c d 5 f d V T L p D 9 S 1 6 N C p I C L 6 G V 3 6 A x 3 8 A 0 F 0 N O j S o d d 1 I U q q d 5 i Z Z 5 5 5 n 3 e e m Z E M V b F s x r o e Y W R 0 b H z C O + m b 8 k / P B I K z c 3 l L b 5 o y z 8 m 6 q p t F S b S 4 q m g 8 Z y u 2 y o u G y c W G p P K C V N / u 7 x d a 3 L Q U X d u 1 2 w Y v N 8 S a p l Q V W b S J S h 9 W g h E W Y 0 6 E h 0 H c B R G 4 k d K D t 9 j D P n T I a K I B D g 0 2 Y R U i L G o l x M F g E F d G h z i T k O L s c x z D R 9 o m Z X H K E I m t 0 1 i j V c l l N V r 3 a 1 q O W q Z T V O o m K c O I s i d 2 x 3 r s k d 2 z F / b x a 6 2 O U 6 P v p U 2 z N N B y o x I 4 W c i + / 6 t T F m r x B 9 y 5 E n S l I C J + h h t / w E U / Q X R X w Y 0 L b 6 c D o q L e k O T k 5 J 6 b k 0 R z D O F J x l o h p a e 3 r 3 9 g c C g 8 P D I 6 F o m O T x Q 8 u + b q P K / b h u 2 W N N X j h r B 4 X g p p 8 J L j c t X U D F 7 U q h u d / W K d u 5 6 w r W 3 Z c T F m r x B 9 y 5 E n S l I C J + h h t / w E U / Q X R X w Y 0 L b 6 c D o q L e k O T k 5 J 6 b k 0 R z D O F J x l o h p a e 3 r 3 9 g c C g 8 P D I 6 F o m O T x Q 8 u + b q P K / b h u 2 W N N X j h r B 4 X g p p 8 J L j c t X U D F 7 U q h u d / W K d u 5 6 w r W 3 Z c / l a t e x i t > B < l a t e x i t s h a 1 _ b a s e 6 4 = " N T g X 2 e o F H J 9 3 G R W g r M s c 8 b f nf 2 s = " > A A A C b H i c h V G 7 S g N B F D 1 Z X z G + 4 q M Q R A g G x S r M S l C x i t h Y W C T G m I h K 2 F 1 H X d w X u 5 N A X P I D t h Y W a q E g I n 6 G j T 9 g k U 8 Q w U b B x s K b z Y K o q H e Y m T N n 7 r l z Z k Z 1 D N 0 T j D U i U l t 7 R 2 d X t D v W 0 9 v X P x A f H Fr 3 7 I q r 8 Y J m G 7 Z b U h W P G 7 r F C 0 I X B i 8 5 L l d M 1 e B F 9 W C p u V + s c t f T b W t N 1 B y + b S p 7 l r 6 r a 4 o g a m O l 7 G + 5 Z m K x X o 4 n W Y o F k f g J 5 B A k E U b W j l 9 j C z u w o a E C E x w W B G E D C j x q m 5 D B 4 B C 3 D Z 8 4 l 5 A e 7 H P U E S N t h b I 4 Z S j E H t C 4 R 6 v N k L V o 3 a z p B W q N T j G o u 6 R M Y J I 9 s B v 2 w u 7 Z L X t k 7 7 / W 8 o M a T S 8 1 m t W W l j v l g a P R / N u / K p N m g f 1 P 1 Z + e B X Y x H 3 j V y b s T M M 1 b a C 1 9 9 f D k J b + w O u l P s U v 2 R P 4 v W I P d 0 Q 2 s 6 q t 2 l e O r p 4 j R B 8 j f n / s n W J 9 J y b O p d C 6 d z K T D r 4 h i D B O Y p v e e Q w b L y K J A 5 5 o 4 x h n O I 8 / S i D Q m j b d S p U i o G c a X k K Y + A L v c j Q o = < / l a t e x i t > L A < l a t e x i t s h a 1 _ b a s e 6 4 = " h b d K q + 5 t E 2 A W z e W / F Z z 9 i B H 3 b d 0 7 O P x S / e n Z w z 7 m A 6 + C v N s B 4 9 9 C a + o r x 2 e N / G J u v D r B r t k L + b 9 i d f Z A N z A r b 9 p N l u f O E a M P k H 8 + d y t Y n 8 7 I s 5 m Z 7 E x 6 a S H 8 i i i S G M M k v f c c l r C C V R T o X A O n u M B l 5 F U a l p L S a D N V i o S a I X w L a e I T S H i M 1 Q = = < / l a t e x i t > (a) < l a t e x i t s h a 1 _ b a s e 6 4 = " Q e z g K B / t 0 d W V C E 5 x 1 8k P r p X F a W M = " > A A A C b H i c h V G 7 S g N B F D 1 Z X z E + E h + F E I R g U L Q J s y K + K s H G 0 i T G B y p h d x 1 1 y b 7 Y n Q Q 0 5 A d s L S z U Q k F E / A w b f 8 A i n y C C T Q Q b C + 9 u F k S D e o e Z O X P m n j t n Z l T H 0 D 3 B W D 0 i t b V 3 d H Z F u 2 M 9 v X 3 9 8 c T A 4 L p n l 1 2 N F z T b s N 1 N V f G 4 o V u 8 I H R h 8 E 3 H 5 Y q p G n x D L S 3 7 + x s V 7 n q 6 b a 2 J I 4 f v m s q B p e / r m i K I 2 p q s 7 r h m S q 1 N F R N p l m F B p F q B H I I 0 w l i 1 E 7 f Y w R 5 s a C j D B I c F Q d i A A o / a N m Q w O M T t o k q c S 0 g P 9 j l q i J G 2 T F m c M h R i S z Q e 0 G o 7 Z C 1 a + z W 9 Q K 3 R K Q Z 1 l 5 Q p j L M nd s c a 7 J H d s 2 f 2 8 W u t a l D D 9 3 J E s 9 r U c q c Y P x n J v / + r M m k W O P x S / e l Z Y B / z g V e d v D s B 4 9 9 C a + o r x 2 e N / G J u v D r B r t k L + b 9 i d f Z A N 7 A q b 9 p N l u f O E a M P k H 8 + d y t Y n 8 7 I s 5 m Z 7 E x 6 a S H 8 i i i S G M M k v f c c l r C C V R T o X B O n u M B l 5 F U a l p L S a D N V i o S a I X w L a e I T S n q M 1 g = = < / l a t e x i t > (b) < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 8 0 X 6 J 5 h R L 8 K 3 M 1 S c 4 B i S i a m n 6 I= " > A A A C b H i c h V G 7 S g N B F D 1 Z 3 / E V H 4 U Q h G B Q t A m z E n x V A R v L P I w P V M L u O O r i v t j d B G L w B 2 w t L N R C Q U T 8 D B t / w C K f I I J N B B s L 7 2 4 W R I N 6 h 5 k 5 c + a e O 2 d m V F v X X I + x e k R q a + / o 7 O r u i f b 2 9 Q 8 M x o a G 1 1 y r 7 H B R 5 J Z u O R u q 4 g p d M 0 X R 0 z x d b N i O U A x V F + v q 4 b K / v 1 4 R j q t Z 5 q p X t c W O o e y b 2 p 7 G F Y + o z e n a t m M k + P F M K Z Z k K R Z E o h X I I U g i j K w V u 8 U 2 d m G B o w w D A i Y 8 w j o U u N S 2 I I P B J m 4 H N e I c Q l q w L 3 C M K G n L l C U o Q y H 2 k M Z 9 W m 2 F r E l r v 6 Y b q D m d o l N 3 S J n A J H t i d 6 z B H t k 9 e 2 Y f v 9 a q B T V 8 L 1 W a 1 a Z W 2 K X B k 7 H C + 7 8 q g 2 Y P B 1 + q P z 1 7 2 M N C 4 F U j 7 3 b A + L f g T X 3 l 6 K x R W M p P 1 q b Y N X s h /1 e s z h 7 o B m b l j d / k R P 4 c U f o A + e d z t 4 K 1 2 Z Q 8 l 0 r n 0 s n M Y v g V 3 Y h j A t P 0 3 v P I Y A V Z F O l c A 6 e 4 w G X k V R q V 4 t J 4 M 1 W K h J o R f A t p 6 h N M f I z X < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " T y V a W / v E 8 h e t Y H p e R v g H s r P 1 B 0 4 8 A E i 4 m e 4 8 Q d c 9 B P E l V R w 4 8 L b 6 Y C o q D c k O T m 5 5 + Y k U R 1 D 9 w R j r Z D U 0 9 v X P z A 4 F B 4 e G R 2 L R M c n t j 2 7 5 m o 8 r 9 m G 7 R Z U x e O G b v G 8 0 I X B C 4 7 L F V M 1 + I 5 a 3 e j s 7 9 S 5 6 + m 2 t S U a D i + b y r 6 l 7 + m a I o g q J Z s l 1 4 w J 9 2 h x Y 6 k S j b M E 8 y P 2 E 8 g e H e z I C r q H W b m z J l 7 7 p y Z k U x V s R 3 G 2 g G h q 7 u n t 6 9 / I D g 4N D w S C o + O 5 W 2 j Z s k 8 J x u q Y W 1 L o s 1 V R e c 5 R 3 F U v m 1 a X N Q k l R e k g z V 3 v 1 D n l q 0 Y + p Z z a P K S J l Z 1 p a L I o k N U O T y R b B a n m + v u U G 4 U L S 2 y 3 t p N l s N R F m d e R H 6 C h A + i 8 G P T C F + h i D 0 Y k F G D B g 4 d D m E V I m x q O 0 i A w S S u h A Z x F i H F 2 + d o I U j a G m V xy h C J P a C x S q s d n 9 V p 7 d a 0 P b V M p 6 j U L V J G E G M P 7 J o 9 s 3 t 2 y x 7 Z + 6 + 1 G l 4 N 1 8 s h z V J H y 8 1 y 6 H g y + / q v S q P Z w f 6 n 6 k / P D i p Y 9 r w q 5 N 3 0 G P c W c k d f P z p 5 z q 5 k Y o 1 Z d s G e y P 8 5 a 7 M 7 u o F e f 5 E v 0 z x z i i B 9 Q O L 7 c / 8 E + W Q 8 s R h P p V P R 1 Z T / F f 2 Y w g z m 6 L 2 X s I o N b C J H 5 x 7 h D N e 4 C b w J 0 k 2 p e L U p I k 4 a i k S p d l d g 4 3 9 y u 6 0 S b U 9 S / w B B y c S E R F X f o C L P + D g J 4 g j i Y u D d 7 e b C I J 3 M j P P P P M + 7 z w z o 9 i 6 5 g r G H l u k 1 r b 2 j s 5 Q V 7 i 7 p 7 e v P z I w u O p a Z U f l B d X S L W d d k V 2 u a y Y v C E 3 o f N 1 2 u G w o O l 9 T 9 r P e / l q F O 6 5 m m S u i a

< l a t e x i t s h a 1 _ b a s e 6 4 = " 1
6 h 1 m 5 s y Z e + 6 c m Z E t T X V c x p 4 C Q l / / w O D Q 8 E h w d G x 8 I h S e n N p 1 z L q t 8 K x i a q a d l y W H a 6 r B s 6 7 q a j x v 2 V z S Z Y 3 n 5 N p G Z z / X 4 L a j m s a O 2 7 R 4 S Z e q h l p R F c k l a i / R K t p 6 R D 5 a L I e j L M m 8

Figure 5 (
Figure 5 (a) shows |⟨ψ(t)|d i h i+n |ψ(t)⟩| 2 = |F i,i+n (t)| 2 as a function of n.The propagating peaks in the figure correspond to those of the Bessel function J n (6Jt) (1 ≤ n ≤ L A − 1) in Eq. (47).Each of them describes a wave packet of entangled doublon-holon pairs induced by the quench.The most dominant peak of J n (6Jt) at n ≃ 6Jt represents pairs emitted at the initial time t = 0.The sub-dominant propagating peaks represent pairs emitted at t > 0. All the peaks propagate with the same velocity v pair = 6J, which coincides with the maximum group velocity of a doublon-holon pair with opposite momenta given by

FIG. 6 .
FIG.6.1st-order RE for the ground state |vac⟩ as a function LA.The crosses represent the results obtained from Eq. (60), where fα are obtained by numerically diagonalizing the matrix M .The squares represent the numerical results of the iTEBD calculation.