Simulating a numerical UV Completion of Quartic Galileons

The Galileon theory is a prototypical effective field theory that incorporates the Vainshtein screening mechanism--a feature that arises in some extensions of General Relativity, such as massive gravity. The Vainshtein effect requires that the theory contain higher order derivative interactions, which results in Galileons, and theories like them, failing to be technically well-posed. While this is not a fundamental issue when the theory is correctly treated as an effective field theory, it nevertheless poses significant practical problems when numerically simulating this model. These problems can be tamed using a number of different approaches: introducing an active low-pass filter and/or constructing a UV completion at the level of the equations of motion, which controls the high momentum modes. These methods have been tested on cubic Galileon interactions, and have been shown to reproduce the correct low-energy behavior. Here we show how the numerical UV-completion method can be applied to quartic Galileon interactions, and present the first simulations of the quartic Galileon model using this technique. We demonstrate that our approach can probe physics in the regime of the effective field theory in which the quartic term dominates, while successfully reproducing the known results for cubic interactions.


I. INTRODUCTION
While General Relativity (GR) has made myriad successful predictions at many different scales, a number of deep puzzles posed by contemporary cosmology have led to great interest in the question of whether the theory might admit robust modifications.The effective field theory formalism has provided a particularly useful way of exploring and categorizing such modifications, while ensuring that GR is recovered in the local regimes where it has been so stunningly successful.Modifications of gravity and general models of dark energy and cosmic acceleration can be characterized by the mechanisms through which they reproduce local scale gravity-i.e. which screening mechanism they employ.Of particular interest in this paper are those theories that use the Vainshtein screening mechanism [1] to suppress deviations from GR by screening out the fifth-force on local scales [2][3][4][5][6][7][8][9][10][11].There exists a range of complicated theories which exhibit the Vainshtein mechanism, but the Galileon [12], a model with a single real scalar field, is the simplest example that encapsulates all the essential features, and is known to arise naturally in the decoupling limit of massive gravity theories [5][6][7][8][9][13][14][15].
Galileons include higher order derivative interaction terms which give rise to the Vainshtein screening mechanism [5][6][7] yet they have second order equations of motion.Away from extremely special configurations, such as spherically symmetric, static sources, Galileons are difficult to describe analytically, precisely because of these nonlinear interactions.Hence numerical investigations are critical for understanding this model.Although such studies are made easier by the second order nature of the equations of motion, the effective theory of the Galileon nevertheless remains not well-posed; it has regimes in which the equations of motion are non-hyperbolic.These regimes arise when the nonlinear interactions are very large and the effective field theory is no longer necessarily under control.While this does not impact the theoretical underpinnings of the model-it is entirely well-behaved when properly interpreted as an effective field theory-in order to successfully describe the dynamics of the Galileon system, numerical simulations need to avoid or mitigate these unstable regions.
The cubic Galileon model was considered in [16] and it was shown that by turning the nonlinear interactions on slowly it was possible to mitigate any potential numerical instabilities, and obtain a power spectrum which matched the general scaling expectations of analytic estimates [35,36,40].This was significantly improved upon in [17] where we considered the numerical evolution of the cubic Galileon using three distinct methods.The first amounted to low-pass filter combined with the slow turn off method of [16].The other two methods involved replacing the Galileon theory with a system of auxiliary higher spin fields that control the high frequency modes-a well-posed (numerical) UV completion model.Motivated by the origins of the Galileon as a massive spin 2 field, these additional spin fields propagate via either hyperbolic or parabolic equations of motion.The parabolic model is similar to systems proposed in [18,19] that have been successfully applied to simulating evolution within scalar-Gauss-Bonnet gravity [20,21], k -essence [22][23][24][25], Horndeski gravity [26], and other extensions to General Relativity [27][28][29].This technique is predicated on the idea that the theory of interest is a truncation of some larger theory that is well-behaved, and the problems arise from the truncation [20].These approaches, often referred to as fixing-the-equations, are based on the Müller-Israel-Stewart formulation [30][31][32][33], in which the additional degrees of freedom obey their own wave equations.
In our work on the cubic Galileon [17], we demonstrated that the UV completion model asymptotically approaches the original cubic Galileon theory in the low energy limit, and argued that our numerical treatment correctly reproduces the dynamics of the Galileon in the IR regime.It is worth emphasizing that the system provides a purely numerical completion, as there is no known Lorentz invariant local and unitary UV completion of the Galileon, and indeed there are suggestions that such a theory does not exist [34].We also compared the three different ways of resolving the physics in the cubic Galileon model with numerical integration techniques.All three models reproduced the same long wavelength physical processes up to expected numerical errors.
In this paper we take the next important step to understanding the Vainshtein screening by including both cubic and quartic Galileon interactions.This is by no means straightforward since it is known that despite being moderately successful for the cubic Galileon [16,35], analytic attempts to describe the power radiated from a rotating system completely fail in the case of the quartic Galileon [36].This occurs because the large multipoles are not sufficiently suppressed.For this reason it has remained unclear if the Vainshtein mechanism is even active for a time-dependent system when the quartic Galileon is active.In the present paper, we extend our successful discussion of the cubic Galileon using a well-posed UV completion of the equations of motion to account for quartic interactions.In a companion paper we similarly extend the low-pass filter method to the quartic case [37].

II. QUARTIC GALILEON
The Galileon is a real scalar field, π, which satisfies the symmetry π → π + b µ x µ + c, for constant parameters b µ and c.The relevant action for the Galileon including both cubic and quartic interactions is [11,12,36]) is where Λ 3 and Λ 4 control the strengths of the nonlinear couplings in the cubic and quartic terms respectively and T is the trace of the stress energy tensor for the matter content of the theory.The slightly unusual choice of normalization is a reflection of how π emerges as the helicity-zero mode in massive gravity theories where it naturally couples to the trace of the stress energy momentum tensor.This gives rise to a classical equation of motion The most important type of nonlinearity is controlled by the relative magnitudes of 1/Λ 3 and 1/Λ 4 .For a spherically symmetric, time-independent source with T = ρ(r) the Galileon system is well studied [36,38,39], and the equation of motion reduces to where E(r) ≡ ∂π/∂r.We can integrate this to obtain a solution of the form where is the mass contained within a finite radius when the density is approximated as a function of r only.Well outside the source we may take M s (r) = M s constant, without considering higher moments of ρ.This solution reveals two important distances associated with the two scales: Λ 3 and Λ 4 .We define the Vainshtein radii of the source as the radial locations at which the different nonlinear interaction terms become important, given by where r * ,3 is the radius at which the cubic term is approximately the same size as the kinetic term, and r * ,4 is that at which the quartic term is approximately the same size at the cubic term [36].The Vainshtein radii define the distances from the source at which the nonlinear interaction terms become important.It is worth noting that for astrophysical sources, and cosmologically-interesting cutoffs, the Vainshtein radii are extremely large distances.
For more complicated sources, the Galileon system has been studied both analytically and numerically (see e.g.[16,35,36,[40][41][42][43][44][45].Numerical challenges arise when simulating realistic sources in this model -see for instance [19,29].Well-posed field theories describe physics at IR energy scales that are decoupled from the UV modes of the theory.However truncated EFT expansions are generically not well-posed.Although not a fundamental issue (for example there is no ambiguity in inferring low energy S-matrix elements) this can create a significant problem for the numerical evolution of the classical system since the UV and IR modes might be coupled to one another [46].Furthermore when the interactions are higher derivative, as is generic, the evolution may transfer energy from the IR into the UV modes leading to unstable solutions.The question surrounding the well-posedness problem, then, is if it is possible to manage the UV modes without affecting the physics occurring in the IR region.
In our previous paper [17], we introduced and verified two classes of techniques which can accomplish this: using a low-pass filter to dampen high frequency modes, and introducing additional fields to implement a UV completion at the level of the equations of motion.Simulations of the cubic Galileon involving each of these techniques produced results that were both stable in evolution and reproduced solutions consistent with analytical expectations.In what follows we extend the second approach to the Galileon system with both cubic and quartic interactions.

A. The numerical UV Completion
Following [17] our goal is to introduce a numerical UV completion which is designed to control the high energy behavior and allow for more stable numerical evolution than the original Galileon system.Note that this is not a true UV completion, i.e we do not require an action, rather this is a device to regulate the numerical evolution.With this in mind, we are allowed to introduce friction terms that mildly break Lorentz invariance to help prevent the instabilities of high k modes [47].Using the same technique as in the cubic Galileon case [17], inspired by how the Galileon can arise as the helicity zero mode of a spin 2 field within massive gravity theories, we introduce an auxiliary massive spin 1 A µ and a massive spin 2 H µν field that satisfy damped harmonic oscillator equations [8,9,48,49].These new fields replace the troublesome higher order derivative terms in the Galileon equation of motion with functions of the H µν fields; thereby allowing us to control the derivative interaction terms that are the origin of the instability problems.We first define auxiliary fields via and We then promote the auxiliary fields to have their own dynamics and, at the same time, trade the higher derivative Galileon interactions for the lower derivative fields so that the UV theory is described by the equations of motion These equations have been constructed such that the friction term ensures that the new propagating degrees of freedom in A µ and H µν decay on a time order of τ .As in the cubic case, we have now turned a single propagating field into a total of fifteen fields, all of which obey a form of the wave equation.Importantly, the fourteen additional degrees of freedom should decay away on a time scale of order τ .The only difference between this quartic Galileon system and the cubic Galileon system considered in [17] is an additional term in the equation of motion for the π field.
The terms on the right-hand side of ( 11) and ( 12) dictate that solutions will asymptote at low energies (k, ω ≪ M ) to the definitions of A µ and H µν respectively, given in (8).Additionally, the definitions of the auxiliary fields are constraints on the boundary conditions of our numerical system.In the low energy limit, it is clear to see that (10) reduces exactly to (2), the equation of motion for the Galileon that we started with.The addition of these fields does not guarantee stability, however it does eliminate any issues that arise from the derivative interaction terms.For another way to write out this numerical UV completion-one which makes the connection to the IR regime more apparent-see [17].

III. NUMERICAL SIMULATIONS
The Galileon effective field theory poses several problems for numerical simulations due to its nonlinear nature, coupled derivative interaction terms, and the many scales of interest present in the system 2 .Having shown in a previous work that the full cubic system can be solved numerically in a way that produces results which agree with analytic expectations, we use here the most stable of our methods to produce the simulations of the quartic Galileon model.These numerical simulations are conducted using GABE, a verified numerical program for simulating scalar fields [51]. 3  To be concrete, we focus on the case when the source is comprised of two rotating Gaussian-shaped mass-energy distributions.This is a challenging example that is closely related to physically relevant sources and also allows for easy comparisons with our previous work with the cubic system as well as with other numerical simulations [16] and analytic expectations [35,40].We parameterize the system through two dimensionless quantities: α ≡ Ωr, which describes the rotational speed of the system, and β ≡ r/r * ,3 , which relates the diameter of the sources to the cubic Vainshtein radius as defined by Eq. (7).To fully constrain the system, we can use Kepler's Law (a reasonable approximation for small velocities) where M s is the total mass of the system, to describe the cubic nonlinear coupling term via We choose β = 0.05 and α = 0.2 as our fiducial model, which translates into a coupling strength κ cubic ≈ 1.7 × 10 5 .

A. Dimensionless Units
For the simulation to be independent of the physical parameters of this model, we rescale both our scalar field and the spatial coordinates to yield dimensionless variables (denoted with the subscript pr), defined by 2 Other work has focused on singularities in the effective metric of perturbations as a cause for numerical issues [50]. 3http://cosmo.kenyon.edu/gabe.html where r is the distance between the centers of the rotating Gaussians, and the rescaling of the spatial derivative clearly follows from this latter definition.These allow us to rewrite the equation of motion given in (2) as where we regularize the sources using J pr J pr = A e −(⃗ r pr + (t)/σpr) where ⃗ r pr ± (t) = (x pr ± cos (Ω pr t pr ) , y pr ± sin (Ω pr t pr ) , z pr ) and the constant are chosen to ensure that the total mass of the system is given by M s = d 3 x ρ = d 3 x T , as expected.Additionally, we define the dimensionless quantities which control the size of the nonlinear terms and the strength of the quartic term relative to the cubic term, respectively.In the limit ξ → 0, this system will reduce to the cubic Galileon model.Additionally, the auxiliary spin fields are rescaled according to We investigate the quartic Galileon system, with these sources, in a variety of regimes; we also define a fiducial case for which the physical and numerical parameters are: box size L = 2.5r * ,3 = 50r, number of points along the box N = 384, and time step size dt = β −1 r/6400 = 0.003125r.

B. Analytic Expectations
The analytic solution for a single Gaussian source is determined by the polynomial equation where we define E(r) ≡ ∂π pr /∂r pr .It is evident that as ξ → 0 we recover the cubic system.For ξ > 0 there will be 3 regions depending on which of the three terms Klein-Gordon, cubic, and quartic dominates.We can solve this analytically and determine the relative importance of each of the three terms for a given value of ξ. Figure 1 depicts the contribution of each of these terms to the equation of motion for the π field for a given value of ξ = 0.6 and a single Gaussian source.It is clear that for ξ = 0.6, the binary system is well inside the region in which the quartic Galileon dominates, namely the Vainshtein radius r * ,4 .For this reason we expect our simulations to be substantively different than those for the pure cubic Galileon ξ = 0.

C. Boundary Conditions
In each simulation, we use outgoing boundary conditions, so that the physical processes inside the box are the only sources to the fields.To isolate the system, we construct a buffer of points inside the lattice at each boundary where the field is no longer propagating according to its equation of motion given in (10) but rather evolves using outgoing boundary conditions.For the scalar π field-or any interaction-less and massless field-this looks like Although this works well for the π field, problems arise for the massive auxiliary fields.As in the cubic case, we minimize those issues by enforcing a parabolic damped constraint equation, given by where C is some parameter that determines the extent to which the auxiliary fields are damped.For our simulations, we set this constant by the decay scale for the auxiliary fields, C = M 2 pr τ , a parameter to which the simulations are largely insensitive below a value of C = τ × 10 4 .These are the same 'damped' versions of the constraint equations used in the cubic system [17], where we determined that directly enforcing the constraint equations on the boundary introduces large discontinuities that caused backreaction into the box.In the low energy limit, we expect that A µ ∼ ∂ µ π and H µν ∼ ∂ µ ∂ ν π, which would result in a zero time derivative-and no backreaction.
As energy increases beyond the low energy limit, these approximations no longer hold, and the time derivatives of the auxiliary fields are no longer trivial.In that case, our goal is to minimize the effect that the discontinuities have on the system, and to enforce that the auxiliary fields behave like derivatives of the π fields, so that discontinuities do not disrupt the physical processes occurring at low energies.The discontinuities appear because of the contribution of the cubic terms to the equation of motion at the edge of the box.At the boundary, we assume that the cubic terms vanish, but as Figure 1 depicts, there is still around a 40 % contribution at the boundary.At first glance, one viable solution to this appears to be to make the box bigger in order to allow the cubic contribution terms to vanish more completely before reaching the boundary.However this significantly reduces the resolution of the sources in the middle of the box.Hence, we enforce damped constraint equations to mitigate backreaction from the discontinuities and to maximize the spatial resolution of the source.

D. Power Calculations
To establish whether the Vainshtein mechanism is active for a time-dependent rotating binary source with the quartic interactions present we can calculate the power emitted per multipole.We calculate this in the same way as described in [17], by setting up a sphere of points, evaluating π, π, and ∂ r π at those points, and then performing a tri-linear interpolation to project the field onto the spherical harmonics.The sphere is set up using the HEALPIX4 standard, defined by a radius of r = 22.5r which is larger than the Vainshtein radius, r v = 20r, but smaller than half the size of the box, L half = 25r.This method provides us with a sphere of evenly spaced points (not necessarily at the lattice point locations) which are used to decompose the field onto spherical harmonics.
In [36], the authors analytically found that the power in higher multipoles are suppressed and that the dominant multipole was the quadrupole.This was confirmed by the numerical simulations of [16,17].By contrast in [35] it was found that including a quartic interaction term greatly alters the behavior of perturbations around a spherically symmetric, time dependent source rendering the approximate perturbative treatment inadequate for calculating the radiated power.Nevertheless some insight can be gained from considering the form of the analytic radiated power which should be valid when there is a large hierarchy between the scales involved which consequently allows for the linear treatment of perturbations from a static, spherically symmetric background solution.In this situation, the radiated power for a given mode depends on m and ℓ, derived in Appendix A of [16], and is given by where we have restricted ourselves to circular orbits in the θ = π/2 plane with equal mass objects, For ℓ = 0, we are constrained to m = 0, which causes the power to vanish, as we expect for the monopole.For ℓ = 1, the two allowed values for m also lead to zero power: ℓ = 1, m = 0 is zero because of the m contribution at leading order, ℓ = 1, m = 1 is zero because the (1 + (−1) m ) term is also zero.For ℓ = 2, the first non vanishing power term, m = 2 is the sole contributor because m must be >0, and it must be even, so 1 + (−1) m ̸ = 0.For ℓ = 3, the constraints on m leave only the m = 2 case as a possibility for non-zero power.However, this is also zero because Y ℓ=3,m=2 (π/2, 0) vanishes.Extrapolating, we expect that the only non-vanishing power will arise in modes that have even ℓ (ℓ = 2n where n ∈ N) with contributions from m modes that are also even and greater than zero.

E. Full Auxiliary Field Method
In our previous work, we described two distinct numerical UV completions for simulating the cubic Galileon model, both of which reproduced the long wavelength behavior expected by analytic analysis.While each of these had their own advantages and shortcomings, we have chosen to employ the 'Full Auxiliary Field' method for simulating the quartic system because it allows nonlinear interactions to turn on sooner and generally requires less spatial resolution.Our chosen UV completion is in program units □ pr H pr µν − where M pr = M r/2 and f 1 (t) is a window function given by that ramps up from zero to unity.Using this system and the relaxed numerical constraints as boundary conditions for the auxiliary fields, we have been able to achieve numerically stable simulations of the quartic Galileon model.
An important check on the numerics is to first ensure that the quartic simulations reproduce the cubic system previously simulated by taking the limit ξ → 0 for a given mass value.Figure 2 depicts the asymptotic approach of the π field profile from the quartic system toward the field profile of the cubic system, taken along a line in the equatorial plane of the system.We argued in [17] that the Full Auxiliary field method reproduces the dynamics of the cubic system in a way that matches analytic expectations-here we validate the current method by comparing it to previous cubic work [17].
In addition, we have investigated the dynamics of this quartic model by examining the multipole power in the system radiated by the π field.Figure 3 depicts the period-averaged power contained in each multipole for a single value of ξ, the parameter controlling the relative strength of the quartic term.The multipoles diminish in power as the multipole number increases and only select even poles (ℓ = 2n, m = 2n where n ∈ N) have appreciable nonzero power, with any power in odd multipoles arising from numerical errors.Although these results match the analytic expectations discussed in Section III D we stress that for comparable mass binary sources for which the orbit lies inside the quartic region, the analytic approach breaks down.Our result is hence nontrivial and match the results of the low-pass filter method in [37].
While the period-averaged power contained in each multipole highlights the differences between the quartic and cubic dominated systems, we can also study the effects of varying M for a fixed quartic strength, ξ. Figure 4 shows how the final quadrupole power depends on the M r parameter.At M r = 0, the only source for π is H 2 µν , but H µν itself is unsourced.Therefore, any oscillations in H µν will decay away and, in the long time limit, π will approximately satisfy the Klein-Gordon equation.For large M , the power converges to a value that is insensitive to increasing M further, indicating that the UV physics is decoupled from the IR regime.Comparing Figure 4 to the analogous one for the cubic only system in [17], we see that the final quadrupole power depends on the M r parameter in a similar way, and we note that the same asymptotic behavior is exhibited with higher resolution runs as was noticed in the previous paper.In other words, in the M r > 7 regime, the simulations approximate the M r → ∞, decoupled limit.We also investigated the M − ξ phase space by examining the effect of increasing the ξ parameter for a single M value on the period averaged quadrupole power.Figure 5 portrays the final quadrupole power for a fixed mass M value.Here, increasing the ξ parameter from 0.0 to 0.7 encodes the progression from the cubic system to non-zero quartic contribution regime, with the quartic term dominating for ξ ≥ 0.6.This highlights how the power is increasing as we turn on the quartic interaction.
In the Full Auxiliary Field Method, issues with the boundary conditions meant that we were only able to simulate the system up to a certain M r parameter, for a low value of ξ, the same limiting value that we found in our previous work with the cubic system [17].Instead of M r = 10, the cutoff for the cubic Galileon, here the cut off is closer to M r ≈ 7.As the quartic term contributes to the equation of motion more and more (i.e. as ξ increases), the cutoff on M r decreases, because the cubic contribution at the edge of the box also grows.For higher mass runs, the fiducial model in the cubic system could stably evolve for several orbits of the source, but eventually power seeped into the higher frequency modes and caused the simulation to crash.At the boundary, we calculate derivatives of the auxiliary fields while assuming that the constraints are satisfied -a good approximation given that the constraints are satisfied exactly and we are far enough away from the source that the π field is Klein-Gordon.Similarly to the cubic case, it seems that these assumptions are violated for M r ≥ 5, because the cubic terms are non-zero at the edge of the box.We have tested the idea that this is a result of the boundary conditions and not the model itself by increasing the spatial resolution of our marginal case M r = 7 to effectively move the boundary further from the source, and found that our simulations evolved longer.In principle, one should be able to probe beyond the range presented here by increasing spatial resolution and the number of points used in the finite derivative equations, although this comes at the cost of computational resources and time (we estimate that on the 20 core machine with 512 GB RAM memory owned rights.Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof.The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof."

FIG. 1 .
FIG.1.Analytic contribution of each of the Klein-Gordon (solid black line), cubic (red dashed), and quartic (blue dot-dashed) terms to the total π field for our fiducial model with α = 0.2, β = 0.05 and ξ = 0.6 illustrating the two distinct Vainshtein screening regimes.

FIG. 2 .FIG. 3 .
FIG. 2. (Left)The π field profile along the x axis after the system has reached stability using a mass value of Mpr = 3 and an increasing nonlinear quartic term, controlled by ξ.The cubic case (ξ = 0) is denoted by an solid black line, the ξ = 0.4 a gray dashed line, and ξ = 0.6 a red dot-dashed line.

FIG. 4 .FIG. 5 .
FIG.4.Late-time quadrupole power emitted by the fiducial system with a constant value of ξ = 0.6 and changing M parameter.The orange dots have 384 3 resolution whereas the blue dots have 512 3 resolution, included so as to demonstrate the asymptotic behavior noticed in the cubic system in[17].