O ct 2 01 9 Constructing spherically symmetric Einstein-Dirac systems with multiple spinors : Ansatz , wormholes and other analytical solutions

In this paper we present a detailed calculation of an Ansatz that allows to obtain spherically symmetric Einstein-Dirac configurations in d-dimensions. We show that this is possible by combining 2 d−2 2 ⌋ Dirac fields, making use of the properties of the angular dependence of the spinors in a spherical background. By applying this Ansatz, we investigate some simple analytical solutions. One of them is a regular wormhole supported by the Dirac fields. Other solutions include a pathological black hole and a naked singularity. We analyze the domain of existence and properties of all these solutions.


I. INTRODUCTION
The study of the Einstein equations coupled to different classes of matter content has received a lot of attention in the recent years, since solutions of these theoretical models could be related to exotic astrophysical systems (composing the dark matter/energy sector) [1]. In higher dimensions these solutions could be of potential interest in the context of supergravity and the AdS/CFT correspondence [2].
There is an ongoing intense exploration of selfgravitating stationary soliton-like solutions composed by different classes of massive fundamental fields. The interest of these settings is because of the contrast with the more standard Einstein-Maxwell theory (and even with the Einstein-Maxwell-scalar theory [3]), where the electro-vac black hole is the only self-gravitating stationary soliton-like solution. But when fundamental fields are considered to be massive the situation changes. It is possible to construct particle like solutions (with regular, stationary space-times, but typically with a harmonic time dependence on the fields).
With scalar fields (in particular, massive and complex), these configurations are well-known to exist and originally obtained in [4,5]. Typically these configurations are known as boson stars [6], and they are considered potential candidates as astrophysical objects [7]. With massive vector fields (known as Proca stars) this was explored in [8] (some results in five dimensions can be found in [9]). Several astrophysical properties of these objects have been analyzed in depth [10][11][12][13].
From a more theoretical point of view, a particularly interesting case is to consider self-gravitating soliton-like solutions of the Einstein-Dirac system, where gravity plays the role of the non-linear interaction that allows for the existence of Dirac solitons in simpler models [14,15]. However, there is an additional challenge for this type of fields. Because of the intrinsic angular momentum of a single spinor field (a preferred direction in space-time), the resulting space-time is forced to rotate in order to accommodate stationary solutions. Such configurations have been very recently constructed in [16], where they were compared with rotating Boson and Proca stars.
Nonetheless, a possible route to enforce that the global solution of the Einstein-Dirac system is actually static and spherically symmetric, is to relax the single spinor condition and consider a collection of Dirac fields. In four dimensions, two fields are enough to cancel the intrinsic angular momentum and realize a global spherically symmetric and non-rotating space-time [17,18] (provided the fields possess a certain harmonic time dependence).
The properties of such multi-Dirac soliton-like solutions have been recently studied in [19], where their properties were compared with similar configurations made of scalar and Proca fields for d ≥ 4. It was shown that some generic features of the solutions actually do not depend qualitatively on the spin of the field, but they are controlled by the dimension of the space-time. These 'Dirac stars' have also been studied in the presence of vector fields [20].
On the other hand, solutions of the Einstein-Dirac system with multiple fermions are of interest in condensed matter [21], where in particular wormhole space-times with two Dirac fields can be used as effective models describing two graphene layers connected by a short nanotube. This model is known as the graphene wormhole [22][23][24][25][26][27][28][29][30][31] The purpose of the present paper is two-fold. First, we want to provide details on how the Ansatz for these multi-Dirac self-gravitating solitons is calculated for arbitrary space-time dimension. Second, we will show that, in addition to the numerical solutions previously obtained in [19], it is possible to obtain a few simple analytical solutions to the equations. One of these solutions is a regular wormhole solution supported by multiple Dirac fields. We will analyze the physical meaning and properties of these configurations.
The paper is organized as follows: in section II we present the general formalism for the Einstein-Dirac system and make an overview of how the Ansatz is built. In section III we explain how to combine the angular dependence of the different Dirac fields in order to get a spherically symmetric stress-energy tensor. In section IV we analyze the effective action and the minimum set of differential equations of this system. In section V we describe several sets of solutions that can be obtained in various particular cases. In section V A we present a regular wormhole solution supported by pairs of Dirac fields. In section V B we present a black hole solution with a pathological behaviour of the Dirac fields at the horizon. In section V C we present a light-like singularity. In section VI we end the paper with a summary and conclusions.

II. OVERVIEW OF THE GENERAL SETTING
We want to construct spherically symmetric solutions of the d-dimensional Einstein-Dirac system. There are many studies on the Dirac equation in higher dimensional, spherically symmetric space-times [32][33][34][35][36][37][38][39][40][41]. An appropriate metric Ansatz is where dΩ 2 d−2 is the line element of the (d − 2)-sphere. In order to build minimally coupled Dirac fields to this metric, we need to specify the vielbein , and for the metric Ansatz (1) we choose where j = 1, ..., d − 2 is an index running over the (d − 2)sphere and ω j d−2 is a vielbein for the (d − 2)-sphere. This allows us to write the Dirac equation, The operator K d−2 is the angular operator of the (d − 2)sphere, given by with Γ d−2 ij being the spin connection of the (d−2)-sphere, e d−2 j being the dual to the vielbein on the (d − 2)-sphere and γ j d−2 = γ t γ r γ j . Because of the spherical symmetry of the metric, the Dirac operator commutes with the angular operator, [D, K d−2 ] = 0. In addition ∂ t is a Killing vector. These two properties allow us to write a spinor with the following Ansatz where Θ κ depends on the angular variables only. The angular part is chosen to fulfill K d−2 Θ κ = κΘ κ , with κ the angular momentum eigenvalue. As we said in the introduction, because of a single spinor having a non-trivial intrinsic angular momentum, it is not possible to construct a compatible solution of the Dirac equation (3) with the metric (1).
A way out is to consider a system of multiple Dirac spinors. If we choose them appropriately, the combination of all of them will have a total stress-energy tensor compatible with the symmetries of the metric (1). To do so, we need 2 ⌊ d−2 2 ⌋ spinors (note ⌊x⌋ means the integer less than or equal to x). These spinors need to have the same radial function, but they differ in their angular parts, with the same (the smallest possible) angular eigenvalue combined incoherently. This means that, written as a formal sum, the spinors combine like where ǫ is the index of each one of the 2 ⌊ d−2 2 ⌋ spinors. In terms of the action, the Einstein-Dirac system with cosmological constant Λ can be written like where α g is the coupling constant between gravity and the spinor fields, and the Lagrangian for the spinor part is then a sum of the form We will show that this leads to a spherically symmetric energy momentum tensor for the spinors, where the T (ǫ) µν is the energy momentum tensor of the e −iωt φ κ ⊗ Θ κ,ǫ spinor.
We will now focus on the (d − 2)-dimensional sphere and the construction of this spherically symmetric configuration.

III. HOW TO COMBINE THE SPINORS
In the following we will make use of expression (5) for each one of the spinors. We will assume that: 1. all the spinors share the same radial dependence; 2. all the spinors share the same temporal dependence, and we will assume it can be written in terms of a phase, introducing the frequency ω; 3. the spinors only differ in the angular part.
In this section we will discuss in detail the properties of this angular part, and how it can be chosen in order to make the stress-energy tensor compatible with spherical symmetry.
A. Peeling the n-sphere Let n denote the dimension of the sphere. Solutions to the Dirac equation on the n-sphere are well-known and called spinor monopole harmonics in the literature [42][43][44]. It is however not trivial to combine these solutions into a field configuration which possesses a spherically symmetric energy momentum tensor. An approach to make this combination more intuitive, is to maximize the commuting Killing vectors on the n-sphere in our coordinate system. For this we will choose angular coordinates in such a way, that the line element on the sphere is given recursively by with and S n = sin θ n , C n = cos θ n . We thus slice off a twosphere from the n-sphere. This is convenient, because it allows us to define the vielbein on the sphere also in a recursive way, meaning ω θn n = dθ n , ω φn n = S n dφ n , ω j n = C n ω j n−2 , with j = 1, . . . , n − 2 being an index running over the (n − 2)-sphere. The spinor covariant derivative on the n-sphere ∇ (n) a is thus with γ Kn n = −γ θn n γ φn n and γ j n−2 = iγ Kn n γ j n and j as before. One can think of γ j n−2 as the γ j n matrices projected down onto the (n − 2)-sphere with γ Kn n governing this projection. The reason we choose this factorization of the γ-matrices will become clear later.
With these choices of line element, vielbein and algebra, we can write the Dirac operator on the sphere (meaning the angular part) as K n = γ a n ∇ (n) a = γ θn n ∂ θn + ∂ θn ln S n C n−2 where the index a runs over the n-sphere and K n−2 = γ k n−2 ∇ (n−2) k is the angular operator for the (n−2)-sphere. The matrices fulfill γ a n , γ b n = −2δ ab , with a, b ∈ {θ n , φ n , j} , (15) γ a n , γ b n = −2δ ab , with a, b ∈ {θ n , φ n , K n } , (16) γ Kn n , γ j n = 0 , γ a n , γ j n−2 = 0 , with a ∈ {θ n , φ n , K n } , where in the above j and k denote indices on the (n − 2)sphere. Equation (15) expresses the Clifford algebra on the nsphere. Equation (16) is the Clifford algebra on the 2sphere. Finally equation (17) is the Clifford algebra on the (n − 2)-sphere, showing us that the projection works correctly and we have sliced off a two sphere from the nsphere. The next equation (18) tells us that the matrix governing the projection onto the (n − 2)-sphere commutes with the matrices of the (n − 2)-sphere. The last equation (19) shows that the projected γ-matrices on the (n− 2)-sphere commute with the γ-matrices on the sliced off two-sphere.
This last relation implies that [K n , K n−2 ] = 0. In addition, since ∂ φn is a Killing vector we also have [∂ φn , K n ] = 0. Even more, for any n, m ∈ N we have that in the tower of angular operators [K n , K m ] = 0 and [K n , ∂ φm ] = 0.

B. Angular solutions of the spinor field
We have rewritten the angular operator into a tower of angular operators. We will study now what this is implying to the angular part of the spinor fields.
Let us study this equation a bit more. It is convenient to define Θ κn,mn = e − θn 2 γ φn n γ Kn n S n C n−2 nΘ n .
Substituting this into the differential equation (21) and multiplying with e θn 2 γ φn n γ Kn n from the left gives the following differential equation forΘ n γ θn At this stage, let us choose as a representation This leads to the coupled first order differential equation with Let us also define p m := m n + 1 2 , p κ := κ n−2 + 1 2 , where F (a, b; c; z) is the hypergeometric function and n κ ≥ 1 is a natural number. For solutions of the differential equation (25) m n is a half integer number, while |κ n−2 | ≥ (n − 2)/2 is an integer number in the case n is even, or a half integer number in the case n is odd.
With these definitions, we can write three solutions of equation (25).
1) The solution in the case K + = 0 is and the angular eigenvalue in this case is where ± κ is a sign choice.
2) In the case of K + = 0 the solution is with m n ≥ 1/2 and κ n−2 ≥ (n − 2)/2. The angular eigenvalue in this case is 3) Lastly in the case of K − = 0 the solution is with m 1 ≤ −1/2 and κ n−2 ≤ −(n − 2)/2. The angular eigenvalue in this case is With these three solutions at hand, let us analyze what are the smallest possible eigenvalues |κ n | we can reach and what are the corresponding angular solutions. An inspection of the solutions allows us to conclude that these are given by the cases K ± = 0 choosing |m n | = 1/2, |κ n−2 | = (n − 2)/2, and by K + = 0 choosing n κ = 1, m n = ±1/2, κ n−2 = ∓(n − 2)/2. This gives as a result |κ n | = n/2.
Note that the previous values of the angular parameters depend only on n and some possible sign choices. In the following we will choose only these minimum values. Hence, since for a given dimension n the only possible choices are the different signs of the angular parameters, for the sake of simplicity we can relabel the angular solution accordingly: Θ κn,mn ≡ Θ (n) sgn(κn),sgn(mn) . Note that the sign of κ n−2 is determined by these sign choices via sgn(κ n−2 ) = −sgn(m n ) sgn(κ n ).
Let us write explicitly what these solutions are. There are four possibilities with |κ n | = n/2: Tracing back our steps we can thus write the angular part of the spinor for the n-sphere with minimal absolute value of the eigenvalue |κ n | = n/2 as with ǫ κ being the sign choice for κ n , and the ǫ j being the sign choices for the m j (j > 1), which we can summarize as a binary vector ǫ (notice that either the even or the odd components of this vector are immaterial for us depending on n), and The sign of m 1 , so ǫ 1 , is fixed by the equation The reason for this is that m 1 plays a double role as an eigenvalue to ∂ φ1 and as the angular eigenvalue κ 1 to the one-sphere (circle).

C. Analyzing the properties of the angular solution on the components of the total stress-energy tensor
Now that we have this set of solutions (35) for the angular part of the spinor, we need to study how it enters into the stress-energy tensor. The stress-energy tensor for a collection of spinors was given in expression (9). From there we can see that it is useful to construct explicitly the matrix elements of the covariant derivative ∇ (n) a multiplied with a matrix Γ, since we will need these objects for the calculation of the total stress-energy tensor.
The first thing to do is to look at the following relations Using this and the inner product table gives the following matrix elements On the other hand, the covariant derivatives are explicitly given by In the case of n being odd we have For the construction of the spherically symmetric stressenergy-tensor we now need the expectation value of ∇ (n) k and γ k n ∇ (n) j with Θ ǫκ,ǫ . Write these as Γ = Θ † ǫκ,ǫ ΓΘ ǫκ,ǫ , for the matrix element of Γ. The following identity proves to be useful, with j ≡ k ≡ n mod 2, n ≥ k < j and α ∈ {θ, φ, K}, After some tedious algebra we find the following ex- where we have defined It is important to note that for l = k and the sum is over all possible sign vectors ǫ.
Hence the non-diagonal terms sum to zero. The non-vanishing sums are in the diagonal parts, which result in D. Combining the spinors We will use the above expressions to construct a field configuration with a spherically symmetric energy momentum tensor.
Fix a sign ǫ κ for a lowest angular eigenvalue κ = ǫ κ d−2 2 of the (d − 2)-sphere. Define 2 ⌊ d−2 2 ⌋ spinor fields as in equation (5), but labeling explicitly all the allowed sign combinations of ǫ. We then combine these spinors in an incoherent superposition so that written here as a formal sum ranging over all possible values for ǫ. Note in the last step of equation (49) we have made use of the fact that the radial and temporal dependence of each individual spinor is the same for all of them.
In the spacetime of the metric given by equation (1) the covariant derivatives are explicitly given by with j being an index on the (d − 2)-sphere and being the covariant derivative on the (d − 2)-sphere and γ j d−2 := γ t γ r γ j being the γ-matrices of the (d−2)-sphere. Using for the sphere the same vielbein as in the previous sections, and after some algebraic manipulations in which one needs to make use of the expressions we have derived in section III C, one arrives at the following energy momentum tensor (note that it is written in vielbein components) where we have used the radial equation for φ κ to simplify some expressions, One can easily see that this tensor is diagonal on the spatial components and thus spherically symmetric. Note however that the tensor has in general a nontrivial t − r component. This means the configuration in general has a radial flux, and will force the configuration to be time dependent. If we want to obtain solutions compatible with the static metric (1), we have to require the radial current to vanish everywhere, meaning φ † κ γ t γ r φ κ = 0 (no-flux).
The only thing left is to choose a particular representation of the remaining γ-matrices and spinor components, The no-flux condition reads The following parametrization incorporates the no-flux condition withφ a complex function and ν a real valued function.
Using this Ansatz and representation in the equation (53) and after some algebra we get the following nonlinear first order system of differential equations forφ and ν The equation for ν forces the frequency ω to be real. Note that the phase ofφ does not vary with r and is not a dynamical quantity. The stress-energy tensor in the vielbein components simplifies into An important quantity we can calculate is the time component of the net current in the vielbein, the Dirac density: We can see that all the components of the stress-energy tensor are proportional to the Dirac density j 0 net .

E. A comment on the time-dependent case
Although we are mainly interested in static metrics, the previous Ansatz can be easily generalized to accommodate the time-dependent case.
In this case the metric functions σ and N also have to depend on time. But it is also necessary to change the Ansatz for the Ψ ǫ to meaning there is no harmonic time dependence in the fields.
With these changes, the stress-energy tensor becomes The equation fulfilled by φ κ is now a partial differential equation, In the following we will consider only a static spacetime, and assume the Dirac fields possess a harmonic time-dependence.

IV. EFFECTIVE ACTION
With the construction we have developed in the previous section, it is possible to simplify the part of the action (7) containing the collection of Dirac fields, where we have defined the effective Lagrangian L eff .
With this the equations of motion read, using j 0 net = 2|φ| 2 , or after using the equation for ν ′ in the equation for N ′ and defining we have Another useful way to write this is usinĝ This means that With this the effective Lagrangian for the spinor part is especially simple The equations of motion are This form is useful for numerical calculations [19]. Finally, if we assume that σ > 0, which we can always do without loss of generality, there is a convenient way to redefine the spinor functions by setting This is convenient, because it makes the form of the field equations a bit more compact, This will be helpful in the next sections. Let us note here that with these definitions the Dirac density is

V. ANALYTICAL SOLUTIONS
Of course a fundamental question that immediately arises is if, for some set of parameters, the previous sys-tem of equations possesses physically meaningful configurations, and what is their interpretation. In the following, we will focus on cases without cosmological constant (Λ = 0).
Soliton-like solutions of this system in several dimensions have been presented in [19]. These solutions, to our knowledge, can only be constructed numerically. The solutions (sometimes called Dirac stars, although they are not expected to have any connection with realistic astrophysical objects) are regular everywhere, and share many features with similar self-gravitating soliton-like configurations found with massive bosonic fields.
In this section we want to present several analytical solutions that the previous system possesses. We will analyze in detail the physical and geometrical properties of these solutions.

A. Multi-Dirac wormhole
Let us specialize to a massless (m = 0) field which does not vary in time (ω = 0). The differential equations (72) simplify in this case to The differential equation for N means that N = 1 − (µ/r) d−3 with µ being a constant. Without loss of generality we can fix positive angular momentum of the fields, ǫ κ = 1, κ = d−2 2 . The solutions for f and g arê with c 0 , c ∆ ∈ R. Notice thatfĝ = −e c0 ∈ R ≤0 . This simplifies the differential equation for σ, which now reads This is easily integrated to be with c σ ∈ R. In total, the solution is parameterized by three real constants c σ , c 0 and c ∆ , in addition to µ, the coupling constant α g and the dimension d. However the parameters satisfy several relations.
Our first requirement is for σ to be a positive function in all of the domain r ∈ [µ, ∞). An analysis of equation (77) reveals that this is only possible if c σ > − 2e c 0 The second requirement is to reach the standard Minkowski metric at infinity. This means that Taking into account the previous condition for σ > 0, this implies the following relation for c 0 Hence we can write and it is easy to see that the g tt component behaves asymptotically like Let us explore the physical meaning of this metric. To simplify the discussion, let us first look at the particular case with c σ = 0. Thus σ(r) = 1/ √ N . The metric is very simple, with N = 1 − (µ/r) d−3 . This looks like the metric of a traversable wormhole [45]. Let us make the following coordinate transformation in this metric, ρ = √ N . This leads to For r ∈ [µ, ∞) we have ρ ∈ [0, 1) mapped such that r = ∞ → 1 = ρ. In this coordinate system, it is possible to extend the above metric (83) to ρ ∈ (−1, 1). Thus the above metric corresponds to a wormhole connecting two asymptotic regions at ρ = ±1. The sphere with minimal surface has radius r = µ, which corresponds to the throat of the wormhole as we will explicitly see later. Note that the Ricci scalar of this metric vanishes, but the Krestchmann scalar is finite. For d = 4, K = 6(1 − ρ 2 ) 6 /µ 4 ; for d = 5, K = 24(1 − ρ 2 ) 4 /µ 4 , etc... An interesting property of the above geometry is that the temporal part of the metric is essentially not curved and thus a test mass can rest at a fixed radius ρ without moving.
We can compute the mass using the standard Komar integral with ξ t being the one form dual to the Killing vector K t = ∂ t and the integral being over a (d − 2)-sphere at infinity. For metric (83) this is trivially zero. Hence the mass of this wormhole vanishes. Now let us choose the more general case with c σ = 0. Changing again the radial coordinate to ρ = √ N we have a metric slightly different than the one in the previous case, Again it is possible to extend the range of ρ from [0, 1) to (−1, 1). We can interpret this geometry also as a wormhole connecting two asymptotic regions at ρ = ±1. The difference now is that we have some non-trivial red-shift between universes. In fact note that g tt → 1 when ρ → 1, but g tt → (1 − 2c σ ) 2 when ρ → −1. This is similar to what happens in the Ellis wormhole [46][47][48][49][50]. Note that solution (85) includes the solution (83) in the limit c σ = 0. Again the Ricci scalar of this solution is zero, but the Kretschmann scalar has a more complicated expression. In four dimensions we have The Kretschmann scalar becomes singular at some radial point if |1 − 1/c σ | < 1. This is actually the case for higher dimensions too. We can prevent the geometry from becoming sick if we choose c σ < 1/2. Note that, if this expression holds, then ρσ > 0 everywhere. We have also assumed that κ > 0, but different sign choices result in equivalent solutions, with some differences in the global signs of the parameters. The spinor functions that support the wormhole geometry are To analyze the behavior of the matter content we look at the density j 0 net , which in this case looks like with the relations (79), (80), and κ = (d − 2)/2. The density j 0 net is in general not zero in any of the asymptotic regions ρ = ±1: From these expressions we can see that the field density, like the Kretschmann scalar, also diverges on the left side when c σ = 1/2, but as long as c σ < 1/2 the expression is regular and positive everywhere. Let us calculate the mass of the wormhole using the Komar integral. We need a time-like Killing vector normalized to one at infinity. For ρ → 1 we can use the due to our normalization. The mass calculated using the dual form ξ + t of the Killing K + t , using the expression is For the other side we cannot use ∂ t , because which is generally not equal to one. Instead, let us define the Killing vector Hence we have lim ρ→−1 g(K − t , K − t ) = 1. We also have to be careful with regard to the vielbein we use. Because e r points towards spatial infinity on the ρ > 0 side, but it points towards the wormhole on the ρ < 0 side. So we have to reorient the vielbein for the mass calculation as well, changing from ω r to −ω r . This introduces a minus sign in the star operator, and the expression of the Komar integral is  with ξ − t the dual form of K − t . This relation indicates that each side measures values for the mass of the wormhole with contrary signs, but also with different absolute values. Note that the mass is finite in both sides, as long as c σ < 1/2 is satisfied.
If we insist in having a positive value for M + , then 0 ≤ c σ < 1/2. In this case the value of M − is always negative. Note also that for the singular solution with c σ = 1/2, the M − mass diverges but the M + mass reaches its maximum possible value.
Note that we can find solutions with the opposite behaviour of the ρ < 0 and ρ > 0 sides if we choose a different sign of κ. In that case the M + mass could diverge while the M − mass would always remain finite.
Keeping with κ positive, notice that for c σ < 1/2 so c ∆ determines how much the ratio of the field amplitude in the asymptotic region differs from the ratio of the absolute values of the masses measured in these regions.
To discuss the role of c ∆ further, let us look at the phase ν(ρ). From equation (68) we have (96) In Figure 2 we show a plot of the phase function ν in d = 4 for various values of c ∆ . As we can see it changes from the boundary value ν(ρ = −1) = π to the boundary value ν(ρ = +1) = 0 at a position determined by c ∆ . For c ∆ = 0 this phase jump happens at ρ = 0, for c ∆ < 0 in the region ρ < 0 and for c ∆ > 0 in the region ρ > 0. This change of phase suggests a relation between the spinors as defined by an asymptotically flat observer on the right side or on the left side. Following [29] let us discuss the spinor field on both asymptotic flat regions. For this let us introduce the observers Alice and Bob. Alice will be the observer living in the asymptotically flat region ρ → 1. Quantities like M + and κ have been defined in the frame and with the vielbein of Alice. The observer Bob lives in ρ → −1. Bob differs from Alice by his choice of time normalization (the temporal Killing vector Bob has to use is given by equation (93) and the radial coordinate r ′ defined on the part ρ < 0 by the equation −ρ = N (r ′ ). Changing to Bobs frame we have the Dirac equation where we used e r . This shows us that so the chirality of the field in Bob frame is the opposite to the chirality in Alice frame. Thus we can write any of the spinors in the incoherent superposition of fields given by the sign choices ǫ in the asymptotic regions as where we should keep in mind that the first part of that product is written in the frame of Alice and the second in the frame of Bob. So one could think of the spinors in the asymptotic regions as anticorrelated entangled pairs completely similiar to the discussion in [29]. The change of chirality is a result of the different orientation of observers at each side of the wormhole, and it has nothing to do with the particular solution we have obtained, so this should be a general feature of Dirac fields in the geometry of wormholes. The only difference with other models is that this wormhole solution is a full back-reacting solution of the Einstein-Dirac equations, so the wormhole geometry is supported by the Dirac spinors.
Let us now explicitly calculate the position of the throat. For this we look at a slice of constant time and keep all angles constant except one. Then we embed this into a two-dimensional metric in a three-dimensional space using a function z(r) The position of the throat can be calculated from dr/dz| r throat = 0, so from Finally, let us comment on the energy conditions of this wormhole. The stress-energy tensor for this particular solution is So there are directions for which the energy conditions (null and weak) are violated. One could also easily see that from the fact that, due to a massless field T a a = 0 but also T tt = 0, so there must be directions for which T ab ξ a ξ b < 0 for time-like ξ a . This is a generic feature of wormholes supported by exotic matter. In addition, in this case we have seen that the density of the Dirac fields does not decay at infinity, meaning that it is also more difficult to interpret the Dirac fields of these solutions.
For instance, these wormholes we have obtained are described in practice by three parameters: µ, c σ and c ∆ (apart from the coupling constant α g ). Essentially, these three parameters are related respectively with the radius of the throat, the mass of the wormhole, and the amplitude of the Dirac fields. In principle, the amplitude can be fixed by imposing an extra normalization condition on the field. For example, we could impose the integral over the density to be equal to one, if we want to connect with the quantum (probabilistic) interpretation of the field. This is not possible in our case since the integral diverges. One could fix this parameter following other criteria, for example, by fixing the value of the density (89) at one of the asymptotical regions, or by choosing a particular relation between the mass and density ratios (95) .
A natural question is to ask if these solutions can be generalized to include massive Dirac fields m = 0 and/or frequency ω. A simple analytical solution doesn't seem to be available, and this suggests that numerical methods may be necessary for the construction of these configurations.
Finally, let us note that the properties of Dirac fields in the background of wormhole geometries have been studied before in the literature [29][30][31]. Solutions with pairs of Dirac fields in the background of a wormhole can be used as effective models describing a short nanotube bridging two different graphene layers. These models are known as graphene wormholes [22][23][24][25][26][27][28]. Such models are constructed in lower dimensions (d = 2 + 1), in the presence of a gauge field and without back-reaction.

B. Schwarzschild black hole with a divergent Dirac flux
Let us now consider the special solution for which either f ≡ 0 or g ≡ 0, also with m = ω = 0. In this case, from equations (74) we can see that the function σ is just a constant, which we choose to be one. The metric is thus simply the d-dimensional Schwarzschild-Tangherlini metric. We will nevertheless have a non-vanishing spinor field in this background. The background is a vacuum black hole, because for the above spinor field configuration the stress-energy tensor vanishes. To have a well-behaved solution at infinity we choose f = 0 in the case of κ < 0 and g = 0 in the case of κ > 0 (recall κ is fixed by the Ansatz construction to κ = ǫ κ d−2 2 ). Define for this with c 1 ∈ R ≥0 being a constant. This solution can actually be reached in the previous solution (75), when taking in the expressions c ∆ = −ǫ κ c 0 + 2ǫ κ ln c 1 and the limit c σ = 1 (c 0 → −∞). The Dirac density of this solution is given by In Figure 3 we show the Dirac density j 0 net as a function of r for several values of the dimension, where we can see that it decays to zero at infinity but diverges at the horizon. Let us analyze this behaviour in more detail.
The asymptotic part for determining the behaviour of the field at spatial infinity is for y ≤ 1 and y ≈ 1. So for r → ∞ we have and thus Close to the horizon j 0 net diverges like From this we can conclude that the integral over the density is not finite: the quantity explodes logarithmically at the horizon. Hence this solution is actually sick at the level of the matter field content, and is not physically reasonable. This is of course expected, since some general results forbid regular black holes with Dirac fields to exist [51,52]. Nonetheless, it is interesting to see that the combination of Dirac fields conspires in such a way so that the effective stress-energy tensor vanishes completely, and hence the geometry (the metric) is not affected at all by the matter configuration. Also it is interesting to see that the problem of the solution can be explicitly tracked to the behaviour of the Dirac fields at the horizon of the black hole.

C. Light-like singularity
Another solution can be obtained when m = 0 but ω = 0. Specializing to solutions with dν dr = 0, we find that σ has to be The equation (66) for σ implies that N has to be a constant too, being With this we can integrate the equation (66) for λ, implying a simple expression for this function, with L ∈ R a length scale defined from the integration constant.
Lastly the equation (66) implies two algebraic relations: This can always be solved for d ≥ 3. With this and thus equation (114) implies the frequency has to be fine-tuned Gathering all the relations, the solution is sin ν = ±1. This means that This fixes ǫ κ ǫ ν = +1. The density of the field goes to zero at infinity, although it is divergent for r → 0. We can simplify the expressions a bit by choosing the length The metric has a very simple form: The Kretschmann scalar for this solution in four dimensions d = 4 is while the curvature scalar vanishes, R = 0. The above singularity is light-like. To see this explicitly, we can calculate the Penrose diagram. The in-and out-going null geodesics obey the equation with r ′ , t ′ ∈ (−1, 1) brings the line element to where v = t + d−1 d−3 ln r and w = t − d−1 d−3 ln r. The singularity is at r → 0, meaning it is at t ′ + r ′ = −1 and t ′ − r ′ = 1. This set defines a light-like surface. In Figure  4 we show the Penrose diagram of metric (125) in the t ′ , r ′ coordinates and showing a few curves with constant t and r. All the radial geodesics of a massive particle begin and end at the singularity r = 0.

VI. CONCLUSIONS
In this paper we have studied the properties of configurations with a collection of Dirac fields, chosen in such a way that the total stress-energy tensor of the matter content is compatible with the spherical symmetry of the metric.
In sections II and III we have given a detailed explanation on how the collection of 2 ⌊ d−2 2 ⌋ Dirac fields can be chosen in order to achieve this symmetry. In order to do this, one proposes the standard separable Ansatz for each indivitual spinor. The radial and temporal dependence of the spinor is assumed to be equal for all fields. We make use of the known solutions for the angular part of the spinor, in particular when the (d − 2)-sphere is factorized as a tower of lower dimensional spheres. Then it is possible to show that, for certain values of the angular momentum of the field, |κ| = (d − 2)/2, the angular dependence of each field combines with the rest in a way so that the total stress-energy tensor is compatible with the spherical symmetry of the background. In order to have static metrics, we have to impose the vanishing of the radial current.
Making use of this Ansatz, we simplified the action and field equations in section IV. With this simplified system we constructed several simple analytical solutions. In section V A we obtained a family of wormholes supported by the Dirac fields. These wormholes connect two asymptotically flat regions, they can have positive, zero or negative mass, and their geometry and matter content are regular everywhere. However the density does not decay to zero at any of the asymptotically flat regions. We analyzed the relation between several quantities at each asymptotically flat region, in particular, the mass, which changes sign and value at each side, and the chirality of the Dirac fields, which also changes sign.
In section V B we found that the Schwarzschild metric can in fact satisfy the Einstein-Dirac field equations with a non-trivial solution for the fermionic fields. The catch is that the matter content becomes sick at the horizon of the black hole (the density diverges), and the fields cannot be normalized in the standard way.
In section V C we have also obtained a light-like naked singularity, where matter fields and geometry become simultaneously sick.
For a study of more general solutions, it is likely that numerical methods are required in order to construct the solutions. For instance, the Ansatz we have presented here in detail was used to obtain numerically regular selfgravitating solitons [19].
A possible continuation of this work is to investigate numerically if it is possible to obtain similar wormhole solutions with finite values of the Dirac mass and/or frequency, or with additional matter fields (like a gauge field, similar to what is used in the graphene wormholes). It would be interesting to explore if in these cases regu-lar solutions still exist, and if the Dirac density can be asymptotically zero. On the other hand, these non-trivial asymptotics of the Dirac fields in the case of the wormhole, suggest that maybe allowing the metric to have other asymptotical behaviour (i.e., allowing for a negative cosmological constant to have an asymptotically AdS space-time) could help regularize the integral of the Dirac density.