Non-linear dynamics of spinning bosonic stars: formation and stability

TL;DR

This study uses numerical simulations to compare the formation and stability of spinning bosonic stars made of scalar and vector fields, revealing distinct stability behaviors and potential implications for astrophysical objects.

Contribution

It provides the first detailed numerical analysis of the non-linear dynamics and stability of spinning bosonic stars, highlighting differences between scalar and Proca (vector) cases.

Findings

Scalar boson stars are unstable and eject angular momentum.

Proca stars are stable and resist large perturbations.

Formation of scalar stars involves transient states and instabilities.

Abstract

We perform numerical evolutions of the fully non-linear Einstein-(complex, massive)Klein-Gordon and Einstein-(complex)Proca systems, to assess the formation and stability of spinning bosonic stars. In the scalar/vector case these are known as boson/Proca stars. Firstly, we consider the formation scenario. Starting with constraint-obeying initial data, describing a dilute, axisymmetric cloud of spinning scalar/Proca field, gravitational collapse towards a spinning star occurs, via gravitational cooling. In the scalar case the formation is transient, even for a non-perturbed initial cloud; a non-axisymmetric instability always develops ejecting all the angular momentum from the scalar star. In the Proca case, by contrast, no instability is observed and the evolutions are compatible with the formation of a spinning Proca star. Secondly, we address the stability of an existing star, a stationary solution of the field equations. In the scalar case, a non-axisymmetric perturbation develops collapsing the star to a spinning black hole. No such instability is found in the Proca case, where the star survives large amplitude perturbations; moreover, some excited Proca stars decay to, and remain as, fundamental states. Our analysis suggests bosonic stars have different stability properties in the scalar/vector case, which we tentatively relate to their toroidal/spheroidal morphology. A parallelism with instabilities of spinning fluid stars is briefly discussed.

Tables

Table 1. Table 1: Physical properties of some illustrative SBSs. The second row identifies if they are scalar ( S 𝑆 S ) or vector/Proca ( P 𝑃 P ). All solutions have m = 1 𝑚 1 m=1 and none has an ergo-region.

Configuration	$1_S$	$2_S$	$1_P$	$2_P$	$3_P$	$4_P$	$5_P$
type (S or P)	S	S	P	P	P	P	P
$n$	0	0	0	0	1	1	1
$\omega$	0.90	0.83	0.95	0.90	0.95	0.90	0.85
$M$	1.119	1.281	0.534	0.726	1.149	1.456	1.564
$J$	1.153	1.338	0.543	0.750	1.171	1.500	1.622

Equations

$${\mathcal{L}}_m=-{\partial }^{\alpha }\varphi {\partial }_{\alpha }\bar{\varphi }-{\mu }^2\varphi \bar{\varphi },{\mathcal{L}}_m=-\frac{{\mathcal{F}}_{\alpha \beta }{\bar{\mathcal{F}}}^{\alpha \beta }}{4}-\frac{{\mu }^2}{2}{\mathcal{A}}_{\alpha }{\bar{\mathcal{A}}}^{\alpha },$$

$$\varphi (t,r,\theta ,\phi )=Ar{e}^{-\frac{r^2}{{\sigma }^2}}\sin \theta e^{i(\phi -\omega t)},$$

$$g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=-{\alpha }^{2}d{t}^2+{\gamma }_{ij}(d{x}^i+{\beta }^{i}dt)(d{x}^j+{\beta }^{j}dt),$$

$${\gamma }_{ij}={\Psi }^4{\tilde{\gamma }}_{ij},$$

$$\Pi =-\frac{1}{\alpha }({\partial }_t-{\beta }^i{\partial }_i)\varphi .$$

$$\varphi (t,r,\theta ,\phi )=R(r)Y_{11}(\theta ,\phi )e^{-i\omega t},$$

$$\Pi =-\frac{i}{\alpha }(\omega +{\beta }^{\phi })\varphi ,$$

$$D_i{E}^i=-{\mu }^2{\mathcal{X}}_{\varphi }.$$

$${\nabla }^{2}V={\mu }^2{\mathcal{X}}_{\varphi }.$$

$${\mathcal{X}}_{\varphi }(t,r,\theta ,\phi )=R(r)Y_{11}(\theta ,\phi )e^{-i\omega t},$$

$$\frac{1}{r}{\partial }_r^2[rV(r)]-\frac{2}{r^2}V(r)={\mu }^{2}R(r).$$

$$\begin{split}V(r)=&\frac{A\,\sigma e^{-\frac{r^{2}}{\sigma^{2}}}}{24r^{2}}\biggl{(}-e^{\frac{r^{2}}{\sigma^{2}}}\biggl{[}2\sqrt{\pi}\sigma^{2}r^{3}+8\sigma^{5}\\ &-2\sqrt{\pi}\sigma^{3}r^{2}\,\rm{Erf}\left(\frac{r}{\sigma}\right)\biggr{]}+8r^{2}\sigma^{3}+8\sigma^{5}\biggr{)}\ .\end{split}$$

$${\partial }_t{\mathcal{X}}_{\varphi }=-{\mathcal{X}}^i{D}_i\alpha -\alpha D_i{\mathcal{X}}^i+{\beta }^{\phi }{\partial }_{\phi }{\mathcal{X}}_{\varphi },$$

$$D_i(\alpha {\mathcal{X}}^i)=i(\omega +{\beta }^{\phi }){\mathcal{X}}_{\varphi },$$

$${\mathcal{X}}^i=\frac{i}{\alpha }(\omega +{\beta }^{\phi })E^i.$$

$${\mathcal{X}}_{\varphi }(r,\theta ,\phi )=$$

$$E^r(r,\theta ,\phi )=$$

$$E^{\theta }(r,\theta ,\phi )=$$

$$E^{\phi }(r,\theta ,\phi )=$$

$${\mathcal{X}}_r(r,\theta ,\phi )=$$

$${\mathcal{X}}_{\theta }(r,\theta ,\phi )=$$

$${\mathcal{X}}_{\phi }(r,\theta ,\phi )=$$

$${\hat{E}}^r(r,\theta )=$$

$$\displaystyle\ 4r^{2}\sigma^{2}+4\sigma^{4}]-\sqrt{\pi}r^{3}\sigma\,\rm{Erf}\left(\frac{r}{\sigma}\right)\biggr{)}\sin{\theta}\ ,$$

$${\hat{E}}^{\theta }(r,\theta )=$$

$$\displaystyle-2\sqrt{\pi}\sigma^{3}r^{2}\,\rm{Erf}\left(\frac{r}{\sigma}\right)\biggr{]}+8r^{2}\sigma^{3}+8\sigma^{5}\biggr{)}\cos{\theta}\ ,$$

$${\hat{E}}^{\phi }(r,\theta )=$$

$$\displaystyle-2\sqrt{\pi}\sigma^{3}r^{2}\,\rm{Erf}\left(\frac{r}{\sigma}\right)\biggr{]}+8r^{2}\sigma^{3}+8\sigma^{5}\biggr{)}\cos{\theta}\ ,$$

$${\hat{\mathcal{X}}}_i(r,\theta )=$$

Full text

Introduction. Recent data from gravitational-wave astronomy Abbott et al. (2018), as well as from electromagnetic VLBI observations near galactic centres Gravity Collaboration et al. (2018); Akiyama et al. (2019), support the black hole (BH) hypothesis: BHs commonly populate the Cosmos, with masses spanning a range of (at least) 10 orders of magnitude. Yet, the elusiveness of the event horizon, the defining property of a BH, rules out an observational “proof” of their existence. Considering, thus, models of BH mimickers is a valuable tool to understand the uniqueness of BH phenomenology.

Within the landscape of BH mimickers, bosonic stars (BSs) are particularly well motivated. They arise in simple and physically sound field theoretical models: complex, massive, bosonic fields (scalar Kaup (1968); Ruffini and Bonazzola (1969) or vector Brito et al. (2016)) minimally coupled to Einstein’s gravity. Dynamically, moreover, static, spherical BSs, are viable; for some range of parameters, the lowest energy stars - the fundamental family (FF) - have a formation mechanism Seidel and Suen (1994); Di Giovanni et al. (2018) and are perturbatively stable Gleiser and Watkins (1989); Lee and Pang (1989); Brito et al. (2016); Sanchis-Gual et al. (2017). The properties and phenomenology of such static BSs have been considered at length (see $e.g.$ the reviews Schunck and Mielke (2003); Liebling and Palenzuela (2012)), including dynamical situations such as orbiting binaries, from which gravitational waveforms have been extracted Palenzuela et al. (2007, 2008); Sanchis-Gual et al. (2019). These studies unveiled a close parallelism in the phenomenology of spherical BSs, regardless of their scalar or vector nature.

Astrophysically, however, rotation is ubiquitous and should, thus, be included in more realistic models of BSs. Both scalar Schunck and Mielke (1998); Yoshida and Eriguchi (1997); Kleihaus et al. (2008) and vector Brito et al. (2016); Herdeiro et al. (2016) axisymmetric, spinning BSs (SBSs) have been constructed and some of their phenomenology has been studied Vincent et al. (2016); Meliani et al. (2015). Yet, their dynamical and stability properties, a key aspect of their physical viability, have remained essentially unexplored - see the discussion in Herdeiro and Radu (2014).

In this letter we describe the dynamical properties of SBSs, obtained from fully non-linear numerical simulations of the corresponding Einstein-matter system. We provide evidence that scalar SBSs in the FF are prone to a non-axisymmetric instability. Thus, such stars are transient states, in a dynamical formation scenario. Assuming an already formed scalar SBS, on the other hand, it collapses into a BH after a non-axially symmetric instability develops. Vector SBSs (aka spinning Proca stars), by contrast, are dynamically robust. In the formation scenario we find no evidence of an instability. In agreement, for already formed vector SBSs we observe that: $(i)$ even large perturbations are dissipated away; and $(ii)$ some stars in excited families decay to the FF where they remain. This suggests scalar/vector SBSs have different dynamical properties and viability, and their toroidal/spheroidal morphology provides a suggestive interpretation.

SBSs as stationary solutions. Scalar and vector BSs, with and without spin, arise as equilibrium states in models with Lagrangian density $\mathcal{L}=R/(16\pi G)+\mathcal{L}_{m}$, where $R$ is the Ricci scalar, $G$ is Newton’s constant and

[TABLE]

describe the scalar and vector cases, respectively. The scalar ($\phi$) and vector ($\mathcal{A}_{\alpha}$) fields are complex-valued, with conjugation denoted by an overbar, both with mass $\mu$. As usual, $\mathcal{F}=d\mathcal{A}$. Henceforth, units with $G=1=c=\mu$ are used.

Scalar SBSs were first constructed numerically in Schunck and Mielke (1998); Yoshida and Eriguchi (1997), as asymptotically flat, stationary and axisymmetric solutions of the above Einstein-Klein-Gordon system. They are a “mass torus” in general relativity - Fig. 1 (left panel). Scalar SBSs form a discrete set of families of continuous solutions. Each family is labelled by two integers: $m$, the azimuthal winding number and $n$, the node (or overtone) number, see $e.g.$ Schunck and Mielke (2003); Kleihaus et al. (2008); Liebling and Palenzuela (2012); Grandclement et al. (2014); Herdeiro and Radu (2015). The FF, which has the lowest energy, has $(m,n)=(1,0)$. Fixing the family, $i.e.$ $(m,n)$, SBS are characterised by their total mass $M$ and angular momentum $J$. They form a one-dimensional set, often labelled by $M$, and oscillation frequency, $\omega$. The (dynamically) most interesting solutions occur in between the Newtonian limit, $\omega\rightarrow 1$ and $M\rightarrow 0$, and the maximal mass solution. The latter occurs for $\omega\rightarrow\omega_{\rm Mmax}$ ($\simeq 0.775$ for the FF) and the ADM mass becomes highest, $M\rightarrow M_{\rm max}$ ($\simeq 1.315$ for the FF). In Table 1 we list the properties of two illustrative scalar SBSs used in the simulations below.

Vector SBSs were first reported as excited states $(n=1)$ in Brito et al. (2016); Herdeiro et al. (2016). The FF was considered in Herdeiro and Radu (2017); Herdeiro et al. (2019). The aforementioned description for scalar SBSs applies, mutatis mutandis. An important distinction, however, is that the energy distribution is now spheroidal, rather than toroidal Herdeiro et al. (2019) - Fig. 1 (right panel). Moreover, for the FF, $\omega_{\rm Mmax}\simeq 0.562$ and $M_{\rm max}\simeq 1.125$ Herdeiro et al. (2019). For the excited family with $(m,n)=(1,1)$, $\omega_{\rm Mmax}\simeq 0.839$ and $M_{\rm max}\simeq 1.568$ Brito et al. (2016). In Table 1 we list the properties of two [three] representative vector SBSs, in the FF [$(m,n)=(1,1)$ family].

Dynamical formation of SBSs. In the spherical case, numerical simulations established that both scalar Seidel and Suen (1994) and vector Di Giovanni et al. (2018) BSs form dynamically, from a spherical “cloud” of dilute scalar or vector field. The cloud collapses due to its self-gravity. The ejection of energetic scalar or vector “particles”, dubbed gravitational cooling, allows the formation of a compact object.

For studying the formation of SBSs, with $m=1$, the Hamiltonian, momentum and (in the vector case) Gauss constraint are solved by appropriately choosing a Gaussian radial dependence for the key variables, together with a non-spherical profile (see Appendix A). For the scalar case the “matter” initial data is:

[TABLE]

where $A$, $\sigma$ are constants and $e^{-i\omega t}$ is the harmonic dependence. Besides this unperturbed initial data, we also evolve perturbed initial data of two types: replacing in (2) $e^{i\varphi}\rightarrow e^{i\varphi}(1+A_{1}\cos(2\varphi))$; or, alternatively, replacing $e^{i\varphi}\rightarrow e^{i\varphi}+A_{2}e^{2i\varphi}$. $A_{1},A_{2}$ are the amplitudes of the perturbations.

Fully non-linear numerical evolutions of the Einstein-matter system using this initial data were carried out with the Einstein Toolkit EinsteinToolkit ; Löffler (2012) - see Appendix B. Two choices of $A$ were considered, both of which yield global data for the scalar cloud $(M_{\rm{Sc}},J_{\rm{Sc}})$ close to that of equilibrium scalar SBS solutions. The first/second choices give $M_{\rm{Sc}}^{(1)}\sim 0.46\sim J_{\rm{Sc}}^{(1)}$, and $M_{\rm{Sc}}^{(2)}\sim 0.89\sim J_{\rm{Sc}}^{(2)}$, respectively. We have run simulations with both perturbed and unperturbed initial data, with $A_{1}=0,0.001,0.01,0.05$ and $A_{2}=0,0.05$. Typically, $\sigma=40$.

All evolutions show the emergence of a non-axisymmetric instability. The time at which the instability kicks in depends on the type and amplitude of the perturbation; but even the lowest mass unperturbed model ($M_{\rm{Sc}}^{(1)}$) exhibits non-axisymmetric features at a sufficiently long time scale ($t\sim{\cal O}(10^{4})$). The instability generically triggers a larger ejection of angular momentum than mass, reshaping the toroidal energy distribution into a spherical one. This suggests that the asymptotic end state of the cloud evolution is either a spherical (non-spinning) scalar BS or even, merely, ejected debris carrying all angular momentum and energy.

As an illustration, Fig. 2 exhibits snapshots of the equatorial plane evolution of the energy density, $\rho_{E}$ (left panels) and angular momentum density, $\rho_{J}$ (middle-left panels), for the unperturbed scalar initial data with mass $M_{\rm{Sc}}^{(2)}$ 111The energy and angular momentum densities are the ones defined from the Komar integrals, as in Cunha et al. (2017); Sanchis-Gual et al. (2019).. Initially, the collapse preserves axial symmetry.

Around $t\sim 4000$, however, the non-axisymmetric instability is visible, producing a fragmentation event: the star splits into a roughly symmetric orbiting binary. The binary is, nonetheless, bound and recollapses to a deformed spinning star, around $t\sim 6500$. This star breaks into two asymmetric pieces, which again recollapse into a spheroidal star with angular momentum. Around $t\sim 10000$, this residual, still evolving, star has $(M,J)$= (0.49, 0.16), evaluated up to $r=30$, and an oscillation frequency $\omega\sim 0.96$. For this $\omega$, the FF static scalar BS has $(M,J)=(0.45,0)$. Thus, this (or a neighbour) static scalar BS appears to be asymptotically approached, after the remaining $J$ is shed away.

Now consider the formation of a vector SBS. The construction of initial data is more complex due to the Gauss constraint Zilhão et al. (2015); Di Giovanni et al. (2018). After a 3+1 splitting of $\mathcal{A}_{\mu}$, the key variables are the scalar and 3-vector potentials together with the electric field. The first of these admits a solution almost identical to (2), but the others are more involved - Appendix A. These initial data can again be perturbed. We have considered a perturbation analogous to the first type considered in the scalar case; the perturbation amplitudes studied were $A_{1}=A_{2}=0,0.05$. Initial data describing a Proca cloud with three different values of global data were used: $M_{\rm{Pc}}^{(1)}\sim 0.46\sim J_{\rm{Pc}}^{(1)}$, $M_{\rm{Pc}}^{(2)}\sim 0.56\sim J_{\rm{Pc}}^{(2)}$ and $M_{\rm{Pc}}^{(3)}\sim 0.77\sim J_{\rm{Pc}}^{(3)}$.

The unperturbed models evolutions are instability-free during the simulations, lasting up to $t\sim 10^{4}$. This is illustrated by the 3rd and 4th columns in Fig. 2 which show snapshots of the time evolution of the unperturbed Proca cloud $M_{\rm{Pc}}^{(2)}$. The gravitational collapse ejects part of the mass and angular momentum, which shows the gravitational cooling mechanism at play. At $t\sim 10^{4}$ the star has $(M,J)$= (0.35, 0.38), evaluated up to $r=60$, and $\omega\sim 0.985$. For this $\omega$, the FF vector SBS has $(M,J)=(0.301,0.303)$. Thus, this (or a neighbour) vector SBS appears to be asymptotically approached. For the perturbed initial Proca clouds, on the other hand, the energy density oscillates strongly. Nonetheless, no sudden loss of angular momentum is observed, which suggests the endpoint is still a vector SBS.

Evolution of equilibrium SBSs. The dichotomy observed in the formation scenario can be further assessed by considering the dynamics of SBSs obtained as equilibrium solutions of the corresponding Einstein-matter system. A perturbative stability analysis of these SBSs, such as the ones in Gleiser and Watkins (1989); Lee and Pang (1989); Brito et al. (2016) for the spherical case, seems challenging. Thus, we resort to non-linear numerical evolutions of the Einstein-matter system, analogue to the ones in the formation scenario, but now starting with the equilibrium solutions as initial data. This generalises the evolutions in Cunha et al. (2017) for non-spinning BSs.

First consider the scalar SBSs. Fig. 3 shows the time evolution of model $2_{S}$. Up to $t\sim 1000$ the star remains essentially undisturbed; then, following the development of a non-axisymmetric perturbation, see upper panels, the star pinches off into two fragments. The resulting binary is gravitationally bound and collapses into a BH at $t\sim 1200$. This is diagnosed by both the appearance of an apparent horizon (AH), whose mass is shown in the main panel, and the vanishing of the lapse function, $\alpha$, as seen in the inset of Fig. 3. A similar evolution is observed for model $1_{S}$. This confirms that scalar SBSs, even in the FF, are prone to a non-axisymmetric instability. Unlike the formation scenario, here the instability leads to a complete gravitational collapse, likely due to the more compact initial data.

The behaviour of the vector SBSs is distinct. FF solutions, such as models $1_{P}$ and $2_{P}$ show no sign of instability, in the absence of large perturbations. They neither disperse away nor collapse to a BH up to $t\sim 4000$, time at which the drift in the Proca field energy and angular momentum for model $1_{P}$ is 2.0% and 2.2% respectively, whereas for model $2_{P}$, the drift is less than 1%. We further tested the dynamical robustness of vector SBSs by perturbing models $1_{P}$ and $2_{P}$ and by considering some excited states, such as model $3_{P}$-$5_{P}$. Fig. 4 exhibits the time evolution of two examples: $(i)$ model $1_{P}$ with a perturbation of the sort considered in the formation scenario for the vector case, and with a sufficiently large amplitude to visibly distort the star (see first panel), and $(ii)$ the excited model $3_{P}$. In the first case, the perturbation, albeit large enough to deform the morphology of the star away from its spheroidal shape, is dissipated away, and the star recovers its shape. In the second case, the excited state structure of the star is manifest in the composite, Saturn-like, structure of its energy distribution Herdeiro et al. (2016). After $t\sim 1000$, the star abruptly loses energy and angular momentum, until $t\sim 3000$ when it asymptotically tends to a new equilibrium configuration. This new configuration has no nodes and it is close to model $2_{P}$. Thus, the star migrates from the excited family to the FF, where it settles, advocating the stability of the latter. Excited models $4_{P}$ and $5_{P}$, on the other hand, collapse to a BH.

Interpretation and further remarks. The contrasting dynamical properties of the scalar/vector SBSs break the phenomenological (qualitative) degeneracy observed between these two types of BSs in the spherical case. It is tempting to attribute this contrast to the different morphology of these stars, as exhibited in Fig. 1. This interpretation is partly supported by the analogy with dynamical instabilities in differentially rotating relativistic (neutron) stars Paschalidis and Stergioulas (2017). In that case the existence of a toroidal shape has been suggested to be a necessary condition for the development of non-axisymmetric corrotational instabilities Saijo et al. (2003); Watts et al. (2005). In fact, the pinching instabilities and fragmentation exhibited in the scalar models above are reminiscent of evolutions of unstable toroidal fluid stars Zink et al. (2006, 2007), but also of other corrotational instabilities observed in toroidal systems such as the Papaloizou-Pringle instability in accretions disks Kiuchi et al. (2011). Preliminary results, moreover, indicate vector SBSs with $m=2$ (which are toroidal) are also unstable. The analogy with corrotational instabilities in relativistic fluid stars will be further explored elsewhere.

The instability of scalar SBSs may explain the inability to find them as endpoints in the evolution of orbiting binaries of (non-spinning) BSs Bezares et al. (2017); Palenzuela et al. (2017). By the same token, however, vector SBSs should form in the equivalent vector scenarios. This suggests revisiting the work in Sanchis-Gual et al. (2019) using constraint-abiding initial data. A related question pertains the impact of matter self-interactions in the dynamics reported herein.

Acknowledgements. We thank Carlos Palenzuela, Darío Núñez, Juan Carlos Degollado, and Sergio Gimeno-Soler, for useful discussions and valuable comments. This work has been supported by the Spanish MINECO (grant AYA2015-66899-C2-1-P), by the Generalitat Valenciana (ACIF/2015/216), by the Fundação para a Ciência e a Tecnologia (FCT) projects PTDC/FIS-OUT/28407/2017, UID/MAT/04106/2019 (CIDMA) and UID/FIS/00099/2013 (CENTRA), by national funds (OE), through FCT, I.P., in the scope of the framework contract foreseen in the numbers 4, 5 and 6 of the article 23, of the Decree-Law 57/2016, of August 29, changed by Law 57/2017, of July 19. This work has further been supported by the European Union’s Horizon 2020 research and innovation (RISE) programmes H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904 and H2020-MSCA-RISE-2017 Grant No. FunFiCO-777740. We would like to acknowledge networking support by the COST Action GWverse CA16104. MZ acknowledges financial support provided by FCT/Portugal through the IF programme, grant IF/00729/2015. PC acknowledges the Ramon y Cajal funding (RYC-2015-19074) supporting his research. Computations have been performed at the Servei d’Informàtica de la Universitat de València, on the “Baltasar Sete-Sois” cluster at IST, and at MareNostrum (via PRACE Tier-0 Grant No. 2016163948).

Appendix A. Initial data and numerical evolutions. To perform numerical relativity evolutions, one must construct an initial configuration for the spacetime and matter fields that satisfies all the Einstein constraint equations. We solve numerically the Hamiltonian and momentum constraints using the extended conformally flatness condition (XCFC) approximation Cordero-Carrion et al. (2009). For fast rotating neutron stars XCFC has been shown to provide accurate space-times within $1\%$ of full GR. In the following we sketch the main steps for our application.

In the $3+1$ spacetime decomposition the metric can be expressed as

[TABLE]

where $\alpha$ is the lapse function, and $\beta^{i}$ is the shift vector. We apply a conformal decomposition of the spatial 3-metric

[TABLE]

where $\Psi$ is the conformal factor. In XCFC, $\widetilde{\gamma}_{ij}$ is approximated by the flat 3-metric, that simplifies notably the constraint equations. The procedure we follow to solve the constraints is the following.

Make an ansatz for the matter fields (which we will detail below for the scalar and vector case). 2. 2.

Set the geometry, as a first step, to be Minkowski spacetime: $\alpha=\Psi=1$ and $\beta^{i}=0$. Then evaluate the matter source terms in the Hamiltonian and momentum constraints. 3. 3.

Having obtained the matter source terms, solve the constraint equations, obtaining an improved (non-Minkowski) value for $\alpha$, $\Psi$ and $\beta^{i}$. 4. 4.

With the updated values for the metric quantities, evaluate again the matter fields and matter source terms. Then, repeat the procedure iteratively.

At each iteration, we evaluate the L2-norm of the difference between the updated and the previous step values. If it is lower than a certain tolerance ($10^{-8}$), we consider it has converged to a constraint-satisfying solution.

Let us now consider in detail the ansatz for both the scalar and Proca field, and the explicit form of the matter source terms.

.1 Scalar case

The evolution formulation adopted is the one described in Okawa et al. (2014), which involves the scalar field $\phi$ and the conjugated momenta $\Pi$ defined as

[TABLE]

We specify the “shape” of the scalar cloud as

[TABLE]

where $Y_{11}(\theta,\varphi)=\sin{\theta}\,e^{i\varphi}$ is the $\ell=m=1$ spherical harmonic and $R(r)=A\,r\,e^{-\frac{r^{2}}{\sigma^{2}}}$. At $t=0$,

[TABLE]

where we use that $\beta^{\varphi}$ is the only non-zero component of the shift vector, a consequence of the axisymmetry invariance of the energy-momentum tensor. Note that the definition of $R(r)$ includes a factor $r$ that cancels out the $1/r^{2}$ factors appearing in the matter source terms, due to the spherical spatial coordinates, and yields a well-defined behaviour in the limit $r\rightarrow 0$.

.2 Proca case

The evolution formulation adopted is the one described in Zilhão et al. (2015), which involves the scalar potential $\mathcal{X}_{\phi}$, the vector potential $\mathcal{X}_{i}$ and the “electric” and “magnetic” fields $E^{i}$ and $B^{i}=\epsilon^{ijk}D_{j}\mathcal{X}_{k}$, where $D_{i}$ is the spatial covariant derivative and $\epsilon^{ijk}$ is the anti-symmetric Levi-Civita symbol. In the Proca case, besides solving the Hamiltonian and momentum constraints, it is necessary to solve the Gauss constraint which reads

[TABLE]

We assume $E^{i}$ is conservative; thus it is the gradient of a potential $V$, $E^{i}=-\nabla^{i}V$. In this way, the Gauss constraint becomes

[TABLE]

Now, the ansatz for the scalar potential is chosen to be

[TABLE]

where $R(r)=A\,r^{2}e^{-\frac{r^{2}}{\sigma^{2}}}$, and we solve equation (9) to obtain $V$ and, thus, $E^{i}$. Due to the fact that the spherical harmonics $Y_{lm}(\theta,\varphi)$ are eigenfunctions of the operator $r^{2}\nabla^{2}$ with eigenvalues $-\ell(\ell+1)$, Eq. (9) can be reduced to an ordinary differential equation which depends only on the radial part and reads, at $t=0$,

[TABLE]

Similarly to the scalar case, the $r^{2}$ factor in $R(r)$ is needed to have a well-defined behaviour in the limit $r\rightarrow 0$. We solve Eq. (11) analytically specifying boundary conditions $\lim_{r\to 0}V(r)=0$ and $\lim_{r\to\infty}V(r)=0$ to find

[TABLE]

From the potential $V$ we can evaluate $E^{i}$ but we still need to specify a shape for the vector potential $\mathcal{X}_{r}$. The XCFC formalism assumes maximal slicing, $K=0$, where $K$ is the trace of the extrinsic curvature. Therefore, the evolution equation for the scalar potential takes the form

[TABLE]

where we have considered that the only non-zero component of the shift vector is $\beta^{\varphi}$, as in the scalar case. Using Eq. (10) we get

[TABLE]

which has the same form as (8) and can be solved similarly, yielding

[TABLE]

Since $\mathcal{X}^{i}$ can be written as the gradient of a potential function, it is an irrotational field. Therefore the “magnetic” field for these initial data is zero by construction.

As a summary, the $t=0$ initial data for the fields read

[TABLE]

where the $r$ and $\theta$ dependence is:

[TABLE]

Appendix B. Code Assessment. For the numerical evolutions we employ the codes introduced in Refs. Cunha et al. (2017); Sanchis-Gual et al. (2019), which make use of the EinsteinToolkit infrastructure Loffler et al. (2012); EinsteinToolkit ; Zilhão and Löffler (2013) with the Carpet package Schnetter et al. (2004); Carpet for mesh-refinement capabilities and AHFinderDirect Thornburg (2004, 1996) for finding apparent horizons. The evolution of the spacetime metric is handled either by Lean Witek and Zilhão - originally presented in Sperhake (2007) for vacuum spacetimes, and now also distributed within the EinsteinToolkit infrastructure - or McLachlan Brown et al. (2009).

The numerical code we have used to evolve the rotating Proca stars has been assessed in Zilhão et al. (2015) for a real Proca field and in Sanchis-Gual et al. (2019) for the complex case. In this section we will discuss the convergence analysis of the rotating Proca star stable solution. We performed the evolution of model $1_{P}$ with three different resolutions, namely with high resolution, corresponding to $\{(96,48,24,12)$, $(1.6,0.8,0.4,0.2)\}$, medium resolution grid structure of $\{(96,48,24,12)$, $(3.2,1.6,0.8,0.4)\}$ and low resolution $\{(96,48,24,12)$, $(6.4,3.2,1.6,0.8)\}$. The first set of numbers indicates the spatial domain of each level and the second set indicates the resolution. In Fig. 5 we show the time evolution of the angular momentum for the different resolutions. The lower the resolution, the larger is the drift at $t=1000$: 12.3% for low resolution, 2.3% for medium resolution and 0.6% for the high resolution. The L1-norm can be rescaled to second order convergence as shown in the inset of Fig 5. For further details concerning the performance of the codes employed, please see Cunha et al. (2017); Witek and Zilhão for the scalar code and Sanchis-Gual et al. (2019) for the Proca one.