Multi-state Dirac stars

In this paper, we construct the multi-state Dirac stars (MSDSs) consisting of two pairs of Dirac fields. The two pairs of Dirac fields are in the ground state and the first excited state, respectively. Each pair consists of two fields with opposite spins, ensuring spherical symmetry of the system. We discuss the solutions of the MSDSs under synchronized and nonsynchronized frequencies. By varying the mass $\tilde{\mu}_1$ of the excited state Dirac field and the frequency $\tilde{\omega}_0$ of the ground state Dirac field, we obtain different types of solutions, including single-branch and double-branch solutions. These two types of solutions do not smoothly transition into each other as the parameters $\tilde{\mu}_1$ and $\tilde{\omega}_0$ continuously change, but undergo a sudden transition when $\tilde{\mu}_1$ ($\tilde{\omega}_0$) is greater than or less than the threshold value of $0.7694$ ($0.733$). Furthermore, we analyze the characteristics of the various MSDSs solutions and analyze the relationship between the ADM mass $M$ of the MSDSs and the synchronized and nonsynchronized frequencies. Subsequently, we calculate the binding energy $E_B$ of the MSDSs and discuss the stability of the solutions. Finally, we discuss the feasibility of simulating the dark matter halos using MSDSs.


I. INTRODUCTION
Recently, there has been rapid development in the field of gravitational wave astronomy, which has provided us with new insights into compact objects such as the black holes (BHs) and the neutron stars (NSs) [1][2][3].The advancements in gravitational wave detection technology have also made it possible to search for exotic compact objects (ECOs) similar to the BHs.One prominent class of ECOs is the bosonic stars, which are particle-like configurations of massive scalar fields [4][5][6][7][8][9] or vector fields [10][11][12][13][14][15] that form under their own gravitational attraction.The repulsive force balancing gravity is provided by the Heisenberg uncertainty principle.The bosonic stars provide a promising framework for studying compact objects, and certain models of the bosonic stars can mimic the BHs [16][17][18][19][20]. Additionally, the bosonic stars are also considered candidates for dark matter [21][22][23][24][25][26].
However, particle-like configurations can also be formed by spin-1/2 fermion fields.For non-gravitational cases, attempts to construct particle-like solutions for the Dirac equation were made as early as the 1930s by Ivanenko [27].Subsequent studies have also been conducted in this regard [28][29][30][31].However, it was not until 1970 that the exact numerical solutions for such particle-like configurations were first studied by Soler [32].When gravitational interactions are considered, numerical calculations become more challenging.In 1999, Finster et al. constructed the exact numerical solutions for the Einstein-Dirac system, which couples spinor fields with Einstein's gravity, for the first time [33].These particle-like configurations, formed by spin-1/2 fermions under their own gravitational attraction, are known as the Dirac stars.Subsequently, research on the Dirac stars has been extended to include charged [34] and gauge field [35] additions, and the existence of the Dirac star solutions has been proven [36,37].Recently, the rotational Dirac star solutions [38] and their charged counterparts [39] have been provided for the first time by Herdeiro et al. Some comparative studies between the bosonic stars and the Dirac stars have been conducted in [40,41].Additionally, various interesting studies on the Einstein-Dirac system have been carried out [42][43][44][45][46][47][48][49][50][51].
In 2010, Bernal et al. constructed the multi-state boson stars (MSBSs) composed of two complex scalar fields in their ground state and the first excited state and analyzed the stability of the solutions [52].Subsequently, the MSBSs were extended to include rotation [53] and self-interactions [54].It is possible that the Dirac field, under its own gravitational attraction, can also form multi-state configurations.In this work, we numerically solve the Einstein-Dirac system and construct spherically symmetric multi-state Dirac stars (MSDSs), where two coexisting states of the Dirac field are present.
The paper is organized as follows.In Sec.II, we introduce the Einstein-Dirac system, which couples four-dimensional Einstein's gravity with two sets of Dirac fields.In Sec.III, we investigate the boundary conditions of the MSDSs.In Sec.IV, we present the numerical results and analyze the solutions of the MSDSs under synchronized and nonsynchronized frequencies.We also discuss the binding energy of the solutions and the problem of the galactic halos.We conclude in Sec.V.

II. THE MODEL SETUP
We consider a system composed of multiple matter fields that are minimally coupled to Einstein's gravity.The matter fields consist of two pairs of Dirac spinor fields, with each pair containing two spinors of opposite spin.This arrangement ensures that the system exhibits spherical symmetry.One pair is in the ground state, while the other pair is in the first excited state.For such a system, the action is given by: where R is the Ricci scalar, G is the gravitational constant,and L 0 and L 1 are the Lagrangians of the spinor fields in the ground state and first excited state, respectively, where Ψ

(k)
n are spinors with mass µ n and n radial nodes, and the index k = 1, 2, corresponds to spinors with opposite spin.The variations of the action (1) with respect to the metric and the field functions yield the Einstein equations and the Dirac equation: where T 0 αβ and T 1 αβ are the energy-momentum tensors of the two sets of spinor fields, Furthermore, it can be seen from Eq. ( 2) and Eq. ( 3) that the Lagrangians of the spinor fields are invariant under a global U (1) transformation Ψ n , where α is an arbitrary constant.As a result, the system possesses a conserved current: Integrating the timelike component of the conserved current over a spacelike hypersurface S yields the Noether charge: To construct spherically symmetric solutions, we choose the metric to be of the following form: where N (r) = 1 − 2m(r)/r.The two pairs of Dirac fields are given by [40]: where the index n also represents the number of radial nodes, and ω n is the frequency of the Dirac field with n radial nodes.We only consider the cases of n = 0, 1 in this paper.
Substituting the above ansatz into the field equations (4-5) yields the following system of ordinary differential equations: And the Noether charges of the system are: III. BOUNDARY CONDITIONS To solve the system of ordinary differential equations obtained in the previous section, it is necessary to impose appropriate boundary conditions.First, for a regular, asymptotically flat spacetime, the metric function should satisfy the following boundary conditions: where the ADM mass M and σ 0 are unknown constants.In addition, the matter field vanishes at infinity: Expanding equations (13)(14) near the origin, we obtain that the field functions satisfy the following condition at the origin:

IV. NUMERICAL RESULTS
In order to facilitate numerical calculations, we employ the following dimensionless quantities: where M P l = 1/ √ G is the Planck mass.For any physical quantity A, we denote the dimensionless quantity under the conditions of ρ = 1/µ 0 and ρ = 1/µ 1 as Ã and A, respectively.
To facilitate computation, we define the radial coordinate x as follows: where the radial coordinate r ∈ [0, ∞), so x ∈ [0, 1].We utilize the finite element method to numerically solve the system of differential equations.The integration region 0 ≤ x ≤ 1 is discretized into 1000 grid points.The Newton-Raphson method is employed as our iterative approach.To ensure the accuracy of the computational results, we enforce a relative error criterion of less than 10 −5 .
To ensure the accuracy of our numerical calculations, it is crucial to verify the numerical precision by validating physical constraints [57,58], in addition to employing the aforementioned numerical analysis methods.In this study, we examined the equivalence between the asymptotic mass and the Komar mass of the numerical solution, and the results consistently showed a discrepancy of less than 10 −5 between these two quantities.
We denote the Dirac stars in the ground state and the first excited state as D 0 and D 1 , respectively, and the multi-state Dirac stars as D 0 D 1 (or MSDSs).The representation of the gamma matrices and the choice of the tetrad in the Dirac equation are the same as in [42]. A.

Synchronized frequency
Through the analysis of the numerical calculations, we observed that the solution of the MSDSs under synchronized frequency depends on the ratio of the masses of the excited state and ground state Dirac fields: µ 1 /µ 0 , which is the dimensionless mass μ1 of the first excited state Dirac field.By varying the value of the mass μ1 , various MSDSs solutions can be obtained.Based on the number of branches in the obtained MSDSs solutions, we categorize them into single-branch and double-branch solutions.For 0.7694 ≤ μ1 < 1, the MSDSs solution corresponds to a single-branch solution, whereas for 0.7573 ≤ μ1 < 0.7694, the MSDSs solution corresponds to a double-branch solution.In the following discussion, we will delve into the characteristics of these two types of solutions.

Single-branch
We first discuss the more general single-branch solution in the multi-field system [42,43,[53][54][55].The characteristic change in the radial profile of the matter field forming the MSDSs as the synchronized frequency ω continuously varies is shown in Fig. 1.The field functions depicted in the figure were obtained under the condition of a fixed mass μ1 = 0.898.The  Furthermore, when the mass μ1 of the excited state Dirac field is small, both endpoints of the orange line are located on the first branch of the black and blue dashed lines.As the mass μ1 decreases to 0.801, the intersection point of the orange line and the blue dashed line is located at the inflection point between the first and second branches of the blue dashed line.As the mass μ1 further decreases to 0.771, the intersection point of the orange line and the black dashed line is located at the inflection point between the first and second branches of the black dashed line.When μ1 decreases to the minimum frequency at which the single-branch solution can exist, 0.7694, the orange line and the black dashed line exhibit a "tangent" form.The two plots at the bottom of Fig. 2 intuitively illustrate the variation of the orange line endpoints.It should be noted that the horizontal axis of the lower-right plot represents ω, not ω, in order to show the changing trend of the intersection points of the orange line and the blue dashed line.

Double-Branch
In addition to the single-branch solution described in the preceding section, the solutions of MSDSs exhibit two branches when the synchronized frequency ω is sufficiently low.In the following, we will first discuss the variation of the matter field functions with respect to the synchronized frequency ω for the double-branch solution.As shown in Fig. 3, we demonstrate the relationship between the radial profile of the matter fields that constitute the MSDSs and the synchronized frequency ω under the condition of a fixed mass μ1 = 0.898.For the first branch field functions of the MSDSs displayed in the left column of the figure, as the synchronized frequency increases, the peak values of the ground state Dirac field functions f0 and g0 also increase, while the peak values of the excited state Dirac field functions f1 and g1 gradually decrease.For the second branch field functions displayed in the right column, as the synchronized frequency decreases, the peak values of the ground state field functions f0 and g0 gradually decrease, while the peak values of the excited state field functions f1 and g1 gradually increase.It is worth noting that the ground state Dirac field disappears at the minimum synchronized frequency of the first and second branches, while the excited state field always exists.
Next, we analyze the characteristics of the ADM mass M of the double-branch solutions   In Fig. 4, the MSDSs solutions exhibit a double-branch structure when 1/µ 0 (μ 1 ) is less than the threshold value of 0.7964.This change in solutions occurs abruptly, and the single-branch solutions do not gradually extend into a second branch as 1/µ 0 decreases.
Extensive numerical calculations were conducted, and even when 1/µ 0 is gradually reduced by an order of magnitude of 10 −5 , the single-branch solutions fail to transition continuously to the double-branch solutions.We refer to this phenomenon as bifurcation.It is worth noting that this phenomenon has not been observed in previous studies of multi-field soliton models [42,43].
The two branches of the double-branch solution are very close to each other, but their intersections with the blue dashed line are clearly different.As 1/µ 0 decreases, the synchronized frequency range of the two branches of the MSDSs gradually decreases.For all double-branch solutions, the intersections of the two branches with the blue dashed line are always in the second branch of the blue dashed line.When 1/µ 0 decreases to 0.7573, the second branch solution is very close to some of the solutions in the first branch, to the extent that the two branches almost coincide.It is noteworthy that at this point, the first branch of the MSDSs appears to be "tangential" to the blue dashed line, which is a similar feature to the single-branch solution of the MSDSs mentioned earlier.

B. Nonsynchronized frequency
In this section, we discuss the nonsynchronized frequency solutions of the MSDSs.To analyze the influence of the parameters on the numerical solutions, we set the masses of the  When the frequency ω0 of the ground state Dirac field is close to 1, the range of the frequency ω1 corresponding to the obtained single-branch solution is very narrow.As the frequency ω0 gradually decreases, the range of the frequency ω1 increases gradually.
Additionally, for a fixed frequency ω0 , as the frequency ω1 increases, the ADM mass of the system decreases gradually.When the ADM mass reaches its maximum value, the frequency ω1 reaches its minimum value, and the ground state Dirac field disappears, causing the MSDSs to degenerate into D 1 .Conversely, when the ADM mass reaches its minimum value, this minimum value is the same as the ADM mass of D 0 when the frequency ω0 takes the values indicated in each plot.This is because when the ADM mass of the MSDSs reaches its minimum value, the frequency ω1 reaches its maximum value, and the excited state Dirac field disappears, resulting in the system degenerating into D 0 .In other words, the minimum ADM mass of MSDSs depends on the frequency ω0 of the ground state Dirac field.

Double-Branch
Next, we discuss the double-branch solutions of the MSDSs under nonsynchronized frequency.In the single-branch solution discussed earlier, the minimum value of the frequency ω0 of the ground state Dirac field is 0.733.Since there is no solution with a frequency less than 0.733 for D 0 , the MSDSs cannot degenerate into D 0 when the frequency ω0 is less than 0.733, and thus the single-branch solution cannot be obtained.Through a series of numer-ical calculations, we found that when the frequency ω0 of the ground state field satisfies 0.6971 ≤ ω0 < 0.733, the double-branch solutions of the MSDSs are obtained.Fig. 7 shows the relationship between the radial profile of the matter fields that constitute the MSDSs and the nonsynchronized frequency ω1 .The left column shows the field functions on the first branch.As the frequency ω1 increases, the peak values of the ground state Dirac field functions f0 and g0 gradually increase, while the changes in the excited state field functions f1 and g1 are relatively small.The right column shows the field functions on the second branch.As the frequency ω1 decreases, the peak values of the ground state Dirac field functions f0 and g0 gradually decrease, while the excited state field functions f1 and g1 gradually increase.It can be observed that when the frequency ω1 of the first and second branch solutions reaches its minimum value, the ground state Dirac field disappears, while the excited state Dirac field does not disappear for any ω1 .
The relationship between the ADM mass M of the nonsynchronized frequency doublebranch solutions of the MSDSs and the frequency ω1 is shown in Fig. 8.The black and blue dashed lines represent the ground state and first excited state solutions of the Dirac star (D 0 and D 1 ), and the orange line represents the double-branch solution of the MSDSs (D 0 D 1 ).Similar to the case of synchronized frequencies mentioned earlier, a bifurcation occurs when the frequency ω0 is below the threshold of 0.733, transforming the singlebranch solution into the double-branch solution.As the frequency ω0 of the ground state Dirac field decreases, the range of existence of the two branches of the double-branch solution gradually decreases with respect to the frequency ω1 .For a fixed ground state field frequency ω0 , when the frequency ω1 reaches its maximum value, the MSDSs do not degenerate into D 0 but transition to another new branch.As the frequency ω1 decreases further, the ADM mass of the system gradually increases, and eventually, the MSDSs transform into D 1 .This characteristic change in the system's mass is also reflected in the profile of the field functions presented in Fig. 7, where the disappearance of the ground state Dirac field functions on the two branches when the frequency ω1 reaches its minimum indicates the transformation of the MSDSs into D 1 .

C. Binding energy
After obtaining various different solutions for the MSDSs, we will analyze the stability of the system from the perspective of binding energy.Consider a MSDS with ADM mass M , where the Noether charge for the ground state Dirac field is denoted as Q 0 , and the Noether charge for the first excited state Dirac field is denoted as Q 1 .The binding energy E B of the system can be expressed as: where the coefficient 2 outside the parentheses on the right-hand side of the equation arises from the fact that in a spherically symmetric MSDS, both the ground state and the first excited state of the Dirac field have two components.
1/μ 0 =0.99We first analyze the binding energy of the MSDSs under synchronized frequency.Fig. 9 shows the relationship between the binding energy E B and the synchronized frequency ω of the MSDSs.The left plot represents the single-branch solutions of the MSDSs.When 1/µ 0 > 0.863, for a fixed value of µ 0 , the binding energy E B of the MSDSs monotonically increases with the synchronized frequency ω, and the binding energy is always less than zero.
When 1/µ 0 ≤ 0.863, for a given µ 0 , the binding energy initially decreases and then increases as the synchronized frequency increases.When 1/µ 0 becomes sufficiently small, the solutions become unstable (e.g., the red curve).Therefore, when 1/µ 0 is sufficiently large, i.e., when the masses µ 0 and µ 1 of the ground state and excited state Dirac fields are sufficiently close, the MSDSs are more stable.The right plot represents the double-branch solutions of the MSDSs.As 1/µ 0 increases, the minimum value of the binding energy gradually decreases, but the binding energies of these double-branch solutions are all greater than zero, indicating that the solutions are unstable.
Next, we consider the binding energy of the MSDSs under nonsynchronized frequency.
Fig. 10 shows the relationship between the binding energy E B and the frequency ω1 for the nonsynchronized frequency solutions of the MSDSs.The left plot represents the singlebranch solutions, where for a fixed ω0 , the binding energy E B monotonically increases with the frequency ω1 and remains negative when the frequency ω0 is sufficiently large.However, for small values of ω0 (e.g., ω0 = 0.733), the MSDSs can undergo a transition from a stable solution to an unstable one as the frequency ω1 increases.The right plot represents the double-branch solutions, where the negative binding energy solutions appear when the frequency ω0 is sufficiently large (e.g., ω0 = 0.732).As the frequency ω0 decreases, the stable solutions in the double-branch solutions gradually disappear, and all solutions eventually become unstable.Therefore, in the case of nonsynchronized frequency, stable solutions exist within the double-branch solutions.

D. Galactic halos as MSDSs
The velocity of stars orbiting the central core of a galaxy remains constant over a large range of distances starting from the galactic center.This phenomenon may be attributed to the presence of a dark matter halo in the outer regions of the galaxy.By using boson stars to simulate the dark matter halo, it is possible to obtain results that are consistent with real observational data [21,52,56].In the following, we will analyze the feasibility of simulating the dark matter halo using MSDSs by computing the velocities of test particles orbiting around them.Considering timelike circular geodesics on the equatorial plane, the rotational velocity of the test particles is given by [52,56]: Substituting Equ.(10) into the expression, we obtain Next, we analyze the rotational curves of the MSDSs.As shown in Fig. 11, the left panel

V. CONCLUSION
In this paper, we investigate the Einstein-Dirac system and construct spherically symmetric multi-state Dirac stars, wherein two coexisting states of the Dirac field are present.
We discuss the field functions, ADM mass, and binding energy of the solutions for the multi-state Dirac star under synchronized and nonsynchronized frequency conditions.Furthermore, we analyze the feasibility of considering the multi-state Dirac star as a candidate for dark matter halos.
In the case of synchronized frequency, we explore different solutions for the MSDSs by varying the ratio of masses between the ground state and excited state Dirac fields (µ 0 /µ 1 ).
Based on the number of solution branches, we classify the obtained numerical results into single-branch solutions and double-branch solutions.For single-branch solutions, the peak values of the ground state and excited state Dirac field functions exhibit monotonic behavior as the synchronized frequency changes.The ADM mass monotonically decreases with increasing synchronized frequency, and at the minimum and maximum synchronized frequencies, the MSDSs degenerate into the first excited state Dirac stars and the ground state Dirac stars, respectively.When μ1 (or 1/µ 0 ) is below the threshold of 0.7694, the single-branch solution undergoes a transition into a double-branch solution.In the case of double-branch solutions, the excited state Dirac field functions persist while the ground state field functions disappear at the minimum synchronized frequency on both branches, leading to the degeneration of the MSDSs into the first excited state Dirac stars.Moreover, as the mass of the ground state field decreases, the range of synchronized frequency values on the two branches of the double-branch solution gradually diminishes.
Next, for the case of nonsynchronized frequency, we set the ratio of masses between the ground state and excited state Dirac fields as µ 0 /µ 1 = 1 and obtain different solutions by varying the frequency of the ground state Dirac field.Similar to the synchronized frequency case, the nonsynchronized frequency solutions of the MSDSs can also be classified into singlebranch solutions and double-branch solutions.In the case of single-branch solutions, the variations of the ground state and excited state field functions with respect to the frequency of the excited state Dirac field exhibit monotonic behavior.The ADM mass decreases as the frequency of the excited state field increases, and its minimum value depends on the frequency of the ground state Dirac field.When ω0 is below the threshold of 0.733, the single-branch solution undergoes a transition into a double-branch solution.For the doublebranch solutions, as the frequency of the ground state field decreases, the minimum ADM mass of the MSDSs gradually increases, the range of nonsynchronized frequency values on the two branches diminishes, and the MSDSs degenerate into the first excited state Dirac stars when the nonsynchronized frequency becomes sufficiently small.
It is worth noting that the MSDSs solutions exhibit certain similarities between the synchronized and nonsynchronized frequency cases.In the synchronized frequency (nonsynchronized frequency) scenario, when μ1 (ω 0 ) falls below the threshold of 0.7694 (0.733), the single-branch solutions undergo an abrupt transition into double-branch solutions.Furthermore, regardless of whether the frequencies are synchronized or not, the single-branch solutions of the MSDSs can always degenerate to the ground state or the first excited state, while the double-branch solutions can only degenerate to the first excited state.Subsequently, we computed the binding energy of various solutions for the MSDSs.For synchronized frequency solutions, stable solutions only exist in the single-branch solutions, and when 1/µ 0 is sufficiently small, all stable solutions within the single-branch solutions vanish.For nonsynchronized frequency solutions, stable solutions exist in both single-branch and double-branch solutions.However, the double-branch solutions become unstable when the frequency of the ground state field reaches a sufficiently low value, whereas the singlebranch solutions maintain stable solutions for any frequency of the ground state field.
Finally, we computed the rotation curves of the MSDSs.It has been observed that MSDSs containing excited state matter fields exhibit a nearly flat region near the velocity peak in their rotation curves, similar to the rotation curves of multi-state boson stars discussed in [52].In our future work, we plan to construct MSDSs with a greater number of matter field nodes.It is possible that rotational curves of models with a higher number of nodes can provide a closer fit to the observed data [56].
excited state Dirac stars (D 1 ), and the orange line represents the MSDSs (D 0 D 1 ).The ADM mass of the system monotonically decreases as the synchronized frequency increases.It can be observed that the upper end of the orange line intersects with the blue dashed line, where the MSDSs degenerate into the D 1 ; the lower end of the orange line intersects with the black dashed line, where the MSDSs degenerate into the D 0 .This degeneration of the MSDSs is manifested in Fig. 1 as the disappearance of the ground state or excited state field functions.

FIG. 3 .
FIG. 3. The matter functions f0 , g0 , f1 and g1 on the first (left column) and second branches (right column) of the MSDSs solution function as functions of x and ω for μ1 = 0.7645.

FIG. 4 .
FIG. 4. The ADM mass M of the MSDSs as a function of the synchronized frequency ω for several values of 1/µ 0 .

FIG. 6 .
FIG. 6.The ADM mass M of the MSDSs as a function of the frequency ω1 for several values of ω0 .

FIG. 7 .
FIG. 7. The matter functions f0 , g0 , f1 and g1 on the first (left column) and second branches (right column) of the MSDSs solution function as functions of x and ω1 for ω0 = 0.721.

FIG. 8 .
FIG. 8.The ADM mass M of the MSDSs as a function of the frequency ω1 for several values of ω0 .

2 FIG. 9 .
FIG. 9.The binding energy E B of the MSDSs as a function of the synchronized frequency ω for several values of 1/µ 0 .

FIG. 11 .
FIG. 11.Left panel: the rotational curves of the MSDSs for several values of μ1 .Right panel: the rotational curves of the Dirac stars at different energy levels.