Resolving modular flow: a toolkit for free fermions

Modular flow is a symmetry of the algebra of observables associated to spacetime regions. Being closely related to entanglement, it has played a key role in recent connections between information theory, QFT and gravity. However, little is known about its action beyond highly symmetric cases. The key idea of this work is to introduce a new formula for modular flows for free chiral fermions in $1+1$ dimensions, working directly from the \textit{resolvent}, a standard technique in complex analysis. We present novel results -- not fixed by conformal symmetry -- for disjoint regions on the plane, cylinder and torus. Depending on temperature and boundary conditions, these display different behaviour ranging from purely local to non-local in relation to the mixing of operators at spacelike separation. We find the modular two-point function, whose analytic structure is in precise agreement with the KMS condition that governs modular evolution. Our ready-to-use formulae may provide new ingredients to explore the connection between spacetime and entanglement.


Introduction
Recent years have witnessed significant progress in our understanding of the role that entanglement, as well as other ideas from quantum information theory, play in the context of high energy physics, including Quantum Field Theory (QFT) and gravity. A remarkable example of this exchange of ideas between different research areas is the Ryu-Takayanagi formula [1] for the entanglement entropy in the AdS/CFT correspondence, which generalises the Bekenstein-Hawking area law for the black hole entropy.
Despite the many contexts in which modular flow appears, there are very few cases where its action is explicitly known. In the general context of QFT, the vacuum modular flow in a Rindler wedge is fixed by Poincaré symmetry alone [42], while conformal symmetry fixes it for diamond shaped geometries [43]. Anything beyond these cases, be it the choice of another state or a different region, depends on the details of the theory under consideration and is largely unknown. Many discussions are concerned with universal properties of modular flows, or deal with highly symmetric configurations where the flow has a geometric (local) interpretation. In a generic case however, we expect to see many forms of non-localities. An initially local operator that is subject to non-local modular flow acquires contributions from spacelike separated regions with increasing modular time. Quantum entanglement of spacetime seems to play a crucial role here. Therefore it is of great interest to obtain further detailed understanding of explicit realizations of modular flow for specific cases.
The example of the chiral fermion is rich enough for understanding non-universal behaviour in detail, but still simple enough for explicit computations. Recently, this has lead to novel results on the modular Hamiltonian and relative entropy for disconnected regions [37-41, 44, 45]. In this paper we go beyond these studies and derive results for the modular flow itself. We provide explicit formulae that may find direct applications in studies of fermionic entanglement.
Let us briefly state our main results. We consider density matrices reduced to an arbitrary set of disjoint intervals V = n [a n , b n ]. Modular flow of a fermion operator localised at y ∈ V is given by the convolution σ t ψ † (y) = V dx ψ † (x)Σ t (x, y) , (1.1) where the kernel is a function of the correlator, Here the reduced propagator G| V is understood as a linear operator acting on smooth functions via convolution. This approach avoids the computation of the modular Hamiltonian and instead directly yields the flow. We determine this kernel, and the explicit formulae for the modular flow are given in (3.15) (for the plane or cylinder (A)), (3.24) (cylinder (P)) and (3.30) (torus), which are illustrated in figures 4, 5, and 6 respectively. Here, the notation P and A refers to periodic and antiperiodic boundary conditions for the fermions on the spatial circle, in other words the Ramond and Neveu-Schwarz sector. As a second novel result, we explicitly compute the modular two-point function, which is given by G mod (x, y; t) = ψ(x) σ t (ψ † (y)) . (1. 3) The final results are, in the same order as above, (4.8), (4.16) and (4.21). A remarkable feature of this approach is that, although σ t (ψ † ) generically involves solving higher-degree polynomials or transcendental equations, the modular two-point function can be determined analitically by direct integration, without the need of solving such equations.
An important concept in our paper will be that of different degrees of locality of the modular flow. This can be understood directly from (1.1). We call a flow completely nonlocal if the kernel Σ t (x, y) is a smooth function of x, y supported on the entire region V , since it mixes operators along the entire region. If however Σ t (x, y) ∼ δ (f (x, y)) for some function f , the integral in (1.1) will localise to a discrete set of isolated contributions, namely the zeroes of f . Generically, these solutions will be non trivial, in the sense that x = y at t = 0. We call these solutions bi-local, since they couple pairs of distinct points. Finally, if there is a solution such that x = y at t = 0, we call it local. We will use this terminology throughout the text. As we shall see, one of the essential features of our operator flow (1.1) is that it changes its locality properties depending on the temperature and the spin boundary conditions. In turn, this manifests itself in the structure of poles and cuts of the modular correlator (1.3).
The paper is organized as follows. In section 2 we specify the objects that we aim to compute: the modular flow of fermion operators and the modular two-point function.
To this effect, we first review some basic notions of Tomita-Takesaki theory, a mathematical framework relevant for local quantum field theories. In section 2.2 we introduce the particular physical system we focus on -the free chiral fermion -together with the necessary tools we will use to study modular flow. The main ingredient is holomorphic functional calculus, including the method of the resolvent. Section 2.3 includes a list of all known fermion resolvents. Section 3 contains the first new results of this paper: We apply the above techniques to find the modular flow of the fermion operator for all cases considered. Section 4 presents our second important result, the modular two-point function. We verify that it obeys all required properties, such as analyticity and the KMS condition. We conclude by a summary and future directions in section 5.

A lightning overview of modular theory
Before we specialize to the free fermion, let us first recall some basic notions of (Tomita- Takesaki where J is antiunitary and ∆ is positive. We denote J and ∆ as the modular conjugation and modular operator, respectively. The importance of these operators stems from Tomita's theorem, which states that i.e. J intertwines R and its commutant R (the set of bounded operators in H that commute with those in R) while the modular flow preserves R. Hence we see that modular flow is a symmetry of R arising just from the algebraic structure. The generator K of this symmetry, in the sense that admit analytic continuations to the strips −1 ≤ (t) ≤ 0 and 0 ≤ (t) ≤ 1, respectively, where they satisfy This is to be compared with the original KMS condition for a state ρ = Z −1 e −βH of inverse temperature β: Denoting time evolution by α t (O) := e itH Oe −itH , we have (2.8) Therefore, σ t behaves like a time evolution with respect to the Hamiltonian −K in a state of inverse temperature −1. In other words: If we only consider operators in R, the vector |Ω behaves like a thermal state ρ R = Z −1 e −K . This coincides with the definition of modular flow in terms of "reduced density matrices" [8]. As a word of caution, we would like to mention here that reduced density matrices do not exist as operators in a genuine quantum field theory (QFT) due to the universal divergence of vacuum entanglement [3,4,8]. However, in the following, we will use the formal analogy to finite dimensional systems (where reduced density matrices do exist) to derive relevant formulae. While not entirely rigorous, this method has proven to be very useful in the past and was confirmed in many cases [41,46].
We now specialize to QFT in flat spacetime: In the Haag-Kastler approach to QFT, a von Neumann algebra R is associated to each (causally complete) region in space time, typically denoted by the same symbol. This algebra can be thought of as consisting of the set of (bounded) operators that have support in R. Since the associated modular flow preserves this algebra and, hence, the associated region, it is tempting to ask, to which extent it has a geometric or physical meaning. Remarkably, this question has an affimative answer in the following scenario: If |Ω is the vacuum state and R a Rindler wedge, then K is nothing but the (approriately scaled) generator of Lorentz boots that preserves this wedge [42]. Furthermore, we can sometimes use additional symmetries (e.g. conformal symmetry) to generalize this geometric action to other regions such as a lightcone, double cones, or even to other states such as the thermal state [43].
Finally, let us elaborate further on the KMS condition in the context of fermionic theories. As an example, we consider eq. (2.7) for the case of two field operators O 1 = ψ(x) and O 2 = ψ † (y). As stated, the two functions defined on (2.6) are analytic when restricted to the lower and upper unit strips respectively. Therefore, it is natural to define the modular two-point function [3,41,47] This function, defined via a different function in each strip, satisfies the following antiperiodicity due to the KMS condition (2.7) allowing to continue G mod to arbitrary non-integer imaginary parts of t. The reason to define (2.9) is because of its relation to the anticommutator. While G mod is analytic in the lower and upper strips by construction, the interesting question regards its regularity properties of along (t) = 0 (and, hence, at all integer imaginary parts due to antiperiodicity). Now it is easy to see that the variation of G mod along the real axis is given by the anticommutator between a field and the modular evolved field: This condition becomes particularly useful in cases where the modular evolution is given by a smearing of the field (this will be precisely the case for gaussian free fermion states) (2.12) with which (2.11) becomes by the canonical anticommutation relation {ψ(x), ψ † (y)} = δ(x − y). Equivalently, we can use the antiperiodicity (2.10) to rewrite this purely in terms of the function on the lower strip, This relation is important because it relates the analytic structure of the modular correlator to the locality properties of the modular flow, via the Kernel Σ t of the operator flow. If the flow under consideration is local, i.e., if Σ t (x, y) ∝ δ(x − y), the right hand side vanishes almost everywhere and we obtain antiperiodicty of the two-point function in imaginary time. On the other hand, this means that every failure of such regularity (e.g. a branch cut of G mod ) is a clear sign of non-locality in the modular flow.
In section 4 we will explicitly determine G mod (t) for the two-dimensional chiral fermion. We will confirm that for the cases where modular flow is local or bi-local (such as a single interval on the plane, antiperiodic cylinder, or torus), the modular correlator restricted to the lower strip yields a function that is analytic everywhere away from isolated simple poles. The salient example where analyticity fails is a single interval on the periodic vacuum on the cylinder, where the flow is completely non-local and G mod has branch cuts.

Modular flows for free Fermions
For the rest of the paper, we restrict to Gaussian states in a Fermionic theory. Also, we are working on a Cauchy slice and consider only subregions V of that slice. In this case, we can formally decompose the modular Hamiltonian as (2.15) where V c is the complement of V in the Cauchy slice. The absence of mixing terms between V and V c reflects the fact that modular flow preserves R and R , which are associated to V and V c respectively. Mathematically, the above Kernels k, k c have to be understood in the distributional sense, i.e. they are only defined as integrated against suitably smooth test functions. For the remained or the text we restrict to k, as the calculation of k c is completely analogous.
To derive an explicit formula for k, we require that the "reduced density matrix" Z −1 e −K reproduces the correct expectation values for operators with support in V . By Wick's theorem it is enough to reproduce the propagator G(x, y) := Ω|ψ(x)ψ † (y)|Ω (2.16) in the subregion x, y ∈ V , which allows to derive the relation [46,48] where G| V is the restriction of G to V . In a similar manner, we arrive at Here and below, we use the compact notation Σ t omitting its space-time dependence, but it should always be kept in mind that it is a linear operator acting on functions. Similarly for other operators such as the resolvent. The problem is thus reduced to computing functions of the (restricted) propagator G| V . Since this is a bounded operator -its spectrum is contained in the interval [0, 1] -we can use functional calculus to write where 1/(λ − G| V ) is the resolvent of G| V and γ denotes that the integral is to be done counter-clockwise along a contour that tightly wraps around the spectrum [0, 1], as shown in figure 1. Eq. (2.20) can easily be seen to be correct in an eigendecomposition of G| V and implies that the resolvent 1/(λ − G| V ), as a function of λ, is analytic in a neighbourhood of [0, 1], but not along the interval itself: If the spectrum of G| V is discrete, we expect a simple pole whenever λ approaches an eigenvalue. For continuous portions of the spectrum, this culminates in a branch cut. In any closed form expression of the resolvent, such a branch cut will only be visible as a pair of branch points and one might be tempted to assume that the precise location of the cut is indeterminate. However, as just discussed, one should always keep in mind that the branch cut will always be situated alongside [0, 1]. Let us compare eq. (2.20) with the general definition of a function of the kernel G| V , given by where E λ is the spectral measure of G| V . Assuming the contributions to eq. (2.20) at 0 and 1 vanish, we obtain which characterizes the spectral measure completely. Note that the requirement of vanishing contributions at 0 and 1 also imposes a regularity constraint on the function f . For the computations in section 3, it will turn our that this contraint is violated and we will have to work with eq. (2.20) directly. To proceed any further, it is necessary to find the resolvent for the state and region under consideration. To this end, we can make the ansatz Notice that while this equation is valid for fermions in arbitrary dimensions, the solutions are only known in two dimensions, which is the focus of the next subsection.

Resolvent for the chiral Fermion
Due to recent developments [39], (2.24) is well understood in the special case where we are dealing with a chiral fermion in one dimension and V = n [a n , b n ] is a finite union of disjoint intervals. This is because there, the propagator is a Cauchy kernel, hence the integral equation can be reduced to yet another complex analysis problem, where the resolvent has a branch cut along V . We omit details here and just state the results in the following cases, classified by the domain/periodicities of G: • No periodicity (the entire complex plane)-the corresponding propagator is given by The solution is (2.27) • One periodicity, taken to be 1 without loss of generality (the complex cylinder) -the propagators are depending on the choice of antiperiodic or periodic boundary conditions, respectively. The corresponding solutions are , (2.31) with the total length L = n (b n − a n ) of V and (2.32) • Two periodicities 1, τ (the complex torus) -here, the propagators are with ν = 2, 3 denoting the periodic-antiperiodic (PA) and antiperiodic-antiperiodic (AA) boundary conditions, respectively. The conventions for Jacobi theta and Dedekind eta functions are the same as in [49]. The solutions of (2.24) are now where h is defined by e 2πh := − 1−λ λ and Note that G (ν) (x, y; τ, Lh) is the propagator of a state with chemical potential Lh, i.e. we have the series representations where the symbol denotes that the sums have to be ordered symmetrically to ensure convergence.

Modular flow of operators
In this section we will compute explicitly the modular flow of the fundamental field, σ t (ψ † ) from (2.18). This is a basic building block that allows to compute the flow of composite operators. As explained in section 2.2, the task reduces to determining the kernel and decomposing the resolvent as in (2.23) we have As it stands, this integral is not completely well defined, as the integrand is both divergent and highly oscillating around the branch points. However, this should come as no surprise: since we know that Σ t (x, y) represents a distribution, we expect the appearance of Dirac delta distributions. Thus the strategy is to regularise the integral, evaluate it, and finally remove the regulator and identify the remaining distributions.
Here the analytic structure of the integrand is crucial. In addition to the cut associated to the resolvent -the last factor in (3.2)f (λ) has introduced another cut, branched over the same endpoints. The latter cut can be freely chosen as long as it does not overlap with the former. For simplicity, we choose it to run along the real complement, In appendix A we provide a rigorous treatment of this integral and evaluate it by residues. Here instead, we proceed with a more straightforward but nevertheless equivalent approach. As explained in appendix A, a standard regularisation consists of avoiding the poles at λ = 0 and λ = 1 by shifting them slightly into the complex plane. As a consequence, the integral with the first term in the square brackets, proportional to 1/λ, vanishes. This can be seen as follows. This term has a branch cut along (−∞, 0) ∪ (1, ∞) and a pole at λ = 0, but everywhere else is holomorphic. In particular, the contribution around λ = 1 vanishes due to the KMS requirement. Since both boundary terms vanish, the integral vanishes. For more details, please refer to appendix A. Thus, we are left with Finally, F λ takes a different form depending on the topology and boundary conditions chosen. We consider them case by case.

Plane
We start with the simplest case. For the vacuum state on the plane (2.26), F takes the form where we have introduced the shorthand notatioñ and throughout the text we uset =t(x, y) and omit the spacetime dependence, which should nevertheless be kept in mind. The notation will become clear shortly, as we will see thatt plays a role closely analogous to modular time t. Also, it is important to realise that here, the propagator G(x, y) has no dependence on λ and therefore can be pulled out of the integral (notice however that this will not hold for the cylinder (R) or torus). Introduced back in (3.3), one obtains where we have defined the integral (see fig. 2) In order to exploit the symmetry and simplify the integral, it is useful to change variable to z = (1 − λ)/λ, which maps the cut along [0, 1] to R + , and R \ [0, 1] to R − , while the image of the contour wraps positively around R + , see fig. 3. To make the branch cuts explicit, we use −π < arg z < π for the cut along R − and 0 < arg z < 2π for the one along R + . The integral (3.7) now reads The discontinuity along R + implies that just above and below the cut we have As mentioned above and discussed in detail in the appendix, the regularisation makes the integrand bounded in a neighbourhood of the origin z = 0, so we can neglect the contribution of the boundary point. Putting everything together, (3.8) can be formally represented as Now we proceed to the regularisation. This integral receives divergent contributions from both z → 0 and z → ∞. Both can be made finite by restricting the integration to z ∈ (e −2πm , e 2πm ) where we eventually take m → ∞. The result is Again, we remind the reader that a rigorous derivation of this via residue analysis is presented in appendix A. This expression is actually familiar. If t +t(x, y) = 0, the fraction is bounded but wildly oscillating, and therefore vanishes when integrated agains regular test functions.
In the vicinity of t +t(x, y) = 0 instead, the fraction diverges as m → ∞. This is the standard Dirichlet kernel representation of the Dirac distribution, so As stated above, it is clear thatt is indeed playing a role somewhat analogous to modular time itself. Putting back all together into (3.6) and replacing (3.5), we learn that the kernel associated to the action of modular flow is (3.13) whose support is given by the solutions of This equation, and its solutions, play a fundamental role in our analysis. It will determine which points are non-locally coupled via modular flow, as well as the magnitude of their coefficients. For a fixed y, we shall call x = x (y) the solutions for x, where the discrete index labels the different intervals in V . The most important property of the function Z is that it increases monotonically from −∞ at each left endpoint a j to +∞ at the right endpoints b j . This guarantees that there will exist one solution to this equation per interval. The action of modular flow finally reads where as before x are the solutions of (3.14). We omit the absolute value in the denominator since Z(x) is monotonically increasing in each interval. In order to gain some intuition, we illustrate these results below with some simple examples.
To conclude here, let us note what happens in the limit t → 0. At zero time, the kernel Σ t in (3.1) must reduce to the identity, localized at x = y. Now the prefactor in (3.15) vanishes linearly with time. Since the propagator has a (unique) simple pole at coincident points, all but the 'local' solution, which obeys x → y at t → 0, vanish in (3.15).
Rindler space. This is the best known explicit case, which obeys a universal formula for the vacuum of any QFT on the Rindler wedge [42]. Physically this corresponds to the standard Unruh effect, where modular evolution is nothing but translations along the worldline of observers with constant acceleration. Here the entangling region is V = R + , which can be obtained by taking the limit b a of the single interval on the plane (2.27), yielding The unique solution to (3.14) is x 1 = e −2πt y, which inserted back into (3.15) leads to the geometric flow σ t ψ † (y) = e −πt ψ † e −2πt y , (3.17) the prefactor being due to the transformation law of a spin 1/2 field under Lorentz boosts.
Multiple intervals on the plane. Here the entagling region is an arbitrary set of disjoint intervals V = ∪ n i=1 (a i , b i ). This was the case solved in the seminal work [37]. However it is important to note how their strategy differs from ours. In [37] the authors first derived the modular Hamiltonian K V = − log ρ V and next used the associated Heisenberg equation This yields a set of coupled differential equations relating the different ψ(x (t)). On the other hand, we computed modular flow directly in terms of the resolvent. This avoids using the modular Hamiltonian itself, and the need for the differential equation, and hands at once the solution.
For completeness, and in order to compare to the new results, we illustrate this case for two intervals (a 1 , b 1 ) ∪ (a 2 , b 2 ). Again, the solutions to (3.14) are essential. In this case, and therefore (3.14) leads to a second degree equation, which can be readily solved but the expression is rather cumbersome. We plot these solutions in Fig.4. This was the first known case of a bi-local or quasi-local modular flow that could be solved analytically [37]. The most important feature is that the flow contains involves two kinds of terms. The local solution lives in the same interval as y, and is continuously connected to it at t = 0. But there are also bi-local terms, one per interval. Also notice that due to chirality, the solutions move towards the left as modular time t evolves, converging to the left endpoints asymptotically, and similarly go to the right endpoints as t → −∞. This simple example illustrates an important point. In general, the configuration on the plane with n intervals involves polynomials of degree n. Thus, solving explicitly for the modular flow of the fermion operator becomes quickly hopeless as we increase n. This becomes even more involved for other states like the cylinder or the torus. Now, since these terms will show up in the modular flow of composite operators, it would seem implausible to find any closed analytic results for those cases. However, as we will show in section 4 with the modular two-point function, there exists a remarkable way to circumvent the need to solve (3.14).

Cylinder
The vacuum state for the fermion on the cylinder possesses two spin sectors, the antiperiodic (A or Neveu-Schwarz) and the periodic (P or Ramond), depending on the boundary conditions we choose for the fermion along the circle. The modular flow for the antiperiodic sector (for any number of intervals) is identical to that of the plane, provided we use the appropriate propagator (2.28) and Z(x) given in (2.32), and therefore we will not elaborate on it. However, this behaviour changes dramatically when we consider the periodic sector.
The periodic sector on the cylinder provides an example of how the present method allows to go beyond previous results in the literature. This case must be considered separately, since the resolvent (2.31) contains an extra term, due to the presence of a zero mode. The first term of (2.31) is of identical form as the one derived in the previous section, namely (3.6), where again S(x, y) is given by (3.12) and the corresponding correlator on the cylinder (P), (2.29). This allows us to decompose the modular evolution as where the extra term, associated to the Ramond sector, can be written in the variable z = (1 − λ)/λ as where again we definedt = Z(x) − Z(y). Using again (3.8), this can be brought into the more convenient form Although the last factor in the integral does not possess a multiplicative branch cut, one can bring it into such a form using the identity where we identify y = (−z) L for L ∈ (0, 1). Then, the same steps leading to (3.12) yield (3.23) As we discuss below, this result is quite remarkable. While the first term -given by (3.13) -produces a local flow, the second one does not and leads to complete non-locality. Indeed, as explained in 2, any contribution to the kernel Σ t (x, y) that is not localised (in the sense of being proportional to δ(x − y)) will produce a discontinuity in the modular two-point function, in accordance with the KMS condition (2.10).

(3.24)
The first term is very similar to the flow on the plane or cylinder (NS) of the previous section, but with the correlator shifted by a constant of 1/2. This factor is not very significant as it can be reabsorbed into the integral; it will show up again in section 4.2. Again, this flow always involves a local (geometric) term, and additional bi-local (quasi-geometric) couplings whenever V contains more than one interval. In particular, the solution for the x is identical to that on the antiperiodic sector, since this depends only on the function Z(x), which does not depend on the periodicity. Notice that, as mentioned above, at t = 0 all terms vanish, except the one involving the propagator. The second piece is another important result of this work. It constitutes an example of a continuously non-local modular evolution. The operator ψ † , initially localised at y at t = 0, receives contributions from the entire interval as modular time evolves. The properties of the coefficient sinh(π(t +t)/L) −1 are again determined by those of Z(x): since this function increases monotonically and diverges at the endpoints, the coefficient vanishes at ∂V , and has an asymptote at the solutions of (3.14). Just as with the bi-local flows of section 3.1, the asymptote moves monotonically from b i at t → −∞ to a i at t → ∞. We illustrate these features in figure 5, for the case of a single interval. The corresponding modular Hamiltonian for the periodic case was first derived in [38]. Although the continuously non-local flow appears to have a different nature than the bi-local quasi-geometric terms, it actually results as the zero temperature limit of bi-locality on the torus, as shown in [39]. We comment on this in the next subsection.
For multiple intervals on the periodic sector, the flow has three components: the local piece, the bi-local term (one per interval as in the plane), and the completely non-local one, which couples a given point to all intervals.

Torus
The resolvent on the torus possesses a qualitatively new behaviour. As shown in [39], the modular Hamiltonian involves a discrete but infinite set of bi-local couplings. The same behaviour is present for the modular flow and the modular correlator as discussed in section 4.3. See also [40] for related work.
The novelty here is that the propagator appearing in F λ (see eq. (2.34)) carries an explicit λ-dependence, and therefore does not factor out of the integral in (3.3), as in the previous cases. Nevertheless, we can still find an analytic solution. In (2.37) we noticed that the propagator appearing in this case can be re-interpreted as the correlator with a chemical potential turned on and has an associated series representation, (2.38) and (2.39) for each spin sector. Using this, and again essentially the same steps as in section 3.2 for the cylinder for t = 0 (dashed-opaqued) and t > 0 (solid blue). The novelty is that in addition, this case contains a completely non-local term, as seen in the second term of (3.24), represented by the dot-dashed line that is proportional the non-local kernel. This vanishes at the endpoints and diverges at the position of the local contribution. where the sum runs over k ∈ Z for ν = 2 and k ∈ Z + 1/2 for ν = 3. Finally, using the periodic/antiperiodic representations of the Dirac distribution, we find that the kernel for the modular operator is given, for ν = 2, 3, by where we have replaced τ = iβ again to emphasize that the argument of the Dirac distribution is strictly real. Here again the argument of the distribution plays a fundamental role. The support of Σ t locates at the roots of For every integer k and interval , (3.28) has a single solution, x ,k = x ,k (y), since again Z(x) is monotonically increasing from −∞ to +∞ within each interval. This is illustrated for a single interval in figure 6. Therefore modular flow connects any given point y ∈ V to an infinite set of other points in V . Thus we can write Finally, the modular flow on the torus is then . (3.30) This result illustrates the 'infinite bi-locality', as was already understood for the modular Hamiltonian [39,40]. As can be seen explicitly from the above equation and depicted in figure 6, the modular flow of a field localised at y receives contributions from a discrete but infinite set of points within V , even for a single interval. As discussed in [39], this structure illuminates the behaviour at zero temperature. Indeed, when β → ∞, the solid curve becomes a straight line with infinite slope, and the red dots in figure 6 form a regular partition of the interval. For ν = 2, this reproduces exactly the definition of the Riemann integral, which leads to the continuous integral for the periodic sector on the cylinder as shown in (3.24).
It is worth emphasizing once more how explicit (3.30) is. Indeed, whereas alternative approaches [40] write the modular flow in terms of a set of infinitely many coupled differential equations, our method yields at once the full solution.

Modular two-point function
In this section we illustrate the power of the tools laid out in section 2 and explicitly calculate the modular two-point function defined in eq. (2.19). Following the distiction of cases from section 2.3 we may work out the different results for the plane, the cylinder with (anti-)periodic boundary conditions and finally the torus. We will explicitly calculate G mod (t) in the lower strip −1 < (t) < 0, and later comment on the continuation to the complex plane.
On a formal level, the modular correlator is already determined by the results of the previous section. All we have to do is take the expectation value of σ t ψ † with another fermionic field, on the global state considered, where the modular flow of the operator σ t (ψ † (y)) has been computed for all mentioned cases in section 3. For instance, in the simplest case of the plane or cylinder (A) this yields where again x k (y; t) are the solutions of (3.14). In practice however, this formula is not very useful. As mentioned above, the problem of finding the roots x k is in general very hard. Moreover, as mentioned in section 2, the modular correlator must be analytic in the strip −1 ≤ (t) ≤ 0 and satisfy the boundary condition (2.14), implied by the KMS condition (2.7). Both properties are obscured in the above formula, because it was derived for t ∈ R and the x k depend implicitly on t. Therefore, one would wish for an alternative expression which makes these properties manifest. Fortunately, such an expression exists. In the next section we compute the modular twopoint function directly from the resolvent. Since all the integrals involved are convergent, this yields a final result for the modular correlator that depends on Z(x) − Z(y) rather than on the roots x k . Furthermore, the resulting expression can explicitly be shown to satisfy the boundary condition (2.14).

Plane
This is a rather trivial case, but the result is a useful building block for the more involved calculations in the other cases. Note that all results hold for arbitrary configurations of intervals. We want to calculate the modular two-point function (2.19) where we used the shorthandt = Z(x)−Z(y) again and γ is a tight counter-clockwise contour around the interval [0, 1] as before. The first term vanishes since the integrand is analytic in the integration region. For the second term we find it convenient to substitute z = (1 − λ)/λ and thus we may write with Γ being a tight counter-clockwise contour around R + . We note that z it has a branch cut on R − while (−z) it has a branch cut on R + which makes it impossible to use the residue theorem. However, we can circumvent this problem with a trick that produces a common branch cut on R + for the whole integrand. This goes as follows.
Assuming −1 < (t) < 0, we can pull tight the integration contour and neglect the contributions at the endpoints 0 and ∞. While (t) < 0 avoids the divergence at zero, (t) > −1 ensures the vanishing of contributions at infinity. Hence, we may write Notice that we went from a complex contour integral to an integral along R + and from (−z) it to z it by acquiring a factor [e −πt − e πt ]. Now we use the inverse logic and go back to the contour integral and from z i(t+t) to (−z) i(t+t) by adding a factor [e −π(t+t) − e π(t+t) ] −1 . Thus, we arrive at .
where we made explicit the above assumption that (t) < 0. The remaining contour integral is solved by closing Γ at infinity with a big circle as illustrated in figure 7. We only pick up the residue at z = −1 which is just one. So the final result for the modular two-point function on the plane is This expression was expected, since it is known already, but it has, to the best of our knowledge, not yet been derived directly from the resolvent. As a first consistency check, we see that for t = 0 we retain the original two-point function. It is also straightforward to check that it satisfies the KMS condition. Since eq. (4.8) is analytic for all t ∈ C, we do not have to care about the direction of any limits and can straighforwardly substitute t → t − i. Then, the sinh in the denominator switches sign and thus the entire expression switches sign. This verifies KMS condition in the form of (2.14) since the right hand side vanishes up to the localized contribution (3.13).  .7) is carried out along the contour Γ displayed here. We can close the contour at infinity with a big circle since those contributions vanish. The only non-analyticity in the integration region is the simple pole at z = −1 for which we obtain the residue.

Cylinder
Now we will take the spatial direction to be periodic with periodicity 1. For the boundary conditions there are two possibilities. In the first case of antiperiodic boundary conditions, there is no more work to be done. Recalling the details given in section 2.3 we see that the only difference to the case of the plane lies in the definitions of the propagator (2.28) and the function Z(x) (2.32). However, both of these definitions did not play a role in the calculation. Hence, the modular two-point function on the cylinder with antiperiodic boundary conditions is already given by the expression in (4.8).
The case of periodic boundary conditions is more complicated insofar as the resolvent (2.31) features an additional term. The corresponding contribution to the modular two-point function is where we used the same definitions for z,t, γ and Γ as in the previous subsection. We may split the expression into two terms after which we recognize the first term to be the same integral as in (4.4) whose solution has been given in (4.8). Let us refer to the second term as ∆G mod . To solve it we make use of the identity (3.22). It follows that Again, the contour integral is of the same type as in (4.4), which leaves us with the integral over s. We find it convenient to substitute w = e πs and arrive at . (4.13) Our further solution is restricted to the special case L = 1 /n, n ∈ N \ {1}. Under this restriction we are able to separate the integrand into partial fractions which results in a sum of simple integrals. To do that we continue with the substitution v = w 1 /n , yielding . (4.14) A decomposition into partial fractions gives .
In the third line we made the negative imaginary part of t explicit as a prescription for how to avoid the pole. Although all terms diverge separately, we can trust the divergences to cancel in total, since the original integral is finite. Importantly, note that the boundary term of the first summand cancels exactly with the first term of δG mod in (4.11). Thus, our final result for the modular two-point function on the cylinder with periodic boundary conditions is sinh πt sinh π t +t + δG It is straightforward to check that δG mod vanishes for t = 0. This is as expected, since we should recover the usual propagator. To see what happens to the KMS condition, we again insert t → t − i. The first term clearly switches sign hence does not contribute to the left hand side of eq. (2.14). For the second term, we can to consider the cases of odd and even n separately. In both cases, we end up with sinh n π t +t , (4.18) reproducing exactly the non-local extra term from (3.24). As far as the authors know, this is the first time in the literature that such a non-local term was explicitly derived in the context of modular two-point functions.

Torus
At last, we will turn to the modular two-point function on the torus, i.e. we will deal with periodicity 1 in the spatial direction and with periodicity τ = iβ in the time direction. Following the labelling in section 2.3, the boundary conditions will be denoted by ν = 2, 3.
To begin with, the same recipe as in the previous subsections can be applied. We take the necessary information about the resolvent and the propagator from the third key point in section 2.3. It strikes us that we can (and will) make use of (3.22) once again. Let us demonstrate the procedure for ν = 2 and then comment on the difference for ν = 3. Combining all things mentioned (again with z,t and Γ defined as before) we have which allows a number of simplifications. The contour integral over z on the right might be an old acquaintance by now. It is of the same type as the one in (4.4), which was solved by (4.8). The sum over k on the left can be recognized as a Dirac comb with periodicity 1.
Putting the changes into place results in which is our final result for the modular two-point function on the torus with PA boundary conditions (ν = 2). To see the slight difference for AA boundary conditions (ν = 3), it suffices to compare eqs. (2.38) and (2.39). The sum is now over k ∈ Z + 1/2 which produces an additional factor (−1) n in the final result. For t = 0, the second factor in (4.21) cancels and we can trace back the steps to G(x, y) as expected. To check for the KMS condition, we replace t → t − i and confirm that the expected sign change is produced by the sinh in the denominator of the second factor, up to the isolated poles which yield local contributions on the right hand side of (2.14).

Analytic structure
As a key feature of two-point functions we deem it instructive to analyze the pole structure of the above results. This will give insights into the non-locality of modular flow, since it causes couplings between multiple points or even entire regions.
We illustrate the analytic structure of G mod (x, y, t) for t ∈ C, for some fixed x, y ∈ V . As anticipated in section 2, the presence of poles and cuts gives information about the locality or not of the modular evolution. First, we know that G mod (x, y, t) possesses simple poles along the real axis, located at the solutions to (3.14) and (3.28) for the plane/cylinder, and the torus respectively. Due to the definition of the modular correlator and the KMS condition (2.10), we also know that the function is antiperiodic in imaginary modular time. The precise values of x, y -and therefore the poles -will not be important for the discussion.
Plane or cylinder (A). The modular two-point function in the plane or on the cylinder with antiperiodic boundary conditions was determined in (4.8). Given fixed x, y, the poles are the values of t which satisfy Notice that, since Z(x) ∈ R for x ∈ V (and similarily for y), the solutions of this equation lie at (t) ∈ Z. In figure 8 we illustrate the structure of the modular two-point function when continued to the complex plane. The sketch can be easily understood with the aim of figure  4. As portrayed there, if we fix y and evolve in t, the values of x that satisfy (4.22) move, and sweep the entirety of V when t varies in (−∞, ∞). In other words, for fixed x, y ∈ V there exists a unique value of t ∈ R that solves (4.22). This is a poles of G mod , depicted as a black dot in figure 8, which gets repeated due to antiperiodicity. for fixed x, y belonging to a single interval on the plane or cylinder(A). The function is analytic everywhere except at its simple poles (black dots), their specific location depending on x, y. Since the modular flow in this case is local, G mod (t) possesses no branch cuts. For multiple intervals, the location of the poles is shifted, but no new poles arise. The KMS condition ensures that G mod is antiperiodic in imaginary time.

Cylinder (P).
On the cylinder with periodic boundary conditions we have the same situation as above plus the contributions coming from the additional term δG mod in (4.17). This is the novel ingredient in this result. As explained in (4.18), this contribution induces a discontinuity along (t) ∈ Z, as depicted in figure 9. At first sight, the appearance of a cut might seem surprising, due to its absence in the antiperiodic sector. However, this will be clarified in brief once we discuss the case on the torus, and consider the limit of zero temperature. : Sketch of G mod (x, y, t) (4.16), for fixed x, y belonging to a single interval on the plane or cylinder (P). The function is analytic in the interior of each strip and antiperiodic amongst strips due to the KMS condition. Just as in the plane, it possesses simple poles (black dots) along (t) ∈ Z. However, the novelty here is that in addition it has branch cuts, indicating that the modular flow is completely non-local as explained around (2.14).
As a remark, notice that there is no cut at t ∈ iZ, since by definition the modular two-point function coincides with the propagator at t = 0.
Torus. To analyze the pole structure of the modular two-point function for the torus with PA boundary conditions, we reconsider our result (4.21). The first factor inside the sum, although it seems obscure, represents the usual propagator G(x, y) and is not important for our discussion. The second factor determines the interesting poles as the solutions to Since we already know the effect that a composition of multiple disjoint intervals has, let us focus on the case of a single interval. Even then there are countably infinitely many solutions for t, namely one for each k, spaced by L/β. Hence, the modular two-point function on the torus (PA) shows a non-locality that consists of infinitely many bi-local couplings. for fixed x, y belonging to a single interval on the torus. The function is analytic everywhere except at simple poles (black dots), located at the solutions to (3.28), which lie on a lattice of spacings i and L/β respectively. In the limit of low temperature β → ∞, the real component of the spacing between poles vanishes. This behaviour is the origin of the branch cut in the case of the cylinder (P), figure 9.
In the low-temperature limit β/L → ∞ we note that the poles move closer to each other until they finally form a branch cut, precisely the situation that we found for the cylinder (P). This is consistent with the interpretation of the cylinder as the low-temperature limit of the torus. By the same line of argument, we expect this infinite set of poles not to appear for the torus with AA boundary conditions. Indeed, the factor (−1) n in the sum of (4.21) causes the poles to cancel each other, in a way completely analogous to the modular Hamiltonian [39].
In order to understand how a cut might originate from an infinite set of simple poles, consider the following example: Here, we produced a function with a branch cut in the interval [a, b] by 'summing' over infinitely many simple poles between a and b. This is precisely what happens in the low temperature limit of the torus, where the lattice of simple poles create a branch cut in the modular two-point function. In turn, this implies that for any fixed modular time t, the modular flow on the cylinder couples a given point y to the entirety of the region V .

Discussion
In this paper we computed the modular flow for the chiral fermion CFT in 1+1 dimensions, for entangling regions consisting of arbitrary sets of disjoint intervals. Working in the framework of functional calculus, we derived two important formulae: 1) The modular flow of the field operators (2.18), 2) the associated modular two-point function (2.19). This was done for the cases of the vacuum and thermal states, both on the infinite line and on the circle, giving an extensive overview of cases that are usually of interest [49]. A central element in our analysis is that we made extensive use of the resolvent method. This technique has allowed us to resolve two main obstacles in the understanding of fermionic entanglement. First, we computed modular flow directly, bypassing the need for the modular Hamiltonian. This is important because although the modular Hamiltonian for the cases of interest was known [37][38][39][40], determining the flow from the Hamiltonian is in general very involved and remained unknown. The second problem one faces is that, even with the knowledge of the operator flow, finding the modular two-point function analytically appears hopeless at first sight. The reason is that, as discussed in section 4, the operator flow generically involves solving transcendental equations whose solutions are unknown. However, we have shown that the resolvent yields directly the modular correlator in closed analytic form, without even the need of considering such equations.
On the circle, we considered both antiperiodic (A) and periodic (P) boundary conditions, yielding closed form expressions for both modular flows, including the extra terms that appear due to a zero mode contribution in the periodic sector. These extra terms appear in the modular two-point function as a branch cut -a feature which, to our knowledge, is novel for exact solutions and which can be traced back to the non-locality of the associated modular flow, as we discussed in section 4.4. On the torus, we considered the ν = 2, 3 sectors, showing that the flow leads to an infinite bi-local set of couplings. Moreover, we illustrated how the analytic structure of the two vacua on the cylinder arises as the low temperature limit of the exact result on the torus.
These solutions display a spectrum of degrees of non-locality. To discuss them, let us restrict to the periodic case: Starting at high temperature, the correlations are dominated by thermal fluctuations. Entanglement is thus suppressed and the reduced density matrix is still a thermal state, albeit with a different temperature β/L due to the introduction of the region size L as an additional scale. As a result, modular flow coincides with regular, local, time evolution. While this is to be expected in fairly general theories, it can also be seen explicitly from our results in eq. (3.28),(3.30) and figure 6: In the limit of small β, the Z-terms vanish, and a field initially localized at y is transported to the point x = y − tβ/L. As we lower the temperature, we see that modular flow adds bi-local couplings between an infinite discrete set of points, one for each value of k in eq. (3.28). As we derived in eq. (3.25), (3.26), these additional couplings can be traced back to the anti-periodicty of the thermal propagator, i.e., the KMS condition in physical time. Finally, as the temperature approaches absolute zero, the discrete couplings condense to a continuum, which we derived in eq. (3.24).
As we saw in this paper, we have explicit results for mixed states in regions that consists of arbitrary sets of disjoint intervals. We cover finite temperature states, as well as that of the perfect mixture of degenerate ground states. It is remarkable how far one can actually go: This method provides not only a recipe to determine modular flow 'in principle', butby virtue of the power of complex analysis -allows to determine modular flow in closed analytic form.
Let us point out some future directions. In addition to naturally continuing the work on the modular theory for free fermions on the torus, our results may give a starting point for computations of more general modular flows. For example, the methods presented in this paper can readily be applied to excited states in the fermion CFT, since the corresponding resolvent is known [38]. We expect to find additional non-local terms there.
Another possible direction for further progress lies in the generalization to higher dimensions, and even to massive (i.e. non-conformal) theories: Since the resolvent formulae that we derived are valid for any free fermionic theory, they can be readily applied if one has the corresponding resolvent at hand, which in turn can always be found by solving an integral equation.
Finally, since the algebra of higher spin fields embeds into the free fermion CFT [50], it will be interesting to see if similar methods can be used to derive modular flows there. We expect this to have applications to higher spin AdS/CFT [51][52][53][54][55] and to bulk reconstruction in this context. A starting point in this direction will be to explore the modular flow of higher spin operators by adapting the techniques presented here.