weight $O(|\boldsymbol{e}_{10}\cup\boldsymbol{e}_{01}|)$. Since the code is
LDPC, it means that this can be done without increasing the error weight too
much.

######  Lemma 28.

For any error $\boldsymbol{e}$, there exists an equivalent error
$\boldsymbol{e}^{\prime}$ such that $\boldsymbol{e}^{\prime}_{01}=0$,
$\boldsymbol{e}^{\prime}_{10}=0$ and

$|\boldsymbol{e}^{\prime}|\leqslant|\boldsymbol{e}|+4\Delta^{2}(|\boldsymbol{e}_{01}|+|\boldsymbol{e}_{10}|).$

###### Proof.

Consider a general error $(\boldsymbol{e}^{X},\boldsymbol{e}^{Z})$. We focus
on the $\sigma_{Z}$-part
$(\boldsymbol{e}_{AB}^{Z},\boldsymbol{e}_{00}^{Z},\boldsymbol{e}_{10}^{Z},\boldsymbol{e}_{01}^{Z},\boldsymbol{e}_{11}^{Z})$.
The $\sigma_{X}$-part is treated identically. Let $U_{0},U_{1}$ be matrices
such that

$U_{0}H_{B}=\boldsymbol{e}_{10}^{Z},\quad U_{1}h_{B}=\boldsymbol{e}_{01}^{Z}.$

We have that $|U_{0}|\leqslant\Delta|\boldsymbol{e}_{10}^{Z}|$ and
$|U_{1}|\leqslant\Delta|\boldsymbol{e}_{01}^{Z}|$. These stabilizers induce
new errors $h_{A}^{T}U_{0}+H_{A}^{T}U_{1}+U_{0}h_{B}^{T}+U_{1}H_{B}^{T}$,
which have a weight at most $2\Delta(|U_{0}|+|U_{1}|)\leqslant
2\Delta^{2}(|\boldsymbol{e}_{10}|+|\boldsymbol{e}_{01}|)$. ∎

Thanks to Lemma 28, we now focus on a $\sigma_{Z}$-error $\boldsymbol{e}$ such
that $\boldsymbol{e}_{01}=0$ and $\boldsymbol{e}_{10}=0$. Its syndrome reads:

$T_{0}=\boldsymbol{e}_{AB}h_{B}^{T}+h_{A}^{T}\boldsymbol{e}_{00},\qquad
T_{1}=\boldsymbol{e}_{AB}H_{B}^{T}+H_{A}^{T}\boldsymbol{e}_{11}.$

We can multiply these syndromes by $g_{A}$ and $G_{A}$ respectively and obtain

$g_{A}T_{0}=g_{A}\boldsymbol{e}_{AB}h_{B}^{T},\qquad
G_{A}T_{1}=G_{A}\boldsymbol{e}_{AB}H_{B}^{T},$

which is exactly the decoding problem for quantum Tanner codes. We note that a
similar strategy was already exploited in [QC21]. The decoder presented in the
previous sections can correct this error provided its weight is small enough
(Theorem 2). More precisely, if $|\boldsymbol{e}_{AB}|\leqslant\kappa
n/\Delta^{4}$, then the decoder outputs an equivalent error
$\hat{\boldsymbol{e}}_{AB}$ for the quantum Tanner code. In other words, there
exists a stabilizer element associated with $(U_{0},U_{1})$ such that

$\hat{\boldsymbol{e}}_{AB}=\boldsymbol{e}_{AB}+h_{A}^{T}U_{0}^{\prime}G_{B}+H_{A}^{T}U_{1}^{\prime}g_{B},$

as described in the previous subsection.

Let us introduce a matrix $n_{A}$ indexed by $C_{0}\times A$ such that
$n_{A}h_{A}^{T}=\mathbbm{1}_{C_{0}}$. We assume here that $h_{A}$ is a full-
rank matrix. We also denote $N_{A}:=n_{A}D_{A}^{-1}$, which implies
$N_{A}H_{A}^{T}=\mathbbm{1}_{C_{1}}$.

The decoder for the lifted product code will simply return the triple
$(\hat{\boldsymbol{e}}_{AB},\hat{\boldsymbol{e}}_{00},\hat{\boldsymbol{e}}_{11})$
with

$\hat{\boldsymbol{e}}_{00}:=n_{A}(T_{0}+\hat{\boldsymbol{e}}_{AB}h_{B}^{T}),\quad\hat{\boldsymbol{e}}_{11}:=N_{A}(T_{1}+\hat{\boldsymbol{e}}_{AB}H_{B}^{T}).$

Straightforward manipulations show that

$\displaystyle\hat{\boldsymbol{e}}_{00}$
$\displaystyle=n_{A}((\boldsymbol{e}_{AB}h_{B}^{T}+h_{A}^{T}\boldsymbol{e}_{00})+(\boldsymbol{e}_{AB}+h_{A}^{T}U_{0}^{\prime}G_{B}+H_{A}^{T}U_{1}^{\prime}g_{B})h_{B}^{T})$
$\displaystyle=\boldsymbol{e}_{00}+U_{0}^{\prime}G_{B}h_{B}^{T}$

and

$\displaystyle\hat{\boldsymbol{e}}_{11}$
$\displaystyle=N_{A}((\boldsymbol{e}_{AB}H_{B}^{T}+H_{A}^{T}\boldsymbol{e}_{11})+(\boldsymbol{e}_{AB}+h_{A}^{T}U_{0}^{\prime}G_{B}+H_{A}^{T}U_{1}^{\prime}g_{B})H_{B}^{T})$
$\displaystyle=\boldsymbol{e}_{11}+U_{1}^{\prime}g_{B}H_{B}^{T}$

which means that the true error and the error returned by the decoder differ
by

$\underbrace{h_{A}^{T}(U_{0}^{\prime}G_{B})+H_{A}^{T}(U_{1}^{\prime}g_{B})}_{A\times
B}+\underbrace{(U_{0}^{\prime}G_{B})h_{B}^{T}}_{C_{0}\times
D_{0}}+\underbrace{(U_{1}^{\prime}g_{B})H_{B}^{T}}_{C_{1}\times D_{1}}.$

In other words, they differ by the stabilizer element
$(U_{0},U_{1})=(U_{0}^{\prime}G_{B},U_{1}^{\prime}g_{B})$ of the lifted
product code, meaning that they are equivalent and that the decoder succeeded.

We have therefore proved the following:

######  Theorem 29.

There exists a constant $\kappa$, depending only on $\delta$, a lower bound
for the minimum distances of both component codes $C_{A}$ and $C_{B}$, such
that for large enough fixed $\Delta$, the above decoding algorithm corrects
all error patterns of weight less than $\kappa n/\Delta^{6}$ for the lifted
product code of length $n=|Q|$.

## References

  * [AE15] Dorit Aharonov and Lior Eldar.  Quantum locally testable codes.  SIAM Journal on Computing, 44(5):1230–1262, 2015.
  * [BE21] Nikolas P Breuckmann and Jens N Eberhardt.  Balanced product quantum codes.  IEEE Transactions on Information Theory, 67(10):6653–6674, 2021\.
  * [BK22] Nouédyn Baspin and Anirudh Krishna.  Connectivity constrains quantum codes.  Quantum, 6:711, May 2022.
  * [BPT10] Sergey Bravyi, David Poulin, and Barbara Terhal.  Tradeoffs for Reliable Quantum Information Storage in 2D Systems.  Phys. Rev. Lett., 104:050503, Feb 2010.
  * [BT09] Sergey Bravyi and Barbara Terhal.  A no-go theorem for a two-dimensional self-correcting quantum memory based on stabilizer codes.  New Journal of Physics, 11(4):043029, apr 2009.
  * [CS96] Robert Calderbank and Peter W Shor.  Good quantum error-correcting codes exist.  Physical Review A, 54:1098–1105, Aug 1996.
  * [Del13] Nicolas Delfosse.  Tradeoffs for reliable quantum information storage in surface codes and color codes.  In 2013 IEEE International Symposium on Information Theory, pages 917–921, 2013.
  * [DEL+21] Irit Dinur, Shai Evra, Ron Livne, Alexander Lubotzky, and Shahar Mozes.  Locally Testable Codes with constant rate, distance, and locality.  arXiv preprint arXiv:2111.04808, 2021.
  * [DHLV22] Irit Dinur, Min-Hsiu Hsieh, Ting-Chun Lin, and Thomas Vidick.  Good quantum LDPC codes with linear time decoders.  arXiv preprint arXiv:2206.07750, 2022.
  * [DN21] Nicolas Delfosse and Naomi H. Nickerson.  Almost-linear time decoding algorithm for topological codes.  Quantum, 5:595, December 2021.
  * [Edm65] Jack Edmonds.  Paths, trees, and flowers.  Canadian Journal of Mathematics, 17:449–467, 1965.
  * [EKZ20] Shai Evra, Tali Kaufman, and Gilles Zémor.  Decodable quantum LDPC codes beyond the square root distance barrier using high dimensional expanders.  In 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), pages 218–227, 2020.
  * [FGL18] Omar Fawzi, Antoine Grospellier, and Anthony Leverrier.  Constant overhead quantum fault-tolerance with quantum expander codes.  In 2018 IEEE 59th Annual Symposium on Foundations of Computer Science (FOCS), pages 743–754, 2018.
  * [Gal62] Robert Gallager.  Low-density parity-check codes.  IRE Transactions on Information Theory, 8(1):21–28, 1962.
  * [GL14] Larry Guth and Alexander Lubotzky.  Quantum error correcting codes and 4-dimensional arithmetic hyperbolic manifolds.  Journal of Mathematical Physics, 55(8):082202, 2014.
  * [Gol21] Oded Goldreich.  On the Locally Testable Code of Dinur et al.( 2021).  In Electron. Colloquium Comput. Complex., volume 28, page 175, 2021\.
  * [Got97] Daniel Gottesman.  Stabilizer codes and quantum error correction.  PhD thesis, California Institute of Technology, 1997.
  * [Got14] Daniel Gottesman.  Fault-tolerant quantum computation with constant overhead.  Quantum Information & Computation, 14(15-16):1338–1372, 2014.
  * [GPT22] Shouzhen Gu, Christopher A Pattison, and Eugene Tang.  An efficient decoder for a linear distance quantum LDPC code.  arXiv preprint arXiv:2206.06557, 2022.
  * [Gur10] Venkatesan Guruswami.  Expander codes and their decoding.  https://www.cs.cmu.edu/~venkatg/teaching/codingtheory/notes/notes8.pdf, 2010\.
  * [Has14] Matthew B Hastings.  Decoding in hyperbolic spaces: quantum LDPC codes with linear rate and efficient error correction.  Quantum Information & Computation, 14(13-14):1187–1202, 2014.
  * [Has17] Matthew B Hastings.  Quantum codes from high-dimensional manifolds.  In 8th Innovations in Theoretical Computer Science Conference (ITCS 2017). Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik, 2017.
  * [HHO20] Matthew B Hastings, Jeongwan Haah, and Ryan O’Donnell.  Fiber Bundle Codes: Breaking the $N^{1/2}\mathrm{polylog}(N)$ Barrier for Quantum LDPC Codes.  arXiv preprint arXiv:2009.03921, 2020.
  * [HL22] Max Hopkins and Ting-Chun Lin.  Explicit Lower Bounds Against $\Omega(n)$-Rounds of Sum-of-Squares.  arXiv preprint arXiv:2204.11469, 2022.
  * [HLW06] Shlomo Hoory, Nathan Linial, and Avi Wigderson.  Expander graphs and their applications.  Bulletin of the American Mathematical Society, 43(4):439–561, 2006\.
  * [Kit03] Alexei Yu Kitaev.  Fault-tolerant quantum computation by anyons.  Annals of Physics, 303(1):2–30, 2003.
  * [KP22] Gleb Kalachev and Pavel Panteleev.  Two-sided robustly testable codes.  arXiv preprint arXiv:2206.09973, 2022.
  * [KT21] Tali Kaufman and Ran J Tessler.  New cosystolic expanders from tensors imply explicit Quantum LDPC codes with $\Omega(\sqrt{n}\log^{k}n)$ distance.  In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, pages 1317–1329, 2021.
  * [LH22] Ting-Chun Lin and Min-Hsiu Hsieh.  Good quantum LDPC codes with linear time decoder from lossless expanders.  arXiv preprint arXiv:2203.03581, 2022.
  * [LLZ22] Anthony Leverrier, Vivien Londe, and Gilles Zémor.  Towards local testability for quantum coding.  Quantum, 6:661, February 2022.
  * [LTZ15] Anthony Leverrier, Jean-Pierre Tillich, and Gilles Zémor.  Quantum expander codes.  In Foundations of Computer Science (FOCS), 2015 IEEE 56th Annual Symposium on, pages 810–824. IEEE, 2015.
  * [LZ22a] Anthony Leverrier and Gilles Zémor.  A parallel decoder for good quantum LDPC codes.  arXiv preprint arXiv:2208.05537, 2022.
  * [LZ22b] Anthony Leverrier and Gilles Zémor.  Quantum Tanner codes.  arXiv preprint arXiv:2202.13641, 2022.
  * [PK20] Pavel Panteleev and Gleb Kalachev.  Quantum LDPC Codes with Almost Linear Minimum Distance.  arXiv preprint arXiv:2012.04068, 2020.
  * [PK21a] Pavel Panteleev and Gleb Kalachev.  Asymptotically good quantum and locally testable classical LDPC codes.  arXiv preprint arXiv:2111.03654, 2021.
  * [PK21b] Pavel Panteleev and Gleb Kalachev.  Degenerate Quantum LDPC Codes With Good Finite Length Performance.  Quantum, 5:585, November 2021.
  * [QC21] Armanda O Quintavalle and Earl T Campbell.  Lifting decoders for classical codes to decoders for quantum codes.  arXiv preprint arXiv:2105.02370, 2021.
  * [SS96] Michael Sipser and Daniel A. Spielman.  Expander codes.  IEEE Transactions on Information Theory, 42(6):1710–1722, 1996\.
  * [Ste96] Andrew M. Steane.  Error correcting codes in quantum theory.  Physical Review Letters, 77:793–797, Jul 1996.
  * [Tan81] R. Tanner.  A recursive approach to low complexity codes.  IEEE Transactions on Information Theory, 27(5):533–547, 1981.
  * [TZ14] Jean-Pierre Tillich and Gilles Zémor.  Quantum LDPC codes with positive rate and minimum distance proportional to the square root of the blocklength.  IEEE Transactions on Information Theory, 60(2):1193–1202, 2014\.


Fairness-Aware Hybrid Precoding for mmWave NOMA Unicast/Multicast Transmissions in Industrial IoT

Luis F. Abanto-Leon, and Gek Hong (Allyson) Sim
Secure Mobile Networking (SEEMOO) Lab, Technische Universität Darmstadt, Germany {labanto<EMAIL_ADDRESS>
Journal of  Class Files, Vol. 14, No. 8, Jan 2019
Shell et al.: Bare Demo of IEEEtran.cls for IEEE Communications Society Journals

fairness, hybrid precoding, unicast, multicast, non-orthogonal multiple access, industrial IoT, mmWave.

[1]
Z. Ding, P. Fan, and H. V. Poor, “Random Beamforming in Millimeter-Wave NOMA
Networks,” IEEE Access, pp. 7667–7681, February 2017.

[2]
Z. Ding, L. Dai, R. Schober, and H. V. Poor, “NOMA Meets Finite Resolution
Analog Beamforming in Massive MIMO and Millimeter-Wave Networks,”
IEEE Communications Letters, vol. 21, no. 8, pp. 1879–1882, August

[3]
3GPP, “3GPP TR 36.859: Study on Downlink Multiuser Superposition Transmission
(MUST) for LTE (Release 13),” 2015.

[4]
H. Nikopour and H. Baligh, “Sparse Code Multiple Access,” in IEEE
PIMRC, September 2013, pp. 332–336.

[5]
A. Loch, C. Cano, G. H. Sim, A. Asadi, and X. Vilajosana, “A Channel
Measurement Campaign for mmWave Communication in Industrial Settings,”
CoRR, vol. abs/1903.10502, 2019. [Online]. Available:

[6]
L. F. Abanto-Leon, M. Hollick, and G. H. Sim, “Hybrid Precoding for
Multi-Group Multicasting in mmWave Systems,” in IEEE GLOBECOM,
December 2019.

[7]
L. F. Abanto-Leon and G. H. Sim, “Learning-based Max-Min Fair Hybrid
Precoding for mmWave Multicasting,” in IEEE ICC, June 2020.

[8]
P. P. Bergmans and T. M. Cover, “Cooperative Broadcasting,” IEEE
Transactions on Information Theory, vol. 20, no. 3, pp. 317–324, May 1974.

[9]
L. Zhang, W. Li, Y. Wu, X. Wang, S. Park, H. M. Kim, J. Lee, P. Angueira, and
J. Montalban, “Layered-Division-Multiplexing: Theory and Practice,”
IEEE Transactions on Broadcasting, vol. 62, pp. 216–232, March 2016.

[10]
D. Kim, F. Khan, C. V. Rensburg, Z. Pi, and S. Yoon, “Superposition of
Broadcast and Unicast in Wireless Cellular Systems,” IEEE
Communications Magazine, vol. 46, no. 7, pp. 110–117, July 2008.

[11]
W. Hao, M. Zeng, Z. Chu, and S. Yang, “Energy-Efficient Power Allocation in
Millimeter Wave Massive MIMO With Non-Orthogonal Multiple Access,”
IEEE Wireless Communications Letters, vol. 6, no. 6, pp. 782–785,
December 2017.

[12]
J. Zhao, O. Simeone, D. Gunduz, and D. Gomez-Barquero, “Non-Orthogonal
Unicast and Broadcast Transmission via Joint Beamforming and LDM in Cellular
Networks,” in IEEE GLOBECOM, December 2016, pp. 1–6.

[13]
J. Wang, H. Xu, B. Zhu, L. Fan, and A. Zhou, “Hybrid Beamforming Design for
mmWave Joint Unicast and Multicast Transmission,” IEEE Communications
Letters, vol. 22, no. 10, pp. 2012–2015, October 2018.

[14]
E. Chen, M. Tao, and Y. Liu, “Joint Base Station Clustering and Beamforming
for Non-Orthogonal Multicast and Unicast Transmission With Backhaul
Constraints,” IEEE Transactions on Wireless Communications, vol. 17,
no. 9, pp. 6265–6279, September 2018.

[15]
Y. Mao, B. Clerckx, and V. O. K. Li, “Rate-Splitting for Multi-Antenna
Non-Orthogonal Unicast and Multicast Transmission,” in IEEE SPAWC,
June 2018, pp. 1–5.

[16]
——, “Rate-Splitting for Multi-Antenna Non-Orthogonal Unicast and
Multicast Transmission: Spectral and Energy Efficiency Analysis,”
IEEE Transactions on Communications, vol. 67, no. 12, pp. 8754–8770,
December 2019.

[17]
T. E. Bogale, L. B. Le, A. Haghighat, and L. Vandendorpe, “On the Number of
RF Chains and Phase Shifters, and Scheduling Design With Hybrid
Analog–Digital Beamforming,” IEEE Transactions on Wireless
Communications, vol. 15, no. 5, pp. 3311–3326, May 2016.

[18]
A. Hjorungnes and D. Gesbert, “Complex-Valued Matrix Differentiation:
Techniques and Key Results,” IEEE Transactions on Signal Processing,
vol. 55, no. 6, pp. 2740–2746, June 2007.


# Rapidly rotating Dirac stars

Vladimir Dzhunushaliev<EMAIL_ADDRESS>Department of Theoretical
and Nuclear Physics, Al-Farabi Kazakh National University, Almaty 050040,
Kazakhstan  Institute of Nuclear Physics, Almaty 050032, Kazakhstan
Academician J. Jeenbaev Institute of Physics of the NAS of the Kyrgyz
Republic, 265 a, Chui Street, Bishkek 720071, Kyrgyzstan  Vladimir Folomeev
<EMAIL_ADDRESS>Institute of Nuclear Physics, Almaty 050032, Kazakhstan
Academician J. Jeenbaev Institute of Physics of the NAS of the Kyrgyz
Republic, 265 a, Chui Street, Bishkek 720071, Kyrgyzstan  Nassurlla Burtebayev
Institute of Nuclear Physics, Almaty 050032, Kazakhstan

###### Abstract

Within general relativity, we construct sequences of rapidly rotating Dirac
stars consisting of a spinor fluid described by an effective equation of
state. We find the physically relevant domain of stable configurations and
calculate their principal characteristics which are completely determined by
the central density of the spinor fluid, the mass of the nonlinear spinor
field, and the velocity of rotation. It is demonstrated that for a certain
choice of the spinor field mass, the main physical characteristics of the
Dirac stars are close to those that are typical of rotating neutron stars.

Spinor fluid, effective equation of state, static and rotating Dirac stars

###### pacs:

04.40.Dg, 04.40.–b

##  I Introduction

Strongly gravitating configurations supported by various types of matter have
been the object of vigorous investigations in different formulations of the
problem. The physical nature of such objects may vary considerably, and is
determined both by a type of matter supporting them and by a theory of gravity
within which their modeling is performed. A typical object of investigations
are neutron stars – compact configurations possessing characteristic masses of
the order of 1-2 solar masses and radii 10-20 km [1]. In studying such
objects, the following aspects determining their physical characteristics are
of great importance: (i) the choice of equation of state of matter creating
such configurations; (ii) the inclusion of rotation and the presence of strong
magnetic fields; (iii) a theory of gravity used in modeling the stars. All
these aspects are intensively studied in the literature by starting, first,
with various physically reasonable assumptions about the nature of neutron
stars and, second, with matching theoretical models to the current
observational data.

However, it is not impossible that other types of compact strongly gravitating
objects do exist in nature. These can be, in particular, hypothetical boson
stars supported by scalar (spin-0) fields. Their investigations are carried
out in the same directions as those of neutron stars. In doing so, depending
on the specific parameters of the scalar field, the physical characteristics
of such objects may vary over a wide range, from those that are typical of
atoms to those which characterize galaxies [2, 3].

Apart from boson stars, a consideration of configurations consisting of fields
with nonzero spins is also possible. Such hypothetical systems may be
supported by massive vector (spin-1) fields (Proca stars [4, 5, 6, 7]) and
massless vector fields (Yang-Mills systems [8, 9]). In the case of fractional-
spin fields it is also possible to obtain various spherically symmetric
systems consisting of gravitating linear [10, 11, 5] and nonlinear spin-1/2
fields [12, 13, 14]. Nonlinear spinor fields are also employed in studying
cylindrically symmetric solutions [15], wormhole solutions [16], and various
cosmological problems (see Refs. [17, 18, 19] and references inside). In turn,
the inclusion of rotation in the system enables one to get axially symmetric
solutions for configurations supported by linear spin-1 and spin-1/2 fields
[20].

In our previous works, we have considered static gravitating configurations
supported by two nonlinear spin-1/2 fields. In doing so, the cases both of
purely Einstein-Dirac systems [21] and of configurations containing also U(1)
Maxwell and Proca fields [22], or SU(2) Yang-Mills and Proca fields [23] have
been investigated. The remarkable feature of all such systems is that, for
some values of the spinor field coupling constant $\lambda$, one can obtain
configurations whose masses are comparable to the Chandrasekhar mass and
effective radii are of the order of kilometers. Such objects may already be
regarded as self-gravitating Dirac stars which are prevented from collapse
under their own gravitational fields due to the Heisenberg uncertainty
principle.

In the next step, in Ref. [24], we have used the fact that for physically
reasonable values of the coupling constant $\lambda$, it is possible to
replace the nonlinear spinor fields by a hydrodynamic fluid described by an
effective equation of state (EOS). This allowed us, within Newtonian gravity,
to consider static compact objects and to study a dynamical process of their
creation through the gravitational collapse of such a fluid.

The next natural step in studying such type of configurations is the inclusion
of rotation in the system; this will be done in the present paper. To do this,
as in Ref. [24], we will employ the effective EOS that describes well
configurations supported by the nonlinear spinor field in the physically
interesting range of system parameters and whose usage does not require a
consideration of a much more complicated problem when the nonlinear spinor
field is directly involved. We will work within Einstein’s general relativity
and study rapidly rotating configurations which are in their own strong
gravitational field.

The paper is organized as follows. In Sec. II, we introduce the effective
equation of state for a spinor fluid. In Sec. III, within general relativity,
we write down the general field equations sourced by the spinor fluid. In Sec.
IV, we solve these equations numerically and obtain spherically and axially
symmetric solutions describing compact static and rapidly rotating
configurations. Finally, in Sec. V, we summarize the results obtained.

##  II Equation of state for a spinor fluid

In Refs. [21, 22, 23] we have studied gravitating systems with a doublet of a
nonlinear spin-1/2 field $\psi$ described by the Lagrangian

$L_{\text{sp}}=\frac{i\hbar
c}{2}\left(\bar{\psi}\gamma^{\mu}\psi_{;\mu}-\bar{\psi}_{;\mu}\gamma^{\mu}\psi\right)-\mu
c^{2}\bar{\psi}\psi+\frac{\lambda}{2}\left(\bar{\psi}\psi\right)^{2},$ (1)

where $\lambda$ is the coupling constant and $\mu$ is the mass of the spinor
field.

In Ref. [21], it was demonstrated that the introduction of the dimensionless
coupling constant
$\bar{\lambda}\sim\left(M_{\text{Pl}}/\mu\right)^{2}\lambda$, where
$M_{\text{Pl}}$ is the Planck mass, in the limit of $|\bar{\lambda}|\gg 1$,
permits one to get configurations whose sizes and masses are comparable to
those typical of neutron stars. It was also shown there that such limiting
configurations can be described by an effective EOS relating the pressure $p$
and the energy density $\varepsilon$ as

$p=\frac{\varepsilon_{0}}{9}\left(1+3\frac{\varepsilon}{\varepsilon_{0}}-\sqrt{1+6\frac{\varepsilon}{\varepsilon_{0}}}\right),$
(2)

where $\varepsilon_{0}=\mu c^{2}/\lambda_{c}^{3}\equiv\mu^{4}c^{5}/\hbar^{3}$
can be regarded as a characteristic energy density of the system [here
$\lambda_{c}=\hbar/(\mu c)$ is the constant, which need not be associated with
the Compton length since we consider a classical theory]. A configuration
described by such an effective EOS can then be regarded as consisting of a
fluid, which will be further called a spinor fluid.

It is worth emphasizing here that in the general case of configurations
supported by spinor fields, there are two different effective pressures – the
radial, $p_{r}$, and tangential, $p_{t}$. Correspondingly, the system is in
general anisotropic [21]. However, making use of the results of Ref. [21], it
can be shown that in the limit $|\bar{\lambda}|\gg 1$ these two pressures
become equal, i.e., the system becomes isotropic. For this reason, we consider
here only one effective pressure $p=p_{r}=p_{t}$ appearing in Eq. (2).

##  III General equations

Within the framework of Einstein’s general relativity, we consider compact
gravitating configurations consisting of a perfect spinor fluid described by
the EOS (2). The corresponding total action for such a system can be
represented in the form [the metric signature is $(+,-,-,-)$]

$S_{\text{tot}}=-\frac{c^{3}}{16\pi G}\int
d^{4}x\sqrt{-g}R+S_{\text{matter}},$ (3)

where $G$ is the Newtonian gravitational constant, $R$ is the scalar
curvature, and $S_{\text{matter}}$ denotes the action of the spinor fluid.

Varying the action (3) with respect to the metric, we derive the Einstein
equations

$E_{\mu}^{\nu}\equiv R_{\mu}^{\nu}-\frac{1}{2}\delta_{\mu}^{\nu}R-\frac{8\pi
G}{c^{4}}T_{\mu}^{\nu}=0$ (4)

with the energy-momentum tensor
$T_{\mu}^{\nu}=\left(\varepsilon+p\right)U_{\mu}U^{\nu}-\delta_{\mu}^{\nu}p,$
where $U_{\mu}$ represents the four-velocity of the fluid.

In order to obtain rotating configurations, we use the Lewis-Papapetrou line
element for a stationary, axially symmetric spacetime with two Killing vector
fields $\xi=\partial_{t},\eta=\partial_{\varphi}$ in a system of spherical
coordinates $\\{t,r,\theta,\varphi\\}$. In these coordinates the metric is
independent of $t$ and $\varphi$, and can be expressed using isotropic
coordinates as [25]

$ds^{2}=f(dx^{0})^{2}-\frac{l}{f}\left[g\left(dr^{2}+r^{2}d\theta^{2}\right)+r^{2}\sin^{2}\theta\left(d\varphi-\frac{\omega}{r}dx^{0}\right)^{2}\right],$
(5)

where the metric functions $f,l,g$, and $\omega$ depend solely on $r$ and
$\theta$, and $x^{0}=ct$. The $z$-axis ($\theta=0$) represents the symmetry
axis of the system. Asymptotically (as $r\to\infty$), the functions $f,l,g\to
1$ and $\omega\to 0$; i.e., the spacetime approaches Minkowski spacetime.

We consider here uniform (or rigid) rotation of the spinor fluid; such an
assumption is well justified, e.g., for most neutron stars [1]. In this case
the four-velocity is of the form

$U^{\mu}=\left(u,0,0,\frac{\Omega}{c}u\right),$ (6)

where $\Omega$ is a constant parameter denoting the angular velocity of the
system.

For numerical calculations, it is convenient to introduce the dimensionless
variables

$x=\sqrt{8\pi}\frac{r}{R_{L}},\quad\left(\tilde{\varepsilon},\tilde{p}\right)=\frac{\left(\varepsilon,p\right)}{\varepsilon_{0}},\quad\tilde{\Omega}=\frac{R_{L}}{\sqrt{8\pi}c}\Omega\quad\text{with}\quad
R_{L}=\lambda_{c}\frac{M_{\text{Pl}}}{\mu}.$ (7)

The parameter $R_{L}$ is the characteristic radius obtained by Landau in
considering compact configurations supported by an ultrarelativistic
degenerate Fermi gas and modeled within the framework of Newtonian gravity.

Next, taking into account the condition for the four-velocity
$U^{\mu}U_{\mu}=1$ and making use of the metric (5) and the expression (6),
the velocity $u$ can be expressed in terms of the metric functions $f,l$, and
$\omega$ as follows [written already in terms of the dimensionless variables
(7)]:

$u^{2}=\frac{1}{f\left[1-x^{2}f^{-2}l\sin^{2}\theta\left(\tilde{\Omega}-\omega/x\right)^{2}\right]}.$
(8)

In turn, the conservation law $\nabla_{\mu}T^{\mu\nu}=0$ yields the
differential equations

$\frac{\partial_{r}p}{\varepsilon+p}=\frac{\partial_{r}u}{u},\quad\frac{\partial_{\theta}p}{\varepsilon+p}=\frac{\partial_{\theta}u}{u}.$

Substituting in these equations the EOS (2) and integrating them, one can find
the expression for the energy density in terms of the four-velocity,

$\tilde{\varepsilon}=\frac{1}{32}\left(3c_{0}^{4}u^{4}-4c_{0}^{2}u^{2}-4\right),$

where $c_{0}$ is an integration constant. In turn, substituting here the
expression (8), one can get

$\tilde{\varepsilon}=\frac{1}{24}\left\\{A\frac{f_{c}f}{f^{2}-l\sin^{2}\theta\left(\omega-\tilde{\Omega}x\right)^{2}}\left[A\frac{f_{c}f}{f^{2}-l\sin^{2}\theta\left(\omega-\tilde{\Omega}x\right)^{2}}-2\right]-3\right\\}.$
(9)

For convenience, we have introduced here a new constant $A$ defined by the
relation $c_{0}^{2}=2/3Af_{c}$, where $f_{c}$ is the central value of the
metric function $f$.

The field equations coming from the Einstein equations (4) can be represented
in the following form [written in terms of the dimensionless variables (7)]:

$\displaystyle
f_{xx}+\frac{l\sin^{2}\theta}{f}\left(\frac{2}{x}\omega-\omega_{x}\right)\omega_{x}+\left(\frac{2}{x}+\frac{l_{x}}{2l}\right)f_{x}-\frac{f_{x}^{2}}{f}+\frac{1}{x^{2}}\left(f_{\theta\theta}-\frac{l\sin^{2}\theta}{f}\omega_{\theta}^{2}+\frac{f_{\theta}l_{\theta}}{2l}-\frac{f_{\theta}^{2}}{f}+\cot\theta
f_{\theta}-\frac{l\sin^{2}\theta\omega^{2}}{f}\right)$
$\displaystyle=\frac{gl}{f}\left\\{\left(\tilde{\varepsilon}+\tilde{p}\right)\left[f^{2}+l\sin^{2}\theta\left(\omega-\tilde{\Omega}x\right)^{2}\right]u^{2}+2f\tilde{p}\right\\},$
(10) $\displaystyle
l_{xx}-\frac{l_{x}^{2}}{2l}+\frac{3l_{x}}{x}+\frac{1}{x^{2}}\left(l_{\theta\theta}-\frac{l_{\theta}^{2}}{2l}+2\cot\theta
l_{\theta}\right)=4\frac{gl^{2}}{f}\tilde{p},$ (11) $\displaystyle
g_{xx}+\frac{3gl\sin^{2}\theta}{f^{2}}\left(\frac{\omega}{x}-\frac{\omega_{x}}{2}\right)\omega_{x}-\frac{g}{l}\left(\frac{2}{x}+\frac{l_{x}}{2l}\right)l_{x}-\frac{g_{x}^{2}}{g}+\frac{g_{x}}{x}+\frac{gf_{x}^{2}}{2f^{2}}$
$\displaystyle+\frac{1}{x^{2}}\left(g_{\theta\theta}-\frac{3gl\sin^{2}\theta\omega_{\theta}^{2}}{2f^{2}}-\frac{gl_{\theta}^{2}}{2l^{2}}-\frac{2\cot\theta
gl_{\theta}}{l}-\frac{g_{\theta}^{2}}{g}+\frac{gf_{\theta}^{2}}{2f^{2}}-\frac{3gl\sin^{2}\theta\omega^{2}}{2f^{2}}\right)$
$\displaystyle=\frac{2g^{2}l}{f^{2}}\left[u^{2}l\sin^{2}\theta\left(\tilde{\varepsilon}+\tilde{p}\right)\left(\omega-\tilde{\Omega}x\right)^{2}-f\tilde{p}\right],$
(12)
$\displaystyle\omega_{xx}+\left(\frac{2}{x}-\frac{2f_{x}}{f}+\frac{3l_{x}}{2l}\right)\omega_{x}-\frac{3\omega
l_{x}}{2xl}+\frac{2\omega
f_{x}}{xf}+\frac{1}{x^{2}}\left[\omega_{\theta\theta}+\left(3\cot\theta-\frac{2f_{\theta}}{f}+\frac{3l_{\theta}}{2l}\right)\omega_{\theta}-2\omega\right]=2u^{2}gl\left(\tilde{\varepsilon}+\tilde{p}\right)\left(\omega-\tilde{\Omega}x\right).$
(13)

(Here the lower indices denote differentiation with respect to the
corresponding coordinate.) These equations are the following combinations of
the components of the Einstein equations (4):
$\left(E_{t}^{t}-E_{r}^{r}-E_{\theta}^{\theta}-E_{\varphi}^{\varphi}+2\,\omega/rE_{\varphi}^{t}\right)=0$,
$\left(E_{r}^{r}+E_{\theta}^{\theta}\right)=0$,
$\left(E_{\varphi}^{\varphi}-E_{r}^{r}-E_{\theta}^{\theta}-\omega/rE_{\varphi}^{t}\right)=0$,
and $E_{\varphi}^{t}=0$.

Apart from the Einstein equations (10)-(13), which are elliptic partial
differential equations, there are also two more gravitational equations
$E^{r}_{\theta}=0$ and $\left(E^{r}_{r}-E^{\theta}_{\theta}\right)=0$, whose
structure is not of Laplace form and which can be regarded as “constraints”.
According to arguments proposed in Ref. [26], to find a self-consistent
solution to Eqs. (10)-(13), it is enough to check that the constraint
$E^{r}_{\theta}=0$ is satisfied on all boundaries; this in turn implies that
it is satisfied throughout the region of integration. Then the constraint
$\left(E^{r}_{r}-E^{\theta}_{\theta}\right)=0$ must only be imposed at a
single point to ensure that it is fulfilled throughout the region of
integration as well. Notice that, as applied to Dirac stars supported by a
linear spinor field, such an approach to obtaining self-consistent solutions
has been used in Ref. [20]. In the present paper, in order to find self-
consistent numerical solutions, we will also check whether these constraints
are fulfilled.

##  IV Numerical solutions

In this section we obtain sequences of static and rotating solutions to Eqs.
(10)-(13) for the whole physically relevant range of system parameters.

###  IV.1 Boundary conditions

We will seek globally regular, asymptotically flat solutions describing
configurations with a finite mass. For such solutions, we impose appropriate
boundary conditions for the metric functions at the origin ($x=0$), at
infinity ($x\to\infty$), on the positive $z$-axis ($\theta=0$), and, making
use of the reflection symmetry with respect to $\theta\to\pi-\theta$, in the
equatorial plane ($\theta=\pi/2$). So we require

$\displaystyle\begin{split}&\left.\frac{\partial f}{\partial
x}\right|_{x=0}=\left.\frac{\partial l}{\partial
x}\right|_{x=0}=0,\left.g\right|_{x=0}=1,\left.\omega\right|_{x=0}=0;\\\
&\left.f\right|_{x=\infty}=\left.g\right|_{x=\infty}=\left.l\right|_{x=\infty}=1,\left.\omega\right|_{x=\infty}=0;\\\
&\left.\frac{\partial
f}{\partial\theta}\right|_{\theta=0,\pi}=\left.\frac{\partial
l}{\partial\theta}\right|_{\theta=0,\pi}=\left.\frac{\partial\omega}{\partial\theta}\right|_{\theta=0,\pi}=0,\left.g\right|_{\theta=0,\pi}=1;\\\
&\left.\frac{\partial
f}{\partial\theta}\right|_{\theta=\pi/2}=\left.\frac{\partial
g}{\partial\theta}\right|_{\theta=\pi/2}=\left.\frac{\partial
l}{\partial\theta}\right|_{\theta=\pi/2}=\left.\frac{\partial\omega}{\partial\theta}\right|_{\theta=\pi/2}=0.\end{split}$
(14)

Notice that the elementary flatness condition implies that on the symmetry
axis the function $g|_{\theta=0,\pi}=1$ [25]; this ensures the absence of a
conical singularity.

###  IV.2 Asymptotic behavior

Since we seek asymptotically flat solutions, we require the metric to approach
the Minkowski metric at spatial infinity, i.e., asymptotically, $f,g,l\to 1$
and $\omega\to 0$. For our purposes, let us write down expressions describing
an asymptotic behavior of the metric functions

$f\approx 1-\frac{2GM}{c^{2}r}+{\cal
O}\left(r^{-2}\right),\quad\omega\approx\frac{2GJ}{c^{3}r^{2}}+{\cal
O}\left(r^{-3}\right).$

These expressions contain the total mass of the star $M$ and the angular
momentum $J$, which can be represented in the form [25]

$M=\frac{c^{2}}{2G}\lim_{r\to\infty}r^{2}\partial_{r}f,\quad
J=\frac{c^{3}}{2G}\lim_{r\to\infty}r^{2}\omega.$ (15)

###  IV.3 Center and surface of the star

Let us now write out expressions for central values of the energy density and
pressure of the star. Using Eqs. (2) and (9) and taking into account that at
the centre $\left.\omega\right|_{x=0}=0$ [see the boundary conditions (14)],
one can obtain

$\tilde{\varepsilon}_{c}=\frac{1}{24}\left[A^{2}-2A-3\right],\quad\tilde{p}_{c}=\frac{1}{9}\left(1+3\tilde{\varepsilon}_{c}-\sqrt{1+6\tilde{\varepsilon}_{c}}\right).$

By choosing the value of the arbitrary constant $A=1+2\sqrt{7}$, we have
$\tilde{\varepsilon}_{c}=1$. In turn, given different values of $A$, one can
construct families of solutions describing configurations with different
physical characteristics (see below).

The boundary of the star is defined as the surface with radius $r_{b}(\theta)$
where $\varepsilon\to 0$ (correspondingly, $p\to 0$ as well). To construct an
embedding diagram that depicts the intrinsic geometry of the stellar surface
[i.e., the slice of spacetime with $t=\text{const.}$ and $r=r_{b}(\theta)$],
one can embed the stellar surface in a flat three-dimensional space [27]. The
metric of the stellar surface induced by the spacetime metric (5) is

$ds_{b}^{2}=\frac{l_{b}}{f_{b}}\left\\{g_{b}\left[\left(\partial_{\theta}r_{b}\right)^{2}+r_{b}^{2}\right]d\theta^{2}+r_{b}^{2}\sin^{2}\theta
d\varphi^{2}\right\\},$

where the index $b$ denotes the boundary of the star, on which the metric
functions are given [for example, $f_{b}\equiv
f\left(r_{b}(\theta),\theta\right)$]. Then, changing to cylindrical
coordinates $\\{\rho,z,\varphi\\}$ for the flat space, one can obtain the
following expressions [27]:

$\rho(\theta)=\sqrt{\frac{l_{b}}{f_{b}}}r_{b}\sin\theta,\quad
z(\theta)=\int_{\theta}^{\pi/2}d\theta^{\prime}\left\\{\frac{g_{b}l_{b}}{f_{b}}\left[\left(\frac{dr_{b}}{d\theta^{\prime}}\right)^{2}+r_{b}^{2}\right]-\left(\frac{d\rho}{d\theta^{\prime}}\right)^{2}\right\\}^{1/2}.$

Making use of them, the equatorial and polar radii of the embedded surface are
then given by

$\displaystyle
R_{e}=\rho\left(\theta=\frac{\pi}{2}\right)=\left.\sqrt{\frac{l_{b}}{f_{b}}}r_{b}\right|_{\theta=\pi/2},$
(16) $\displaystyle
R_{p}=z(\theta=0)=\int_{0}^{\pi/2}d\theta^{\prime}\left\\{\frac{g_{b}l_{b}}{f_{b}}\left[\left(\frac{dr_{b}}{d\theta^{\prime}}\right)^{2}+r_{b}^{2}\right]-\left(\frac{d\rho}{d\theta^{\prime}}\right)^{2}\right\\}^{1/2}.$
(17)

Consider now the question of a limiting angular velocity of rotation (or the
Keplerian limit) of the star under consideration. This limit is reached when
the star rotates so fast that the angular velocity $\Omega$ of the fluid
approaches the angular velocity $\Omega_{p}$ of a free particle located at the
equator of the stellar surface. As a result, a fluid element cannot already be
confined on the surface of the star and this ultimately leads to losing mass.

The corresponding expression for $\Omega_{p}$ can be found from the geodesic
equation for the angular velocity of the particle [28]

$\left(\frac{\Omega_{p}}{c}-\frac{\omega_{b}}{r}\right)^{2}-2a_{p}\left(\frac{\Omega_{p}}{c}-\frac{\omega_{b}}{r}\right)+b_{p}=0,$
(18)

where $\omega_{b}\equiv\omega(r_{b},\pi/2)$ and

$a_{p}=\left.\frac{fl\left(r\partial_{r}\omega-\omega\right)}{r\left[f\left(2l+r\partial_{r}l\right)-rl\partial_{r}f\right]}\right|_{r_{b},\pi/2},\quad
b_{p}=\left.-\frac{f^{2}\partial_{r}f}{r\left[f\left(2l+r\partial_{r}l\right)-rl\partial_{r}f\right]}\right|_{r_{b},\pi/2}.$

Solving Eq. (18) for $\Omega_{p}$, one can obtain

$\frac{\Omega_{p}}{c}=\frac{\omega_{b}}{r}+a_{p}+\sqrt{a_{p}^{2}-b_{p}}.$

The Keplerian angular velocity $\Omega_{K}$ is defined as
$\Omega_{K}=\Omega_{p}$. A necessary condition for the existence of a
stationary rotating configuration implies that the equatorial velocity of
stellar matter is smaller than the Keplerian velocity of a particle moving
along the $\varphi$-direction in the equatorial plane. The Keplerian velocity
corresponds to the Keplerian angular velocity $\Omega_{K}$, also called the
mass-shedding angular velocity.

###  IV.4 Numerical approach

The set of four coupled nonlinear elliptic partial differential equations
(10)-(13) for the unknown functions $f,g,l$, and $\omega$ will be solved
numerically subject to the boundary conditions (14). Keeping in mind that we
will seek even parity solutions which are symmetrical about the equatorial
plane $\theta=\pi/2$, numerical computations will be carried out only in the
region $0\leq\theta\leq\pi/2$. In doing so, for numerical calculations, it is
convenient to introduce new compactified coordinate

$\bar{x}=\frac{x}{1+x},$ (19)

the use of which permits one to map the infinite region $[0,\infty)$ to the
finite interval $[0,1]$.

All results of numerical calculations for axially symmetric systems exhibited
below have been obtained using the package FIDISOL [29] where the numerical
procedure based on the Newton-Raphson method is employed. This method provides
an iterative procedure for getting an exact solution starting from an
approximate solution (an initial guess). As such initial guess, we use the
solutions for nonrotating (i.e., static) spherically symmetric configurations
which can be obtained using Eqs. (10)-(13) with $g=1$ and $\Omega,\omega=0$.
In this case these equations reduce to the set of two ordinary differential
equations for the metric functions $f$ and $l$ which have been solved using
the NDSolve routine from Mathematica.

In turn, solutions to the partial differential equations (10)-(13) are sought
on a grid of $401\times 101$ points, covering the region of integration
$0\leq\bar{x}\leq 1$ [given by the compactified coordinate from Eq. (19)] and
$0\leq\theta\leq\pi/2$. In performing calculations, we have kept track of the
behavior of the energy density at every point in space $(r,\theta)$ in order
to determine the position of the boundary of the star $r_{b}(\theta)$ where
$\varepsilon\to 0$. The case of $\varepsilon<0$ corresponds to the fact that
the meshpoint is outside the star; in this case the right-hand sides of Eqs.
(10)-(13) are taken to be zero (i.e., the vacuum Einstein equations are
solved).

###  IV.5 Static and rotating solutions

Using the approach described in the previous subsection, let us obtain
sequences of solutions describing Dirac stars. The strategy of their obtaining
is as follows. First, as seed configurations, we obtain spherically symmetric,
nonrotating configurations, starting the solution in the neighborhood of the
center where the metric functions behave as

$f\approx f_{c}+\frac{1}{2}f_{2}x^{2}+\ldots,\quad l\approx
l_{c}+\frac{1}{2}l_{2}x^{2}+\ldots$

In doing so, to ensure the asymptotic flatness of the spacetime, it is
necessary to take strictly defined central values $f_{c}$ and $l_{c}$
(eigenvalue problem) which ensure $f,l\to 1$ at infinity. Then, by varying the
integration constant $A$ appearing in the expression for the energy density
(9), one can find the family of static configurations with $\omega=\Omega=0$
possessing different masses and radii of the spinor fluid, which are
parameterized by $A$.

Figure 1:  Left panel: the mass-radius relation for Dirac stars in the
physically relevant domain for different values of the dimensionless angular
velocity $\tilde{\Omega}$ (designated by the numbers near the curves). For
$\mu=0.3~{}\text{GeV}$, the value $\tilde{\Omega}=0.01$ corresponds to a
frequency of $f\approx 90\,\text{Hz}$. Right panel: the mass-central energy
density relation for the systems of the left panel [the energy density is
expressed in units of $\rho_{0}$ from Eq. (21)].

Next, using the static systems as an initial guess, we gradually increase the
velocity of rotation ($\Omega\neq 0$) and keep track of changes of
characteristics of the stars under consideration. The results of calculations
are given in Fig. 1, where the left panel shows the mass-radius relation for
the systems with different velocities of rotation. To plot these dependences,
we have used the following dimensional quantities [obtained using Eqs. (7) and
(15)]:

$r_{b}(\theta)\approx
0.48\frac{\bar{x}_{b}}{1-\bar{x}_{b}}\left(\frac{\text{GeV}}{\mu}\right)^{2}\,\text{km},\quad
M\approx 0.1623\lim_{\bar{x}\to 1}\bar{x}^{2}\frac{\partial
f}{\partial\bar{x}}M_{\odot}\left(\frac{\text{GeV}}{\mu}\right)^{2}.$ (20)

To find the equatorial radius $R_{e}$, this $r_{b}(\theta)$ is substituted in
Eq. (16). In turn, the spin frequency of the configurations under
consideration is

$f\equiv\frac{\Omega}{2\pi}\approx
10^{5}\tilde{\Omega}\left(\frac{\mu}{\text{GeV}}\right)^{2}\,\text{Hz}.$

The domain of Dirac stars shown in Fig. 1 is bounded according to the
following physical reasons:

  * (i)

The left boundary of the domain corresponds to the set of static stars whose
mass increases monotonically with decreasing the radius up to its maximum.
This maximum corresponds to a point that separates stable configurations
(located to the right of the maximum) and unstable ones (located to the left
of the maximum, not shown here).

  * (ii)

The right boundary of the domain corresponds to the set of configurations
rotating at the Keplerian limit, $\Omega=\Omega_{K}$. On this boundary the
mass also increases monotonically with decreasing the radius up to its maximum
that separates stable and unstable systems.

  * (iii)

The two aforementioned boundaries are connected at the extreme points by the
so-called secular instability line (shown by the dashed line) which is the
locus of the points corresponding to maxima of masses of the configurations
rotating with different fixed values of the angular velocity $\tilde{\Omega}$.
Configurations located to the left of this boundary are unstable (see below in
Sec. IV.6).

Technically, the process of obtaining the curve corresponding to the
configurations rotating at the Keplerian limit consists in successive
increasing the value of $\Omega$ up to the point where
$\Omega\approx\Omega_{p}$. Since the numerical calculations do not permit us
to continue the solution to the point $\Omega=\Omega_{p}$, the physical
parameters of the systems rotating at the Keplerian limit have been obtained
by their extrapolation as functions of $\Omega$ in the limit
$\Omega\to\Omega_{p}$ [28].

Notice that the above behavior of the mass-radius curves of the Dirac stars
under consideration is similar to the behavior of analogous curves typical of
neutron stars (see, e.g., Refs. [28, 30]). In particular, by choosing the mass
of the spinor field $\mu=0.3~{}\text{GeV}$, total masses and sizes of the
Dirac stars are comparable to characteristics typical of neutron stars.

Figure 2: The mass $M$ vs. the angular momentum $J$ (left panel) and the
moment of inertia $I$ vs. the mass $M$ (right panel) are shown for Dirac stars
in the physically relevant domain for different values of the dimensionless
angular velocity $\tilde{\Omega}$ (designated by the numbers near the curves).

###  IV.6 Stability

Let us now address the question of stability of the systems under
consideration. In doing so, we will use the turning-point method of Ref. [31]
which provides a simple criterion to locate points of secular instability
along a sequence of equilibrium configurations. This method does not require
the use of perturbation analysis but instead deals with a consideration of
one- or many-parameter families of solutions.

Namely, in the case of static systems one can estimate stability by
considering the behavior of a mass-central density curve. The maximum mass of
a stable configuration coincides with the first maximum of a sequence of
configurations with increasing central density; this takes place at the first
point where $\partial M/\partial\varepsilon_{c}=0$ (the turning point). As the
central density increases further (after it has passed the turning point),
there occurs instability that evolves on a secular timescale, i.e., it is not
dynamical [31].

When the rotation is turned on, apart from the central density, one more extra
parameter is involved; this can be the ratio of the coordinate radii
$r_{p}/r_{e}$, the angular velocity $\Omega$, the angular momentum $J$, the
gravitational mass $M$. In this case it is always possible to construct a set
of rotating configurations by fixing a value of the second parameter and
letting the central energy density vary in a range which is constrained by
stability limits [30]. As shown in Ref. [31], the turning-point method can be
applied to uniformly rotating stars. For a sequence of configurations with a
constant angular velocity $\Omega$, the turning point of the sequence with
increasing central density serves as a boundary between secularly stable and
unstable configurations; that is, a secular axisymmetric instability occurs at
the point where $\left.\partial
M(\varepsilon_{c},\Omega)/\partial\varepsilon_{c}\right|_{\Omega=\text{const.}}=0$.
Then a curve connecting all the maxima (i.e., the turning points obtained for
different fixed values of $\Omega$) bounds the stability domain (cf. the
secular instability line shown by the dashed curve in the left panel of Fig.
1). The point of intersection of such a limiting curve with the Keplerian
sequence corresponds to a configuration rotating with the largest possible
angular velocity.

It is worth mentioning here that, as was pointed out in Ref. [32], the turning
point criterion is only a sufficient condition for secular instability of
rotating stars, and an analysis of the dynamical stability reveals that
dynamical instability occurs for smaller central densities. Nevertheless, the
turning point method employed here, while being technically more simple,
ensures sufficient accuracy in determining a stability domain.

To illustrate the above for the systems under investigation, consider the
mass-central energy density relation. For greater clarity, instead of the
energy density $\varepsilon$, we will use the mass density $\rho$ and
introduce the characteristic mass density

$\rho_{0}\equiv\frac{\varepsilon_{0}}{c^{2}}\approx 2.33\times
10^{17}\left(\frac{\mu}{\text{GeV}}\right)^{4}\,\text{g cm}^{-3}$ (21)

[the expression for $\varepsilon_{0}$ can be found after Eq. (2)]. The
corresponding dependences of the mass of the systems under consideration on
the central density are shown in the right panel of Fig. 1. The behavior of
these dependences is qualitatively similar to the behavior of the relations
obtained for neutron stars: the central density is maximum for static Dirac
stars in the stability limit and decreases with increasing angular velocity
while moving along the secular instability line (cf. Refs. [28, 30]). The
maximally rotating stable configurations are located on the Keplerian (or
mass-shedding) sequence; here the gravitational force equals the centrifugal
force at the star’s equator, and a faster rotation would inevitably lead to
the expulsion of mass from the star.

To conclude this subsection, notice that if as an example one takes
$\mu=0.3~{}\text{GeV}$, the characteristic central density of the systems
under consideration near the maximum of the mass is of the order of
$10^{15}\,\text{g cm}^{-3}$; as in the case of masses and sizes (see the end
of the previous subsection), this is the value that is also typical of neutron
stars.

###  IV.7 Angular momentum and moment of inertia

We turn now to the angular momentum $J$ of rotating Dirac stars. In the left
panel of Fig. 2, we exhibit the dependences of the total mass $M$ on the
angular momentum for different velocities of rotation. To plot these graphs,
we have used the expression for the mass from Eq. (20) and the dimensional
quantity [obtained using Eqs. (7) and (15)]

$J\approx 5.27\times 10^{-2}\lim_{\bar{x}\to
1}\left(\frac{\bar{x}}{1-\bar{x}}\right)^{2}\omega\,\frac{GM_{\odot}^{2}}{c}\left(\frac{\text{GeV}}{\mu}\right)^{4}.$

It is seen from the figure that (i) the secular instability line forms the
upper limit for the mass, while the lower limit corresponds to the Keplerian
sequence; and (ii) the relationship between the star’s total mass and the
angular momentum for the Keplerian sequence is practically linear.

For the sake of comparison, consider the angular momentum of Kerr black holes.
For them, the reduced dimensionless angular momentum

$\frac{a}{M}=\frac{cJ}{GM^{2}}$

is limited by the maximum value attained by the extremal black holes,
$\left(a/M\right)_{\text{max}}=1$. For neutron stars supported by different
EOSs $\left(a/M\right)_{\text{max}}\approx 0.7$ [30], i.e., it is considerably
lower than the Kerr value. Also, the quantity $a/M$ varies only slightly with
the mass. In turn, in the case of Dirac stars supported by the EOS (2), we
observe a similar behavior, with the maximum value
$\left(a/M\right)_{\text{max}}\approx 0.78$.

Compute now the moment of inertia of the star, which is one of the most
relevant properties in analysis of neutron stars (pulsars). The moment of
inertia can be evaluated as

$I=\frac{J}{\Omega}.$

The right panel of Fig. 2 shows $I$ as a function of the star’s mass for
several fixed values of the angular velocity $\tilde{\Omega}$, as well as for
the Keplerian sequence. On comparing these dependences with the results
obtained for neutron stars modeled by different EOSs [30], one can observe
that there are some qualitative differences in the behavior of the curves,
related obviously to the use of another EOS in the form (2). In turn, by
choosing as an example $\mu=0.3~{}\text{GeV}$, one can see that $I\sim
10^{45}\text{g cm}^{2}$; this is a typical value for neutron stars.

Figure 3:  Left panel: the ratio of the polar to the equatorial radius,
$R_{p}/R_{e}$, as a function of the total mass $M$ of Dirac stars in the
physically relevant domain for different values of the dimensionless angular
velocity $\tilde{\Omega}$ (designated by the numbers near the curves). Right
panel: the isometric embedding of the surface of a star close to the secular
instability line (shown by the dashed curve) and of a star close to the
Keplerian limit (shown by the solid curve); for both surfaces
$\tilde{\Omega}=0.05$.  Figure 4:  Isosurfaces of constant energy density
$\varepsilon/c^{2}$ [in units of $\rho_{0}$ from Eq. (21)] for the
configurations shown in the right panel of Fig. 3. The plots are made in a
meridional plane $\varphi=\text{const.}$ spanned by the coordinates
$X=r\sin\theta$ and $Z=r\cos\theta$. The corresponding values of the energy
density are designated by the numbers near the curves.

###  IV.8 Deformation of the stars

The presence of rotation results in a deviation of the shape of configurations
from spherical symmetry due to the centrifugal forces which deform the stars.
The flattening of a spinning star can be described by the ratio of the polar
radius, $R_{p}$ to the equatorial radius, $R_{e}$, defined by Eqs. (17) and
(16), respectively. The left panel of Fig. 3 depicts this ratio. One can see a
rapid increase of the flattening as the angular velocity approaches the
Keplerian limit $\Omega_{K}$. Another effect of the centrifugal forces is
that, at a given $M$, the central density decreases with increasing $\Omega$
(see the right panel of Fig. 1).

To visualize the above effects, in the right panel of Fig. 3, we show the
geometry of the surface of the Dirac stars using the isometric embeddings of
the surface for two configurations: one for the star that is close to the
secular instability line and one for the system close to the Keplerian limit.
It is seen that for the stars rotating close to the Keplerian limit, there
occurs a cusp at the star’s surface in the equatorial plane.

One more physical quantity of interest is the distribution of the energy
density $\varepsilon$. In Fig. 4, we depict contours of constant energy
density. For illustrative purposes, we employ the same Dirac star solutions as
those used for constructing the isometric embeddings given in the right panel
of Fig. 3. It is clearly seen from the left panel of Fig. 4 that the
centrifugal forces produce the strongest deformation in the equatorial layer
of the outer part of the configuration rotating with the angular velocity
close to the Keplerian limit (shown by the shaded region).

##  V Conclusion

In the present paper, we have been continuing our investigations of the Dirac
stars began in Refs. [21, 22, 23, 24]. Our main purpose here is to study the
influence of rapid rotation on the characteristics of the Dirac stars. In
order to describe the matter of such stars, we have taken advantage of the
fact that, for some values of the spinor field coupling constant $\lambda$
appearing in the Lagrangian (1), it is possible to introduce an effective
equation of state in the form (2) that adequately describes the system under
investigation in the physically interesting range of system parameters. For
the case of uniform rotation considered here, the EOS (2) permits to integrate
analytically the differential equations coming from the conservation law and
to derive expressions for the energy density and pressure in terms of metric
functions. These expressions have been used as a source of a gravitational
field modeled within Einstein’s general relativity. We have started out with a
consideration of static equilibrium solutions describing spherically symmetric
Dirac stars. The inclusion of rotation in such systems results in the fact
that the configurations become axisymmetric and, as the velocity of rotation
approaches the Keplerian limit, considerably flattened.

Summarizing the results obtained,

  * (i)

We have found regular solutions to the Einstein equations with a source in the
form of a perfect spinor fluid describing sequences of static and uniformly
rotating Dirac stars.

  * (ii)

For such configurations, the gravitational mass, equatorial and polar radii,
angular momentum, and moment of inertia have been calculated.

  * (iii)

The mass-radius relations have been constructed, using which we have
established the domain of stability which is bounded by the Keplerian limit
and secular instability line.

  * (iv)

The typical distributions of the energy density inside the star and shapes of
the objects obtained are exhibited.

  * (v)

It is demonstrated that for the choice of the spinor field mass $\mu\sim
0.3~{}\text{GeV}$ the physical characteristics of the Dirac stars (total
masses, radii, densities, velocities of rotation, angular momentum, and moment
of inertia) are close to those that are typical of rotating neutron stars.

## Acknowledgments

We would like to thank L. Rezzolla for discussions and useful suggestions.
This research has been funded by the Science Committee of the Ministry of
Education and Science of the Republic of Kazakhstan (Grant No. BR10965191
“Complex research in nuclear and radiation physics, high-energy physics and
cosmology for development of the competitive technologies”). We are also
grateful to the Research Group Linkage Programme of the Alexander von Humboldt
Foundation for the support of this research.

## References

  * [1] P. Haensel, A. Y. Potekhin, and D. G. Yakovlev, Neutron Stars 1: Equation of State and Structure (Springer, New York, 2007).
  * [2] F. E. Schunck and E. W. Mielke, General relativistic boson stars, Classical Quantum Gravity 20, R301 (2003).
  * [3] S. L. Liebling and C. Palenzuela, Dynamical boson stars, Living Rev. Relativity 15, 6 (2012); 20, 5 (2017).
  * [4] R. Brito, V. Cardoso, C. A. R. Herdeiro, and E. Radu, Proca stars: Gravitating Bose-Einstein condensates of massive spin 1 particles, Phys. Lett. B 752, 291 (2016).
  * [5] C. A. R. Herdeiro, A. M. Pombo, and E. Radu, Asymptotically flat scalar, Dirac and Proca stars: Discrete vs. continuous families of solutions, Phys. Lett. B 773, 654 (2017).
  * [6] N. Sanchis-Gual, F. Di Giovanni, M. Zilhão, C. Herdeiro, P. Cerdá-Durán, J. A. Font, and E. Radu, Non-Linear Dynamics of Spinning Bosonic Stars: Formation and Stability, Phys. Rev. Lett. 123, 221101 (2019).
  * [7] V. Dzhunushaliev and V. Folomeev, Axially symmetric Proca-Higgs boson stars, Phys. Rev. D 104, 104024 (2021).
  * [8] R. Bartnik and J. Mckinnon, Particle-Like Solutions of the Einstein Yang-Mills Equations, Phys. Rev. Lett. 61, 141 (1988).
  * [9] M. S. Volkov and D. V. Gal’tsov, Gravitating non-Abelian solitons and black holes with Yang-Mills fields, Phys. Rep. 319, 1 (1999).
  * [10] F. Finster, J. Smoller, and S. T. Yau, Particle-like solutions of the Einstein-Dirac equations, Phys. Rev. D 59, 104020 (1999).
  * [11] F. Finster, J. Smoller, and S. T. Yau, Particle-like solutions of the Einstein-Dirac-Maxwell equations, Phys. Lett. A 259, 431 (1999).
  * [12] V. G. Krechet and I. V. Sinilshchikova, Self-gravitating nonlinear spinor field in stationary spaces with spherical symmetry, Russ. Phys. J. 57, 870 (2014).
  * [13] V. Adanhounme, A. Adomou, F. P. Codo, and M. N. Hounkonnou, Nonlinear spinor field equations in gravitational theory: Spherical symmetric soliton-like solutions, J. Mod. Phys. 3, 935 (2012).
  * [14] K. A. Bronnikov, Y. P. Rybakov, and B. Saha, Spinor fields in spherical symmetry: Einstein-Dirac and other space-times, Eur. Phys. J. Plus 135, 124 (2020).
  * [15] K. A. Bronnikov, E. N. Chudaeva, and G. N. Shikin, Self-gravitating string-like configurations of nonlinear spinor fields, Gen. Relativ. Gravit. 36, 1537 (2004).
  * [16] K. A. Bronnikov and J. P. S. Lemos, Cylindrical wormholes, Phys. Rev. D 79, 104019 (2009).
  * [17] M. O. Ribas, F. P. Devecchi, and G. M. Kremer, Fermionic cosmologies with Yukawa type interactions, Europhys. Lett. 93, 19002 (2011).
  * [18] M. O. Ribas, F. P. Devecchi, and G. M. Kremer, Cosmological model with fermion and tachyon fields interacting via Yukawa-type potential, Mod. Phys. Lett. A 31, 1650039 (2016).
  * [19] B. Saha, Nonlinear spinor field in isotropic space-time and dark energy models, Eur. Phys. J. Plus 131, 242 (2016).
  * [20] C. Herdeiro, I. Perapechka, E. Radu, and Y. Shnir, Asymptotically flat spinning scalar, Dirac and Proca stars, Phys. Lett. B 797, 134845 (2019).
  * [21] V. Dzhunushaliev and V. Folomeev, Dirac stars supported by nonlinear spinor fields, Phys. Rev. D 99, 084030 (2019).
  * [22] V. Dzhunushaliev and V. Folomeev, Dirac star in the presence of Maxwell and Proca fields, Phys. Rev. D 99, 104066 (2019).
  * [23] V. Dzhunushaliev and V. Folomeev, Dirac star with SU(2) Yang-Mills and Proca fields, Phys. Rev. D 101, 024023 (2020).
  * [24] V. Dzhunushaliev and V. Folomeev, Static and collapsing configurations supported by the spinor fluid, Phys. Rev. D 103, 123017 (2021).
  * [25] B. Kleihaus, J. Kunz, and M. List, Rotating boson stars and Q-balls, Phys. Rev. D 72, 064002 (2005).
  * [26] T. Wiseman, Static axisymmetric vacuum solutions and nonuniform black strings, Classical Quantum Gravity 20, 1137 (2003).
  * [27] J. L. Friedman, L. Parker, and J. R. Ipser, Rapidly rotating neutron star models, Astrophys. J. 304, 115 (1986).
  * [28] B. Kleihaus, J. Kunz, S. Mojica, and M. Zagermann, Rapidly rotating neutron stars in dilatonic Einstein-Gauss-Bonnet theory, Phys. Rev. D 93, 064077 (2016).
  * [29] W. Schönauer and R. Weiß, Efficient vectorizable PDE solvers, J. Comput. Appl. Math. 27, 279 (1989).
  * [30] F. Cipolletta, C. Cherubini, S. Filippi, J. A. Rueda, and R. Ruffini, Fast rotating neutron stars with realistic nuclear matter equation of state, Phys. Rev. D 92, 023007 (2015).
  * [31] J. L. Friedman, J. R. Ipser, and R. D. Sorkin, Turning point method for axisymmetric stability of rotating relativistic stars, Astrophys. J. 325, 722 (1988).
  * [32] K. Takami, L. Rezzolla, and S. Yoshida, A quasi-radial stability criterion for rotating relativistic stars, Mon. Not. R. Astron. Soc. 416, L1 (2011).