Emergence of curved light-cones in a class of inhomogeneous Luttinger liquids

The light-cone spreading of entanglement and correlation is a fundamental and ubiquitous feature of homogeneous extended quantum systems. Here we point out that a class of inhomogenous Luttinger liquids (those with a uniform Luttinger parameter $K$) at low energy display the universal phenomenon of curved light cones: gapless excitations propagate along the geodesics of the metric $ds^2=dx^2+v(x)^2 d\tau^2$, with $v(x)$ being the calculable spatial dependent velocity induced by the inhomogeneity. We confirm our findings with explicit analytic and numerical calculations both in- and out-of-equilibrium for a Tonks-Girardeau gas in a harmonic potential and in lattice systems with artificially tuned hamiltonian density.


Introduction
Light-cone propagation of signals in condensed matter and cold atom physics has attracted a lot of attention in the past decade. Contrarily to what happens in relativistic high-energy physics, where no signal can propagate faster than the speed of light, in non-relativistic quantum systems a maximal speed of propagation does not necessarily exist. A very profound result of mathematical physics is that excitations in discrete models with a finite-dimensional local Hilbert space and with local interactions have indeed a bounded velocity v max [1]. This classic result by Lieb and Robinson has had many ramifications over the years [2]. Roughly, it implies that when these systems are perturbed at a given position x 0 , after a time t the signal generated by the perturbation cannot propagate farther than x 0 ± v max t. The existence of a maximum speed of propagation has deep consequences in the non-equilibrium dynamics of these systems, resulting in the so called light-cone spreading of correlations [3][4][5][6][7] and entanglement [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]28] e.g. following a global or a local quantum quench, as reviewed e.g. in [34,35]. Such light-cone spreading has been also experimentally observed [36][37][38]. Furthermore, the presence of a velocity bound also affects the transport properties of these systems [39][40][41]. Apart from the Lieb-Robinson rigorous mathematical bound, condensed matter systems sometimes display the same behavior as relativistic models of high-energy physics: this happens at quantum critical points with dynamic exponent z = 1, where emergent massless quasi-particles propagate at a critical (sound) velocity v which plays the same role as the speed of light in high-energy physics. This relativistic behavior emerges only when probing the low-energy (large distance) universal features of these systems. In particular in one spatial dimension, these systems mainly fall in the Luttinger liquid paradigm, as pointed out long ago by Haldane [42]. A wide class of one-dimensional systems such as spin chains, onedimensional superfluids, electrons in one-dimensional wires, or clouds of ultra-cold atoms in tightly elongated traps are indeed well described by Luttinger liquid theory.
In the limit of large distances and time separations, all Luttinger liquids in equilibrium display power law singularities in time-dependent correlators (such as the Green function) on the light cone. This is a straightforward manifestation of the light-cone propagation of signal at low energy. However, this is more cleanly devised in non-equilibrium situations such as a local quench following a cut-and-glue protocol [21][22][23]28]. Quite remarkably, since a local quench only probes the low-energy part of the many-body spectrum, the light-cone spreading of entanglement and correlations can be observed in general Luttinger liquids, even in those systems that have unbounded velocities and do not satisfy Lieb-Robinson [28].
The purpose of this paper is to investigate light-cone propagation in certain inhomogeneous quantum critical systems in 1d. We focus on two particular microscopic realizations of an inhomogeneous Luttinger liquid: impenetrable bosons in a non-uniform trap potential, and (Right) in an inhomogeneous system, the velocity v of gapless excitations becomes a function of position, and, as a result, the propagation of gapless excitation is along curved light-cones. The goal of this paper is to confirm the emergence of such curved light cones in some concrete microscopic models.
spin chains with position-dependent couplings between neighboring spins. These two systems share a crucial property, on which our entire analysis is relying: the Luttinger parameter K is constant in those systems. This will be explained in more details shortly; it is important, however, to emphasize that there exist other systems, also in the Luttinger liquid universality class, which would not have a constant Luttinger parameter K in inhomogeneous situations, to which our analysis does not apply (the most remarkable example being the 1d delta Bose gas in a non-uniform trap potential, away from the Tonks-Girardeau limit, see e.g. Ref. [29]). Our calculations are based on the recently introduced framework to describe gapless inhomogeneous systems by means of conformal field theories (CFT) in a curved background [43], see also [44][45][46]. Our main conclusion will be that, in those systems that are inhomogeneous, but have a constant Luttinger parameter K, the velocity of gapless excitations v, analogous to the speed of light, is position-dependent, v → v(x). This has a dramatic consequence on light-cone propagation: the light-cones become curved, as illustrated in Fig. 1. Our approach allows to quantitatively characterize the curvature of the light-cone, i.e. to determine analytically the explicit space dependence of the sound velocity. Similar curved light-cones have also been observed in time-dependent quenches [30][31][32][33], where the velocity depends on time but not position. In the present work we focus on space-dependent velocities instead.
The paper is organized as follows. In Section 2 we show how the inhomogeneous Luttinger liquid naturally appears when describing a gapless 1d quantum system and we introduce a few models with constant Luttinger parameter K, as well as all other assumptions exploited in our calculation. In Section 3 we apply the framework of CFT in curved background to the Tonks-Girardeau gas in a harmonic trap both in equilibrium and for a local quench. In Section 4 we present our results for lattice models with inhomogeneous hamiltonian density. Finally in Section 5 we draw our conclusions and discuss some open problems.

Inhomogeneous Luttinger liquid with constant Luttinger parameter K
The Luttinger liquid is an effective description of systems of particle or spins that is valid on distance-and time-scales that are large compared to microscopic ones (for instance, the lattice spacing, or the interparticle distance). Here we quickly introduce the formalism that we will apply later to some specific microscopic models. Only the aspects of Luttinger liquid physics that are strictly necessary to this paper will be discussed; moreover, our treatment is slightly non-standard, because our goal is to emphasize the role of the inhomogeneity of the model. For a more standard introduction and thorough discussion of Luttinger liquid physics, we refer the reader to the textbooks [47] and [48].

The height field
Imagine that we have a microscopic model of mobile particles on the real line, at positions x 1 , x 2 , . . . , x N , and that we are interested in describing correlations of the density operator ρ(x) ≡ j δ(x − x j ). A configuration of particles can be represented by the height function h(x), defined such that The factor 2π is introduced for notational convenience. Notice that this relation defines h(x) only up to a constant shift h(x) → h(x) + const. If we know how to calculate correlation functions of the height field, then we know how to calculate correlation functions asymptotically in the microscopic model. This is the basic idea of both Luttinger liquid theory and the Coulomb gas approach to 2d statistical mechanics [49]. The functions ρ(x) and h(x) represent the configuration of the system at a given instant, say t = 0. The time dependence is easily written in terms of the Hamiltonian H of the microscopic model or in imaginary time For instance, the two-point function ψ 0 | h(x, τ )h(x , τ ) |ψ 0 , in the ground state |ψ 0 of H, is defined as the (imaginary-) time-ordered expression We have a height field h(x, τ ) that lives in 2d spacetime; the goal is now to write an action S[h] for this field, that can be used to calculate correlation functions.
2.2 Coarse-grained action for the height field, and emergence of the gaussian free field We now assume that the Hamiltonian H of the microscopic model is local, and we adopt a coarse-grained perspective. Because the model is local and invariant under constant shifts h(x, τ ) → h(x, τ ) + const, the fluctuations of h must be captured by a local action of the form One assumes for simplicity that there is a unique classical configuration h cl. that minimizes the action (5), up to constant shifts of h. Then, substituting h → h + h cl. , one can expand around h cl. , which gives + higher order terms. (6) All higher order terms have scaling dimension larger than 2, and are irrelevant in two dimensions. This is the fundamental reason for the universality of the Luttinger liquid model: under very general assumptions, the underlying field theory turns out to be a free, gaussian, field theory, with the quadratic action (6). There are only three adjustable parameters in this action, namely the three independent entries of the hessian ∇ 2 L at the point h = h cl. . It is convenient to interpret the inverse of the hessian as a metric tensor Then the action (6) becomes, with (x 1 , x 2 ) = (x, τ ), where 1 K ≡ 4π det(∇ 2 L). In principle, the value of K may depend on position x, however in this paper, we will assume that it does not. Below, we will discuss in more details in which type of microscopic models this assumption is valid. The case of non-constant K is also very interesting (and it has appeared, in Hamiltonian form, for some particular choices of the function K(x) in Refs. [29,50,51], but in general it is more complicated, and it will be studied elsewhere).
The key point, which allows to solve the Luttinger liquid model exactly, is that the action (7) is nothing but a gaussian free field in the conformal class of the metric [g]; indeed the action (7) is invariant under Weyl transformations g → e 2Φ g, and thus depends on the conformal class of g, noted [g], rather than g itself. Solving the model then boils down to finding its propagator h(x)h(x ) , which is an easy exercise in conformal mapping and potential theory. We will come back to this below, but first, let us introduce and discuss the two additional assumptions that we make throughout this paper.

Our assumptions
Throughout this paper, we will be dealing with microscopic systems in the Luttinger liquid universality class, that satisfy the following three assumptions.
Assumptions in this paper: 1. the parameter K is uniform, i.e. it does not depend on position (x, τ ) in spacetime; consequently, large-scale correlations can be obtained from the action (7); 2. we consider small perturbations around the ground state of H; consequently, the entries of the metric tensor g do not depend on time τ (but they do depend on position x); 3. the system is time-reversal invariant; consequently the component g xτ of the metric tensor (which transforms as g xτ → −g xτ under time-reversal symmetry) is zero.
The first one, namely the fact that K is uniform, has already been mentioned, and is further discussed in the next section. For now, let us focus on the consequences of assumptions 2 and 3. Because of Weyl invariance of (7) only the two ratios g xτ /g xx and g τ τ /g xx play a role in the model, Then, as a consequence of assumptions 2 and 3, we have where we defined v 2 (x) ≡ g τ τ (x)/g xx (x). Eq. (9) is very simple, yet it is the most profound equation of the paper. We will see that gapless excitations automatically propagate along the geodesics of this metric, which lead to the appearance of the curved light-cones illustrated in Fig. 1. It is worth emphasizing that our assumptions 2 and 3 are merely technical ones; while they simplify our task considerably in the following, they are in principle not crucial. One could easily relax them without affecting our main conclusion, which is the emergence of curved light-cones. One particular example where the entries of the metric g ij are all non-zero and depend both on position and time is the one of the quench from a domain-wall initial state, which was treated in Refs. [43,44]. Throughout the paper, it will be useful to use the coordinatẽ where x L is some reference point, that we will take to be the position of the left boundary of the system. Physically,x is the time that is needed by a signal that propagates at the local velocity v(x) to travel from the boundary at position x L to position x. The metric (9) then takes the form where e σ = v(x) and (z,z) = (x + iτ,x − iτ ). In these coordinates, the action (7) takes the standard form

Discussion: why K is uniform for a certain subset of Luttinger liquids
We have explained that we restrict our analysis to the case of systems with constant Luttinger parameter K; let us now emphasize that, even though this is a strong restriction, there are still various interesting microscopic models that satisfy this assumption.
Models that map to free fermions: K = 1. The first important class of microscopic models with constant K is the ones that map to free fermions. Relevant examples include models of hard-core boson [52] or hard-core anyons [53], either in the continuum or on the lattice and spin-chains such as the XX model. Let us briefly recall why, in these models, one necessarily has K = 1, ruling out the possibility of a varying Luttinger parameter K.
Since the model has underlying non-interacting fermionic degrees of freedom, the theory described by the action (12) must contain an operator that creates a fermionic excitation, that evolves freely. Since we are in an effective theory with no mass, the free fermion excitation must be right-or left-propagating, with velocity ±v (or ±v(x) in the inhomogeneous case with position-dependent velocity). To construct an operator that creates left-and right-moving excitations, one proceeds as follows. The propagator derived from the action (12) satisfies the Equivalently, one can think of h(z,z) as being the sum of chiral and anti-chiral components which have the following propagators for consistency of (14) with (13). The operators that create chiral excitations are then of the form e iαϕ(z) and where the exponential is normal ordered. These operators create one particle, so they must have charge one, by definition. This means that the coefficient of their operator product expansion with the density operator ρ(z,z) = must be i 2π and −i 2π respectively. This implies Similarly, by looking at the operator product expansion of e iαϕ(z) and e −iαϕ(z) with the stresstensor, one gets their spin, which is equal to But these operators are fermions, so they must have spin ± 1 2 . Thus we arrive at K = 1 (for free fermions) (20) as claimed.
Models with artificially varying hamiltonian density. Models with underlying free fermionic degrees of freedom are not the only ones that possess a fixed parameter K. If a Luttinger liquid is SU(2) symmetric, then it necessarily has K = 1/2, for reasons that are similar to the ones we just explained for the free fermion case. Namely, SU(2) symmetry implies that the massless left-and right-excitations can be separately organized in SU (2) multiplets. For this to be possible, one must find chiral operators that act as the generators S ± , S z of the SU(2) algebra, up to central extension. These three current operators are of the form e ±i ϕ √ K , and they must have scaling dimension 1 (or, equivalently, spin 1), and this fixes K = 1/2. Thus, we could think of investigating inhomogeneous microscopic models with SU(2) symmetry and emergent Luttinger liquid physics. It turns out that such microscopic models are very constrained. To see this, let us start from the homogeneous antiferromagnetic Heisenberg spin-1/2 chain with S = (S x , S y , S z ). It is well-known that the large-distance, low-energy, behavior of this model is captured by the Luttinger liquid model with K = 1/2. Now we ask: how can we make this model inhomogeneous, while preserving SU(2) symmetry, such that the large-scale behavior will be captured by an inhomogeneous Luttinger liquid with constant K = 1/2?
The simplest answer is obviously: vary the coupling J with position, J → J(x). If J j varies slowly enough with position x, then we will indeed have a model that is locally like a homogeneous Luttinger liquid, but since the velocity v is proportional to the coupling J, it will naturally vary with position. Another possibility would be to include position-dependent next-nearest neighbor terms S x · S x+2 . Such terms would also affect the local velocity; in addition, they could also introduce perturbations by operators that are marginally relevant or irrelevant, depending on the sign of the coupling (in the marginally relevant case, this would lead to dimerization, so it would destroy the long-range Luttinger liquid physics in which we are interested).
Thus, we arrive at the conclusion that, imposing SU(2) symmetry in order to protect the value of the Luttinger parameter K, the simplest possible inhomogeneous deformation of the Heisenberg Hamiltonian is to vary the coupling J with position.
A similar conclusion holds for supersymmetric models 1 . The superysmmetric spin chain of Refs. [54][55][56] falls in the Luttinger liquid class with K = 1/3. Lattice supersymmetry is reflected in the continuous field theory in the fact that, at this particular value of K, there are two chiral operators e ±i √ 3ϕ (and anti-chiral ones e ±i √ 3ϕ ) with conformal dimension 3/2, that generate an N = (2, 2) supersymmetry algebra; see Ref. [54] for a detailed study of this model. If one tries to find an inhomogeneous deformation of this model while keeping supersymmetry, such that the value of K remains protected, then one finds that the simplest deformation allowed is also just to vary the amplitude of the Hamiltonian density with position.
Hence, we are naturally lead to the following class of inhomogeneous models, where we start from a homogeneous hamiltonian, and simply vary its hamiltonian density with position. Yet, such deformations can be made for any hamiltonian, and we do not need to refer to either SU(2) symmetry or supersymmetry, or any other symmetry, any longer (although, as we just discussed, such symmetries may be a strong motivation to look into this specific class of Hamiltonians). Making the coupling strength position-dependent is a very simple modification that leads to many other examples of inhomogeneous Luttinger liquids with constant parameter K. Namely, given a homogeneous Hamiltonian H with hamiltonian density h x , which is known to be described at large distance by the (homogeneous) Luttinger liquid model with velocity v and parameter K, we can artificially vary the hamiltonian density for some slowly-varying function f (x) > 0, and this will give an inhomogeneous Luttinger liquid with position-dependent velocity v(x) = f (x)v. We will study this class of models in more detail in section 4. The rainbow chain [57][58][59] is a prototypical example of these models with artificially varying Hamiltonian density, which indeed has been already studied with the curved spacetime approach [45]. Models with modulated hamiltonian density have also appeared in Refs. [68][69][70][71][72].

A concrete microscopic example: the Tonks-Girardeau gas in a harmonic trap
For simplicity, we start with a very well-known exactly solvable model: impenetrable bosons in a harmonic trap, defined by the Hamiltonian with c → +∞, and we calculate the boson propagator Ψ † (x, t)Ψ(x , t ) , at different times. We compare the prediction from the effective Luttinger liquid description, valid at large scales, to the exact solution.

Prediction from inhomogeneous Luttinger liquid approach
In the ground state, the average density is given by Wigner's semi-circle law for a large number of particles N 1. The position-dependent velocity of low-energy excita- and the metric of the effective field theory that describes those low-energy excitations is ds 2 = dx 2 + v(x) 2 dτ 2 . As explained in section 2.3, it is more convenient to work with the stretched spatial coordinatex which represents the time needed by a low-energy excitation emitted at the left boundary of the system to arrive at point x. The time needed to travel from the left to the right boundary isL = du v(u) = π. The one-particle density matrix in imaginary time was calculated recently by Y. Brun and one of us in Ref. [46]. The expression at different times was not the primary focus of that reference, but it appeared as an intermediate step in the derivation of the main result, see Eq. (35) in [46]: where w = e iπ(x+iτ )/L and |C| 2 is a numerical constant whose exact value was obtained in Refs. [60,61]: |C| 2 = G 4 (3/2)/ √ 2π 0.521414. After a few simplifications, and injecting the explicit formula forx for the harmonic trap, this becomes which is a non-trivial generalization, for τ = τ , of the equal time correlation in a harmonic trap [62,63]. Now let us perform the Wick rotation τ → it. Without loss of generality, we fix t = 0, and write x 0 instead of x . This gives where is one of the two positions reached by the signal emitted from (x 0 , 0) at time t -one position is reached by the signal propagating to the right, the other by the signal propagating to the left. Straightforwardly, an expression for the different time correlation of hard-core anyons can be obtained generalizing the recent result for equal times [64].
Color density plots as a function of both position and time are shown in Fig. 2. The microscopic and CFT results are shown side by side, and display the same curved lightcones, given by (31). The calculation in the microscopic model is explained in the next section.

Numerical check
The Tonks-Girardeau gas can be mapped onto a system of free fermions. The evaluation of correlations functions in such systems boils down to the computation of determinants [65][66][67] (or Pfaffians if there are pairing terms, which is not the case here). Such a general observation is a consequence of Wick's theorem [66], which applies for any correlation in a fermionic state, even at finite temperature.
Here we use this to derive exact determinant formulas for the time-dependent one-particle density matrix, which we express as a (N + 1) × (N + 1) determinant for a system of N particles. The only slight difficulty lies in the fact that the operator that creates a particle at position x, Ψ † (x) is bosonic. The problem can be circumvented by using the second quantized formalism, and attaching a Jordan-Wigner string to a fermionic operator, as is routinely done to diagonalize spin chains such as the XX or Ising chain. More precisely, we write where ψ † (x) is fermionic operator, that satisfies the canonical anticommutation relations In such a language the ground state wave function simply reads 2 where |0 is the fermion vacuum annihilated by any ψ(x). The operators satisfy the commutation relations {ψ k , ψ † l } = R u k (x)u l (x) dx = δ kl . The u k are eigenfunctions of the single-particle Schrödinger operator − 1 2 ∂ 2 x +V (x). When the potential is harmonic the single particle states are known explicitly where H p (x) is the p-th Hermite polynomial. The associated single particle energies are determined from [H, ψ † k ] = (k)ψ † k , with (k) = k − 1/2. Let us now write down the timedependent one-particle density matrix as a fermionic correlator. Using operators in Heisenberg representation Ψ † (x, t) = e iHt Ψ † (x)e −iHt and (33) the density matrix reads
x t 2π to rewrite the above correlator as Now we finally make use of Wick's theorem: where the averages are taken in the fermion vacuum. The matrix elements are given by The determinant formula (40,41) generalizes the determinant results of Refs. [60,62] to unequal times t = t . We also note that the formula holds for arbitrary potential, provided the u k (x) solve the corresponding single-particle Schrödinger equation.
While an asymptotic evaluation of this formula in the limit N → ∞ is tedious and would probably require advanced techniques (see e.g. [73] for the uniform case), numerical evaluation for large but finite N is straightforward. In practice the time-consuming step is to evaluate (44), which we rewrite as with The α kp (x) may be evaluated to arbitrary precision using a gaussian quadrature method, and the infinite sum in (45) truncated to some value p max much larger than the Fermi level p = N . Two comparisons between the CFT and exact lattice formula for a finite number of particles are shown in Figs. 2 and 3. In Fig. 3, we fix the two position x, x to some values and look at the evolution with time. The agreement is very good, and improves further when increasing the number of particles.

Cut-and-glue protocol at the center of the harmonic trap: light-cones and evolution of entanglement entropy
Another way of providing compelling evidence for light-cone propagation of signals is to calculate the evolution of the entanglement entropy after a local quench. Here we focus on the "cut-and-glue" protocol [21], which consists in preparing the system in the ground state of the Hamiltonian where we imagine that V cut (x) is a Dirac delta V cut = V 0 δ(x) at the center of the trap, with very large amplitude V 0 → +∞, such that no tunneling is allowed between the left and right half-systems. Then, at t > 0, the infinite barrier V 0 δ(x) is switched off, and we let the system evolve with the original hamiltonian H. This "cut-and-glue" protocol is standard in the literature, and is viewed as a small, localized, perturbation of the system at position x 0 = 0. It has been largely studied with CFT methods in the homogeneous situation [21,23]. Following the notations of the previous section, we anticipate that the positions will play an important role. The system is inhomogeneous, with a spacetime metric ds 2 = dx 2 + v(x) 2 dτ 2 . In stretched coordinatesx (right), the metric is ds 2 = e 2σ [dx 2 + dτ 2 ] with e σ(x) = v(x), so it can be made flat with the Weyl transformation g µν → e −2σ g µν . This is the trick on which we rely to adapt the calculation of the one-point function T n (x, τ ) from Ref. [23] to this inhomogeneous problem.
While it is possible to calculate the expectation values of local operators after this local quench [21,74], here we focus on the time evolution of the entanglement entropy. The calculation of the entanglement entropy follows from Ref. [23]-see Eqs. (34)- (39) in that reference-. We need to evaluate the one-point function of a certain primary operator (called twist operator [24][25][26][27], coming from the implementation of the replica trick in CFT [75]), noted T n (x, τ ), in the euclidean space-time geometry shown in Fig. 3.3, in the metric ds 2 = dx 2 + v(x) 2 dτ 2 = e 2σ [dx 2 + dτ 2 ] with e σ(x) = v(x). The first step is to relate this one-point function in curved space to the one in flat space, with a Weyl transformation [43] where ∆ n = c 12 (n − 1/n) is the scaling dimension of the twist field T n [75]. Then we use a result of Ref. [23]: the formula for T n (x, t) flat , obtained after performing the Wick rotation τ → it, is given by Eq. (38) in that reference. In our notations, it reads Here , which has the dimension of a time, comes from the definition of the spacetime geometry shown in Fig. 3.3, and it plays the role of a UV cutoff. By dimensional analysis, we expect it to scale as ∼ 1 ρ v in the center of the trap, which gives ∼ 1/N . When we analyze the behavior of that expression when → 0, we find that we need to distinguish two cases: Now let us relate this to the entanglement entropy. In CFT, the Renyi entropy is given by , where a 0 is some UV length scale. In the Tonks-Girardeau gas, the only microscopic length scale is the inverse particle density 1/ ρ(x) , so we must have a 0 ∝ 1/ ρ(x) up to some numerical constant. It is important to emphasize that, contrary to the homogeneous case, the UV length scale a 0 is actually position-dependent, as noted in Ref. [43]. The Renyi entropy is thus with T n (x, τ ) flat given by Eq. (51). Finally, by specifyingx = arccos (53) where x ± (t) was given in Eq. (48), and where the additive numerical constants are independent of N , the number of particles in the trap. This clearly illustrates the light-cone effect: the entropy stays constant while x is outside the interval [x − (t), x + (t)], and it starts growing as soon as the signal reaches the position x. Again, what is interesting here, compared to the well-known homogeneous case, is that the light-cone is curved: the propagation speed is not constant, but it varies with position. We now turn to numerical checks of formula (53).
Numerical checks. We compare here the CFT result to the microscopic model, where the entanglement entropy may be computed using a variation of the correlation matrix techniques [76][77][78]. We are dealing here with particle conserving Hamiltonians, and states that are obtained by filling some Fermi sea, |Ω = N k=1 χ † k |0 , where N is the number of particles and |0 is the Fermion vacuum. For such states the von Neumann entropy of some subsystem A is given by the explicit formula where C A = (C kl ) 1≤k,l≤N , is the overlap matrix with elements [77,78] The χ A † k are obtained by real-space projection of χ † k onto subsystem A: This result may be used to compute the entanglement entropy after the cut and glue quench, in the harmonic potential. Denote as before by ψ † k = dx u k (x)ψ † (x) the fermionic modes that diagonalize the final hamiltonian, [H, ψ † k ] = (k − 1/2)ψ † k , and take an initial state of the form Thus, the problem of computing the entanglement entropy is reduced to computing the N ×N overlap matrix C A (t), with matrix elements and then evaluating (54), which can be done by simple diagonalization. In the case of the Harmonic potential the single particle wave functions before and after the quench may be expressed in terms of Hermite polynomials, which leads after some algebra to The kernel K is known as the Mehler kernel, and may be written as which means the overlap matrix may be efficiently computed by performing gaussian quadrature on (59) and (60). The results are shown in Fig. 5. It is evident that the agreement between numerics and CFT prediction is excellent, but there are also some deviations which are easily understood in terms of supersonic propagation of some particles as already similarly pointed out in [28]. Indeed in the Tonks-Giradeau gas, there is no bound on the speed of propagation of information and so, even if subdominant, there are quasi-particles traveling with a velocity larger than the Luttinger liquid one v(x). This is evident from the fact that the entanglement entropy slowly and slightly deviates from the initial value before than predicted by CFT. Anyhow, this is clearly a small subdominant effect and the light-cone dynamics is clearly visible in Fig. 5.

Models with inhomogeneous hamiltonian density
Let us consider a critical spin chain Hamiltonian with open boundaries, of the form whereĥ x is the Hamiltonian density (an operator acting on a finite number of neighboring sites). We denote by v f the speed of propagation of low-energy excitations. For free fermionic system this coincides with the Fermi speed; for more complicated models this needs to be determined a priori. It is always possible to consider the following deformation where f (u) is some smooth non-negative function of u. In particular, the so-called sinesquare deformation f (u) = sin 2 π(u+1) 2 has been the subject of several investigations in the past few years [68][69][70][71][72], as the ground state in that case mimicks that of a chain with periodic boundary conditions. Here we look at generic examples of such deformations, in a sense which we explain in the following subsection. For finite but large these automatically give models with spatially-dependent velocity v(x) = f (x/ )v f , that also varies slowly at the lattice scale, and where K is the same as in the original homogeneous chain.
The purpose of this section is to investigate the curved light cone physics that emerges from a non-uniform velocity. Our main example will be an inhomogeneous XXZ spin chain, with Hamiltonian densityĥ and |∆| ≤ 1 (gapless regime).

Dynamical correlations
Perhaps the simplest illustration of our ideas is provided by the following dynamical spin-spin correlation function where S 3 x e −iH[f ]t , and the average is taken in the ground state of the inhomogeneous chain. In a homogeneous setting such correlations are well described by CFT, and the lightcones given by the simple equation The red thick curves are the CFT predictions (68) for the lightcone. As can be seen they match very well the numerical data.
In the XXZ spin chain the speed v f is known exactly from the Bethe-ansatz solution [73] v f = π sin γ 2γ , and the lightcone is simply given byx which is the main result of this section. We have assumed here for simplicity that the stretched coordinate (67) is always well defined. This means we need to exclude deformations that vanish too fast -typically as (x − x c ) or faster -for a certain discrete set of points x c . In that case, the integral defining (67) could become divergent if one of the critical points x c lies in the interval [x 0 , x]. Physically this simply means that any signal would take an infinite time to go from x 0 to x. We find it convenient to avoid such situations; in particular, this rules out the sine-square deformation of Refs. [68][69][70][71].
A numerical check of formula (68) can be performed as follows. First, find the ground state |ψ of H[f ] using DMRG. Then, apply S 3 x to the ground state, and perform the unitary time evolution. Since applying S 3 x amounts to a local perturbation, the entropy generated after   Figure 7: Entanglement entropy after the cut-and-glue quench, and comparison with the CFT formula (71). As in Fig. 6 we chose f (u) = 1/(1 + 3u 2 ), which implies˜ = 2 , to make the chain inhomogeneous. Left: Cut in the center, x 0 = 0, with excellent agreement. Right: Cut at x 0 = /2, with very good agreement, but slight deviations. Note that the entropy stays constant until time t =x 0 = 5 /8, yet another illustration of the curved light-cones.
such a "local quench" only grows logarithmically with time, which means relatively large times may be accessed. Then, the overlap of the time evolved state |ψ(t) = e iHt S 3 x |ψ with S 3 x 0 |ψ may be easily computed. In practice, we used the ITensor C++ library [79] and were able to access times of the order of t 100 for system sizes 2 = 128, which was sufficient to achieve the desired separation of scales, and clearly see the lightcones. The results are shown in Fig. 6 for two different values of ∆ in the gapless regime. As it can be seen the agreement is excellent.

Entanglement entropy after a local quench
Another illustration is provided by the entanglement entropy following the cut and glue quench. Similar to what was done in section. 3.3, we consider the ground state of the Hamiltonian The entropy in the flat (or homogeneous) setting has been derived in [23], which yields S(x 0 , t, ) = c 6 log sin 2 πv f t 2˜ − sin 2 πx 0 2˜ + cst.
A numerical check of this formula is shown in Fig. 7 for the inhomogeneous XX chain (∆ = 0), where very large systems and late times may be achieved, and v f = 1. The agreement between the inhomogeneous CFT predictions and the numerics is extremely good and a clear quantitative evidence for the curved light-cone is shown. We must anyhow mention that at a more careful look, very tiny deviations between the two curves can be observed when the measure is not done in the middle (x 0 = 0, see right panel), similar to what we observed in section 3.3. They do not seem to disappear when increasing the system size, and we found them also in the uniform case. Their origin is unclear, possibly related to issues with the analytic continuation from imaginary to real time, whose subtleties are still not fully understood at present [22,23]. This effect will be studied elsewhere. Nevertheless, we stress that the analytical prediction for the curved light-cone (no signal until time t =x 0 ), which is the main message of the present paper, is verified to high accuracy, and should become exact in the limit of infinite system size.

Conclusions
In this paper we have shown that inhomogeneous Luttinger liquids with a uniform Luttinger parameter K at low energy display the universal phenomenon of curved light cones. Gapless excitations propagate along the geodesics of the metric ds 2 = dx 2 + v(x) 2 dτ 2 . Our approach, based on CFT in curved background [43], allows to quantitatively characterize the curvature of the light-cone and so to determine analytically the explicit space dependence of the sound velocity. We have checked our findings with explicit analytic and numerical calculations both in-and out-of-equilibrium for a Tonks-Girardeau gas in a harmonic potential and in systems with artificially tuned hamiltonian density. We stress however, that there is a large class of inhomogeneous systems which have uniform Luttinger parameter. This includes models with some particular internal symmetry, such as all free fermionic models, model with SU(2) symmetry or supersymmetric ones, as discussed in the text.
The main still open problem is how to generalise our approach to systems with nonuniform Luttinger parameter, such as 1d Lieb-Liniger model [80] and Fermi gases [81] in a trap potential. For those systems we still lack a complete analytic and quantitative understanding of the propagation of gapless excitations. There seems to be only a very few references that have attacked this problem: the particular case of a step-function K(x) was studied in Refs. [50,51], while partial results where obtained for correlations of the trapped Lieb-Liniger gas in [29]. We believe that this is a fascinating problem on which many physically relevant results remain to be discovered. Work in this direction is in progress and those systems with non-uniform parameter K will be studied elsewhere.