# Boson stars and their relatives in semiclassical gravity

Miguel Alcubierre  Instituto de Ciencias Nucleares, Universidad Nacional
Autónoma de México, Circuito Exterior C.U., A.P. 70-543, México D.F. 04510,
México  Juan Barranco  Departamento de Física, División de Ciencias e
Ingenierías, Campus León, Universidad de Guanajuato, León 37150, México
Argelia Bernal  Departamento de Física, División de Ciencias e Ingenierías,
Campus León, Universidad de Guanajuato, León 37150, México  Juan Carlos
Degollado  Instituto de Ciencias Físicas, Universidad Nacional Autónoma de
México, Apdo. Postal 48-3, 62251, Cuernavaca, Morelos, México  Alberto Diez-
Tejedor  Departamento de Física, División de Ciencias e Ingenierías, Campus
León, Universidad de Guanajuato, León 37150, México  Miguel Megevand
Instituto de Física Enrique Gaviola, CONICET. Ciudad Universitaria, 5000
Córdoba, Argentina  Darío Núñez  Instituto de Ciencias Nucleares, Universidad
Nacional Autónoma de México, Circuito Exterior C.U., A.P. 70-543, México D.F.
04510, México  Olivier Sarbach  Instituto de Física y Matemáticas, Universidad
Michoacana de San Nicolás de Hidalgo, Edificio C-3, Ciudad Universitaria,
58040 Morelia, Michoacán, México

###### Abstract

We construct boson star configurations in quantum field theory using the
semiclassical gravity approximation. Restricting our attention to the static
case, we show that the semiclassical Einstein-Klein-Gordon system for a single
real quantum scalar field whose state describes the excitation of $N$
identical particles, each one corresponding to a given energy level, can be
reduced to the Einstein-Klein-Gordon system for $N$ complex classical scalar
fields. Particular consideration is given to the spherically symmetric static
scenario, where energy levels are labeled by quantum numbers $n$, $\ell$ and
$m$. When all particles are accommodated in the ground state $n=\ell=m=0$, one
recovers the standard static boson star solutions, that can be excited if
$n\neq 0$. On the other hand, for the case where all particles have fixed
radial and total angular momentum numbers $n$ and $\ell$, with $\ell\neq 0$,
but are homogeneously distributed with respect to their magnetic number $m$,
one obtains the $\ell$-boson stars, whereas when $\ell=m=0$ and $n$ takes
multiple values, the multi-state boson star solutions are obtained. Further
generalizations of these configurations are presented, including the
multi-$\ell$ multi-state boson stars, that constitute the most general
solutions to the $N$-particle, static, spherically symmetric, semiclassical
real Einstein-Klein-Gordon system, in which the total number of particles is
definite. In spite of the fact that the same spacetime configurations also
appear in multi-field classical theories, in semiclassical gravity they arise
naturally as the quantum fluctuations associated with the state of a single
field describing a many-body system. Our results could have potential impact
on direct detection experiments in the context of ultralight scalar
field/fuzzy dark matter candidates.

###### pacs:

95.30.Sf, 95.35.+d, 04.62.+v, 98.80.Jk

##  I Introduction

Boson stars are exotic objects made of bosons in which the gravitational force
that pulls matter together is counterbalanced by the dispersive nature of a
scalar field. They were first proposed in the late 1960s by Kaup Kaup68  and
Ruffini and Bonazzola Ruffini:1969qy , and since then they have been actively
studied for more than half a century Colpi:1986ye ; Friedberg87 ;
Gleiser:1988rq ; Lee:1988av ; Bizon:2000 ; Guzman:2004wj ; Chavanis:2011cz
—see Jetzer:1991jr ; Schunck:2003kk ; Liebling:2012fv ; Zhang:2018slz ;
Visinelli:2021uve  for reviews on boson stars and Lee:1991ax ; Rajaraman ;
Manton2007 ; EJWeinberg ; Shnir  for references on other soliton solutions. At
present, boson stars remain largely theoretical, although they have been
employed to describe dark compact objects Torres:2000dw ; Guzman:2005bs ;
Cardoso:2019rvt ; Hogan:1988mp ; Kolb:1993zz ; Palenzuela:2017kcg  and
galactic halo cores Sin:1992bg ; Schive:2014dra ; Schive:2014hza ;
Schwabe:2016rze ; Veltmaat:2016rxo ; Du:2016aik ; Gonzalez-Morales:2016yaf  in
models of axion Peccei77 ; Preskill83 ; Abbott83 ; Dine83 ; Hertzberg08 ;
Sikivie08  and axion-like Hu:2000ke ; Matos:2000ss ; Arvanitaki:2009fg ;
Suarez:2013iw ; Marsh:2015xka ; Hui:2016ltb ; Niemeyer:2019aqm ; Luis2019 ;
Ferreira:2020fam  particles.

More pragmatically, a boson star is a regular, localized solution to the
classical Einstein-Klein-Gordon (EKG) system. Nevertheless, nature is quantum
at a fundamental level, and as such these objects must also allow an
interpretation in quantum field theory. The purpose of this paper is to
construct boson star configurations in semiclassical gravity, to catalog their
spectrum of spherically symmetric equilibrium solutions, and to compare them
with those of the classical theory. Previous attempts to construct
semiclassical boson stars have been carried out in Ruffini:1969qy ; Ho:2002vz
; Matos2007 ; Bernal:2009zy ; UrenaLopez:2010ur ; Barranco:2010ib . See also
Refs. Berczi:2020nqy ; Berczi:2021hdh ; Guenther:2020kro  for an analysis of
the semiclassical gravitational collapse of quantum matter and
arrechea2021semiclassical  for a recent study on fluid stars in semiclassical
gravity.

The semiclassical theory of gravity is an effective description of
gravitational phenomena that deals with gravitons at tree level and with
matter fields at one loop Verdaguer2020 ; semiclassical ; Jordan:1986ug ;
Calzetta:1986ey ; Paz:1990jg ; Campos:1993ug . The resulting equations of
motion are those of quantum field theory on curved spacetimes coupled to the
semiclassical Einstein equations, where the expectation value of the stress
energy-momentum tensor operator acts as the source term on its right-hand
side. In this article, we introduce a well-defined operational program that
deals with free quantum fields acting as sources of stationary spacetimes. Our
program relies on the semiclassical self-consistent configurations proposed in
Ref. DiezTejedor2012 , and it can be summarized into three steps:

i) Consider a stationary, globally hyperbolic background spacetime on which
the free quantum fields are defined. The assumption of stationarity allows one
to introduce a preferred space of “positive-norm” solutions of the matter
field equations, hence a preferred vacuum state. This in turn provides a well-
defined theory for the quantum fields, which we describe in terms of a Fock
space representation. In particular, field operators can be written (formally)
as linear combinations of creation and annihilation operators, wherein the
“coefficients” are mode functions $f_{I}(x)$ that solve the classical complex
field equations.

ii) Compute the expectation value $\langle\hat{T}_{\mu\nu}\rangle$ of the
stress energy-momentum tensor operator with respect to a given state in the
Fock space. In order to do so we need a regularization and renormalization
prescription that removes the ill-defined ultraviolet behavior of the theory,
leading to sensible finite outcomes. To achieve this, in this work we impose
normal ordering. More sophisticated approaches include e.g. adiabatic
subtraction Parker1974  and Pauli-Villars renormalization Pauli1940 ;
Weinberg:2010wq ; Armendariz-Picon:2019csc , although we expect the
differences between such methods and ours to be suppressed in the limit of
large occupation numbers, as we consider in our configurations, which we also
assume to be far from the Planck scale. More generally, we can also compute a
statistical average by tracing $\hat{T}_{\mu\nu}$ with a density operator.
This offers the interesting possibility of considering, for instance, thermal
states with a given temperature.

iii) Solve the semiclassical Einstein equations $G_{\mu\nu}=8\pi
G\langle\hat{T}_{\mu\nu}\rangle$ sourced by the expectation value (or
statistical average) of the (renormalized) stress energy-momentum tensor. This
step takes into account the backreaction of the quantum fields on the
classical geometry.

Of course, one of the main difficulties of this approach is that the exact
spacetime geometry is not known a priori in step i, and it has to be
constructed in a self-consistent way together with the other two steps. For
the purpose of illustration we will restrict our attention to the case in
which matter consists of a single, free, minimally coupled real scalar field,
although more involved situations can also be explored, including quantum
fields of higher rank. (The case of a complex scalar field is analyzed in an
appendix.) In order to simplify the analysis we shall further concentrate on
static, spherically symmetric configurations, but a generalization of our
formalism to describe stationary and axisymmetric rotating objects should be
possible. In particular, as we show, in the static case the semiclassical EKG
equations can be reduced to a system of self-gravitating classical complex
scalar fields with harmonic time-dependency of the form $e^{-i\omega_{I}t}$,
which leads to a non-linear multi-eigenvalue problem for the frequencies
$\omega_{I}$. These “classical” fields arise from the mode functions
$f_{I}(x)$ that appear in the decomposition of the field operator and
represent the “wave functions” of individual particles in first quantization.
Such problems were treated in the Newtonian limit long ago; see, for example
the seminal work by Lieb Lieb77  or Ref. KavianMuschler2015  for more recent
work in this direction.

In the static case, the resulting semiclassical solutions can be interpreted
as describing equilibrium selfgravitating objects made of bosons.
Specifically, we construct general boson star configurations in spherical
symmetry, for which the number of particles in the different energy levels is
definite. These objects interpolate between standard boson stars Kaup68 ;
Ruffini:1969qy ; Colpi:1986ye ; Friedberg87 ; Gleiser:1988rq ; Lee:1988av ;
Guzman:2004wj ; Bizon:2000 ; Chavanis:2011cz ; Jetzer:1991jr ; Schunck:2003kk
; Liebling:2012fv ; Zhang:2018slz ; Visinelli:2021uve , whose particles all
lie in the lowest possible energy configuration, to more general situations
where the particles are accommodated in states with higher energy and angular
momentum, which include $\ell$-boson stars Alcubierre:2018ahf ;
Alcubierre:2019qnh ; Alcubierre:2021mvs ; Alcubierre:2021psa ;
Jaramillo:2022zwg  and multi-state boson stars Matos2007 ; Bernal:2009zy ;
UrenaLopez:2010ur , as well as new configurations obtained in this article:
multi-$\ell$ multi-state boson stars. As we show, these constitute the most
general solutions to the static, spherically symmetric, semiclassical real EKG
system for which the total number of particles is definite (we shall refer to
these configurations as $N$-particle systems in this paper), and encompasses
the previous boson star solutions reported in the literature. A family tree of
these solutions is provided in Table 2, where we show how they are connected
to each other. A relevant question is whether or not such configurations in
which particles populate not only the ground state, but also higher energy
levels, are stable. For instance, it has been found that boson stars with
$\ell=0$ including only excited states are unstable Lee:1988av . However, as
discussed in Matos2007 ; Bernal:2009zy ; UrenaLopez:2010ur , a possible
mechanism of stabilization is to have a suitable combination of particles in
the ground and higher energy states.111See also DiGiovanni:2021vlu  for a
stabilization mechanism that includes another type of matter. These multi-
state configurations arise naturally within our semiclassical description. See
jaramillo19 ; nambo21  for $\ell$-boson stars in the Newtonian limit and also
Vicens:2018kdk ; Lin:2018whl  for related configurations which include
particles in the excited states.

Before we proceed with the construction of these objects, some words are
needed regarding the relation between the classical and quantum descriptions
of boson stars. Classical fields emerge from quantum theories in the limit
when the quantum fluctuations become negligible. This limit is manifest, for
instance, in the case of coherent states, that saturate the quantum
uncertainty principle and lead to field configurations in which quantum
fluctuations are reduced to their minimum. However, in the real scalar field
theory coherent states are not compatible with a static spacetime geometry, as
we prove later. This is not surprising, given that for a real scalar field
there are no static configurations in the classical theory Seidel91 . This can
be traced back to the properties of the classical limit itself. The existence
of soliton solutions relies on the presence of conserved charges
Coleman:1985ki ; Lee:1991ax  which allow localized field configurations whose
energy per unit charge is less than in any other solutions, including those in
which all the charge is radiated to infinity. In the classical field theory
for a real scalar field there is no such charge, which explains the absence of
static solutions associated with the coherent states.

The situation is different in the quantum theory, where the conserved charge
is the particle number that remains constant if the configuration is
static.222As we show in appendix A, if the spacetime is static, the particle
number operator commutes with the Hamiltonian of the system, which is the
generator of the time translations. Interestingly, despite the fact that the
$N$-particle solutions describe many-body systems and appear as a result of
quantum fluctuations (and, consequently, lie beyond the classicality of boson
stars based on path integral arguments discussed in, e.g., reference
Herdeiro:2022gzp ), they still have a counterpart in multi-field classical
theories. This is due to the relation that exists between the semiclassical
EKG system describing a single real scalar field in a quantum state with a
definite number of particles $N$ and the EKG system for $N$ complex classical
scalar fields. In this way, we show that boson stars and their relatives
(i.e., the multi-$\ell$ multi-state boson stars) can be understood within our
program by invoking a single real quantum scalar field without the need of
postulating the existence of a fixed number of independent complex classical
scalar fields, like, for example, the number $2\ell+1$ in the construction of
$\ell$-boson stars as originally required in Alcubierre:2018ahf .
Nevertheless, the interpretation of the solutions is different in the
classical and the quantum limits, and the difference can be found in the role
that the mode functions $f_{I}(x)$ play in the different regimes of the
theory. On one side, in an $N$-particle state the resulting complex fields
$f_{I}(x)$ are understood as the wave functions of the individual particles in
first quantization, and they represent the many-body Hartree approximation
Hartree1928  of a system of $N$ particles that live in the mean gravitational
field that they produce, where the $N$-particle wave function is just the
product of one-particle wave functions Guth:2014hsa . On the other hand, if
the state of the quantum field is coherent, the mode functions play the role
of a classical field excitation
$\sum_{I}[\alpha_{I}f_{I}(x)+\alpha_{I}^{*}f_{I}^{*}(x)]$, and a description
in terms of particles is not appropriate in this case, in the same way that a
description in terms of photons is not suitable in classical electrodynamics.
Table 1 sketches the connection between the classical and the quantum regimes.
Bearing in mind the different regimes of the theory may be relevant for
potential direct detection experiments, such as those carried out in
Graham2013 ; Carney:2019cio ; Donohue:2021jbv ; Tsai:2021lly . References
Perez:2005gh ; Sudarsky:2009za ; Aguirre:2015mva ; Sikivie:2016enz ;
Hertzberg:2016tal ; Allali:2020shm  delve on the discussion of the classical
and the quantum regimes of a scalar field in different cosmological and
astrophysical situations.

This paper is organized as follows. In section II we review the main
ingredients of quantum field theory on curved spaces and semiclassical
gravity. In section III we focus on the static case and next, in section IV we
further specialize to the static spherically symmetric situation. This leads
to the main theoretical result of this article, which is summarized in the
semiclassical EKG system of Eqs. (40), (44a), and (44c). Remarkably, as a
consequence of the semiclassical approach the resulting system of equations
includes as particular case the system for $N$ classical complex fields. This
constitutes the starting point for the subsequent analysis of this paper.
Numerical solutions presenting new configurations which arise naturally in our
formalism, including multi-$\ell$, multi-state, and multi-$\ell$ multi-state
boson stars are presented in section V. Conclusions are drawn in section VI
and technical aspects of our calculations are included in appendix A. In
appendix B we introduce the static, spherically symmetric, semiclassical
complex EKG system which, in addition to the $N$-particle configurations,
allows solutions sourced by coherent states that are not static.

Our conventions are as follows. We use the mostly plus signature convention
for the spacetime metric, $(-,+,+,+)$, and to simplify the notation, we work
in terms of natural units for which $\hbar=c=1$. Numerical results are
obtained using Planck units, where in addition we set $G=1$.

Table 1: Minimal ingredients to construct boson stars and their relatives in the different regimes of a quantum scalar field theory in the semiclassical gravity approximation. In this paper we concentrate on static configurations. When the quantum field is in a coherent state, the mode functions are related with the excitation of a classical field. In contrast, when the quantum field is in a $N$-particle state, the mode functions are related to the particle wave functions describing a many-body system. Regime |  | State |  | Boson stars |  | Multi-$\ell$ multi-state boson stars
---|---|---|---|---|---|---
classical field excitation |  | coherent state |  | one complex field |  |  $N$ complex fields
many-body system |  |  $N$-particle state |  | one real/complex field |  | one real/complex field

##  II Theoretical framework

The quantization of a free field on a curved, globally hyperbolic spacetime is
well understood. This program was initiated by Parker in the late 1960s and
developed further by Fulling, Ford and Wald, among others (see e.g.
Birrell1984 ; Wald1994 ; Mukhanov2007 ; Parker2009  for relevant textbooks on
this subject and references to the aforementioned original work). If in
addition the quantum fields act as a source of the spacetime metric, the
semiclassical theory of gravity Verdaguer2020 ; semiclassical  provides an
effective description that combines the quantum nature of matter with the
classical behavior that gravity exhibits at macroscopic scales.333The regime
of applicability of semiclassical gravity is an open question, due mainly to
the fact that we do not have access to a complete, satisfactory theory of
quantum gravity. See the discussion in Ref. Flanagan:1996gw , section IIA (and
references therein), for a critical examination of the origin of the
semiclassical equations. In this section we review the main ingredients of
this construction.

###  II.1 Quantum spin-0 fields in curved spaces

For the following, we consider a globally hyperbolic spacetime
$(\mathcal{M},ds^{2})$ which is foliated by three-dimensional Cauchy
hypersurfaces $\Sigma_{t}$. In terms of the standard $3+1$ decomposition the
spacetime metric is written as

$ds^{2}=-(\alpha^{2}-\beta_{j}\beta^{j})dt^{2}+2\beta_{i}dtdx^{i}+\gamma_{ij}dx^{i}dx^{j}.$
(1)

Here, $\alpha(x)$ is the lapse function, $\beta^{j}(x)$ is the shift vector,
and $\gamma_{ij}(x)$ is the induced metric on $\Sigma_{t}$, with
$x=(t,\vec{x})$ denoting a generic point in the spacetime manifold. Latin
indices $i,j,k,\ldots$ take natural values in the range from $1$ to $3$ and
are raised and lowered with the three-metric $\gamma_{ij}$, e.g.,
$\beta_{i}=\gamma_{ij}\beta^{j}$.

At the classical level, a real free massive scalar field satisfies the Klein-
Gordon equation

$(\Box-{m_{0}}^{2})\,\phi=0,$ (2)

where $\Box:=g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$ is the curved d’Alembertian
operator in four dimensions, $g^{\mu\nu}$ is the inverse of the spacetime
metric, and $\nabla_{\mu}$ is the covariant derivative with respect to this
metric. The parameter $m_{0}$, which we assume to be positive, denotes the
inverse Compton length of the field (that plays the role of the rest mass of
the particles in the quantum theory), and for simplicity a minimal coupling
with gravity has been considered.

For the quantization of the field $\phi$, one extends the space of real
classical solutions to the space of complex-valued classical solutions of Eq.
(2). Let us call this space $X$ in the following. Given two such solutions
$\phi_{1},\phi_{2}\in X$, one introduces the four-current vector field (with
$\nabla^{\mu}:=g^{\mu\nu}\nabla_{\nu}$)

$j^{\mu}(\phi_{1},\phi_{2}):=-i\left[\phi_{1}(\nabla^{\mu}\phi_{2}^{*})-(\nabla^{\mu}\phi_{1})\phi_{2}^{*}\right],$
(3)

which, by virtue of Eq. (2), is divergence-free ($\nabla_{\mu}j^{\mu}=0$) and
satisfies the symmetries
$j^{\mu}(\phi_{1},\phi_{2})=[j^{\mu}(\phi_{2},\phi_{1})]^{*}=-j^{\mu}(\phi_{2}^{*},\phi_{1}^{*})$.
Here and in the following, $\phi_{2}^{*}(x)$ denotes the complex conjugate of
$\phi_{2}(x)$. The four-current (3) gives rise to an inner product on $X$,
defined as

$(\phi_{1},\phi_{2}):=\int_{\Sigma_{t}}j^{\mu}n_{\mu}d\gamma=-i\int_{\Sigma_{t}}[\phi_{1}(\pounds_{n}\phi_{2}^{*})-(\pounds_{n}\phi_{1})\phi_{2}^{*}]d\gamma.$
(4)

In this equation $(n_{\mu})=(-\alpha,0,0,0)$ is the future-directed time-like
unit normal covector field to the Cauchy hypersurfaces $\Sigma_{t}$,
$\pounds_{n}\phi_{2}=n^{\mu}\nabla_{\mu}\phi_{2}$ refers to the Lie derivative
of $\phi_{2}$ with respect to the corresponding vector field
$n=g^{\mu\nu}n_{\mu}\partial_{\nu}$, and
$d\gamma=\sqrt{\det(\gamma_{ij})}\,d^{3}x$ denotes the volume element on this
hypersurface. As long as the space $X$ is restricted to those solutions of Eq.
(2) which decay sufficiently fast at spatial infinity, the inner product (4)
does not depend on the choice of the Cauchy hypersurface. Note also that by
construction the inner product (4) is linear in its first argument and
inherits the symmetries of the four-current, such that
$(\phi_{1},\phi_{2})=(\phi_{2},\phi_{1})^{*}=-(\phi_{2}^{*},\phi_{1}^{*})$.
However, it fails to be positive definite. Indeed, for
$\phi_{1}=\phi_{2}=:\phi\in X$,

$(\phi,\phi)=2\,\mbox{Im}\int_{\Sigma_{t}}\phi(\pounds_{n}\phi^{*})d\gamma$
(5)

may assume any real (positive or negative) value, since the restrictions of
the functions $\phi$ and $\pounds_{n}\phi$ on $\Sigma_{t}$ represent the
Cauchy data for Eq. (2), which is free.

At the quantum level, the scalar field and its conjugate momentum
$\pi(x):=\sqrt{\det(\gamma_{ij})}\pounds_{n}\phi(x)$ are promoted to self-
adjoint field operators $\hat{\phi}(x)$ and
$\hat{\pi}(x)=\sqrt{\det(\gamma_{ij})}\pounds_{n}\hat{\phi}(x)$ acting on an
abstract Hilbert space $\mathscr{H}$. These operators satisfy the standard
equal time commutation relations

$[\hat{\phi}(t,\vec{x}),\hat{\pi}(t,\vec{y})]=i\delta^{(3)}(\vec{x}-\vec{y}),\qquad[\hat{\phi}(t,\vec{x}),\hat{\phi}(t,\vec{y})]=[\hat{\pi}(t,\vec{x}),\hat{\pi}(t,\vec{y})]=0,$
(6)

for all $(t,\vec{x}),(t,\vec{y})\in\Sigma_{t}$ and $\hat{\phi}(x)$ satisfying
the Klein-Gordon equation (2). Note that we follow a canonical quantization
scheme, and we are working in the Heisenberg representation, where the
evolution is codified in the operators and the state-vectors remain
independent of time. Instead of the field operator $\hat{\phi}(x)$, which is
really an operator-valued distribution on ${\cal M}$, it is sometimes
convenient to work with its “smeared-out” versions, given by the operators

$\hat{a}(f):=(\hat{\phi},f),\quad f\in X,$ (7)

with $(\cdot,\cdot)$ the same inner product as in Eq. (4), such that the
definition is again independent of the choice of the Cauchy surface and
$\hat{a}(f)$ is constant in time. Because $\hat{\phi}=\hat{\phi}^{\dagger}$ is
self-adjoint, it follows that
$\hat{a}^{\dagger}(f)=-(\hat{\phi},f^{*})=-\hat{a}(f^{*})$ and the commutation
relations (6) imply

$[\hat{a}(f),\hat{a}(g)^{\dagger}]=(g,f),\qquad[\hat{a}(f),\hat{a}(g)]=-(g^{*},f),$
(8)

for all $f,g\in X$.

For the following, we make the important assumption that the space of complex-
valued, classical solutions $X$ of the Klein-Gordon equation (2) can be
decomposed in the form Jacobson:2003vx

$X=X_{+}\oplus X_{+}^{*},$ (9)

with the subspace $X_{+}$ consisting of solutions with positive norm (that is,
$(f,f)>0$ for all $f\in X_{+}$ with $f\neq 0$) and its complex conjugate
$X_{+}^{*}$ being orthogonal to it, such that $(f,g^{*})=0$ for all $f,g\in
X_{+}$.444Given the properties of the inner product, the elements of
$X_{+}^{*}$ have negative norm; however, this will not be relevant for what
follows. For a detailed discussion on the validity and uniqueness of this
decomposition we refer the reader to Ashtekar:1975zn . For a generic spacetime
manifold $({\cal M},ds^{2})$, it does not seem clear if a split of this kind
exists and is unique; however, for static or stationary spacetimes, i.e.,
those admitting a globally-defined timelike Killing vector field, the
decomposition of $X$ exists, and, moreover, the “energy-requirement” of Ref.
Ashtekar:1975zn  selects a preferred one. Under these assumptions, the vacuum
state (which in general depends on the choice of the decomposition) can be
characterized as the state $|0\rangle\in\mathscr{H}$ for which $\langle
0|0\rangle=1$ and $\hat{a}(f)|0\rangle=0$ for all $f\in X_{+}$ in the
preferred decomposition. The particular case for which $({\cal M},ds^{2})$ is
static will be reviewed in the next section. In this case, the natural choice
for $X_{+}$ satisfying the energy-requirement can be constructed directly from
the space of “positive-frequency” solutions of the Klein-Gordon equation (2).

To proceed, it is convenient to work with an orthonormal set of basis
functions $f_{1},f_{2},\ldots\in X_{+}$, such that
$(f_{I},f_{J})=\delta_{IJ}$, which are usually refereed to as the mode
functions, and to introduce the corresponding creation and annihilation
operators

$\hat{a}_{I}:=\hat{a}(f_{I}),\qquad\hat{a}_{I}^{\dagger}:=\hat{a}^{\dagger}(f_{I}),$
(10)

which, by virtue of Eq. (8) and the decomposition (9) fulfill the commutation
relations

$[\hat{a}_{I},\hat{a}^{\dagger}_{J}]=\delta_{IJ},\qquad[\hat{a}_{I},\hat{a}_{J}]=0.$
(11)

In terms of the operators $\hat{a}_{I}$ and $\hat{a}_{I}^{\dagger}$, the field
operator $\hat{\phi}(x)$ can be decomposed as

$\hat{\phi}(x)=\sum_{I}\left[\hat{a}_{I}f_{I}(x)+\hat{a}^{\dagger}_{I}f^{*}_{I}(x)\right],$
(12)

which allows one to disentangle the field properties, codified in the
spacetime functions $f_{I}(x)$ and $f_{I}^{*}(x)$, from the time-independent
quantum operators $\hat{a}_{I}$ and $\hat{a}_{I}^{\dagger}$. Notice that this
decomposition is not unique, and any choice of the orthonormal set of basis
functions works equally well.

The Hilbert space $\mathscr{H}$ can now be constructed (à la Fock) by
successive applications of creation operators on the vacuum state. A generic
element in the base of the Fock construction can be written in the form

$|N_{1},N_{2},\ldots\rangle=\frac{(\hat{a}^{\dagger}_{1})^{N_{1}}}{\sqrt{N_{1}!}}\frac{(\hat{a}^{\dagger}_{2})^{N_{2}}}{\sqrt{N_{2}!}}\ldots|0\rangle,$
(13)

with $N_{1},N_{2},\ldots$ nonnegative integer numbers such that
$\sum_{I}N_{I}$ is finite. Using the commutation relations (11), one easily
shows that the basis vectors (13) are normalized and mutually orthogonal and
that

$\displaystyle\hat{a}_{K}^{\dagger}|N_{1},\ldots N_{K},\ldots\rangle$
$\displaystyle=$
$\displaystyle\sqrt{N_{K}+1}|N_{1},\ldots,N_{K}+1,\ldots\rangle,$ (14a)
$\displaystyle\hat{a}_{K}|N_{1},\ldots,N_{K},\ldots\rangle$ $\displaystyle=$
$\displaystyle\sqrt{N_{K}}|N_{1},\ldots,N_{K}-1,\ldots\rangle.$ (14b)

Furthermore, the states (13) are eigenvectors of the particle number operator
$\hat{N}_{K}:=\hat{a}^{\dagger}_{K}\hat{a}_{K}$, such that
$\hat{N}_{K}|N_{1},\ldots,N_{K},\ldots\rangle=N_{K}|N_{1},\ldots,N_{K},\ldots\rangle$,
with $N_{K}$ representing the number of particles in the K-th mode, and thus
the states (13) describe a system of $N=\sum_{I}N_{I}$ identical quantum
particles, with $N_{1}$ of them in the one-particle state corresponding to
mode $1$, $N_{2}$ of them in the state corresponding to mode $2$, and so on.
Note that the expectation value of the field operator $\hat{\phi}(x)$ vanishes
when evaluated on a state with a definite number of particles, i.e., $\langle
N_{1},N_{2},\ldots|\hat{\phi}(x)|N_{1},N_{2},\ldots\rangle=0$, although this
does not imply that the expectation value of the stress energy-momentum tensor
also vanishes, as we will see later. Because the creation operators
$\hat{a}_{I}^{\dagger}$ commute with each other, these states are totally
symmetric and thus the particles satisfy the Bose-Einstein statistics and
describe bosons.

An arbitrary (pure) state in the Hilbert space $\mathscr{H}$ can be expressed
as a linear combination of the elements in the Fock construction,

$|\psi\rangle=\sum_{N_{1},N_{2},\ldots=0}^{\infty}(C_{N_{1}N_{2}\ldots})|N_{1},N_{2},\ldots\rangle,$
(15)

with $C_{N_{1}N_{2}\ldots}$ arbitrary complex numbers such that
$\sum_{N_{1},N_{2},\ldots=0}^{\infty}|C_{N_{1}N_{2}\ldots}|^{2}=1$. A case of
particular interest consists of the coherent states, defined as those elements
of $\mathscr{H}$ that saturate the quantum uncertainty principle and most
closely resemble a classical field excitation; see, e.g., page 97 of Ref.
Sakurai94  for a brief description of the coherent states in the context of
the single-particle quantum harmonic oscillator and Zhang1990 ; Sanders2012
for comprehensive reviews. In the context of a field theory, they are usually
referred as Glauber states Glauber63  and are defined as the eigenstates of
the (non-Hermitian) annihilation operators $\hat{a}_{I}$,

$\hat{a}_{K}|\alpha_{1},\alpha_{2},\ldots\rangle=\alpha_{K}|\alpha_{1},\alpha_{2},\ldots\rangle,$
(16)

with $\alpha_{K}$ being in general complex numbers. Note that, contrary to
what happens for the states with a definite number of particles, the
expectation value of the scalar field does not vanish when evaluated on a
coherent state, where we obtain
$\langle\alpha_{1},\alpha_{2},\ldots|\hat{\phi}(x)|\alpha_{1},\alpha_{2},\ldots\rangle=\sum_{I}[\alpha_{I}f_{I}(x)+\alpha_{I}^{*}f_{I}^{*}(x)]$,
which is a solution to the classical Klein-Gordon equation (2). This is not
surprising and actually is a consequence of Ehrenfest’s theorem and the fact
that we are dealing with a linear theory, so the expectation value of the
field operator $\hat{\phi}(x)$ always satisfies the classical equations of
motion.

###  II.2 Semiclassical gravity

So far, we have ignored the backreaction of the quantum fields, and we have
assumed that the spacetime background is given a priori. However, according to
general relativity, the spacetime metric is determined dynamically by the
distribution of matter through Einstein’s field equations, for which a
(classical) stress energy-momentum tensor is required. One possibility to
address this problem is to follow an effective field theory approach where,
starting from the generating functional $Z[J,T^{\mu\nu}]$ (with $J$ and
$T^{\mu\nu}$ external sources of $\phi$ and $g_{\mu\nu}$), one expands the
effective action $\Gamma[\phi,g]$ at tree level in gravitons and one loop in
matter fields Jordan:1986ug ; Calzetta:1986ey ; Paz:1990jg ; Campos:1993ug
(see also Ref. Armendariz-Picon:2020tkc  for a derivation of the quantum
corrected equations of motion of the metric in terms of an analysis of
graviton fluctuations). The resulting theory is known as semiclassical
gravity, where, in addition to the quantum field theory summarized in section
II.1, one enforces Einstein’s equations sourced by the expectation value of
the stress energy-momentum tensor Verdaguer2020 ; semiclassical ,

$G_{\mu\nu}=8\pi G\langle\hat{T}_{\mu\nu}\rangle.$ (17)

Here, $G_{\mu\nu}$ is the Einstein tensor, and
$\langle\hat{T}_{\mu\nu}\rangle=\langle\psi|\hat{T}_{\mu\nu}|\psi\rangle$
denotes the expectation value of the stress energy-momentum operator when
evaluated on an arbitrary state $|\psi\rangle$ in $\mathscr{H}$. Further
details on how to solve this problem based on the notion of semiclassical
self-consistent configurations introduced in Ref. DiezTejedor2012  will be
given below. For recent rigorous results on the initial-value problem for
semiclassical gravity, see, for instance, references Juarez-Aubry:2019jon ;
bjuarez2021a ; bjuarez2021b ; Juarez-Aubry:2022qdp .

For the case of a real free massive scalar field, the operator associated with
the stress energy-momentum tensor takes the form

$\hat{T}_{\mu\nu}=(\nabla_{\mu}\hat{\phi})(\nabla_{\nu}\hat{\phi})-\frac{1}{2}g_{\mu\nu}\left[(\nabla_{\alpha}\hat{\phi})(\nabla^{\alpha}\hat{\phi})+{m_{0}}^{2}\hat{\phi}\hat{\phi}\right].$
(18)

It is important to notice that this quantity is quadratic in field operators
and that it contains products of these operators evaluated at the same
spacetime point. Since $\hat{\phi}$ is a distribution, these products are not
mathematically well-defined, and this problem manifests itself as divergences
when computing the right-hand side of Eq. (17). Some regularization and
renormalization prescription is needed in order to subtract the ill-defined
ultraviolet behavior from the expectation value of higher order operators,
providing sensible finite results. On the one hand, this requires the
introduction of counterterms into the effective action, in such a way that the
divergences that appear in the free theory are absorbed into the cosmological
constant, Newton’s gravitational constant, and the coupling constants
accompanying quadratic curvature scalars such as $R^{2}$ and
$R_{\mu\nu}R^{\mu\nu}$ semiclassical  (which are expected to be suppressed in
the low energy regime and we do not include here).555The observed value of the
cosmological constant is so small that its relevance at local scales is
negligible and for that reason we will not include this term in our analysis
either. On the other hand, this also leads to a finite contribution to the
expectation value of the stress energy-momentum tensor originating from the
structure of the vacuum itself, that even if interesting in its own right,
will not be explored in more detail in the present paper (this contribution is
expected to be suppressed for large occupation numbers, and this is what we
assume in the following.) In practice, this corresponds to assuming normal
(Wick) ordering and writing, e.g.,
$:\\!\hat{a}_{I}\hat{a}_{I}^{\dagger}\\!:\;=\hat{a}_{I}^{\dagger}\hat{a}_{I}$
in our expressions, moving all the creation operators to the left.

Introducing the field decomposition (12) in terms of the creation and
annihilation operators into the expression for the stress energy-momentum
tensor (18), one obtains

$\hat{T}_{\mu\nu}=\frac{1}{2}\sum_{I,J}\left[\hat{a}_{I}\hat{a}_{J}T_{\mu\nu}(f_{I},f_{J})+\hat{a}^{\dagger}_{I}\hat{a}_{J}T_{\mu\nu}(f_{I}^{*},f_{J})+\textrm{H.c.}\right].$
(19)

As usual, H.c. stands for Hermitian conjugation, and to abbreviate the
notation we have defined

$T_{\mu\nu}(f_{I},f_{J}):=(\nabla_{\mu}f_{I})(\nabla_{\nu}f_{J})+(\nabla_{\nu}f_{I})(\nabla_{\mu}f_{J})-g_{\mu\nu}\left[(\nabla_{\alpha}f_{I})(\nabla^{\alpha}f_{J})+{m_{0}}^{2}f_{I}f_{J}\right],$
(20)

such that $T_{\mu\nu}(f_{I},f_{I}^{*})$ is the stress energy-momentum tensor
corresponding to a classical complex scalar field of amplitude $f_{I}(x)$.
With the normal order we have imposed, Eq. (19) provides sensible results. In
particular, for coherent states such as (16) one has
$\langle\hat{a}_{I}\hat{a}_{J}\rangle=\alpha_{I}\alpha_{J}$ and
$\langle\hat{a}_{I}^{\dagger}\hat{a}_{J}\rangle=\alpha_{I}^{*}\alpha_{J}$, and
the expectation value of the stress energy-momentum tensor operator takes the
same form as its classical counterpart with
$\phi_{\textrm{cl}}(x)=\langle\hat{\phi}(x)\rangle=\sum_{I}[\alpha_{I}f_{I}(x)+\alpha_{I}^{*}f_{I}^{*}(x)]$,
that is,

$\langle\alpha_{1},\alpha_{2},\ldots|\hat{T}_{\mu\nu}|\alpha_{1},\alpha_{2},\ldots\rangle=\frac{1}{2}T_{\mu\nu}\left(\phi_{\textrm{cl}},\phi_{\textrm{cl}}\right),$
(21)

where the factor $1/2$ on the right hand side of Eq. (21) is due to the
difference in the definition of the stress energy-momentum tensor of a real
and a complex field [cf. Eqs. (19) and (55)]. For a state with a definite
number of particles of the form (13) we have, however,
$\langle\hat{a}_{I}\hat{a}_{J}\rangle=0$ and
$\langle\hat{a}_{I}^{\dagger}\hat{a}_{J}\rangle=N_{I}\delta_{IJ}$, and the
expectation value of the stress energy-momentum tensor reduces to

$\langle
N_{1},N_{2},\ldots|\hat{T}_{\mu\nu}|N_{1},N_{2},\ldots\rangle=\sum\limits_{I}N_{I}T_{\mu\nu}(f_{I},f_{I}^{*}),$
(22)

which is also finite and equal to the weighted sum over the stress energy-
momentum tensors $T_{\mu\nu}(f_{I},f_{I}^{*})$ associated with each mode
function $f_{I}(x)$. Note that there is no analog of Eq. (22) in the classical
real scalar field theory; this is because the eigenstates (13) of the particle
number operator satisfy $\langle\hat{\phi}(x)\rangle=0$, and in this case
quantum fluctuations
$[\langle\hat{\phi}^{2}(x)\rangle-\langle\hat{\phi}(x)\rangle^{2}]^{1/2}$
source the entire stress energy-momentum tensor (22). Note also the different
purpose that the mode functions $f_{I}(x)$ serve in Eqs. (21) and (22);
whereas in the former expression they are associated with the excitations of a
classical field, in the latter, they represent the wave functions of the
quantum particles, which can be also interpreted as $N$ equal classical
complex independent fields.

###  II.3 Statistical ensembles

Up to now, we have restricted our attention to pure states, corresponding to
rays in Hilbert space. More generally, one may consider a statistical ensemble
described by a density operator $\hat{\rho}$, that is, a self-adjoint
nonnegative operator $\hat{\rho}=\hat{\rho}^{\dagger}\geq 0$ of unit trace
$\mbox{Tr}(\hat{\rho})=1$. For the particular case in which this operator is
diagonal with respect to the eigenvectors (13) of the particle number
operator, $\hat{\rho}$ has the representation

$\hat{\rho}=\sum\limits_{N_{1},N_{2},\ldots}(p_{N_{1}N_{2}\ldots})|N_{1},N_{2},\ldots\rangle\langle
N_{1},N_{2},\ldots|,$ (23)

with the probabilities $0\leq p_{N_{1}N_{2}\ldots}\leq 1$ satisfying
$\sum_{N_{1},N_{2},\ldots}(p_{N_{1}N_{2}\ldots})=1$. Although this does not
describe the most general situation, it is sufficient to describe, e.g.,
equilibrium systems at constant temperature $T$, in which case the
probabilities $p_{N_{1}N_{2}\ldots}$ are subject to the Bose-Einstein thermal
equilibrium distribution. A more detailed analysis of such thermal
configurations lies beyond the scope of this article and will be studied in
future work.

When dealing with mixed states, one needs to replace the expectation value
that appears in the semiclassical Einstein equations (17) with the
statistically averaged stress energy-momentum tensor
$\langle\hat{T}_{\mu\nu}\rangle_{\textrm{stat}}=\mbox{Tr}(\hat{\rho}\hat{T}_{\mu\nu})$,
with Tr denoting the trace. For a statistical ensemble of the form (23), it
follows that

$\mbox{Tr}(\hat{\rho}\hat{T}_{\mu\nu})=\sum\limits_{I}\langle
N_{I}\rangle_{\textrm{stat}}T_{\mu\nu}(f_{I},f_{I}^{*}).$ (24)

This has the same form as the right-hand side of Eq. (22), with $N_{I}$
replaced with its statistical average

$\langle
N_{I}\rangle_{\textrm{stat}}:=\sum\limits_{N_{1},N_{2},\ldots}(p_{N_{1}N_{2}\ldots})N_{I},$
(25)

and the previous result (22) is recovered by choosing all the
$p_{N_{1}N_{2}\ldots}$’s equal to zero except for
$p_{0,\ldots,N_{I},0,\ldots}=1$.

###  II.4 Semiclassical self-consistent configurations

In order to address a problem in semiclassical gravity it is convenient to
introduce the notion of semiclassical self-consistent configurations
DiezTejedor2012  (see also references Juarez-Aubry:2022qdp ; Canate:2018wtx ).
A semiclassical self-consistent configuration
$\\{\mathcal{M},ds^{2};\hat{\phi}(x),\hat{\pi}(x),\mathscr{H};|\psi\rangle\in\mathscr{H}\\}$
consists of: a) a spacetime manifold $\mathcal{M}$ equipped with a metric
$ds^{2}$, b) a quantum field theory $\hat{\phi}(x)$, $\hat{\pi}(x)$ with the
Hilbert space $\mathscr{H}$ defined on this fixed classical background
geometry, and c) a state $|\psi\rangle$ in $\mathscr{H}$ such that the Klein-
Gordon equation (2) and the semiclassical Einstein equations (17) are
satisfied simultaneously at every point in the spacetime. This is a non-
trivial task; in order to construct the Hilbert space $\mathscr{H}$, we need
to determine the subspace $X_{+}$ of positive norm solutions $f_{I}(x)$ of the
Klein-Gordon equation (2), and these solutions depend on the spacetime
background which is obtained by solving the semiclassical Einstein equations
(17), so that both the metric field and the quantum state need to be
determined in a self-consistent way. In the case when the quantum theory
consists of a real scalar field, we will say that a semiclassical self-
consistent configuration constitutes a solution to the semiclassical real EKG
theory (2) and (17).

Having said this, we have identified three scenarios for which the expectation
value $\langle\hat{T}_{\mu\nu}\rangle$ of the stress energy-momentum tensor
operator has a special structure: i) coherent states (16), for which
$\langle\hat{T}_{\mu\nu}\rangle$ is equal to the stress energy-momentum tensor
of the corresponding classical solution
$\langle\hat{\phi}(x)\rangle=\sum_{I}[\alpha_{I}f_{I}(x)+\alpha_{I}^{*}f_{I}^{*}(x)]$,
saturating the quantum uncertainty principle; ii) states with a definite
number of particles (13), for which $\langle\hat{T}_{\mu\nu}\rangle$ is
sourced by quantum fluctuations and represents a weighted sum over the
classical stress energy-momentum tensors associated with the complex fields
$f_{I}(x)$; and iii) statistical ensembles described by a density operator
$\hat{\rho}$ of the form (23), for which the statistical average
$\mbox{Tr}(\hat{\rho}\hat{T}_{\mu\nu})$ yields again a weighted sum of
classical stress energy-momentum tensors. In the first case, the semiclassical
system is identical to the classical EKG system for the single, real, free,
minimally coupled scalar field
$\sum_{I}[\alpha_{I}f_{I}(x)+\alpha_{I}^{*}f_{I}^{*}(x)]$. In the second and
third cases the semiclassical equations are equivalent to the classical EKG
system for a family of non-interacting, complex, free, minimally coupled
scalar fields $f_{I}(x)$ which need to form an orthonormal set of basis
functions of the subspace $X_{+}$ of positive norm solutions of the Klein-
Gordon equation. Note the different roles that the mode functions $f_{I}(x)$
play: in scenario i, they combine into a single real field, whereas in
scenarios ii and iii, they all constitute independent complex fields. In this
paper we concentrate mainly on scenarios ii and iii; however, we also discuss
scenario i for the case of a complex scalar field in appendix B.

Up to this point in the presentation, we have intended to provide the reader
with a general perspective of the problem of self-gravitating boson systems in
the semiclassical theory. In the remainder of this article, we focus on static
configurations, in which case there is a well-defined way of performing the
decomposition (9).

##  III Static case

In addition to being globally hyperbolic, we now assume the spacetime $({\cal
M},ds^{2})$ to be static, which implies that there exists a preferred
foliation ${\cal M}=\mathbb{R}\times\Sigma$ of the spacetime manifold such
that the metric has the form

$ds^{2}=-\alpha^{2}(\vec{x})dt^{2}+\gamma_{ij}(\vec{x})dx^{i}dx^{j},$ (26)

i.e., the shift vector is zero, $\beta^{i}=0$, and the lapse function
$\alpha>0$ and the induced three-metric $\gamma_{ij}$ only depend on the
spatial coordinates $\vec{x}$ on $\Sigma$. Introducing this ansatz into the
Klein-Gordon equation (2), we obtain

$\partial_{t}^{2}\phi-\alpha D^{i}\left(\alpha
D_{i}\phi\right)+\alpha^{2}{m_{0}}^{2}\phi=0,$ (27)

where $\partial_{t}:=\partial/\partial t$ is the partial derivative with
respect to the time coordinate, and $D_{i}$ denotes the covariant derivative
operator associated with the induced metric $\gamma_{ij}$. This equation
contains no crossed terms of the form $\partial_{t}D_{i}$ and suggests the
following ansatz for the basis functions,

$f_{I}(t,\vec{x})=\frac{1}{\sqrt{2\omega_{I}}}e^{-i\omega_{I}t}u_{I}(\vec{x}),$
(28)

with $\omega_{I}>0$ and $u_{I}$ a complex-valued function666Due to the fact
that the operator $H$ is real, we could in fact assume that the functions
$u_{I}$ are real-valued. However, for later convenience (see Sec. IV), we
shall only assume that the complex conjugate $u_{I}^{*}$ of $u_{I}$ is
proportional to another member of the same basis, which we call
$u_{I^{\prime}}$. Note that $\omega_{I}=\omega_{I^{\prime}}$. of the spatial
coordinates $\vec{x}$ only, and where the factor $1/\sqrt{2\omega_{I}}$ has
been introduced for future convenience. In this case, the Klein-Gordon
equation (27) leads to the following eigenvalue problem for the square of the
frequency ${\omega_{I}}^{2}$:

$Hu_{I}:=-\alpha D^{j}\left(\alpha
D_{j}u_{I}\right)+\alpha^{2}{m_{0}}^{2}u_{I}={\omega_{I}}^{2}u_{I}.$ (29)

The linear operator $H$ is formally self-adjoint on the Hilbert space $Y$ of
square-integrable functions $u:\Sigma\to\mathbb{C}$, with scalar product

$\langle
u_{1},u_{2}\rangle:=\int_{\Sigma}u_{1}^{*}(\vec{x})u_{2}(\vec{x})\frac{d\gamma}{\alpha(\vec{x})},\quad
u_{1},u_{2}\in Y.$ (30)

Indeed, one can check that, for a suitable definition of the domain $D(H)$ of
the operator incorporating appropriate regularity and fall-off conditions, we
have $\langle u_{1},Hu_{2}\rangle=\langle Hu_{1},u_{2}\rangle$ for all
$u_{1},u_{2}\in D(H)$. Furthermore,

$\langle
u,Hu\rangle=\int_{\Sigma}\left(|Du(\vec{x})|^{2}+{m_{0}}^{2}|u(\vec{x})|^{2}\right)\alpha(\vec{x})d\gamma,$
(31)

which is strictly positive for all $u\in D(H)$ different from zero. Hence, $H$
is a symmetric positive operator, and since it commutes with complex
conjugation, Neumann’s theorem (see Theorem X.3 in ReedSimonVol2 ) implies
that it possesses a positive self-adjoint extension. This offers the
possibility of studying the eigenvalue problem (29) using the powerful tools
of spectral theory for self-adjoint operators ReedSimonVol1 ; ReedSimonVol4 ;
Much2021 .

In the following, we shall assume $H$ has a discrete spectrum with
corresponding eigenvalues $0<\omega_{1}^{2}\leq\omega_{2}^{2}\leq\ldots$, and
associated eigenfunctions $u_{1},u_{2},\ldots$, which can be chosen such that

$\langle u_{I},u_{J}\rangle=\delta_{IJ},\quad I,J=1,2,\ldots.$ (32)

Each of these eigenfunctions gives rise to a (complex-valued) solution of the
Klein-Gordon equation of the form (28), which together with its complex
conjugate solution $f^{*}_{I}(t,\vec{x})$ can easily be verified to satisfy
the following properties:

$(f_{I},f_{J})=-(f_{I}^{*},f_{J}^{*})=\delta_{IJ},\qquad(f_{I},f_{J}^{*})=0,$
(33)

where $(\cdot,\cdot)$ denotes the inner product defined in Eq. (4).

If $H$ has a pure discrete spectrum, then the functions $u_{I}$ provide an
orthonormal basis for $Y$, and the functions $f_{I}$ and $f_{I}^{*}$ defined
by Eq. (28) provide a basis of complex-valued classical solutions of the
Klein-Gordon equation, spanning the spaces of “positive-frequency” and
“negative-frequency” solutions, respectively, which give rise to the spaces
$X_{+}$ and $X_{+}^{*}$ of the decomposition (9). Notice that this choice of
the decomposition makes essential use of the staticity of the spacetime, i.e.
the existence of the globally-defined hypersurface-orthogonal timelike Killing
vector field $\zeta=\partial_{t}$, where the functions $f_{I}$ spanning the
space $X_{+}$ have the property of being eigenstates
$\pounds_{\zeta}f_{I}=-i\omega_{I}f_{I}$ of $\zeta=\partial_{t}$ with
eigenvalues $-i\omega_{I}$, $\omega_{I}>0$ (see the “energy-requirement” of
Ref. Ashtekar:1975zn  for details). If $H$ has a discrete and continuous
spectrum (as will be the case for the boson star solutions discussed in the
next section), the eigenfunctions $f_{I}$ are incomplete; however, they may be
completed by considering “generalized” eigenfunctions lying outside the
Hilbert space $Y$, as it is usually done when dealing, for example, with free
particles in Minkowski space. Alternatively, one can also consider “cutting
off” the spatial domain $\Sigma$ by replacing it with a compact subdomain
$\Sigma_{R}\subset\Sigma$ with a smooth outer boundary $\partial\Sigma_{R}$
with large areal radius $R$, and by solving the eigenvalue problem (29) on
$\Sigma_{R}$ with homogeneous Dirichlet conditions for $u_{I}$ on
$\partial\Sigma_{R}$. One then obtains a pure discrete spectrum at the cost of
introducing the cut-off parameter $R$. However, as long as $R$ is much larger
than the size of the configuration, one would expect boundary effects to be
negligible. Coming back to the case of free particles in Minkowski space, this
is what one usually does when introducing a fictitious box of periodic
boundary conditions and taking the limit of infinite volume at the end of the
calculation. For the configurations that we construct in this paper, only the
discrete spectrum will be excited.

For a static configuration as described by Eqs. (26) and (28), the projections
of the semiclassical Einstein equations (17) normal and tangential to the
hypersurface $\Sigma$ reduce to the system [cf. Eqs. (2.4.10) and (2.5.4) in
Ref. Alcubierre2008 ]

$\displaystyle R^{(3)}$ $\displaystyle=$ $\displaystyle 16\pi G\rho,$ (34a)
$\displaystyle R^{(3)}_{ij}-\frac{1}{\alpha}D_{i}D_{j}\alpha$ $\displaystyle=$
$\displaystyle 4\pi G\left[\gamma_{ij}(\rho-S)+2S_{ij}\right],$ (34b)

where $R^{(3)}_{ij}$ and $R^{(3)}:=\gamma^{ij}R^{(3)}_{ij}$ refer to the
three-dimensional Ricci tensor and Ricci scalar with respect to $\gamma_{ij}$,
and

$\displaystyle\rho$ $\displaystyle:=$ $\displaystyle
n^{\mu}n^{\nu}\langle\hat{T}_{\mu\nu}\rangle,$ (35a) $\displaystyle S_{ij}$
$\displaystyle:=$
$\displaystyle(\delta_{i}{}^{\mu}+n_{i}n^{\mu})(\delta_{j}{}^{\nu}+n_{j}n^{\nu})\langle\hat{T}_{\mu\nu}\rangle$
(35b) are the expectation value of the energy density and the spatial stress
tensor as measured by the so-called Eulerian observers (those moving along the
normal direction to the spatial hypersurfaces), with $S:=\gamma^{ij}S_{ij}$.
For self-consistency with the staticity property, $\rho$ and $S_{ij}$ also
need to be time-independent, and the momentum flux given by
$j_{i}:=(\delta_{i}{}^{\mu}+n_{i}n^{\mu})n^{\nu}\langle\hat{T}_{\mu\nu}\rangle,$
(35c)

must vanish. For an arbitrary state in $\mathscr{H}$, the energy density, the
momentum flux, and the spatial stress tensor can be expressed in the form (51)
of appendix A. As shown in this appendix, $\rho$, $j_{i}$ and $S_{ij}$ are
time-independent as long as $\langle\hat{a}_{I}\hat{a}_{J}\rangle=0$ for all
$I,J$, and $\langle\hat{a}^{\dagger}_{I}\hat{a}_{J}\rangle=0$ whenever
$\omega_{I}\neq\omega_{J}$. These conditions cannot be fulfilled for a non-
trivial coherent state like in Eq. (16), where
$\langle\hat{a}_{I}\hat{a}_{J}\rangle=\alpha_{I}\alpha_{J}$ is different from
zero at least for some values of $I$ and $J$. However, for a state with a
definite number of particles as in Eq. (13), it follows that
$\langle\hat{a}_{I}\hat{a}_{J}\rangle=0$ and
$\langle\hat{a}^{\dagger}_{I}\hat{a}_{J}\rangle=N_{I}\delta_{IJ}$, and we
obtain

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\sum_{I}\frac{N_{I}}{2\omega_{I}}\left[|Du_{I}|^{2}+\left(\frac{{\omega_{I}}^{2}}{\alpha^{2}}+{m_{0}}^{2}\right)|u_{I}|^{2}\right],$
(36a) $\displaystyle j_{k}$ $\displaystyle=$
$\displaystyle\sum_{I}\frac{N_{I}}{2}\frac{i}{\alpha}\left[(D_{k}u_{I})u_{I}^{*}-u_{I}(D_{k}u_{I}^{*})\right],$
(36b) $\displaystyle S_{ij}$ $\displaystyle=$
$\displaystyle\sum_{I}\frac{N_{I}}{2\omega_{I}}\left\\{(D_{i}u_{I})(D_{j}u_{I}^{*})+(D_{j}u_{I})(D_{i}u_{I}^{*})-\gamma_{ij}\left[|D_{i}u_{I}|^{2}-\left(\frac{{\omega_{I}}^{2}}{\alpha^{2}}-{m_{0}}^{2}\right)|u_{I}|^{2}\right]\right\\},$
(36c)

where we have abbreviated
$|Du_{I}|^{2}:=\gamma^{ij}(D_{i}u_{J})(D_{j}u_{J}^{*})$, and where we recall
that the functions $u_{I}$ are subject to the orthogonality condition $\langle
u_{I},u_{J}\rangle=\delta_{IJ}$. If all the $u_{I}$’s are chosen to be real-
valued, the momentum flux obviously vanishes. More generally, if $u_{I}$ is
complex-valued, following the convention in footnote 6, $j_{k}=0$ follows
provided that $N_{I}=N_{I^{\prime}}$. See appendix A for further information
regarding these conditions and their necessity in the context of static and
stationary states.

Taking into account Eqs. (36), and imposing $j_{k}=0$, the Hamiltonian
constraint (34a) and the trace and traceless parts of Eq. (34b) yield

$\displaystyle R^{(3)}$ $\displaystyle=$ $\displaystyle 8\pi
G\sum\limits_{I}\frac{N_{I}}{\omega_{I}}\left[|Du_{I}|^{2}+\left(\frac{{\omega_{I}}^{2}}{\alpha^{2}}+{m_{0}}^{2}\right)|u_{I}|^{2}\right],$
(37a) $\displaystyle\frac{D^{j}D_{j}\alpha}{\alpha}$ $\displaystyle=$
$\displaystyle 8\pi
G\sum\limits_{I}\frac{N_{I}}{\omega_{I}}\left[\left(\frac{{\omega_{I}}^{2}}{\alpha^{2}}-\frac{{m_{0}}^{2}}{2}\right)|u_{I}|^{2}\right],$
(37b)
$\displaystyle\left[R^{(3)}_{ij}-\frac{1}{\alpha}D_{i}D_{j}\alpha\right]^{\rm
tf}$ $\displaystyle=$ $\displaystyle 4\pi
G\sum\limits_{I}\frac{N_{I}}{\omega_{I}}\left[(D_{i}u_{I})(D_{j}u_{I}^{*})+(D_{j}u_{I})(D_{i}u_{I}^{*})\right]^{\rm
tf},$ (37c)

where the superscript “tf” refers to the trace-free part with respect to
$\gamma_{ij}$, i.e. $(A_{ij})^{\rm
tf}:=A_{ij}-\frac{1}{3}\gamma_{ij}(\gamma^{mn}A_{mn})$. These equations,
together with the Klein-Gordon equation (29), constitute a nonlinear multi-
eigenvalue problem for the frequencies $\omega_{I}$ describing a system of
$N=\sum_{I}N_{I}$ identical quantum particles in self-gravitating equilibrium.
Note that Eqs. (37) are also applicable to systems that are described in terms
of a statistical ensemble of the form (23); in this case, $N_{I}$ needs to be
replaced with its statistical average $\langle N_{I}\rangle_{\textrm{stat}}$,
see Eq. (25).

In the next section, we further specialize these equations to the static,
spherically symmetric case, and we show that for states of definite number of
particles, as well as for statistical ensembles of the form (23), the
resulting equations give rise to the $\ell$-boson star configurations
constructed in Alcubierre:2018ahf  (see also Alcubierre:2019qnh ;
Alcubierre:2021mvs ), and even to more general solutions, a few examples of
which are constructed numerically in section V.

##  IV Static spherically symmetric configurations

We now further specialize to a static, spherically symmetric spacetime, for
which the three-metric can be expressed in the form:

$\gamma_{ij}dx^{i}dx^{j}=\gamma^{2}dr^{2}+r^{2}d\Omega^{2},\qquad\gamma=\left(1-\frac{2GM}{r}\right)^{-1/2},$
(38)

where $M$ denotes the Misner-Sharp mass function and
$d\Omega^{2}=d\vartheta^{2}+\sin^{2}\vartheta d\varphi^{2}$ is the standard
line-element on the unit two-sphere $S^{2}$. Furthermore, in these
coordinates, the lapse $\alpha$, the function $\gamma$ and the Misner-Sharp
mass $M$ only depend on the areal radius coordinate $r$. Because of the
spherical symmetry, the mode solutions of the Klein-Gordon equation (29) can
be assumed to be of the form

$u_{I}(\vec{x})=v_{n\ell}(r)Y^{\ell m}(\vartheta,\varphi),\quad I=(n\ell m),$
(39)

with $Y^{\ell m}$ denoting the standard spherical harmonics, and where no sum
in the total angular momentum number $\ell$ is considered. Note that
$u_{I}(\vec{x})^{*}=(-1)^{m}u_{I^{\prime}}(\vec{x})$ with
$I^{\prime}=(n,\ell,-m)$, such that the property assumed in footnote 6 is
satisfied. Since the magnetic number $m$ does not appear explicitly in the
radial differential equation

$-\frac{\alpha}{\gamma r^{2}}\left(\frac{\alpha
r^{2}}{\gamma}v^{\prime}_{n\ell}\right)^{\prime}+\alpha^{2}\left[\frac{\ell(\ell+1)}{r^{2}}+{m_{0}}^{2}\right]v_{n\ell}={(\omega_{n\ell})}^{2}v_{n\ell}$
(40)

that is obtained from Eq. (29) with the ansatz (39), the radial functions
$v_{n\ell}(r)$ can be chosen to be independent of this number. Therefore, the
eigenvalue problem (29) reduces to finding (for each $\ell=0,1,2,\ldots$) a
set of suitable radial basis functions $v_{n\ell}$ solving Eq. (40). Using the
orthonormality property of the spherical harmonics, $\int_{S^{2}}Y^{\ell
m}Y^{\ell^{\prime}m^{\prime}*}d\Omega=\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}$,
the orthogonality condition $\langle u_{I},u_{J}\rangle=\delta_{IJ}$ reduces
to

$\int_{0}^{\infty}v_{n\ell}(r)v^{*}_{n^{\prime}\ell}(r)\frac{\gamma(r)}{\alpha(r)}\>r^{2}dr=\delta_{nn^{\prime}}.$
(41)

Assuming the functions $\alpha$ and $\gamma$ are regular at the center $r=0$,
such that they have local expansions of the form
$\alpha(r)=\alpha_{0}+\alpha_{2}r^{2}+\ldots$ and
$\gamma(r)=1+\gamma_{2}r^{2}+\ldots$, one can show that the local solution
that is finite at $r=0$ has the form $v_{n\ell}(r)\sim r^{\ell}$; see Ref.
Alcubierre:2018ahf . Likewise, as $r\to\infty$, we impose that the metric
functions $\alpha$ and $\gamma$ converge to 1, and that $v_{n\ell}(r)$ are
bounded, which implies that they have the form $v_{n\ell}(r)\sim
e^{-\sqrt{{m_{0}}^{2}-{(\omega_{n\ell})}^{2}}r}$, with
$0<\omega_{n\ell}<m_{0}$. A further restriction arises from the identity
$\langle u_{I},Hu_{I}\rangle={\omega_{I}}^{2}$, which yields

$\int_{\Sigma}\left\\{|Du_{I}|^{2}+\left[{m_{0}}^{2}-\frac{{\omega_{I}}^{2}}{\alpha^{2}}\right]|u_{I}|^{2}\right\\}\alpha
d\gamma=0,$ (42)

and shows that ${m_{0}}^{2}-{\omega_{I}}^{2}/\alpha^{2}$ cannot be positive
everywhere, since otherwise it would follow from Eq. (42) that $u_{I}=0$.

The functions $\alpha$ and $\gamma$ must be determined by solving the static
semiclassical Einstein field equations (37), where for consistency the right-
hand side must be a spherically symmetric tensor.777A different approach to
achieving this property was considered in Mielke:1980sa , where the stress
energy-momentum tensor is averaged over the spheres in order to get rid of the
angular dependency. Our approach requires no such averaging. This is clearly
the case if only the ground state is populated, i.e., if $N_{I}=0$ for all
$I\neq(000)$, which gives rise to the standard boson star equations. More
generally, one can demand that $N_{n\ell m}=0$ for all $\ell>0$, meaning that
all the particles have zero angular momentum but may nevertheless be in
excited energy states. This gives rise to the multi-state boson star equations
solved in Matos2007 ; Bernal:2009zy ; UrenaLopez:2010ur . Following the same
arguments as in appendix A of Ref. Alcubierre:2018ahf , in order for the
expectation value of the stress energy-momentum tensor to be spherically
symmetric, it is in fact sufficient to choose $N_{n\ell m}$ independent of the
magnetic number $m$, such that

$N_{n,\ell,-\ell}=N_{n,\ell,-(\ell-1)}=\ldots=N_{n,\ell,(\ell-1)}=N_{n,\ell,\ell},$
(43)

which implies that the total angular momentum vanishes, even if the individual
particles posses angular momentum. Note that this choice also guarantees that
the momentum flux is zero, as required for staticity. In other words, the
excitation numbers $N_{n\ell m}$ are functions of the energy levels $n$ and
the total angular momentum $\ell$, but not of the magnetic quantum number $m$.
This is rather similar to the case of a kinetic gas, in which a one-particle
distribution function depending only on the energy and the total angular
momentum gives rise to a static, spherically symmetric configuration (see, for
instance, section 5.1 in Andreasson-LivRev ).

Assuming the validity of condition (43) and the spherically symmetric ansatz
(38), equations (37) reduce to

$\displaystyle\frac{2GM^{\prime}}{r^{2}}$ $\displaystyle=$
$\displaystyle\sum\limits_{n\ell}\frac{\kappa_{\ell}N_{n\ell
m}}{\omega_{n\ell}}\left[\frac{|v_{n\ell}^{\prime}|^{2}}{\gamma^{2}}+\left(\frac{{(\omega_{n\ell})}^{2}}{\alpha^{2}}+{m_{0}}^{2}+\frac{\ell(\ell+1)}{r^{2}}\right)|v_{n\ell}|^{2}\right],$
(44a) $\displaystyle\frac{1}{\alpha\gamma
r^{2}}\left(\frac{r^{2}\alpha^{\prime}}{\gamma}\right)^{\prime}$
$\displaystyle=$ $\displaystyle\sum\limits_{n\ell}\frac{\kappa_{\ell}N_{n\ell
m}}{\omega_{n\ell}}\left[\left(2\frac{{(\omega_{n\ell})}^{2}}{\alpha^{2}}-{m_{0}}^{2}\right)|v_{n\ell}|^{2}\right],$
(44b) $\displaystyle\frac{(\alpha\gamma)^{\prime}}{r\alpha\gamma^{3}}$
$\displaystyle=$ $\displaystyle\sum\limits_{n\ell}\frac{\kappa_{\ell}N_{n\ell
m}}{\omega_{n\ell}}\left[\frac{|v_{n\ell}^{\prime}|^{2}}{\gamma^{2}}+\frac{{(\omega_{n\ell})}^{2}}{\alpha^{2}}|v_{n\ell}|^{2}\right],$
(44c)

with $\kappa_{\ell}:=(2\ell+1)G$, and where the last identity was obtained by
contracting the angular components of Eq. (37c) with the metric of the unit
two-sphere, $\hat{g}_{AB}$, and making use of Eq. (44a). In addition, we have
also used the identities $\sum_{m=-\ell}^{\ell}Y^{\ell m}Y^{\ell
m*}=\frac{1}{4\pi}(2\ell+1)$ and
$\sum_{m=-\ell}^{\ell}(\hat{\nabla}_{A}Y^{\ell m})(\hat{\nabla}^{A}Y^{\ell
m*})=\frac{1}{4\pi}\ell(\ell+1)(2\ell+1)$, where $\hat{\nabla}_{A}$ makes
reference to the covariant derivative with respect to $\hat{g}_{AB}$ (see
appendix A in Alcubierre:2018ahf  for details). The full system of reduced
static, spherically symmetric, semiclassical real EKG equations consists of
Eqs. (40) and (44), where, due to the twice contracted Bianchi identities, Eq.
(44b) can be omitted. Further, the eigenfunctions $v_{n\ell}$ should satisfy
the normalization condition (41), although we can also do without this
equation if we absorb the occupation numbers $N_{n\ell m}$ in the radial
functions $v_{n\ell}$, as described in the next section. Once the functions
$v_{n\ell}$ are known, the quantum field $\hat{\phi}(x)$ can be reconstructed
using Eqs. (12), (28), and (39), which gives

$\hat{\phi}(x)=\sum\limits_{n\ell
m}\frac{1}{\sqrt{2\omega_{n\ell}}}\left[\hat{a}_{I}e^{-i\omega_{n\ell}t}v_{n\ell}(r)Y^{\ell
m}(\vartheta,\varphi)+\textrm{H.c.}\right].$ (45)

The $\ell$-boson star configurations we have discussed in Refs.
Alcubierre:2018ahf ; Alcubierre:2019qnh ; Alcubierre:2021mvs  are obtained by
solving a particular case of this system, in which all the $N_{n\ell m}$’s
vanish except the ones for $n=0$ and some specific value of $\ell$. In this
case, after absorbing the factor $N_{0\ell m}/\omega_{0\ell}$ into the
amplitude of $v_{0\ell}$, the system of Eqs. (44a), (44c), and (40) reduces
precisely to the system (7a), (7b), and (7c) of Alcubierre:2018ahf . However,
in contrast to the purely classical description in Alcubierre:2018ahf ;
Alcubierre:2019qnh ; Alcubierre:2021mvs , which requires precisely $2\ell+1$
complex scalar fields, the semiclassical interpretation of the $\ell$-boson
stars becomes much more natural: they correspond to a particular excitation of
a single real quantum spin zero field that describes a selfgravitating system
of $(2\ell+1)N_{0\ell m}$ identical quantum particles of definite energy
$E=\omega_{0\ell}$ and angular momentum $L=\sqrt{\ell(\ell+1)}$ (both
evaluated in natural units). Furthermore, as is evident from the equations
above, there are many other possible configurations that can be constructed in
this way, involving excitations of different energy levels $n$ and different
total angular momentum numbers $\ell$. We summarize these more general
solutions and their sub-families, as well as a few references to corresponding
Newtonian configurations, in Table 2. Numerical examples of some of these more
general configurations are constructed in the next section.

Table 2:  Classification of the solutions to the static, spherically
symmetric, semiclassical real EKG system. They represent self-gravitating
equilibrium configurations of a definite number of identical quantum
particles. These are all the cases obtained when combining two options for the
radial quantum number $n$ and three options for the total angular momentum
number $\ell$. The two options for $n$ are: i) multiple values and ii) one
value; while the three options for $\ell$ are: i) multiple values, ii) one
value, and iii) fixed value $\ell=0$. Note that some of these solutions are
particular cases of others. In order to show this hierarchy more clearly, we
enclose with a bracket solutions that are included in a more general solution,
which is indicated with an arrow.

|  |  |  |  |  | Name | $n$ | $\ell$ | Relativistic | Newtonian
---|---|---|---|---|---|---|---|---|---|---
|  |  |  |  |  | Multi-$\ell$ multi-state boson star | $n_{1},n_{2},\dotsc,n_{p}$ | $\ell_{1},\ell_{2},\dotsc,\ell_{q}$ | Sec. V | $\ldots$
|  |  |  |  |  | Multi-state $\ell$-boson star | $n_{1},n_{2},\dotsc,n_{p}$ | $\ell_{1}$ | Sec. V | $\ldots$
|  |  |  |  |  | Multi-state boson star | $n_{1},n_{2},\dotsc,n_{p}$ | $0$ | Matos2007 ; Bernal:2009zy ; UrenaLopez:2010ur  | KavianMuschler2015
|  |  |  |  |  | Boson star | $n_{1}$ | $0$ | Kaup68 ; Ruffini:1969qy ; Colpi:1986ye ; Friedberg87 ; Gleiser:1988rq ; Lee:1988av ; Guzman:2004wj  | Lieb77 ; MorozPenroseTod1998
|  |  |  |  |  |  $\ell$-Boson star | $n_{1}$ | $\ell_{1}$ | Alcubierre:2018ahf ; Alcubierre:2019qnh ; Alcubierre:2021mvs ; Alcubierre:2021psa  | jaramillo19 ; nambo21 ; GU2020
|  |  |  |  |  | Multi-$\ell$ boson star | $n_{1}$ | $\ell_{1},\ell_{2},\dotsc,\ell_{q}$ | Sec. V | GU2020

##  V Numerical solutions: a few examples

In this section, we present numerical solutions to the static, spherically
symmetric, semiclassical real EKG system described by Eqs. (40), (44a), and
(44c). These solutions complement the mathematical analysis of the previous
section. Specifically, we obtain three particular solutions as representative
examples of the three types of solutions that have not been presented so far
in the literature (see Table 2), that is, a multi-$\ell$ boson star, a multi-
state $\ell$-boson star, and a multi-$\ell$ multi-state boson star.

To proceed, from this point onward, in addition to $\hbar=c=1$, we also set
$G=1$, such that all quantities are dimensionless and measured in Planck
units, although the solutions can be rescaled arbitrarily using a symmetry
transformation, as we explain later. In practice, we solve the semiclassical
real EKG system (40), (44a), and (44c) expressed in the form

$\displaystyle\psi_{n\ell}^{\prime\prime}$ $\displaystyle=$
$\displaystyle-\left[\gamma^{2}+1-(2\ell+1)r^{2}\gamma^{2}\left(\frac{\ell(\ell+1)}{r^{2}}+{m_{0}}^{2}\right){(\psi_{n\ell})}^{2}\right]\frac{\psi_{n\ell}^{\prime}}{r}-\left(\frac{{(\omega_{n\ell})}^{2}}{\alpha^{2}}-\frac{\ell(\ell+1)}{r^{2}}-{m_{0}}^{2}\right)\gamma^{2}\psi_{n\ell},$
(46a) $\displaystyle\gamma^{\prime}$ $\displaystyle=$
$\displaystyle\sum_{n\ell}\frac{2\ell+1}{2}r\gamma\left[\left(\frac{{(\omega_{n\ell})}^{2}}{\alpha^{2}}+\frac{\ell(\ell+1)}{r^{2}}+{m_{0}}^{2}\right)\gamma^{2}{(\psi_{n\ell})}^{2}+{(\psi_{n\ell}^{\prime})}^{2}\right]-\left(\frac{\gamma^{2}-1}{2r}\right)\gamma,$
(46b) $\displaystyle\alpha^{\prime}$ $\displaystyle=$
$\displaystyle\sum_{n\ell}\frac{2\ell+1}{2}r\alpha\left[\left(\frac{{(\omega_{n\ell})}^{2}}{\alpha^{2}}-\frac{\ell(\ell+1)}{r^{2}}-{m_{0}}^{2}\right)\gamma^{2}{(\psi_{n\ell})}^{2}+{(\psi_{n\ell}^{\prime})}^{2}\right]+\left(\frac{\gamma^{2}-1}{2r}\right)\alpha,$
(46c)

where for convenience we have introduced the rescaled fields

$\psi_{n\ell}=\sqrt{\frac{N_{n\ell m}}{\omega_{n\ell}}}v_{n\ell}.$ (47)

Finally, one can choose an arbitrary value for the mass $m_{0}$, since
solutions for a different value can then be obtained by a simple rescaling;
see Eq. (50) below. In particular, we set $m_{0}=1$ for the numerical
integrations, but present the results in an $m_{0}$-independent form. Notice
that from the normalization condition in Eq. (41) one can read off the number
of particles in the different states using expression

$N_{n\ell
m}=\omega_{n\ell}\int_{0}^{\infty}{(\psi_{n\ell})}^{2}\frac{\gamma}{\alpha}r^{2}dr.$
(48)

The choice of appropriate boundary conditions must guarantee that the boson
star solutions are regular and asymptotically flat, and additionally that they
possess finite total energy and finite energy density everywhere. Demanding
regularity at the origin, i.e.,

$\displaystyle\psi_{n\ell}(r)$ $\displaystyle=$
$\displaystyle\frac{\psi_{n\ell}^{0}}{2\ell+1}\>r^{\ell},$ (49a)
$\displaystyle\psi_{n\ell}^{\prime}(r)$ $\displaystyle=$
$\displaystyle\frac{\ell\psi_{n\ell}^{0}}{2\ell+1}\>r^{\ell-1},$ (49b)
$\displaystyle\alpha(r)$ $\displaystyle=$ $\displaystyle 1,$ (49c)
$\displaystyle\gamma(r)$ $\displaystyle=$ $\displaystyle 1,$ (49d)

when $r\to 0$, and a vanishing field at infinity, one obtains a nonlinear
multiple-eigenvalue problem for the different mode-frequencies
$\omega_{n\ell}$. Here, $\psi_{n\ell}^{0}$ are some arbitrary positive
constants related to the number of particles in the different energy levels,
and with no loss of generality, we have fixed the value of the lapse function
at the origin to 1. Notice that, since the system of equations is invariant
under $(\alpha,\omega_{n\ell})\mapsto\lambda(\alpha,\omega_{n\ell})$, with
some positive arbitrary constant $\lambda$, one can always rescale the value
of the lapse function in such a way that $\alpha(r\to\infty)=1$, as we do
later.

The integration of the system is performed numerically using a shooting
algorithm to find the frequencies $\omega_{n\ell}$. To proceed, one integrates
the system of Eqs. (46) outward, starting from the initial conditions in Eqs.
(49) at a point very close to the origin (we used $r_{0}=5\times 10^{-6}$),
and search for the values of the frequencies $\omega_{n\ell}$ to match the
asymptotic behavior of the mode functions until the shooting parameter
converges to the desired accuracy. As we already mentioned, for simplicity we
have assumed $m_{0}=1$, although the solutions can be rescaled to an arbitrary
value of the mass parameter using the invariance of the system under the
transformation

$m_{0}\mapsto\lambda
m_{0},\quad\omega_{n\ell}\mapsto\lambda\omega_{n\ell},\quad
r\mapsto\lambda^{-1}r,$ (50)

with the functions $\psi_{n\ell}$, $\alpha$ and $\gamma$ unchanged. Under this
transformation, the occupation numbers (48) change according to $N_{n\ell
m}\mapsto\lambda^{-2}N_{n\ell m}$. Note that, for instance, $\lambda\sim
10^{-50}$ in the case of an ultralight axion dark matter particle of mass
$m_{0}\sim 10^{-22}\,$eV.

In Figure 1, we present the results of our numerical solutions. In the left
column —subfigures LABEL:sub@figa— we present a multi-$\ell$ boson star
solution. The configuration displayed is characterized by the quantum numbers
$n=0$ and $\ell=0,\,1,\,2$.888Note that, as in our previous works, $n$ refers
to the number of nodes of the radial function in the interval $0<r<\infty$.
Hence, it should be compared to the “radial quantum number” in the theory of
the hydrogen atom rather than the “principal quantum number”. Unlike the
$\ell$-boson stars introduced in Alcubierre:2018ahf , these new solutions do
not require all the individual fields to have the same amplitude, and, still,
the resulting spacetime is spherically symmetric. In the central column
—subfigures LABEL:sub@figb— we present a multi-state $\ell$-boson star
solution. The configuration displayed is characterized by the quantum numbers
$n=0,\,1,\,2$ and $\ell=1$. In the right column —subfigures LABEL:sub@figc— we
present a multi-$\ell$ multi-state boson star solution. The configuration
displayed is characterized by the quantum numbers $n=1,\,3$ and
$\ell=0,\,1,\,2$. In all three cases the lapse function $\alpha$ and the
metric component $\gamma$ of each configuration are shown in the top panels,
the radial profiles for the different wave functions $\psi_{n\ell}$ that are
excited in the configuration are shown in the middle panels, and the total
energy density, as well as the density of the individual energy levels are
shown in the bottom panel. In Table 3, we report the main numbers associated
with these configurations.

(a) Multi-$\ell$ boson star  (b) Multi-state $\ell$-boson star  (c)
Multi-$\ell$ multi-state boson star

Figure 1: As examples we describe here three particular solutions, each one of a different type, as indicated in the subfigures. For all cases we show the metric components, $\alpha(r)$ and $\gamma(r)$ (top panels); the state components $\psi_{n\ell}(r)$ (middle panels); and the energy density $\rho(r)$, as well as the individual components $\rho_{n\ell}(r)$ (bottom panels). Note that $\rho$ can be obtained simply by summing its components (including all degenerate states in $m$). Thus, we have $\rho=\rho_{00}+3\rho_{01}+5\rho_{02}$ in LABEL:sub@figa, $\rho=3\rho_{01}+3\rho_{11}+3\rho_{21}$ in LABEL:sub@figb, and $\rho=\rho_{30}+3\rho_{11}+5\rho_{12}$ in LABEL:sub@figc. We show more properties of these solutions in Table 3. Finally, we note that the actual numerical integration domain in $r$ is much larger than what is shown (for clarity) in the figures. Table 3: Quantum numbers, amplitudes, eigenfrequencies, particle numbers in each state, and total particle number for the configurations presented in Fig. 1. Note that for the case of an ultralight axion $m_{0}\sim 10^{-50}$ in natural units; hence in this example the particle numbers are large for the configurations presented. Name | $n$ | $\ell$ | $\psi^{0}_{n\ell}/{m_{0}}^{\ell}$ | $\omega_{n\ell}/m_{0}$ | ${m_{0}}^{2}N_{nl}$ | ${m_{0}}^{2}N$ | Fig.
---|---|---|---|---|---|---|---
| $0$ | $0$ | $0.1$ | $0.5278$ | $0.0195$ |  |
Multi-$\ell$ boson star | $0$ | $1$ | $0.2$ | $0.6453$ | $0.0243$ | 0.8134 |  1LABEL:sub@figa
| $0$ | $2$ | $0.4$ | $0.7736$ | $0.1442$ |  |
| $0$ | $1$ | $0.3$ | $0.7438$ | $0.0289$ |  |
Multi-state $\ell$-boson star | $1$ | $1$ | $0.6$ | $0.8235$ | $0.1150$ | 1.3884 |  1LABEL:sub@figb
| $2$ | $1$ | $0.9$ | $0.8792$ | $0.3189$ |  |
| $3$ | $0$ | $0.1$ | $0.8679$ | $0.0133$ |  |
Multi-$\ell$ multi-state boson star | $1$ | $1$ | $0.3$ | $0.7497$ | $0.0439$ | 1.2745 |  1LABEL:sub@figc
| $1$ | $2$ | $0.5$ | $0.8247$ | $0.2259$ |  |

In this section we presented just three particular solutions as illustrative
examples of the three new types of boson stars described in this work.
However, going into further analysis of such solutions is beyond the scope of
this article. Instead, we expect to carry out a more detailed study in a
separate work.

##  VI Discussion and conclusions

We have shown that static, spherically symmetric boson star configurations
Kaup68 ; Ruffini:1969qy ; Colpi:1986ye ; Friedberg87 ; Gleiser:1988rq ;
Lee:1988av ; Guzman:2004wj ; Bizon:2000 ; Chavanis:2011cz ; Jetzer:1991jr ;
Schunck:2003kk ; Liebling:2012fv ; Zhang:2018slz ; Visinelli:2021uve  and many
of their generalizations (including $\ell$-boson stars Alcubierre:2018ahf ;
Alcubierre:2019qnh ; Alcubierre:2021mvs ; Alcubierre:2021psa ;
Jaramillo:2022zwg  and multi-state boson stars Matos2007 ; Bernal:2009zy ;
UrenaLopez:2010ur ), arise naturally within the semiclassical gravity
approximation in quantum field theory. Furthermore, we have found new possible
generalizations, namely the multi-$\ell$ multi-state boson stars, that
represent the most general solutions to the $N$-particle, static, spherically
symmetric, semiclassical real EKG system, in which the total number of
particles is definite.

Our approach is based on the expansion of a single, real, free quantum scalar
field in terms of a linear combination of creation and annihilation operators.
We then construct the Hilbert space by successive applications of creation
operators on the vacuum state, that we assume is well defined (which is
guaranteed for the static configurations considered in this article). Taking
particular Fock states with a definite number of particles, the expectation
value of the renormalized stress energy-momentum tensor operator takes the
same form as its counterpart in a classical theory with $N$ fields. Each of
these $N$ fields accounts for one excitation mode of the quantum field, and
corresponds to the quantum particles of the system. The number of particles
contained in each mode is then fixed by the number of classical fields with
the same quantum numbers and can be chosen as a free parameter. In this way,
standard boson stars correspond to the population of only one mode (the ground
state) of the quantum field, for which $n=\ell=m=0$. Other self-gravitating
scalar configurations such as $\ell$-boson stars and multi-state boson stars
are naturally related to them and are just different manifestations of the
quantum fluctuations of the same scalar field with different modes populated.

Regarding $\ell$-boson stars in particular, we would like to highlight an
important difference between the classical and semiclassical interpretations:
the construction of the classical solutions described in our previous works
might be considered somehow “artificial”, in the sense that they must consist
of a particular combination of $2\ell+1$ independent, complex scalar fields
having the exact same radial amplitude in order to constitute as a whole a
spherically symmetric matter distribution. On the other hand, within the
semiclassical gravity approximation, the same $\ell$-boson star configurations
are obtained more naturally, starting from a single real quantum scalar field
in a state that describes a distribution of particles populating a given
energy level and containing angular momentum. Apart from providing a more
natural explanation for their existence, the semiclassical description of
$\ell$-boson stars should have other interesting applications. For example, a
relevant problem is analyzing the dynamical stability of these configurations
in semiclassical gravity. This is clearly a much more difficult problem that
requires generalizing our framework for spacetimes being perturbed off a
static one.

Further, in the realm of an ultralight scalar field/fuzzy dark matter
component, the diversity of self-gravitating structures that might model the
dark matter halos is greatly enhanced in the semiclassical theory. Standard
$\ell=0$ boson stars have been considered as a possible explanation for the
dark matter distribution in small galaxies, such as dwarf spheroidals
Gonzalez-Morales:2016yaf . However, it is clear that these configurations
alone are not sufficient to describe the properties observed in galaxies
spanning different size scales. On the one hand, the radius of stable $\ell=0$
boson stars decreases when their total mass increases, in stark contrast to
what is being observed in real galaxies. On the other hand, their mass density
decreases exponentially with the radius at large distances, and they cannot
accommodate the flatness of the rotation curves observed in spiral galaxies
(although they can still describe the core of galaxies like the Milky Way
DeMartino:2018zkx ). The multi-$\ell$ multi-state configurations offer the
possibility to account for the desired behavior of the dark matter
distribution in large galaxies beyond the core region. Likewise, several
examples of anomalous halo systems recently reported in Atlas3D  could
potentially be explained by means of scalar field configurations with
nontrivial values of $n$, $\ell$, and $m$. The connection between the
multi-$\ell$ multi-state boson star solutions and the configurations obtained
in numerical simulations Schive:2014dra ; Schive:2014hza ; Schwabe:2016rze ;
Veltmaat:2016rxo ; Du:2016aik  is left as an open question.

With regard to dark compact objects, the state of maximal compactness of
$\ell$-boson stars has been identified to grow with $\ell$, as discussed in
Alcubierre:2021psa . This could lead to potential observational signatures in
the gravitational wave spectra resulting from the merger of these objects
Palenzuela:2017kcg .

In this article we have only presented some representative cases by solving
numerically the spherically symmetric, static semiclassical system of
equations for certain quantum numbers; however a more exhaustive analysis is
required in order to compare our solutions with the astrophysical data. In
conclusion, the interpretation of boson stars within the semiclassical gravity
approximation considerably increases the possibility for boson fields to
describe more realistic astrophysical systems.

###### Acknowledgements.

A.D.T. acknowledges conversations with Cristian Armendariz-Picon and Daniel
Sudarsky. O.S. wishes to thank stimulating discussions with Robert Oeckl,
Jorge Pullin, Armando Roque, and Thomas Zannias. This work was supported in
part by the CONACYT Network Projects No. 376127 “Sombras, lentes y ondas
gravitatorias generadas por objetos compactos astrofísicos”, No. 304001
“Estudio de campos escalares con aplicaciones en cosmología y astrofísica”,
No. 1406300 “Explorando los confines de las teorías relativistas de la
gravitación y sus consecuencias”, and No. 286897 “Materia oscura:
Implicaciones de sus propiedades fundamentales en las observaciones
astrofísicas y cosmológicas”, by DGAPA-UNAM grant IN105920, by the European
Union’s Horizon 2020 research and innovation (RISE) program H2020-MSCA-
RISE-2017 Grant No. FunFiCO-777740 and by a CIC grant to Universidad
Michoacana de San Nicolás de Hidalgo.

##  Appendix A Stress energy-momentum tensor

In this appendix we summarize the results for the energy density (35a),
momentum flux (35c) and spatial stress tensor (35b) of an arbitrary quantum
state on a static spacetime that are used in section III. These expressions
are

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\sum_{I,J}\frac{1}{4\sqrt{\omega_{I}\omega_{J}}}\left\\{\langle\hat{a}_{I}\hat{a}_{J}\rangle\left[\left(-\frac{\omega_{I}\omega_{J}}{\alpha^{2}}+{m_{0}}^{2}\right)u_{I}u_{J}+(D_{k}u_{I})(D^{k}u_{J})\right]e^{-i(\omega_{I}+\omega_{J})t}\right.$
(51a) $\displaystyle\hskip
54.06006pt\left.+\langle\hat{a}^{\dagger}_{I}\hat{a}_{J}\rangle\left[\left(+\frac{\omega_{I}\omega_{J}}{\alpha^{2}}+{m_{0}}^{2}\right)u_{I}^{*}u_{J}+(D_{k}u_{I}^{*})(D^{k}u_{J})\right]e^{+i(\omega_{I}-\omega_{J})t}+c.c.\right\\},$
$\displaystyle j_{k}$ $\displaystyle=$
$\displaystyle\sum_{I,J}\frac{1}{4\sqrt{\omega_{I}\omega_{J}}}\frac{i}{\alpha}\left\\{\langle\hat{a}_{I}\hat{a}_{J}\rangle\left[-\omega_{J}(D_{k}u_{I})u_{J}-\omega_{I}u_{I}(D_{k}u_{J})\right]e^{-i(\omega_{I}+\omega_{J})t}\right.$
(51b) $\displaystyle\hskip
64.01869pt\left.+\langle\hat{a}^{\dagger}_{I}\hat{a}_{J}\rangle\left[-\omega_{J}(D_{k}u_{I}^{*})u_{J}+\omega_{I}u_{I}^{*}(D_{k}u_{J})\right]e^{+i(\omega_{I}-\omega_{J})t}-c.c.\right\\},$
and $\displaystyle S_{ij}$ $\displaystyle=$
$\displaystyle\sum_{I,J}\frac{1}{4\sqrt{\omega_{I}\omega_{J}}}$ (51c)
$\displaystyle\times$
$\displaystyle\left\\{\langle\hat{a}_{I}\hat{a}_{J}\rangle\left[(D_{i}u_{I})(D_{j}u_{J})+(D_{j}u_{I})(D_{i}u_{J})-\gamma_{ij}\left(\left(+\frac{\omega_{I}\omega_{J}}{\alpha^{2}}+{m_{0}}^{2}\right)u_{I}u_{J}+(D_{k}u_{I})(D^{k}u_{J})\right)\right]e^{-i(\omega_{I}+\omega_{J})t}\right.$
$\displaystyle\left.+\langle\hat{a}^{\dagger}_{I}\hat{a}_{J}\rangle\left[(D_{i}u_{I}^{*})(D_{j}u_{J})+(D_{j}u_{I}^{*})(D_{i}u_{J})-\gamma_{ij}\left(\left(-\frac{\omega_{I}\omega_{J}}{\alpha^{2}}+{m_{0}}^{2}\right)u_{I}^{*}u_{J}+(D_{k}u_{I}^{*})(D^{k}u_{J})\right)\right]e^{+i(\omega_{I}-\omega_{J})t}\right.$
$\displaystyle\left.+c.c.\right\\},$

where $c.c.$ stands for complex conjugation.

In this article, we are interested in static configurations, that is in states
satisfying the property that $\rho$ and $S_{ij}$ are time-independent and
$j_{k}=0$, such that the expectation value of the stress energy-momentum
tensor is compatible, through the semiclassical Einstein equations, with the
static spacetime metric (26). It is clear from Eqs. (51) that $\rho$, $j_{k}$,
and $S_{ij}$ are time-independent if $\langle\hat{a}_{I}\hat{a}_{J}\rangle=0$
for all $I,J$ and $\langle\hat{a}_{I}^{\dagger}\hat{a}_{J}\rangle=0$ whenever
$\omega_{I}\neq\omega_{J}$. Furthermore, assuming the convention in footnote
6, the additional condition
$\langle\hat{N}_{I}\rangle=\langle\hat{N}_{I^{\prime}}\rangle$ implies that
$j_{k}=0$.

The question of whether these conditions are also necessary for a static
configuration is more subtle. To proceed, we introduce the notion of
stationary states, defined as states in the Schrödinger picture which are
invariant under time translations. In quantum mechanics time translations are
represented by unitary operators of the form
$\hat{\mathcal{U}}[t_{0},t]=\exp[-i\hat{H}(t-t_{0})]$, where the Hamiltonian
$\hat{H}$ is the generator of the transformations. If a state $|\psi\rangle$
is invariant under time evolution, then
$\hat{\mathcal{U}}[t_{0},t]|\psi\rangle=|\psi\rangle$, apart from a possible
phase that does not contribute to physical observables. Hence, stationary
states are eigenvectors of the Hamiltonian operator. If, in addition, the
condition $\langle\hat{N}_{I}\rangle=\langle\hat{N}_{I^{\prime}}\rangle$ is
satisfied, the state is called static.

As we show now, static states give rise to static configurations. To prove
this, we first observe that the Hamiltonian can be defined as

$\hat{H}:=\int_{\Sigma_{t}}\hat{T}_{\mu\nu}k^{\mu}n^{\nu}d\gamma,$ (52)

with $n_{\nu}$ and $d\gamma$ given as in Eq. (4), and $k=\partial_{t}$ the
timelike Killing vector field associated with the static symmetry. Taking
$\Sigma_{t}$ to be a $t=const.$ hypersurface, such that $k^{\nu}=\alpha
n^{\nu}$, and using Eqs. (19), (26), (28), (29), (30), and (32) and
integration by parts, we obtain

$\hat{H}=\sum\limits_{I}\hat{N}_{I}\omega_{I}.$ (53)

In particular, it follows that each of the states $|N_{1},N_{2},\ldots\rangle$
is an eigenfunction of $\hat{H}$ with energy $E=\sum_{I}N_{I}\omega_{I}$. Any
stationary state can be written as a superposition of eigenfunctions of the
form $|N_{1},N_{2},\ldots\rangle$ which have the same definite value of the
energy $E$ (note that if $E$ is degenerate, this superposition might contain
several such eigenfunctions.) Since $\omega_{I}>0$ for each $I$, it follows
that $\hat{a}_{I}\hat{a}_{J}|N_{1},N_{2},\ldots\rangle$ is either zero, or is
an eigenfunction of $\hat{H}$ with energy smaller than $E$. This implies that
$\langle\hat{a}_{I}\hat{a}_{J}\rangle=0$ for all $I,J$. Similarly,
$\hat{a}_{I}|N_{1},N_{2},\ldots\rangle$ is either zero, or is an eigenfunction
of $\hat{H}$ with energy $E-\omega_{I}$, which implies that
$\langle\hat{a}_{I}^{\dagger}\hat{a}_{J}\rangle=0$ if
$\omega_{I}\neq\omega_{J}$. These conditions imply that $\rho$, $j_{k}$, and
$S_{ij}$ are time-independent. Therefore, if in addition the state is static,
we conclude that the configuration is static.

An interesting question is whether or not static configurations can only be
sourced by static states. For a real scalar field, we have not been able to
obtain a particular realization of a non-static state that is associated with
a static configuration. However, in the next section we provide such a
counterexample in the complex scalar field theory.

##  Appendix B Complex scalar fields

In this appendix, we generalize the static, spherically symmetric,
semiclassical real EKG system to the case of a complex scalar field. The main
interest of this extension is that, as we are going to demonstrate below (and
in stark contrast to the real case), when the field is complex there are
coherent states which are compatible with a static spacetime background.

The construction of the semiclassical, complex EKG system is very similar to
that presented in the main text, so for the sake of simplicity we will only
stress the main differences with respect to the real case. When expressed in
terms of the creation $(\hat{a}_{I}^{\dagger}$, $\hat{b}_{I}^{\dagger})$ and
annihilation $(\hat{a}_{I}$, $\hat{b}_{I})$ operators, a complex scalar field
takes the form

$\hat{\phi}(x)=\sum_{I}\left[\hat{a}_{I}f_{I}(x)+\hat{b}^{\dagger}_{I}f^{*}_{I}(x)\right],$
(54)

where $f_{1},f_{2},\ldots$ is an orthonormal set of basis functions of the
subspace $X_{+}$, and the creation and annihilation operators satisfy the
standard commutation relations
$[\hat{a}_{I},\hat{a}_{J}^{\dagger}]=[\hat{b}_{I},\hat{b}_{J}^{\dagger}]=\delta_{IJ}$,
with all other possible commutators vanishing. The Hilbert space $\mathscr{H}$
can be constructed by successive application of creation operators on the
vacuum state. The elements of this construction
$|N^{a}_{1},\ldots,N_{K}^{a},\ldots,N^{b}_{1},\ldots,N_{K}^{b},\ldots\rangle$
are the eigenstates of the particle and antiparticle number operators
$\hat{N}_{I}^{a}=\hat{a}_{I}^{\dagger}\hat{a}_{I}$ and
$\hat{N}_{I}^{b}=\hat{b}_{I}^{\dagger}\hat{b}_{I}$, respectively. As in the
real case, the coherent states are defined as the eigenstates of the
annihilation operators $\hat{a}_{I}$ and $\hat{b}_{I}$. Note that the main
difference with respect to the real field is that the complex theory contains
particles and antiparticles, although they are indistinguishable unless we
introduce interactions with other fields that break the degeneracy, which is
not the case here.

The stress energy-momentum tensor operator associated with a complex scalar
field takes the form

$\hat{T}_{\mu\nu}=(\nabla_{\mu}\hat{\phi})(\nabla_{\nu}\hat{\phi})^{\dagger}+(\nabla_{\nu}\hat{\phi})(\nabla_{\mu}\hat{\phi})^{\dagger}-g_{\mu\nu}\left[(\nabla_{\alpha}\hat{\phi})(\nabla^{\alpha}\hat{\phi})^{\dagger}+{m_{0}}^{2}\hat{\phi}\hat{\phi}^{\dagger}\right].$
(55)

Introducing the decomposition (54) into (55), we obtain

$\hat{T}_{\mu\nu}=\sum_{I,J}\left[\hat{a}_{I}\hat{a}^{\dagger}_{J}T_{\mu\nu}(f_{I},f_{J}^{*})+\hat{a}_{I}\hat{b}_{J}T_{\mu\nu}(f_{I},f_{J})+\hat{b}_{I}^{\dagger}\hat{a}^{\dagger}_{J}T_{\mu\nu}(f_{I}^{*},f_{J}^{*})+\hat{b}_{I}^{\dagger}\hat{b}_{J}T_{\mu\nu}(f_{I}^{*},f_{J})\right],$
(56)

with $T_{\mu\nu}(f_{I},f_{J})$ as in Eq. (20).

If the spacetime is static, the energy density (35a), momentum flux (35c), and
spatial stress tensor (35b) take the form

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\sum_{I,J}\frac{1}{2\sqrt{\omega_{I}\omega_{J}}}\left\\{\langle\hat{a}_{I}^{\dagger}\hat{a}_{J}\rangle\left[\left(+\frac{\omega_{I}\omega_{J}}{\alpha^{2}}+{m_{0}}^{2}\right)u_{I}^{*}u_{J}+(D_{k}u_{I})(D^{k}u^{*}_{J})\right]e^{+i(\omega_{I}-\omega_{J})t}\right.$
(57a) $\displaystyle\hskip
55.48277pt\left.+\langle\hat{a}_{I}\hat{b}_{J}\rangle\left[\left(-\frac{\omega_{I}\omega_{J}}{\alpha^{2}}+{m_{0}}^{2}\right)u_{I}u_{J}+(D_{k}u_{I})(D^{k}u_{J})\right]e^{-i(\omega_{I}+\omega_{J})t}\right.$
$\displaystyle\hskip
55.48277pt\left.+\langle\hat{b}^{\dagger}_{I}\hat{a}_{J}^{\dagger}\rangle\left[\left(-\frac{\omega_{I}\omega_{J}}{\alpha^{2}}+{m_{0}}^{2}\right)u_{I}^{*}u^{*}_{J}+(D_{k}u_{I}^{*})(D^{k}u^{*}_{J})\right]e^{+i(\omega_{I}+\omega_{J})t}\right.$
$\displaystyle\hskip
55.48277pt\left.+\langle\hat{b}^{\dagger}_{I}\hat{b}_{J}\rangle\left[\left(+\frac{\omega_{I}\omega_{J}}{\alpha^{2}}+{m_{0}}^{2}\right)u_{I}^{*}u_{J}+(D_{k}u_{I}^{*})(D^{k}u_{J})\right]e^{+i(\omega_{I}-\omega_{J})t}\right\\},$
$\displaystyle j_{k}$ $\displaystyle=$
$\displaystyle\sum_{I,J}\frac{1}{2\sqrt{\omega_{I}\omega_{J}}}\frac{i}{\alpha}\left\\{\langle\hat{a}^{\dagger}_{I}\hat{a}_{J}\rangle\left[-\omega_{J}(D_{k}u_{I}^{*})u_{J}+\omega_{I}u_{I}^{*}(D_{k}u_{J})\right]e^{-i(\omega_{I}-\omega_{J})t}\right.$
(57b) $\displaystyle\hskip
64.01869pt\left.+\langle\hat{a}_{I}\hat{b}_{J}\rangle\left[-\omega_{J}(D_{k}u_{I})u_{J}-\omega_{I}u_{I}(D_{k}u_{J})\right]e^{-i(\omega_{I}+\omega_{J})t}\right.$
$\displaystyle\hskip
64.01869pt\left.+\langle\hat{b}^{\dagger}_{I}\hat{a}_{J}^{\dagger}\rangle\left[+\omega_{J}(D_{k}u_{I}^{*})u_{J}^{*}+\omega_{I}u_{I}^{*}(D_{k}u_{J}^{*})\right]e^{+i(\omega_{I}+\omega_{J})t}\right.$
$\displaystyle\hskip
64.01869pt\left.+\langle\hat{b}^{\dagger}_{I}\hat{b}_{J}\rangle\left[-\omega_{J}(D_{k}u_{I}^{*})u_{J}+\omega_{I}u_{I}^{*}(D_{k}u_{J})\right]e^{+i(\omega_{I}-\omega_{J})t}\right\\},$
and $\displaystyle S_{ij}$ $\displaystyle=$
$\displaystyle\sum_{I,J}\frac{1}{2\sqrt{\omega_{I}\omega_{J}}}$
$\displaystyle\times$
$\displaystyle\left\\{\langle\hat{a}^{\dagger}_{I}\hat{a}_{J}\rangle\left[(D_{i}u_{I}^{*})(D_{j}u_{J})+(D_{j}u_{I}^{*})(D_{i}u_{J})-\gamma_{ij}\left(\left(-\frac{\omega_{I}\omega_{J}}{\alpha^{2}}+{m_{0}}^{2}\right)u_{I}^{*}u_{J}+(D_{k}u_{I}^{*})(D^{k}u_{J})\right)\right]e^{+i(\omega_{I}-\omega_{J})t}\right.$
$\displaystyle\left.+\langle\hat{a}_{I}\hat{b}_{J}\rangle\left[(D_{i}u_{I})(D_{j}u_{J})+(D_{j}u_{I})(D_{i}u_{J})-\gamma_{ij}\left(\left(+\frac{\omega_{I}\omega_{J}}{\alpha^{2}}+{m_{0}}^{2}\right)u_{I}u_{J}+(D_{k}u_{I})(D^{k}u_{J})\right)\right]e^{-i(\omega_{I}+\omega_{J})t}\right.$
$\displaystyle\left.+\langle\hat{b}^{\dagger}_{I}\hat{a}_{J}^{\dagger}\rangle\left[(D_{i}u_{I}^{*})(D_{j}u_{J}^{*})+(D_{j}u_{I}^{*})(D_{i}u_{J}^{*})-\gamma_{ij}\left(\left(+\frac{\omega_{I}\omega_{J}}{\alpha^{2}}+{m_{0}}^{2}\right)u_{I}^{*}u_{J}^{*}+(D_{k}u_{I}^{*})(D^{k}u_{J}^{*})\right)\right]e^{+i(\omega_{I}+\omega_{J})t}\right.$
$\displaystyle\left.+\langle\hat{b}^{\dagger}_{I}\hat{b}_{J}\rangle\left[(D_{i}u_{I}^{*})(D_{j}u_{J})+(D_{j}u_{I}^{*})(D_{i}u_{J})-\gamma_{ij}\left(\left(+\frac{\omega_{I}\omega_{J}}{\alpha^{2}}+{m_{0}}^{2}\right)u_{I}^{*}u_{J}+(D_{k}u_{I}^{*})(D^{k}u_{J})\right)\right]e^{-i(\omega_{I}-\omega_{J})t}\right\\},$

where we have used Wick’s ordering to obtain these expressions. Comparing them
with those of Eq. (51) we find that, apart from an overall factor of 2, they
differ in some of the expectation values which are quadratic in creation and
annihilation operators. Of course, if we choose a state with a definite number
of particles, for which
$\langle\hat{a}^{\dagger}_{I}\hat{a}_{J}\rangle=N_{I}\delta_{IJ}$ and
$\langle\hat{a}_{I}\hat{b}_{J}\rangle=\langle\hat{a}^{\dagger}_{I}\hat{b}^{\dagger}_{J}\rangle=\langle\hat{b}^{\dagger}_{I}\hat{b}_{J}\rangle=0$,
Eqs. (51) and (57) coincide, and we recover the relations in Eqs. (36). The
same holds true for a state with a definite number of antiparticles.
Furthermore, it is even possible to have a state with a definite number of
particles and antiparticles simultaneously; in that case, we only have to
replace $N_{I}$ by $N_{I}^{a}+N_{I}^{b}$ in Eqs. (36). This is because the
stress energy-momentum tensor sources gravity, and this tensor does not
differentiate between particles and antiparticles. However, if the field is
complex, it is still possible to guarantee a static source even if the system
is not in an eigenstate of the Hamiltonian operator. This is what happens for
coherent states that involve a single mode $I$ of the particle or antiparticle
sector. For these states
$\langle\hat{a}^{\dagger}_{I}\hat{a}_{J}\rangle=|\alpha_{I}|^{2}\delta_{IJ}$,
$\langle\hat{a}_{I}\hat{b}_{J}\rangle=\langle\hat{b}^{\dagger}_{I}\hat{a}_{J}^{\dagger}\rangle=\langle\hat{b}^{\dagger}_{I}\hat{b}_{J}\rangle=0$
(where in this case we have assumed that we only have particles), and it is
easy to convince oneself that we recover Eqs. (36) with $N_{I}$ replaced by
$|\alpha_{I}|^{2}$ and no sum over $I$. Note that this is not possible in the
real case, where the coefficients
$\langle\hat{a}_{I}\hat{a}_{J}\rangle=(\alpha_{i})^{2}\delta_{IJ}$ and
$\langle\hat{a}^{\dagger}_{I}\hat{a}_{J}^{\dagger}\rangle=(\alpha_{I}^{*})^{2}\delta_{IJ}$
of a coherent state give rise to a time dependency of the source terms that is
not compatible with a static metric. The main difference of these solutions
with respect to those with a definite number of particles is that the
expectation value of the scalar field evolves nontrivially in time.

## References

  * (1) D.J. Kaup, “Klein-Gordon geon,” Phys. Rev. 172, 1331 (1968)
  * (2) R. Ruffini and S. Bonazzola, “Systems of selfgravitating particles in general relativity and the concept of an equation of state,” Phys. Rev. 187, 1767-1783 (1969)
  * (3) M. Colpi, S.L. Shapiro and I. Wasserman, “Boson stars: Gravitational equilibria of selfinteracting scalar fields,” Phys. Rev. Lett. 57, 2485-2488 (1986)
  * (4) R. Friedberg, T.D. Lee and Y. Pang, “Mini-soliton stars,” Phys. Rev. D 35, 3640-3657 (1987)
  * (5) M. Gleiser, “Stability of boson stars,” Phys. Rev. D 38, 2376 (1988)
  * (6) T.D. Lee and Y. Pang, “Stability of mini-boson stars,” Nucl. Phys. B 315, 477 (1989)
  * (7) P. Bizon and A. Wasserman, “On existence of mini-boson stars,” Commun. Math. Phys. 215, 357 (2000) [arXiv:gr-qc/0002034]
  * (8) F.S. Guzman and L.A. Ureña-López, “Evolution of the Schrödinger-Newton system for a selfgravitating scalar field,” Phys. Rev. D 69, 124033 (2004) [arXiv:gr-qc/0404014]
  * (9) P.H. Chavanis and T. Harko, “Bose-Einstein condensate general relativistic stars,” Phys. Rev. D 86, 064011 (2012) [arXiv:1108.3986 [astro-ph.SR]
  * (10) P. Jetzer, “Boson stars,” Phys. Rept. 220, 163-227 (1992)
  * (11) F.E. Schunck and E.W. Mielke, “General relativistic boson stars,” Class. Quant. Grav. 20, R301-R356 (2003) [arXiv:0801.0307 [astro-ph]]
  * (12) S.L. Liebling and C. Palenzuela, “Dynamical boson stars,” Living Rev. Rel. 15, 6 (2012) [arXiv:1202.5809 [gr-qc]]
  * (13) E. Braaten and H. Zhang, “Axion stars,” Symmetry 12, 25 (2019) [arXiv:1810.11473]
  * (14) L. Visinelli, “Boson stars and oscillatons: A review,” Int. J. Mod. Phys. D 30, 2130006 (2021) [arXiv:2109.05481 [gr-qc]]
  * (15) T.D. Lee and Y. Pang, “Nontopological solitons,” Phys. Rept. 221, 251-350 (1992)
  * (16) R. Rajaraman, Solitons and Instantons: An Introduction to Solitons and Instantons in Quantum Field Theory, North-Holland (1987)
  * (17) N. Manton and P. Sutcliffe, Topological Solitons, Cambridge Monographs on Mathematical Physics. Cambridge University Press (2007) 508pp
  * (18) E.J. Weinberg, Classical Solutions in Quantum Field Theory: Solitons and Instantons in High Energy Physics, Cambridge Monographs on Mathematical Physics. Cambridge University Press (2012) 342pp
  * (19) Y.M. Shnir, Topological and Non-Topological Solitons in Scalar Field Theories, Cambridge Monographs on Mathematical Physics. Cambridge University Press (2018) 278pp
  * (20) C.J. Hogan and M.J. Rees, “Axion miniclusters,” Phys. Lett. B 205, 228-230 (1988)
  * (21) E.W. Kolb and I.I. Tkachev, “Axion miniclusters and Bose stars,” Phys. Rev. Lett. 71, 3051-3054 (1993) [arXiv:9303313 [hep-ph]]
  * (22) D.F. Torres, S. Capozziello and G. Lambiase, “A Supermassive scalar star at the galactic center?,” Phys. Rev. D 62, 104012 (2000) [arXiv:104012 [astro-ph]]
  * (23) F.G. Guzman, “Accretion disc onto boson stars: A Way to supplant black holes candidates,” Phys. Rev. D 73, 021501(R) (2006) [arXiv:0512081 [gr-qc]]
  * (24) C. Palenzuela, P. Pani, M. Bezares, V. Cardoso, L. Lehner and S. Liebling, “Gravitational wave signatures of highly compact boson star binaries,” Phys. Rev. D 96, 104058 (2017) [arXiv:1710.09432 [gr-qc]]
  * (25) V. Cardoso and P. Pani, “Testing the nature of dark compact objects: a status report,” Living Rev. Rel. 22, 4 (2019) [arXiv:1904.05363 [gr-qc]]
  * (26) S-J. Sin, “Late time cosmological phase transition and galactic halo as Bose liquid,” Phys. Rev. D 50, 3650-3654 (1994) [arXiv:9205208 [hep-ph]]
  * (27) H-Y. Schive, T. Chiueh and T. Broadhurst, “Cosmic structure as the quantum interference of a coherent dark wave,” Nature Phys. 10, 496-499 (2014) [arXiv:1406.6586 [astro-ph.GA]]
  * (28) H.Y. Schive, M.H. Liao, T.P. Woo, S.K. Wong, T. Chiueh, T. Broadhurst and W.Y.P. Hwang, “Understanding the core-halo relation of quantum wave dark matter, $\psi$DM, from 3D simulations,” Phys. Rev. Lett. 113, 261302 (2014) [arXiv:1407.7762]
  * (29) B. Schwabe, J.C. Niemeyer and J.F. Engels, “Simulations of solitonic core mergers in ultra-light axion dark matter cosmologies,” Phys. Rev. D 94, 043513 (2016) [arXiv:1606.05151]
  * (30) J. Veltmaat and J.C. Niemeyer, “Cosmological particle-in-cell simulations with ultra-light axion dark matter,” Phys. Rev. D 94, 123523 (2016) [arXiv:1608.00802]
  * (31) X. Du, C. Behrens, J.C. Niemeyer and B. Schwabe, “The core-halo mass relation of ultra-light axion dark matter from merger history,” Phys. Rev. D 95, 043519 (2017) [arXiv:1609.09414]
  * (32) A.X. Gonzalez-Morales, D.J.E. Marsh, J. Peñarrubia and L.A. Ureña-López, “Unbiased constraints on ultralight axion mass from dwarf spheroidal galaxies,” Mon. Not. Roy. Astron. Soc. 472, 1346–1360 (2017) [arXiv:1609.05856 [astro-ph.CO]]
  * (33) R.D. Peccei and H.R. Quinn, “CP conservation in the presence of pseudoparticles,” Phys. Rev. Lett. 38, 1440 (1977)
  * (34) J. Preskill, M.B. Wise and F. Wilczek, “Cosmology of the invisible axion,” Phys. Lett. B 120, 127 (1983)
  * (35) L.F. Abbott and P. Sikivie, “A cosmological bound on the invisible axion,” Phys. Lett. B 120, 133 (1983)
  * (36) M. Dine and W. Fischler, “The not so harmless axion,” Phys. Lett. B 120, 137 (1983)
  * (37) M.P Hertzberg, M. Tegmark and F. Wilczek, “Axion cosmology and the energy scale of inflation,” Phys. Rev. D 78, 083507 (2008)
  * (38) P. Sikivie, “Axion Cosmology,” Lect. Notes Phys. 741, 19-50 (2008) [astro-ph/0610440]
  * (39) W. Hu, R. Barkana and A. Gruzinov,“Cold and fuzzy dark matter,” Phys. Rev. Lett. 85, 1158-1161 (2000) [arXiv:astro-ph/0003365].
  * (40) T. Matos and L. A. Ureña-López,“A Further analysis of a cosmological model of quintessence and scalar dark matter,” Phys. Rev. D 63, 063506 (2001) [arXiv:astro-ph/0006024]
  * (41) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper and J. March-Russell, “String axiverse,” Phys. Rev. D 81, 123530 (2010) [arXiv:0905.4720 [hep-th]]
  * (42) A. Suárez, V.H. Robles and T. Matos, “A review on the scalar field/Bose-Einstein condensate dark matter model,” Astrophys. Space Sci. Proc. 38, 107–142 (2014) [arXiv:1302.0903 [astro-ph.CO]]
  * (43) D.J.E. Marsh, “Axion cosmology,” Phys. Rept. 643, 1-79 (2016) [arXiv:1510.07633 [astro-ph.CO]]
  * (44) L. Hui, J.P. Ostriker, S. Tremaine, E. Witten, “Ultralight scalars as cosmological dark matter,” Phys. Rev. D 95, 043541 (2017) [arXiv:1610.08297 [astro-ph.CO]]
  * (45) J.C. Niemeyer, “Small-scale structure of fuzzy and axion-like dark matter,” Prog. Part. Nucl. Phys. 113, 103787 (2020) [arXiv:1912.07064]
  * (46) L.A. Ureña-López, “Brief review on scalar field dark matter models,” Front. Astron. Space Sci. 12, (2019) [https://doi.org/10.3389/fspas.2019.00047]
  * (47) E.G.M. Ferreira, “Ultra-light dark matter,” Astron. Astrophys. Rev. 29, 7 (2021) [arXiv:2005.03254 [astro-ph.CO]]
  * (48) J. Ho, F. C. Khanna, C. H. Lee, “Boson stars with self-interacting quantum scalar fields,” (2003) [arXiv:gr-qc/0207073]
  * (49) T. Matos and L.A. Ureña-López, “Flat rotation curves in scalar field galaxy halos,” Gen. Rel. Grav. 39, 1279-1286 (2007)
  * (50) A. Bernal, J. Barranco, D. Alic and C. Palenzuela, “Multi-state boson stars,” Phys. Rev. D 81, 044031 (2010) [arXiv:0908.2435 [gr-qc]]
  * (51) L.A. Uneña-López and A. Bernal, “Bosonic gas as a galactic dark matter halo,” Phys. Rev. D 82, 123535 (2010) [arXiv:1008.1231 [gr-qc]]
  * (52) J. Barranco and A. Bernal, “Self-gravitating system made of axions,” Phys. Rev. D 83, 043525 (2011) [arXiv:1001.1769 [astro-ph.CO]]
  * (53) B. Berczi, P.M. Saffin and S.Y. Zhou, “Gravitational collapse with quantum fields,” Phys. Rev. D 104, L041703 (2021) [arXiv:2010.10142 [gr-qc]]
  * (54) B. Berczi, P.M. Saffin and S.Y. Zhou, “Gravitational collapse of quantum fields and Choptuik scaling,” JHEP 02, 183 (2022) [arXiv:2111.11400 [hep-th]]
  * (55) J. Guenther, C. Hoelbling and L. Varnhorst, “Semiclassical gravitational collapse of a radially symmetric massless scalar quantum field,” Phys. Rev. D 105, 105010 (2022) [arXiv:2010.13215 [gr-qc]]
  * (56) J. Arrechea, C. Barceló, R. Carballo-Rubio and L. J. Garay, “Semiclassical relativistic stars,” Sci. Rep. 12, 15958 (2021) [arXiv:2110.15808 [gr-qc]]
  * (57) B.-L. Hu and E. Verdaguer, Semiclassical and Stochastic Gravity, Cambridge University Press (2020) 599pp
  * (58) See e.g. Birrelll and Davies Birrell1984  Chapter 6 “Stress-tensor renormalization”; Wald Wald1994  Section 4.6 “The stress-energy tensor”; Mukhanov and Winitzki Mukhanov2007  Part 2 “Path integrals and vacuum polarization”; and Parker and Toms Parker2009  Chapter 3 “Expectation values quadratic in fields”. See also Ref. Burgess:2003jk  for a general discussion on quantum gravity at accessible scales.
  * (59) N. D. Birrell and P. C. W. Davies, Quantum Fields in Curved Space, Cambridge University Press (1984) 352pp
  * (60) R. M. Wald, Quantum Field Theory in Curved Spacetime and Black Hole Thermodynamics, University Of Chicago Press (1994) 220pp
  * (61) V. Mukhanov and S. Winitzki, Introduction to Quantum Effects in Gravity, Cambridge University Press (2007) 284pp
  * (62) L. Parker and D. Toms, Quantum Field Theory in Curved Spacetime: Quantized Fields and Gravity, Cambridge University Press (2009) 472pp
  * (63) C.P. Burgess, “Quantum gravity in everyday life: General relativity as an effective field theory,” Living Rev. Rel. 7, 5–56 (2004) [arXiv:gr-qc/0311082]
  * (64) R.D. Jordan, “Effective field equations for expectation values,” Phys. Rev. D 33, 444-454 (1986)
  * (65) E. Calzetta and B.L. Hu, “Closed time path functional formalism in curved space-time: Application to cosmological back reaction problems,” Phys. Rev. D 35, 495 (1987)
  * (66) J.P. Paz, “Anisotropy dissipation in the early universe: Finite temperature effects reexamined,” Phys. Rev. D 41, 1054-1066 (1990)
  * (67) A. Campos and E. Verdaguer, “Semiclassical equations for weakly inhomogeneous cosmologies,” Phys. Rev. D 49, 1861-1880 (1994)
  * (68) A. Diez-Tejedor and D. Sudarsky, “Towards a formal description of the collapse approach to the inflationary origin of the seeds of cosmic structure”, JCAP, 07, 045 (2012)
  * (69) L. Parker and S.A. Fulling, “Adiabatic regularization of the energy momentum tensor of a quantized field in homogeneous spaces,” Phys. Rev. D 9, 341-354 (1974)
  * (70) W. Pauli and F. Villars, “On the invariant regularization in relativistic quantum theory,” Rev. Mod. Phys. 21, 434-444 (1949)
  * (71) S. Weinberg, “Ultraviolet divergences in cosmological correlations,” Phys. Rev. D 83, 063508 (2011) [arXiv:1011.1630 [hep-th]]
  * (72) C. Armendariz-Picon, “On the expected production of gravitational waves during preheating,” JCAP 08, 012 (2019) [arXiv:1905.05724 [astro-ph.CO]].
  * (73) E.H. Lieb, “Existence and Uniqueness of the Minimizing Solution of Choquard’s Nonlinear Equation”, Studies in Applied Mathematics 57 (1977) 93–105.
  * (74) O. Kavian and S. Mischler, “A global approach to the Schrödinger-Poisson system: An existence result in the case of infinitely many states”, J. Math. Pures Appl. 104 (2015) 942-964.
  * (75) M. Alcubierre, J. Barranco, A. Bernal, J.C. Degollado, A. Diez-Tejedor, M. Megevand, D. Núñez and O. Sarbach, “$\ell$-Boson stars,” Class. Quant. Grav. 35, 19LT01 (2018) [arXiv:1805.11488 [gr-qc]]
  * (76) M. Alcubierre, J. Barranco, A. Bernal, J.C. Degollado, A. Diez-Tejedor, M. Megevand, D. Núñez and O. Sarbach, “Dynamical evolutions of $\ell$-boson stars in spherical symmetry,” Class. Quant. Grav. 36, 215013 (2019) [arXiv:1906.08959 [gr-qc]]
  * (77) M. Alcubierre, J. Barranco, A. Bernal, J.C. Degollado, A. Diez-Tejedor, M. Megevand, D. Núñez and O. Sarbach, “On the linear stability of $\ell$-boson stars with respect to radial perturbations,” Class. Quant. Grav. 38, 174001 (2021) [arXiv:2103.15012 [gr-qc]]
  * (78) M. Alcubierre, J. Barranco, A. Bernal, J. C. Degollado, A. Diez-Tejedor, V. Jaramillo, M. Megevand, D. Núñez, O. Sarbach, “Extreme $\ell$-boson stars,” Class. Quant. Grav. 39, 094001 (2022) [arXiv:2112.04529 [gr-qc]]
  * (79) V. Jaramillo, N. Sanchis-Gual, J. Barranco, A. Bernal, J. .C. Degollado, C. Herdeiro, M. Megevand, D. Núñez, “Head-on collisions of $\ell$-boson stars,” Phys. Rev. D 105, 104057 (2022) [arXiv:2202.00696 [gr-qc]]
  * (80) V. Jaramillo Pérez, Límite de campo débil para el campo escalar autogravitante, Master’s Thesis, Universidad Nacional Autónoma de México (2019).
  * (81) F. S. Guzman and L. A. Ureña-López, “Gravitational atoms: General framework for the construction of multistate axially symmetric solutions of the Schrödinger-Poisson system,” Physical Review D 101, 8 (2020) [arXiv:1912.10585 [astro-ph.GA]]
  * (82) E. Chávez Nambo, Sobre la existencia de estrellas de bosones newtonianas con momento angular en simetría esférica, Master’s Thesis, Universidad Michoacana de San Nicolás de Hidalgo (2021)
  * (83) J. Vicens, J. Salvado and J. Miralda-Escude, “Bosonic dark matter halos: excited states and relaxation in the potential of the ground state,” (2018) [arXiv:1802.10513 [astro-ph.CO]]
  * (84) S-C. Lin, H-Y. Schive, S-K. Wong and T. Chiueh, “Self-consistent construction of virialized wave dark matter halos,” Phys. Rev. D 97, 103523 (2018) [arXiv:1801.02320 [astro-ph.CO]]
  * (85) E. Seidel and W.M. Suen, “Oscillating soliton stars,” Phys. Rev. Lett. 66, 1659 (1991)
  * (86) S.R. Coleman, “Q-balls,” Nucl. Phys. B 262, 263 (1985) [Addendum: Nucl. Phys. B 269, 744 (1986)]
  * (87) F. Di Giovanni, S. Fakhry, N. Sanchis-Gual, J. C. Degollado and J. A. Font, “A stabilization mechanism for excited fermion–boson stars,” Class. Quant. Grav. 38, (2021) no.19, 194001 [arXiv:2105.00530 [gr-qc]]
  * (88) C. Herdeiro and E. Radu, “On the classicality of bosonic stars,” Int. J. Mod. Phys. D 31, 2242022 (2022) [arXiv:2205.05395 [gr-qc]]
  * (89) D.R. Hartree, “The wave mechanics of an atom with a non-coulomb central field,” Math. Proc. Camb. Philos. Soc. 24, 111 (1928)
  * (90) A.H. Guth, M.P. Hertzberg and C. Prescod-Weinstein, “Do dark matter axions form a condensate with long-range correlation?,” Phys. Rev. D 92, 103513 (2015) [arXiv:1412.5930 [astro-ph.CO]]
  * (91) P.W. Graham and R. Surjeet, “New observables for direct detection of axion dark matter,” Phys. Rev. D 88, 035023 (2013) [arXiv:1306.6088 [hep-ph]]
  * (92) D. Carney, A. Hook, Z. Liu, J.M. Taylor and Y. Zhao, “Ultralight dark matter detection with mechanical quantum sensors,” New J. Phys. 23, 023041 (2021) [arXiv:1908.04797 [hep-ph]]
  * (93) C.M. Donohue, S. Gardner and W. Korsch, “LC circuits for the direct detection of ultralight dark matter candidates,” (2021) [arXiv:2109.08163 [hep-ph]]
  * (94) Y.-D. Tsai, J. Eby, and M.S. Safronova, Marianna S, “SpaceQ – Direct detection of ultralight dark matter with space quantum sensors,” (2021) [arXiv:2112.07674 [hep-ph]]
  * (95) A. Perez, H. Sahlmann and D. Sudarsky, “On the quantum origin of the seeds of cosmic structure,” Class. Quant. Grav. 23, 2317–2354 (2006) [arXiv:gr-qc/0508100]
  * (96) D. Sudarsky, “Shortcomings in the understanding of why cosmological perturbations look classical,” Int. J. Mod. Phys. D 20, 509–552 (2011) [arXiv:gr-qc/0906.0315]
  * (97) A. Aguirre and A. Diez-Tejedor, “Asymmetric condensed dark matter,” JCAP 04, 019 (2016) [arXiv:1502.07354 [astro-ph.CO]]
  * (98) P. Sikivie and E.M. Todarello, “Duration of classicality in highly degenerate interacting Bosonic systems,” Phys. Lett. B 770, 331–334 (2017) [arXiv:1607.00949 [hep-ph]]
  * (99) M.P. Hertzberg, “Quantum and classical behavior in interacting bosonic systems,” JCAP 11, 037 (2016) [arXiv:1609.01342 [hep-ph]]
  * (100) I.J. Allali and M.P. Hertzberg, “Decoherence from General Relativity,” Phys. Rev. D 103, 104053 (2021) [arXiv:2012.12903 [gr-qc]]
  * (101) E.E. Flanagan and R.M. Wald, “Does back reaction enforce the averaged null energy condition in semiclassical gravity?” Phys. Rev. D 54, 6233-6283 (1996) [gr-qc/9602052]
  * (102) T. Jacobson, “Introduction to quantum fields in curved space-time and the Hawking effect,” contribution to School on Quantum Gravity, pages 39–89 (2003) [gr-qc/0308048 [gr-qc]]
  * (103) A. Ashtekar and A. Magnon, “Quantum fields in curved space-times,” Proc. Roy. Soc. Lond. A, 346, 375–394, 1975.
  * (104) J.J. Sakurai, Modern Quantum Mechanics, Addison-Wesley Publishing Company (1994) 500pp
  * (105) W.M. Zhang, D.H. Feng and R. Gilmore, “Coherent states: Theory and some applications,” Rev. Mod. Phys. 62, 867 (1990)
  * (106) B.C. Sanders, “Review of entangled coherent states,” J. Phys. A: Math. Theor. 45, 244002 (2012)
  * (107) R.J. Glauber, “Coherent and incoherent states of radiation field,” Phys. Rev. 131, 2766-2788 (1963)
  * (108) C. Armendariz-Picon, “On the expected backreaction during preheating,” JCAP 05, 035 (2020) [arXiv:2003.01542 [gr-qc]]
  * (109) B. A. Juarez-Aubry, T. Miramontes and D. Sudarsky, “Semiclassical theories as initial value problems,” J. Math. Phys. 61, 032301 (2020) [arXiv:1907.09960 [math-ph]]
  * (110) B. A. Juarez-Aubry, “Semiclassical gravity in static spacetimes as a constrained initial value problem,” Annales Henri Poincaré 23, 4, 1451-1487 (2022) [arXiv:2011.05947 [gr-qc]]
  * (111) B. A. Juarez-Aubry and S. L. Modak, “Semiclassical gravity with a conformally covariant field in globally hyperbolic spacetimes,” J. Math. Phys. 63, 9, 092303 (2022) [arXiv:2110.01719 [math-ph]]
  * (112) B. A. Juarez-Aubry, B. Kay, T. Miramontes and D. Sudarsky, “On the initial value problem for semiclassical gravity without and with quantum state collapses,” (2022) [arXiv:2205.11671 [gr-qc]]
  * (113) P. Canate, E. Ramirez and D. Sudarsky, “Semiclassical self consistent treatment of the emergence of seeds of cosmic structure. The second order construction,” JCAP 08, 043 (2018); Erratum: JCAP 10, E01 (2018) [arXiv:1802.02238 [gr-qc]]
  * (114) M. Reed and B. Simon, Methods of Modern Mathematical Physics, Vol. II: Fourier Analysis, Self-Adjointness, Academic Press, San Diego (1980)
  * (115) M. Reed and B. Simon, Methods of Modern Mathematical Physics, Vol. I: Functional Analysis, Academic Press, San Diego (1980)
  * (116) M. Reed and B. Simon, Methods of Modern Mathematical Physics, Vol. IV: Analysis of Operators, Academic Press, San Diego (1980)
  * (117) A. Much and R. Oeckl, “Self-Adjointness in Klein-Gordon theory on globally hyperbolic spacetimes,” Math. Phys. Anal. Geom. 24, 5 (2021) [arXiv:1804.07782 [math-ph]]
  * (118) M. Alcubierre, Introduction to 3+1 Numerical Relativity, Oxford University Press (2008) 444pp
  * (119) E.W. Mielke and R. Scherzer, “Geon type solutions of the nonlinear Heisenberg-Klein-Gordon equation,” Phys. Rev. D 24, 2111 (1981)
  * (120) H. Andréasson, “The Einstein-Vlasov system/kinetic theory,” Living Rev. Rel. 14, 4 (2011) [arXiv:1106.1367 [gr-qc]]
  * (121) I.M. Moroz, R. Penrose and P. Tod, “Spherically-symmetric solutions of the Schrödinger-Newton equations”, Class. and Quantum Grav. 15, (1998), 2733–2742.
  * (122) I. De Martino, T. Broadhurst, S.-H.H. Tye, T. Chiueh, and H.-Y. Schive, “Dynamical evidence of a solitonic core of $10^{9}M_{\odot}$ in the Milky Way,” Phys. Dark Univ. 28, 100503 (2020) [arXiv:1807.08153 [astro-ph.GA]]
  * (123) Pierre-Alan Duc et. al. “The Atlas 3D project XXIX. The new look of early-type galaxies and surrounding fields disclosed by extremely deep optical images,” MNRAS 446, 120 (2015) [arXiv:1410.0981 [astro-ph.GA]]