Scalar Radiation with a Quartic Galileon

The class of Galileon scalar fields theories encapsulate the Vainshtein screening mechanism which is characteristic of a large range of infrared modified theories of gravity. Such theories can lead to testable departures from General Relativity through fifth forces and new scalar modes of gravitational radiation. However, the inherent non-linearity of the Vainshtein mechanism has limited analytic attempts to describe Galileon theories with both cubic and quartic interactions. To improve on this, we perform direct numerical simulations of the quartic Galileon model for a rotating binary source and infer the power spectrum of given multipoles. To tame numerical instabilities we utilize a low-pass filter, extending previous work on the cubic Galileon. Our findings show that the multipole expansion is well-defined and under control. Moreover, our results confirm that despite being a non-linear scalar, the dominant Galileon radiation is quadrupole, and we find a new scaling behaviour deep inside the Vainshtein region.

The simplest example of the Vainshtein mechanism is provided by the cubic Galileon [11] which originally emerged as the decoupling limit of the Dvali-Gabadadze-Porrati model [39].This model is sufficiently simple to be reasonably acquiescent to analytic approximations [40,41].In a recent work, three distinct numerical approaches were utilized to describe the scalar radiation emitted by a binary system, such as a binary pulsar [42].One of these approaches built on previous works [43] by directly simulating the cubic Galileon using a lowpass filter to suppress any numerical instabilities.The other approaches utilized an extended system of equations of motion, in effect a well-posed UV completion, designed to reproduce the Galileon theory at low energies.In the present work we will extend the low-pass filter approach to the quartic Galileon, and in a companion work [44] the UV completion method will similarly be utilized for the quartic Galileon theories.What makes the quartic Galileon particularly challenging is that unlike in the cubic case, analytic attempts to describe scalar radiation using the one-body effective approach fail [41].This failure is a reflection of the fact that on linearizing around spherically symmetric backgrounds, the usual centrifugal repulsion which suppresses power in high order multipoles is lost, leading to a superficially divergent power spectrum.This simply indicates the failure of perturbative analytic approaches.It is worth noting that this problem is ameliorated by having a larger cubic Galileon term, and this is mirrored by our numerical simulations which are more stable for large cubic interactions.
In the present work we show that just as for the cubic case, the Vainshtein mechanism is successfully realized for these time-dependent solution in the quartic Galileon models (in the presence of a cubic interaction), something which could not be inferred from any approximate analytic treatment.In common with the pure cubic Galileon, we find that the dominant scalar radiation is quadrupole.However, when the source radiates in the region in which the quartic interaction dominates, the scaling of the radiated power with both frequency and multipole number is distinct from the cubic case demonstrating that the quartic Galileon screening while qualitatively similar to the cubic is quantitatively distinct.The results in this paper closely parallel those in [44] where the quartic Galileon system is treated by embedding it in a 'UV completion' that renders the dynamics well-posed.

The Model
The Galileon is a scalar field theory which exhibits a non-linearly realized 'Galileon' symmetry π → π+v µ x µ [18].Among the class of interactions that are invariant under this symmetry are a finite set of interactions that lead to second order equations of motion.In 3 + 1 dimensions these are the cubic, quartic and quintic Galileon terms.In the present work we shall be concerned with the generic quartic Galileon meaning that we shall include both cubic and quartic interactions and consider a coupling to the trace of the stress energy of matter The precise form of the normalization is characteristic of how the Galileon arises as the helicity-zero mode in massive gravity theories [17,45].If we assume spherical symmetry and a static point-source, we can integrate (2.1) to obtain an ordinary differential equation for E ≡ dπ/dr [35,41,46] We will use a system consisting of two gaussian 'masses' orbiting around their center of mass 3) The locations of each mass, at a given time, are ⃗ r(t) = (x ± r cos (Ω p t) /2, y ± r sin (Ω p t) /2, z) and the constant A is chosen so that the total mass of the system is M s = d 3 x ρ = d 3 x T and r is the diameter of the circular orbit.We constrain the system to be Keplerian, for an orbital frequency Ω p , which is to say we will ignore relativistic corrections to the orbit.The quantity Ω p r is an important quantity that parameterizes the power in scalar radiation from the Klein-Gordon system, (2.5) 3 Numerical Implementation The Galileon system has been studied both analytically and numerically (see e.g.[23, 27, 36-38, 40, 41, 43, 47]).At least two of the numerical implementations use modified versions of GABE using a standard, CPU implementation [43] as well as a GPU implementation [42].
In [43] the authors show that turning on the sources and the non-linear interactions slowly at the beginning of the simulation can avoid the issue of a singluar effective metric [41].This method was made more robust in the GPU implementation, one of the methods described in [42], where efficient Fast Fourier Transforms (FFTs) made it possible to calculate mixed spatial derivatives using spectral methods and cut-off high frequency modes.In the present work our goal is to extend these methods to the quartic Galileon.
All of the simulations presented here use grid sizes of N 3 = 384 3 points where the length of each side is L = 50r.These parameters are chosen as to be comparable to the previous literature.We use a set of dimensionless units to rescale spacetime (we can choose physical units such that M Pl = 1) and the Galileon field, In these dimensionless units, we can re-write the quartic equation of motion, (2.1), in a pictorial form that allows us to separate the non-linear interactions O n from the strength of these interactions κ n and the numerical turn-on functions, f n , where O 2 is the usual Klein-Gordon term and we have two non-linear operators and The source is taken to be where we take σ pr = σ 2/r = 1/3.The cubic interaction is parameterized by the dimensionless quantity whereas the quartic interaction is parameterized by If we only consider the cubic term, there is a clear connection between the size of κ 3 and the Vainstein radius, In the presence of a quartic term, we have two important radii that differentiate between the regions where the quartic term, the cubic term, or the Klein-Gordon term dominates the left hand side of (2.2).While we can set an expectation of where the second of these boundaries is based on (3.9), we can numerically solve (2.2) to find both.
In a particular EFT realization, there might be a clear reason to expect that the two interactions scales are related Λ 4 ∼ Λ 3 , however generically they can be independent.In order to truly probe the Vainshtein screening effect of the quartic Galileon we will chose a value of κ 4 which is large enough for the quartic term to dominate in the region of the source.To use the parameterization of [44], we set where ξ = 0.6 is a fiducial value that realizes our goal.Fig. 1 demonstrates how this choice puts the boundary between the region of quartic dominance and the region of cubic dominate at a value of r/r ≈ 1.67 -this value is substantially far away from the sources which orbit at r/r = 0.5 with a width of σ/r = 1/6.Larger values of ξ should, in principle, be possible however require significantly more computational time due to the need to slowly turn on the interactions.
In the present work we need to be more careful than in previous work [42,43] when turning on the sources and interactions.These choices are not unique, but they allow us to get all the simulations presented here to a stable state, where the code can run until its endtime.We turn on the source using In the fiducial case, Ω p r = 0.2 so the system orbits about every t/r = 30.
and the two interactions using These are smoother and slower functions than are needed for pure cubic simulations and are shown graphically in 2. In addition to the simulations needing longer program time to initialize, we needed to reduce the timestep from the choice made in [42] to ∆t = 3∆x/625 = 0.0048 ∆x.
We use the same CUDA-accelerated version of GABE as described in section III a. of [42].In this numerical scheme, we employ spectral methods to calculate spatial derivatives of π and φ; because the derivatives exist in Fourier space, we apply a low-pass filter, on φ at the end of every timestep.As in [42], the cutoff scale is chosen to be the equivalent 1-dimensional Nyquist frequency, k 1DN = dk N/2 = πN/2, where dk = 2π/L.Again, this is not a unique choice, but one that allows for our simulations to remain stable until late times.
One major challenge that has been apparent in numerical simulations of the Galileon theories is the treatment of out-going wave boundary conditions [42].In the present work we continue to use which is strictly speaking only valid for massless Klein-Gordon fields.Although in principle this should be suitable when applied to an interacting Galileon at large distances, in practice it is hard to perform simulations for large hierarchies, meaning that in practice the radii at which we evaluate the boundary conditions is not significantly larger than the Vainshtein radii.Given this, the interaction terms have not truly switched off and the above boundary conditions can lead to waves reflecting at the boundary of the simulation and potentially generating instabilities that terminate the runs too soon.In the alternative UV completion method considered for the cubic Galileon in [42] and the quartic Galileon [44] a second problem arises which is that it is difficult to implement out-going boundary conditions even for a linear massive scalar.We leave to future work a potential solution to these problems.

Results
To illustrate the efficacy of our numerical scheme, we consider the fiducial model with parameters Ω p r = 0.2, r/r v = 0.05 together with a choice of ξ = 0.6.For these values, our dimensionless couplings are κ 3 = 1.70 × 10 5 and κ 4 = 1.35 × 10 9 .We find that when running the simulations, a large ξ (and hence large quartic Galileon) necessitates a longer turn on time for the simulations, making the runs inefficient.The alternative approach considered in [44] has the advantage that it may be implemented without needing to slowly turn on the interactions.
While there is no analytic solution to the time-dependent problem, we can still look to validate the code, by comparing the time-averaged field profile to that of the (semi-)analytic solution to (2.2).Fig. 3 shows a comparison of rE(r), as calculated by taking the positive x-axis in the simulation to be the r-direction and averaging each point along that axis over the final two periods of the simulation.We show both the entire range of the positive x-axis as well as a close-up near the center.We can see that the solution is nearly identical to the spherically symmetric profile-and much different from the cubic-only profile.Deviations from the solution are expected at small r/r ≲ 1 where the source is not well described as a point source.We can now turn to calculating the power emitted by the system.The outgoing energy flux is given by [43] Note that this is not the non-linear expression for the power, but rather that for perturbations around a spherically symmetric background.At large distances where the background monopole dominates this is expected to be a sufficiently good approximation to the true radiated power.We evaluate the relation (4.1) at a radius where the Klein-Gordon term dominates, r/r = 22.5, which is halfway between the Vainshtein radius and the closest edge of the box.Again, following the procedure first described in [42], we evaluate dπ/dr and π on a set of points on the sphere defined by the HEALPIX2 standard using a tri-linear interpolation.We can then use the efficient Healpy [48] routines to decompose this power onto the spherical harmonics.When we report the power, we further perform a rolling time-average over one orbital period.Fig. 4 shows how the first three moments, namely the monopole, dipole and quadrupole, behave as a function of time.We note that the dipole power is zero to machine precision throughout the simulation.The initially large monopole is an artifact of the way the interactions and source are turned on for which energy is not conserved and a large monopole artificially appears.We find that time-averaged quadrupole power remains the dominant mode, as in Fig. 4 in the presence of the quartic interaction.(a) Power versus multipole.Blue dots denote the cubic, κ4 = 0, Galileon whereas the black dots the fiducial quartic Galileon system.In both cases, Ωp r = 0.2.The dashed red line is the best-fit to the low-ℓ multipoles for the cubic system and the dashed red line is the best-fit for the quartic system.The quadrupole power, P 2 , for our fiducial simulation, Ω p r = 0.2, κ 4 = 0.6 (red, solid) compared to the quadrupole power in a cubic Galileon, κ 4 = 0 simulation (black, dotted).Note the increase in power as the quartic interactions are turned on.Fig. 5 shows how the power in the quadrupole differs in the quartic case from the purely cubic Galileon.Interestingly for our fiducial values the power with the quartic interaction included for the same value of cubic terms is actually larger, although not by orders of magnitude.Another major difference that arises when the quartic term is added is the modification of the spectrum of power as distributed among the modes.This is illustrated in Fig. 6(a).For cubic case the best-fit shows an ℓ-dependence of P ∝ ℓ −6.4756 whereas the quartic case shows a best-fit of P ∝ ℓ −9.414645 as in Fig. 6(a).
Finally, we can now turn our attention to parameterizing how the power depends on Ω p r.To isolate this, we calculate the ratio of the quadrupole power from our simulations to the Klein-Gordon expectation, (2.5), as can be seen in Fig. 6(b).A simple power-law fit of the scaling of the power with frequency gives which varies significantly from the cubic-only case, where [43] where the latter was consistent with analytic approximations [40,41].Unfortunately we do not at present have a semi-analytic understanding of this distinct scaling, however it is worth noting that already for the static solutions the scaling of π with r in the region where the quartic dominates is quite different from the pure cubic theory.

Conclusions
In this work we have been able to successfully simulate a generic quartic Galileon, and use this to determine the power radiated into scalar radiation for a rotating binary source.This is of great interest in models of modified gravity where additional scalar degrees of freedom can couple to the trace of the stress energy with gravitational strength, but are at the same time screened by the Vainshtein mechanism.Understanding the amount of power for a generic binary system such as a binary pulsar is important in order to put observational constraints on such models.The present work is a step in this direction for modified gravity models whose decoupling limit is well described by a quartic Galileon, as is common in both soft and hard massive gravity theories.
First and foremost our results show that the Vainshtein mechanism is fully active in this time-dependent situation, as previously found in simulations of the cubic Galileon [42,43].In particular the time averaged field configuration matches well analytic expectations of the screened solution.This is a non-trivial result since attempts to provide an approximate semianalytic treatment fail [40,41].We confirm that despite the highly non-linear nature of the system, the dominant scalar radiation is quadrupole, with the next most significant mode ℓ = 4 being typically several orders of magnitude smaller.Although these results parallel the cubic case, we find that the quartic Galileon leads to a qualitatively similar but quantitatively different scaling of the power with orbital velocity and multipole number.
Although the particular numerical scheme we have used here is successful, it was necessary to turn each interaction and source on slowly to tame any potential numerical instabilities which may arise due to the fact that the Galileon system is not strictly well-posed.Furthermore the larger the quartic coupling parameter ξ the slower the rate of turn on needed to avoid any instabilities (we were able to simulate ξ ≤ 0.6 with smaller ξ being considerably easier).This unfortunately renders simulating large hierarchies ξ ≫ 1 too costly in time at present.A solution to the turn on problem is to use the UV completion method proposed in [42] which is considered for the quartic case in [44].This latter method avoids the particular instabilities associated with the system of equations not being well-posed.Both the present simulation and that of [44] have difficulty giving a proper treatment of the radial boundary conditions and we suspect that a better treatment of the boundary conditions will improve the stability of typical runs.

Figure 1 .
Figure1.The relative contributions of the different terms on the left hand side of the spherically symmetric equation of motion (2.2).We show the relative contributions of the quartic (blue), cubic (red) and Klein-Gordon (black) contributions for our choice of ξ = 0.6 (solid) and what would happen in the case of a larger, ξ = 0.95 (dashed).In both cases we set the size of κ 3 such that the Vainstein radius is approximately r/r = 20.

Figure 3 .
Figure 3.The profile of rE averaged over two periods at the end of a simulation with Ω p r = 0.2 (black dots).This is compared to the static spherically symmetric solution, (2.2) (red, solid) and a cubic-only static spherically symmetric solution (blue, dotted) and a purely Klein-Gordon static spherically symmetric solution (grey, dashed).The left panel shows the full range of 0 < r < 50r resolved in our simulations and the right panel shows just the range 0 < r < 20r to more clearly show disagreement between the spherically symmetric expectations and the simulation near the source.

Figure 4 .
Figure 4. Power as a function of time for the fiducial, Ω p r = 0.2 and ξ = 0.6 model.The curves show the period-averaged power in the monopole (black), dipole (blue) and quadrupole (red) as a function of time.
Power versus orbital velocity.The dots show the the late-time qudrapole power versus Ωp r for a set of quartic Gaileon simulations.The dotted red line is a best-fit line.

Figure 6 .Figure 5 .
Figure 6.The period-averaged power at the end of a set of simulations with (a) multipoles and (b) varying Ω p r.