Modular Hamiltonians for the massless Dirac field in the presence of a boundary

We study the modular Hamiltonians of an interval for the massless Dirac fermion on the half-line. The most general boundary conditions ensuring the global energy conservation lead to consider two phases, where either the vector or the axial symmetry is preserved. In these two phases we derive the corresponding modular Hamiltonian in explicit form. Its density involves a bi-local term localised in two points of the interval, one conjugate to the other. The associated modular flows are also established. Depending on the phase, they mix fields with different chirality or charge that follow different modular trajectories. Accordingly, the modular flow preserves either the vector or the axial symmetry. We compute the two-point correlation functions along the modular flow and show that they satisfy the Kubo-Martin-Schwinger condition in both phases. The entanglement entropies are also derived.


Introduction
Entanglement quantifiers are crucial quantities to explore quantum field theories, quantum gravity models, condensed matter systems and quantum information theory.
The geometric entanglement between complementary spatial regions is studied by considering a quantum system in a state described by its density matrix ρ and a spatial bipartition A∪B of the entire space. Assuming that the Hilbert space can be factorised as H = H A ⊗H B , the reduced density matrix ρ A ≡ Tr H B ρ ≡ Tr B ρ of the region A is obtained by tracing out the degrees of freedom corresponding to the complementary region. The reduced density matrix ρ A is both hermitean and positive semidefinite, hence it can be written as ρ A ∝ e −K A , where the proportionality constant guarantees the normalisation condition Tr A ρ A = 1. The hermitean operator K A is the modular Hamiltonian (also known as entanglement Hamiltonian) of the region A [1]. The spectrum of the modular Hamiltonian provides interesting quantities like e.g. the entanglement entropy.
The modular Hamiltonian K A allows to introduce a family of unitary operators U (τ ) = e −iτ K A , parameterised by τ ∈ R, that generates a flow O(τ ) ≡ U (τ ) O U (−τ ) for any operator O localised in A, which is called modular flow of the operator O [2]. The modular flow of O satisfies Tr A ρ A O = Tr A ρ A O(τ ) . The modular flow provides the intrinsic internal dynamics induced by the reduced density matrix of a subsystem.
The modular Hamiltonian is known analytically in terms of the fundamental fields of the model in few cases. A seminal result due to Bisognano and Wichmann [3,4] states that the modular Hamiltonian of half-space x > 0 for a Lorentz invariant quantum field theory in its vacuum is given by the boost generator in the x-direction. Other analytic expressions of modular Hamiltonians have been found for conformal field theories by combining this crucial result with the conformal symmetry. For a conformal field theory in generic spacetime dimension and in its ground state, this analysis provides the modular Hamiltonian of a spherical region [5][6][7]. In 1 + 1 dimensions, where the conformal symmetry has infinitely many generators, the Bisognano-Wichmann modular Hamiltonian and the conformal transformations allow to find the modular Hamiltonians in some other cases of physical interest [8,9], including time-dependent scenarios [9] and spatially inhomogeneous systems [10]. All these modular Hamiltonians are local because the corresponding densities are local operators.
Very few analytic expressions of modular Hamiltonians are available in the literature that cannot be found through the result of Bisognano and Wichmann and the conformal symmetry, even in conformal field theories. An important example is the modular Hamiltonian of the 1 + 1 dimensional massless Dirac fermion in the ground state for the union of disjoint intervals on the infinite line. By employing the modular Hamiltonian on the lattice obtained by Peschel [11], this modular Hamiltonian has been written by Casini and Huerta [12], who found also the corresponding modular flow for the Dirac field. This modular Hamiltonian is very interesting because, beside a local term, it contains also a bi-local operator that induces a mixing along the modular flow between the field evaluated in two different points. The correlators of the Dirac field along the modular flow generated by this modular Hamiltonian have been obtained in [13], verifying that they satisfy the Kubo-Martin-Schwinger (KMS) condition [1]. The validity of this condition guarantees the uniqueness of the modular flow [2]. Further analyses of this modular Hamiltonian and of its modular flow are reported in [14,15].
Physical boundaries heavily influence the entanglement of complementary spatial regions. For instance, in two-dimensional boundary conformal field theories [16][17][18], the Affleck-Ludwig boundary entropy [19] occurs in the entanglement entropy [20]. The role of a physical boundary has been studied also through algebraic quantum field theory methods [21][22][23].
A class of modular Hamiltonians in two-dimensional conformal field theories can be found through boundary conformal field theory techniques [9] and this approach allows to explain various features of some entanglement spectra [9,10,[24][25][26][27][28]. The modular Hamiltonians of more complicated configurations can be found by employing techniques based on the specific model. The typical example is the massless Dirac fermion, where the modular Hamiltonians of an arbitrary number of disjoint intervals on the infinite line when the entire system is in its ground state [12] and of an interval on the circle when the system is in a thermal state [15,29,30] have been found. Other interesting related studies are [31][32][33].
In this manuscript we consider the massless Dirac fermion in its ground state on the halfline x > 0. Imposing boundary conditions at x = 0 that guarantee the energy conservation, two inequivalent models (phases) are allowed. These two phases are characterised by different conservation laws. In particular, either the charge or the helicity is preserved but not both of them. Instead, for the massless Dirac field on the line both these symmetries are preserved. Furthermore, on the half-line the two components of the massless Dirac field are coupled through the boundary condition.
In this model, we study the modular Hamiltonians K A of an interval A on the half-line. Analytic expressions for K A in terms of the components of the Dirac field are obtained in both phases. Beside the expected local term, also a bi-local term occurs which involves fields evaluated in two conjugate points within the interval. This bi-local term breaks either the vector or the axial symmetry. This is the main difference with respect to the bi-local term in the modular Hamiltonian of two disjoint equal intervals on the line for the massless Dirac fermion in the ground state found in [12], which preserves both these symmetries.
The modular flows of the Dirac field generated by these modular Hamiltonians are also obtained. Depending on the phase, the modular flow mixes fields with either different chirality or different charge evaluated in conjugate points. The characteristic symmetry of each phase is preserved along the corresponding modular flow. We find the two-point correlators of the fields along the modular flow that satisfy the KMS condition.
The outline of the manuscript is as follows. In Sec. 2 we discuss the massless Dirac fermion on the half-line and its correlation functions. In Sec. 3 we derive the modular Hamiltonians K A of an interval in the two phases. The corresponding entanglement entropies are computed in Sec 4. The modular flows of the components of the Dirac field and their correlators are discussed in Sec. 5 and in Sec. 6 respectively. Some interesting limiting regimes are considered in Sec. 7. In Sec 8 we study K A and its flow in the spacetime. The results are summarized in Sec. 9. In the Appendices A, B and C we provide some technical details, the derivations of some expressions reported in the main text and further analyses, including some results for the massless Dirac fermion in spacetimes that are invariant under spatial translations.

Dirac fermions on the half-line
In this section we summarise the main properties of massless Dirac fermions on the half-line R + ≡ [0, ∞), fixing also the notation adopted throughout the manuscript.

Bulk dynamics
The massless Dirac field ψ(t, x) is the following doublet made by the two complex fields In the bulk of the half-line, it satisfies the scale invariant equation of motion Denoting the Hermitean conjugation through the asterisk, the energy-momentum tensor for the Dirac field (2.1) can be written in terms of the following components 1 The equation of motion (2.2) implies the local energy conservation of the energy-momentum tensor The bulk dynamics is invariant both under the vector phase transformation ψ r (t, x) −→ e i θv ψ r (t, x) θ v ∈ [0, 2π) (2.7) and under the axial phase transformation ψ r (t, x) −→ e i(−1) r θa ψ r (t, x) θ a ∈ [0, 2π) .
Denoting by j t (t, x) and j x (t, x) the components of the current corresponding to the vector U (1) symmetry and by k t (t, x) and k x (t, x) the components of the current corresponding to the axial U (1) symmetry, we have that These currents are locally conserved when the equation of motion (2.2) holds, namely and they describe the electric and helical transport in the system. 1 The remaining components can be obtained from the symmetry and tracelessness conditions.

Boundary conditions
In order to fully determine the dynamics on the half-line R + , the boundary condition at x = 0 must be imposed. This choice deeply influences the symmetry content of the model.
We adopt the most general boundary conditions that ensure global energy conservation. This means that the vanishing of the energy flow through the boundary is imposed, namely [16][17][18] T xt (t, 0) = 0 t ∈ R . (2.11) This boundary condition can be satisfied in two ways: either which provide a scale invariant coupling of fields with different chirality at x = 0.
We remark that α v and α a parametrise all self-adjoint extensions [34] of the Hamiltonian which is obtained by rewriting the equation of motion (2.2) in the form The boundary conditions (2.12) and (2.13) lead to different conservation laws for the system. In particular, the boundary condition (2.12) preserves the vector symmetry and breaks the axial one, while the opposite holds in the model where (2.13) is imposed. Thus, on the halfline the basic physical requirement of energy conservation (i.e. the self-adjointness of the Hamiltonian) implies that either the vector symmetry or the axial symmetry is necessarily broken. Instead, for the massless Dirac fermion on the infinite line both these symmetries are preserved. At quantum level this feature provides two different phases for the massless Dirac fermion on the half-line, which are characterised by either the conservation of the charge or by the conservation of the helicity (but not both of them) [35]. Throughout this manuscript we will refer to these two inequivalent models respectively as vector phase and axial phase.

Quantization
The quantum fields ψ r (t, x) in (2.1) satisfy the equation of motion (2.2), the following the equal-time anti-commutation relations and also a boundary condition that is either (2.12) or (2.13). These fields can be described in terms of two mutually anticommuting algebras A + and B + generated respectively by which satisfy the canonical anticommutation relations.

Vector phase
In the vector phase, where the boundary condition (2.12) is imposed. We denote the components of the Dirac field by can which be written as [36] Notice that each component λ r depends on the angle α v and involves the same generators of A + and B + , which is not the case of the Dirac fermion on the line 2 .
For the vacuum expectation values of the fields in (2.20) and (2.21) in the Fock representation of A + and B + , one finds where we have introduced and In this phase the axial symmetry (2.8) is broken, hence the nontrivial correlation functions (2.26) and (2.27) that mix different helicities are allowed.
We find it convenient to collect the correlation functions (2.24), (2.25), (2.26) and (2.27) at equal times t 1 = t 2 ≡ t into the following matrix On the infinite line R, two copies of A+ and B+ are needed to write the components of the Dirac field as It is worth discussing the meaning of the phase factor e iαv in (2.12) and the possibility to absorb it through the field redefinition which leaves invariant both the equations of motion (2.2) and the equal-time canonical relations (2.16) and (2.17). This leads us to describe the boundary condition at x → ∞.
The boundary condition at x → ∞ can be studied by introducing the following family of states where Ω is the Fock vacuum and D(R + ) is the space of C ∞ functions with compact support in R + . It is well known that the Fourier transform is a smooth function that is rapidly decreasing as |y| → ∞. We take h ∈ D(R + ) such that Within the family of one-particle states (2.32) defined by (2.32) and (2.34), one can evaluate the following expectation values However, considering the limit along the components of the light cone, given by t = −x and t = x, one finds lim x→∞ x=−t Comparing (2.12) and (2.37), we conclude that the phase e iαv cannot be absorbed in λ 2 as in (2.31) at x = 0 and at x = ∞ simultaneously.
Finally we observe that the angle α v has a simple physical interpretation in the context of scattering theory: the boundary induces a non-trivial one-body scattering matrix, which describes the particle reflection at x = 0. The angle α v defines the scattering phase shift of an outgoing particle with respect to an incoming one.

Axial phase
The axial phase implements the boundary condition (2.13). Denoting by χ r (t, x) the components of the Dirac field in this phase throughout the manuscript, we have that these fields can be decomposed through the generators of A + and B + introduced in Sec. 2.3.1. These decompositions read which depend on the angle α a . The corresponding correlation functions on the vacuum are Thus, the matrices (2.30) and (2.44) coincide, except for the angles α v and α a .
The discussion made in Sec. 2.3.1 about the meaning of the phase e iαa in (2.13) and the possibility to reabsorb it through a redefinition of the fields (2.45) can be easily adapted also to this phase, arriving at the same conclusion.
Since in the vector phase both components of the doublet λ(t, x) defined in (2.1) transform in the same way under the vector phase transformations (2.7), in the axial phase we find it convenient to introduce the following doublet whose components transform in the same way under the axial phase transformations (2.8).
At this point, in order to treat both phases in a unified way, we introduce the doublet with λ and χ given by (2.19) and (2.45) respectively. With this notation the boundary conditions (2.12) and (2.13) take the form

Modular Hamiltonians of an interval on the half-line
The massless Dirac fermion is described by a quadratic field theory, hence the modular Hamiltonian of an interval A = [a, b] ⊂ R + on the half-line for a massless Dirac fermion in its ground state can be written in the following quadratic form where : · · · : denotes the normal product in the oscillator algebras A + and B + and ψ(t, x) is the two-components field in (2.46).
The kernel H A (x, y) in (3.1) is the 2 × 2 matrix given by [11,37,38] where I is the identity matrix and C A is the reduced correlation functions matrix, obtained by restricting the correlation functions matrix (given by (2.30) in the vector phase and by (2.44) in the axial phase) to the interval A.
In order to obtain an explicit expression for H A , the spectral problem associated to the reduced correlation functions matrix C A must be solved. This means that we have to find the eigenvalues σ s and the eigenfunctions Φ s, where s and p are two parameters specified below.

The spectral problem
In order to solve the spectral problem (3.3) for the massless Dirac field on the half-line when the subsystem is the interval A = [a, b] ⊂ R + , let us first consider the auxiliary spectral problem corresponding to the massless Dirac field in its ground state on the infinite line and where the bipartition of the line is given by the two disjoint equal intervals A sym ≡ [−b, −a] ∪ [a, b] ⊂ R and by its complement on the line. This auxiliary spectral problem reads where C(x − y) is the distribution defined in (2.29). The solution of this spectral problem can be found by specialising to A sym the solution of the spectral problem corresponding to an arbitrary number of disjoint intervals of generic lengths on the line, which has been found in [12] by employing [39]. The explicit form of the spectral data {σ s , φ s,p | s ∈ R, p = 1, 2} for (3.4) has been reported in the Appendix A (see (A.1) and (A.2)).
The spectral problem (3.3) can be solved through the auxiliary spectral problem (3.4) by first observing that the latter one can be rewritten as follows b a C(x − y) φ s,p (y) dy + b a C(x + y) φ s,p (−y) dy = σ s φ s,p (x) x ∈ A sym .
Since A sym is symmetric under reflection with respect to the origin at x = 0, we can conclude that (3.5) holds also for x → −x. This observation provides another identity that can be combined with (3.5). This leads to write the solution of (3.3) in terms of φ s,p and α as follows We stress that the function φ s,p in (3.6) is the solution of the spectral problem (3.4) in A sym . The completeness and orthonormality of the system {φ s,p (x) : x ∈ A sym , s ∈ R, p = 1, 2} implies that (3.6) form a complete set of orthonormal eigenfunctions in A, namely where the overline indicates the complex conjugation.

The kernel H A
The solution of the spectral problem (3.3) leads us to write the reduced correlation functions matrix through its spectral representation The relation (3.2) tells us that C A and H A share the same eigenfunctions and that, since the eigenvalues σ s of C A are (A.1), the eigenvalues of H A are given by −2πs = log(1/σ s − 1) with s ∈ R. Thus, the spectral representation of the kernel H A reads Expressing the eigenfunctions (3.6) through the explicit form of φ s,p given in (A.2), we find that plays an important role throughout our analysis, and m(x, y) is defined in terms of the functions m p (x) in (A.3) as follows The integration in (3.9) can be performed by observing that where η x , η y ∈ {−1, +1}. The support of the Dirac delta function in (3.13) is given by the zeros of the following function (3.14) A crucial role is played by the pointx conjugate to x defined as The integration in (3.9) of the diagonal elements of the matrix (3.10) provides Dirac delta functions localised where z w (x, y) = 0 and z w (−x, −y) = 0. These equations are solved by y = x and y = −x, but only the former solution is allowed because the latter one does not belong to R + when x ∈ A. Thus, for these integrals we have to use (3.16) which tells us that the diagonal elements in (3.10) lead to local terms in the modular Hamiltonian. Instead, the off-diagonal elements of (3.10) give Dirac delta functions localised on the zeros of z w (x, −y) and z w (−x, y). They are y =x and y = −x, but only the former solution is allowed because −x / ∈ R + when x ∈ A. Thus, for the integration in (3.9) of the off-diagonal elements of (3.10) we need which tells us that these terms give origin to bi-local terms in the modular Hamiltonian.
The above discussion suggests to separate the diagonal terms and the off-diagonal terms in (3.10). This leads to write the kernel (3.9) as follows where H loc A and H bi-loc A give origin respectively to the local terms and the bi-local terms in the modular Hamiltonian. More explicitly, the matrices in the r.h.s. of (3.18) are and where we have introduced and that satisfies the relation D(−x, y) = −D(x, −y).
The decomposition (3.18) leads to write the modular Hamiltonian (3.1) as follows where K loc A is a local operator, while K bi-loc A is a bi-local operator, which are discussed in Sec. 3.2.1 and Sec. 3.2.2 respectively.
Setting either α = α v or α = α a in the above expressions, we obtain the corresponding results respectively for the vector phase and for the axial phase.

Local term
The local term in (3.23) reads where the kernel H loc A (x, y) has been defined in (3.19). The Dirac delta functions occurring in H loc A (x, y) guarantee that (3.24) is a local operator because it can be written as an integral over A of fields evaluated at the same point. This integral can be found by first plugging (3.19) into (3.24) and then integrating by parts. Exploiting the fact that [∂ x M + (x, y) − ∂ y M + (x, y)] | x=y = 0, we find that (3.24) becomes where the operator T tt (t, x) is the normal ordered version of the energy density (2.4), namely (3.26) and the weight function reads (3.28)

Bi-local term
The bi-local term in the decomposition (3.23) of the modular Hamiltonian (3.1) reads A (x, y) ψ(0, y) : dx dy (3.29) where the kernel H bi-loc A (x, y) is defined in (3.20).
The operator (3.29) is bi-local because the Dirac delta functions in H bi-loc A (x, y) allow to write it as an integral over A of fields evaluated at different (conjugated) points. This integral can be computed by plugging (3.20) into (3.29) first and then integrating by parts. Since M − (x, ab/x) = 0, we find that (3.29) becomes where we have introduced the following bi-local operator and the weight function is which satisfies

Modular Hamiltonians in the vector and axial phases
The explicit expressions of the modular Hamiltonians (3.23) for the interval A ⊂ R + on the half-line when the Dirac field is in the ground state can be written by first using (3.25) and (3.30) and then performing the substitutions (ψ, α) → (λ, α v ) for the vector phase and (ψ, α) → (χ, α a ) for the axial phase into the operators (3.26) and (3.31), where λ and χ have been defined in (2.19) and in (2.45) respectively.
In the local term (3.25), taking into account that the fermion fields anticommute under the normal product, we observe that T tt in (3.26) has the same form when expressed in terms of λ or χ. This is a consequence of the fact that T tt is invariant under both the vector and the axial transformations given in (2.7) and (2.8) respectively.
In the bi-local term (3.30), the bi-local operator (3.31) in the vector phase reads while in the axial phase it becomes These expressions depend explicitly on the angle determining the boundary condition at x = 0.
It is instructive to compare these modular Hamiltonians to the modular Hamiltonian of two disjoint equal intervals A sym on the line, obtained as a special case of the modular Hamiltonian of the union of a generic number of disjoint intervals on the line found in [12]. Also the modular Hamiltonian of A sym ⊂ R can be written as the sum of a local term and a bi-local one (3.37) In these expressions T tt (t, x) is the normal ordered energy density (3.26), the weight functions β loc (x) and β bi-loc (x) are (3.27) and (3.32) respectively and the bi-local operator T bi-loc (t, x, y) occurring in the bi-local term is defined as follows where we remark that the fields in this operator are given by (2.22) and (2.23).
In the local term (3.25), the weight function and the functional dependence on the fermion fields in the integrand are the same occurring in the local term K loc Asym of the modular Hamiltonian of A sym on the line, that is the first expression in (3.37). However, we stress that the fermionic fields on the half-line are given by (2.20), (2.21), (2.38) and (2.39), which depend explicitly on the angles α v and α a characterising the boundary condition at x = 0.
As for the bi-local terms, comparing (3.30) and the second expression in (3.37), one observes that, while the same weight function occurs in the integrands, the corresponding bi-local operators are very different. Indeed, we have that (a) the operator (3.38) is invariant under both vector (2.7) and axial (2.8) transformations, while this is not the case for (3.34) and (3.35) which preserve separately only the vector and axial symmetry; (b) the point conjugate to x in the integrand of K bi-loc Asym is −x, which belongs to the opposite interval in A sym with respect to x, while in the integrand of (3.30) both x and its conjugate pointx belong to the interval A; hence also the self-conjugate point x = √ a b (where x =x) occurs; (c) the bi-local operators (3.34) and (3.35) depend explicitly on the boundary conditions through the angular parameter characterising the corresponding phase.

Entanglement entropies
The Rényi entropies are defined through the moments of the reduced density matrix Tr A ρ n A for integers n 2 as They provide the entanglement entropy S A by means of the following analytic continuation of the integer parameter n In a two dimensional conformal field theory, the moments of the reduced density matrix have been computed as the correlation functions of the branch point twist fields T n [20].
When the interval A = [0, b] is adjacent to the boundary of the half-line and the entire system is in the ground state, Trρ n A = T n (b) ∝ (2b/ ) −∆n is the one-point function of the twist field, where is the ultraviolet cutoff and ∆ n is the conformal dimension of the twist field given by which is proportional to the central charge c of the model.
When the interval A = [a, b] on the half-line is not adjacent to the boundary, Trρ n A can be found as the two-point function of the twist fields. Combining the Rényi entropies of two disjoint equal intervals on the line [40][41][42][43] with the method of the images, one obtains where r is the cross ratio of the endpoints of the interval and of their images, while c n is a constant. The function F n (r) depends on the full operator content of the boundary conformal field theory, hence it encodes also the conformal boundary state, in a highly non-trivial way. These one-point and two-point correlation functions of twist fields allow to construct the following ultraviolet finite combinations .

(4.5)
A relevant quantity encoding some properties of the boundary is the Affleck-Ludwig boundary entropy log(g) [19]. Considering a conformal field theory in its ground state either on the half-line or on the line, the Affleck-Ludwig boundary entropy can be found by combining the Rényi entropies S The massless Dirac field is a conformal field theory with c = 1. In this model, the Rényi entropies and the entanglement entropy can be evaluated also through the reduced correlation functions matrix respectively as [37,38] S (n) and where I is the 2 × 2 identity matrix, γ is a generic 2 × 2 matrix and g(γ) ≡ −∂ n g n (γ)| n=1 .
By employing the spectral representation of C A in (4.8), for the entanglement entropies S

(n)
A with integer n 1, which are given by S (1) A ≡ S A and (4.7) for n 2, we find where σ s is (A.1), g n (σ) is defined in (4.8) for n = 1 and in (4.7) for n 2, and A ≡ (a+ , b− ) with → 0 + has been introduced to regularise the integral. By using (A.2) and (3.6), we find where the explicit expressions of m p (x) are given in (A.3). Notice that (4.10) is independent of both the parameters α and s. The independence of s leads to the factorisation of the two integrals in (4.9); hence we can evaluate them separately. By using +∞ −∞ g n (σ s ) ds = π(n+1) 12 n ; for the entanglement entropies we obtain (4.11) We remark that this expression is independent of α.
Plugging (4.4) into (4.1) and comparing the resulting expression with (4.11), we find that F n (r) = 1 (4.12) identically in both the phases and for any choice of the boundary condition parameter. A similar simplification has been observed also in [12] for the massless Dirac fermion in its ground state on the line when the subsystem is the union of disjoint intervals, and further explored in [46][47][48]. Lattice results in the XX chain with open boundary conditions for the bipartition that we are considering have been discussed in [49].
The entanglement entropies for an interval A = [0, b] adjacent to the boundary are given by (4.9) From (4.11) and (4.13), we find that the UV finite combination (4.5) in this model reads (4.14) In order to evaluate the Affleck-Ludwig boundary entropy for the massless Dirac field through (4.6), we have to consider the entanglement entropies of an interval A = [a, b] of length = b − a on the line, when the entire system is in its ground state. For the massless Dirac field, the underlying spectral problem is solved by the eigenvalues (A.1) and by the following eigenfunctions [12,37] The corresponding entanglement entropies are and the factor 2 occurs because the two components of the Dirac field give the same contribution. By using (4.13) and (4.16), we can specialise (4.6) to our case, finding g = 1, which has been obtained also in [19,50] through different methods.

Modular flows of the Dirac field
The reduced density matrix ρ A ∝ e −K A generates the one-parameter family of unitary operators given by {ρ iτ A : τ ∈ R} that defines an automorphism on the operator algebra known as modular flow [1].
The modular flow of the two components of the Dirac field are defined as where the initial configuration ψ r (x) in the right hand side is defined by (2.46) at t = 0. The fields (5.1) can be found by solving the following differential equations The modular Hamiltonian K A can be decomposed as in (3.23); hence the r.h.s. of (5.2) is the sum of two terms determined by the local operator K loc A in (3.25) and by the bi-local operator K bi-loc A in (3.30). Both these contributions can be evaluated by means of the equaltime anticommutators (2.16) and (2.17). As for the term provided by K loc A in (3.25), by using the expression of T tt in (3.26) one obtains where the differential operator B loc is defined in terms of the weight function (3.27) as follows The contribution of the bi-local operator K bi-loc A in (3.30) to the r.h.s. of (5.2) can be found by exploiting the expression (3.31). The result reads where the weight function (3.32) occurs andỹ is the point conjugate to y.
The two partial differential equations in (5.2) are coupled because both the components of the Dirac field occur in (5.5). In order to study the resulting system of partial differential equations, we find it convenient to introduce the following doublet of fields By using (5.3) and (5.5) with y = x and y =x, combined with (5.2), in terms of the doublet (5.6) we obtain the following system of four partial differential equations d dτ where the 4 × 4 block diagonal matrix within the square brackets is a differential operator expressed through the 2 × 2 matrix differential operator defined as in terms of the differential operator B loc (x) introduced in (5.4). The solution of the system of partial differential equations in (5.7) can be found by solving with the initial condition This system can be solved by first performing a suitable transformation that decouples the two equations in (5.9). Each of the two resulting decoupled equations defines the flow of a one-parameter abelian group that can be determined through a standard technique. All the details of this procedure are reported in Appendix B.
x < l a t e x i t s h a 1 _ b a s e 6 4 = "       The solution of (5.9) and (5.10) is one of the main results of this manuscript. It reads and ξ(τ, x) can be defined in terms of the function w(x) introduced in (3.11) as follows (5.13) This function describes the modular evolution (parameterised by τ ∈ R) of any point x ∈ A. It satisfies ≡ξ(τ, x) . (5.14) Notice that, since x = √ ab is self-conjugate under (3.15), the second expression in (5.14) implies that ξ(τ, √ ab ) ξ(−τ, √ ab ) = ab for any τ ∈ R.
In Fig. 1, we show the arguments of the fields in the r.h.s.'s of (5.11) in the half-plane parameterised by (ξ, τ ), by employing (5.13). We consider a point x ∈ (a, b) at τ = 0, whose conjugate point isx. The blue curves represent the arguments ξ(τ, x) and ab/ξ(τ, x) occurring in the first equation of (5.11), while the red curves correspond to ξ(−τ, x) and ab/ξ(−τ, x), which appear in the second equation of (5.11). At τ = 0, the solid and the dashed curves pass through x andx respectively. A solid curve and the corresponding dashed one having the same colour intersect at one of the points whose coordinates are (ξ, τ ) = ( √ ab , ±τ 0 ), depending on the position of x with respect to √ ab, where τ 0 is obtained by solving ξ(τ 0 , x) = √ ab. The value of τ 0 can be determined from (B.7), (C.13) and (3.11), We find it instructive to write the explicit expressions of the solution (5.11) in the two inequivalent phases by using (2.46). In the vector phase, the modular flow of the massless Dirac field is . (5.15) This modular flow mixes fields with different chiralities, but the electric charge is preserved. The mixing is non-local because it involves a field in ξ and another one in its conjugate point ab/ξ.
In the axial phase, the modular flow reads . (5.16) This bi-local modular flow mixes fields with different electric charge, but the chirality is preserved, contrary to (5.15).

Correlation functions along the modular flow
The correlation functions of the fermion fields λ r (τ, x) and χ r (τ, x), evolving trough the modular flow generated by the modular Hamiltonian K A , describe the quantum fluctuations along the modular "time" τ . They can be obtained by first employing either (5.15) or (5.16) and then by writing the initial data in the Fock representation of the algebras A + and B + , i.e. using either (2.20) and (2.21) or (2.38) and (2.39) for the initial field configurations.
The derivation significantly simplifies by adopting the identity w(ξ(τ, x)) = 2π τ + w(x) written in the exponential form instead of the cumbersome explicit expression (5.13).
We find it convenient to introduce the following distribution where the first factor in the r.h.s. is regular at x = y; indeed, from (3.11) we have .
We remark that (6.1) satisfies By using (3.11), we find that for the interval A ⊂ R + the distribution (6.1) becomes The distribution (6.1) represents the modular counterpart of (2.29). Indeed, analogously to the conventional time evolution, where all the correlators (2.24), (2.25), (2.26) and (2.27) are written in terms of the distribution (2.29), we find that all the non-vanishing correlators of the fields along the modular flow can be expressed in terms of (6.1).
After some algebra, we find that the eight non-vanishing two-point functions in the vector phase can be written as follows where we have adopted the notation The vacuum expectation values (6.5), (6.6), (6.7) and (6.8) satisfy some basic properties. First, the Hilbert space structure of the theory provides the following relations Second, they obey the modular equations of motion following from (5.9) when x 1 = x 2 . For instance, we have that By exploiting (6.5) and (6.7), this equation is equivalent to the following one satisfied by (6.1) (in the limit ε → 0) Notice that the partial differential equation (6.12) involves both x andx, keeping trace of the bi-local character of the modular evolution. We remark that equations similar to (6.12) hold for all correlation functions in (6.5), (6.6), (6.7) and (6.8). Their validity provides a valuable consistency check of the whole construction.
A fundamental property satisfied by the correlation functions (6.5), (6.6), (6.7) and (6.8) is the Kubo-Martin-Schwinger (KMS) condition [1] where r 1 , r 2 ∈ {1, 2}, whose validity is a consequence of the first equality in (6.3). From these KMS relations, one infers that the expectation values (6.5), (6.6), (6.7) and (6.8) behave like thermal correlators with inverse temperature β = 1 in our units. We remark that the KMS condition is a distinguishing property of the modular group (see Theorem 1.2 in chapter VIII of [2]); hence their validity provides a strong check of our results.
Notice that the two properties in (6.3) hold for a distribution W(τ ; x, y) of the form where w(x) and g(x, y) are real functions and g(y, x) = g(x, y). This class of distributions includes the ones occurring in the modular correlators of the massless Dirac field for various interesting cases like the bipartition of the infinite line given by a generic number of disjoint intervals when the entire system is in its ground state [13,15] and the ones discussed in the Appendix B.2.1.
The above analysis can be adapted to write the corresponding quantities for the axial phase. In the axial phase, the correlators of the fields along the modular flow read The crucial difference with respect to the vector phase is that the mixed correlation functions (6.19) and (6.20) are not invariant under the vector phase transformation (2.7).

Special bipartitions
In this section we consider some limiting regimes for the position of the interval A = [a, b] on the half-line. In all these limits, the modular flow of the Dirac field becomes local.

Interval at large distance from the boundary
The first case that we find worth considering is an interval of length at large distance from the boundary. This limit can be performed by first setting b = a + , x = a + v with v ∈ [0, ] and then sending a → ∞. In this limit, for the conjugate point we havex = a +ṽ + O(1/a), whereṽ ≡ − v, and the function (3.11) simplifies to The weight functions (3.27) and (3.32) in the modular Hamiltonian become respectively where one recognises that β 0 (v) is the weight function occurring in the modular Hamiltonian of an interval of length in the infinite line (see (B.22)) [5,7].
In order to get the modular flow of the Dirac field, one first observes that, in this limit, (5.13) becomes This allows to write the modular evolutions of the fields in this regime by specifying (B.17) to this case. The result reads where we used that β 0 (ζ) β 0 (v) = ∂ v ζ. The two-point functions of these fields is obtained by first observing that, in this limiting regime, (6.4) becomes (we also introduce y = a + z) and then employing this result in the correlators of Sec. 6 expressed through this function.

Interval adjacent to the boundary
The case of an interval A = [0, b] at the beginning of the half-line can be studied by taking the limit a → 0 in the above results.
In this limit, the function (3.11) becomes and for the weight functions (3.27) and (3.32) we find Thus, in this limit the bi-local operator in the modular Hamiltonian (3.23) vanishes and the remaining local term becomes the modular Hamiltonian found in [9] for an interval adjacent to the boundary of the half line, specialised to the model that we are considering.
The modular flow of the Dirac field when the interval is adjacent to the boundary can be found by specifying the analysis of the Appendix B.2 to this case. In particular, β 0 (x) is the weight function obtained in (7.7). By using (B.7), the corresponding ξ(τ, x) can be written in terms of w(x) in (7.6), where we remark that x ∈ (0, b), hence the inverse function w −1 (y) is defined when y 0. This result reads Notice that the limit a → 0 of (5.13) gives (7.8) when τ τ 0 and ξ = 0 when τ τ 0 .
Alternatively, we can take the limit a → 0 in the modular flow (5.11), observing that the mixing terms vanish. The modular flow of the Dirac field can be written through β 0 (x) in (7.7) as follows .
The correlators of these fields can be found by first taking a → 0 in (6.4), that gives and then employing this result into the expressions in Sec. 6 written in terms of this function.

Semi-infinite line separated from the boundary
Another interesting bipartition of the half-line to consider is the one where the subsystem A is the semi-infinite line separated from the boundary of the half-line at x = 0, namely

T q M A U W d i H A z g C D 8 6 g A l d Q h R o Q E P A A T / D s K O
f R e X F e p 6 0 Z Z z Z T g D 9 w 3 n 4 A e g m T + g = = < / l a t e x i t > ⌧ < l a t e x i t s h a 1 _ b a s e 6 4 = " j 0 i Y l x h 4 I 7 d t j S R 7 a m w 4 L k 6 s n l Y = " > A A A B 6 3 i c b Z D L S g M x F I Y z 9 V b r r e p S k G A R X J U Z E d S V h W 5 c t m I v 0 A 4 l k 2 b a 0 C Q z J G e E U r p 0 6 8 a F I m 5 9 B x / B t T u f w Z c w 0 3 a h r T 8 E P v 7 / H H L O C W L B D b j u l 5 N Z W l 5 Z X c u u 5 z Y 2 t 7 Z 3 8 r t 7 d R M l m r I a j U S k m w E x T H D F a s B B s G a s G Z G B Y I 1 g U E 7 z x h 3 T h k f q F o Y x 8 y X p K R 5 y S i C 1 2 k C S T r 7 g F t 2 J 8 C J 4 M y h c f V S / 7 w / f b y q d / G e 7 G 9 F E M g V U E G N a n h u D P y I a O B V s n G s n h s W E D k i P t S w q I p n x R 5 N Z x / j Y O l 0 c R t o + B X j i / u 4 Y E W n M U A a 2 U h L o m / k s N f / L W g m E F / 6 I q z g B p u j 0 o z A R G C K c L o 6 7 X D M K Y m i B U M 3 t r J j 2 i S Y U 7 H l y 9 g j e / M q L U D 8 t e m f F y 6 p b K J X R V F l 0 g I 7 Q C f L Q O S q h a 1 R B N U R R H z 2 g J / T s S O f R e X F e p 6 U Z Z 9 a z j / 7 I e f s B c U K S Q Q = = < / l a t e x i t >  F a s B B s G a s G Z G B Y I 1 g U E 7 z x h 3 T h k f q F o Y x 8 y X p K R 5 y S i C 1 2 k C S T r 7 g F t 2 J 8 C J 4 M y h c f V S / 7 w / f b y q d / G e 7 G 9 F E M g V U E G N a n h u D P y I a O B V s n G s n h s W E D k i P t S w q I p n x R 5 N Z x / j Y O l 0 c R t o + B X j i / u 4 Y E W n M U A a 2 U h L o m / k s N f / L W g m E F / 6 I q z g B p u j 0 o z A R G C K c L o 6 7 X D M K Y m i B U M 3 t r J j 2 i S Y U 7 H l y 9 g j e / M q L U D 8 t e m f F y 6 p b K J X R V F l 0 g I 7 Q C f L Q O S q h a 1 R B N U R R H z 2 g J / T s S O f R e X F e p 6 U Z Z 9 a z j / 7 I e f s B c U K S Q Q = = < / l a t e x i t > In this limiting regime, the function (3.11) becomes w(x) = log x − a x + a x > a (7.11) and for the weight functions (3.27) and (3.32), occurring in the local and in the bi-local term of the modular Hamiltonian (3.23), we find respectively (see also [51]) We remark that, despite the fact that β bi-loc (x) is non-vanishing in this regime, the modular Hamiltonian becomes local becausex → ∞ in this limit and the fields {ψ r (x) : r = 1, 2} occurring in the bi-local term (3.30) vanish (see (2.36)).
The arguments of the fields provided by the r.h.s. of (B.17) give the curves shown in the right panel of Fig. 2. The blue and the red curves are ξ(τ, x) and ξ(−τ, x) respectively and they intersect at (ξ, τ ) = (x, 0). Furthermore, they diverge as τ → τ 0 or τ → −τ 0 respectively, where 2πτ 0 = |w(x)|, with w(x) given by (7.11). The occurrence of τ 0 is due to the presence of the boundary. Indeed, for the semi-infinite line in the line [3,4], one finds that ξ(τ, x) = x e 2πτ , that diverges as τ → +∞ (see also the corresponding comment in the Appendix B.2.1). Thus, the internal dynamics of a semi-infinite line depends also on its complement.
The correlators of the fields along this modular flows can be written by employing the results discussed in Sec. 6 and observing that the function W (τ ; x, y) in this regime can be obtained by plugging (7.11) into (6.1), that gives Combining the modular Hamiltonian of this limiting regime and the one obtained in the Sec. 7.2 for the interval adjacent to the boundary, we can write the full modular Hamiltonian K A∪B [1] associated to the bipartition A ∪ B of the half-line where A = [0, ] is an interval of length adjacent to the boundary and B = [ , +∞) is its complement. It reads where 1 A and 1 B denote the identity operators on A and B respectively, while are the modular Hamiltonians respectively of the interval A = [0, ] adjacent to the boundary and of its complement, that is a semi-infinite line at distance from the boundary.
Let us remark that, considering the system on the line, in its ground state and the bipartition A ∪ B of the line where the finite subsystem is the interval A = (− , ), the full modular Hamiltonian is (7.15) where K A and K B are very similar to the operators reported in (7.16) 3 . Indeed, only the integration domains are different, which are A = (− , ) for K A and B = (−∞, − ) ∪ ( , +∞) for K B , as expected.

Modular evolution in the spacetime
The modular evolution whose initial field configurations (at τ = 0) are either (2.20) and (2.21) in the vector phase or (2.38) and (2.39) in the axial phase at t = 0 has been described in Sec. 5. It is worth investigating the most general modular evolution where also the physical time t is involved.
For the massless Dirac field, the physical time can be included by exploiting the fact that each component ψ r of the Dirac field depends only on one of the light-cone coordinates In these coordinates and on the half line, ψ r satisfy the following anti-commutation relations Here we consider two intervals u ± ∈ [a , b], which parameterise the grey diamond D shown in Fig. 3. By employing the spectral problem discussed in Sec. 3.1 in the form and repeating the analysis described in Sec. 3, we obtain the modular Hamiltonian K A in the coordinates u ± . It can be decomposed again as where the local term is while the bi-local term can be written by introducing theũ ± ≡ ab/u ± conjugate to u ± as follows which contains the bi-local hermitean operator defined as follows The weight functions β loc (z) and β bi-loc (z) are defined in (3.27) and (3.32) respectively. Notice that, restricting to the slice defined by constant t = 0, one finds that (8.6) and (8.8) become (3.25) and (3.30) respectively.
The modular flow is defined by where u ± ∈ [a, b] and, according to (2.46), the initial configurations ψ 1 (u + ) and

Z F B f c = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 6 G A y C I I R d Q e J J A l 4 8 R j Q m s F n C 7 K Q 3 G T I 7 u 8 z M C i H k E 7 y I K M G r X + T N v 3 H y O G h i Q U N R 1 U 1 3 V 5 g K r o 3 r f j u 5 l d W 1 9 Y 3 8 Z m F r e 2 d 3 r 7 h / 8 K i T T D G s s 0 Q k q h l S j Y J L r B t u B D Z T h T Q O B T b C / s 3 E b z y h 0 j y R D 2 a Q Y h D T r u Q R Z 9 R Y 6 T 5 r n 7 e L J b f s T k G W i T c n p e r x + N U X F V p r F 7 9 a n Y R l M U r D B N X a 9 9 z U B E O q D G c C R 4 V W p j G l r E + 7 6 F s q a Y w 6 G E 5 P H Z F T q 3 R I l C h b 0 p C p + n t i S G O t B 3 F o O 2 N q e n r R m 4 j / e X 5 m o q t g y G W a G Z R s t i j K B D E J m f x N O l w h M 2 J g C W W K 2 1 s J 6 1 F F m b H p F G w I 3 u L L y + T x o u x d l t 0 7 m 8 Y 1 z J C H I z i B M / C g A l W 4 h R r U g U E X n u E N 3 h 3 h v D h j 5 2 P W m n P m M 4 f w B 8 7 n D + X u k H s = < / l a t e x i t > u + < l a t e x i t s h a 1 _ b a s e 6 4 = " Y F O l J x P Q s m d D W s j N D a c 8 u U P H E G s = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 6 G A y C F 8 O u I P E k A S 8 e I x o T 2 C x h d t K b D J m d X W Z m h R D y C V 5 E l O D V L / L m 3 z h 5 H D S x o K G o 6 q a 7 K 0 w F 1 8 Z 1 v 5 3 c y u r a + k Z + s 7 C 1 v b O 7 V 9 w / e N R J p h j W W S I S 1 Q y p R s E l 1 g 0 3 A p u p Q h q H A h t h / 2 b i N 5 5 Q a Z 7 I B z N I M Y h p V / K I M 2 q s d J + 1 z 9 v F k l t 2 p y D L x J u T U v V 4 / O q L C q 2 1 i 1 + t T s K y G K V h g m r t e 2 5 q g i F V h j O B o 0 I r 0 5 h S 1 q d d 9 C 2 V N E Y d D K e n j s i p V T o k S p Q t a c h U / T 0 x p L H W g z i 0 n T E 1 P b 3 o T c T / P D 8 z 0 V U w 5 D L N D E o 2 W x R l g p i E T P 4 m H a 6 Q G T G w h D L F 7 a 2 E 9 a i i z N h 0 C j Y E b / H l Z f J 4 U f Y u y + 6 d T e M a Z s j D E Z z A G X h Q g S r c Q g 3 q w K A L z / A G 7 4 5 w X p y x 8 z F r z T n z m U P 4 A + f z B + j 2 k H 0 = < / l a t e x i t >
u < l a t e x i t s h a 1 _ b a s e 6 4 = " M B k U 1 g W j j N S / X S o u g H V x 8 b e 3 E + A = " > A A A B 6 H i c b Z D J S g N B E I Z 7 4 h b j F p e b l 8 Y g e A o z g u j N g A c 9 J m A W S E L o 6 d Q k b X p 6 h u 4 a I Q 5 5 A i 8 e F P H q A / g U P o E 3 j 7 6 J n e W g 0 R 8 a P v 6 / i q 4 q P 5 b C o O t + O p m F x a X l l e x q b m 1 9 Y 3 M r v 7 1 T M 1 G i O V R 5 J C P d 8 J k B K R R U U a C E R q y B h b 6 E u j + 4 G O f 1 W 9 B G R O o a h z G 0 Q 9 Z T I h C c o b U q 2 M k X 3 K I 7 E f 0 L 3 g w K 5 + 9 3 X 5 d v e 2 m 5 k / 9 o d S O e h K C Q S 2 Z M 0 3 N j b K d M o + A S R r l W Y i B m f M B 6 0 L S o W A i m n U 4 G H d F D 6 3 R p E G n 7 F N K J + 7 M j Z a E x w 9 C 3 l S H D v p n P x u Z / W T P B 4 K y d C h U n C I p P P w o S S T G i 4 6 1 p V 2 j g K I c W G N f C z k p 5 n 2 n G 0 d 4 m Z 4 / g z a / 8 F 2 r H R e + k 6 F b c Q q l E p s q S f X J A j o h H T k m J X J E y q R J O g N y T R / L k 3 D g P z r P z M i 3 N O L O e X f J L z u s 3 G o G Q 2 A = = < / l a t e x i t > the mixed anti-commutators (8.3) vanish for u ± , v ± ∈ [a, b], we can analyse (8.10) and (8.11) and solve them exactly as in Sec. 5, finding where P (ξ; u ± ) has been defined in (5.12). At t = 0 the expressions in (5.11) are recovered.
The solution (8.12) for the modular flow of the Dirac field allows us to write the coordinates of the generic point along the modular trajectory in the spacetime where p ±,0 denote the light-cone coordinates of the initial point (at τ = 0) of the modular trajectory. In Fig. 3 we show some modular trajectories obtained through (8.13). A filled marker denotes an initial points p ±,0 = u ±,0 , while the empty version of the same marker denotes the conjugate initial point such that p ±,0 = ab/u ±,0 . These are the initial points of two conjugate modular trajectories (a solid curve and a dashed curve in Fig. 3) that determine the modular evolution of the field in the spacetime. In Fig. 3 three pairs of modular trajectories are shown. Notice that the modular trajectory passing through the point whose light-cone coordinates are u ± = √ ab coincides with its conjugate trajectory (see the solid black curve and the yellow dashed curve in Fig. 3). Changing the initial point, the resulting pairs of modular trajectories span the entire grey diamond D in Fig. 3. Any modular trajectory begins in the lower vertex of D as τ → −∞ and ends into the upper vertex of D as τ → +∞.
The correlation functions along the modular flow in the light cone coordinates can be obtained from (8.12). By adopting the notation u r± ∈ [a, b] with r ∈ {1, 2} for the two points, in the vector phase we find In the axial phase, this analysis leads to

Conclusions
In this manuscript we have studied the modular Hamiltonians of an interval A = [a, b] ⊂ R + for the boundary conformal field theory defined by the massless Dirac fermion on the half-line, in both the inequivalent phases that can be introduced by imposing the global conservation of energy, where either the charge or the helicity is preserved. These modular Hamiltonians have been obtained in terms of the components of the Dirac field as the sum of the local term (3.25), where the energy density (3.26) occurs, and the bi-local terms given by (3.30) and (3.31). Each bi-local term involves fields that must be evaluated in two conjugate points x andx = ab/x within the interval A. Moreover, the modular Hamiltonians preserve the global symmetry characterising the phase, that is either the vector phase transformation (2.7) or the axial phase transformation (2.8). In this respect we recall that the bi-local terms in the modular Hamiltonians for the massless Dirac field defined on translationally invariant spaces [12,[29][30][31][32] preserve both vector and axial symmetry.
Another main result of this manuscript are the modular flows of the Dirac field, given by (5.15) in the vector phase and by (5.16) in the axial phase. These flows mix two modular trajectories that start at conjugate points at τ = 0 (see Fig. 1) of fields having different chirality in the vector phase and different charge in the axial phase. By employing these results, we have found the correlation functions of the fields along the modular flow, which are (6.5)-(6.8) in the vector phase and (6.17)- (6.20) in the axial phase. These correlators satisfy the KMS condition, which characterises the modular flow [2]. Furthermore, these correlation functions satisfy the modular equation of motion determined by the modular flow (see e.g. (6.12)). The above results, first obtained at t = 0 have been also extended in the diamond region shown in Fig. 3.
As for the entanglement entropies of the interval on the half-line, the expression (4.11) has been found and, by considering the special case where the interval is adjacent to the boundary, the Affleck-Ludwig boundary entropy g = 1 has been obtained, in agreement with [19,50].
In the model considered in this manuscript only reflection occurs at the boundary. In a companion manuscript [52], the results obtained here are employed to study the modular Hamiltonians of equal disjoint intervals on the line placed symmetrically with respect to a defect that allows both reflection and transmission.
It would be interesting to find the modular Hamiltonians of bipartitions involving many disjoint intervals for other conformal field theory models (see [40-43, 47, 53] for the entanglement entropies). In [54][55][56][57] some entanglement hamiltonians for free quantum field theories have been obtained as the continuum limit of the corresponding entanglement hamiltonians of the lattice model [11,37,38,58]. The contours of the entanglement entropies are also interesting related quantities to consider [59,60]. It would be instructive to study also the modular Hamiltonians found in this manuscript through lattice calculations. The modular Hamiltonians for free quantum field theories in higher dimensions and in the presence of boundaries are also natural operators to explore (see e.g. [61] for the entanglement entropies). It is important to study also the modular Hamiltonians in quantum field theories where the conformal symmetry does not occur [54,57]. They could provide important insights to understand some features of the renormalisation group flows, even in the cases where the conformal symmetry is broken only by boundary terms [62][63][64]. Finally, let us mention that it is important to address the above open problem also in the context of the AdS/BCFT correspondence [65][66][67][68][69][70][71][72].

A Spectral problem for two equal intervals on the line
In this appendix we report the explicit form of the solution of the spectral problem (3.4) for the union A sym ≡ [−b, −a] ∪ [a, b] ⊂ R of two disjoint equal intervals on the line [12].
The eigenvalues in (3.4) can be written in terms of the real parameter s as follows The eigenfunctions in (3.4) are where the function w(x) occurring in the phase is defined in (3.11), the functions m p (x) determining the amplitude are 3) and Θ A (x) can be written in terms of the Heaviside step function Θ as follows This expression is non vanishing only in A sym . In particular, Θ A (x) = +1 when x ∈ (−b, −a) and Θ A (x) = −1 when x ∈ (a, b).

B Modular flow from a shift operator
In this Appendix we discuss the solution of the partial differential equation that allows to determine the modular flows reported in Sec. 5. In the Appendix B.1 the solution of the general case is derived and in the Appendix B.2 we describe its application to some cases characterised by local modular Hamiltonians.

B.1 General case
The general form of the partial differential equation underlying the modular flows discussed in this manuscript reads where ψ(x) corresponds to the given initial configuration of the field.
In order to solve (B.1), first one identifies a proper multiplicative factor that simplifies the form of the equation. In particular, let us redefine the field ψ(t, x) as where A(x) satisfies the ordinary differential equation V A + Y A = 0, which can be also written as d dx log A = − Y /V , whose solution reads Plugging (B.2) into (B.1), the partial differential equation for Ψ(t, x) becomes where Ψ(x) provides the known initial configuration. When Y (x) vanishes identically, we have that A(x) = 1 identically.
The unique solution of the partial differential equation (B.4) can be written in terms of the shift operator e t V (x) ∂x as follows [73] where we have introduced and which satisfies ξ(0, x) = x.
Finally, from (B.2), (B.3) and (B.4), we can construct the solution of (B.1) as follows where in the last expression we can employ (B.3) to get A(ξ(t, x)) = A(x 0 ) e −γ(ξ(t,x);x 0 ) . The final expression for the solution of (B.1) can be written as which is independent of x 0 . Since ξ(0, x) = x, it is straightforward to check that the solution (B.10) satisfies the initial condition ψ(0, x) = ψ(x), as required in (B.1).

B.2 Modular flows from local modular Hamiltonians
In various cases of physical interest [5,[7][8][9], the modular Hamiltonian of the interval A = [a, b] at t = 0 for the Dirac field reads where T tt (t, x) is the energy density and β 0 (x) is a weight function that characterises the underlying case.
The modular flow of the Dirac field generated by (B.12) is given by (5.1). It can be found by solving the following partial differential equation with a given initial configuration for the field This differential equation corresponds to the special case of (B.1) where V (x) and Y (x) are obtained from weight function β 0 (x) occurring in (B.12) as follows This implies that the modular flow ψ(τ, x) can be found by specialising to this case the general solution (B.10) derived in the Appendix B.1. The function w(x) in (B.6) and therefore also ξ(τ, x) in (B.7) cannot be written without the explicit expression of β 0 (x). Indeed, in this case the equations in (B.8) hold for ξ(τ, x) with V (x) given in (B.14).
Specialising the function γ(x; x 0 ) in (B.3) to the case defined by (B.14), we obtain hence the function multiplying the field in (B.10) becomes Thus, the solution (B.10) in the special case defined in (B.14) reads where (B.8) has been used. Since ξ(0, x) = x, the initial condition ψ(0, x) = ψ(x) is satisfied.
The solution (B.17) tells us that the modular flows generated by the modular Hamiltonians (B.12) are local because they do not mix fields localised in different points.

B.2.1 Examples with translation invariance
We find worth enumerating some interesting configurations where the bipartition involves a single interval A and whose modular Hamiltonian takes the local form (B.12). In these examples the underlying system is invariant under spatial translations.
Let us consider a configuration characterised by its weight function β 0 (x). This weight function allows to find w(x) through (B.6) and (B.14). Then, the function ξ(τ, x) is constructed through (B.7) and (B.14). The modular flow for the massless Dirac field can be found by specialising the expression (B.17) to the case of interest.
In the following cases, where the underlying system is translations invariant, the solution (B.17) can be employed to find the correlators of the components ψ r (with r ∈ {1, 2}) of the Dirac field along the modular flow More explicitly, in the following cases one finds that this correlator takes the form (6.16) where we have introduced the antisymmetric function G(x, y) ≡ g(x, y)/[e w(x) − e w(y) ]. We observe that (B.19) satisfies the following partial differential equation under the condition that G(x, y) is a solution of (we recall that w (x) = 1/β 0 (x)) In the following cases, it turns out that B w (x, y) is independent of the underlying bipartition and that G(x, y) is proportional to the two-point correlator of the entire system with x, y ∈ A.
The first example that we consider is given by a system in its ground state on the line, where the bipartition is defined by a single interval A = [a, b] ∈ R and by its complement [5][6][7]. The modular Hamiltonian takes the local form (B.12) with (see also (4.15)) The function w(x) provides ξ(τ, x) as follows For any x ∈ (a, b), we have that ξ(τ, x) → a as τ → −∞ and ξ(τ, x) → b as τ → +∞ in (B.23). Instead, ξ(τ, a) = a and ξ(τ, b) = b for any τ ∈ R.
The correlators along the modular flow are given (B.19) with It is worth considering the limiting regime of this first example where the subsystem becomes the semi-infinite line x > 0 (i.e. a = 0 and b → +∞), which corresponds to the case studied by Bisognano and Wichmann [3,4]. In this limit, from (B.22), we have β 0 (x) = x with x > 0 and the expression (B.23) reduces to ξ(τ, x) = x e 2πτ with x > 0, which are the dilations of the semi-infinite line parameterised by τ .
The expression (B.12) describes also the modular Hamiltonian of an interval A = [a, b] on a circle of finite length L, when the entire system is in its ground state. We remark that we consider anti-periodic boundary conditions. Non local terms can occur in the modular Hamiltonian if other boundary conditions are imposed [31,32,74]. In this case, which has been already studied in [8,9], the weight function β 0 (x) and the corresponding w(x) read respectively The expression of w(x) allows to write ξ(τ, x) as The last example is given by an interval A = [a, b] ⊂ R on the infinite line and a system at finite temperature 1/β. In this case the modular Hamiltonian has the local form (B.12) with [8,9] This expression for w(x) leads to the corresponding ξ(τ, x), which can be written as ξ(τ, x) = β 2π log e π(b+a)/β + e 2πb/β e w(x)+2πτ e π(b−a)/β + e w(x)+2πτ (B. 29) and to the corresponding correlators along the modular flow, which are (B.19) with Notice that the first case can be obtained either from the second one in the limit L → ∞ or from the last one in the limit β → ∞.

B.2.2 Modular flow of a primary field in CFT
A conformal field theory in two space-time dimensions contains the symmetric, conserved and traceless energy momentum tensor {T µν (t, x) : µ, ν = t, x}. Let us consider the chiral right and left-moving combinations of the energy-momentum tensor. Introducing A = [a, b], we are interested in the flows of a primary field of the theory generated by the modular Hamiltonians Consider, for instance, any right-moving primary field φ(v + ) with conformal dimension h. By employing the transformation law [75] and assuming that β 0 (a) = β 0 (b) = 0, we find which generalises (B.13) to an arbitrary conformal dimension h. Notice that the differential equation (B.34) determining the modular flow of φ is given by (B.1) with Specifying the solution (B.10) to this case, one obtains The modular evolution of a left-moving primary field of dimension h is obtained from (B.36) by replacing v + with v − .
We remark that the explicit form of the weight function β 0 (x) has not been specified in the above discussion. For the special case of a two dimensional conformal field theory in its ground state and the bipartition defined by the single interval A = [a, b], the modular trajectory ξ is given by (B.23) and the result (B.36) has been already reported in [5,6,15].

C Modular flows from bi-local modular Hamiltonians
In this Appendix we describe the derivation of the modular flows of the massless Dirac field generated by two different bi-local modular Hamiltonians. In the Appendix C.1 the case of the interval in the half-line is considered, obtaining (5.11), while in Appendix C.2 we discuss the case given by two disjoint equal intervals on the line, that has been solved in [12]. In the Appendix C.3 we report a partial differential equation satisfied by the correlators of the Dirac field along the modular flow generated by the modular Hamiltonian of two disjoint intervals in a generic configuration on the line.

C.1 Interval on the half-line
In the following we discuss the derivation of (5.11) as the solution of the system of the two coupled partial differential equations (5.9) with the initial condition (5.10). By using (3.28) to find −β loc (x)∂x = β loc (x) ∂ x , the system (5.9) can be written as which provides the coupling between the two equations.
In order to write the system (C.1) in a simpler form, we redefine the fields in Ψ(τ, x) as in terms of the function c(x) satisfying the condition where c 0 is a non vanishing constant. By plugging (C.3) into (C.1) and exploiting (C.4), one finds thatΨ(τ, x) must solve the following system which is independent of the constant c 0 and the constant matrix J α can be written in terms of the Pauli matrices σ 1 and σ 2 as follows The unitary matrix J α , that satisfies also J −1 α = −J α , provides the dependence on α in (C.6). This matrix can be diagonalised through the unitary matrix U α , namely We find it convenient to introduce the following linear combinations of fields Indeed, by employing (C.9) and (C.10) into (C.6), we observe that the partial differential equations for the fieldsμ 1 (τ, x) andμ 2 (τ,x) are decoupled. They read where the second equation tells us thatx is a convenient spatial variable forμ 2 .
Since one of the equation in (C.11) is obtained from the other one by exchanging (τ, x) ↔ (−τ,x) andμ 1 ↔μ 2 , the solution of this system can be constructed from the solution of a single partial differential equation given by This partial differential equation belongs to the class defined by (B.1), whose solution has been derived in the Appendix B.1. In particular, (C.12) corresponds to the case defined by Following the discussion reported in the Appendix B.1, one first introduces ξ(τ, x) from (B.7), that in this case becomes ξ(τ, x) = w −1 w(x) + 2π τ (C.14) where w(x) is given by (3.11) and its explicit expression of ξ(τ, x) has been reported in (5.13). Then, one specialises the integral in (B.3) to the case defined by (C.13), finding b(ξ) β loc (ξ) dξ = ω(ξ) + const ω(ξ) ≡ arctan ξ/ √ ab .

C.2 Two disjoint equal intervals on the line
In order to make this manuscript self-consistent, in the following we find it worth deriving the modular flow of the massless Dirac field in the ground state and the correlators of the resulting field when the subsystem is given by two disjoint equal intervals in the infinite line A sym ≡ [−b, −a] ∪ [a, b] ⊂ R, by adapting to this case the procedure discussed in the Appendix C.1 for the interval in the half-line. This analysis corresponds to a special case of the results obtained in [12,13] for the union of a generic number of disjoint intervals with arbitrary lengths.
The modular flow of the Dirac field is the solution of the following system of partial differential equations where the 2 × 2 matrix differential operator B sym (x) can be expressed in terms of B loc defined in (5.4) as follows . (C.29) The block diagonal structure in (C.28) implies that we can focus on the doublet given by Ψ(τ, x) ≡ ψ(τ, x) ψ(τ, −x) (C.30) that must solve the following system of partial differential equations .

C.3 Two disjoint intervals on the line: Modular equation of motion
The modular Hamiltonian of a subregion made by the union of a generic number of disjoint intervals on the line for the massless Dirac field and the corresponding modular flow have been studied by Casini and Huerta in [12], while the correlators the Dirac field along the modular flow have been found in [13].
The modular Hamiltonian of two disjoint intervals A ≡ [a 1 , b 1 ] ∪ [a 2 , b 2 ] on the line can be written as the sum K A = K loc A + K bi-loc A , where the local term K loc A and the bi-local term K bi-loc A are defined respectively as [12] K loc A = 2π where T tt (t, x) is the local operator (3.26), while T bi-loc (t, x, x c ) is the bi-local operator (3.38). The weight functions in (C.57) can be written as follows and x c (x) is the point conjugate to x ∈ A satisfying the condition w(x c (x)) = w(x), namely (C. 60) We remark that x and its conjugate point x c (x) belong to different intervals in A.
The correlators of the components of the Dirac field along the modular flow found in [13,15]  We remark that (C.63) satisfies the properties given in (6.3) and this implies that the modular correlators obey the KMS condition [13].
We remark that the occurrence of the last term in the l.h.s. of (C.64) is due to the fact that the modular Hamiltonian contains also a bi-local term. Indeed, for a single interval (a, b) on the line, the modular correlators (which are (C.61,C.62) with w(x) given by (B.22)) satisfy a modular equation of motion like (C.64) with β bi-loc (x) = 0 identically and β loc (x) = β 0 (x) defined in (B.22).