A Well-Posed UV Completion for Simulating Scalar Galileons

The Galileon scalar field theory is a prototypical example of an effective field theory that exhibits the Vainshtein screening mechanism, which is incorporated into many extensions to Einstein gravity. The Galileon describes the helicity zero mode of gravitational radiation, the presence of which has significant implications for predictions of gravitational waves from orbiting objects, and for tests of gravity sensitive to additional polarizations. Because of the derivative nature of their interactions, Galileons are superficially not well-posed as effective field theories. Although this property is properly understood merely as an artifact of the effective field theory truncation, and is not theoretically worrisome, at the practical level it nevertheless renders numerical simulation highly problematic. Notwithstanding, previous numerical approaches have successfully evolved the system for reasonable initial data by slowly turning on the interactions. We present here two alternative approaches to improving numerical stability in Galileon numerical simulations. One of these is a minor modification of previous approaches, which introduces a low pass filter that amounts to imposing a UV cutoff together with a relaxation method of turning on interactions. The second approach amounts to constructing a (numerical) UV completion for which the dynamics of the high momentum modes is under control, and for which it is unnecessary to slowly turn on nonlinear interactions. We show that numerical simulations of the UV theory successfully reproduce the correct Galileon dynamics at low energies, consistent with the low-pass filter method and with previous numerical simulations.


I. INTRODUCTION
In the current era of high precision cosmology and gravitational wave physics, there is a significant interest in understanding how this precise data can be used to test our theories of gravity and low energy particle physics, and to search for potential new physics beyond the standard models. In order to do this we need to have a detailed understanding of the predictions of extensions to standard theories. Effective field theories are a natural tool to describe these corrections, and in the last decade many interesting effective field theories have been developed that extend Einstein gravity while successfully reproducing its successful predictions. A large class of theoretical models that have been developed as alternative descriptions of dark energy and late time acceleration are those that incorporate the Vainshtein screening mechanism [1] (see [2] for a review). This mechanism is built into massive theories of gravity [3][4][5][6][7][8][9][10] and allows them to be made consistent with solar system tests of gravity by screening the would-be fifth forces that are propagated by the additional helicity-zero mode of the massive graviton. While the complete theory of massive gravity is quite complex [10], these essential features are captured by a simplified scalar field theory, the Galileon [11], which incorporates the nonlinear interactions that are responsible for the screening mechanism [9,12]. * gerhardinger1@kenyon.edu giblinj@kenyon.edu a.tolley@imperial.ac.uk trodden@physics.upenn.edu The nonlinear nature of the Vainshtein screening mechanism means that it is difficult to describe analytically beyond very special configurations, and the traditional post-Newtonian or post-Minkowskian formalisms fail to adequately describe the essential features. Given this, numerical progress is crucial, but is unfortunately hampered by the fact that the Galileon effective theory is not wellposed, and despite admitting second order equations of motion it does have regimes in which the equations are no longer hyperbolic. However, such regions arise only when the interactions are taken to be large, for which it is unclear that the effective field theory is under control. Thus, while not a fundamental theoretical problem, successful numerical approaches need to incorporate a mechanism through which to avoid these dangerous regions. For example, in [13] the cubic Galileon was successfully simulated by slowly turning on the non-linear interactions to avoid instabilities.
In this paper we propose a different scheme which replaces the original Galileon theory with a well-posed (numerical) UV completion by means of the introduction of auxiliary higher spin fields. This approach is similar in spirit to that proposed in [14,15], based on the Müller-Israel-Stewart formulation [16][17][18][19] 1 . The key difference is that here the additional spin 1 and spin 2 fields will 1 An alternative numerical scheme is to construct the solution perturbatively in the EFT corrections as in [20,21]. Issues with secular growth of such a perturbative expansion can potentially be resummed in the manner proposed in [22]. We do not consider these approaches as the Vainshtein screening region is necessarily non-perturbative in the leading EFT derivative interactions. arXiv:2205.05697v2 [hep-th] 15 May 2022 be given propagating (hyperbolic) equations in a manner which is closely motivated by the massive spin 2 origin of the Galileon. The cubic Galileon arises consistently as the leading terms in the low energy effective theory of our proposed UV completion, and so we anticipate that a successful numerical treatment will correctly reproduce the dynamics of the Galileon at long wavelengths. We stress that the UV completion here is a numerical one, since no Lorentz invariant local and unitary UV completion of the Galileon is known, and there are now strong arguments suggesting that one does not exist [23]. 2 Since our goal here is to render the low energy theory numerically well defined, we are not constrained by the need to find a local Lorentz invariant action, and thus we propose only a local UV extension of the equations of motion which reproduces the Galileon at low energies. Our approach should render the theory well-behaved when simulating any type of physics, but we will focus mostly on solutions relevant to radiation generated by orbiting binary objects, since this is where the most immediately interesting numerical applications are likely to be.

II. CUBIC GALILEON
We begin by looking at the action for the cubic Galileon (see e.g. [27,28]) where T is the trace of the stress-energy tensor for the matter content. Note that we use a non-standard choice of normalization and coupling that is consistent with how the Galileon degree of freedom arises as the helicity-zero mode in theories of gravity where the graviton is effectively massive. This action yields a classical equation of motion, which makes manifest the Galileon symmetry π → π + c + v µ x µ , and in which the nonlinearity is parameterized by the strength of the coupling, 1/3Λ 3 . In the case where the source is spherically-symmetric and time-independent, with an associated source mass M , there exist static, analytic solutions to eq. (2), 2 These arguments hinge crucially on locality and Lorentz invariance. If the UV theory is mildly non-local then there may be no problem [24]. In addition fractons are a Lorentz violating realization of Galileons [25,26] where E ≡ ∂π/∂r, and where we have defined the Vainshtein radius, The Vainshtein radius sets the distance from the center of mass of the source at which the nonlinear interactions of the Galileon become important. A key feature of the Vainshtein screening mechanism is that this distance is astrophysically large, meaning that most dynamical systems, such as binary pulsars, lie well inside their own Vainshtein radius. For this reason, all linear or perturbative approaches fail to describe the physics of these systems.
This system was studied numerically in [13], where it was shown that eq. (1) could be simulated using numerical tools [29]. This work showed that the scalar power radiated by this system followed anticipated scaling relationships. While this was an important proof of concept, the numerical challenges of simulating this system for more realistic hierarchies fundamentally arise from the physical system itself. When setting up numerical simulations, it is the normal practice to chose numerical parameters to place the physics of greatest interest 'well inside' the box. However, in situations where there are many physical scales, this becomes a more difficult problem. Nevertheless, if the system is formally well-posed, then any UV dynamics of the system remains decoupled from the IR physics of interest. When the system is not manifestly well-posed (as happens, for example, when including higher-derivative operators that inevitably arise from quantum corrections within all EFTs (for a detailed discussion see [30])), then this decoupling of the UV and IR modes is no longer guaranteed. Although effective field theorists have well-established analytic techniques for handling this behavior, it can cause serious problems for numerical implementations. In particular, in our system, if these modes become populated, they run the risk of violating the assumptions of the effective field theory and are analytically unstable.
The question, then, is whether we can regulate these higher-frequency modes using either numerical techniques to dampen them or by finding a physical UV-completion. In this work we compare these two techniques by developing examples of both, and studying whether they are stable, and provide solutions that are consistent with a full-numerical solution to eq. (2).

A. The UV-completion
In standard effective field theories, higher dimension operators can be understood as arising from integrating out high energy (UV) physics, either via tree level or loop level effects. For example, when there exists a 'tree level' UV completion, this means that it is possible to find an action for a well defined classical UV theory for which explicitly solving the equations of motion for the heavy fields in terms of the light fields as a derivative expansion, and substituting back in the action will result in the action for the desired low energy effective field theory. Since the would-be UV completion is valid at arbitrary high energy scales, we would expect it to be well posed. Indeed if the UV completion is Lorentz invariant, we would expect the characteristics of the UV theory to match the Lorentz lightcone, which is to say that the front velocity of propagating modes should be luminal.
In practice we are rarely lucky enough to know the UV theory and in many cases it may be possible to argue that one does not exist, at least satisfying familiar principles. In the particular case of the massless Galileon [31] or massive/weakly-broken Galileon [23,32], there are now well established arguments from positivity bounds that appear to rule out a standard local Lorentz invariant UV completion. It should be stressed, however, that there are implicit assumptions in these arguments which are not required of a UV completion (the UV completion may for example be mildly non-local [24]), and so this does not rule out the Galileon playing a role as an interesting effective field theory. In particular Lorentz violating Galileons emerge in the context of fractons [25,26]. They also seem to play a special role in scattering amplitude methods both for Lorentz invariant theories [33,34] and non-relativistic theories [35,36].
In the present context, our goal is not to find a UV completion satisfying all the principles of unitarity and analyticity, but rather the more modest goal of a completion with high energy behavior that is numerically more stable than that of the initial system (2). Given this, we do not require an action, and at the price of a mild breaking of Lorentz invariance can introduce friction terms to tame unphysical modes. Our proposed method is motivated by how the Galileon arises in massive gravity theories as the helicity zero mode of a massive spin 2 field. In particular given a spin-2 field H µν , the helicity-zero part of it is encoded in H µν ∼ ∂ µ ∂ ν π. Indeed, in massive theories of gravity, this enters explicitly via a dynamical gauge transformation 10,37,38]. With this in mind, we introduce an auxiliary massive spin-1 field A µ and an auxiliary massive spin-2 field H µν that satisfy damped hyperbolic sourced equations. The problematic derivative terms in the Galileon equation of motion are replaced by interactions build out of algebraic functions of the massive spin-2 field. Thus, the UV theory is defined by The presence of the friction terms, parametrized by τ −1 , ensures that the homogenous spin-1 and spin-2 mode solutions of (6) and (7) decay in a time of order τ . The sources on the RHS of (6) and (7) are introduced to ensure that the particular solutions asmptote at low energies k, ω M to and Assuming the approximate validity of (8) and (9), then it is simple to see that (5) reduces to (2), which ensures a faithful UV extension. It is apparent that the UV completion (5),(6), (7) has conventional second order equations of motion with characteristics at high energy determined by the Minkowski lightcone. While not a guarantee of stability of the system, this removes the particular problems associated with the derivative interactions present in the Galileon equations of motion (2). This comes at the cost of replacing the original single field system with a system of 15 dynamical fields. Crucially though, the additional degrees of freedom, even if initially excited, decay away over a time scale τ . There is an alternative way to write the UV completion that makes its connection with the IR theory more transparent. Assuming the homogenous modes of H µν and A µ are set to zero initially, then we may solve for them directly via where D ret (x, y) is the retarded Green's function satisfying with solution is the conventional retarded Green's function for a Klein-Gordon field of mass squared Combining these relations we have where vanishes for x 0 − y 0 < 0, given that the poles lie in the upper half k 0 plane. Substituting (13) in (5) yields a causal integro-differential equation for a single degree of freedom π. At low energies |k| M 2 it is apparent that which shows that the leading term in the EFT expansion reproduces the original system (2).

III. NUMERICAL SIMULATIONS
From a practical standpoint, the full system described by eq. (2) is difficult to study numerically due to the facts that: (1) the system is highly non-linear and derivatively coupled, (2) the system has a number of relevant scales that need to be simultaneously resolved and (3) the effective metric for perturbations can become singular [39]. Nonetheless, it has been shown that the full system can be simulated numerically, with results that are consistent with analytic estimates.
To work with the system numerically, we define dimensionless variables, x µ = x µ prr /2, π pr = π r/M s to rewrite (2) as is a dimensionless parameter that sets the size of the nonlinear terms and J pr is the source. Throughout our work here, we will focus on the system described in [13] (which is a numerical implementation of the binary system studied in [27,40]) and, as we have commented, will also focus on the dynamics of an orbiting binary system. To parameterize this system, we generally use two dimensionless quantities: β ≡r/r v , which relates the diameter of the orbit (roughly the size of the source) to the Vainshtein radius, and α ≡ Ωr, which parameterizes the rotational speed of the system. Roughly speaking, β sets the overall mass of the binary system and α sets the strength of the non-linear effects. With these, Kepler's Law tells us that which fully constrains the system. In practical terms, this sets For our fiducial model here, we take α = 0.2, β = 0.05, which leads to a value of κ ≈ 1.70 × 10 5 .
We expect [13,27,28,40] the system (a Galileon with a cubic interaction) to emit radiation in the quadrupole, with power given by which we can express in a dimensionless way as The analytic expression (20) is computed using the outgoing power for the perturbations of the field about a spherically symmetric background which accounts for the Vainshtein screening due to the monopole. This power is obtained from integrating t π 0r , where t π µν is the stressenergy tensor for the Galileon perturbations (see, e.g. [13]) Provided the power is computed by integrating over a sphere much larger than the Vainshtein radius, this should provide a reliable estimate of the nonlinear power.
In practice, we calculate the power by defining a sphere of radius r * = r, where = 22.5, which is somewhat larger than r v = 20r but less than half the size of the box, L/2 = 50r. Unlike the analysis of [13], we choose to evaluate the radial flux on a set of points defined by the HEALPIX 3 standard. The values of π,π and ∂ r π are calculated at all points over this sphere, even though they are not grid-points, by doing a tri-linear interpolation. Using this process allows us to (i) have assurances that the points are approximately equallyweighted when integrating over the sphere, and (ii) use efficient methods, provided by Healpy, to decompose the fields onto spherical harmonics.
We use different software for each of the sections below; however all are based on GABE [29]-a verified numerical tool for studying scalar fields. While the numerical methods (and hardware) will vary from case to case, we will use the same fiducial physical system and numerical parameters, such as box-size, L = 2.5r v = 50r and number of points along each side, N = 384. In each of the simulations there will be a buffer of points-we generally take this buffer to the 6 grid-points nearest to any boundaryaround the boundary in which the field will not evolve according to the eq. (2), but rather will evolve according to outward-going wave boundary conditions. For the π field, using an assumption that the non-linear terms are negligible at the boundary, this meanṡ While this assumption works very well for massless, Klein-Gordon scalar fields, it remains one of the greatest challenges to successfully simulating our systems.
In each simulation, we take the source to be two, rotating gaussian sources, J = A e −( r pr where r pr ± (t) = (x pr ± cos (Ω pr t pr ) , y pr ± sin (Ω pr t pr ) , z pr ) and the constant is set such so that the total mass of the system is M s = d 3 x ρ = d 3 T . For our fiducial model, we will take L = 50r, N = 384 3 and will use a 'standard' timestep dt fid = β −1r /6400 ≈ 0.003125r.

A. An active low-pass filter
To recover and go beyond the analysis of [13], we begin by attempting to simulate the fully-nonlinear system (2), using just one degree of freedom.
In this system, terms that involve products of secondderivatives are a particular numerical challenge. In these cases, the accuracy of finite-differencing stencils (finiteapproximations to calculating derivatives) is one of the main roadblocks for accurately evolving the field. Accuracy can be gained by increasing the number of nearbyneighbors used to calculate these derivatives; however, there is a substantial run-time cost to that strategy. Therefore, rather than using finite-derivative stencils, we use a spectral method, in which we take a Fourier transform of the scalar field, π, and its time-derivative, π, at each step, and we then calculate the first-and second-derivatives of π andπ in momentum space before performing a set of inverse-Fourier transforms to recover the configuration-space derivatives. This process gives excellent approximations to the derivatives of the field away from the boundary. Luckily, we only need to evolve the boundaries using eq. (23) and we employ 2nd-order, inward finite-differencing stencils to calculate ∂ r π in that region.
This process can be computationally expensive on a CPU, so we employ the GPU-accelerated version of GABE. This version, written in CUDA, maintains all the same structures of the original software, but is written such that the field evolution occurs on a GPU. This acceleration is particularly useful for taking Fourier transforms and is ideally suited for our task. In addition the GPU-accelerated version of GABE also uses a 4th order Runge-Kutta integrator, which, in principle allows us to use slightly larger time-steps as compared to GABE.
However, the greatest benefit to using the GPU accelerated version of GABE is the ability to quickly apply a low-pass filter on the field (or on its derivatives) and thereby to actively remove any high-frequency modes. In practice, we found that applying a low-pass filter, toπ at the end of every Runge-Kutta step, as well as to the Fourier-transform of mixed-spatial derivatives, gave excellent stability without the need to apply additional filters. This filter is designed to cut-off power in modes larger than the one-dimensional Nyquist Frequency to k 1DN = dk N/2 = πN/L, where dk = 2π/L is the standard unit for discrete Fourier Transforms, .
In order to achieve stability, however, one needs to employ slow 'turn-on' strategies like those used in [13]. In this scheme, pr π pr +f 3 (t pr )κ ( pr π pr ) 2 − (∂ pr µ ∂ pr ν π pr ) 2 = f 1 (t pr )J pr (27) where f 1 and f 3 are window functions that start at zero 'ramp-up' to unity, where the choice of smoothing parameter, 0.015 is chosen as a reasonable numerical parameter and 5.25 = 350 × 0.015. Figure 1 shows the quadrupole power emitted in this system as a function of time. The ratio of the late-time quadrupole power in the full system compared to a Klein-Gordon system is 1.81, which is consistent with the values seen in [13].
In Figure 1 (and in following figures) for simplicity we plot the quantityP which is related to the analytic expression (20) for the power radiated bȳ which is useful when comparing multiple methods.

B. Full Auxiliary Field Method
The second approach is to define auxiliary fields, A µ ≡ ∂ µ π and H µν ≡ (∂ µ A µ + ∂ ν A µ ) /2, for which the classical equations of motion describing their interactions are (5),(6), (7). When converting these equations to program units, only a single window function, f 1 (t pr ) is now needed, where M pr = Mr/2.
These massive auxiliary fields cannot use the same massless outgoing wave boundary conditions that we described earlier for the π field. Instead, we enforce the constraint equations for A µ and H µν given in (8) and (9) when the waves reach the buffer (defined as N/64 where N is the size of the box). Using these relaxed constraints, we have been able to achieve numerical stability regardless of when the source 'turned-on'.
In addition to the fiducial tests, we also make a single comparison to a larger, N 3 = 512 3 simulation for Mr = 10. This run will be important as a comparison where we keep the grid-spacing, dx = L/N , constant therefore moving the boundary away from the source without changing the range of high-frequency modes in our system. This run is of particularly importance in diagnosing the limitations of Auxiliary fields as a numerical scheme.
We can compare the results of this system to the full system in a couple of different ways. As a first comparison, we look at the profile of the π field along a line in the equatorial plane of the binary system (take here to be the x-axis). Figure 2 shows excellent agreement between the active low-pass filter and the Auxiliary field methods, particularly at N 3 = 384 3 , as well as agreement with the larger, N 3 = 512 3 simulation.
Next, we can look to the multipole power radiated in the system by the π-field. Figure 3 shows the period- averaged power in the quadrupole for a fairly low, Mr ≈ 0.8, approximately Klein-Gordon simulation as well as the largest mass, Mr = 10, that reached equilibrium. Figure 3 also shows the parametric dependence of the final, quadrupolar power versus Mr. The progression from near-Klein-Gordon to approaching the full, nonlinear system occurs as M transitions from a small to large number compared to one. Below we comment on the limitations of our numerical system to go to higher values of M , however, we anticipate that the trend shown on the right panel of Figure 3 would continue until it matches Figure 1.
In addition to considering the period-averaged power, we can also look for consistency in the power spectra of the π field and its derivative. Figure 4 compares the dimensionless power spectra of the π field, as well as its time-derivative,π. These plots show exceptional consistency between our different methods as a mode-by-mode comparison. This figure also shows the effect of the non- linear terms on the system; the Klein-Gordon (or near Klein-Gordon) simulations have significantly more power on smaller scales which is suppressed as the nonlinear terms become important.
One of the issues we encountered while simulating this model was that the code would crash as we increased the mass of the Auxiliary fields, M . For our simulations, long-term stability became intractable around Mr ≈ 10. For the specific borderline case of Mr = 10, our fiducial model was able to achieve stability for many orbits of the system; however, after some time high-frequency modes are excited and the code becomes unstable. This instability does not seem to arise from a problem with the dynamics of the system, rather, it emerges as a consequence of our boundary conditions. In the boundary, we calculate the derivatives of the Auxiliary fields assuming that the constraints are satisfied and eq. (23). This is a good approximation if (1) we are sufficiently far away from the source such that the π field is Klein-Gordon and (2) the constraints are satisfied exactly. For values of Mr > 5 we seem to violate these assumptions. To demonstrate, we look at our marginal, Mr = 10, case, and test whether the instability is a consequence of numerical instability (by reducing the time-step) or a result of the boundary conditions (by keeping dx the same, but increasing resolution to send the boundary further away from the source). Figure 5 shows that the simulations are not stabilized by increasing time resolution (which would indicate that we're not numerically resolving the problem well); however, the system remains stable for much longer if the boundary is moved away from the source.

C. Restricted Auxiliary Field Method
In addition to the above described UV completion, we also numerically explore a partial UV completion which is obtained from the system (5),(6),(7) by taking the scaling limit M → ∞ for fixedτ In this limit, the equation of motion for the π field remains the same, however those for the additional fields can be reduced to second order equations of motion for the ten auxiliary fields, H µν , given by or more explicitlÿ This restricted system is similar in spirit to the approach taken in [14,15] based on the Müller-Israel-Stewart formulation [16][17][18][19] which has recently been successfully applied to effective field theories of gravity in [41] (for related work on cubic Horndeski theories see [42]). In this approximation, as with the fully UV complete system, we employed outward-going boundary conditions on the π field. In contrast to the full system of Auxilliary fields, however, we did not need to enforce any boundary conditions on the H µν fields since these restricted Auxiliary fields are not propagating degrees of freedom and neither the equations of motion for π nor for H µν depend on derivatives of H µν . Given this, we only need to define H µν in the bulk.
To compare it to the first system, we simulate this system using numerical parameters comparable to the largest stable value of Mr = 8.22-calculatingτ from eq. (35). Figure 6 shows a comparison of the periodaveraged quadruple power.

IV. DISCUSSION
Effective field theories inevitably involve derivative interactions, the effects of which can have important and interesting implications in a number of settings, particularly in gravitational physics and cosmology. While it is well-understood how to analytically deal with the subtleties of solving the resulting equations of motions, significant problems can arise in numerical implementations. This fact has seriously hampered progress in understanding the detailed predictions of large classes of theories that have received much recent attention.
In this paper we have developed, compared, and contrasted three ways of dealing with this problem in numerical implementations of such theories. The first approach is to employ a low-pass filter to tame the UV modes. The second approach is to construct an example of a "UV-completion" of the equations of motion, involving auxiliary fields that constitute new propagating degrees of freedom. The effect of these fields is to render the full system of equations formally well-posed (the system is hyperbolic for all degrees of freedom, and the characteristic speed is unity for all modes), but also to ensure that the IR behavior lies in the same universality class as the original set of equations. The third approach is a restricted "UV-completion", also using auxiliary fields but without introducing new propagating degrees of freedom.
The key point here is that we posit equations of motion that, while remaining second order, now involve damping terms to again tame the UV behavior. Explicitly, we have simulated an orbiting two body system, and determined the power spectrum of scalar radiation of relevance for example to binary pulsars in common examples of modified theories of gravity. We have demonstrated that for the same initial data, all three methods reproduce the same long wavelength physics with the expected errors of the numerical simulations. In the case of the low pass filter, both the source and interactions need to be turned on slowly in order to maintain numerical stability. On the other hand, both of the UV completions are found to be under better control, allowing the interactions to be turned on at the initial timestep. These results parallel those of [43] and more recently [44,45], which consider UV completions of theories with kinetic screening along the lines of [46][47][48].
One remaining technical issue that prevents us from treating large hierarchies of scale (large M ) is that our treatment of the boundary conditions for the massive degrees of freedom is in tension with the damping of the bulk degrees of freedom. This problem arises because of a known issue with imposing boundary conditions in real space for massive fields (see, for example, [49]). A better treatment of boundary conditions should remove this issue.
Our hope is that the techniques described in this paper will be of direct use to those wishing to simulate generic effective field theories, including known difficult examples such as Galileons, massive gravity, and the effects of higher-curvature corrections in gravity.