# Rotating Boson Stars Using Finite Differences
and Global Newton Methods

Santiago Ontañón<EMAIL_ADDRESS>Miguel Alcubierre
<EMAIL_ADDRESS>Instituto de Ciencias Nucleares, Universidad
Nacional Autónoma de México, A.P. 70-543, México D.F. 04510, México.

(August 30, 2024)

###### Abstract

We study Rotating Boson Star initial data for Numerical Relativity as
previously considered by Yoshida and Eriguchi [1], Lai [2], and Grandclement,
Somé and Gourgoulhon [3]. We use a 3 + 1 decomposition as presented by
Gourgoulhon [4] and Alcubierre [5], adapted to an axisymmetric quasi-isotropic
spacetime with added regularization at the axis following work by Ruíz,
Alcubierre and Núñez [6] and Torres [7]. The Einstein-Klein-Gordon equations
result in a system of six coupled, elliptic, nonlinear equations with an added
unknown for the scalar field’s frequency $\omega$. Utilizing a Cartesian two-
dimensional grid, finite differences, Global Newton Methods adapted from
Deuflhard [8], the sparse direct linear solver PARDISO [9, 10], and properly
constraining all variables generates data sets for rotation azimuthal integers
$l\in[0,6]$. Our numerical implementation, published in GitHub [11], is shown
to correctly converge both with respect to the resolution size and boundary
extension (fourth order and third order, respectively). Thus, global
parameters such as the Komar masses and angular momenta can be precisely
calculated to characterize these spacetimes. Furthermore, analyzing the full
family at fixed rotation integer produces maximum masses and minimum
frequencies. These coincide with previous results in literature for
$l\in[0,2]$ as in [1, 2, 3, 12] and are new for $l>2$. In particular, the
study of high-amplitude and localized scalar fields in axial symmetry is
revealed to be only possible by adding the sixth regularization variable.

###### pacs:

04.20.Ex, 04.25.Dm, 95.30.Sf

††preprint: IOP/CQG

##  I Introduction

Boson stars are localized configurations of scalar fields first studied by
Kaup [13], and Ruffini and Bonazzola [14]. They arise as a natural coupling of
a smooth classical complex scalar field with Einstein’s theory of general
relativity [15]. Whereas Wheeler’s original idea for a “gravitational atom”,
the geon [16], was proven to be unstable, for a scalar field the gravitational
collapse can be avoided by the dispersive nature of the Klein-Gordon equation.
Although Derrick’s theorem [17, 18] shows that no regular, static, and
localized scalar fields are stable in three-dimensional space, boson stars
solve this conundrum by demoting staticity to stationarity, and assuming a
harmonic time behavior for the complex scalar field of the form

$\Phi(\mathbf{r},t)=\phi(\mathbf{r})\,e^{-i\omega t}\,,$ (I.1)

where $\omega$ is called the scalar field’s frequency that arises as an
eigenvalue when impossing the stationarity condition. Thus, boson stars are
stationary and compact bundles of a complex scalar field, and given some
hypothetical scalar field they may serve different astrophysical purposes, for
example as mimics of neutron stars. Boson stars have also been used to model
dark matter halos by Schunck and Torres [19] and Ureña-López and Bernal [20],
and also as black-hole mimickers by Guzmán and Rueda-Becerril [21] and
Barranco and Bernal [22]. On the question of their actual astrophysical
existence, one must first consider the type of scalar field one has in mind.
Indeed, for a free Klein-Gordon field a simple application of the uncertainty
principle states that for a boson star with constituents of mass $m$ and
localized within a Schwarzschild radius, the maximum mass must be of the order
of

$M_{\text{max}}\sim\frac{\hbar c}{2Gm}=\frac{m_{P}^{2}}{2m}\,,$ (I.2)

where $m_{P}\equiv\sqrt{\hbar c/m}$ is the Planck mass. For a scalar such as
the Higgs boson, the resulting star mass is only $2\times 10^{9}$ kg, which is
a tiny fraction of a Solar Mass. On the other hand, for an interacting field,
such as the quartic interaction studied by Colpi et al. [23], it has been
shown that this limit can be increased almost arbitrarily to astronomically
expected scales. However, throughout this paper we will only consider the case
of a free Klein-Gordon field.

The simplest solution for a boson star is a spherical boson star that
considers a spherically symmetric spacetime and scalar field such that
$\phi(\mathbf{r})=\phi(r)$. This problem is then one-dimensional and has been
studied in the literature multiple times, including numerical results by Kaup
[13] and Friedberg et al. [24], using a shooting method that solves for an
appropriate frequency such that the field and metric quantities decay
correctly at spatial infinity. This solution is constrained by setting the
field’s value at the origin, $\phi_{0}\equiv\phi(r=0)$, and solving the
eigenvalue problem in order to find the frequency $\omega$, which results in a
family of solutions parametrized by the value of $\phi_{0}$. These solutions
have been shown to have a maximum mass of (as reported by Lai [2] and
Grandclement et al. [3])

$M_{\text{max}}^{l=0}\approx 0.633\>\frac{m_{P}^{2}}{m}\,,$ (I.3)

in agreement with the uncertainty principle limit mentioned above (equation
(I.2)).

Here we will consider the case of rotating boson stars which are a natural
extension of the problem into an axisymmetric spacetime. These objects
consider an additional harmonic decomposition in the rotation angle $\varphi$,

$\Phi(\mathbf{r},\,t)=\phi(\mathbf{r})\,e^{-i(\omega t-l\varphi)}\,,$ (I.4)

where $l$ is an integer known as the azimuthal rotation number. As expected,
rotating boson stars have non-zero angular momentum and greater complexities
than their spherical siblings (notice that this is quite different from the
$\ell$-boson stars recently studied by Alcubierre et al. [25], which are
particular combinations of several scalar fields that result in spherical
objects with non-zero total angular momentum).

Rotating boson stars were first studied by Silveira and Sousa in Newtonian
theory [26]. The first general relativity numerical results were obtained by
Yoshida and Eriguchi for the $l=1$ case [1], for which they found a maximum
mass of $M_{\text{max}}^{l=1}\approx 1.31\,m_{P}^{2}/m$. Further analysis was
done by Lai in his PhD thesis [2] where he extended results to the $l=2$ case,
and also found a maximum mass for $l=1$, although it differed from the result
obtained by Yoshida and Eriguchi. Mielke and Schunck [19] found solutions all
throughout $l\in[1,10]$ (and also for $l=500$), however their results were
limited to extremely small amplitude near-Newtonian cases. Finally,
Grandclement, Somé and Gourgoulhon [3] found rotating boson stars for
$l\in[1,4]$ which confirmed Yoshida and Eriguchi’s $l=1$ results, and also
estimated the maximum mass for $l=2$. More background and information on boson
stars can be found in Liebling and Palenzuela’s review [12], which covers the
widest breadth of topics (whereas Mielke’s review considers rotating boson
stars exclusively [27]).

For this paper we have considered free-field boson stars with rotation
azimuthal number $l\in[1,6]$. Differing from Lai’s and Grandclement et al.’s
approach, we have chosen to use a non-compactified spacetime where the grid is
a simple uniform grid with finite difference discretization. However, in
reading both aforementioned works, the issue of the field’s regularity at the
axis was posited as a reason for numerical failure at high amplitudes. Thus,
by looking into previous works by Ruíz, Alcubierre, and Núñez [6], the
regularization of the axis has been a major component in our analysis.

This article is the organized as follows. In Section II we present the field
equations and the deduction of the system of five nonlinear elliptic equations
that describe rotating boson stars. Section III is key in presenting our
boundary conditions, since an inadequate constraint will lead to trivial
Minkowski spacetime or excited boson stars states. Section IV introduces a
sixth variable that guarantees regularity at the axis even for high
amplitudes, and will prove critical for extending previous results. Section V
details how masses and angular momenta were calculated given our grid
limitations. Section VI presents our numerical algorithm for obtaining initial
data for a highly nonlinear problem where initial guesses are hard to
construct. Having presented both the theory and numerical algorithms, Section
VII gives and explanation and overview of the numerical code, ROTBSON
specially written for rotating boson star initial data generation. Section
VIII presents ROTBOSON’s expected fourth-order convergence and gives estimates
for relative errors in global quantities. Finally, our results are presented
in Section IX, and Section X gives some concluding remarks on these results.

A final comment on conventions and units. Throughout this work we use
geometrized units where $c=G=\hbar=1$, unless stated otherwise. In these
units, the Planck mass is given by $m_{P}=\sqrt{\hbar c/G}=1$. This implies
that our coordinate lengths will have dimensions directly inverse to the
field’s mass, $(1/m)$. In the results section we recover Planck’s reduced
constant to present them as is usual in the literature, such as in [3].
Lastly, although all calculations are invariant to this, our metric signature
is the usual convention for general relativity $(-,\,+,\,+,\,+)$.

##  II Field Equations

We use the 3+1 general relativity formalism as presented by York [28],
Alcubierre [5], and Gourgoulhon [4], where the spacetime metric is written as

$\differential s^{2}=g_{\mu\nu}\,\differential x^{\mu}\,\differential
x^{\nu}=(-\alpha^{2}+\gamma_{ij}\,\beta^{i}\,\beta^{j})\,\differential
t^{2}+2\gamma_{ij}\,\beta^{j}\,\differential t\,\differential
x^{i}+\gamma_{ij}\,\differential x^{i}\,\differential x^{j}\,.$ (II.1)

Here $\alpha$ is the lapse function, $\beta^{i}$ is the shift vector, and
$\gamma_{ij}$ is the spatial metric induced on hypersurfaces of constant
coordinate time $\Sigma_{t}$. Furthermore, we consider a stationary and
axisymmetric spacetime. In an adapted coordinate system, we can generate these
symmetries with the stationary generator
$(\partialderivative*{t})^{\mu}=\alpha\,n^{\mu}+\beta^{\mu}$ (where $n^{\mu}$
is the normal vector to the spatial hypersurfaces $\Sigma_{t}$) and the
axisymmetric generator $(\partialderivative*{\varphi})^{\mu}$. An extensive
study of these symmetries and formalism can be found in [4]. From this same
reference, we use quasi-isotropic coordinates $(t,\,\rho,\,z,\,\varphi)$, also
called Lewis-Papapetrou coordinates [29, 30]. Notice that we are using
cylindrical coordinates $(\rho,\,z)$ instead of polar $(r,\,\theta)$. However,
this is perfectly equivalent with the usual transformation
$\rho=r\,\sin\theta$ and $z=r\,\cos\theta$. In this adapted coordinate system
the shift vector can be written as

$\beta^{i}=\Omega(\rho,\,z)\,(\partialderivative*{\varphi})^{i}\,,$ (II.2)

whereas the spatial metric takes the simple form

$\differential l^{2}=\gamma_{ij}\,\differential x^{i}\,\differential x^{j}\\\
=A(\rho,\,z)\,(\differential\rho^{2}+\differential
z^{2})+\rho^{2}\,H(\rho,\,z)\,\differential\varphi^{2}\,.$ (II.3)

Thus, the spacetime geometry is reduced to four functions
$(\alpha,\,\Omega,\,A,\,H)$ of the two coordinates $(\rho,\,z)$. This is not
all, since the Einstein equation will couple geometry to matter via the
stress-energy tensor $T_{\mu\nu}$. Our scalar field will be written as the
standard ansatz with harmonic time and angular dependence [2, 27]

$\Phi(t,\,\rho,\,z,\,\varphi)=\phi(\rho,\,z)\,e^{-i(\omega t-l\varphi)}\,.$
(II.4)

Above, $l$ must be an integer for single-valuedness at $\varphi=0,\,2\pi$, and
is known as the azimuthal rotational number (sometimes referred to as $k$ or
$m$ in other works).

Thus, we need five equations for the five unknown variables
$\\{\alpha,\,\Omega,\,A,\,H,\,\phi\\}$. These are given via:

  1. 1.

The Lewis-Papapetrou or quasi-isotropic coordinates imply that $K$, the trace
of the extrinsic curvature $K_{ij}=-\tfrac{1}{2}\,\pounds_{n}\gamma_{ij}$, is
zero, which results in the maximal slicing condition,

$D^{2}\,\alpha-\alpha\,\left(4\pi\,(u+S)+K_{ij}K^{ij}\right)=0\,,$ (II.5)

  2. 2.

The Hamiltonian constraint,

$\mathcal{H}={}^{(3)}R+K^{2}-K_{ij}-16\pi\,u=0\,,$ (II.6)

  3. 3.

The momentum constraint for $i=\varphi$,

$\mathcal{M}^{i}=D_{j}\,K^{ij}-D^{i}\,K-8\pi\,j^{i}=0\,,$ (II.7)

  4. 4.

The stationarity condition $\partial_{t}K_{ij}=0$ for $i=j=\varphi$,

$\partialderivative{K_{ij}}{t}=\pounds_{\beta}K_{ij}-D_{i}D_{j}\,\alpha+\alpha\,\left({}^{(3)}R_{ij}+K\,K_{ij}-2K_{ik}\,K^{k}\,_{j}+4\pi\,\left((S-u)\,\gamma_{ij}-2S_{ij}\right)\right)=0\,,$
(II.8)

  5. 5.

The Klein-Gordon equation

$\left(\nabla^{2}-m^{2}\right)\,\Phi=0\,.$ (II.9)

Above, $D_{i}$ is the covariant derivative associated to the spatial metric
$\gamma_{ij}$, ${}^{(3)}R_{ij}$ its Ricci tensor, and $u,\,j^{i},\,S_{ij}$ are
respectively the energy density, momentum density and stress tensor as seen by
the Eulerian observers, defined as the following projections of $T_{\mu\nu}$
[5]

$\displaystyle u\;=$ $\displaystyle\,n^{\mu}n^{\nu}\,T_{\mu\nu}\,,$ (II.10)
$\displaystyle j^{i}\;=$
$\displaystyle\,-\gamma^{i\mu}\,n^{\nu}\,T_{\mu\nu}\,,$ (II.11) $\displaystyle
S_{ij}=$
$\displaystyle\,\gamma^{\mu}\,_{i}\,\gamma^{\nu}\,_{j}\,T_{\mu\nu}\,,$ (II.12)
$\displaystyle S\;=$ $\displaystyle\,\gamma^{ij}\,S_{ij}\,,$ (II.13)

and $m$ is the scalar field’s mass parameter.

To obtain the final equations as presented below, linear combinations of these
five equations must be taken. The resulting system of equations is

$\displaystyle\begin{split}f_{\alpha}\equiv&\,\left(\partialderivative[2]{\rho}+\partialderivative[2]{z}+\frac{1}{\rho}\,\partialderivative{\rho}\right)\,\alpha+\frac{1}{2H}\,\partial\alpha\,\cdot\,\partial
H-\frac{\rho^{2}H}{2\alpha}\,\partial\Omega\,\cdot\,\partial\Omega\\\
&-4\pi\,A\,\left(\frac{2(\omega+l\Omega)^{2}}{\alpha}-m^{2}\,\alpha\right)\,\phi^{2}=0\,,\end{split}$
(II.14)
$\displaystyle\begin{split}f_{\Omega}\equiv&\,\left(\partialderivative[2]{\rho}+\partialderivative[2]{z}+\frac{3}{\rho}\,\partialderivative{\rho}\right)\,\Omega+\frac{3}{2H}\,\partial\Omega\,\cdot\,\partial
H-\frac{1}{\alpha}\,\partial\Omega\,\cdot\,\partial\alpha-16\pi\,\frac{A}{H}\,l\,(\omega+l\Omega)\,\left(\frac{\phi}{\rho}\right)^{2}=0\,,\end{split}$
(II.15)
$\displaystyle\begin{split}f_{A}\equiv&\,\left(\partialderivative[2]{\rho}+\partialderivative[2]{z}\right)\,A-\frac{1}{A}\,\partial
A\,\cdot\,\partial A+\left[-\frac{1}{\alpha H}\,\partial
H\,\cdot\,\partial\alpha-\frac{\rho^{2}H}{2\alpha^{2}}\,\partial\Omega\,\cdot\,\partial\Omega-\frac{2}{\rho\,\alpha}\,\partial_{\rho}\alpha\right.\\\
&\left.+8\pi\,\left\\{\left(\frac{\phi}{\rho}\right)^{2}\,\left(\rho^{2}\,\frac{(\omega+l\Omega)^{2}}{\alpha^{2}}\,A-l^{2}\,\frac{H}{A}\right)+\partial\phi\,\cdot\,\partial\phi\right\\}\right]\,A=0\,,\end{split}$
(II.16)
$\displaystyle\begin{split}f_{H}\equiv&\,\left(\partialderivative[2]{\rho}+\partialderivative[2]{z}+\frac{2}{\rho}\,\partialderivative{\rho}\right)\,H-\frac{1}{2H}\,\partial
H\,\cdot\,\partial H+\frac{1}{\alpha}\,\partial H\,\cdot\,\partial\alpha\\\
&\quad+\left(\frac{\rho^{2}H}{\alpha^{2}}\,\partial\Omega\,\cdot\,\partial\Omega+\frac{2}{\rho\,\alpha}\,\partial_{\rho}\alpha+8\pi\,A\,\left(\frac{\phi}{\rho}\right)^{2}\left(\rho^{2}m^{2}+\frac{2l^{2}}{H}\right)\right)\,H=0\,,\end{split}$
(II.17)
$\displaystyle\begin{split}f_{\phi}\equiv&\,\left(\partialderivative[2]{\rho}+\partialderivative[2]{z}+\frac{1}{\rho}\,\partialderivative{\rho}-\frac{l^{2}}{\rho^{2}}\right)\,\phi+\frac{1}{\alpha}\,\partial\alpha\,\cdot\,\partial\phi+\frac{1}{2H}\,\partial
H\,\cdot\,\partial\phi\\\
&+\left(A\,\left(\frac{(\omega+l\Omega)^{2}}{\alpha^{2}}-m^{2}\right)-\frac{l^{2}}{H}\,\left(\frac{A-H}{\rho^{2}}\right)\right)\,\phi=0\,.\end{split}$
(II.18)

These equations use the short-hand defined in [3, 4]:

$\partial u\,\cdot\,\partial
v\equiv\partialderivative{u}{\rho}\,\partialderivative{v}{\rho}+\partialderivative{u}{z}\,\partialderivative{v}{z}\,.$
(II.19)

Equations (II.14) - (II.18) are not quite the final equations we will be
using. In fact, $\\{\alpha,\,\Omega,\,A,\,H,\,\phi\\}$ are also not the actual
variables used in our work. Instead we will work with

$\displaystyle\aleph\equiv$ $\displaystyle\,\log\alpha\,,$ (II.20)
$\displaystyle a\equiv$ $\displaystyle\,\frac{1}{2}\,\log A\,,$ (II.21)
$\displaystyle h\equiv$ $\displaystyle\,\frac{1}{2}\,\log H\,,$ (II.22)
$\displaystyle\psi\equiv$ $\displaystyle\,\frac{\phi}{\rho^{l}}\,.$ (II.23)

The first three variables are introduced to guarantee that $\alpha,\,H,\,A>0$
(as should be from metric-positivity). The final variable comes from examining
regularity conditions in Eq. (II.18). More precisely, the first four terms
that make up the elliptic operator in this equation imply that, if everything
else in the equation is regular at the rotation axis $\rho=0$, $\phi$ must
then vanish as $\rho^{l}$ on this same axis. This might not seem true,
considering the term $(A-H)/\rho^{2}$, however, this term is shown to be
regular at the axis via local-flatness, as will be further discussed in
Section IV below. The immediate implication of this decomposition is that the
term $\phi/\rho=\rho^{l-1}\,\psi$ is regular for $l\geq 1$, which is true for
all rotating boson stars. Using this new $\psi$ variable, its resulting
equation is

$\displaystyle\begin{split}f_{\psi}\equiv&\,\left(\partialderivative[2]{\rho}+\partialderivative[2]{z}+\frac{(2l+1)}{\rho}\,\partialderivative{\rho}\right)\,\psi+\frac{1}{\alpha}\,\partial\alpha\,\cdot\,\partial\psi+\frac{1}{2H}\,\partial
H\,\cdot\,\partial\psi\\\
&+\left(\frac{l}{\rho}\,\left(\frac{\partial_{\rho}\alpha}{\alpha}+\frac{\partial_{\rho}H}{2H}\right)+A\,\left(\frac{(\omega+l\Omega)^{2}}{\alpha^{2}}-m^{2}\right)-\frac{l^{2}}{H}\,\left(\frac{A-H}{\rho^{2}}\right)\right)\,\psi=0\,.\end{split}$
(II.24)

Notice also that using $\psi$ is convenient since otherwise the field’s first
$l$ derivatives on the $\rho$ direction would need to vanish, something very
difficult to guarantee numerically. This was reported to be an issue in
Yoshida and Eriguchi’s first work [1].

We will not rewrite equations (II.14) - (II.17) and (II.24) in terms of
$\\{\aleph,\,\Omega,\,a,\,h,\,\psi\\}$ since the change of variables is quite
simple. Thus, we have a set of five nonlinear elliptic equations to solve for
the initial data. Unfortunately, this system has a trivial solution given by
$\alpha=A=H=1$, $\Omega=\phi=0$, which corresponds to Minkowski spacetime. We
obviously do not wish to obtain this spacetime configuration, so a first step
is to constrain solutions to have $|\phi|>0$, or, equivalently $|\psi|>0$. In
fact, we can be more specific, since we are searching for ground-state
solutions where the scalar field has no “nodes” (save at the rotation axis)
and thus we will ask for $\psi>0$. This constraint will be further developed
in the next section.

##  III Boundary Conditions

First a comment on our grid structure: we will be solving this system of
equations with a two-dimensional cartesian $(\rho,\,z)$ uniform grid. We could
use a more complex non-uniform or adaptive grid, but for simplicity’s sake, we
have opted to use the easier approach where the grid spaces are constant and
given by $\Delta\rho$ and $\Delta z$ in their respective directions. As is
common for Numerical Relativity, a staggered grid is used with the goal of
avoiding possible problematic divisions by zero on the $\rho$ axis, and also
to allow us to impose equatorial symmetry in a simple way. Thus, our discrete
coordinates are given by

$(\rho_{i},\,z_{j})=\left(\Delta\rho\,(i-g+\tfrac{1}{2}),\,\Delta
z\,(j-g+\tfrac{1}{2})\right)\,,$ (III.1)

where $g$ is the number of ghost zones necessary for finite differences, and
$i$, $j$ are integers in the intervals:

$i\in[0,\,N_{\rho}+2g-1]\,,\quad j\in[0,\,N_{z}+2g-1]\,.$ (III.2)

This means that $N_{\rho},\,N_{z}$ are the number of interior points in each
respective direction, and we have two ghost zones boundary bands with the
first interior point being $(i,\,j)=(g,\,g)$ and the last interior point
corresponding to $(i,\,j)=(g+N_{\rho}-1,\,g+N_{z}-1)$.

Notice that our system of equations also has equatorial symmetry. Thus, the
most optimal solution algorithm only solves for a single quadrant of the
$(\rho,\,z)$ coordinates, namely, the positive quadrant. The remaining
quadrants can be obtained by parity transformations. It is easy to see that
invariance under the independent symmetries $\rho\to-\rho$ and $z\to-z$
implies that all our functions are even functions of both variables. This is a
boundary condition for our grid, specifically for our left and bottom strips.
The mathematical expression is

$u(i,\,j)=u(2g-1-i,\,j)\,,\,i\in[0,\,g-1]\,,\quad
u(i,\,j)=u(i,\,2g-1-j)\,,\,j\in[0,\,g-1]\,,$ (III.3)

for $u\in\\{\aleph,\,\Omega,\,a,\,h,\,\psi\\}$.

We also need boundary conditions for the external boundaries $i\geq
N_{\rho}+g,\,j\geq N_{z}+g$. Asymptotic flatness has the simple condition
$\alpha,\,A,\,H=1$ and $\Omega=\psi=0$ at spatial infinity. However, we do not
have a compactified spacetime and must use another condition instead.
Enforcing the previous condition as a Dirichlet approximation at finite
distance is not ideal, rather we use the Robin boundary conditions [31, 32]
that come from expanding a grid variable $u$ asymptotically for large $r$ as:

$u(r,\,\theta)=u_{\infty}+\frac{C_{n}}{r^{n}}+\mathcal{O}(r^{-(n+1)},\,\theta)\,,$
(III.4)

where $u_{\infty}$ is the value or limit of $u$ at spatial infinity, $C_{n}$
is a constant independent of the direction $\theta$, $n$ a positive integer,
and further terms can depend on the direction $\theta$ but are of order
$r^{-(n+1)}$ or greater.

The actual Robin boundary condition is written as a mix of Dirichlet and
Neumann conditions which is a natural consequence of the above Equation
(III.4):

$r_{\text{bdy}}\,\partialderivative{u}{r}\,(r_{\text{bdy}},\,\theta)+n\,(u(r_{\text{bdy}},\,\theta)-u_{\infty})=\mathcal{O}(r_{\text{bdy}}^{-(n+1)},\,\theta)\,.$
(III.5)

Above, $r_{\text{bdy}}$ stands for an external boundary point such that
$r_{\text{bdy}}\gg 1$. The numerical implementation sets the right-hand side
of (III.5) to zero and, for our two dimensional grid, we have
$r\,\partialderivative*{u}{r}=\rho\,\partialderivative*{u}{\rho}+z\,\partialderivative*{u}{z}$.

Notice that the unknown constant $C_{n}$ has dropped out from our boundary
condition. This is not so for the integer $n$, which must be known a priori.
This integer can be calculated from expanding the system of equations around
spatial infinity and determining the first nonzero power of $1/r$ for each
variable. This is tedious since we first have to rewrite the equations in
$(r,\,\theta)$ coordinates, but quite straightforward. In this way we find
that $n=1$ for $\aleph,\,a,\,h$, and $n=3$ for $\Omega$. In fact, it can also
be shown that the constant $C_{n}$ must have the same value for $a$ and $h$,
and minus this value for $\aleph$. This will be no surprise when we talk about
global quantities in Section V. Also, the Dirichlet condition at spatial
infinity means that all variables $\\{\aleph,\,\Omega,\,a,\,h\\}$ have
$u_{\infty}=0$.

We have yet to deal with the scalar field. Its decay is in fact not a power of
$1/r$, but instead exponential. This can be seen by expanding equation (II.18)
at spatial infinity and keeping only the lowest order in $1/r$, i.e.
$\alpha=A=H=1$ and $\Omega=0$. This leads to

$\left(\partialderivative[2]{r}+\frac{2}{r}\,\partialderivative{r}+\frac{1}{r^{2}}\,\partialderivative[2]{\theta}+\frac{1}{r^{2}\,\tan\theta}\,\partialderivative{\theta}-\frac{l^{2}}{r^{2}\sin^{2}\theta}-(m^{2}-\omega^{2})\right)\,\phi\approx
0\,.$ (III.6)

We recognize the Laplacian operator in cylindrical coordinates above. By
separation of variables, one can show that the appropriate solution of (III.6)
is

$\phi=C\,P^{l}_{l}(\cos\theta)\,k_{l}(\sqrt{m^{2}-\omega^{2}}\,r)\,,$ (III.7)

where $P^{l}_{l}(\cos\theta)$ is the associated Legendre polynomial,
$k_{l}(\sqrt{m^{2}-\omega^{2}}\,r)$ is the modified spherical Bessel function,
and $C$ is a constant. This solution is the appropriate one since it goes to
zero at spatial infinity, and vanishes as $\sin^{l}\theta$ on the rotation
axis. Indeed, from [33], given the asymptotic form of the Bessel function
$k_{l}(x)\sim e^{-x}/x$, and a direct calculation of
$P^{l}_{l}(\cos\theta)\sim\sin^{l}\theta$, one finds,

$\phi\to C\,\sin^{l}\theta\,\frac{e^{-\sqrt{m^{2}-\omega^{2}}\,r}}{r}\,,$
(III.8)

or, in terms of our variable $\psi=\rho^{-l}\,\phi$,

$\psi\to C\,\frac{e^{-\sqrt{m^{2}-\omega^{2}}\,r}}{r^{l+1}}\,.$ (III.9)

Observe that this asymptotic behavior requires that (reinserting Planck’s
Constant)

$0<\omega<\frac{m}{\hbar}\,,$ (III.10)

as otherwise there would not be an appropriate exponential decay of the scalar
field. Thus, the differential equation for the boundary condition for $\psi$
is

$r_{\text{bdy}}\,\partialderivative{\psi}{r}\,(r_{\text{bdy}},\,\theta)+\left((l+1)+r_{\text{bdy}}\,\sqrt{m^{2}-\omega^{2}}\right)\,\psi(r_{\text{bdy}},\,\theta)=\mathcal{O}(r_{\text{bdy}}^{-(l+2)},\,\theta)\,.$
(III.11)

In conclusion, the above boundary conditions have given us a full list of the
constraints our variables must satisfy:

  1. 1.

The metric requires $\alpha,\,A,\,H$ positive, and decaying as $1/r$
asymptotically.

  2. 2.

Gourgoulhon [4] shows that in terms of our notation we must have $\Omega<0$.
Also, $\Omega$ decays as $1/r^{3}$.

  3. 3.

The scalar field auxiliary variable $\psi$ must always be positive to
guarantee that we have the ground-state, and must decay as
$e^{-\sqrt{m^{2}-\omega^{2}}\,r}/r^{l+1}$. Nevertheless, to avoid the trivial
Minkowski solution, $\psi$ must still be nonzero at the boundary, however
small.

  4. 4.

Finally, the field’s frequency is bounded by the scalar field mass parameter.

Condition 1 is satisfied by taking logarithmic variables as previously
mentioned. Conditions 2 and 3 could be guaranteed by taking another variable
such as $\log(-\Omega)$ and $\log(\psi\,e^{\sqrt{m^{2}-\omega^{2}}\,r})$.
However, this leads to Robin boundary conditions where the value at spatial
infinity is infinite. Therefore, we do not constrain these variables in
practice, but violations of the aforementioned constraints would show that we
are solving the wrong problem. Lastly, Condition 4 is affirmed by working with
a new sigmoid variable that guarantees it

$\xi\equiv\tanh^{-1}\,\left(2\,\frac{\omega}{m}-1\right)\,\quad\Rightarrow\quad\omega=\frac{m}{2}\,\left(1+\tanh(\xi)\right)\,.$
(III.12)

##  IV Regularization

Before proceeding with technical and algorithmic details, we should carefully
examine that all terms in our equations are not only analytically regular at
both axes, but also, that we can write them in a fashion that is less likely
to cause numerical problems. Terms such as
$(1/\rho)\,\partialderivative*{u}{\rho}$ (for a given
$u\in\\{\aleph,\,\Omega,\,a,\,h,\psi\\}$) are regular since all grid functions
are even about the rotation axis. Therefore, $\partialderivative*{u}{\rho}$
behaves as $\mathcal{O}(\rho)$ at this same axis, and its combined division by
$\rho$ is an even function. Nevertheless, it is recommended to write these
terms with the order of operations explicit in numerical code, i.e., instead
of writing in source code

$\frac{l}{\rho}\,\left(\partialderivative{\aleph}{\rho}+\partialderivative{h}{\rho}\right)\,\psi\,,$

write

$l\,\left(\left(\frac{1}{\rho}\,\partialderivative{\aleph}{\rho}\right)+\left(\frac{1}{\rho}\,\partialderivative{h}{\rho}\right)\right)\,\psi\,.$

As previously stated, terms such as $\phi/\rho$ are regular since the scalar
field vanishes at least as $\rho$ on the axis. With this in mind, it is easy
to verify that all terms in our equations are regular, with the possible
exception of a single term in Equation (II.24) proportional to:

$\lambda\equiv\frac{A-H}{\rho^{2}}\,.$ (IV.1)

An initial approximation might seem to indicate that this should not be an
issue, given that we are using a staggered grid. Thus, no actual division by
zero is ever done. Indeed, the smallest value of the $\rho$ coordinate is,
from (III.1), $\rho=\rho_{g}=\Delta\rho/2$. Thus, we should not see
regularization problems for low resolutions, i.e., for $\Delta\rho$ relatively
big (empirically, $\Delta\rho\gtrsim 0.1$). However, as almost all problems in
this field, a greater resolution is best to reduce truncation errors, and it
is also necessary to closely examine small regions of spacetime, specially as
the scalar field amplitude is increased.

Furthermore, as will be seen in Section IX, as we consider scalar fields with
greater amplitude, the field concentrates closer to the rotation axis and
develops a sharp maximum, i.e., the $\rho$ coordinate location of the field’s
maximum becomes smaller and smaller. Thus, we need greater resolution to
properly resolve and characterize this spike. As we increase the resolution,
$\Delta\rho$ decreases and we start seeing regularization problems in equation
(II.24) coming from the term (IV.1).

It is worthwhile to point out that these regularization issues come from our
curvilinear coordinates. In a full 3D Cartesian scheme, $1/\rho$ or $1/r$
terms are absent. However, curvilinear coordinates have an obvious advantage
when dealing with symmetric spacetimes such as axisymmetry where the Killing
vectors are fully adapted. This means that a 3D calculation is reduced to 2D
(for our axysimmetric work), making it far less computationally expensive.

For an exact solution the term (IV.1) is in fact perfectly regular. Indeed, a
simple argument of local flatness shows that near the rotation axis, we can
expand $A$ and $H$ as

$A(\rho,\,z)=\gamma_{\rho\rho}=C(0,\,z)+\frac{\rho^{2}}{2}\,A_{2}(0,\,z)+\mathcal{O}(\rho^{4})\,,\quad
H(\rho,\,z)=\frac{\gamma_{\varphi\varphi}}{\rho^{2}}=C(0,\,z)+\frac{\rho^{2}}{2}\,H_{2}(0,\,z)+\mathcal{O}(\rho^{4})\,,$
(IV.2)

i.e., from local flatness the value at the rotation axis, $C(0,\,z)$ must be
equal for both metric functions, and thus $\lambda$ is a regular and even
function at the rotation axis,

$\lambda(\rho,\,z)=\frac{(A_{2}(0,\,z)-H_{2}(0,\,z))}{2}+\frac{\rho^{2}}{2}\,\lambda_{2}(0,\,z)+\mathcal{O}(\rho^{4})\,.$
(IV.3)

This clearly implies that the term in (II.24), $-(l^{2}\,\lambda/H)\,\psi$, is
analytically regular at the axis. But, unfortunately, numerical methods
introduce truncation and round-off errors that may cause that this exact
cancellation will not hold numerically, and an increase in resolution will
introduce loss of significance if we divide $(A-H)$ by a small quantity. To
deal with this potential catastrophe, we follow previous work by Milton,
Alcubierre and Núñez [6], and Torres [7], whence the quantity $\lambda$ is
promoted to a new independent variable. This is expensive in terms of
computational resources since we are introducing, intuitively, at least $20\%$
more work. Also, $\lambda$ must have its own elliptic equation and boundary
conditions, otherwise the problem is indeterminate. Obtaining this equation is
tedious work and nontrivial. A first guess is to take equations $f_{A}$ and
$f_{H}$, (II.16) and (II.17), subtract them and divide them by $\rho^{2}$.
Afterwards, we construct an elliptic operator for $\lambda$, i.e.,
$\partialderivative*[2]{\lambda}{\rho}+\partialderivative*[2]{\lambda}{z}$,
and regularize all terms. For example, terms that initially appear irregular
can be combined into a regular expressions,

$\frac{1}{\rho^{2}}\,\left(\partialderivative[2]{\aleph}{\rho}-\frac{1}{\rho}\,\partialderivative{\aleph}{\rho}\right)=\frac{1}{\rho}\,\partialderivative{\rho}\,\left(\frac{1}{\rho}\,\partialderivative{\aleph}{\rho}\right)\,.$
(IV.4)

Above we have the derivative of an even function, since we previously noted
that $(1/\rho)\,\partialderivative*{\aleph}{\rho}$ was regular and even.
Therefore, its derivative is odd, and its division by $\rho$ turns out to be
an overall regular expression.

Unfortunately, the combination $(f_{A}-f_{H})/\rho^{2}$ turns out to be
inappropriate to reduce via this procedure in order to generate a regular
equation. Happily, however, we have a large amount of equations to work with.
Previously we specified the five equations that give rise to our system of
equations. But there are in fact many more equations. Some are trivial, such
as the momentum constraints $M^{i}$ for $i=\rho,\,z$, but stationarity
$\partial_{t}K_{ij}=0$ yields more equations for $i,j\neq\varphi$. The proper
combination to obtain a regular equation for $\lambda$ turns out to be:

$f_{\lambda}\equiv-\frac{2A}{\rho^{2}\alpha}\,\left(\partial_{t}K_{\rho\rho}-\frac{1}{\rho^{2}}\,\partial_{t}K_{\varphi\varphi}\right)=\partialderivative[2]{\lambda}{\rho}+\partialderivative[2]{\lambda}{z}+\frac{3}{\rho}\,\partialderivative{\lambda}{\rho}+\ldots=0\,.$
(IV.5)

This is a rather large equation so we will not write it explicitly here (the
full equation is given in Appendix A, Equation (A.1)). For our discussion, we
only want to notice the elliptic operator in Equation (IV.5).

Before continuing, we must specify boundary conditions for $\lambda$. We have
already determined that $\lambda$ is an even function of $\rho$. By examining
the full equation (A.1), we can also deduce that it is even about the $z$
axis, which is fortunate since otherwise we would need to solve another
quadrant of $(\rho,\,z)$ space. For its external boundary condition, now that
we have the asymptotic behavior of the other variables, a similar expansion at
spatial infinity shows that $\lambda$ must decay as $1/r^{4}$, i.e.

$r_{\text{bdy}}\,\partialderivative{\lambda}{r}\,(r_{\text{bdy}},\,\theta)+4\lambda(r_{\text{bdy}},\,\theta)=\mathcal{O}(r_{\text{bdy}}^{-5},\,\theta)\,.$
(IV.6)

##  V Global Quantities

It is of great interest to characterize our solutions with a set of parameters
that globally describe the spacetime. This can be done via conserved
quantities which for General Relativity are associated with spacetime
symmetries. Possibly the most general concept is the Arnowitt-Deser-Misner
(ADM) mass which only requires that the spacetime be asymptotically flat. From
Wald [34], this mass is defined as

$M_{\text{ADM}}=\lim_{\mathscr{S}\to\infty}\,\frac{1}{16\pi}\,\oint_{\mathscr{S}}\,\left(\mathcal{D}^{j}\,\gamma_{ij}-\mathcal{D}_{i}\,(f^{kl}\gamma_{kl})\right)\,s^{i}\,\sqrt{q}\,\differential\theta\,\differential\varphi\,,$
(V.1)

where a fiducial metric $f_{ij}$ is introduced on $\Sigma_{t}$ which must be
Euclidian, i.e. flat. Since the integral must be calculated over 2-spheres
$\mathscr{S}$, the obvious choice is
$f_{ij}=\text{diag}(1,\,r^{2},\,r^{2}\,\sin^{2}\theta)$. In the above
expression furthermore $\mathcal{D}_{i}$ is $f_{ij}$’s Levi-Civita connection,
$s^{i}$ the normal vector induced on the 2-spheres, and $\sqrt{q}$ the volume
element with respect to the physical metric $\gamma_{ij}$. Thus, one must take
special care to distiguish the quantities that depend on the physical metric,
$s^{i}$ and $\sqrt{q}$, from those that depend on the flat fiducial metric,
$f^{ij}$ and $\mathcal{D}_{i}$. With this in mind, the resulting expression
for our particular metric (II.3) is [4]

$M_{\text{ADM}}=\lim_{r_{\text{bdy}}\to\infty}\,-\frac{1}{8}\,\int\limits_{0}^{\pi}\,\differential\theta\,\left(\partialderivative{r}\,(A+H)+\frac{H-A}{r}\right)\,r^{2}\,\sin\theta\,.$
(V.2)

It is well known that the ADM mass converges very slowly with radius. Even for
a Schwarzschild black-hole with mass parameter $M$ (in isotropic coordinates
so that $A=H=(1+M/2r)^{4}$) the expression turns out to be

$M_{\text{ADM}}=\lim_{r\to\infty}M\,\left(1+\frac{M}{2r}\right)^{3}=M+\mathcal{O}(r^{-1})\,.$
(V.3)

In fact, to have a relative error of less than $1\%$ in this expression, we
must have the external boundary at $r_{\infty}>150M$, which can be very
computationally expensive. Fortunately, alternatives have been derived. In
Schwarzschild’s case, a pseudo-mass has been introduced by Alcubierre et al.
[35]

$M_{\text{PS}}=\left(\frac{\tilde{A}}{16\pi}\right)^{1/2}\,\left(1-\frac{\left(\derivative*{\tilde{A}}{r}\right)^{2}}{16\pi\,\tilde{\gamma}_{rr}\,\tilde{A}}\right)\,,$
(V.4)

where $\tilde{A}$ is the area of a 2-sphere, $r$ is the radial coordinate, and
$\tilde{\gamma}_{rr}$ is the average of $\gamma_{rr}$ over the 2-sphere. This
expression is ideal for the case of an areal radial coordinate, where
$\tilde{A}=4\pi\,r^{2}$, so that

$M_{\text{PS}}=\frac{r}{2}\,\left(1-\frac{1}{\gamma_{rr}}\right)\,.$ (V.5)

For Schwarzschild’s case in the areal coordinates we have
$\gamma_{rr}=(1-2M/r)^{-1}$, and thus we can easily see that the pseudo-
Schwarzschild mass gives an exact expression, $M_{\text{PS}}=M$, independent
of our radial coordinate $r$. This certainly is an improvement from the
$\mathcal{O}(1/r)$ convergence in the ADM case. However, notice that it will
not always be possible to work with an areal radial coordinate, especially for
spacetimes that are not spherically-symmetric. This is certainly not possible
in our case, since $A\neq H$ for rotating boson stars. But even so, we will
later show that in the limit in which $A\approx H$ this expression yields a
better result than the ADM mass at finite radius.

Better suited to our problem is the Komar mass [34]. Here the only added
assumptions are that the spacetime has a Killing vector $\xi^{\mu}$ associated
to stationary symmetry, and that the Ricci tensor
$R_{\mu\nu}=8\pi\,(T_{\mu\nu}-g_{\mu\nu}T/2)$ is zero beyond a 2-sphere
$\mathscr{S}$. The expression for the Komar mass is

$M_{\text{Komar}}=\frac{1}{8\pi}\,\oint_{\mathscr{S}}\,\nabla^{\mu}\xi^{\nu}\,\differential
S_{\mu\nu}\,,$ (V.6)

where $\differential S_{\mu\nu}$ is the 2-form normal to the surface
$\mathscr{S}$, i.e. induced by $n^{\mu}$ and $s^{\mu}$, the normal vector to
the hypersurface $\Sigma_{t}$ and the normal vector to $\mathscr{S}$, so that
$\differential
S_{\mu\nu}=(s_{\mu}n_{\nu}-s_{\nu}n_{\mu})\,\sqrt{q}\,\differential\theta\,\differential\varphi$.

Our problem does not actually satisfy the fact that $T_{\mu\nu}=0$ beyond any
2-sphere $\mathscr{S}$. However, this necessity can be relaxed as long as this
tensor goes sufficiently rapid to zero at spatial infinity [34]. Indeed, the
stress energy tensor is quadratic in the scalar field, and we determined
previously that the scalar must decay exponentially, so that at sufficiently
large $r$, the field is dominated by this decay,

$T_{\mu\nu}\to e^{-2\sqrt{m^{2}-\omega}\,r}\,.$ (V.7)

Therefore, an expression can be calculated in terms of our metric variables as
long as we take the limit at spatial infinity, $r_{\text{bdy}}\to\infty$. We
find:

$M_{\text{Komar}}=\lim_{r_{\text{bdy}}\to\infty}\,\frac{1}{2}\int\limits_{0}^{\pi}\,\differential\theta\,\left(\partialderivative{\alpha}{r}-\frac{Hr^{2}\sin^{2}\theta}{2\alpha}\,\Omega\,\partialderivative{\Omega}{r}\right)\,\sqrt{H}\,r^{2}\sin\theta\,.$
(V.8)

In addition, another expression can be obtained by using Stokes’ theorem and
Einstein’s equations to convert the surface integral into a volume integral.

$M_{\text{Komar}}=2\,\int_{\Sigma}\,\left(T_{\mu\nu}-\frac{T}{2}\,g_{\mu\nu}\right)\,n^{\mu}\xi^{\nu}\,\differential
V\,.$ (V.9)

Above, $\Sigma$ is a spatial hypersurface such that its boundary is the
2-sphere $\mathscr{S}$. Like in the previous case, even though $R_{\mu\nu}\neq
0$, if we take the limit at spatial infinity, the expression can hold. In this
case, we just integrate over the complete hypersurface $\Sigma_{t}$. There is,
however, another caveat to this expression: it is clear that for a black-hole
like Schwarzschild’s, $R_{\mu\nu}=T_{\mu\nu}=0$, yet the mass is certainly not
zero. The issue here is the presence of a spacetime singularity which does not
allow an application of Stokes’ theorem. Thus, (V.9) must not be used in the
presence of singularities. However, all our functions are regular and smooth
at all our hypersurfaces by construction, so we can use the Komar mass in
terms of a volume integral [4], which in our case reduces to:

$M_{\text{Komar}}=\lim_{r_{\text{bdy}}\to\infty}\,2\pi\,\int\limits_{0}^{r_{\text{bdy}}}\,\differential
r\int\limits_{0}^{\pi}\,\differential\theta\,\left(\frac{2\omega}{\alpha}\,(\omega+l\Omega)-\alpha\,m^{2}\right)\,\phi^{2}\,A\,\sqrt{H}\,r^{2}\sin\theta\,.$
(V.10)

This last expression makes explicit the fact that the integral must converge
exponentially, since all terms are quadratic in the scalar field. This
certainly is a very different behavior than that of the ADM mass. Also, having
two expressions for the same quantity can be very useful for numerical
convergence-monitoring.

Nevertheless, a priori there is no connection between the ADM mass and the
Komar mass, but it has been shown that they are identical if the Killing
vector $\xi^{\mu}$ is orthogonal to $\Sigma_{t}$ at spatial infinity. This is
true in our case since, as we approach spatial infinity we have

$n^{\mu}=\frac{1}{\alpha}\,\left((\partialderivative*{t})^{\mu}-\beta^{\mu}\right)\to(\partialderivative*{t})^{\mu}=\xi^{\mu}\,.$
(V.11)

In conclusion, if we obtain a rapidly convergent Komar mass, this is
equivalent to calculating the ADM mass. The clear advantage is that the Komar
mass may be calculated within a smaller computational domain, which is ideal
for our numerical implementation. In other words, as long as the scalar field
is appropriately confined within our computational domain, we will be able to
calculate spacetime’s global mass parameter.

Our spacetime possesses another symmetry, namely axisymmetry under the Killing
vector $\chi^{\mu}=(\partialderivative*{\varphi})^{\mu}$. An ADM and Komar
analogue quantity may be derived: the angular momentum. Once again, the ADM
formalism converges too slowly and a Komar expression is preferred [34, 4].
For our metric the final expression becomes:

$\displaystyle\begin{split}J_{\text{Komar}}=&\,\lim_{\mathscr{S}\to\infty}\,\frac{1}{16\pi}\,\oint\,\frac{H\,r^{2}\,\sin^{2}\theta}{\alpha\,\sqrt{A}}\,\partialderivative{\Omega}{r}\,\differential
A=\lim_{r_{\text{bdy}}\to\infty}\,\frac{1}{8}\,\int\limits_{0}^{\pi}\,\differential\theta\,\frac{H^{3/2}\,r^{4}\,\sin^{3}\theta}{\alpha}\,\partialderivative{\Omega}{r}\,,\end{split}$
(V.12)
$\displaystyle\begin{split}=&\,l\,\int_{\Sigma_{t}}\,\frac{(\omega+l\,\Omega)}{\alpha}\,\phi^{2}\,\differential
V=\lim_{r_{\text{bdy}}\to\infty}\,2\pi\,l\,\int\limits_{0}^{r_{\text{bdy}}}\,\differential
r\,\int\limits_{0}^{\pi}\,\differential\theta\,\frac{(\omega+l\Omega)}{\alpha}\,\phi^{2}\,A\,\sqrt{H}\,r^{2}\sin\theta\,,\end{split}$
(V.13)

where we have written both expressions in terms of surface and volume
integrals. Analogously to the previous expression (V.10), the volume integral
makes explicit the fact that the integral will converge exponentially as
$r_{\text{bdy}}\to\infty$.

##  VI Global Newton Methods

Initially, our system of equations for
$\\{f_{\aleph},\,f_{\Omega},\,f_{a},\,f_{h},\,f_{\psi},\,f_{\lambda}\\}$
appears to consist of six elliptic, nonlinear, coupled equations. However,
this can be misleading. When we discretize each variable
$\\{\aleph,\,\Omega,\,a,\,h,\,\psi,\,\lambda\\}$ on our grid
$(\rho_{i},\,z_{j})$ we obtain six equations for each grid point. For example,
for an illustrative $800$ interior points in each direction, we have $6\times
800\times 800\approx 4\times 10^{6}$ equations, and each function
$\\{\aleph(\rho_{i},\,z_{j}),\,\Omega(\rho_{i},\,z_{j}),\,a(\rho_{i},\,z_{j}),\,h(\rho_{i},\,z_{j}),\,\psi(\rho_{i},\,z_{j}),\,\lambda(\rho_{i},\,z_{j})\\}$
is a true variable we must solve for. Adding to the complexity is that we also
must calculate the field’s frequency, $\omega$, as a sort of eigenvalue. It is
easy to see that we do not have enough equations by adding $\omega$ (actually
$\xi$ as defined in Equation (III.12) above) as an unknown. The solution is to
remove another variable, i.e., constrain the value of $\psi$ at a certain grid
point, which also solves the previous requirement that $\psi>0$ in order to
avoid the trivial Minkowski solution. Of course, we can fix $\psi$ at any
point in theory, however, it is best to constrain it near the axis and origin,
i.e., the first interior point $\psi(\Delta\rho/2,\,\Delta z/2)$. Otherwise,
we might be constraining it at a point where it is already too small given its
exponential decay. Due to the axial and equatorial symmetries, this is
equivalent to constraining the field’s value at the origin:

$\psi_{0}\equiv\psi(0,\,0)\,.$ (VI.1)

Now, suppose that for each $\psi_{0}$ (and each rotation number $l$) there
exists a unique $\omega$ for the ground-state (the converse is not true, i.e.,
for a single $\omega$ there can be multiple $\psi_{0}$’s, as will be seen
later in Section IX), then the system can be solved with $\omega$ as an
unknown and the above constraint $\psi(0,\,0)=\psi_{0}$.

This is can be done by a Newton-Raphson method [36, 37]. The following
discussion on Newton methods and their globalization closely follows Deuflhard
[8]; other globalization strategies have also been developed, for example, by
Bank and Rose [38]. Specifically, we will use a finite difference
discretization [39] and calculate the Jacobian matrix with respect to the
discrete variables and $\omega$ (since it is also an unknown variable).
Newton’s method is a very powerful tool used to solve nonlinear equations,
specifically when we have a system of equations from a domain
$D\subset\mathbb{R}^{n}$ into a codomain space $U\subset\mathbb{R}^{n}$, i.e.,
for $x\in\mathbb{R}^{n}$, we seek the solution of the $n$ equations
represented as a vector in $\mathbb{R}^{n}$

$F(x)=0\,.$ (VI.2)

The standard Newton method relies on an initial guess $x^{0}$ and a basic
linearization of the equation:

$F(x^{k}+\Delta x^{k})=F(x^{k})+F^{\prime}(x^{k})\,\Delta
x^{k}+\mathcal{O}((\Delta x^{k})^{2})\,,$ (VI.3)

where $F^{\prime}(x^{k})$ is the $n\times n$ Jacobian matrix at $x^{k}$. Thus,
Newton affirms that a “better” estimate is constructed iteratively as

$x^{k+1}\equiv x^{k}+\Delta x^{k}\,,\quad\Delta
x^{k}=-F^{\prime}(x^{k})^{-1}\,F(x^{k})\,.$ (VI.4)

The (famous) advantage of this method is that it converges quadratically [40],
however it can be very sensitive to the initial guess $x^{0}$. The “more
nonlinear” the problem, this sensitivity is exacerbated. Unfortunately, for a
general problem, we might really have no idea as to what is a “good initial
guess”. Therefore, we would like to globalize Newton’s method, so that in the
best case scenario any type of initial guess leads to the true solution. This
is immediately limited by the fact that a system of equations might have
multiple (or none at all) solutions. However, if we constrain our solution
appropriately via boundary conditions, we may recover uniqueness. Thus,
although a formal proof can be very cumbersome, uniqueness is critical and we
will normally “build-in” constraints into the very same variables. For
example, in our problem we have guaranteed $m^{2}>\omega^{2}$ by introducing
the variable $\xi$ (III.12), and we are also constraining the metric variables
to be positive by taking their logarithms. Also, as previously discussed, in
looking for the scalar field’s ground state, we should enforce $\psi>0$.

Our first step in this discussion is having the means to determine if a
problem requires a globalization, or if the classical Newton algorithm will
suffice. In other words, we must determine if a problem is “highly nonlinear”.
First, we notice that the classical existence, uniqueness and convergence
theorems [40, 41] for Newton’s method rely on having a priori that
$F^{\prime}(x)^{-1}$ exists and is bounded across the domain $D$,

$\norm{F^{\prime}(x)^{-1}}\leq\beta<\infty\,,\quad x\in D\,.$ (VI.5)

However, this is computationally problematic. First of all, notice that this
norm is taken over matrices. Also, $\beta$ is realistically near impossible to
calculate for nontrivial problems. The best we can do is to sample across the
domain

$\norm{F^{\prime}(x^{0})^{-1}}\leq\beta_{0}\,,\quad x^{0}\in D\,.$ (VI.6)

But even this is still expensive as we have to explicitly calculate the
inverse $F^{\prime}(x^{0})^{-1}$ which can be numerically very expensive. The
solution is to use another condition instead of (VI.5), namely a Lipschitz
condition

$\norm{F^{\prime}(x)-F^{\prime}(y)}\leq\gamma\,\norm{x-y}\,,\quad x,\,y\in
D\,.$ (VI.7)

This leads to the Newton-Kantorovich theorem [40] that asserts existence and
uniqueness of a solution $x^{*}$, and quadratic convergence for initial guess
$x^{0}$ that is in a neighborhood characterized by the Kantorovich quantity

$h_{0}\equiv\norm{x^{*}-x^{0}}\,\beta_{0}\,\gamma<\frac{1}{2}\,,$ (VI.8)

where $\beta_{0}$ and $\gamma$ were defined in (VI.6) and (VI.7),
respectively. However, this does not seem to solve the problem, since we must
still somehow calculate the Lipschitz quantity, $\gamma$, and this again
implies norms over matrices. Fortunately, a reformulation of the Newton-
Kantorovich theorem [8] exists which only requires vector norms and another
Lipschitz constant, $\epsilon$:

$\norm{F^{\prime}(x)^{-1}\,\left(F^{\prime}(x)-F^{\prime}(y)\right)\,\left(x-y\right)}\leq\epsilon\,\norm{x-y}^{2}\,,\quad
x,\,y\in D\,.$ (VI.9)

Furthermore, with this condition, the Kantorovich quantity is written simply
as

$h_{0}\equiv\norm{x^{*}-x^{0}}\,\epsilon<2\,.$ (VI.10)

In other words, we have unique quadratic convergence for a neighborhood of
radius [8]

$\rho_{0}\equiv\norm{x^{*}-x^{0}}=\frac{2}{\epsilon}\,.$ (VI.11)

As previously stated, we now have only vector norms in both sides of (VI.9),
which is computationally advantageous. Still, calculating the Lipschitz
constant $\epsilon$ seems to be a hopeless task based on random sampling.
However, condition (VI.9) has an enormous structural advantage in that it is
affine covariant, i.e., it is invariant under the transformation of the
original system of $n$ equations, $F(x)=0$ into

$G(x)=A\,F(x)=0\,,$ (VI.12)

where $A\in\text{GL}(n)$ is a nonsingular $n\times n$ matrix. This unique
property allows a clever construction of a scalar function $g(x,\,y,\,z)$ such
that

$\epsilon=\text{sup}\,g(x,\,y,\,z)\,,\quad x,\,y,\,z\in D\,.$ (VI.13)

Once this function, $g$, is obtained, we can calculate a computational local
estimate

$\left[\epsilon_{0}\right]=g(x^{1},\,x^{0},\,x^{0})\leq\epsilon\,.$ (VI.14)

A straightforward example (although not the actual estimate used in
algorithms) is [8]

$\left[\epsilon_{0}\right]=\frac{\norm{F^{\prime}(x^{0})^{-1}\,\left(F^{\prime}(x^{1})-F^{\prime}(x^{0})\right)}}{\norm{x^{1}-x^{0}}}\,.$
(VI.15)

There will obviously be an error $\epsilon-\left[\epsilon_{0}\right]$, however
an efficient adaptive Newton algorithm can be constructed in such a way that
successive iterations catch more binary digits of the true $\epsilon$, via the
bit counting lemmas as presented by Deuflhard [42, 8]. In summary, Newton’s
method, the vector norm Lipschitz condition (VI.9), and affine covariance
gives us a way to calculate an approximate $\left[\epsilon\right]$. Via the
convergence radius (VI.11), we can now estimate whether we have an initial
guess “sufficiently close” to the solution, $x^{*}$, i.e., we characterize a
problem as highly nonlinear [8] if the Newton update is such that

$\norm{\Delta x^{0}}=x^{1}-x^{0}\gg\frac{2}{\left[\epsilon_{0}\right]}\,.$
(VI.16)

In such a case, Newton’s method will not have guaranteed convergence and
requires a globalization. This can only be constructed using additional
structure on the equations $F(x)$. A globalization of Newton’s method should
deal with (within reason) “bad” initial guesses and should merge with the
classical Newton method when sufficiently near to the solution $x^{*}$ as
determined via Equation (VI.11) since we wish to make use of Newton method’s
quadratic convergence.

There exists multiple globalization strategies, however the most intuitive is
probably that of steepest descent [43] where the iterates $x^{k}$ successively
approach the solution point whose direction is calculated via the gradient of
the residual level function

$T(x)\equiv\frac{1}{2}\,\norm{F(x)}^{2}=\frac{1}{2}\,F(x)^{T}\,F(x)\,.$
(VI.17)

This is a concave function, since $T(x)=0$ if and only if $x=x^{*}$ and
$T(x)>0$ if and only if $x\neq x^{*}$. Thus, a global solution can easily be
found via a monotonicity criterion

$T(x^{k+1})<T(x^{k})\quad\iff\quad\norm{x^{k+1}-x^{*}}<\norm{x^{k}-x^{*}}\,.$
(VI.18)

A formal iterative method is thus

$\displaystyle\Delta x^{k}\equiv$
$\displaystyle\,\gradient{T(x^{k})}=-F^{\prime}(x^{k})^{T}\,F(x^{k})\,,$
(VI.19) $\displaystyle x^{k+1}\equiv$ $\displaystyle\,x^{k}+s_{k}\,\Delta
x^{k}\,,$ (VI.20) $\displaystyle s^{k}>$ $\displaystyle\,0\,,$ (VI.21)

where $s^{k}$ is a steplength parameter, which is chosen to satisfy the
monotonicity criterion (VI.18) by taking the successive steps
$i=0,\,1,\,\ldots$ until it is true that

$T(x^{k}+s^{i}_{k}\,\Delta x^{k})<T(x^{k})\,,$ (VI.22)

via a steplength strategy

$s^{i+1}_{k}\equiv\kappa\,s^{i}_{k}\,,\quad\kappa<1\,,\quad s^{0}_{k}=1\,,$
(VI.23)

where $\kappa$ is usually chosen to be $\tfrac{1}{2}$. This method can be
shown to have linear convergence [8] for arbitrarily “bad” initial guesses
$x^{0}$ and it is also clear that as the steplength goes to one, we recover
Newton’s classical method. However, this convergence, although linear, can be
also arbitrarily slow. Worse still is the fact that we can find “pseudo-
convergence” characterized by small

$\norm{F^{\prime}(x)^{T}\,F(x)}\ll 1\,.$ (VI.24)

This discussion on the method of steepest descent is not fruitless since it
serves to introduce the analogous and more powerful extensions where we do not
utilize the level function (VI.17) but a general level function to fully
utilize affine covariance [8]

$T(x\,|A)\equiv\frac{1}{2}\,\norm{A\,F(x)}^{2}\,,\quad A\in\text{GL}(n)\,.$
(VI.25)

Using this general level function, it is easily shown that its minimization
still points “downhill” independently of the matrix $A$ (notice that $A=I$
gives the steepest descent method).

Several choices of matrices are possible to implement different methods,
including the trust-region models of Levenberg-Marquardt [44, 45]. The choice
of $A=F^{\prime}(x^{k})^{-1}$ gives the error oriented descent where the
iterates now are required to satisfy the natural monotonicity test [42, 8]

$\norm{\overline{\Delta x}^{k+1}}<\norm{\Delta x^{k}}\,,$ (VI.26)

where we introduce the simplified Newton correction

$\overline{\Delta x}^{k+1}\equiv F^{\prime}(x^{k})^{-1}\,F(x^{k+1})\,.$
(VI.27)

Furthermore, we now interpret $s_{k}$ as a damping factor

$x^{k+1}=x^{k}+s_{k}\,\Delta x^{k+1}\,,$ (VI.28)

that must be chosen via a damping strategy that obeys the criteria (VI.26) and
goes to one once we are within the local convergence radius (VI.11). This
damping strategy is given by Theorem 3.12 of Deuflhard [8] and states that
given a Lipschitz condition (VI.9) and a Kantorovich quantity

$h_{k}\equiv\norm{\Delta x^{k}}\,\epsilon\,,$ (VI.29)

the natural monotonicity criteria (VI.26) is satisfied and the optimal choice
is

$s_{k}\equiv\min(1,\,1/h_{k})\,.$ (VI.30)

Deulfhard also gives a proof for global convergence in Theorem 3.14 of his
book [8]. Our only remaining limitation is that this theoretical damping
strategy cannot be implemented directly and only via the computational
available estimates $[\epsilon_{k}]$ introduced in Equation (VI.14) and
$[h_{k}]=[\epsilon_{k}]\,\norm{\Delta x^{k}}$. The supremum property implies
that

$[s_{k}]\equiv\min(1,\,1/[h_{k}])\geq s_{k}\,,$ (VI.31)

so that the estimated damping factor may be “too large”. This conundrum is
addressed with a prediction and correction strategies. In essence, these
strategies use trial iterates and a user-input of an initial damping factor
$s_{0}$ to guarantee the monotonicity criteria. These are also termination
criteria in the algorithm since the damping factor may become arbitrarily
small for singular problems. More details on these criteria and their
optimization is presented in Deuflhard [8] where the complete implementation
of these concepts results in the error oriented global Newton method, called
NLEQ-ERR [46, 8], that serves as our solution algorithm presented in Appendix
B.

##  VII Code Overview

All previous functionalities are implemented in a numerical code specially
tailored to this problem, ROTBOSON, openly available at a GitHub repository
[11]. ROTBOSON is a finite difference discretization code spanning over 20,000
lines of C++ OpenMP [47] code influenced by other cousins OllinAxis [48] and
OllinBrill [49] developed at the ICN, UNAM. Its design philosophy centers on
being able to run on personal or professional-grade machines and as such it is
scalable and optimized. Its main usage is detailed in the GitHub repository
but is briefly summarized in the interest of the following discussion.

First, the user must compile the code with either the GNU or Intel C++
compiler. An installation of the freely available Intel Math Kernel Library
software [50] is also necessary since it used for the sparse direct solver and
other vector operations. Once compiled, the user must provide a parameter file
with a list of relevant parameters such as the grid resolutions
$\Delta\rho,\,\Delta z$; the number of interior points $N_{\rho},\,N_{z}$; the
finite differentiation order (2 or 4); the scalar field’s rotation number $l$,
an initial guess for the scalar frequency $\omega_{0}$; and other I/O options
to read or interpolate previously obtained datasets. A selection of parameters
inside a typical parameter file is presented below.

⬇

1# GRID

2dr = 0.015625 # Step size in rho (r) direction.

3dz = 0.015625 # Step size in z direction.

4NrInterior = 512 # Number of interior points in r direction.

5NzInterior = 512 # Number of interior points in z direction.

6order = 4 # Finite difference discretization order 2|4.

7

8# SCALAR FIELD PROPERTIES

9l = 4 # Rotation number.

10m = 1.0 # Field mass parameter.

11

12# INITIAL DATA I/O.

13readInitialData = 2 # Select no initial data 0, direct read 1, interpolation
2.

14# File names for variables.

15log_alpha_i = "./log_alpha_i.asc"

16beta_i = "./beta_i.asc"

17log_h_i = "./log_h_i.asc"

18log_a_i = "./log_a_i.asc"

19psi_i = "./psi_i.asc"

20lambda_i = "./lambda_i.asc"

21w_i = "./w_i.asc"

22

23# FIXED VARIABLE.

24fixedPhi = 1 # Select if scalar field is fixed and if so, where.

25fixedPhiR = 2

26fixedPhiZ = 2

27

28# SOLVER PARAMETERS

29solverType = 1 # Select Newton solver: 0 classical, 1 NLEQ-ERR, 2 NLEQ-RES.

30localSolver = 1 # Enable local solver.

31epsilon = 1.0E-10 # Convergence exit tolerance.

32maxNewtonIter = 20 # Maximum number of iterations.

33lambda0 = 1.0E-03 # Initial damping factor.

34lambdaMin = 1.0E-06 # Minimum damping factor.

35useLowRank = 1 # Use Low-Rank optimization 1, otherwise use CGS
preconditioning.

36

37# Parameter file continues...

A critical element of this code is its sparse direct solver (the theory behind
these solvers is beyond the scope of this work but may be consulted in
Tewarson [51], and Gupta and Kumar [52], among others). Indeed, about 50% of
its lines consist of generating Jacobian matrices at different discretization
orders (this was done automatically via Python and Mathematica [53] scripts).
Once the matrices are written in sparse CSR format they are used in Intel
MKL’s included software, PARDISO [54, 9, 10]. Since Newton’s method solves
iteratively, the solver can be optimized greatly given the fact that succesive
Jacobians have the same sparse structure, thus an $LU$ sparse factorization
scheme only requires a single analysis phase which can then be saved for other
executions with entirely different parameters (such as rotation number and
step size) but identical sparse structure. For more details on PARDISO’s
optimization for problems in Numerical Relativity, see Ontanon [49].

In terms of memory use and scalability, ROTBOSON has been tested succesfully
at a cluster node running on up to 54 processing units and with a grid size of
approximately $1000\times 1000$ points, ROTBOSON uses about 120 Gb of memory.
For a personal machine, adequate results can be obtained over 4 processing
units and a grid of $400\times 400$ which uses instead 12 Gb of RAM memory.
For now, execution is limited to a single node since the sparse solver loses
significant speed when forced to run under a message passing interface, such
as MPI [55], over multiple nodes.

##  VIII Error Indicators

Before reporting any results, we must verify that our numerical implementation
converges appropriately. Since we are using finite differences, we should
expect that they introduce truncation errors. For all the following results,
finite differences were implemented to fourth order (although second order is
also supported for faster performance and testing). Furthermore, the exit
tolerance for Newton’s algorithm is always set to a relative error
$\epsilon<10^{-10}$, i.e., in terms of the language of Section VI, by taking
2-norms, one of these conditions must be true to exit and guarantee a solution
$x^{*}$

$\frac{\norm{\Delta
x^{k}}}{\norm{x^{k}}}<10^{-10}\,,\quad\text{or}\quad\,\frac{\norm{\overline{\Delta
x^{k}}}}{\norm{x^{k}}}<10^{-10}\,,$

To check for fourth order convergence we vary the resolutions
$\Delta\rho=\Delta z$, but we must fix the external boundary $r_{\text{bdy}}$
in all solutions. For example, we consider eight different configurations with
interior points

$N_{\rho}=N_{z}\in\\{128,\,256,\,384,\,512,\,640,\,768,\,896,\,1024\\}\,,$

and $\Delta\rho=\Delta z$ chosen such that in all cases we have
$r_{\text{bdy}}=128\,(\hbar/m)$ according to Equation (III.1). Furthermore,
since we are fixing the position of the boundary, this is equivalent to
constricting $\omega$, and in our tests we take $\omega=0.7\,(m/\hbar)$. Now,
given these eight solutions, we must compare quantities to test for fourth
order convergence. A possibility is to interpolate the coarse resolutions into
the finest and make a point-wise comparison. However, to avoid compounding
interpolation errors and to simplify, we can make convergence tests with
global quantities, namely $M_{\text{Komar}}$ and $J_{\text{Komar}}$ in their
surface and volume integral expressions (Equations (V.8), (V.10), (V.12),
(V.13)). Another error indicator is the scalar field’s value at the origin
$\psi_{0}$, obtained by interpolation at each solution.

Using these indicators and varying resolutions, we compare, i.e. subtract in
absolute value, two quantities at neighboring resolutions. For example, for
two resolutions $\Delta\rho_{0}>\Delta\rho_{1}$ and their reported masses
$M(\Delta\rho_{0}),\,M(\Delta\rho_{1})$ we define the relative error as

$\epsilon(\Delta\rho_{1})\equiv\epsilon(M,\Delta\rho_{0},\Delta\rho_{1})=\absolutevalue{M(\Delta\rho_{1})-M(\Delta\rho_{0})}\,,$
(VIII.1)

Our convergence results are shown as $\log_{10}$-$\log_{10}$ plots in Figures
VIII.1 and VIII.2. In Figure VIII.1 an “S” or “V” denotes whether the
expression was calculated as a surface or a volume integral. Both Figures
include a reference straight line for $\mathcal{O}((\Delta\rho)^{4})$ fourth
order convergence. From the Figures it is clear that we have the expected
fourth order convergence. In the case of $\psi_{0}$, the error stagnates at
about $\varepsilon\sim 10^{-11}$, which serves to indicate that point-wise it
has reached a limit due to machine round-off error instead.

Figure VIII.1: Komar mass and angular momentum differences between resolutions
as error indicators for $l=6$, $\omega=0.7\,(m/\hbar)$, at fixed
$r_{\text{bdy}}=128\,(\hbar/m)$ and varying resolution $\Delta\rho$.
Figure VIII.2: Scalar field value at origin $\psi_{0}$ difference between
resolutions as error indicator for $l=6$, $\omega=0.7\,(m/\hbar)$, at fixed
$r_{\text{bdy}}=128\,(\hbar/m)$ and varying resolution $\Delta\rho$.

Now that we have determined fourth-order convergence, we can find an error
approximation via Richardson extrapolation [56] or our global quantities.
Thus, at the finest resolution with $\Delta\rho=0.125\,(\hbar/m)$ and
$N_{\rho}=1024$, we find for this example

$M_{\text{Komar}}=(5.550074\pm 0.000002)\,(m_{P}^{2}/m)\,\quad\text{and}\quad
J_{\text{Komar}}=(37.1475\pm 0.0001)\,(m_{P}^{2}/m)^{2}\,,$

which correspond to relative errors of $4\times 10^{-5}\,\%$ and $3\times
10^{-4}\,\%$, respectively.

Further analysis and testing is prudent, specifically with respect to varying
the position of the boundary. In this case, we expect that error arises from
the Robin-type boundary conditions such as Equation (III.5) where the right-
hand side was not zero, but rather $\mathcal{O}(r_{\text{bdy}}^{-(n+1)})$ for
an $n$-type decay. Since the lapse and metrics have $n=1$, an initial estimate
might lead us to expect that there is second order convergence with respect to
$r_{\text{bdy}}$. To verify this, other convergence tests are done now with
fixed resolution $\Delta\rho=\Delta z=0.125\,(\hbar/m)$ and varying
$r_{\text{bdy}}$ (which is equivalent to varying the number of interior
points). Once again, global quantities are used as error monitors and are
presented in Figure VIII.3 as $\log_{10}$-$\log_{10}$ plots.

Figure VIII.3: Komar mass and angular momentum as error indicators for $l=6$,
$\omega=0.7\,(m/\hbar)$, for fixed resolution $\Delta\rho=0.125\,(\hbar/m)$
and varying boundary position $r_{\text{bdy}}$.

In this case, we can conclude that the errors due to the boundary still
converge. By adding the $\mathcal{O}(r_{\text{bdy}}^{-3})$, the results
indicate that we have third order convergence instead of the initially
expected second order. However, given the complexities of this analysis, it
would be best to not read too much into these results and instead just
conclude that the error is indeed decreasing as we place the boundary further
away, even though the magnitude of the boundary error is at least an order of
magnitude greater than the resolution error (in this example, the boundary
error is about 30 times greater than the resolution error). It is clear that
further analysis and studies are necessary for boundary problems and
approximations in Numerical Relativity, but this is outside our current scope.

##  IX Results for $l\in[0,\,6]$

Having verified that the numerical implementation converges correctly, we now
present our results for $l\in[0,\,6]$, which are summarized in Figure IX.1.
All solutions therein have a maximum $\omega=0.9\,(m/\hbar)$. Most solutions
are done with grids of $800\times 800$ interior points and $\Delta\rho=\Delta
z=0.04\,(\hbar/m)$. However, for $l=1,2$ a finer resolution of
$\Delta\rho=\Delta z=0.02\,(\hbar/m)$ is necessary to examine the high-
amplitude behavior. Conversely, for $l=5,6$ a coarser resolution of
$\Delta\rho=\Delta z=0.08\,(\hbar/m)$ is required at low amplitudes
$\omega>0.8\,(m/\hbar)$ for which the field’s width increases, and in order to
properly fit it in our computational domain (given our computational
resources) the boundary is placed farther away and the resolution is
decreased. The changes in resolution have been done seamlessly, i.e. where the
difference between the global parameters is less than the truncation error.

Figure IX.1: Rotating boson stars’ characterization for $l\in[0,\,6]$.
Included global parameters are the Komar mass, the Komar angular momentum, the
field’s frequency $\omega$, the field’s maximum $\max(\phi)$, and the radial
coordinate of this maximum.

Notice that the spherically symmetric case is also included for reference as
$l=0$. In this case the solution is found using an analogous method, but it is
worth pointing out that they are not simply obtained by setting $l=0$ in our
elliptic equations. For more information on our results for spherical boson
stars see [57]. Spherical boson stars [12] trivially have zero angular
momentum, $J=0$, and the field’s maximum is always located at the origin [3].

Some important maximums and minimums are presented in Table IX.1. Given our
error analysis, and considering the worst-case scenario where the error’s main
contributors are the boundary conditions, we can still report five significant
figures for masses and frequencies, and four significant figures for angular
momenta. These critical points were obtained by interpolating a fourth-degree
spline into the data points and calculating the derivative’s zeros. These
results are in agreement for $l=0,1,2$ with those given in references [3, 12,
2, 19, 27]. From Table IX.1 we also can report new “turning-points” (i.e. the
minimum value of $\omega$) and maximum masses which, to our knowledge, are not
yet present in the literature.

Boson Star $l$ | Maximum $M$ ($m_{\text{P}}^{2}/m$) | Maximum $J$ $(m_{\text{P}}^{2}/m)^{2}$ | Minimum $\omega$ $(m/\hbar)$
---|---|---|---
0 | 0.63300 | 0 | 0.76754
1 | 1.3155 | 1.382 | 0.64561
2 | 2.2159 | 4.810 | 0.51657
3 | 3.5287 | 12.49 | 0.44339
4 | 5.0590 | 25.83 | 0.40756
5 | 6.6681 | 44.63 | 0.38819
6 | 8.2824 | 69.02 | 0.37391
Table IX.1: Critical points for boson stars $l\in[0,\,6]$.

Out of the hundreds of computed data sets, we turn now to examine a couple of
illustrative examples for the $l=6$ case at three different frequencies
$\omega\approx\\{0.8\,(m/\hbar),\,0.6\,(m/\hbar),\,0.37\,(m/\hbar)\\}$. These
correspond, respectively, to “low”, “medium”, and “high” scalar field
amplitudes. The latter case is at the turning point for $l=6$ reported in
Table IX.1. These results are presented as contour plots for the physical
variables $\\{\alpha,\,\Omega,\,H,\,A,\,\phi,\,\lambda\\}$ in Figures IX.2,
IX.3, IX.4. Although, as previously stated, equatorial and axial symmetry
imply that our computational domain can be restricted to the quadrant
$\rho,\,z>0$, the figures are plotted in the full space by reflecting about
the corresponding axes. A three-dimensional picture can be obtained or
imagined by rotating the $z$ axis. Notice that the color bar scale is not
equal between figures.

The contours show that, as we decrease $\omega$, the field’s spatial extension
shrinks and at the same time that its amplitude grows. This behavior is also
manifest in the physical variables, namely the metric coefficients
$\Omega,\,H,\,A,\,\lambda$, with an opposite behavior for the lapse (i.e. the
lapse decreases as the other metric functions grow). Notice that in the “high-
amplitude” case the central value of the lapse is below $3\times 10^{-2}$.

Considering our three bosons stars with different amplitudes, the results
might seem to indicate variable trends as a function of decreasing $\omega$.
However, caution is necessary since $\omega$ is not monotonic as a function of
$\psi_{0}$. Rather, as previously discussed the solutions are obtained by
varying $\psi_{0}$. Therefore, an improved picture is obtained in Figure IX.5
that shows how the variables behave as we vary this parameter (which has
dimensions $(\hbar/m)^{l}$ given that $\phi$ is dimensionless). This behavior
is also studied for the $l=1$ case in Figure IX.6. These two figures show the
true trends of the physical variables. As $\psi_{0}$ increases, the lapse’s
minimum decreases, whereas $(-\Omega),\,A,\,H,\,(-\lambda)$ increase. The
field’s frequency behavior is not monotonic, and its critical points lead to
the turning points in the $M$ vs. $\omega$ diagrams.

The field’s localization behavior is also examined for fixed
$\omega=0.8\,(m/\hbar)$ and varying rotation number $l$ in Figure IX.7.
Clearly, only at $l=0$ the maximum may be present at the origin. This figure
also shows that an increase in rotation number brings a decrease in the
field’s amplitude and another increase in terms of the field’s extension.

As previously discussed, we prefer the Komar mass instead of the ADM mass
because of the latter’s slow convergence. Indeed, Figure IX.8 shows how the
ADM mass from Equation (V.1) converges too slowly at
$r_{\text{bdy}}=32\,(\hbar/m)$. This is heavily contrasted with the Komar
equivalents which converge exponentially. Notice that the difference between
the ADM and Komar masses is about 50% at $r_{\text{bdy}}$, making the former
useless for our results. An analogous results is presented in Figure IX.9 for
angular momenta where the Komar expressions Equations (V.12) and (V.13) are
compared to the asymptotic form

$J_{\text{Asymptotic}}=\lim_{r_{\text{bdy}}\to\infty}\,\frac{1}{8}\,\int\limits_{0}^{\pi}\,\differential\theta\,r^{4}\,\sin^{3}\theta\,\partialderivative{\Omega}{r}\,.$
(IX.1)

The final result of this section is Figure IX.10, which demonstrates the
importance of regularization for a “high-amplitude” case such as $l=6$,
$\omega=0.37\,(m/\hbar)$. The yellow curve represents the regularization
variable $\lambda$, whereas the blue curve is the expression $(A-H)/\rho^{2}$
calculated using the two other variables $A,\,H$. Notice that the vertical
scale is logarithmic so that at the origin we have $\lambda\approx 2\times
10^{2}\,(m/\hbar)^{2}$, whereas $(A-H)/\rho^{2}$, though indeterminate at the
origin, at the first interior point is about $8\times 10^{4}\,(m/\hbar)^{2}$
which is an enormous difference of two orders of magnitude. Beyond $r\gtrsim
1\,(\hbar/m)$ both expressions match and can be used interchangeably. However,
$r\ll 1\,(\hbar/m)$ is precisely the region where the scalar field is
localized and of non-negligible amplitude. Thus, due to the
$-l^{2}\,\lambda\,\psi/H$ term in Equation (II.24), it is clear that a
solution scheme with regularization will have substantial advantages over one
without regularization. In our experience, we found that without
regularization our solver cannot converge at high amplitudes precisely for
this reason. Therefore, thanks to regularization, it has been possible to
further explore rotating boson stars solutions at higher amplitudes.

Figure IX.2: Rotating boson star with $l=6$ for “low-amplitude” $\omega\approx
0.8\,(m/\hbar)$.

Figure IX.3: Rotating boson star with $l=6$ for “medium-amplitude”
$\omega\approx 0.6\,(m/\hbar)$.

Figure IX.4: Rotating boson star with $l=6$ for “high-amplitude”
$\omega\approx 0.4\,(m/\hbar)$.

Figure IX.5: Rotating boson stars with $l=6$. We show the behaviour of the
physical variables as functions of $\psi_{0}$ for the range
$\psi_{0}\in[1.0\times 10^{-8},\,3\times 10^{-3}]\,(\hbar/m)^{6}$, at
resolutions $\Delta\rho=\Delta z=0.04\,(\hbar/m)$ and
$r_{\text{bdy}}=32\,(\hbar/m)$.

Figure IX.6: Rotating boson stars with $l=1$. We show the bahavious of the
physical variables as functions of $\psi_{0}$ for the range
$\psi_{0}\in[1.0\times 10^{-2},\,2.0]\,(\hbar/m)$, at resolutions
$\Delta\rho=\Delta z=\\{0.02,\,0.04\\}\,(\hbar/m)$ and
$r_{\text{bdy}}=32\,(\hbar/m)$.

Figure IX.7: Boson stars’ scalar field profiles at the equator for fixed
$\omega=0.8\,(m/\hbar)$ and $l\in[0,6]$.

Figure IX.8: Masses comparison and convergence along radial coordinate $r$ for
Komar surface, Komar volume, ADM and pseudo-Schwarzschild expressions.
Figure IX.9: Angular momenta comparison and convergence along radial
coordinate $r$ for Komar surface, Komar volume and the asymptotic expression.

Figure IX.10: Comparison between the regularization variable $\lambda$ and its
expression in terms of metric variables $(A-H)/\rho^{2}$. Both expressions are
plotted in absolute value using a logarithmic scale. Notice the considerable
differences near the origin, and how they match beyond $r\gg 1\,(\hbar/m)$.

##  X Summary and Conclusions

In this paper we have studied solutions for rotating boon stars, i.e. a
stationary and axisymmetric ground-state of the Klein-Gordon complex scalar
field coupled to Einstein’s equations. The ADM equations initially result in a
system of five coupled and highly-nonlinear equations with an added eigenvalue
problem for the field’s time frequency, $\omega$. Guided by boundary
conditions at both spatial infinity and the $\rho,\,z$ axes, we propose
variables that mostly automatically satisfy the system’s physical constraints.
Given the curvilinear formulation of our variables, terms with powers of
$1/\rho$ appear throughout the equations. While most these terms are regular
under parity conditions, a single term in the scalar field’s equations
proportional to $\lambda=(A-H)/r^{2}$ is not immediately regular. Using local
flatness this term is shown analytically regular, but $\lambda$ must be
promoted numerically to a new independent variable to guarantee regularity at
the rotation axis: thus a sixth elliptic equation is introduced in conjunction
with its proper boundary conditions.

For solving our system of equations we have developed an open-source numerical
code, ROTBOSON, based on finite difference fourth-order discretization and a
Cartesian grid in the cylindrical coordinates $(\rho,\,z)$. This set up is
deal for curvilinear axisymmetric spacetimes where the Killing vectors
$\xi^{a}=(\partialderivative*{t})^{a}$ and
$\chi^{a}=(\partialderivative*{\varphi})^{a}$ are adapted to the coordinates.
Due to the high non-linearity and lack of simple initial guesses to input into
a regular Newton-Raphson method, an alternative Global Newton Method using
error oriented affine covariant Lipschitz conditions is implemented to solve
for an initial seed for each family of rotating boson stars. From this seed,
the entire branch is computed by varying the scalar field’s value at the
origin $\psi_{0}$.

We present data sets for rotation numbers $l\in[0,6]$. These families are
characterized and studied using global parameters such as Komar masses and
angular momenta. The Komar quantities are chosen instead of the ADM
alternatives because of their exponential convergence, which makes them ideal
for a non-compactified grid where spatial infinity is an approximation at
$r_{\text{bdy}}$. Analyzing the results at fixed rotation integer we determine
maximum masses and minimum frequencies. In particular,
$M_{\text{max}}^{l=0}=0.633\,m_{P}^{2}/m$;
$M_{\text{max}}^{l=1}=1.316\,m_{P}^{2}/m$;
$M_{\text{max}}^{l=2}=2.216\,m_{P}^{2}/m$. These latter masses coincide with
previous results in literature found by Yoshida and Eriguchi [1], and
Grandclement, Somé and Gourgoulhon [3]. Most notably, we also obtain maximum
masses and turning-point frequencies for $l>4$ which are new results in the
literature to our knowledge. Given the fact that as the scalar field’s
amplitude increases, it also approaches the rotation axis, high resolutions
are needed to properly resolve and study the field. The new maximum masses for
$l>2$ directly correspond to these limit cases, and their proper study is
revealed to be only possible by adding the sixth regularization variable, thus
proving again the necessity of adequate regularization in curvilinear
formulations at high resolution.

Future work will concentrate in studying these datasets (including the new
high-amplitude regularized members) regarding their stability properties both
under 3D and 2D axisymmetric perturbations. Also, given the high versatility
of Global Newton Methods, an extension to neighboring problems such as
interacting scalar fields and scalar hair is also possible.

###### Acknowledgements.

This work was partially supported by CONACyT Network Projects No. 376127
“Sombras, lentes y ondas gravitatorias generadas por objetos compactos
astrofísicos”, and No. 304001 “Estudio de campos escalares con aplicaciones en
cosmología y astrofísica”, and by a CONACyT National Graduate Grant.

## Appendixes

##  Appendix A Complete Regularization Equation

Here we present the full equation for the regularization variable $\lambda$
that was obtained following the discussion in Section IV:

$\displaystyle\begin{split}f_{\lambda}\equiv&\,\partial_{\rho}^{2}\lambda+\partial_{z}^{2}\lambda+\frac{3}{\rho}\,\partial_{\rho}\lambda\\\
&\,+\frac{1}{2H}\,\left(\partial_{z}H\,\partial_{z}\lambda-\partial_{\rho}H\,\partial_{\rho}\lambda\right)-\frac{2}{A}\,\partial\lambda\cdot\partial
H-\frac{\rho^{2}}{A}\,\partial\lambda\cdot\partial\lambda+\frac{1}{\alpha}\,\left(\partial_{z}\alpha\,\partial\lambda-\partial_{\rho}\alpha\,\partial_{\rho}\lambda\right)\\\
&\,-\frac{4\lambda^{2}}{A}-\frac{4\lambda\,\rho}{A}\,\partial_{\rho}\lambda-\frac{2\lambda}{\alpha}\,\left(\frac{\partial_{\rho}\alpha}{\rho}\right)+\left(\frac{\partial_{\rho}H}{\rho}\right)\,\left(-\frac{4}{A}+\frac{1}{H}\right)\,\lambda-\left(\frac{\partial_{\rho}H}{\rho}\right)^{2}\,\left(\frac{1}{A}+\frac{1}{2H}\right)\\\
&\,+\frac{(\partial_{z}H)^{2}}{A\,H}\,\lambda+\frac{(\partial_{\rho}H)^{2}}{2H^{2}}\,\lambda-\frac{2}{\alpha}\,\left(\frac{\partial_{\rho}H}{\rho}\right)\,\left(\frac{\partial_{\rho}\alpha}{\rho}\right)-\frac{H^{2}}{\alpha^{2}}\,\partial\Omega\cdot\partial\Omega-\frac{A\,H}{\alpha^{2}}\,(\partial_{\rho}\Omega)^{2}\\\
&\,+\frac{2\lambda}{\alpha}\,\partial_{\rho}^{2}\alpha+\frac{\lambda}{H}\,\partial_{\rho}^{2}H+\left[\frac{1}{\rho}\,\partial_{\rho}\,\left(\frac{\partial_{\rho}H}{\rho}\right)\right]+\frac{2H}{\alpha}\,\left[\frac{1}{\rho}\,\partial_{\rho}\,\left(\frac{\partial_{\rho}\alpha}{\rho}\right)\right]\\\
&\,+8\pi\,(\rho^{l-1})^{2}\,A\,\left(\rho^{2}\,m^{2}\,\lambda+2\,\left(\frac{\partial_{\rho}\psi}{\rho}\right)\,\left(2l\,\psi+\rho\,\partial_{\rho}\psi\right)\right)=0\,.\end{split}$
(A.1)

Notice how all terms are written in a regularized manner. For example,
parentheses explicitly contain the overall-even terms
$(1/\rho)\,\partial_{\rho}u$, and the combination

$\frac{1}{\rho^{2}}\,\left(\partial_{\rho}^{2}-\frac{\partial_{\rho}u}{\rho}\right)=\left[\frac{1}{\rho}\,\partial_{\rho}\,\left(\frac{\partial_{\rho}u}{\rho}\right)\right]\,,$
(A.2)

is also regular by differentiating the even function
$((1/\rho)\,\partial_{\rho}u)$ and then dividing by $\rho$: thus resulting in
an even and regular expansion.

##  Appendix B NLEQ-ERR Algorithm

The Global Newton Method, NLEQ-ERR, is taken from Defulhard [8].

Inputs: Error accuracy $\epsilon$ sufficiently above machine precision,
initial iterate $x^{0}$, right-hand side $F(x^{0})$, and damping factor
$s_{0}$.

Output: Newton solution $x^{*}$ for a convergent algorithm.

for  _$k=0,\,1,\,\ldots$_ do

1\. Step $k$: Evaluate Jacobian matrix $F^{\prime}(x^{k})$. Solve linear
system

$F^{\prime}(x^{k})\,\Delta x^{k}=-F(x^{k})\,.$

Convergence test: if  _$\left\|\Delta x^{k}\right\| <\epsilon$_ then

Stop. Solution found $x^{\star}=x^{k}+\Delta x^{k}$.

end if

for  _$k >0$_ do

Compute a prediction value for the damping factor
$\mu_{k}=\frac{\left\|\Delta x^{k-1}\right\|\cdot\left\|\overline{\Delta
x^{k}}\right\|}{\left\|\overline{\Delta x^{k}}-\Delta
x^{k}\right\|\cdot\left\|\Delta x^{k}\right\|}\,\cdot\,s_{k-1}\,,\quad
s_{k}=\min\left(1,\,\mu_{k}\right)\,.$

end for

Regularity test: if  _$s_{k} <s_{\min}$_ then

Stop. Convergence failure.

end if

2. Compute the trial iterate $x^{k+1}=x^{k}+s_{k}\,\Delta x^{k}$ and evaluate $F(x^{k+1})$. Solve the linear system (“old” Jacobian, “new” right-hand side):

$F^{\prime}(x^{k})\,\overline{\Delta x^{k+1}}=-F(x^{k+1})\,.$

3\.  Compute the monitoring quantities

$\Theta_{k}=\frac{\left\|\overline{\Delta x^{k+1}}\right\|}{\left\|\Delta
x^{k}\right\|}\,,\quad\mu^{\prime}_{k}=\frac{\tfrac{1}{2}\,\left\|\Delta
x^{k}\right\|\cdot s_{k}^{2}}{\left\|\overline{\Delta
x^{k+1}}-(1-s_{k})\,\Delta x^{k}\right\|}\,.$

if  _$\Theta_{k}\geq 1$_ then

Replace $s_{k}$ by
$s_{k}^{\prime}=\min\left(\mu^{\prime}_{k},\,\tfrac{1}{2}\,s_{k}\right)$. Go
to Regularity test.

else

Let $s_{k}^{\prime}=\min(1,\,\mu_{k}^{\prime})$.

end if

if  _$s_{k}^{\prime}=s_{k}=1$_ then

if  _$\left\|\overline{\Delta x^{k+1}}\right\| <\epsilon$_ then

Stop. Solution found $u^{*}=u^{k+1}+\overline{\Delta x^{k+1}}$

else if  _$s_{k}^{\prime} >4\,s_{k}$_ then

Replace $s_{k}$ by $s_{k}^{\prime}$ and go to 2.

else

Accept $x^{k+1}$ as new iterate. Continue loop with $k=k+1$

end if

end for

Algorithm 1 NLEQ-ERR

The algorithm as written in our numerical code ROTBOSON has additional
features to avoid getting stuck within infinite non-convergent iterations and
also uses a local solver, QNERR [8], once the current iteration is determined
to be sufficiently close to the final solution, as was detailed in Section VI.

## References

  * Yoshida and Eriguchi [1997] S. Yoshida and Y. Eriguchi, Phys. Rev. D (1997).
  * Lai [2004] C. Lai,  _A numerical study of boson stars_ , Ph.D. thesis, The University of British Columbia (2004).
  * Grandclément  _et al._ [2014] P. Grandclément, C. Somé, and E. Gourgoulhon, Physical Review D (2014).
  * Gourgoulhon [2010] E. Gourgoulhon, An introduction to the theory of rotating relativistic stars (2010), arXiv:1003.5015 .
  * Alcubierre [2012] M. Alcubierre,  _Introduction to 3 + 1 numerical relativity_ (Oxford University Press, 2012).
  * Ruiz  _et al._ [2008] M. Ruiz, M. Alcubierre, and D. Núñez, General Relativity and Gravitation (2008).
  * Torres [2012] J. M. Torres, AIP Conference Proceedings (2012).
  * Deuflhard [2011] P. Deuflhard,  _Newton methods for nonlinear problems affine invariance and adaptive algorithms_ (Springer, 2011).
  * Petra  _et al._ [2014a] C. G. Petra, O. Schenk, and M. Anitescu, Computing in Science Engineering (2014a).
  * Petra  _et al._ [2014b] C. G. Petra, O. Schenk, M. Lubin, and K. Gäertner, SIAM Journal on Scientific Computing (2014b).
  * Ontanon [2021a] S. Ontanon, ROTBOSON: Rotating boson stars initial data for numerical relativity, https://github.com/sontanon/ROTBOSON (2021a).
  * Liebling and Palenzuela [2017] S. L. Liebling and C. Palenzuela, Living Reviews in Relativity (2017).
  * Kaup [1968] D. J. Kaup, Phys. Rev. (1968).
  * Ruffini and Bonzzola [1969] R. Ruffini and S. Bonzzola, Phys. Rev. (1969).
  * Einstein [1915] A. Einstein, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften (Berlin), Seite 844-847. (1915).
  * Wheeler [1955] J. A. Wheeler, Phys. Rev. (1955).
  * Derrick [1964] G. H. Derrick, Journal of Mathematical Physics (1964).
  * Diez-Tejedor and Gonzalez-Morales [2013] A. Diez-Tejedor and A. X. Gonzalez-Morales, Phys. Rev. D (2013).
  * Schunck and Mielke [2003] F. E. Schunck and E. W. Mielke, Classical and Quantum Gravity (2003).
  * Ureña-López and Bernal [2010] L. A. Ureña-López and A. Bernal, Phys. Rev. D (2010).
  * Guzmán and Rueda-Becerril [2009] F. S. Guzmán and J. M. Rueda-Becerril, Phys. Rev. D (2009).
  * Barranco and Bernal [2011] J. Barranco and A. Bernal, AIP Conference Proceedings (2011).
  * Colpi  _et al._ [1986] M. Colpi, S. L. Shapiro, and I. Wasserman, Phys. Rev. Lett. (1986).
  * Friedberg  _et al._ [1987] R. Friedberg, T. D. Lee, and Y. Pang, Phys. Rev. D (1987).
  * Alcubierre  _et al._ [2018] M. Alcubierre, J. Barranco, A. Bernal, J. C. Degollado, A. Diez-Tejedor, M. Megevand, D. Núñez, and O. Sarbach, Classical and Quantum Gravity (2018).
  * Silveira and de Sousa [1995] V. Silveira and C. M. G. de Sousa, Phys. Rev. D (1995).
  * Mielke [2016] E. W. Mielke, Rotating boson stars, in  _At the Frontier of Spacetime: Scalar-Tensor Theory, Bells Inequality, Machs Principle, Exotic Smoothness_, edited by T. Asselmeyer-Maluga (Springer International Publishing, 2016) pp. 115–131.
  * Smarr and York [1978] L. Smarr and J. W. York, Phys. Rev. D (1978).
  * Lewis [1932] T. Lewis, Proceedings of the Royal Society A (1932).
  * Papapetrou [1945] A. Papapetrou, Proceedings of the Roya Irish Academy. Section A: Mathematical and Physical Sciences (1945).
  * Gustafson [1998] K. Gustafson, Contemporary Mathematics (1998).
  * Lüthen  _et al._ [2018] N. Lüthen, M. Juntunen, and R. Stenberg, Numerische Mathematik (2018).
  * Arfken and Weber [2005] G. Arfken and H. Weber,  _Mathematical Methods for Physicists_, Mathematical Methods for Physicists (Elsevier, 2005).
  * Wald [1984] R. M. Wald,  _General relativity_ (Chicago Univ. Press, Chicago, IL, 1984).
  * Alcubierre  _et al._ [2000] M. Alcubierre, S. Brandt, B. Brügmann, C. Gundlach, J. Massó, E. Seidel, and P. Walker, Classical and Quantum Gravity (2000).
  * Kollerstrom [1992] N. Kollerstrom, The British Journal for the History of Science (1992).
  * Ypma [1984] T. J. Ypma, SIAM Journal on Numerical Analysis (1984).
  * Bank and Rose [1981] R. E. Bank and D. J. Rose, Numerische Mathematik (1981).
  * LeVeque [2007] R. LeVeque,  _Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems_, Other Titles in Applied Mathematics (Society for Industrial and Applied Mathematics, 2007).
  * Kantorovich and Akilov [1982] L. V. Kantorovich and G. P. Akilov,  _Functional analysis_ (Pergamon, 1982).
  * Ortega and Rheinboldt [2000] J. M. Ortega and W. C. Rheinboldt,  _Iterative solution of nonlinear equations in several variables_ (Society for Industrial and Applied Mathematics, 2000).
  * Deuflhard [1975] P. Deuflhard, Lecture Notes in Mathematics Optimization and Optimal Control (1975).
  * Cauchy [1847] A.-L. Cauchy, Oeuvres complètes (1847).
  * Levenberg [1944] K. Levenberg, Quarterly of Applied Mathematics (1944).
  * Marquardt [1963] D. W. Marquardt, Journal of the Society for Industrial and Applied Mathematics (1963).
  * Nowak and Weimann [1991] U. Nowak and L. Weimann,  _A Family of Newton Codes for Systems of Highly Nonlinear Equations_, Tech. Rep. (Zuse Institute Berlin, 1991).
  * Dagum and Menon [1998] L. Dagum and R. Menon, IEEE Comput. Sci. Eng. (1998).
  * Torres [2016] J. M. Torres,  _Dinámica de materia cargada en relatividad general_ , Ph.D. thesis, Universidad Nacional Autónoma de México, Instituto de Ciencias Nucleares (2016).
  * Ontanon [2018] S. Ontanon, Resolvedor elíptico para relatividad numérica y su aplicación en ondas de Brill, Bachelor’s Thesis (2018).
  * Corporation [2018] I. Corporation, Intel mkl pardiso - parallel direct sparse solver interface (2018), online.
  * Tewarson [1973] R. P. Tewarson, ed.,  _Sparse Matrices_, Mathematics in Science and Engineering (Elsevier Science, 1973).
  * Gupta and Kumar [1995] A. Gupta and V. Kumar, in  _Proceedings of the IEEE/ACM SC95 Conference_ (1995).
  * Inc. [2020] W. R. Inc., Mathematica, Version 12.2 (2020), champaign, IL, 2020.
  * Schenk and Gärtner [2020] O. Schenk and K. Gärtner,  _PARDISO User Guide_, PARDISO Project (2020).
  * Forum [1994] M. P. Forum,  _MPI: A Message-Passing Interface Standard_, Tech. Rep. (University of Tennessee, Knoxville, TN, USA, 1994).
  * Richardson [1911] L. F. Richardson, Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character (1911).
  * Ontanon [2021b] S. Ontanon, SPHBOSON: Spherical boson stars initial data for numerical relativity, https://github.com/sontanon/SPHBOSON (2021b).