R\'enyi entropy and negativity for massless Dirac fermions at conformal interfaces and junctions

We investigate the ground state of a (1+1)-dimensional conformal field theory built with $M$ species of massless free Dirac fermions coupled at one boundary point via a conformal junction/interface. Each CFT represents a wire of finite length $L$. We develop a systematic strategy to compute the R\'enyi entropies for a generic bipartition between the wires and the entanglement negativity between two non-complementary sets of wires. Both these entanglement measures turn out to grow logarithmically with $L$ with an exactly calculated universal prefactor depending on the details of the junction and of the bipartition. These analytic predictions are tested numerically for junctions of free Fermi gases, finding perfect agreement.


Introduction
The entanglement content of extended quantum systems has been investigated in the last two decades in many different contexts, ranging from condensed matter [1][2][3][4] to high energy and black hole physics [5][6][7][8][9]. The most successful way to quantify the many-body entanglement is via the Rényi entropies: given a system in a pure state |Φ and a bipartition A ∪ B, the subsystem A is described by the reduced density matrix ρ A = Tr B |Φ Φ|, and the the Rényi entropies are S n (A) = 1 1 − n log Tr (ρ n A ) . (1.1) For n → 1 this definition gives the von Neumann entropy S(A) = −Tr (ρ A log ρ A ), often called just entanglement entropy. One of the most important use of the Rényi entropies has been the characterisation of critical one-dimensional systems: the distinctive feature is the logarithmic divergence of the entanglement entropy with the (sub)system size and conformal field theory (CFT) provides universal predictions for the prefactor of such logarithm. For example, the vacuum Rényi entanglement entropy of an interval A of length embedded in an infinite system is given by [10][11][12] S n (A) = c 6 1 + 1 n log ε + . . . , (1.2) where c is the central charge and ε is a ultraviolet (UV) cutoff. This remarkable scaling behaviour is altered at leading order by the presence of a boundaries [11]. For conformally invariant boundary conditions (bc's), the entanglement entropy can be studied via boundary CFT [13][14][15], a framework that already found a large number of applications in condensed matter and particle physics, such as quantum impurity problems [16], the multi-channel Kondo problem [17], D-brane physics [18] etc. For a finite size CFT of length 2L with conformal invariant bc's at the two edges, the Rényi entanglement entropy between the half-chain A = [0, L] and the other half is [11,12] S n (A) = c 12 1 + 1 n log L ε + . . . , (1.3) up to finite terms that depend on the bc's. At leading order in L/ε → ∞, Eqs. (1.2) and (1.3) differ by a factor 2. This is heuristically understood because the two geometries differ by the number of entangling points and, in general, one expects the entanglement entropy to be proportional to the size of the boundary of the subsystem. In both geometries mentioned above the origin of the entanglement relies on the presence of completely transmissive entangling points, resulting in some degree of quantum coherence among the subsystem and its complement. Conversely, when the entangling points are completely reflective because of some additional boundary conditions, the subsystems decouple, and the entanglement entropy between them vanishes. A natural generalisation of the above scenarios regards the intermediate setting in which the entangling points are partially transmitting and reflecting [19]. In the literature, such special situations are known as permeable interfaces, defects, or impurities (and indeed we will refer to them using all these equivalent names). A crucial result is that for free massless theories the defect is marginal [20] and so can alter the leading behaviour of the entanglement entropy. Conversely, interactions make the defect either relevant or irrelevant [20] ending up asymptotically in a completely reflective or transmitting situation, respectively, as shown also by the scaling of the entanglement entropy itself [21,22]. For free theories in the presence of a conformal interface (i.e., scale invariant) Sakai and Satoh exploited boundary CFT to show that the scaling of the entanglement entropy, (1.3) for n = 1, is modified as [23] S(A) = c eff ( √ T ) 6 log L ε + . . . , (1.4) where T is a parameter which represents the transmission probability and c eff ( √ T ), dubbed as effective central charge, is a monotonic function of its argument satisfying c eff (0) = 0, c eff (1) = c. (1.5) Ref. [23] focuses on the free massless boson, but Eq. (1.4) with a different c eff ( √ T ) has been subsequently derived also for free massless fermions both by means of CFT [24][25][26], and explicitly solving microscopic models [27][28][29][30][31] in the same universality class. While the scaling in Eq. (1.4) is expected to be a generic feature of conformal invariant (1+1)-dimensional systems, the explicit functional form of the effective central charge depends both on the theory and the details of the interface, eventually encoded in a interface operator (or, equivalently, in a boundary state as explained in [24]). We mention that a class of completely transmissive interfaces, dubbed topological interfaces, has been also considered in the literature [32][33][34][35][36][37][38]. While their effective central charge is always c, and they could be erroneously considered trivial, the O(1) terms shrugged off in Eq. (1.4) still contains important information about the boundary conditions, strictly related to the boundary entropy of Affleck-Ludwig [39].
The permeable interface between two CFTs can be generalised to a junction of M wires. The resulting geometry is depicted in Fig. 1 in which the junction is fully characterised by a scattering matrix S between the wires. Imposing that this matrix S preserves conformal invariance, one finds consistency conditions that have been studied and solved for a large number of physical configurations [40][41][42][43][44][45][46][47][48]. The bipartite entanglement in these conformal (or star) junctions has been studied in Refs. [29,[49][50][51] but focusing on the entanglement between a single wire and the remaining M − 1 ones. A unifying framework to compute the entanglement of a generic bipartition among the wires of the junctions is still missing.
The conformal junction is also a very obvious setup for the study of multipartite entanglement because it is made of several wires and it is very natural to wonder about the entanglement between a subset of them, not only two complementary subsystems. In this respect, the first configuration that comes to mind is the tripartition in A, B, C with M A , M B , M C wires each, as depicted in Fig. 1. To study the entanglement of this tripartition, one can integrate out the M C wires in C, to get the reduced density matrix ρ A∪B . Then the entanglement between A and B with the mixed density matrix ρ A∪B is measured by the negativity [52,53] where T B denotes the partial transposition with respect to the degrees of freedom in B and || · || stands for the trace norm. The negativity in the presence of a defect has been computed for a bipartite geometry with M = 2 [54], exploiting its relation with the 1/2-Rényi entropy for the bipartition of a pure state, but for a genuinely tripartite geometry at a junction there are no results yet.
The main goal of this manuscript is to provide a general framework to deal with the entanglement through permeable junctions of M (1 + 1)-dimensional free-fermion CFT. Following Refs. [19,23], the strategy is to constrain the form of the general boundary state in a folded theory. Then, being the theory free, we can reduce the problem to the computation of a charged partition function in the presence of this boundary state. This approach also allows us to compute the negativity in a tripartite geometry by properly implementing a partial transpose operation for free fermions.
The paper is organised as follows. In section 2 we review the folding trick which turns the problem of constructing conformal interfaces into the one of building boundary states. We review the construction of fermionic boundary states and we compute the partition functions in the junction geometry. Using this result and the replica trick, we obtain the entanglement entropy analytically for a generic bipartition between wires. In Section 3, we combine the previous formalism with the replica trick for the fermionic negativity. This allows us to obtain an analytic prediction that we benchmark against numerical computations in Section 4. In the same section, we also describe an alternative technique for the computation of the entanglement of a fermion gas on a star graph modelling the junction of interest. We draw our conclusions in section 5 and we relegate some technical material about our computations in the Appendix A.
2 CFT approach: description of the method and the application to the Rényi entropy In this section we present the CFT approach for the evaluation of the entanglement in permeable junctions of (1+1)-dimensional free-fermion CFTs, following closely Refs. [19,23,24]. As a first application, we employ this method to compute the Rényi entropies between an arbitrary number of wires at the junction. Let us consider M wires of length L, each of them described by a CFT denoted by In Euclidean space-time, the junction looks like a booklet (with each page corresponding to one CFT) bound along the imaginary axis at x = 0, see the left panel of Fig. 2. As custom in this kind of systems, we are going to work in the folded picture in which the system is represented as a single CFT i.e. the world-sheet is a single infinite strip of width L (we are dealing with a finite size quantum system at zero temperature, so the space-time coordinate w satisfies Re(w) ∈ [0, L]) Figure 2: The folding procedure. The junction in Euclidean spacetime is the booklet with the CFTs bound in x = 0 (left panel). The folding consists in merging together the M CFTs in a worldsheet being a single infinite strip with an appropriate boundary state |S at x = 0 (middle panel). To compute the entanglement, we cut the system for |w| < ε and w > |L| (dashed lines) and map the M-CFT onto a rectangle of size log L ε × π (right) with Re(z) ∈ [log ε, log L] and Im(z) ∈ [−π/2, π/2].
where M copies of the CFT live. See the middle panel of Fig. 2 for a pictorial representation. The joining between the distinct wires is specified by the boundary conditions along the lines We require that the boundary condition at Re(w) = L decouples the replicas, and can be thus described by a boundary state |B factorised as with |B j being a boundary state of CFT j . Instead, we assume that the boundary conditions at Re(w) = 0 (describing the defect/junction), in general couple explicitly distinct wires. We denote by |S the associated boundary state in M-CFT. In the remainder of the manuscript, the precise details of the boundary state |B appearing at Re(w) = L would not matter and so we do not specify more about it. The physical motivation is that, as long as it decouples the wires, we do not expect that its features affect (at least at leading order) the correlation properties among distinct wires. In contrast, this is not the case for the boundary state |S , and for this reason we have to be very careful about its characterisation. In order to have under control ultraviolet and infrared divergences in the entanglement entropy, a standard trick [10,[55][56][57][58] consists in cutting the theory for |w| < ε and |w| > L (see Fig. 2, middle panel). The cut strip can then be mapped into a rectangle by the conformal transformation z = log w. (2.5) The semicircles |w| = ε and |w| = L are mapped respectively onto the segments The defect line at Re(w) = 0 is split into the two lines Im(z) = ±π/2. This mapping is shown in the right panel of Fig. 2. The partition function in this geometry can be written as where π is the height of the rectangle (see Fig. 2), while the hamiltonian is (up to the Casimir energy which does not play any role in our discussion) with L 0 ,L 0 being generators of the Virasoro algebra of CFT 1 ⊗ . . . CFT M . So far, everything is general and no assumption on the bulk theory or the boundary state |S has been made yet. However, the knowledge of |S is required to evaluate the partition function (and, by replicas, the entanglement). From now on, we thus restrict the analysis to massless free fermions for which we can provide a precise characterisation for the boundary state |S .

Boundary states for free-fermions
In this section, we first review the construction of boundary states for a theory of many species of massless Majorana fermions [19,51]. Then we discuss the straightforward generalisation to Dirac fermions, obtained through a doubling of the degrees of freedom [59].
We consider M species of Majorana fermions. This CFT has central charge c = M/2 and it is described in terms of the left/right chiral fermionic fields ψ j ,ψ j , j = 1, . . . , M. (2.9) In radial quantisation [59], restricting the analysis to the Neveu-Schwarz (NS) sector, one can decompose the fermionic fields in their Laurent modes (In the Ramond sector, k would be integer and the discussion would be slightly more involved due to the presence of a zero mode for k = 0.) Within this convention, the creation/annihilation operators of a fermion of the j-th species in the mode k (k > 0) are ψ j ∓k . The number k is (proportional to) the momentum of the particle. More precisely, one can show that the commutation relations between the fermionic fields and the Virasoro operators The effect of the scattering matrix S at the junction (as in in Fig. 1) is nothing but a consistency condition for the state boundary state |S reading where, hereafter, repeated indices are summed over. It has been shown [25,60], that in order to preserve conformal invariance at the boundary, S must be orthogonal In particular the possible k-dependence of the scattering matrix is ruled out by scale invariance. The solution for |S of Eq. (2.12) is simply with |0 being the vacuum of the theory. Notice that the different values of k are decoupled, a fact that will simplify the forthcoming computations. Nevertheless, in general, different species of particles are coupled, due to the possible occurrence of non-diagonal terms in the matrix S. Those terms represent physically the amplitudes of transmission between different wires and cause the entanglement among them.
We now consider a theory of M free Dirac fermions (having central charge c = M ), for which the associated fields are where Ψ and Ψ † represent the particles/antiparticles respectively. This theory is equivalent to a theory with 2M Majorana fermions, and so the previous derivation is valid also in this case. The number of degrees of freedom is doubled and one should take an orthogonal real 2M × 2M scattering matrix S ∈ O(2M ). However, if we further impose that the global U (1) is preserved by the boundary conditions, there are additional constraints on the scattering matrix. This requirement corresponds to the property that a left/right particle can be produced from the vacuum (through the boundary state) together with its right/left antiparticle only. Requiring that this symmetry is preserved by the boundary conditions, we end up into a complex unitary scattering matrix S ∈ U (M ), (2.17) that constrains the boundary state |S as withS being the matrix complex conjugated to S. The solution of such constraint is A property of the state |S in Eq. (2.19) is that it contains two decoupled contributions, depending on the right and left moving particles. We will use this property to simplify the computations in the following sections.

Rényi entropies for a generic bipartition between wires
We describe how to compute the n-th Rényi entropy of a subset made up of M A ≤ M wires via the replica trick. Given a subsystem A of a generic QFT, the Rényi entropies (1.1) of integer order n can be obtained in a replicated theory with n copies of the QFT, i.e. in QFT ⊗n , which are cyclically joined along A by a branch-cut connecting the i-th and the (i + 1)-th replica [11]. The moments of the reduced density matrices can be then written in terms of a ratio of partition functions as [11,12] Tr where Z n is the partition function of the replicated theory while Z n 1 is just the partition function of a single replica raised to the n-th power.
In the case of the bulk free Dirac fermion, the partition function Z n can be further factorised using the replica diagonalisation as, e.g., shown in [61]. Within this method, the replicated partition function Z n becomes the product of n single-replica U (1) charged partition functions with flux α p = 2πp n (p = − n−1 2 , . . . , n−1 2 ), each of them denoted by Z 1 (α p ). The flux e iαp is inserted along the branch-cut of Z n that can be rewritten as In Ref. [61] this factorisation is derived by writing Z n as a 2n-dimensional Gaussian integral, whose diagonalisation leads to the product of n two-dimensional Gaussian integrals. Plugging Eq. (2.21) for Z n into Eq. (2.20) one has Tr (ρ n A ) = p (Z 1 (α p )/Z 1 ) in which the ratio Z 1 (α p )/Z 1 can be expressed as the vacuum expectation value of the operator associated to the action of the U (1) symmetry restricted to A, namely Here Q A is the charge operator which counts the difference between particles and antiparticles in the subsystem A, while |0 is the vacuum of the theory. Notice that these charged partition sums are the same appearing in the calculation of the symmetry resolved entanglement [62][63][64][65]. While Ref. [61] and most of the subsequent literature focus on the ground state of the system in the absence of boundaries, the same considerations apply more generically and in particular to the case of interest here. The reason is that the boundary state of interest (2.19) is Gaussian (it is an exponential of a bilinear of fermions) and thus the functional measure is Gaussian too: in other words, the theory is free both in the bulk and at the boundary (we stress that this property cannot be assumed a priori for a quadratic bulk theory, because there are interactions at the boundary that spoil the Gaussianity of the state, see for example [66]). Hence, in our specific case, we start from the theory M-CFT = CFT 1 ⊗ . . . CFT M and we replicate it n times, ending up with M-CFT ⊗n . Then, to compute Z n we perform a diagonalisation in replica space and end up with the product of n charged partition functions which are given by Eq. (2.7) with the insertion of the appropriate flux, i.e. (with our normalisation Here the modular parameter has been introduced for later convenience. Our goal then becomes the computation of in the limit q → 1, corresponding to L ε → ∞. For this purpose, we firstly decompose S, which is a unitary M × M matrix, in a block diagonal form Here to shorthand the species of A and B respectively. In this way, the charge operator Q A is where the summation over the index a is understood. We consider the contribution to the partition function (2.25) coming from the single Laurent mode k, which requires the evaluation of The commutation relations among q L 0 +L 0 e iαQ A and the fermionic fields are and they can be easily derived from the momentum/charge of the Laurent modes. Using q L 0 +L 0 e iαQ A |0 = |0 and the commutation relations (2.30), we get The last expression can be evaluated (see Eq. (A.1) in the Appendix), and we get Using the unitarity of S, SS † = S † S = 1, one can show that (see Eq. (A.10) in the Appendix), where the proportionality constant is an unimportant α-independent prefactor (see the appendix). Putting all the pieces together and taking into account the contribution coming from exchanging Ψ ↔ Ψ † , we find the analytic expression of the U (1) charged partition function in Eq. (2.25) According to Eq. (2.21), the n-sheeted partition function Z n can be written finally as as the contribution coming from the generic eigenvalue T a ∈ Spec(1 − S † AA S AA ), one has where T a can be interpreted as generalised effective transmission probabilities. Plugging this relation in the definition of the Rényi entropy in Eq. (1.1), one gets being the Rényi entropy associated to each T a . For the sake of completeness, we provide the explicit result for the partition functions and for the entanglement entropies in the relevant limit L ε → ∞. Since the total entropy is just given by the sum of M A independent contributions with effective transmission T a it is sufficient to write only one term. For convenience, we also define a parameter α , being a function of α and the effective transmission T a , satisfying 2 cos α = 2(1 − T a + T a cos α). (2.40) The infinite product appearing in Eq. (2.34) which gives the U (1) partition function is explicitly evaluated in Appendix (A.3), obtaining (2.41) In the limit q → 1, the leading term of the partition function gives which matches the one in Ref. [29] where also the analytic analytical continuation to n → 1 can be found and it is not repeated here. We stress that the major advance in this section compared to the existing literature [29,51] has been to understand how the elements of S AA combine (via the eigenvalues of (1 − S † AA S AA )) to give the entanglement entropy of more than one wire, while previous studies focused on a single one. This result is also preparatory to the calculation of the negativity reported in the following section.

CFT approach: Fermionic Negativity
In this section we apply the CFT formalism to the calculation of the negativity between two subsets of wires of a conformal junction. In particular, we consider here the partial time-reversal negativity (often just called fermionic negativity), which is a more suitable entanglement measure for fermionic systems in mixed states (see [67][68][69][70][71][72][73][74][75]). We will proceed via the evaluation of the Rényi negativity for even n = n e and then we will study the replica limit n e → 1.

Rényi negativities
We consider the conformal junction of Fig. 1 with the subsystems A, B, and C formed by three sets of wires. We are interested in the negativity between A and B. The definition in Eq. (1.6) is not well-suited to study the entanglement properties in the context of fermionic systems (see [67][68][69]). To circumvent this issue, the partial time-reversal transformation of the reduced density matrix has been introduced [67] and here we use the same symbols ρ T B AB and E, as for the standard partial transpose operation, having in mind that we refer always and only to the fermionic negativity. The replica approach to the negativity [76,77] starts from the computation of the moments of the partial transpose reduced density matrix that can be written in terms of a ratio of partition functions as This path integral representation of the moments is similar to the one of the reduced density matrix in Eq. (2.20), but here Z n is replaced byẐ n . The latter is the partition function in the n-sheeted Riemann surface built in such a way to implement the partial time reversal transposition in the subsystem B (see Refs. [76,77] for more details on the partial transpose and [67,68] for the fermion case). The negativity is finally obtained as [76,77] i.e. by taking the analytic continuation from the even sequence of replicas, n = n e .Ẑ n can be further factorised using the replica diagonalisation, such that it becomes the product of n single-replica U (1) charged partition functions, similarly to what has been done for Z n in the previous section, but with some differences. Let us focus on even n = n e , which is the only necessary object to compute the negativity (instead for the negativity spectrum, i.e. the spectrum of ρ T B AB also odd values of n matters [68,78], as well as for other entanglement witnesses [79,80]). The needed charged partitionẐ 1 (α) has twisting phases equal to e iα in A and e i(π−α) in B, i.e. it reads [67,72] The operator e iαQ A implements the U (1) symmetry restricted to A, while e −i(α−π)Q B inverts the flux (α → −α) and it introduces an additional phases −1 along B, which is the combined net effect of the partial transpose operation (or, equivalently, partial time reversal) on fermionic systems. The final result of this approach is that the n e -th Rényi negativity can be computed as In the presence of boundaries, the U (1) charged partition function straightforwardly becomeŝ Eq. (3.4) with (3.5) holds for a generic tripartite fermionic system with no assumption. To proceed for the calculation of the wire junction we restrict to the following specific situation: • M A , M B , M C = 1, so that the total number of wires is M = 3.
• S is not only unitary but also Hermitian, which means that S 2 = 1 and its eigenvalues can be just ±1. For some physical systems (including the Schrodinger junction in the next section), the hermiticity of the S matrix is a necessary condition for physical consistency. Hence this is not at all a very restrictive assumption.
With these working assumptions it is possible to obtain nice analytic results in a rather compact form. More general expressions (e.g. for more wires) can also be obtained, but at the price of more cumbersome computations and less intelligible final results without any major physical insight. It is clear from Eq. (3.5) that for M = 3 the key object to be evaluated is where we used the commutation relations between q L 0 +L 0 e iαQ A −i(α−π)Q B and the fields (see Eq. (2.30)), and the formula (A.10), which provides the vacuum expectation value. The determinant of the 3 × 3 matrix appearing in (3.6) can be evaluated directly, but it is useful to discuss first the constraints due to the unitarity of S. We define the matrix O as and we verify the following properties: • O is unitary (OO † = 1) and its eigenvalues are phases; • det(O) = 1 and the product of the eigenvalues is 1; • O = SO † S and the spectrum of O is thus invariant under complex conjugation, a feature that relies on our hermiticity assumption S = S † .
These properties imply that the spectrum of O has to take this form withα a real parameter depending on S and α, defined by the following property which is a consequence of Eq. (3.8). The determinant det(1 + q 2k O) can be thus computed taking the product over the eigenvalues of O as follows (3.10) Evaluating Tr(O) and using again the unitarity of S, we can rewrite Eq. (3.10) as det 1 + q 2k O = (1 + q 2k )(1 + 2 cosα q 2k + q 4k ), (3.11) and the explicit expression ofα as a function of the S matrix is (3.12) Notice thatα does only depend on the diagonal entries of the matrix S and on the flux α. Putting all the pieces together, we express the partition functionẐ 1 (α) aŝ (3.13) We find the same formal structure of the partition function which appeared for the Rényi entropies in Eq. (2.34), up to the replacement α →α . Analogously to Eq. (2.42), for L/ 1, we have withα given by Eq. We conclude this subsection by providing few simple consistency checks forα in some limits.
If the wire C is decoupled from the other two, then S 2 CC = 1 and the transmission probability between A and B is 1 − S 2 AA = 1 − S 2 BB . In that case 2 cosα = 2S 2 AA − 2(1 − S 2 AA ) cos(2α).

(3.16)
This value ofα is the same one would obtain for α in the U (1) partition function of A and B in the presence of a flux −e i2α inserted along A. The reason is that in this limit the system becomes invariant under the U (1) symmetry generated by Q A + Q B and, thanks to one recognises an equivalence with the insertion of e i(2α−π)Q A . Finally, we consider α = π/2 because it corresponds to the evaluation of the 2-Rényi negativity (indeed, from Eq. (3.4), E 2 = log(Ẑ 1 (π/2)Ẑ 1 (−π/2)/Ẑ 2 1 (0)) = 2 logẐ 1 (π/2)/Ẑ 1 (0)). In this case 2 cosα = 2S 2 CC , (3.18) and there is no explicit dependence on S AA , S BB . Now this parameter is the same α one would get in the presence of a flux e iπ/2 along C only. In this way, we reproduced the general identity [76] Tr which is the well known relation between the 2-Rényi negativity and the 2-Rényi entropy.

Analytic continuation
The representation of E ne as a sum, appearing in Eq. (3.15), gives an expression valid only when n e is an even natural number. To proceed to the calculation of the negativity, we should provide its analytic continuation for n e being a generic number in the complex plane. To this goal, the strategy we device is the following: • Using the identities of Appendix A.3, the U (1) partition functionẐ 1 (α) in Eq. (3.14) can be expressed through an integral representation in the limit q → 1 • The sum over the value of fluxes (3.15) can now be performed inside the integral. Through some simple trigonometric identities, this leads to an analytic continuation of the integrand.
• The final result is an integral, which represents our analytic continuation.
As aforementioned, using the trigonometric identities studied in the appendix A.4, we find with a, b, c being the following functions of the S matrix and t Using the integral representation for the product over the k modes in the limit q → 1 reported in A.3, we get Eq. (3.23) is the desired analytic continuation. At this point, the negativity is simply obtained by taking n e → 1: (3.25) Let us conclude this section by providing some useful cross-checks of our result. In the limit in which S 2 CC = 1, the wire C decouples and one recovers the result for two wires, which is given by the Rényi entropy with n = 1/2 (A and B now form a pure state) [54]. In this case, Eq. (3.23) simplifies as (3.26) Performing the change of variables t = e 2x , one can recover the result (in the integral form) for the Rényi 1/2 obtained in [29] or, by solving the integral, the result for the negativity between two CFTs in [54], with s = 1 − S 2 AA .

Schroedinger junction
In this section, we describe a fermion gas on a star graph modelling a junction made up of M wires of length L, joined together through a single defect. We introduce a slightly different framework (compared to the existing ones in the literature) which allows us to efficiently perform exact numerical computations also for the negativity. The CFT predictions of the previous sections are checked against these exact numerical results.

Correlation functions
Let us consider a star graph like the one in Fig. 1 where now on each wire there is a gas made of N spinless fermions. The M wires are decoupled everywhere but in the vertex of the graph and their mixing is described by a non-trivial scattering matrix. We consider the ground state of such system with N particles. The same system has been studied in Ref. [29] by the overlap matrix approach [49,82] which is the starting point of our analysis. Each point of the junction is parametrised by a pair where j is the index identifying the wire and x > 0 the spatial coordinate along the wire. The bulk hamiltonian of the system is with Ψ j , Ψ † j being the fermionic fields associated to the j-th wire (also called Schroedinger field, from which the name Schroedinger junction). We consider a scattering matrix S ij , i, j = 1, . . . , M, (4.3) describing the defect at x = 0, which has to be hermitian and unitary [29,43] S = S † , SS † = 1.

(4.4)
The most general boundary condition along the junction is where Ψ = {Ψ j } j=1,...,M , λ is an arbitrary real parameter with the dimension of mass. To fully specify the problem, we also need to impose boundary conditions at the external edges of each wire that generically take the form where µ i are again real parameters with the dimension of mass. In order to simplify the treatment, it is possible to diagonalise S via a unitary transformation U and its eigenvalues are just ±1. It is custom [43,44] to introduce a set of unphysical fields {ϕ j (x)} so that in terms of these ϕ j (x), the boundary conditions decouple as where the µ i 's are linear function of the µ i 's whose form is irrelevant. For the junction to be scale-invariant, we require each of the dimensionful parameters η i and µ i to be either 0 or ∞. The choice corresponds to either Neumann (∂ x ϕ i = 0) or Dirichlet (ϕ i = 0) boundary conditions at x = 0 and x = L. We impose Dirichlet boundary conditions (µ i = ∞) at x = L for all wires. Conversely, the values of η i being 0 or ∞ depends on the diagonalisation of the S matrix, see [29,43,44]. Hence, the unphysical fields ϕ j (x) have Dirichlet bc's at x = L and either (Neumann or Dirichlet) at x = 0, so that it is natural to use the short-hand notation N D Neumann-Dirichlet, DD Dirichlet-Dirichlet, (4.9) to refer to the two possibbilities. For these two possible boundary conditions, the singleparticle wavefunctions are

(4.10)
We work in the ground state with fixed particle number N for each wire, so that the correlation function is  Going back to the physical fields {Ψ j } j , linear algebra straightforwardly gives The matrices 1±S 2 are the projectors over the eigenspaces of S with eigenvalues ±1 respectively. In fact, given any eigenvector v ± of S satisfying Sv ± = ±v ± ,, by inspection it holds (4.14)

Finite-dimensional representation of the correlation function
The correlation functions (4.13) are continuous kernel of the spatial variables. While it is possible to work directly with such kernels (as done, e.g., in Refs. [83,84]), it is more convenient to work with a finite-dimensional representation of such correlation. In this subsection we derive a representation which is particularly useful for numerical applications and it is equivalent to the overlap matrix approach [49]. The main result can be read in Eq. (4.28).
Hereafter, we set L = 1 without loss of generality. We start noticing that C ij (x, y) can be thought as an operator acting on the Hilbert space C which is not orthonormal. Indeed, since the single-particle eigenfunctions are normalised in both ND and DD sectors, one has e n , e n = dx φ DD (n, x)φ DD (n , x) = δ n,n , n, n = 1, . . . , N, (4.16) and similarly if n, n = 1 + N, . . . , 2N . Instead, for n ≤ N and n ≥ N + 1, their scalar product Q n,n is Q n,n ≡ e n , e n = 2 . (4.17) Since the basis is not orthonormal, we have to be careful to correctly give a matrix representations of C N D (x, y), C DD (x, y), and thus of C ij (x, y). We introduce the dual space H * 0 , as the space of linear functional on H 0 , and the dual basis defined by e * n (e n ) = δ nn . To avoid confusion we use another symbol for the bra associated to e n , which is denoted by 20) and it is defined by e † n (v) ≡ e n , v , ∀v ∈ H 0 . (4.21) We stress that e * n = e † n , since e † n (e n ) = Q n,n , e * n (e n ) = δ nn , n ≤ N, n ≥ N + 1, which is a consequence of the non-orthonormality of the basis. The projectors C N D and C DD , seen as operators of H 0 and belonging to the space have the following expression which means that their associated 2N × 2N matrices are In conclusion, we represented the correlation function C ij (x, y) in Eq. (4.13) as a (2M N, 2M N ) matrix acting on as follows with Q being a N × N matrix defined by Eq. (4.17).

The Rényi entropy between two arbitrary sets of wires
A useful auxiliary quantity for the computation of the entanglement entropy and negativity is the matrix Γ = 1 − 2C (sometimes referred to as covariance matrix). Using Eq. (4.13) and the finite-dimensional representation of the correlation matrix in the previous section, we can express it as When Γ refers to the entire system, given that S 2 = 1, it is zero in the "wire space". However, Eq. (4.29) has a very convenient form for the restriction of Γ to a subsystem A made of M A wires because of the tensor product structure in internal and spatial coordinates. It is in fact enough to replace S → S AA , i.e. the projected S-matrix to obtain Γ AA , the covariance matrix restricted to the subsystem of interest: with S AA being a M A × M A matrix. In general 1 − S 2 AA ≥ 0, and so 1 − Γ 2 AA is positive semidefinite. We recall that the matrix representation for (C N D − C DD ) 2 is Since 1 − QQ † and 1 − Q † Q have the same spectrum, the spectrum of (C N D − C DD ) 2 is obtained by two copies of the spectrum of 1 − QQ † . To start, let us recover the results [29] for the case when A is a single wire. First, taking M = 2 and the completely transmissive S-matrix the correlation function of a single wire A satisfies and so the eigenvalues of 1 − Γ 2 AA are just rescaled by a factor |s| 2 , in agreement with what known [29] from the overlap matrix. In the case of A being one of the M wires in a junction, in the above equation is enough to replace s with the transmission coefficient of A.
Once the covariance matrix Γ AA is known, the entanglement Rényi entropy is [85] S n (A) = 1 1 − n Tr log 1 + Γ AA 2 Using the basic property about the spectrum of a tensor product we get that the Rényi entropy between A and the complementary wires B is where the T a 's are the eigenvalues of 1−S 2 AA , which play the role of a transmission probability, and S n,Ta is the Rényi entropy of a single wire with transmission probability T a . This analytic results for the microscopic model perfectly match Eq. (2.43) in CFT. We also notice that the relation S n (A) = S n (B) comes from the fact that 1 − S 2 AA and 1 − S 2 BB have the same non-zero spectrum (i.e., the same eigenvalues up to the vanishing ones).
To conclude this subsection, we present a numerical test for the validity of the CFT result for the logarithmic scaling of the Rényi entropy. Since the case of A consisting of a single wire has been discussed and tested in Ref. [29], we focus here on a four-wire junction and the subsystem A consisting of two wires. The S matrix is chosen of the form (4.39) The numerical results are reported in Fig. 3 finding a perfect agreement with CFT.

Entanglement negativity
We now consider a tripartition A∪B ∪C, where A (B) contains M A (M B ) wires, and we study the entanglement negativity between A and B. This amounts to project the scattering-matrix S over a subset of rows/columns belonging to A ∪ B. In particular, we denote from which we construct the matrix [67,69]  The latter matrix Γ × A∪B is the crucial object to write the Rényi negativities E ne which indeed are [67] E ne ≡ log Tr (|ρ A∪B | ne ) = The above equation is valid for arbitrary real n e (i.e. also for a non-even integer) and so the negativity E is obtained just by taking n e = 1.
Eq. (4.45) gives the Rényi negativities in terms of the correlation matrices that, once numerically evaluated, provides a test of the CFT results for the coefficient of the logarithm obtained in Section 3. For the numerical evaluation, we focus on a three-wire junction and on the two-parameter family of scattering matrices given by (4.46) We select as subsystems A and B the first two wires and compute numerically the Rényi negativity for several values of s, θ, and N . In Fig. 4 we reported the coefficient of the logarithm obtained as follows. We fixed s = 0.75 and we selected some values of θ; for each value of (s, θ), we calculated numerically the negativity, for several values of N up to 200. We fitted the obtained numerical results with a log N + b 0 + b 1 N −1 . Fig. 4 finally reports the best fit of a as a function of theta and compares it to the corresponding analytic result in Eq. (3.23), finding perfect agreement To conclude, in Fig. 5 we report the N dependence of the Rényi negativity and we benchmark the prefactor of the logarithmic term in Eq. (3.23) for different pairs (s, θ) and different replica indices n e .

Conclusions
In this manuscript we investigated the entanglement entropy and negativity for the ground state of M species of free massless Dirac fermions coupled at one boundary point via a conformal interface/junction. We consider a bipartion into two complementary sets with M A and M B wires each. We generalised the CFT approach (introduced in Ref. [23] for M A = M B = 1) to evaluate the Rényi entanglement entropies for arbitrary M A . In terms of the wire length L, the resulting Rényi entropy is where C n (S AA ) has been explicitly calculated and turns out to depend only on the scattering matrix reduced to the subsystem A, denoted by S AA . One interesting aspect of our result is that the S n (A) is the sum of M A independent contributions which have the form of Rényi entropies of a single wire with a transmission probability given by the eigenvalues of 1−|S AA | 2 (see Eq. (2.39)). We obtained the same result also for a microscopical model of a free Fermi gas on a star junction. To do so, we generalised the approach of Ref. [29] for M A = 1 to the case of a subsystem A made of more wires. Such generalisation also provided numerical tests for the correctness of the CFT prediction.
We then moved to the study of the entanglement between two edges embedded in a multiterminal junction, considering the fermionic negativity E of Ref. [67]. To this aim, we adapted the CFT formalism above to the computation of the Rényi negativities between two wires of a tripartite geometry. As for the entanglement entropies, we found that the negativity grows logarithmically with the system size L, with a prefactor that depends on the details of the junction. Also the negativity has been explicitly constructed and analysed for a free Fermi gas on a star junction, finding results fully compatible with the CFT.
We stress that our results apply to some different physical situations, too. For example, our predictions are expected to hold also for lattice models of free fermions, in particular for M tight binding chains of length L joined at a single common vertex (as, e.g., done in [30] for M = 2). Furthermore, the logarithmic prefactors we obtained for entropy and negativity should appear also in the study of a star junction of M infinite CFT, but when the subsystems consist of segments of length starting from the interface.
We conclude the manuscript discussing few outlooks. The main focus of this work has been the free-fermion CFT, but our formalism can be easily adapted to free complex boson too, to study, e.g., the entanglement entropy and negativity across junctions of harmonic chains. The same formalism can be further applied to the study of some out-of-equilibrium protocols for free CFTs in the presence of defects (see, e.g., [22,31,54,86]). Another interesting open problem is the generalisation to other CFTs such as the compact boson, WZW models (see [35] for the interfaces of WZW), or minimal models. In all those cases we do not expect any kind of replica diagonalisation, due to the lack of Gaussian measures, but still one could employ the replica construction to investigate the negativity, as well as other entanglement measures.

A Useful identities
In this appendix we report some technical details about the calculations appearing in Sections 2 and 3.

A.1 Expectation value of Gaussian operators
We want to prove the following identity for the expectation value of Gaussian operators Here ψ j ∓k is the creation/annihilation operator of a Majorana fermion in the k-th left mode of the j-th species (among the M ones), whileψ j ∓k is the corresponding right mover. O jj and O jj are M × M matrices and the sum over j, j is implicit.
Before proving Eq. (A.1) in the most general case, we highlight simple cases in which it holds. Let us suppose there is just one species of fermions, so that we can suppress the indices j, j and the matrices O and O become numbers. Using that the annihilation operators annihilate the vacuum, that the fermionic operators square to zero, and applying Wick theorem, one gets Similarly, when O and O commute, one can diagonalise them simultaneously and apply the previous consideration to show Eq. (A.1).
We provide a general proof of Eq. (A.1) using Gaussian integrals over Grassmann variables (whose basic properties can be found on [87]). We start by representing the Gaussian operators as integrals over Grassmann variables. In particular, for each j-th Dirac fermionic field we associate a pair of Grassmann variables η j ,η j , and we express the Gaussian operator where the sums over j and j are implicit. Similarly, we introduce a pair of Grassmann variables θ j ,θ j to each species j and we express exp O jj ψ j The last step is the evaluation of the Gaussian integral over the 4M Grassmann variables η j ,η j , θ j ,θ j . We introduce a 4M -dimensional vector Θ of Grassmann variables as follows and the Gaussian integral in (A.5) as withÕ being the following 4M × 4M matrix The integral over Θ in Eq. (A.7) gives the Pfaffian of the matrixÕ [87], which we write as This completes the proof of Eq. (A.1), which is the main result of this section.

A.2 A useful determinant
In this section we show that for a complex unitary M ×M matrix S having the block structure (2.26) the following relation holds Before proceeding with the proof, we notice that the α dependence in the rhs above is related only to the non-zero eigenvalues of 1−S † AA S AA . Despite the explicit dependence of 1−S † AA S AA on A, this fact does not lead to any asymmetry between A and B. Since we have already shown that 1−S † AA S AA , 1−S † BB S BB have the same non-zero spectrum, we get det 1 + q 2k O ∝ det 1 + q 2k (2S † AA S AA + 2(1 − S † AA S AA ) cos α) + q 4k , (A.15) The α-independent proportionality constant has to be a power of (1 + q 2k ), which comes form the possible presence of zero eigenvalues of 1 − S † AA S AA (or 1 − S † BB S BB ). We match this constant by power counting. More precisely, since det 1 + q 2k O is a polynomial in q 2k of order M and det 1 + q 2k (2S † AA S AA + 2(1 − S † AA S AA ) cos α) + q 4k (A. 16) is a polynomial in q 2k of order 2M A , the right power of (1 + q 2k ) which matches the proportionality constant has to be (1 + q 2k ) M −2M A . This concludes the proof of Eq. (A.10).

A.3 Jacobi Theta functions and Dilogarithm
Here we review some properties of the Jacobi Theta functions and the dilogarithm, useful for the evaluation of the partition functions. We consider the following infinite product representation for the Jacobi theta function θ 3 (z, q) [81] θ 3 (z, q) = ∞ m=1 1 − q 2m 1 + 2 cos(2πz)q 2m−1 + q 4m−2 . (A.17) We want to evaluate θ 3 (z, q)/θ 3 (0, q) in the limit q → 1 − . To do so, we first take its logarithm, which turns the infinite product representation into a sum, i.e. In the previous computation, the integral representation of the dilogarithm function [81] Li 2 (z) = −

A.4 Trigonometric identities
We consider some useful algebraic identities which will be applied for the analytical continuation of the Rényi negativity. The first identity, which holds for any n ∈ N is We conclude this appendix with a straightforward consequence of (A.24), which is the evaluation of the product ne−1 2 p=− ne−1 2 x + cos 2πp n e y + cos 2 2πp n e z .
(A. 28) Since the term in the parenthesis is a second order polynomial in cos 2πp ne , it can be factorised in the two roots and Eq. (A.24) can be used for the computation of each of the two products in which it splits.