Proca in an Expanding Universe

The superradiant growth of massive vector fields in rotating black hole spacetimes has garnered significant attention in recent literature. However, the majority of these studies overlook the influence of a cosmological constant, which likely constitutes the primary energy content of our universe. In this paper, we extend recent research by incorporating a cosmological constant into the Einstein+Proca system and numerically evolving the resulting equations of motion. Utilizing the newly released GRBoondi numerical relativity code, designed specifically for the numerical evolution of (generalized) Proca fields, we discover that parameters causing a growing instability in the $\Lambda = 0$ scenario transition to a decaying state when $\Lambda>0$. This results in a more intriguing phenomenology. These simulations pave the way for future full Einstein+Proca simulations to explore the secular decay of the resultant cloud from gravitational emission.


I. INTRODUCTION
As an effective field theory, General Relativity performs exceptionally well at intermediate scales.However, it faces significant challenges at both extremely small and vast cosmological scales, necessitating extensions to the theory [1].Its problematic reconciliation with quantum mechanics and the inconsistencies observed in cosmology indicate the need for new physics.Two of the most pressing questions in contemporary cosmological research pertain to the seemingly invisible components of the universe, whose existence is inferred from compelling observational data [2][3][4].These unknown components are colloqiually dubbed 'dark matter' and 'dark energy', the latter being responsible for the accelerated expansion of the universe.The former is predicted to exist from observations of galactic rotation curves and gravitational lensing measurements [5][6][7][8][9][10][11][12][13][14].The nature of these components remains elusive, however there are many experiments and theoretical analyses today that are looking to uncover these invisible sectors of our model of the universe (see for example [15][16][17][18][19][20][21] and references therein).
The dark photon in particular is a well-motivated candidate for dark matter and it has several possible production mechanisms, including a misalignment mechanism [22], arising naturally in certain string theories [23], tachyonic couplings to misaligned axions [24], and quantum flucuations during the inflationary epoch [25].Constraints on the couplings of such particles to the Standard Model come from equivalence principle tests, such as of the Eöt-Wash group [26,27] and Lunar Laser Ranging groups [28,29].Near future constraints will come from GW measurements, for example from black hole superradiance [30][31][32][33].The dark photon does not need to coincide with the photon of the Standard Model.Observational data strongly suggest that it must be a completely new particle beyond the Standard Model.This dark photon could exhibit intriguing interactions, including derivative self-couplings [34].This could lead to fascinating cosmological applications beyond just explaining dark matter [35,36].
Regarding dark photon applications as dark matter, superradiance will soon offer stringent tests on the mass of massive vector fields, thereby constraining the allowable parameter space for dark matter candidates.Due to this, several theoretical analyses have been carried out to understand the evolutionary behavior of a massive vector field surrounding black holes [30,32,[37][38][39][40][41].It has been shown that the Proca field exhibits a very rapid instability phase before saturation, afterwhich it experiences a secular decay due to gravitational radiation emission.The instability phase can be very rapid, typically on the order of seconds.On the other hand, the secular decay due to gravitational radiation emission can last from months to several billion years depending on the mass of the vector field.Both the fast and slow gravitational radiation decay can be used to constrain the vector field mass using upcoming gravitational wave (GW) antennas, such as the LISA mission [42,43].
It is thus important to develop accurate models of the evolution of the Proca cloud surrounding black holes.Previous studies have focused on the evolution of the Proca cloud surrounding a Schwarzschild or Kerr black hole, with no influence from the cosmic expansion [30,32,[39][40][41][44][45][46][47][48][49][50][51].This is normally a justifiable approximation since recent cosmological observations have shown that the cosmological constant is quite small [2][3][4].Nonetheless, the effects of a cosmological constant on the growth rate of the Proca cloud is still required to complete the understanding of Proca superradiance in our universe.It is thus the goal of this study to further the understanding of Proca in a Kerr-de Sitter background.We opt to neglect the backreaction of the Proca cloud on the background since the magnitude of the stress-energy tensor of the Proca field is negligible.This is a common assumption for the initial phase of superradiant instability studies since the cloud begins as a perturbation on the spacetime.As long as the Proca field does not constitute a dominant contribution and its interactions remain weak, the backreaction effect is expected to be minimal.Going beyond the initial phase to studies of the resulting gravitational radiation would require full mutual evolution of the matter and metric variables, hence this is beyond the scope of this study and is relegated to future work.
The paper is divided into 4 sections.Section II takes care of outlining the basic theoretical background, including the equations of motion and the particular coordinate system.Section III then discusses the results of this study.Finally, Section IV concludes and discusses future work.An appendix expands on the various formulas utilized in the main text.

II. THEORY
The starting point is the Einstein-Hilbert-Proca theory, where the Einstein-Hilbert action is augmented by the Proca action where g µν is the metric, g its determinant, R is the Ricci scalar, Λ is the cosmological constant, and A µ is the Proca 4-vector.Variation of S[g, A] with respect to these two fields yields the equations of motion where G ρσ is the Einstein tensor and T µν is the stress energy tensor of the Proca field.We assume the backreaction of the Proca field on the spacetime is negligible, as is usually the case, meaning we can set T µν = 0 without significant loss of accuracy.This implies that the field equations reduce to The solution of the gravity sector describing a spinning black hole is called the Kerr-de Sitter (KdS) solution.This solution is in fact a special case of the more general Plebanski-Demianski family of metrics [52,53].We choose to work in the Kerr-Schild form of the solution, which takes the form where g 0,µν is the background de Sitter metric and K µ is a null vector (with respect to both g and g 0 ).In Kerr-Schild coordinates (t, r, θ, ϕ), the de Sitter background metric takes the form where we've defined and H = 2M r ρ 2 .The ∆ r definition becomes important for the analysis of the black hole horizons.The poles of ∆ r correspond to the poles of the metric in Boyer-Lindquist coordinates [54].In a similar fashion, the null vector is defined via Prior to plugging this metric into the numerical solver, we need to understand the basic causal structure of the background spacetime on which the Proca cloud will evolve.This amounts to determining the location of the horizons, which follows by solving the quartic polynomial ∆ r = 0.The existence of the horizons can be determined by analyzing the discriminant of the quartic polynomial, which will be denoted by Q.The positivity of Q guarantees that either all four roots are simultaneously either real or complex.This means we only have to verify one of the roots is real to ensure the other 3 are as well.The discriminant is easily solved for and takes the form Fig. 1 shows a contour plot of Q, with the spin and cosmological constant rescaled by the black hole mass.The region Q > 0 shows the allowed parameters for the background.The shaded region denotes the disallowed region.Beyond the allowed region, one or more horizons will disappear, yielding a naked singularity.Another interesting feature is the existence of a minimum spin for a certain range of values for the cosmological constant.For spacetimes satisfying Λ ≥ 1/9, the black hole is required to possess spin in order for the Kerr-de Sitter black hole to exist.That is to say, for cosmological constant values greater than 1/9, static black holes in the form of Eq. 6 do not exist.
We label the three positive roots as r − , r + , and r Λ , which denote the inner, outer, and cosmological horizons, respectively.The fourth root is negative and corresponds to a 'horizon' inside the singularity at r = 0.In Kerr-Schild coordinates, besides the singularity at r = 0, the poles in the metric disappear, including r Λ .However, a new pole appears at rΛ = 3 Λ .We find that r Λ ≤ rΛ throughout the allowed parameter space.In the limit of small cosmological constant, A peculiar property of the causal structure is the existence of a maximum cosmological constant which permits the existence of a black hole.At the maximum cosmological constant value Λ max , the outer and cosmological horizons merge, see for example Fig. 2. Additionally, dimensionless spin values greater than unity are allowed without producing a naked singularity, as long as the cosmological constant is greater than zero.Past the value of Λ where the horizons merge, no Kerr-de Sitter black hole is possible.
On the other hand, the solution in the Proca sector of Eqs. 4 is unclear, since analytic studies have yet to be performed.However, based on knowledge from the Proca on Kerr solutions, we can expect the cloud to undergo a superradiant instability due to the existence of the black hole horizon.This implies that the energy of the cloud will be of the form where ω is the complex frequency of the Proca field.The factor of 2 comes from the fact that the stress-energy tensor of the Proca field is quadratic in the Proca field itself.Since the frequency is of the form ω = ω r + iω i , the imaginary part implies an exponential evolution superimposed over an oscillatory one.The imaginary frequency is what is solved for in this paper.

A. 3+1 decomposition
Towards a numerical solution of the Proca system, the next step is to decompose the field equations into a form pertinent for numerical computations.We follow the standard procedure, which is to decompose the spacetime via a foliation into a series of three dimensional time-like hypersurfaces.Since the spacetime admits a time-like killing vector, we can choose the Kerr-Schild time coordinate as the function that defines the foliation leaves.Hence, the (r, θ, ϕ) coordinates become coordinates on the hypersurfaces.We thus define our 3 + 1 decomposition as [55] ) where g and γ are the determinants of the full and spatial metrics, respectively, α is the lapse function, β i is the shift vector, γ ij is the spatial metric, and Latin indices range from one to three.In the coordinate system of Eq. 7, these take the form where we've defined Γ = 2M r∆ θ + ρ 2 ΘΛ r .As our numerical solver computes the time evolution using cubic cells, we transform the (r, θ, ϕ) coordinates to a Cartesian-like coordinate system defined by The shift, spatial metric, and all derivatives are then transformed using the resulting Jacobian matrix.See App.A for more details.The last step is to decompose the Proca equations Eq. 5 under the foliation.A standard calculation yields where Z is an auxiliary field introduced to damp violations of the Proca constraint with a tuning parameter κ [32,56,57], and n µ is the time-like normal to the spatial hypersurfaces.

III. NUMERICAL RESULTS
For the numerical computation of the Proca field, we use the recently released GRBoondi software [58].GRBoondi is an offshoot of the GRChombo [59] framework that is catered towards numerical solutions of Generalized Proca theories [34].GRBoondi allows for the numerical evolution of Generalized Proca theories on arbitrary fixed backgrounds and so is an excellent choice for the system under consideration here, even if we concentrate here solely on the mass term.Generalizing our study to the more intriguing case of self-interacting Generalized Proca fields is straightforward with GRBoondi.
GRBoondi solves the Proca equations Eq. 23-26 using an adaptive mesh refinement (AMR) grid with time evolution determined using a Runge-Kutta fourth-order time integrator and fourth-order accurate finite difference stencils for the spatial derivatives.We use a box width of L = 60M with N = 192 grid points across each edge of the computational box.We use 4 refinements levels at a 2 : 1 refinement ratio, resulting in a resolution of the finest level of dx f ine = 0.01953M .To prevent boundary effects from contaminating the simulation, we use Sommerfeld-outgoing radiation boundary conditions, which allows oscillations to exit the simulation region with minimal reflections due to finite-size effects.This is especially important since we introduced an auxiliary field which dampens violations of the constraint equation and the evolution equation for Z is a generalized telegraph equation.This implies that not only are the values of Z damped, but also propagate at the speed of light.Hence, the outgoing radiation boundary conditions are vital for ensuring violations of the constraint equation propagate outside the computational domain.
To understand the effect of a cosmological constant on the dynamical evolution of the superradiant Proca cloud, we perform a various number of simulations with parameters that yield the highest growth rates.We choose three different values of the cosmological constant, Λ = 5 • 10 −6 , 10 −4 , 10 −3 .Higher values of the cosmological constant are more difficult to simulate numerically as the cosmological horizon quickly becomes small.We reserve probing this region of the parameter space to future studies, which will likely entail a new coordinate system.Additionally, we fix the black hole spin to χ = 0.99.We sample the Proca mass at 6 different values, µ = (0.35, 0.4, 0.45, 0.5, 0.6, 0.7).In addition to the main simulations, we also perform a convergence study to ensure our choice of resolution produces accurate data, which we discuss in sec.III A.
For initial data, we take a Gaussian profile with width determined by analytic approximation studies [41], r 0 = 1 M 2 µ .The initial data is then A x = A γ e − r r 0 , where A is some pre-determined amplitude which we take to be A = 0.1, and all other variables are chosen to be zero.
The data from our simulations is available in table I.The cosmological constant has been rescaled to Λ = Λ M 2 , where Λ is the unscaled parameter.Plots of the normalized total energy as a function of time are shown in fig. 3. A curious new feature seems to arise in the case of decaying modes, namely a secondary scale.For example, in the case of µ = 0.6 and Λ = 5 * 10 −6 , the decay rate of the total energy slows at around t = 4000M .Similar features can be seen in the µ = 0.7/Λ = 10 −3 plot as well as µ = 0.6/Λ = 10 −4 .Whether this is a numerical artifact or a real emerging scale is uncertain.Analytic studies will need to be performed to determine the true nature of this emerging scale.Due to the variations in the µ = 0.35/Λ = 10 −3 , such a feature cannot be determined.It should be noted that the rapid oscillations and cessation of decay in the µ = 0.35/Λ = 10 −3 plot is likely purely numerical errors.The energy of the Proca cloud becomes incredibly small, it likely runs in to a precision floor and the simulation becomes numerically meaningless past t ∼ 2000M .

A. Convergence tests
To ensure the reliability and accuracy of our simulations, we conduct two types of convergence tests.These tests are crucial for validating the fidelity of the simulated data and confirming that it accurately reflects the underlying physical system.The first type of test we perform is the grid resolution test.By running our simulations at multiple grid resolutions, we can verify that the results converge to a stable solution as the grid is refined.This process helps us identify and minimize numerical artifacts that might arise from discretization errors.Specifically, we monitor key physical quantities and ensure that their values stabilize with increasing grid resolution, indicating that our results are not dependent on the grid size.The second type of test is the time step refinement test.Here, we run our simulations with progressively smaller time steps to ensure that the temporal evolution of the system is accurately captured.By comparing results obtained with different time steps, we can check for stability and convergence in the time integration scheme.This test is particularly important for dynamic simulations where rapid changes in the system need to be accurately resolved.In order to maintain the Courant-Friedrichs-Lewy condition, the temporal timestep decreases proportionally with the spatial timestep.Thus, the two checks are unified in our resolution tests.Below, we provide detailed insights into the methodologies employed for conducting these two tests and assess the reliability of our simulations based on their outcomes.

Resolution Test
The first convergence test is a resolution test.This ensures that finite-differencing effects have a diminishing effect on the data as the resolution is turned up.In other words, the resulting data should converge to stable values as the resolution is increased.After some point, increasing the resolution further should have little effect on the data.This ensures that the resolution we have chosen for our simulations is satisfactory enough to produce good data.
The results of this convergence test are as follow: we repeat two of the simulations with parameters (µ, Λ) = (0.4, 5 * 10 −6 ) and (0.4, 10 −3 ) with a resolution on the finest level of dx f ine = 0.0167M .Since the Courant-Friedrichs-Lewy condition is automatically adjusted for in GRBoondi, the corresponding time resolution is dt f ine = 0.0033, an increase in resolution from the main data, which utilized dt f ine = 0.0039.It was found that the instability rate for

Analytic Derivatives Test
The second convergence test is a self-check test on the metric derivatives.Since the procedure in transforming the metric variables and their derivatives from Kerr-Schild coordinates to a Cartesian-like coordinate system is involved, we run a check to make sure the numerically computed derivatives converge to the analytically computed ones from Eq. A3-A4.We follow the procedure of [60] for carrying out this convergence test.The test procedure is as follows: where the convergence factor c(t) is defined as and ∆ i denotes the i'th resolution, ϵ ∆i = U (t, x) − U ∆i (t, x), U (t, x) denotes the variable computed using the exact analytic derivatives, and U ∆i (t, x) denotes the variable computed using the numerically computed derivatives at a resolution of ∆ i .The resolutions are typically chosen to be twice the previous resolution.In other words, ∆ i = 2 i ∆ 1 .This implies that for an nth-order finite differencing scheme, Since GRBoondi uses fourth-order finite differencing stencils for the spatial derivatives, we expect c(t) = 16 for each variable.We perform the convergence test for two different resolutions, which differ by a multiple of two per the previous discussion.We find that the minimum convergence factor across all grid variables was 14.35, in fairly good agreement with the expected convergence factor of 16 for a fourth-order finite differencing routine.

IV. CONCLUSION
This study investigates the dynamics of a Proca field surrounding a spinning black hole within an expanding universe, a novel approach that advances our understanding of superradiant vector fields interacting with cosmological constants-a crucial yet overlooked aspect in prior research.Leveraging the advanced capabilities of the recently developed GRBoondi software, we conducted a series of simulations aimed at quantifying the growth rates under diverse scenarios involving different Proca masses and cosmological constants.Furthermore, we conducted a rigorous resolution analysis to validate the accuracy and reliability of our simulation data, revealing robust agreement with expected theoretical predictions.
There are a few limitations of the current study.Firstly, the effects of backreaction were ignored, however this approximation is suitable for this analysis as the energy densities were small and gravitational radiation was not desired.Should gravitational emission from the resulting Proca cloud be desired, full numerical computations will need to be performed.Due to GRBoondi outputting checkpoint files in the same format as GRChombo, using GRChombo to perform full computations with initial data from the simulations here will be straightforward.Additionally, many parts of the parameter space were left unsampled, due to the numerical complexities.Primarily, simulating larger cosmological constants is difficult due to the cosmological horizon rapidly approaching the outer horizon of the black hole.Additionally, simulating lower spins increases the simulation time considerably, hence we only focused on a single value of the black hole spin, χ = 0.99.A more complete analysis would likely require a different coordinate choice that penetrates the cosmological horizon.
Future studies to be performed include theoretical analyses of various coordinate systems in the search of ones that can penetrate the outer and cosmological horizons.Additionally, future studies will turn to full evolution of the Einstein+Proca system to compute the emitted gravitational radiation and apply the results to gravitational observatory forecasting.An intriguing progression would involve extending beyond the Proca field's mass term to incorporate derivative self-interactions inherent in Generalized Proca theories.This expansion would enhance our exploration of the field's dynamics, encompassing interactions that go beyond simple mass considerations and delve into the complexities introduced by derivative couplings within these theories.

Appendix A: Coordinate Transformations
In transforming metric variables from Kerr-Schild to the Cartesian-like coordinate system Eq.20-22, we follow the standard procedure for coordinate transformations.That is β i cart (x, y, z) = Λ i j β j KS (r, θ, ϕ) (A1) where Λ i j = dX i dR j , X = (x, y, z), and R = (r, θ, ϕ).The derivatives of the metric variables are then computed straight forwardly using the following rule for a rank-r contravariant tensor where a tilde repesents the quantity in the base coordinate system and non-tilde in the new coordinate system.For a covariant rank-r tensor, we find a similar rule For example, the spatial metric in the new Cartesian-like coordinate system can be computed as The rest of the derivatives of the metric variables follow similarly.The Jacobian matrix of the transformation Eq. 20-22 is simple to compute.Since at each grid point, we are given the (x, y, z) coordinates, we represent the Jacobian matrix in terms of the (x, y, z) variables.Hence, it is given by where r 2 = x 2 + y 2 + z 2 and ρ 2 = x 2 + y 2 .Due to the division by the two radii, GRBoondi sets a minimum value for the radii of 10 −12 .For more details on the numerical implementation, see [61].

FIG. 1 :
FIG. 1: Plot of the values of the discriminant Q.The shaded region shows the unallowed parameters for the existence of a Kerr-de Sitter black hole.The unshaded region are the allowable parameters.

√ 3 1M 2 9 16 + 3 √ 3 8
Fig.1shows a contour plot of Q, with the spin and cosmological constant rescaled by the black hole mass.The region Q > 0 shows the allowed parameters for the background.The shaded region denotes the disallowed region.Beyond the allowed region, one or more horizons will disappear, yielding a naked singularity.Fig. 1 hence tells us the allowed values of the black hole parameters that can be used in the simulations.An absolute maximum value of the cosmological constant is Λ max = 16 45+26 √ 3 1 M 2 and corresponding absolute maximum value of the black hole spin of a max = 9 16 + 3 √ 3 8 M .This point corresponds to the cusp of the non-shaded region on the upper-right quadrant of Fig. 1.Another interesting feature is the existence of a minimum spin for a certain range of values for the cosmological constant.For spacetimes satisfying Λ ≥ 1/9, the black hole is required to possess spin in order for the Kerr-de Sitter black hole to exist.That is to say, for cosmological constant values greater than 1/9, static black holes in the form of Eq. 6 do not exist.We label the three positive roots as r − , r + , and r Λ , which denote the inner, outer, and cosmological horizons, respectively.The fourth root is negative and corresponds to a 'horizon' inside the singularity at r = 0.In Kerr-Schild coordinates, besides the singularity at r = 0, the poles in the metric disappear, including r Λ .However, a new

FIG. 2 :
FIG.2: Evolution of the cosmological and outer horizons, eventually merging.A value of a = 0.5 is chosen as a representative value.Past the value of Λ where the horizons merge, no Kerr-de Sitter black hole is possible.

FIG. 3 :FIG. 4 :
FIG.3: Growth of the energy of the cloud over simulation time.Energy values are normalized to their value at t = 250M .The energy data is fitted in logarithmic space to a linear function, which captures the exponential characteristics of the instability.The fitting of the data starts at t = 250M .The green line is the numerical data and the blue straight line is the fit function.

TABLE I :
All available simulation data.