Conformal Field Theory for Inhomogeneous One-dimensional Quantum Systems: the Example of Non-Interacting Fermi Gases

Conformal field theory (CFT) has been extremely successful in describing large-scale universal effects in one-dimensional (1D) systems at quantum critical points. Unfortunately, its applicability in condensed matter physics has been limited to situations in which the bulk is uniform because CFT describes low-energy excitations around some energy scale, taken to be constant throughout the system. However, in many experimental contexts, such as quantum gases in trapping potentials and in several out-of-equilibrium situations, systems are strongly inhomogeneous. We show here that the powerful CFT methods can be extended to deal with such 1D situations, providing a few concrete examples for non-interacting Fermi gases. The system's inhomogeneity enters the field theory action through parameters that vary with position; in particular, the metric itself varies, resulting in a CFT in curved space. This approach allows us to derive exact formulas for entanglement entropies which were not known by other means.


Introduction
Low-dimensional quantum systems are a formidable arena for the study of many-body physics: in one or two spatial dimensions (1D or 2D), the effects of strong correlations and interactions are enhanced and lead to dramatic effects. Celebrated examples from condensed matter physics include such diverse cases as the fractionalization of charge and emergence of topological order in the quantum Hall effect, high-T c superconductivity, or the breakdown of Landau's Fermi liquid theory in 1D, replaced by the Luttinger liquid paradigm [1]. In the past decade, breakthroughs in the field of optically trapped ultra-cold atomic gases [2] have lead to a new generation of quantum experiments that allow to directly observe fundamental phenomena such as quantum phase transitions [3] and coherent quantum dynamics [4,5] in low-dimensional systems, including 1D gases [6,7,8,9,10,11]. These revolutionary experiments are an ideal playground for the interplay with theory, as they allow to directly realize, in the laboratory, ideal setups that were previously regarded only as oversimplified thought experiments.
On the theory side, many exact results in 1D can be obtained by a blend of methods that comprises lattice integrability [12,13] and non-perturbative field theory approaches, in particular 2D (1+1D) conformal field theory (CFT) [14,15] and integrable field theory [16]. CFT has been incredibly successful at making exact universal predictions for 1D condensed matter systems at a quantum critical point; these include the Kondo effect and other quantum impurity problems [17], and the many insights on quantum quenches [18] as well as universal characterization of entanglement at quantum criticality [19,20,21]. Entanglement entropies, in particular, are currently in the limelight, as they have become experimentally measurable both in- [22] and out-of-equilibrium [23], opening the route to a direct comparison between experimental data and many exact analytical results obtained by CFT. Indeed, entanglement entropies are usually difficult to compute in microscopic models [24,25], but their scaling limit is obtained within the powerful CFT approach by solving elementary exercises on Riemann surfaces [21].
There is a caveat in the CFT approach to 1D physics though: since it describes lowenergy excitations around some fixed energy scale (e.g. the Fermi energy), CFT does not accommodate strong variations of that scale throughout the system. This rules out a priori the possibility of tackling inhomogeneous systems, in which the relevant energy scale varies. This caveat is important, as inhomogeneous systems are the rule rather than the exception in the realm of quantum experiments: quantum gases at equilibrium always lie in trapping potentials (often harmonic) and therefore usually come with a non-uniform density; this is also the case of many out-of-equilibrium situations, such as trap releases.
In this paper, we make one step forward. We focus on the example of the free Fermi gas, in a few illustrative in-and out-of-equilibrium inhomogeneous situations. The free Fermi gas is technically simpler than interacting models, and yet it allows to draw interesting lessons that will hold more generally. We find that the varying energy scale is taken into account rather naturally in the effective field theory, in the form of varying parameters in the action. Interestingly, the metric is one such parameter, so one generically ends up with a CFT in curved 2D space. These conclusions are very general and hold under the reasonable assumption of separation of scales: there must exist an intermediate scale which is simultaneously large compared to the microscopic scale (the inter-particle distance ∼ ρ −1 , where ρ is the density), but small compared to the scale on which physical quantities vary macroscopically (of order ρ|∂ x ρ| −1 ). Indeed, at the intermediate scale , the system is well described by continuous fields because ρ −1 , and is locally homogeneous because |∂ x ρ|/ρ 1, so one knows that it corresponds to a standard (i.e. flat-space, translation-invariant) field theory. From there, it is clear that unravelling the global theory for the inhomogeneous system is a problem of geometric nature: it is about understanding the global geometric data (e.g. metric tensor, coupling constants, gauge fields, etc.) that enter the action. This is the program we illustrate with the few simple examples below. We demonstrate the power of this formalism by providing new exact asymptotic formulas for entanglement entropies.
2 The one-dimensional Fermi gas 2.1 Translation-invariant Fermi gas and the euclidean 2D Dirac action Let us start by considering a free Fermi gas in 1D in the absence of an external potential, i.e. V (x) = 0. For the reader's convenience, we briefly review the well known relation between this (homogeneous) 1D system and the CFT of massless Dirac fermions in 2D euclidean space-time. The question is: what is the proper field theory framework that captures the behavior of long-range correlations of arbitrary local observables φ 1 (x 1 , y 1 ) . . . φ n (x n , y n ) ? Here and below, the y-coordinate is imaginary time.
The starting point to answer this question is the ground-state propagator where ε(k) = 2 k 2 2m − µ is the dispersion relation and k F = 1 √ 2mµ is the Fermi momentum. Its large distance behavior is obtained by linearizing the dispersion relation around the two Fermi points k ± These two terms coincide with the right(R)-/left(L)-components of a massless Dirac fermion in 2D euclidean spacetime, ψ † R,L (x, y)ψ R,L (0, 0) = 1 x±iv F y , that derive from the action (with z = x + iv F y,z = x − iv F y): The action (4) is a holomorphic function of z. The phase factors in (3) may be incorporated into (4) with a chiral gauge transformation ψ R,L → e ±iα(x) ψ R,L . Here we do not keep track of these phase factors, as they are unimportant for our purposes, and simply discard them; these aspects are discussed in appendix A.
From Wick's theorem, it then follows that the large-scale behavior of arbitrary multi-point correlations of the original fermions can be obtained from correlators in the massless Dirac theory.

Harmonic trap and the euclidean Dirac action in curved 2d space
Consider the Fermi gas (1) in a harmonic potential V (x) = mω 2 x 2 /2. We ask: what is the underlying Dirac action? The system is now inhomogeneous, so there should be parameters in the effective action that vary with position. But what are these parameters in (4)? In order to understand this, one needs to find the density profile first. The latter is obtained from the exact solution of the microscopic problem in the thermodynamic limit, which in this case is extremely simple. The single-particle eigenstates are just those of the harmonic oscillator, and the many-body ground state is obtained by simply filling up all eigenstates with negative energies. In the thermodynamic limit (µ ω), the density profile follows the Wigner semicircle law where v F (x) = ε (k F (x)), δϕ = k F (x)δx + iv(k F (x))δy, and δϕ is its complex conjugate. Like in (3), we would like to view the terms (δx ± iv F (x)δy) −1 as the R-/L-components of a massless Dirac field. Of course, in the neighborhood of (x, y), one can always do that. But the real question is: is there a consistent Dirac theory defined globally on the entire domain (x, y) ∈ [−L, L] × R, such that its propagator has the required local behavior everywhere? If the reader is familiar with quantum field theories in curved background, they will probably have guessed that Eq. (6) is, in fact, related to the propagator of the massless Dirac fermion in curved space. The action of the latter theory in 2D is written in isothermal coordinates (z,z) and in a fixed frame (see appendix A for all details). The underlying Riemannian metric is To connect this theory to Eq. (6), notice that, in the coordinate z, the propagator behaves locally as Thus, to prove that Eq. (6) is the propagator of a Dirac fermion in a curved metric, it is sufficient to exhibit a map (x, y) → z(x, y) such that for some real-valued function σ(x, y). This equation can be solved easily. First, notice that it is equivalent to e σ ∂ x z = 1, e σ ∂ y z = iv F . Writing that ∂ x ∂ y z = ∂ y ∂ x z, we find a constraint on σ: Looking for a solution that is independent of y, we can set e σ(x) = v F (x), which implies ∂ x z = 1/v F (x) and ∂ y z = 1. The solution is straightforward: up to an additive constant, we find the complex coordinate system (z,z), which lives on the infinite strip [− π 2 , π 2 ] × R: This fixes the underlying geometry of the problem.

Application to the entanglement entropy
To illustrate the power of this approach, we exhibit new exact results for the entanglement entropies of the Fermi gas in external trapping potentials. Such calculations are, in general, difficult. But within the framework we just developed, they boil down to elementary manipulations of complex analytic functions. The Renyi entanglement entropies of a subsystem A are defined as where ρ A is the reduced density matrix of A and n is an arbitrary real number. For n → 1, S n reduces to the von Neumann entropy of the subsystem which is usually referred to as entanglement entropy.
The main property we will use in the following is that the Renyi entanglement entropies are related to the expectation values of the twist fields T n [21,26,27] which under conformal mapping share the same transformation properties of primaries with dimension with c the central charge of the CFT (c = 1 for the free Fermi gas).

Entanglement entropy of the Fermi gas in a Harmonic trap
Let us now apply Eqs. (8)- (11) to the problem of calculating the entanglement entropy in a harmonic trap.
We start with the case of a bipartition A ∪ B consisting of two semi-infinite systems, . In a homogeneous system, the Renyi entanglement entropy is [21] S n (x) is a UV cutoff, sometimes dismissed in homogeneous systems, because it simply appears in the form of a non-universal constant offset. In inhomogeneous situations, however, it is crucial to have a closer look at this cutoff: since the energy scale changes throughout the system, why shouldn't vary as well? And, if (x) varies, then it must obviously affect the dependence of S n on x. And indeed, there is a good reason why should vary with position: the continuous Fermi gas is locally galilean invariant, and the only relevant microscopic scale is the inverse density 1/ρ(x), or equivalently k −1 F (x). So the UV cutoff must simply be proportional to that scale: (x) = 0 /k F (x), for some dimensionless constant 0 .
Coming back to the harmonic potential, we now make use of the coordinate system (z,z), with z(x, y) = arcsin (x/L) + iy. The conformal mapping g(z) = e i(π/2+z) maps the z−strip to the upper half plane. To evaluate T n , we first perform a Weyl transformation e 2σ dzdz → dzdz, which changes T n into e σ∆n T n . Next, we notice that, under the conformal mapping z → g(z) = e i(z+ π 2 ) from the strip onto the upper half-plane, T n (z,z) becomes dg(z) dz ∆n T n (g(z), g(z)) uhp . The latter factor, which is the one-point function in the upper half-plane, is equal to (Im g(z)) −∆n . Putting everything together, Hence, using Eq. (14) with (x) = 0 /k F (x), we finally have for the entanglement entropy up to an additive constant and subleading corrections, which we systematically drop from now on. This gives A more complicated bipartition A ∪ B that can be considered in our framework is The calculation is straightforward, but rather cumbersome and so we report it in appendix B.1. The final result, setting Figure 1: An example: a Fermi gas trapped in a double-well potential. To obtain non trivial results, we scaled the potential with the number N of particles so as to have a finite density inside each well. Our analytic prediction for the entanglement entropy is checked against numerical data for n = 1: the agreement is excellent, and improves further when increasing N . Inset: in phase space, the quasi-classical orbitals are equipotentials enclosing an area that is an integer multiple of 2π.
which is a highly non-trivial generalization of recent results (for x 2 = −x 1 ) obtained by means of random matrix theory [28]. We checked the validity of this formula against exact finite size computations for lattice and continuous Fermi gases using the approaches of Refs. [29,30].

Fermi gas in an arbitrary external potential
The generalization to arbitrary V (x) is very simple, even though the single particle problem is not always exactly solvable for a general potential. Indeed, we are always interested in the thermodynamic limit, where the single particle states that matter are those up in the spectrum, and for which the semi-classical approximation becomes exact. Thus, focusing on the thermodynamic limit (µ ω), we proceed as follows: semi-classically, the singleparticle eigenstates correspond to the equipotentials (x, k) that satisfy the Bohr-Sommerfeld quantization rule. To get the many-body ground state, one simply fills up all eigenstates with negative energies. The equipotentials can be locally parametrized as (x, ±k F (x)), where The position-dependent Fermi velocity is v F (x) = k F (x). The metric that underlies the problem is ds 2 = dx 2 + v F (x) 2 dy 2 , one can obtain a complex coordinate system on the worldsheet just like before; the result for a general potential V (x) reads Notice that the real part of z(x, y) is the time that a massless excitation takes to arrive at point x, starting from some reference point x 0 with Re z(x 0 , y) = 0. The coordinate z always lives on an infinite strip [τ 1 , τ 2 ] × R, whose width τ 2 − τ 1 depends on V (x) and µ, and which can be conformally mapped onto the upper half-plane by z → g(z) = e iπ z−τ 1 τ 2 −τ 1 . Correlation functions are then evaluated exactly as in the harmonic case above. Formula (16) holds, and can be used to derive new exact results for the entanglement entropy. For instance, considering V (x) ∝ |x| p one recovers the results from the so-called trap size scaling [31].
One example that, to the best of our knowledge, cannot be solved with other means is that of a double well potential In Fig. 1 we compare some exact finite size calculations for this potential with our novel prediction. The figure shows that (apart from well-known finite size oscillations [32]) exact numerical results match perfectly our CFT prediction.

A non-equilibrium situation
We now show how the above framework can be adapted to deal with out-of-equilibrium situations. The most natural problem to attack is the Fermi gas (1) released from a harmonic trap. However, this problem may be solved by other methods, and it is known, for instance, that various observables obey a dynamical symmetry [33], including the entanglement entropy [30,34,35]. This symmetry relates time-dependent quantities to their time-independent, ground-state, counterpart, simply by rescaling the coordinate x to x/ √ 1 + ω 2 t 2 . Since the dynamical symmetry allows to deal with this problem in an efficient way, it is not the most illustrative example for our purposes.
Instead, we turn to a lattice gas, released from a semi-infinite box: the initial state is such that all sites x < 0 are filled, and all sites x ≥ 0 are empty (this is also known as a quench from a domain wall initial state [36,37,38,39,40,41]). At t > 0, one lets the system evolve with the Hamiltonian In Fourier space the Hamiltonian is (up to an additive constant) with dispersion relation ε(k) = − cos k. . The density at coordinate (x, y) is defined as (25). Green colors correspond to a density close to one, blue to a density close to zero. In both cases the system is said to be "frozen": observables do not fluctuate at all. Intermediate colors correspond to finite densities, and have non trivial fluctuations. The fluctuating region is a disk of radius R, with density given by (25).
The relevant regime for an effective field theory description is that of large distances and late times, in a way such that the ratio x/t is kept finite. In this limit the density profile at time t is [36] ρ(x, t) = 1 π arccos x t .
Again, the question we want to answer is: what is the effective theory that captures long-range correlators φ 1 (x 1 , t 1 ) . . . φ n (x n , t n ) in this inhomogeneous system? We expect that it should be a (lorentzian) Dirac theory in curved 1+1D spacetime. There are technical issues, however, associated with the lorentzian formulation of the problem-for instance, the metric would be degenerate, ds 2 = (dx− x t dt) 2 , and there would be no clear distinction between right-and leftmovers-, so we chose to look at the problem in imaginary time, as routinely done in quench problems in CFT [18]. In this imaginary time approach to quantum quenches, the initial state becomes a boundary condition on the two sides of an infinite strip of width 2R in imaginary time direction y [42,43,18]. Real-time correlators are recovered by first performing a Wick rotation y → it, and then sending R → 0. We focus on correlators φ 1 (x 1 , y 1 ) . . . φ n (x n , y n ) R , where y j ∈ [−R, R] (see Fig. 2 for the application to the domain wall quench). For example, the imaginary-time density profile is [44] ρ(x, y) which gives back the real-time profile (24) after performing the Wick rotation y → it and taking the limit R → 0. The density is different from zero or one only inside the disc x 2 + y 2 < R 2 ; thus, there is a phase separation phenomenon, known as arctic circle [45]. This is shown in Fig. 2. In Ref. [44], the field theory that describes long-range correlations inside the disc was unraveled: it is a Dirac theory in curved space, with euclidean metric Once we know this, it is again a straightforward exercise in CFT to compute correlation functions. As an illustration, we calculate again the entanglement entropy, for the bipartition A ∪ B with A = [−∞, x], B = [x, +∞] at times t > x (without loss of generality we assume x > 0). This can be done by performing elementary manipulations similar to those of Sec. 3.1 The main difference with the previous calculation is that, due to the presence of a finite lattice spacing, Eq. (16) needs to be modified. Namely, the position-dependent cut-off is no longer simply proportional to k F . In the homogeneous problem, it is known from the exact lattice solution [46] that the cut-off enters the formula for the Renyi entropies as sin(k F ), instead of k F in the continuous gas. As a consequence of separation of scales, in the inhomogeneous case, we thus need to replace the local cut-off k F (x, y) in Eq. (16) by sin(k F (x, y)). This leads to the following formula for the Renyi entropies in imaginary time, Here the complex coordinate z = z(x, y) must be compatible with the conformal structure of the metric (26), i.e. one must have ds 2 = e 2σ dzdz, for some function e σ . Both z and e σ are given in Ref. [44], The function (x, y) → z maps the interior of the disc x 2 + y 2 < R 2 , namely the interior of the critical region in the (x, y) plane, onto the infinite strip [0, π] × R. The strip itself can be sent to the upper half plane by the conformal map g(z) = e iz . Finally, recalling that [44] k F (x, y) = Re[z(x, y)] = arccos x elementary algebraic manipulations and the use of (27) lead to the following expression for the Renyi entropies in imaginary time The analytic continuation to real time is obtained by first substituting y → it, and then sending R → 0. This gives a formula which was guessed from numerics in Ref. [48], and was calling for an analytic derivation. We just provided this derivation, which crucially relies on the metric (26) that underlies the whole problem. In Figure 3, we report a comparison of this prediction with numerical data and the agreement is perfect, up to the usual finite-size effects. The entanglement entropy for other bipartitions can be calculated as well, but the resulting formulas are more complicated and are therefore deferred to the appendix B.2. We mention that it is also possible to study different dispersion relations, but the results are rather technical; they are reported in appendix C.

Conclusion
Inhomogeneous 1D quantum systems are difficult to tackle and this is motivating an enormous activity in order to provide exact results in some regimes, as for example the recently developed integrable hydrodynamics [49,50] (anticipated in [47]) which may have important ramifications into transport in 1D systems [51,52] such as quantum wires. To shed some light on these timely problems, we have shown here, with a few free fermion examples, that CFT methods may be extended to attack this class of problems. The key assumption of our work is separation of scales: the system is locally homogeneous (but only locally), which is also the working limit of all approaches constructed so far [47,49,50]. In this regime, one can write a field theory action with varying parameters. In the examples tackled in this paper, we found that the metric tensor should vary, leading to CFT in curved space. In particular, this new approach allows us to compute in a simple manner the entanglement entropies of these inhomogeneous systems (both in-and out-of-equilibrium) which are otherwise very difficult (in most cases impossible) to obtain with other methods. It is important to stress that the background metric and the inhomogeneous cut-off are non-universal functions and must be viewed as inputs of the formalism. They should be obtained a priori by a proper microscopic computation.
We believe anyway that the results of this paper should open the door to several new developments. For example, by following the CFT approach of Ref. [53], our method could be used not only to determine the entanglement entropy, but also the entanglement negativity which is a proper measure of entanglement in mixed states [54] and it is not yet known how to compute it exactly for free fermionic systems. Another interesting development would be to recover by random matrix techniques [55,28] the results we obtained for the entanglement entropy in a harmonic potential, possibly unveiling new structures of the random matrices. Having exact results also for arbitrary trapping potential could also help understanding whether these general cases could be tackled by random matrix techniques.
From a physical point of view, the main open problem is how to describe inhomogeneous interacting 1D systems, most importantly Luttinger liquids such as Heisenberg spin chains and Lieb-Liniger gases. In these cases, one expects also the coupling constant (or Luttinger parameter K) to vary with position in spacetime. Fixing such parameters is a challenging problem. We hope that this paper will stimulate activity in that direction.

A The Dirac action in curved space-time
Here we explain the form of the Dirac action in a curved space-time, that was used in the main text. The most generic form of the Dirac action in a curved two-dimensional eulidean space-time is given by where we chose γ 1 = −σ 2 , γ 2 = σ 1 , γ 5 = σ 3 andΨ = Ψ † γ 2 . The σ α are the usual Pauli matrices. In this representation the two components of the spinor Ψ are the chiral components ψ R and ψ L ; the function A µ are the vector and axial gauge field associated with the U (1) gauge symmetry and the (anomalous) chiral symmetry. The matrix e µ a is called a tetrad: it maps the tangent space of the manifold into R 2 while preserving the inner product. It satisfies g µν e µ a e ν a = δ ab , where g µν is the metric tensor and δ ab the flat metric. In two dimensions, it is convenient to use complex coordinates z = x 1 + ix 2 andz = x 1 − ix 2 . The metric is conformally flat and off-diagonal, the line element thus reads ds 2 = e 2σ dz dz. The tetrad is diagonal if one chooses complex coordinates both in the tangent space and in R 2 . Its only non-vanishing components are complex conjugated and have with modulus e −σ . Using finally √ g = e 2σ we can then rewrite the Dirac action in complex coordinates as where A (v/a) andĀ (v/a) are the complex z,z components of the vector and axial potential. The rotation of angle θ in (33) between the tangent space of the manifold and the flat euclidean space does not alter any of the conformal maps discussed in the main text. This holds because the twist fields needed to compute the entropy are spinless [27]. The extra phase factors in the fermionic propagators may be restored by performing a U (1) gauge transformation and a chiral transformation. The former acts on the chiral fermions as ψ R/L → e iϕ ψ R/L , and on the gauge field as The entropies of a finite interval [x 1 , x 2 ] are slightly more complicated to calculate than the case of the semi-infinite line (cf. Eq. (17)). Following Ref. [21], the Renyi entanglement entropy of order n is related to the two-point function of the twist field T n on the desired geometry. Here we work out the generalization to inhomogeneous systems and in particular for the case of fermions in a harmonic trap. In this case, the two-point function can be related to the one of the upper half plane by the combined action of a Weil transformation and the conformal mapping g(z) = e i(z+π/2) from the strip [− π 2 , π 2 ] × R to the upper half-plane. This allows us to write such two-point function in the z coordinate as The fieldT n is the conjugated twist field [27,21]. The calculation of two-point functions of twist fields in the upper half plane is, in general, a very challenging problem (indeed by images trick, it can be turned into a four point function in the complex plane which have been considered e.g. in Refs. [56,57]). However, in the case of the massless free fermion field-theory, this two-point function in the half plane simplifies considerably and our desired object boils down to [58,30] (up to unimportant multiplicative constants) T n (z 1 )T n (z 2 ) = e −σ(z 1 ) dg(z 1 ) dz 1 e −σ(z 2 ) dg(z 2 ) dz 2 |g * (z 1 ) − g(z 2 )| 2 Im(g(z 1 ))Im(g(z 2 ))|g(z 1 ) − g(z 2 )| 2 ∆n . (35) In our setup z 1 = arcsin(x 1 /L) and z 2 = arcsin(x 1 /L). Using again g(z) = e i(z+π/2) and (11), we obtain where S n (x i ) is given in Eq. (17). This is the result reported in the main text.

B.2 The domain-wall quench
For the case of a finite interval, the Renyi entanglement entropy are easily obtained by using Eq. (35) and the conformal mapping from the strip [0, π]×R to the upper half-plane g(z) = e iz . Setting y 1 = y 2 ≡ y (i.e. equal imaginary times) and using (28), we get S n (x 1 , x 2 , y) = S n (x 1 , y) + S n (x 2 , y) + n + 1 12n ln . (37) Analytically continuing the result (37) to real time and using the rescaled variables This formula is valid in the regime t > |x 1 |, |x 2 |. Notice that this is very similar to Eq. (36), and indeed it has the same dependence of the rescaled variables ζ i . However, the leading dimensional term in one case scales like ln t while in the other as ln L 2 (times (n + 1)/(6n) in both cases). This is somehow reminiscent of a similar anomalous scaling found in local quenches [59] and it is unclear whether there is a connection between the two.
C The domain wall quench for a more general dispersion relation C.1 Imaginary time treatment Here we come back to the entropy of a single interval (−∞; x] following a domain wall quench, but with a more general dispersion relation. We consider Hamiltonian (23) with a dispersion relation of the form which contains only odd Fourier modes. The reason for choosing this special form is mainly technical. As we shall see it makes explicit computations much simpler due to the symmetries ε(−k) = ε(k) and ε(π − k) = −ε(k).
In real time the stationary phase equation governing the long time dynamics is [41] v and may be used to compute systematically all correlation functions. Now, due to the form (39), if k s (x, t) is a solution, then so is π − k s (x, t). In the following we will assume that there are at most two solutions k s (x, t) and π − k s (x, t) to the previous equation. Under these assumptions, the density profile is As is emphasized in the main text and above, it is convenient to think of this problem in imaginary time. The stationary phase equation for the imaginary time problem reads [44] the Hilbert transform of ε(k). In case there are only a finite number of Fourier modes in the dispersion, solving (42) amounts to solving algebraic equations in e ik , which can be done in principle systematically. Now we make the extra assumption that there are only two solutions k = z(x, y) and k = −z * (x, y). The new variable z lives in an infinite strip 0 ≤ Re z ≤ π. The metric then reads The conformal distance to the boundary may be obtained by mapping the infinite strip to the upper half plane. Such a mapping is once again provided by g(z) = e iz , so we obtain Img(z) = e σ(x,y) sin (Re z(x, y)) .
The reasoning to get the entanglement entropy should be clear at this point: first we get κ and Q from the above set of equations, and then we use them to compute the conformal distance to be plugged in the appropriate equation for the entropy. For example in the case of the standard cosine dispersion relation, using this procedure, we find the already known result d(x, y) = (R 2 − x 2 − y 2 )/ R 2 − y 2 , which can then be continued to real time y = it, R → 0. It is important to stress that the concept of conformal distance only makes sense in imaginary time: in principle only after computing d(x, y) in imaginary time we are allowed to make the analytic continuation to real time.
In practice, one can avoid finding the general solution of the above system of equations. Indeed, because of their analytic structure, we can perform the analytic continuation at the level of Eqs. (47) and (48), relaxing the requirement that κ and Q are real. Plugging y = it and taking the limit R → 0, we find that Q = iπ/2 is a trivial solution to (48), due to the absence of even Fourier modes. The variable κ is then the solution of x t = n≥1 (−1) n+1 cos[(2n + 1)κ] = v(κ + π/2), and so we get κ = k s − π/2. The analytic continuation of the conformal distance becomes lim R→0 d(x, it) = t dv dk k=ks cos k s = tv (k s ) cos k s .
Putting back the contribution from the UV cutoff, sin k F = cos k s , we finally obtain S n (x, t) = 1 12 1 + 1 n ln t v (k s ) cos 2 k s .
Let us briefly comment on the case with more than two solutions to the stationary phase equation (40). With the form (39) of the dispersion they always come in pairs. We label them as k (i) s and π − k (i) s for i = 1, . . . , p. p may be interpreted as a number of Fermi seas. As before we compute the entanglement entropy using field theory, and add to this the contribution coming from the position-dependent UV cutoff. The generalization of our field-theoretical framework to account for several Fermi seas is straightforward: we simply introduce p different species of Dirac Fermions, one for each Fermi sea (pair of stationary modes). Being non-interacting particles, all the Dirac species contribute independently to the entanglement entropy. Obtaining the position-dependent cutoff is more tricky, but can nevertheless be done using exact equilibrium results obtained for several Fermi seas 1 . We refer to Ref. [60] for the details. The lattice cutoff takes now the more complicated form The final result for the entropy reads

C.2 An explicit example
Let us now work out an example where the formula (51) may be computed explicitely. We consider a dispersion of the form ε(k) = − cos k + α 9 cos(3k) , 0 < α < 1.
The restriction on α ensures that there are at most two solutions to the stationary phase equation, k s (x, t) and π − k s (x, t). k s may be obtained by solving a cubic equation:  : von Neumann entanglement entropy as a function of x/t after a quench from the domain wall initial state, for dispersion relation ε(k) = − cos k + α cos 3k and α = 1/10. he numerical simulations are performed using a finite system with 4096 sites. Data for t = 512 are reported (black circles) and compared to our analytical prediction (red lines): the agreement is nearly perfect. Inset: corresponding density profile ρ(x, t) after the quench.
The density profile in the inhomogeneous region |x/t| ≤ 1 + α/3 is then given by and the entropy is obtained by plugging (55) in (51). In Fig. 4 we report the numerical data for the von Neumann entanglement entropy for the dispersion relation with α = 0.1: this shows perfect agreement with our prediction.