Systematic strong coupling expansion for out-of-equilibrium dynamics in the Lieb-Liniger model

We consider the time evolution of local observables after an interaction quench in the repulsive Lieb-Liniger model. The system is initialized in the ground state for vanishing interaction and then time-evolved with the Lieb-Liniger Hamiltonian for large, finite interacting strength $c$. We employ the Quench Action approach to express the full time evolution of local observables in terms of sums over energy eigenstates and then derive the leading terms of a $1/c$ expansion for several one and two-point functions as a function of time $t>0$ after the quantum quench. We observe delicate cancellations of contributions to the spectral sums that depend on the details of the choice of representative state in the Quench Action approach and our final results are independent of this choice. Our results provide a highly non-trivial confirmation of the typicality assumptions underlying the Quench Action approach.


Introduction
The non-equilibrium dynamics in isolated many-particle quantum systems has attracted a great deal of attention over the last decade [1][2][3][4][5][6]. These developments were driven by the ability to realize almost isolated many-particle quantum systems using trapped, ultra-cold atoms and investigate their time evolution when driven out of equilibrium in exquisite detail, see e.g. Refs [7][8][9][10][11][12][13][14][15][16][17]. It was realized early on that conservation laws play a crucial role in the late time relaxational behaviour of isolated systems [8,18]. This implies in particular that in the thermodynamic limit integrable systems with extensive numbers of conservation laws will typically relax to non-thermal stationary states . The full time evolution of local observables in integrable models is equally interesting, but significantly harder to determine. Early work focused on rational conformal field theories [5,40,41] and non-interacting models [19,[42][43][44]. The low density regime after weak quantum quenches has be analyzed by means of linked-cluster expansions [42,[45][46][47][48] and semiclassical methods [49][50][51]. Arguably the method of choice for studying the time evolution of local operators in interacting integrable models is the so-called Quench-Action approach [23,52]. To date it mostly has been applied to determine and characterize the stationary state [26][27][28][29]32,[34][35][36][37]. Exceptions are Refs [46], [53] and [54], which address respectively the asymptotic late-time regimes after quenches to the sine-Gordon, Lieb-Liniger and transverse field Ising models respectively. According to the Quench-Action approach the expectation values of local operators after a quantum quench from an initial state |Ψ are given by Here L denotes the system size, the so-called representative state |Φ s is a simultaneous eigenstate of the Hamiltonian and of the (quasi)local [33] conservation laws I (n) of the theory under consideration, such that it correctly reproduces the extensive parts of the expectation values of the I (n) in the initial state The structure of (1) is similar to that of response functions in equilibrium and provides a spectral representation in terms of (normalized) energy eigenstates |n by writing Ψ|O (t) |Φ s = n Ψ|n n |O (0)| Φ s e it(En−Es) .
In practice the Quench Action approach faces two challenges: • It requires knowledge of the overlaps Ψ|n between the initial state and energy eigenstates. This is known as the "initial state problem". To date such overlaps are known for a number of specific examples only [55][56][57][58][59][60][61][62], but many of these are physically interesting.
• Determining the time evolution requires carrying out spectral sums like (3). Given that these generally involve an exponentially (in system size) large number of terms this is a formidable challenge.
In this work we focus on the second of these problems, namely how to extract the time dependence of local observables after a quantum quench from the spectral representation. We consider the case of a quantum quench to the repulsive Lieb-Liniger model, and bring to bear strong-coupling expansion methods we recently developed in the context of equilibrium response functions [63].
In the following we set = 2m = 1, impose periodic boundary conditions and restrict ourselves to the repulsive case c > 0. For later convenience we define the local operators of interest, namely the density operator at position x and the interaction potential The Lieb-Liniger model is solvable by the Bethe ansatz [64][65][66]. Its eigenfunctions can be parametrized by N rapidity variables λ 1 , ..., λ N that on a ring of radius L satisfy a set of quantization conditions known as "Bethe equations" 1 π arctan λ k − λ j c , k = 1, . . . , N.
Here I k are integer if N is odd and half-odd integer if N is even. The corresponding eigenstate |λ λ λ can be written as |λ λ λ = B(λ 1 )...B(λ N )|0 , where B(λ) is a creation operator acting on a particular reference state |0 . The eigenvalues of the Hamiltonian and other conserved quantities are expressed in terms of the rapidities as well. For example the energy E(λ λ λ) and momentum P (λ λ λ) read For c > 0 all the solutions λ i to the Bethe equations are real [66].

Quench protocol and observables of interest
Following [28] we consider the following quantum quench protocol. We assume that the system is prepared in the Bose-Einstein condensate (BEC) ground state for N particles in the absence of interactions At t = 0 we then suddenly turn on the interactions, so that for t > 0 the time evolution of the system |Ψ(t) is governed by the Hamiltonian (4) Our aim is to determine the full time evolution of a number of different observables after the quench in the framework of the systematic 1/c-expansion developed in [63]. We have considered the following one and two-point functions: • One-point function of the interaction potential • Density-density correlation function • Steady-state expectation value of the two-point function of the interaction potential Here we have defined σ 2 (x, τ ) = e iHτ σ 2 (x)e −iHτ . We note that we use a different notation for the time difference τ to avoid confusion with the time t according to which the system evolves after the quench. The analogous two-point function for the density operator was derived in [63] up to order 1/c 2 .

Summary of results
As the derivations of our results are quite technical we start by presenting our final answers and discuss their physical implications. All correlators are expressed in terms of distribution functions of particles ρ(λ) and holes ρ h (λ) defined as follows [28] where τ = D c and with I the modified Bessel function. The particle density D is related to ρ(λ) by 2.1 Relaxation of the one-point function σ 2 (0) t Our final result for the time evolution of the interaction potential σ 2 (0) after the quench, valid at all finite times t > 0 and expanded in 1/c up to and including order O(c −4 ) is The steady-state value σ 2 (0) ∞ has been previously calculated in [28]. From (18) the latetime asymptotics can be straightforwardly extracted with a saddle point approximation Here ρ (0) denotes the second derivative of ρ(λ) evaluated at λ = 0. The asymptotic t −3 dependence is in agreement with a previous conjecture [53]. However, our results show that this regime is reached only at rather late times when the expectation value is already negligibly small. This is shown in Figure 1, where we plot Our 1/c-expansion provides us with the first few terms of an expression of the form g 2 (t) = The dotted red line is the result for the leading asymptotics ∝ t −3 . The inset shows the same quantities on a logarithmic scale. ∞ n=2 γ −n a n (t), where γ = c D and the functions a n (t) incorporate non-perturbative summations of certain terms at all orders in 1/c . In order to assess the parameter range in which the series may be convergent we consider the ratios r n (t) = an(t) a 2 (t) 1 n−2 , n = 3, 4.
In Figure 2 we plots these ratios as functions of t for c = 3 and D = 0. 16. We see that both 5 10 15 2 4 6 8 t r 3,4 (t) Figure 2: r 3 (resp. r 4 ) as a function of t, in light (resp. dark) green. These ratios give an estimate of the smallest value of γ for which the series is convergent.
ratios grow at short times, indicating that the series is not likely to be uniformly convergent near t = 0. Moreover, it follows from the fact that g 2 (0) = 1 while g 2 (t) = O(c −2 ) for all t > 0 that a resummation of the series is required to capture limit t → 0. For t 1 the results for r 3,4 (t) suggest that the series could be convergent for γ 4. As a comparison, we recall that the series in γ for the ground state energy density is convergent for γ > 4.527 [67].

Relaxation of the two-point function σ(x)σ(0) t
We find that the leading contributions in the 1/c-expansion of the density-density correlation function can be cast in the form where and Here − denotes a principal value integral defined as The limit c → ∞ of (22) was previously computed in [28]. The density-density correlator (22) is shown in Figs 3 and 4. We see that for the chosen parameters D = 1 and c = 10 the effects of the O(c −1 ) term are clearly visible and significantly modify the c = ∞ result. In particular the oscillatory behaviour as a function of distance for short times becomes more pronounced for smaller values of c. Perhaps the most striking feature of Fig. 3 is the apparent absence of any light cone effect [41]. This can be understood by noting that (i) our initial state has an infinite correlation length and any light cone like feature would therefore be weak; (ii) the local Hilbert space is infinite dimensional and the dispersion relation of elementary excitations unbounded. Hence the Lieb-Robinson bound [68] does not apply and "superluminal" effects [69] are allowed.
An alternative representation of (22) more suitable for numerical evaluations and an analysis of the x → 0 and t → 0 limits are presented in Appendix D.
The large x and t asymptotics of (22) at fixed ratio α = x 4t can be determined by a stationary phase approximation, which results in with α = x 4t .

c in the stationary state
We discussed how to determine the non-equal time density-density correlation function in an arbitrary energy eigenstate described by a root density ρ(λ) in our previous work [63]. The results in this Section are thus valid for a generic root density ρ, the steady state one (15) being a particular case. Applying the same method to the connected dynamical two-point function of σ 2 (x) gives the following result where we have defined The result for the connected two-point function in the stationary state reached at late times after the quench is obtained by substituting the particle and hole densities (15) into (27) and (28). The leading asymptotic behaviour for large x and τ with α = x 2τ kept fixed can be obtained by a stationary phase approximation These results can be compared with predictions of Generalized Hydrodynamics [70]. According to these the leading large time and distance asymptotics of connected correlations between two local observables O 1 and O 2 is where V O (λ) is the so-called hydrodynamic projection of the operator O, and v eff (λ) the effective velocity associated with the macro-state defined by the particle and hole densities ρ(λ) and ρ h (λ). The hydrodynamic projection V σ 2 (λ) of σ 2 has been determined in [70] Here h n (λ) = λ n−1 , g µ (λ) = µ−λ (λ−µ) 2 +c 2 , and the dressing operation h dr is defined by We find that the asymptotics (29) agrees with this GHD prediction at leading order in 1/c. The dynamical two-point function of σ 2 (x) is related to the Drude weight D and the Onsager coefficient L by In contrast to the density-density correlator the two-point function of σ 2 (x) is expected to exhibit diffusive behaviour, i.e. have a non-vanishing Onsager coefficient L = 0. Our expression for the two point function translates into the following results for D and L This shows that higher orders in the 1/c-expansion are required to determine the Onsager coefficient. We note that this specific result holds only for root densities ρ that decay sufficiently fast at infinity. In the case of the steady state root density (15), because of the slow decay of the density, the next terms in the 1/c expansion should be re-summed to yield a convergent integral.

Quench Action approach and 1/c expansion
In this Section we discuss the implementation of a 1/c expansion of the Quench Action approach [23], which we will then apply to several observables of interest in the remainder of the paper.
Here we have assumed that Ψ BEC |Ψ BEC = 1. The Quench Action approach [23] posits that one of the two sums in (35) is completely dominated by states around a saddle point characterized by a certain root distribution ρ s that is fixed by the overlaps. This allows one to rewrite (35) in the form i.e. a generalized micro-canonical average [23,71] over a set S L of microstates corresponding to the root density ρ s . Employing typicality ideas the micro-canonical average is then replaced by the expectation value with respect to a single "representative state" |λ λ λ [23] lim We note that this last step implies that in the thermodynamic limit (37) depends on the representative state |λ λ λ only through its root density ρ(λ).

"Initial data" for the quench protocol of interest
To be of practical use the representation (37) requires closed-form expressions for the overlaps Ψ BEC |λ λ λ . For our quench protocol an efficient representation for the overlaps was derived in [27,28]. Importantly, the overlaps are non-zero only for "pair" states, i.e. states whose rapidities are of the form −λ N/2 , ..., −λ 1 , λ 1 , ..., λ N/2 with 0 < λ j ∀j. We will denote a set of positive λ j by λ λ λ > 0. We will use the notationλ λ λ = (−λ λ λ) ∪ λ λ λ for such sets of rapidities. The overlaps are then given by where G ± (λ λ λ) are (N/2) × (N/2) matrices of the form For our quench protocol the saddle point root distribution was determined in Ref. [28] and is given in (15).

The 1/c expansion
Our objective is to combine the Quench Action approach to non-equilibrium dynamics (36) with a strong coupling expansion around c = ∞. A detailed exposition of the 1/c expansion technique for dynamical correlation functions in equilibrium has been given in [63]. In the following we recall the key steps of the method and then extend it to the out-of-equilibrium case.
In order to facilitate the 1/c-expansion of the form factors and Bethe equations we first fix an arbitrary, large Λ > 0 that will be sent to ∞ at the end of the calculation. We then select an arbitrary averaging state λ λ λ by fixing its Bethe numbers I I I, impose the constraint that ∀i, |λ i | < Λ, and define the following overlap-weighted spectral sum The overlap-weighted form factor can then be expanded in powers of 1/c at fixed L, I I I where J J J denotes the Bethe numbers of µ µ µ. We also expand the argument of the phase but do not expand the phase e it(E(λ λ λ)−E(µ µ µ)) itself in powers of 1/c. The truncation of the resulting series at a given order O(c −m ) defines the m-th term of our expansion. Once this truncation has been done, the thermodynamic limit and (if necessary) the average in (36) can be performed. By construction, the result depends only on the root density ρ of the fixed averaging state λ λ λ. Finally one would like to take the limit Λ → ∞. As we will see, the thermodynamic limit of the quantity O [λ λ λ],Λ t at finite Λ > 0 involves integrals of the form The limit Λ → ∞ of these integrals for n > 0 only exists in a distribution sense, i.e. their integral with any smooth function of x, t has a well-defined limit when Λ → ∞. The resulting limits are denoted by I n (t, x) and have been worked out in [63] for n = 0, 1, 2 An equivalent representation is The process described above provides closed-form expressions at order O(c −m ) for the quantities O Finally, in order to obtain the out-of-equilibrium time evolution (36) this result needs to be evaluated for the saddle point root density ρ describing the quench protocol of interest.
4 Calculation of the one-point function σ 2 (0) t In this Section we apply the Quench Action approach combined with a 1/c expansion to compute the one-point function σ 2 (0) t .

The form factors
In order to evaluate the expression (36), one requires a closed-form expression for the form factors of σ 2 between energy eigenstates. In the case of interest, because of the structure of the non-vanishing overlaps with the initial state |Ψ BEC , the states entering (36) have a pair structure and will be denoteed |λ λ λ and |μ μ µ . Hence they have same (vanishing) momentum.

1/c expansion and particle-hole excitations
Employing a saddle-point argument in (37) shows that in the limit t → ∞ we have We use this and the pair structure of the states entering (37) to rewrite (37) as We now analyze this expression in terms of a 1/c-expansion [63]. In the c → ∞ limit G ± become the identity matrix and the ratio of overlaps takes a simple form Next we turn to the 1/c-expansion of the form factor. It is convenient to introduce some shorthand notations where The rapidities {λ i } and {µ j } are solutions to the Bethe equations (7) with Bethe numbers I j and J j respectively. The 1/c-expansion of the rapidity differences µ i − λ i is given by The 1/c-expansion of V ± j is computed by writing and then Taylor expanding the exponential and the logarithms. For a pair state this gives Combining (56) and (54) we obtain that the large c limit of the matrix U is given by To evaluate the determinant appearing in the form factor we use that for an invertible matrix A and two vectors u, v we have which implies that det Let us now introduce i.e. the number of Bethe numbers associated with the rapiditiesλ λ λ that are distinct from the Bethe numbers corresponding to the rapiditiesμ μ µ. Using (54) we find Putting everything together it follows that This establishes that the 1/c-expansion of the spectral sum (50) corresponds to an expansion in the number of particle-hole excitations. Since ν has to be even because of the pair structure of the states, the leading order term forμ μ µ =λ λ λ is obtained for ν = 2, i.e. two particle-hole excitations and is of order O(c −2 ). The next terms involve four particle-hole excitations and contribute only at order O(c −8 ). Since our goal is to compute the relaxation dynamics up to order c −4 , we can restrict our analysis to two particle-hole excitations.

Two particle-hole excitations
We now fix the rapidities λ λ λ > 0 of the representative state and denote its Bethe numbers by {I j }. We then consider µ µ µ > 0 such that the corresponding Bethe numbers J j are equal to I j except for J a = I a + n .
The usual exclusion principles in the Bethe ansatz impose that n = 0, J a > 0 and ∀i = 1, ..., N , J a = I i . The Bethe state |μ µ µ constructed in this way is a pair state that corresponds to a two particle-hole excitation over the representative state |λ λ λ . Taking into account only such states in the spectral sum (50) provides a 1/c-expansion up to and including O(c −4 ). Taking the difference between Bethe equations for the roots µ i and λ i and using the pair structure we obtain the following expansion for the positive Bethe roots with i = a Here we have introduced the convenient notation We next turn to the 1/c-expansion of the matrix U . We choose λ p = −λ a , so that the first term in U jk is O(c −1 ) except for j = a. This gives where We then employ the following identity obtained from (58) det j,k to obtain Here we have defined Using (56) we then obtain the following result for the 1/c-expansion of f (z) Noting that we finally arrive at the following expression for the determinant appearing in the form factor The expansion of the remaining terms in the form factor is more straightforward. We find Putting everything together we obtain The expansion of the ratio of the normalized overlaps is similarly straightforward Our final result for the 1/c-expansion of the summand in (36) is then This is a regular function of λ a and n and in the thermodynamic limit the sums over λ a and n can therefore be turned into integrals We refer the reader to Section 3.3 for the → 0 limit. Importantly (83) depends on the representative state only via the particle and hole densities. This shows that the typicality assumption underlying (37) indeed holds, at least to the order of the 1/c-expansion we are working in.
The usual typicality arguments suggest that in the thermodynamic this quantity will depend on the representative state only via its particle and hole densities. If this holds true then the generalized micro-canonical average in (84) can be dropped and We will see below that this is indeed the case due to rather delicate cancellations of contributions that depend on details of the representative state. The form factors entering (85) are given by [72][73][74][75][76][77] µ µ µ|σ (0) |λ λ λ λ λ λ |λ λ λ µ µ µ| µ µ µ where

Structure of the contributing "excited states"
In order to determine the order O(c −1 ) in our 1/c-expansion of (85) we need to know which "excitations" ν ν ν andμ μ µ will contribute to the spectral sums. Let us first remark that the limiting value taken by σ(x)σ(0) t when t → ∞ is obtained whenμ μ µ =λ λ λ in (85), as written in (49). In order to investigate the relaxation dynamics we will thus assume from now onμ μ µ =λ λ λ. We recall from [63] that the density form-factor for a one particle-hole excitation with rapidities µ µ µ above a state with rapidities λ λ λ takes the following form at order O(c −1 ) µ µ µ|σ (0) |λ λ λ λ λ λ |λ λ λ µ µ µ| µ µ µ The product of signs in this formula arises because we chose an Algebraic Bethe Ansatz description of the eigenstates (8) that is symmetric in the rapidities λ 1 , ..., λ N , in contrast to the Coordinate Bethe Ansatz description which is antisymmetric. In normalized form and for zero-momentum states, the two are related by a factor i<j sgn (λ i − λ j ) times a phase independent of λ i 's [78]. For a two particle-hole excitation where the Bethe numbers of µ µ µ are the same as the ones of λ λ λ except for I a , I b , we have [63] µ µ µ|σ (0) |λ λ λ λ λ λ |λ λ λ µ µ µ| µ µ µ Form factors with a higher number of particle-hole excitations are suppressed by at least a factor c −2 and we will ignore them in the following. We are now in a position to identify the dominant "excitations" contributing to the spectral representation at large c.
(i) "Type I" configurations contributing at O(c 0 ) and higher Because of the pair structure of bothλ λ λ andμ μ µ the leading order of the 1/c-expansion is obtained with states corresponding to a two particle-hole excitationμ μ µ aboveλ λ λ such that the Bethe numbers I a , −I a ofλ λ λ are replaced by J a , −J a inμ μ µ. We will assume this structure to be satisfied in the following.
Then the intermediate state ν ν ν that provides the leading O(c 0 ) contribution is obtained by imposing that it is a one particle-hole excitation above bothλ λ λ andμ μ µ. This implies that the Bethe numbers of ν ν ν have to be the same as those ofλ λ λ with the exception of I a or −I a , which is replaced by either J a or −J a . These contributions give the full result in the c → ∞ limit, which correspond to a quench directly from the BEC to the Tonks-Girardeau gas [44]. However, they also incorporate c −1 corrections due to subleading terms in the form factors.
(ii) "Type II" configurations contributing at O(c −1 ) At order O(c −1 ) contribution arise from other terms in the spectral sum as well. One class of terms corresponds to the case where ν ν ν is equal toλ λ λ (μ μ µ) and corresponds to a two particle-hole excitations aboveμ μ µ (λ λ λ). In this case one of the two form factors in (85) reduces to the expectation value of σ which equals the density D, while other form factor is of order O(c −1 ) since it involves states related by two particle-hole excitations. Closer inspection of (85) reveals that these contributions cancel Here we have used that since σ is a conserved quantity we have This leaves one remaining source for O(c −1 ) contributions, namely when the ν ν ν correspond to a one particle-hole excitation aboveλ λ λ (μ μ µ) and a two particle-hole excitation aboveμ μ µ (λ λ λ). As we will see below these terms give non-vanishing contributions to the spectral sum.
Lemma 1. Let f (λ, µ, ν) be a regular function that grows sufficiently slowly at infinity. Then in the thermodynamic limit we obtain Here we have defined Ω(λ λ λ) ≡ lim We stress that Ω(λ λ λ) is a quantity that in the thermodynamic limit depends on the choice of representative state not only through the root density ρ(λ). A proof of Lemma 1 is given in Appendix B.
Lemma 2. Let f (λ, µ, ν) be a regular function that grows sufficiently slowly at infinity. Then in the thermodynamic limit A proof of Lemma 1 is given in Appendix C.

First sum Σ 1 (x, t)
Writing out the various constraints in the summations explicitly we have We note that since the states have a pair structure, N is even and the Bethe numbers are half-odd integers, and so neither µ a or λ i + 2πn L can vanish in the denominators. The last two sums are two-dimensional sums with a prefactor 1/L 3 but only simple poles and hence vanish in the thermodynamic limit. The remaining two sums in (104) can be carried out using Lemma 2 and Lemma 1 respectively. This gives whereρ(λ) denotes the Hilbert transform of ρ(λ) defined bỹ

Second sum Σ 2 (x , t)
Writing out the constraints on the various summations explicitly we have The third sum is a two-dimensional sum with a prefactor 1/L 3 and no double poles and hence vanishes in the thermodynamic limit. The first two sums can be turned into principal value integrals, which gives

Third sum Σ 3 (x , t)
The third sum is a two-dimensional sum with a prefactor 1/L 3 , so can contribute in the thermodynamic limit only if there is a double pole. It follows that By writing out the constraint explicitly we have We see that the sum over λ a can be turned into an integral, while the remaining sums can be respectively carried out explicitly and expressed in terms of Ω(λ λ λ) (102). This gives

Result
Putting everything together, we obtain the following result for the contribution of two one particle-hole excitations to the spectral sum (85) We stress that C 1,1 λ λ λ (x, t) depends on the representative state λ λ λ not only through the root density ρ, but via the quantity Ω(λ λ λ) (102) as well.

Contributions arising from type II configurations
Let us denote by C 2,1 λ λ λ (x, t) the sum of contributions of type-II configurations to (85), i.e. configurations where ν ν ν corresponds to a one particle-hole excitation aboveλ λ λ (μ μ µ) and a two particle-hole excitation aboveμ μ µ (λ λ λ) respectively . There are altogether four cases: (i) The Bethe numbers of ν ν ν are those ofλ λ λ except for the replacement of I a or −I a by K a .
Denoting the corresponding root by ν a we have the following restrictions: ∀i, µ a = λ i ; ∀i, ν a = λ i ; ν a = µ a , −µ a .
(ii) The Bethe numbers of ν ν ν are those ofμ μ µ with only J a or −J a replaced by a K a . Denoting ν a the corresponding root, we have the restrictions ∀i, µ a = λ i and ∀i, ν a = λ i and ν a = µ a , −µ a . Cases (i) and (ii) are sketched in Figure 6.   Case (i) can be accounted for by always changing λ a for ν a , but allowing λ a to range between −∞ and ∞. One can also allow µ a to range between −∞ and ∞ by introducing a combinatorial factor 1 2 . In case (ii) the same holds true with λ a and µ a interchanged. Case (iii) can be accounted for by always changing λ b for µ a , but allowing both λ b and µ a to range between −∞ and ∞. One can also allow λ a to range between −∞ and ∞ by introducing a combinatorial factor 1 2 . In case (iv) the same holds true with λ a and µ a interchanged. In cases (i) and (ii) the product of all the signs appearing in (89) and (90) give a factor − sgn (λ a µ a ). In cases (iii) and (iv) they give a factor sgn (λ a µ a ). It follows that in these four cases we have Ψ BEC |μ μ µ λ λ λ|σ(0)|ν ν ν ν ν ν|σ(0)|μ μ µ Ψ BEC |λ λ λ μ μ µ|μ μ µ ν ν ν|ν ν ν e 2it(E(λ λ λ)−E(µ µ µ))+ix P (ν ν ν) In order to proceed it is convenient to decompose the rational functions in (113) using Using (113) and (114) we can express the sum of all type-II contributions to (85) in the form where In all four contributions Σ j (x , t) the respective first term only involves simple poles and therefore can be straightforwardly expressed in terms of principal value integrals in the thermodynamic limit. The other terms involve two simple poles and require a more elaborate treatment.

First term Σ 1 (x , t)
The contribution to Σ 1 (x , t) involving two simple poles is of the form where Resolving all the constraints, we have at leading order in 1/c The first two contributions can be computed by first summing over m, and then summing over n and λ i respectively, which involves one-dimensional sums with only a single simple pole. In the thermodynamic limit they can be readily turned into principal value integrals. The fifth and sixth terms are double sums with a factor L −3 and hence are completely dominated by their respective double poles. They yield The third term is of a very similar structure to Lemma 2 (103) and can be treated analogously. We write In the thermodynamic limit this becomes The two principal values can be brought under a single principal value as in (163), cf. Appendix A.2. Finally the fourth term in (119) can be calculated using Lemma 1 (101).
The first terms on the right-hand side can all be computed by performing successive onedimensional sums with only a single simple pole, which allows them to be turned into principal value integrals in the thermodynamic limit. The last term involves a two-dimensional sum with a factor L −3 and a summand featuring only simple poles. Hence it vanishes in the thermodynamic limit. We conclude that

Third term Σ 3 (x , t)
This contribution is straightforward to deal with. After writing the sum over µ a as the difference of a sum over vacancies and holes the sums over λ a,b can be factorized and will involve only single simple poles. It then follows that

Fourth term Σ 4 (x , t)
The contribution to Σ 4 (x , t) involving two simple poles is of the form where The first sum can be straightforwardly turned into a principal value integral and the second sum can be carried out using Lemma 1 (101). This gives

Result for all contributions arising from type II configurations
The combined contribution of all Σ n (x , t) can be brought into a simpler form by using that (i) the root distribution is even; (ii) at leading order in 1/c we can write and (iii) for x = 0 we have in a distribution sense This allows us to combine the contributions of the terms in Σ n (x , t) involving only a single simple pole into the following expression Our final result for the thermodynamic limit of all contributions to (85) arising from type II configurations is then where F 2 (λ, µ, ν; x ) is the function defined in (24).

Cancellation of the representative state dependence
Once both contributions (112) and (134) to the spectral sum are summed up, we observe that the dependence on the representative state through the quantity Ω(λ λ λ) exactly vanishes! This non-trivial cancellation suggests that the typicality assumption underlying (37) is indeed correct, even though the partial contributions do carry an additional dependence on the chosen representative state.
To arrive at the expression (22) written in the introduction, we sum up (112) and (134) and use that at leading order in 1/c 6 Calculation of the two-point function σ 2 (x, τ )σ 2 (0, 0) ∞ in the steady state We saw in (50) that the expectation value of an observable O t after the quench converges when t → ∞ to O ∞ given in (49). This limit value is thus expressed as an equilibrium expectation value of O in a representative state corresponding to the steady-state root density ρ that is fixed by the quench protocol. An interesting question is then how to characterize the physical properties of this steady state through its response functions. The dynamical correlation function of an observable O in an energy eigenstate |λ λ λ has a spectral representation in a basis of (unnormalized) energy eigenstates |µ µ µ of the form λ λ λ |λ λ λ µ µ µ| µ µ µ e iτ (E(λ λ λ)−E(µ µ µ))+ix(P (µ µ µ)−P (λ λ λ)) .

1/c expansion and particle-hole excitations
Let us again follow the same reasoning as in the previous sections, and investigate the leading behaviour of the form factor when c → ∞, for generic λ λ λ, µ µ µ satisfying the Bethe equations. The simple relation (137) allows us to directly use the results of [63] for the density correlations.

One particle-hole excitations
We now consider a one-particle-hole excitation above λ λ λ, namely a state µ µ µ such that all its Bethe numbers are those of λ λ λ except for I a which is replaced by I a + n. This results in constraints on the Bethe numbers n = 0 , ∀i = a, I a + n = I i .
We then can use the results of [63] because of the simple relation (137), which always holds since the momenta between the two states involved are necessarily different. We obtain | λ λ λ|σ 2 (0)|µ µ µ | 2 λ λ λ|λ λ λ µ µ µ|µ µ µ = 16 Interestingly, the a priori leading order O(c −2 ) contribution vanishes. As a result the one and two particle-hole excitations contribute at the same order in 1/c. Since (141) does not have poles the corresponding contribution to the spectral sum (136) is straightforward to compute and gives the first line of (27).

Two particle-hole excitations
Assuming that the momenta of the two states are different, i.e. that n = −m, one can again use (137) and [63] to obtain | λ λ λ|σ 2 (0)|µ µ µ | 2 λ λ λ|λ λ λ µ µ µ|µ µ µ = 16 This expression has no singularities and the corresponding contribution to the spectral sum is straightforwardly expressed as an integral in the thermodynamic limit. This gives the second line of (27). When n = −m, i.e. when the momenta of the two states are identical, we obtain from and that there are no singularities in n. As in this case there are only three sums we conclude that such contributions vanish in the thermodynamic limit.

Summary and Conclusions
In this work we have combined the Quench Action approach with our recently developed 1/c-expansion method for form factor sums in the Lieb-Liniger model to analyze a number of different observables after a quantum quench starting in the ground state of a non-interacting Bose gas. To the best of our knowledge our work is the first to obtain analytic results for quench dynamics in an interacting integrable theory beyond the asymptotic late-time regime.
Our work also uncovered a novel aspect regarding the application of typicality ideas to the analysis of quantum quenches in integrable models. We observed that carrying out partial summations of the spectral sums in the Quench Action approach can lead to results that violate the underlying typicality assumption and depend on details of the particular representative state selected. In the case at hand this dependence arises from the singular behaviour of overlaps at zero rapidity. But remarkably, we observe that this representative-state dependence cancels out between different types of particle-hole excitations at the order in 1/c of our calculation, yielding a significant check of typicality in an out-of-equilibrium setting. However, we are able to construct ad hoc initial states in a free theory for which these cancellations do not occur. This results in a failure of typicality, but this failure is weak in the sense that the problematic representative states are rare and can be avoided through a regularization procedure. A brief discussion of these findings is given in Appendix E.
Our work raises a number of interesting questions that should be investigated further. First, it is important to work out higher orders in the 1/c-expansion both for dynamical response functions and in the quench context. In particular, conjectured extensions of GHD predict that the two-point functions of σ 2 (x) will exhibit diffusive behaviour [70]. This is not seen in the leading order of the 1/c-expansion worked out here, but supposedly will appear at the next order. Second, it should be explored how to define truncations of the spectral sum that would be finite in the thermodynamic limit (not divergent and not exponentially small) for finite c. Indeed, the spectral sum truncation induced by the 1/c expansion generically exhibits terms polynomial in the system size that cross-cancel between different numbers of particle-hole excitations. Third, it would be very interesting to apply our strong coupling expansion method to dynamical correlations in other models like the Heisenberg XXZ chain [79][80][81][82]. These typically will involve bound states, and an important question is how to extend the strong coupling expansion in order to take their contributions into account. Fourth, it would be interesting to extend the analysis presented above to quantum quenches starting in inhomogeneous initial states [83,84]. Finally, we think it is important to arrive at a more complete understanding of the scope and limitations for applying typicality ideas to the calculation of dynamical correlations in and out of equilibrium.

A Principal value integrals
In this appendix we present details on principal value integrals used in the main text and the proofs of Lemma 1 and 2.

A.1 Double principal values
Given a function F (λ, µ, ν), we define its integral with successive double principal value as where the − symbols appearing in the right-hand side of this expression denote single principal values defined in (25). As shown in [63], the following relations hold and The integral with simultaneous double principal value is defined by As shown in [63], it is related to the integral with successive double principal value through the Poincaré-Bertrand-like formula A.2 Proof of equation (163) Using the identity (149) we obtain and In (163) the sum of these two quantities appears. The latter can be brought under a single simultaneous principal value because the excluded regions of the integral are identical (which is not the case of the successive principal values). Hence Using (146) we arrive at (163).

B Proof of Lemma 1 (101)
We start by adding the condition λ j = −λ i Finally we employ (149) to arrive at Eq (101).

C Proof of Lemma 2 (103)
We start by rewriting the multiple sum of interest as The first and second terms on the right-hand side can be turned into principal part integrals in the thermodynamic limit by first summing over n and then over λ i . The third sum, although two-dimensional with a prefactor 1/L 3 , is not negligible in the thermodynamic limit since it involves a double pole in λ i . Its thermodynamic in fact depends on the representative state λ λ λ through the quantity Ω(λ λ λ) defined in (102).
The two principal values can be brought under a single principal value according to the following relation, proved in Appendix A.2 This gives the desired result D Further results on σ 2 (x)σ 2 (0) In this appendix we collect a number of additional results on the two-point function after our interaction quench (22).

D.1 Alternative expression for (22)
In this section we present an alternative expression for the two-point function after the quench (22), that is particularly useful for numerical purposes. It is based on the observation that the first terms in the 1/c-expansion of the steady state root density (15) take the simple form which allows one to carry out some of the integrals in (22). To that end we introduce We then find and H t (λ) = − e −2itµ 2 λ − µ dµ .
The integrals in (168) can be deduced by differentiating this with respect to x. A straightforward calculation then shows that at t = 0 we indeed recover the two-point function in the BEC state at order O(c −2 ) D.3 Consistency check II: x → 0 limit Since for x = 0 we have σ(x)σ(0) = ψ † (x)ψ † (0)ψ(x)ψ(0) , the correlation function σ(x)σ(0) t should approach σ 2 (0) t = O(c −2 ) in the x → 0 limit. Using (171) at x = 0 and simplifying (22) by exploiting that the root density ρ(λ) is even we find after some calculations that indeed D.4 Some remarks on the limit x, t → 0 As we have noted in the main text the limit t → 0 of our result for σ 2 (0) t does not recover the correct result for the expectation value of σ 2 in the BEC initial state, D 2 . On the other hand, we have just shown that the limit x → 0 of σ(x)σ(0) t=0 does reduce to D 2 . On a technical level it can be traced back to properties of the integral which vanishes if one first takes the limit x → 0 and then t → 0, but gives a finite result if one takes first t → 0 and then x → 0.

E Typicality and Quench Action method
In this Appendix we present an ad hoc initial state in a free theory for which the Quench Action spectral sum for the out-of-equilibrium dynamics is representative state dependent. We consider a simple tight-binding Hamiltonian on a ring where a † j , a j are fermionic creation and annihilation operators satisfying canonical anticommutation relations {a j , a † k } = δ j,k . The Hamiltonian is straightforwardly diagonalized by a canonical transformation to Bogoliubov fermions in momentum space where k n = 2πn L and {b p , b † k } = δ p,k . We denote the Bogoliubov vacuum state by |0 . We now consider a quantum quench where the system is initialized in a Gaussian state parametrized by a fixed arbitrary function K(p) For our purposes it is sufficient to focus on the Green's function G(n, t) = I(t)|a n+1 a 1 |I(t) .
Let us now try to recover this with the Quench Action approach. The normalized overlaps of the initial state with an eigenstate |λ λ λ = k∈λ λ λ b † −k b † k |0 are λ λ λ|I = k∈λ λ λ iK(k) L/2−1 n=1 from which one finds the root density characterizing the non-equilibrium steady state reached at late times after the quench ρ(k) = 1 2π .
So far we have closely followed the discussion in [23]. However, let us now consider the following singular behaviour with m ≥ 1 an integer, and define a representative state |λ λ λ by replacing k 0 ∈ λ λ λ by k 0 . By construction |λ λ λ is a micro-state that for any choice of k 0 , k 0 corresponds to the macro-state with particle density ρ in the thermodynamic limit, and in particular the extensive parts of all local conservation laws are the same for |λ λ λ and |λ λ λ . Let us choose k 0 finite in the thermodynamic limit, and k 0 = O(L −1 ). We observe that λ λ λ|(a n+1 a 1 )(t)|I λ λ λ|I = λ λ λ |(a n+1 a 1 )(t)|I λ λ λ |I + i L K(k 0 )e 8it sin 2 (k 0 /2) e ik 0 n − K(k 0 )e 8it sin 2 (k 0 /2) e ik 0 n .
This shows that the two choices of representative state lead to different results in the thermodynamic limit, which generally does not even exist as K(k 0 ) ∝ L m . This shows that for this particular initial state a naive application of typicality ideas fails. However, a few comments are in order. First, since ρ h (k) ∼ k 2m at small k, the smallest hole in a representative state λ λ λ is typically of order L −1/(2m+1) , and in this case the additional terms are negligible indeed. Hence for such "typical" states, typicality ideas can be applied. This fact is confirmed numerically by observing that when one averages (185) over representative states, one indeed recovers (181). Second, in the problem at hand one can slightly change the initial state by imposing for example K(k) = K(δ) for k < δ for a fixed small δ. With this "regularisation" one obtains (187), which is now well-behaved and allows for the limit δ → 0 to be taken. In this limit one recovers the expected result (181).