  Direct detection of dark energy or modified gravity may finally be within
reach due to ultrasensitive instrumentation such as atom interferometry capable
of detecting incredibly small scale accelerations. Forecasts, constraints and
measurement bounds can now too perhaps be estimated from accurate numerical
simulations of the fifth force and its Laplacian field at solar system scales.
The cubic Galileon gravity scalar field model (CGG), which derives from the DGP
braneworld model, describes modified gravity incorporating a Vainshtein
screening mechanism. The nonlinear derivative interactions in the CGG equation
suppress the field near regions of high density, thereby restoring general
relativity (GR) while far from such regions, field enhancement is comparable to
GR and the equation is dominated by a linear term. This feature of the
governing PDE poses some numerical challenges for computation of the scalar
potential, force and Laplacian fields even under stationary conditions. Here we
present a numerical method based on finite differences for solution of the
static CGG scalar field for a 2D axisymmetric Sun-Earth system and a 3D
Cartesian Sun-Earth-Moon system. The method relies on gradient descent of an
integrated residual based on the normal attractive branch of the CGG equation.
The algorithm is shown to be stable, accurate and rapidly convergent toward the
global minimum state. We hope this numerical study, which can easily be
extended to include smaller bodies such as detection satellites, will prove
useful to future measurement of modified gravity force fields at solar system
scales.
