THE MASS OF THE MILKY WAY AND M31 USING THE METHOD OF LEAST ACTION

Estimating the total masses of galaxies, our own in particular, is a continuing challenge of precision cosmology. Part of the challenge lies in the unknown extent of the dark matter halos within which they are presumably embedded: while the measurement of galaxy rotation curves from coherent stellar motions allows the mass within the visible radius to be inferred, the total mass of the associated dark matter halos predicted in the standard model of cosmology, whose physical extent is not known, is more difficult to estimate. To probe the total effective gravitational mass, the analysis must include the effect on the peculiar motions of nearby galaxies.

The measurement of total galaxy masses from their relative motions was pioneered by Kahn & Woltjer (1959). Their "timing argument" (TA) method, which assumes purely radial infall, indicated a total mass for the MW+M31 system of about 3 × 10¹² M_☉—a lower bound, since the possibility of transverse motions is excluded. The total mass of the Local Group (LG) can also be computed from the velocity dispersion of its various member galaxies, assuming that it is in virial equilibrium and that its velocity ellipsoid is isotropic: Courteau & van den Bergh (1999), using this method, found an LG mass of (2.3 ± 0.6) × 10¹² M_☉. More recent applications of the TA tend to suggest higher masses. Li & White (2008) confirmed that the TA method used on mock galaxies drawn from the Millennium Run (Springel et al. 2005) systematically underestimates the true mass, and revised the TA method to predict an LG mass of (5.27 ± 0.5) × 10¹² M_☉. Van der Marel et al. (2012) used the full proper motion of M31 to improve the TA method, estimating an LG mass of about 5 × 10¹² M_☉, which is somewhat higher than the combined prediction of (3.17  ±  0.57) × 10¹² M_☉ from a Bayesian combination of estimates from different methods. The TA can also be used with the proper motion of Leo I (Sohn et al. 2013), assuming it is gravitationally bound to the Milky Way (MW), to estimate the mass of the MW alone. Boylan-Kolchin et al. (2013) combine the TA and other methods in estimating a virial mass for the MW of 1.6 × 10¹² M_☉ with a 90% confidence interval of [1.0–2.4] × 10¹² M_☉.

Mass estimation methods using the TA are based on the analysis of single-galaxy interactions with the MW. We show in this paper that the numerical action method (NAM), by taking into account the peculiar motions of a large subset of LG satellites, effectively breaks the mass degeneracy in the two-body TA and identifies separate ranges of likely masses for the two principal actors in the LG. The method avoids the TA assumption that galaxies are gravitationally bound, makes no assumptions about virialization of the LG, and is sensitive to more widely diffused concentrations of dark matter that could remain undetected using other methods. NAM takes as input the cosmological parameters H₀ and Ω₀, and assumes that linear theory correctly describes velocities at early times and that galaxies and their progenitors back in time can be approximated as simple paths representing the center-of-mass motions of their associated dark matter halos.

Earlier papers on NAM (including Peebles 1989, which introduced the method; Peebles 1995; Peebles et al. 2001, 2011; Phelps 2002; Peebles & Tully 2013) developed the method in the context of the LG, focusing on the analysis of individual solutions or small ensembles. In what follows, we explore the behavior of several thousand independent solutions in order to identify the combinations of MW and M31 masses that yield the best fit to the observational constraints. Section 1 introduces our implementation of NAM. In Section 2, we apply it to the LG, and in Section 3 we test NAM predictions in the LG with mock catalogs drawn from the Millennium Run. A brief discussion including an assessment of future prospects is found in Section 4.

2.1. Method and Approximations

Our version of NAM is based on the improved algorithm described in Peebles et al. (2011). We have extended it to include an optional partial canonical transformation of coordinates between the radial distance and the radial velocity (first used in Peebles et al. 2001 and more fully described in Phelps 2002), allowing either the galaxy distance or redshift to be chosen as the boundary condition in the action integration at z = 0, in addition to the two components of the observed galaxy position in the plane of the sky. This effectively doubles the solution space that can be explored. The Appendix gives details of the coordinate transformation.

We model galaxies as constant-density spheres with a MW radius of 100 kpc and all other galaxy radii scaling proportionately with the cube root of their masses. From the start of the computation time at a = 0.1 until a = 1, galaxy paths are interpreted as the mean motions of the coalescing systems of baryonic matter and their associated dark matter halos. As shown in Peebles et al. (2011), the initial velocities are by construction consistent with linear perturbation theory. This keeps the galaxies well separated from each other at early times, and so passages of galaxies through each other's (nonphysical) cutoff radii are rare.

Galaxy flight paths are reconstructed from initial randomized straight line trial orbits by successive iterations in the direction of the first and second derivatives of the action until a stationary point is reached. Solutions are verified by comparing them to the solutions to the equations of motion from the same initial timestep in a leapfrog approximation. With a target of 10⁻¹¹ in the sum of squares of the gradient of the action, deviations between NAM-generated flight paths and the leapfrog approximation are at most a few kpc at the final timestep. In the present work, we follow Peebles & Tully (2013) in excluding several close satellites to the MW, but we use a smaller number of time steps (30, versus 500 in Peebles et al. 2011; Peebles & Tully 2013), which is still sufficient to produce tightly bending half-orbits. This significantly reduces computation time and maximizes the number of solutions generated.

NAM solutions are non-unique owing to the mixed boundary conditions: velocities are constrained at the initial time (see discussion in Appendix A.1.4), while some combination of distances and velocities are fixed at the final time. As a consequence, different choices of initial trial orbits with the same observational constraints will in general yield a different set of galaxy paths, each of which are solutions to the equations of motion. However, as there are more constraints than those required for a solution (it is sufficient to specify either distances or redshifts), a χ² measure of fit is defined for each observable,

$$\begin{equation} \chi ^2 = \sum _{i} \left({{\rm model}_i - {\rm cat}_i \over \sigma _i}\right)^2, \end{equation} \tag{ 1 }$$

from which a best-fit solution can be selected from an ensemble.

Apart from the initial trial orbits, solutions may also be sensitive functions of the input parameters, particularly for galaxies in close proximity to each other. These parameters are fixed in the formal NAM computation even though there is an observational uncertainty associated with each. Since we wish to explore the widest possible range of reasonable physical configurations, we adopt the method, similar to Peebles et al. (2011) and Peebles & Tully (2013), of defining a χ² measure as a sum over all constraints (distance, redshift, angular position, mass, and, where available, velocity transverse to the line of sight), fixing the observational constraints as their given values plus a random error within the observational uncertainty that differs for each trial solution. Our per-particle χ² is defined as follows

$$\begin{eqnarray} \chi ^2_{{\rm tot}} & = & \sum _{i}\left(\chi _{i,d}^2 +\chi _{i,cz}^2 + \chi _{i,\theta }^2 + \chi _{i,\phi }^2\right. \nonumber \\ & &\left. + \chi _{i,{\rm mass}}^2 + \chi _{i,{\rm vtrv}}^2 + \chi _{i,v0}^2\right). \end{eqnarray} \tag{ 2 }$$

We relax the initial trial orbits to a solution to the equations of motion using NAM, and then relax each NAM solution to a minimum in χ². By holding some of the observational constraints relatively fixed, we can repeat this procedure for a large number of NAM solutions to explore an n-dimensional space of solutions whose minimum in χ² gives the most likely values of the desired n constraints. In our case, we are interested in the masses of the principal actors in the LG and so we are looking for a minimum within the two-dimensional space defined by the masses of MW and M31. Further details of our minimization approach, which is fully general and can be applied to any desired combination of observational parameters, are as follows.

For each solution we assign masses to MW and M31 within a range of 0.5 and 6 × 10¹² M_☉ in a Gaussian random distribution centered around the nominal masses given in the catalog. This produces higher quality solutions by focusing less attention on combinations of masses, which experience has shown are less likely to yield good fits to the constraints. Random errors are added to all other observables (distances, redshifts, angular positions, and masses), with standard deviations taken generically to be 10% of the distance for distances, 5 km s⁻¹ for redshifts, 60% of the nominal catalog masses (giving these quantities the widest possible latitude), half a degree for angular positions, and the published uncertainties for proper motions (as summarized in Peebles & Tully 2013; these are available for M31, LMC, M33, IC10, and LeoI). Our relative leniency in angular positions, which we allow owing to the approximate nature of the model, also accounts for the physical possibility that galaxies are offset from the center of mass of their respective dark matter halos. Since the main objective of this study is investigating the masses of MW and M31, and good agreement with M31 constraints is therefore a priority, standard deviations in distances, redshifts, and positions for galaxies beyond 1.5 Mpc are doubled and those for M31 are halved. Standard deviations in MW and M31 masses are reduced by a factor of 20 relative to their initial randomized guesses, essentially fixing them in the χ² relaxation step but allowing them to shift by a small amount if a significantly better solution is found at a slightly different mass. As large initial velocities are permitted in the solutions (see Appendix A.1.4), a contribution to χ² from the magnitude of the initial velocity (the final term in Equation (2) above) is additionally assigned with a standard deviation of 40 km s⁻¹. This effectively suppresses implausibly large initial velocities in the search for the best solutions.

We find that NAM solutions are efficiently generated by generating a solution for the two dominant galaxies first and then adding the other galaxies, one at a time, in descending order of mass. With each NAM solution, we then take each galaxy in turn in the same order, using Powell's method to relax the individual galaxy distance, redshift, and angular position to a minimum in its χ², since derivatives in χ² with respect to these quantities cannot be reliably computed. Discontinuities in χ² may be encountered where a galaxy jumps to a qualitatively different orbit in the descent to the target in the action; these are allowed so long as it leads to an improvement in χ². If after relaxing with Powell's method χ² remains above a given threshold (100 yielded a good balance between time to solution and quality of solution), then we recast the orbit for this galaxy and find another solution. If after 50 attempts χ² is still above the threshold, then we switch from distance to redshift boundary conditions for this galaxy and allow another 25 attempts. If χ² still remains above the threshold, we take the solution corresponding to the lowest χ² found thus far and move to the next galaxy. In practice, between 25% and 50% of the galaxies in any given NAM solution will have been fitted to the velocity boundary condition and the rest to the distance boundary condition. Once χ² has been minimized in this way, we then jointly relax the galaxy masses, also using Powell's method, to minimize χ² still further (we found that relaxing the masses separately from the other quantities gave more rapid convergence to a minimum). This method of relaxation between the various observational quantities typically reduces χ² by an order of magnitude from its initial value; further improvements to this ad hoc procedure are no doubt possible.

We generate four thousand independent solutions using this method, varying the mass of MW and M31 in each solution as described above. Parallelizing the code to run on 12 supercomputing nodes using 2.40 GHz six-core Xeon processors, each solution typically takes 10 or 20 s for catalogs up to a dozen or two particles. From an ensemble of solutions, we then plot a Gaussian-smoothed χ² map against the MW and M31 masses, dividing the map into 24 × 24 equal bins and keeping the best χ² in each bin. The smoothing is desirable because the solution space can feature many local minima, each corresponding to a qualitatively different configuration of orbits.

Our LG catalog, based on Peebles & Tully (2013), is listed in Table 1. The omission of the smallest, tightly bound satellites of MW and M31 facilitates direct comparison to the simulations, which lack comparable range in mass, and maximizes the number of solutions that can be explored. The four actors listed immediately after MW and M31 are intended to approximately model the influence of the significant mass concentrations residing immediately beyond the LG. While inclusion of all LG actors may hold the promise of producing even better dynamical constraints on the mass of MW and M31, it is offset by the greatly increased effort required to produce acceptable solutions.

χ² maps for 4000 solutions in four different scenarios, all with H₀ = 67 and Ω₀ = 0.27, are shown in Figure 1. At the upper left are the contours in χ² generated from a simplified catalog consisting of only MW and M31, to check the consistency of our implementation of NAM with the TA. As expected from the classic TA, we find a well-defined constraint on the sum m_MW + m_M31 but are unable to resolve the individual masses. The bare TA sum indicated by NAM, (m_LG = 6 ± 1) × 10¹² M_☉ at 95% confidence, is consistent with Li & White (2008), who calibrated the TA against galaxy pairs drawn from the Millennium simulation to find, at 95% confidence, 1.9 × 10¹² M_☉ < M_LG < 1.0 × 10¹³ M_☉, with a median likelihood estimate of 5.7 × 10¹² M_☉ (the confidence interval in the latter is larger since it takes cosmic variance into account).

Zoom In Zoom Out Reset image size

In Figure 1, in the upper right, we show the results from a reduced version of our catalog which includes the LG actors but excludes the four external groups. The additional dynamical actors have broken the degeneracy in the TA, giving independent masses of 2.5 ± 1.5 × 10¹² M_☉ for the MW and 3.5 ± 1.0 × 10¹² M_☉ for M31. With the addition of the four external groups (Figure 1, lower left), the best mass for the MW increases to 3.5 ± 1.0 × 10¹² M_☉. This is consistent at the lower end with previous TA measurements of the total LG mass and the individual MW mass. When the transverse velocity constraints on M31, LMC, M33, IC10, and LeoI are added (Figure 1, lower right), the confidence intervals are broadened and the best-fit mass for MW decreases slightly, to 3.0 ± 1.5 × 10¹² M_☉, reflecting the fact that lower masses for MW are correlated to lower transverse velocities for M31 and other nearby galaxies.

As a check on the result for the LG, we use publicly available data from the Millennium Run (Springel et al. 2005; Lemson & the Virgo Consortium 2006; De Lucia & Blaizot 2007) and follow Li & White (2008) in generating mock LG catalogs satisfying the following conditions: we select type 0 or 1 galaxies with rotation velocities in the range 200 < V_max < 250 km s⁻¹ and a bulge-to-total luminosity ratio in the range 1.2–2.5. We then sub-select close pairs with comoving separation between 0.5 and 1 Mpc, negative relative peculiar velocity, and no massive companion with V_max > 150 km s⁻¹ within 2.5 Mpc from the center-of-mass. This led to the initial identification of about 100 Local Group candidates. Since the simulation overproduces satellite galaxies relative to observations (the "missing satellite problem"; see, e.g., Bullock 2012), and since the resulting dynamical complexity poses a special difficulty for NAM reconstructions, we additionally exclude catalogs where the gravitational acceleration on the mock-MW due to satellites is greater by more than a factor of six than what we expect from the observed distribution (assuming the distance and nominal masses listed in our catalog). This limited our set of mock catalogs to 32, with reasonably similar dynamics to what we expect to be the case with the LG, although in every case but two the dynamical complexity as defined above is greater than in the LG. To speed up computation time, within each of the 32 selected mock catalogs we include in the NAM computation only those galaxies and satellites within the distance to mock-M31, and all other galaxies out to 7 Mpc which produce an acceleration at least 5% that of M31. This reduced the number of particles to those that are most relevant to the dynamics of the principal actors – from as few as 9 to as many as 34, with an average across the catalogs of 19.

Contours in χ² for different halo masses of mock-MW and mock-M31 are shown for the four dynamically simplest catalogs in Figure 2. The predicted values for their masses are consistent with the true galaxy + halo masses from the simulation, within the 95% confidence range. We may expect the NAM prediction to be somewhat higher if it is sensitive to dark matter associated with the two principal actors but beyond their respective cutoff radii (at an overdensity of 200) that define their masses. A detailed comparison of the underlying dark matter distribution with the halo mass predictions may indicate whether this sensitivity is present. For the other 28 mock catalogs the mass predictions for both actors were likewise consistent with the 95% confidence intervals in all but three cases. The NAM-based mass predictions from the mock catalogs are thus accurate within statistical expectations.

Zoom In Zoom Out Reset image size

We have shown that the NAM offers a promising method of constraining the individual masses of the principal actors in the LG using an approach that investigates the behavior of χ² across a large ensemble of solutions. It is complementary to earlier work focusing on the detailed dynamical analysis of individual solutions.

The relatively large masses of MW and M31 suggested in this first-of-its-kind implementation of NAM are a point of interest. They may be due to the method's sensitivity to both the extended dark matter halo as well as purely dark concentrations between galaxies, since these otherwise invisible concentrations should leave a signature in galaxy proper motions. In this case, our result would suggest the presence of mass concentrations larger than previously suspected. Another possibility is that the result is biased by poor reconstruction of satellite paths, since solutions placing M31 at exactly the observed distance and redshift also tend to place the dwarf galaxies in the near vicinity of the MW at distances greater than the catalog distance. Related to the above point, our observational constraints assume a value of 220 km s⁻¹ for the circular velocity of the sun around the MW center; a higher value would reduce the approach velocity of M31, reducing the predicted mass.

Further improvements to this method are suggested by the above observations, and can be expected from a number of directions. On the computational side, increases in the efficiency of the solution finding algorithm will permit a larger number of time steps to be used and potentially allow for more complex orbits, improving in particular the reconstruction of nearby satellites. It will also permit a larger complement of dwarf and satellite galaxies to be included. Adding full three-dimensional proper motions of a larger number of nearby galaxies will certainly further constrain the likely masses. A comparison of the underlying dark matter distribution with the galaxy halos in the simulations will confirm whether NAM is potentially sensitive to extended distributions of dark matter beyond the nominal halo radii.

We thank Brent Tully for helpful suggestions that improved the presentation of the paper. This research was supported by the I-CORE Program of the Planning and Budgeting Committee, the Israeli Science Foundation (grants No. 1829/12 and No. 203/09), the German-Israeli Foundation for Research and Development, and the Asher Space Research Institute. V. Desjacques acknowledges support by the Swiss National Science Foundation.

Sections 1 and 2 are from Peebles et al. (2011). Section 3 is new to this study and shows how the action can be modified to accommodate redshift boundary conditions.

A.1. Review of the Theory

A.1.1. Equations of Motion

In a cosmologically flat universe, the expansion parameter satisfies

$$\begin{equation} {\dot{a}^2\over a^2} = {H_o^2\Omega \over a^3} + (1 - \Omega)H_o^2,\qquad {\ddot{a}\over a} = -{H_o^2\Omega \over 2a^3} + (1 - \Omega)H_o^2, \end{equation} \tag{ A1 }$$

with the present value a_o = 1. The equations of motion in physical length units are

$$\begin{equation} {d^2r_{i,k}\over dt^2} = \sum _{j\not= i} {Gm_j(r_{j,k}-r_{i,k})\over |r_i-r_j|^3} +(1-\Omega)H_o^2r_{i,k}. \end{equation} \tag{ A2 }$$

Changing variables to the comoving coordinates x_{i, k} = r_{i, k}/a(t) used here brings Equation (A2) to

$$\begin{equation} {d\over dt}a^2{d\,x_{i,k}\over dt}= {1\over a}\left[ \sum _{j} Gm_j{(x_{j,k}-x_{i,k})\over |x_i-x_j|^3} +{1\over 2}\Omega H_o^2x_{i,k}\right]. \end{equation} \tag{ A3 }$$

This is derived from the action

$$\begin{eqnarray} S & = & \int _0^{t_{o}}dt\left[\sum _i {m_i a^2 \dot{x}_i^2\over 2} + {1\over a}\left(\sum _{j\not= i}{Gm_im_j\over |x_i-x_j|} + {1\over 4}\sum _i m_i\Omega H_o^2 x_i^2\right)\right] \end{eqnarray} \tag{ A4 }$$

when the present positions are fixed, δx_i(t_o) = 0, and initial conditions satisfy

$$\begin{equation} a^2\dot{x}_i\rightarrow 0\quad \hbox{ at } a(t)\rightarrow 0. \end{equation} \tag{ A5 }$$

A.1.2. Discrete Representation

In a discrete representation, the coordinates are x_{i, k, n}, where i labels the particles, k = 1, 2, 3 the Cartesian coordinates, and 1 ⩽ n ⩽ n_x + 1 the time steps. The present positions $x_{i,k,n_x+1}$ are fixed and given. The relevant derivatives of the action are

$$\begin{equation} S_{i,k,n}={\partial S\over \partial x_{i,k,n}}, \quad S_{i,k,n;j,k^{\prime },n^{\prime }}={\partial ^2 S\over \partial x_{i,k,n}\partial x_{j,k^{\prime },n^{\prime }}}, \quad 1\le n,n^{\prime }\le n_x. \end{equation} \tag{ A6 }$$

If S is close to quadratic in the x_{i, k, n}, then position shifts δx_{i, k, n} to a solution at an extremum of S satisfy

$$\begin{equation} S_{i,k,n} + \sum _{j,k^{\prime },n^{\prime }}S_{i,k,n;j,k^{\prime },n^{\prime }}\delta x_{j,k^{\prime },n^{\prime }} =0. \end{equation} \tag{ A7 }$$

If the x_{i, k, n} are not close to a solution, then S is not close to quadratic in the x_{i, k, n}, but experience shows that coordinate shifts in the direction of δx_{i, k, n} move toward a solution.

Approximate the action (A2) as

$$\begin{eqnarray} &&\scriptsize{S = \sum _{i,k,n=1,n_x}{m_i\over 2}{(x_{i,k,n+1}- x_{i,k,n})^2\over (a_{n+1}-a_n)} \dot{a}_{n+1/2}a_{n+1/2}^2 + \sum _{i,j,n=1,n_x}{t_{n+1/2}-t_{n-1/2}\over a_n} \left[\sum _{j<i}{Gm_im_j\over |x_{i,n}-x_{j,n}|} + {1\over 4} \sum _i m_i\Omega H_o^2 x_{i,n}^2\right].} \end{eqnarray} \tag{ A8 }$$

The times t_{n ± 1/2} interpolate between the time steps at n and n ± 1 in leapfrog fashion. The approximation to the kinetic energy in Equation (A8) is motivated by linear perturbation theory, where dx/da is nearly independent of time, so (x_{n + 1} − x_n)/(a_{n + 1} − a_n) is a good approximation to dx/da at a_{n + 1/2}. The earliest time at which positions are computed is at a₁ > 0. The leapfrog back in time from a₁ is to a_1/2 = 0 = t_1/2. Recall that present positions at $a_{n_x+1}=1$ are given at $x_{i,k,n_x+1}$.

The derivative of the action with respect to the coordinates x_{i, k, n} for 1 ⩽ n ⩽ n_x, gives

$$\begin{eqnarray} &&\scriptsize{S_{i,k,n} = -{a_{n+1/2}^2\dot{a}_{n+1/2}\over a_{n+1} - a_n}(x_{i,k,n+1}-x_{i,k,n}) + {t_{n+1/2}-t_{n-1/2}\over a_n} \left[\sum _{j\not=i}Gm_j{x_{j,k,n}-x_{i,k,n}\over |x_{i,n}-x_{j,n}|^3} + {1\over 2}\Omega H_o^2 x_{i,k,n}\right]+ {a_{n-1/2}^2\dot{a}_{n-1/2}\over a_{n} - a_{n-1}}(x_{i,k,n}-x_{i,k,n-1}).}\nonumber\\ \end{eqnarray} \tag{ A9 }$$

When S_{i, k, n} = 0, this is a discrete approximation to the equation of motion (A3). The common factor m_i has been dropped to reduce clutter, which means $S_{i,k,n;j,k^{\prime },n^{\prime }}\not= S_{j,k^{\prime },n^{\prime };i,k,n}$. (The asymmetry is in the gravity term. There still is the symmetry $S_{i,k,n;i,k^{\prime },n^{\prime }} = S_{i,k^{\prime },n^{\prime };i,k,n}$.)

To simplify Equation (A9) and its derivatives wrt x_{i, k, n}, let

$$\begin{equation} F^+_n= {a_{n+1/2}^2\dot{a}_{n+1/2}\over a_{n+1} - a_n}, \quad F^-_n= {a_{n-1/2}^2\dot{a}_{n-1/2}\over a_{n} - a_{n-1}} = F^+_{n-1}, \quad {dt_n\over a_n} = {t_{n+1/2}-t_{n-1/2}\over a_n}. \end{equation} \tag{ A10 }$$

Note that

$$\begin{equation} F^-_1 = 0=F^+_0 \end{equation} \tag{ A11 }$$

follows from $a^2\dot{a}\rightarrow 0$ at a → 0. Also, write the acceleration (apart from the factor a²) of particle i due to the other particles $j\not= i$ as

$$\begin{equation} g_{i,k,n} = \sum _{j\not= i} Gm_j{x_{j,k,n}-x_{i,k,n}\over |x_{i,n}-x_{j,n}|^3}. \end{equation} \tag{ A12 }$$

All this notation brings Equation (A9) to

$$\begin{equation} S_{i,k,n} = -F^+_n(x_{i,k,n+1}-x_{i,k,n}) + F^-_n(x_{i,k,n}-x_{i,k,n-1}) + {dt_n\over a_n}\left[g_{i,k,n} + {1\over 2}\Omega H_o^2 x_{i,k,n}\right]. \end{equation} \tag{ A13 }$$

We need the derivatives of the action with respect to the positions of the particles. Let the derivative of the acceleration of particle i wrt the position of particle $j\not= i$ be

$$\begin{eqnarray} {\cal G}_{i,k,n;j,k^{\prime }} &=& {\partial g_{i,k,n}\over \partial x_{j,k^{\prime },n}}\nonumber\\ &=& Gm_j\left({\delta _{k,k^{\prime }}\over |x_{i,n} - x_{j,n}|^3} - 3{(x_{j,k,n}-x_{i,k,n})(x_{j,k^{\prime },n}-x_{i,k^{\prime },n}) \over |x_{i,n} - x_{j,n}|^5 }\right). \end{eqnarray} \tag{ A14 }$$

The derivative of the acceleration of particle i with respect to its own position is

$$\begin{equation} {\cal G}_{i,k,n;i,k^{\prime }} = {\partial g_{i,k,n}\over \partial x_{i,k^{\prime },n}} = -\sum _{j\not= i} {\cal G}_{i,k,n;j,k^{\prime }}. \end{equation} \tag{ A15 }$$

So the nonzero derivatives of Equation (A13) with respect to the coordinates are

$$\begin{eqnarray} S_{i,k,n;i,k,n+1} &=& - F^+_n, \qquad S_{i,k,n;i,k,n-1} = - F^-_n,\nonumber\\ S_{i,k,n;j,k^{\prime },n} &=& {dt_n\over a_n} {\cal G}_{i,k,n;j,k^{\prime }}, \ j\not= i, \nonumber\\ S_{i,k,n;i,k^{\prime },n}&=& \left(F^+_n + F^-_n\right)\delta _{k,k^{\prime }} + {dt_n\over a_n}\left[ {\cal G}_{i,k,n;i,k^{\prime }} +{1\over 2}\Omega H_o^2\delta _{k,k^{\prime }} \right]. \end{eqnarray} \tag{ A16 }$$

The goal is to use these second derivatives of the action to drive the first derivatives to zero at a stationary point, S_{i, k, n} = 0.

A.1.3. Leapfrog Integration Forward in Time

It is worth recording that when the action is at a stationary point, S_{i, k, n} = 0, a solution of Equation (A13) is equivalent to a standard leapfrog numerical integration of the equation of motion (A3) forward in time. In this leapfrog, the positions and velocities are computed at interleaved time steps as (in the notation in Equation (A10)),

$$\begin{equation} a_{n+1/2}^2\dot{x}_{i,k,n+1/2} = a_{n-1/2}^2\dot{x}_{i,k,n-1/2} + {dt_n\over a_n}\left[\sum _{j\not= i}g_{i,k,n;j} + {1\over 2}\Omega H_o^2 x_{i,k,n}\right], \end{equation} \tag{ A17 }$$

and

$$\begin{equation} x_{i,k,n+1} = x_{i,k,n} + \dot{x}_{i,k,n+1/2}(t_{n+1}-t_n)= x_{i,k,n} + \dot{x}_{i,k,n+1/2}(a_{n+1}-a_n)/\dot{a}_{n+1/2}, \end{equation} \tag{ A18 }$$

which in the notation in Equation (A10) is

$$\begin{equation} x_{i,k,n+1} = x_{i,k,n} + a_{n+1/2}^2\dot{x}_{i,k,n+1/2}/F^+_n. \end{equation} \tag{ A19 }$$

Equations (A17) and (A19) are equivalent to Equation (A13) at S_{i, k, n} = 0.

The difference from a conventional leapfrog integration is the boundary conditions, which here are the present positions and a condition on the initial velocities, as follows.

A.1.4. Initial Conditions

The representation of the mass distribution in the early universe by galaxy-size particles is certainly crude, but is a reasonably useful approximation that motivates the following consideration.

In linear perturbation theory for a continuous pressureless fluid, the unwanted decaying mode has a peculiar velocity that is decreasing as $v=a\dot{x}\propto 1/a(t)$. Equation (A5) formally eliminates this decaying mode. In the wanted growing mode, the coordinate position of a particle is changing with time as

$$\begin{equation} x(t) - x(0)\propto t^{2/3}\propto a(t). \end{equation} \tag{ A20 }$$

This implies that in the wanted growing mode, the left hand side of Equation (A3) is constant, meaning a²dx_{i, k}/dt ≃  constant × t, where t is the time measured from a = 0. This motivates approximating Equation (A3) at the time t_3/2 intermediate between the first two time steps in the leapfrog, a₁ and a₂, as

$$\begin{equation} a_{3/2}^2{dx_{i,k,3/2}\over dt} \simeq a_{3/2}^2\dot{a}_{3/2}{x_{i,k,2} - x_{i,k,1}\over a_2 - a_1}= {t_{3/2}\over a_1}\left[ \sum _{j} Gm_j{(x_{j,k,1}-x_{i,k,1})\over |x_i-x_j|^3} +{1\over 2}\Omega H_o^2x_{i,k,1}\right]. \end{equation} \tag{ A21 }$$

Recall that the half time step earlier than a₁ is at a_1/2 = 0 = t_1/2 = 0, where $a^2\dot{x}$ is supposed to vanish. Equation (A21) agrees with Equations (A11) and (A13) at S_{i, k, 1} = 0.

It should be noted that a₁ may be much larger than a₂ − a₁, meaning that the numerical solution commences at modest redshift with small time steps. However, t_3/2 is still the time from a_1/2 = 0 to the time midway between a₁ and a₂.

The prescription in Equation (A21) allows large peculiar velocities at high redshift. This is not inconsistent with Equation (A5); it corresponds to a large primeval departure from homogeneity. It does mean that one must select solutions that are judged to have realistic initial peculiar velocities (at a_3/2).

To summarize, the usual initial conditions—position and velocity—in a leapfrog integration of the equation of motion forward in time are replaced by the present position, $x_{i,k,n_x+1}$, and the relation in Equation (A21) between the two earliest positions, x_{i, k, 1} and x_{i, k, 2}, at times t₁ and t₂. If Equation (A20) is a good approximation, then this is equivalent to specifying the initial velocity, for then Equation (A21) determines $\dot{x}_{3/2}$. That is, the boundary conditions for a solution S_{i, k, n} = 0 are the present position and the time-variation of the initial velocity.

A.2. Method of Solution

A.2.1. Single Orbit Adjustment

Since we're adjusting only the orbit of particle i, we drop the label i and write the first derivative of S as

$$\begin{equation} {\partial S\over \partial x_{k,n}} = S_{k,n} = -F^+_n(x_{k,n+1}-x_{k,n}) + F^-_n(x_{k,n}-x_{k,n-1}) + {dt_n\over a_n}\left(\sum _{j\not= i}g_{k,n;j} + {\Omega H_o^2\over 2} x_{k,n}\right), \end{equation} \tag{ A22 }$$

and write the equation to be solved as

$$\begin{equation} S_{k,n} + \sum _{k^{\prime },n^{\prime }} S_{k,n,k^{\prime },n^{\prime }}\delta x_{k^{\prime },n^{\prime }} = 0, \qquad {\partial ^2 S\over \partial x_{k,n}\partial x_{k^{\prime },n^{\prime }}} = S_{k,n;k^{\prime },n^{\prime }} = S_{k^{\prime },n^{\prime };k,n}. \end{equation} \tag{ A23 }$$

The nonzero second derivatives are

$$\begin{equation} S_{k,n;k,n+1} = - F^+_n, \end{equation} \tag{ A24 }$$

$$\begin{equation} S_{k,n,k^{\prime },n} = \left(F^+_n + F^-_n\right)\delta _{k,k^{\prime }} + {dt_n\over a_n}\left(\sum _{j\not= i}{\cal G}_{k,n;jk^{\prime }} + {\Omega H_o^2\over 2}\delta _{k,k^{\prime }}\right). \end{equation} \tag{ A25 }$$

Since

$$\begin{equation} S_{k,n;k^{\prime },n^{\prime }} = 0 \hbox{ unless } n^{\prime }=n \hbox{ or else } n^{\prime } = n\pm 1 \hbox{ and } k^{\prime } = k, \end{equation} \tag{ A26 }$$

Equation (A23) is

$$\begin{equation} S_{k,n} + S_{k,n;k, n+1}\delta x_{k,n+1}+ \sum _{k^{\prime }}S_{k,n;k^{\prime },n}\delta x_{k^{\prime },n} + S_{k,n;k,n-1}\delta x_{k,n-1}= 0. \end{equation} \tag{ A27 }$$

Set n → n − 1 in this equation and rearrange it to

$$\begin{equation} \delta x_{k,n} = - {S_{k,n-1}+ \sum _{k^{\prime }} S_{k,n-1;k^{\prime },n-1}\delta x_{k^{\prime },n-1} + S_{k,n-1;k,n-2}\delta x_{k,n-2} \over S_{k,n-1;k,n} }. \end{equation} \tag{ A28 }$$

This gives δx_{k, n} in terms of δx_{k, n − 1} and δx_{k, n − 2}. On iterating, we get the form

$$\begin{equation} \delta x_{k,n} = A_{k,n} + \sum _{k^{\prime \prime }} B_{k,n;k^{\prime \prime }}\delta x_{k^{\prime \prime },1}. \end{equation} \tag{ A29 }$$

At n = 1, this is just

$$\begin{equation} A_{k,1} = 0,\qquad B_{k,1;k^{\prime \prime }} = \delta _{k,k^{\prime \prime }}. \end{equation} \tag{ A30 }$$

At the second time step from the start, n = 2, Equation (A28) is

$$\begin{equation} \delta x_{k,2} = - \left [S_{k,1}+ \sum _{k^{\prime }} S_{k,1;k^{\prime },1}\delta x_{k^{\prime },1} \right ]\left/\right. S_{k,1;k,2}, \end{equation} \tag{ A31 }$$

because there is no δx_{k, 0}. Comparing this with Equation (A29), we see that

$$\begin{equation} A_{k,2} = - S_{k,1}/S_{k,1;k,2}, \qquad B_{k,2;k^{\prime }} = - S_{k,1;k^{\prime },1}/S_{k,1;k,2}. \end{equation} \tag{ A32 }$$

At n ⩾ 3, the result of substituting the form (A29) into Equation (A28) is

$$\begin{eqnarray} &&\scriptsize{\delta x_{k,n} = - \left [S_{k,n-1} + \sum _{k^{\prime }} S_{k,n-1;k^{\prime },n-1}(A_{k^{\prime },n-1} + \sum _{k^{\prime \prime }} B_{k^{\prime },n-1;k^{\prime \prime }}\delta x_{k^{\prime \prime },1}) + S_{k,n-1;k,n-2} \left(A_{k,n-2} + \sum _{k^{\prime \prime }} B_{k,n-2;k^{\prime \prime }}\delta x_{k^{\prime \prime },1}\right)\right ]\left/\right. S_{k,n-1;k,n}.} \nonumber\\ \end{eqnarray} \tag{ A33 }$$

So at 3 ⩽ n ⩽ n_x

$$\begin{eqnarray} A_{k,n} &= - {S_{k,n-1} + \sum _{k^{\prime }} S_{k,n-1;k^{\prime },n-1}A_{k^{\prime },n-1} + S_{k,n-1;k,n-2}A_{k,n-2}\over S_{k,n-1;k,n}},\nonumber\\ B_{k,n;k^{\prime \prime }} &= - {\sum _{k^{\prime }} S_{k,n-1;k^{\prime },n-1}B_{k^{\prime },n-1;k^{\prime \prime }} + S_{k,n-1;k,n-2}B_{k,n-2;k^{\prime \prime }}\over S_{k,n-1;k,n}}. \nonumber\\ \end{eqnarray} \tag{ A34 }$$

Equation (A34), with Equations (A30) and (A32), gives the A_{k, n} and $B_{k,n;k^{\prime \prime }}$, for all the time steps in terms of the input derivatives of the action. We use these to find the position shifts δx_{k, n} in Equation (A29). All that remains before this can be done is to compute the δx_{k, 1}, which we find by setting n = n_x in Equation (A27) and recalling that $\delta x_{k,n_x+1}= 0$:

$$\begin{eqnarray} 0 &=& S_{k,n_x} + \sum _{k^{\prime }}S_{k,n_x;k^{\prime },n_x}\delta x_{k^{\prime },n_x} + S_{k,n_x;k,n_x-1}\delta x_{k,n_x-1}\nonumber\\ &=& S_{k,n_x} + \sum _{k^{\prime }}S_{k,n_x;k^{\prime },n_x} \left [ A_{k^{\prime },n_x} + \sum _{k^{\prime \prime }} B_{k^{\prime },n_x;k^{\prime \prime }}\delta x_{k^{\prime \prime },1}\right ] + S_{k,n_x;k,n_x-1}\left [ A_{k,n_x-1} + \sum _{k^{\prime \prime }} B_{k,n_x-1;k^{\prime \prime }}\delta x_{k^{\prime \prime },1}\right ]. \end{eqnarray} \tag{ A35 }$$

So write this as

$$\begin{eqnarray} 0 &=& T_k + \sum _{k^{\prime }}T_{k,k^{\prime \prime }}\delta x_{k^{\prime \prime },1},\nonumber\\ T_k &=& S_{k,n_x} + \sum _{k^{\prime }}S_{k,n_x;k^{\prime },n_x}A_{k^{\prime },n_x} + S_{k,n_x;k,n_x-1}A_{k,n_x-1}, \nonumber\\ T_{k,k^{\prime \prime }} &=& \sum _{k^{\prime }}S_{k,n_x;k^{\prime },n_x}B_{k^{\prime },n_x;k^{\prime \prime }} + S_{k,n_x;k,n_x-1} B_{k,n_x-1;k^{\prime \prime }}, \end{eqnarray} \tag{ A36 }$$

solve this 3 × 3 set of equations for the δx_{k, 1}, and then get the rest of the δx_{k, n} from Equation (A29).

A.3. Redshift Boundary Condition

Solutions to the equations of motion in the method outlined in Appendix A.2 above satisfy the constraint that the distances at the present epoch are fixed. However, since redshifts are known more accurately than distances, and since we may want to explore a larger space of solutions, it may desirable to recast the problem with fixed redshifts, with the distances at the present epoch emerging as predictions (note that the choice of boundary condition can be made particle-by-particle). This can be accomplished through a partial transformation of coordinates that exchanges radial distances with radial velocities while leaving the angular position coordinates unchanged. Details are given in Phelps (2002); the procedure is summarized below.

The change to radial velocity coordinates is carried out in the Hamiltonian frame through a canonical transformation of the conjugate variables (positions and momenta). The generating function of the transformation, which is added to the action outside the integral, is

$$\begin{equation} {-}p_{i,r,nx+1} q_{i,r,nx+1}, \end{equation} \tag{ A37 }$$

where q_r and p_r are the conjugate distance and momentum in the radial direction relative to the reference galaxy. With the addition of the generating function, the problem becomes equivalent to one expressed with a new set of conjugate coordinates Q_r and P_r, where Q_r is the radial momentum and P_r is the radial distance. In these coordinates, the boundary term in the variational derivative PδQ vanishes at z = 0 when the angular positions and the radial velocity vanish. However, since the transformed Hamiltonian cannot be expressed analytically (the gravitational term cannot be written out in terms of Q), the computation must be carried out in the original coordinate system, expressed in the Lagrangian frame, where the new boundary condition must be imposed by hand to recover the correct equations of motion. That constraint takes the form of additional terms in the action, as follows.

A.3.1. Modification to the Action

The modified action is

$$\begin{eqnarray} S &=& {1 \over 2} H_0 x_{i,r,nx+1}^2 - cz x_{i,r,nx+1} + x_{i,r,nx+1} \boldsymbol {v}_{mw,nx+1} \cdot \hat{r}_{i,nx+1} +\sum _{i,k,n=1,n_x}{1\over 2}{(x_{i,k,n+1}- x_{i,k,n})^2\over (a_{n+1}-a_n)} \dot{a}_{n+1/2}a_{n+1/2}^2 \nonumber \\ &&+\sum _{i,j,n=1,n_x}{t_{n+1/2}-t_{n-1/2}\over a_n} \left[\sum _{j<i}{Gm_j\over |x_{i,n}-x_{j,n}|} + {1\over 4} \sum _i \Omega H_o^2 x_{i,n}^2\right]. \end{eqnarray} \tag{ A38 }$$

The gradient at the final time step is

$$\begin{eqnarray} &&\scriptsize{S_{i,k,n} = H_0 x_{i,r,nx+1} - cz + \boldsymbol {v}_{mw,nx+1} \cdot \hat{r}_{i,nx+1} -F^+_n(x_{i,k,n+1}-x_{i,k,n}) + F^-_n(x_{i,k,n}-x_{i,k,n-1}) + {dt_n\over a_n}\left[g_{i,k,n} + {1\over 2}\Omega H_o^2 x_{i,k,n}\right].}\nonumber\\ \end{eqnarray} \tag{ A39 }$$

The modified second derivatives of the action at the final two timesteps are

$$\begin{equation} S_{i,k,nx;i,k,nx+1} = - {\hat{x}_{i,k,nx+1} \over dt_{nx+1}}. \end{equation} \tag{ A40 }$$

Note that a new term in the action $S_{k,n_x+1;k,n_x+1}$ is also present, but it is not used in the computation that follows.

Since the boundary condition is now velocity-limited, $\delta x_{k,n_x+1}\not= 0$, and we must add a new term in Equation (A37) arising from the final two timesteps:

$$\begin{eqnarray} 0 &=& S_{k,n_x} + \sum _{k^{\prime }}S_{k,n_x;k^{\prime },n_x}\delta x_{k^{\prime },n_x} + S_{k,n_x;k,n_x-1}\delta x_{k,n_x-1} + S_{k,n_x;k,n_x+1}\delta x_{k,n_x+1} \nonumber \\ &=& S_{k,n_x} + \sum _{k^{\prime }}S_{k,n_x;k^{\prime },n_x} \left [ A_{k^{\prime },n_x} + \sum _{k^{\prime \prime }} B_{k^{\prime },n_x;k^{\prime \prime }}\delta x_{k^{\prime \prime },1}\right ]+ S_{k,n_x;k,n_x-1}\left [ A_{k,n_x-1} + \sum _{k^{\prime \prime }} B_{k,n_x-1;k^{\prime \prime }}\delta x_{k^{\prime \prime },1}\right ] \nonumber\\ &&+ S_{k,n_x;k,n_x+1} \left [ A_{k,n_x+1} + \sum _{k^{\prime \prime }} B_{k,n_x+1;k^{\prime \prime }}\delta x_{k^{\prime \prime },1}\right ]. \end{eqnarray} \tag{ A41 }$$

This is now written as

$$\begin{eqnarray} &&0 = T_k + \sum _{k^{\prime }}T_{k,k^{\prime \prime }}\delta x_{k^{\prime \prime },1}, \nonumber \\ && T_k = S_{k,n_x} + \sum _{k^{\prime }}S_{k,n_x;k^{\prime },n_x}A_{k^{\prime },n_x} + S_{k,n_x;k,n_x-1}A_{k,n_x-1}, + S_{k,n_x;k,n_x+1}A_{k,n_x+1},\nonumber \\ && T_{k,k^{\prime \prime }} = \sum _{k^{\prime }}S_{k,n_x;k^{\prime },n_x}B_{k^{\prime },n_x;k^{\prime \prime }} + S_{k,n_x;k,n_x-1} B_{k,n_x-1;k^{\prime \prime }}, + S_{k,n_x;k,n_x+1}B_{k,n_x+1;k^{\prime \prime }}. \end{eqnarray} \tag{ A42 }$$

As before, this 3 × 3 set of equations is solved for the δx_{k, 1}, and we get the rest of the δx_{k, n} from Equation (A29). However, we will find below that the boundary conditions will enable us to replace the δx_{nx + 1} above with functions of δx_nx.

A.3.2. Modification to the Position Shifts

The change of boundary conditions implies a relationship between the position shifts δx_{k, nx + 1} and δx_{k, nx} that will further modify the new terms above. The relationship is simplest to see in the one-dimensional case (setting y = z = 0), where

$$\begin{eqnarray} cz &=& v_{nx+1} + H_0 x_{nx+1} - v_{mw,nx+1}. \nonumber \\ &=& {x_{nx+1} - a_{nx} x_{nx} \over dt_{nx+1/2}} + H_0 (x_{nx+1} - x_{mw,nx+1}) - v_{mw,nx+1}, \end{eqnarray} \tag{ A43 }$$

where we continue to drop the particle subscript unless referring to the reference galaxy, mw, since we are working on one particle at a time while holding the rest of the catalog fixed.

After the variation,

$$\begin{eqnarray} cz^{\prime } &=& {x_{nx+1} + \delta x_{nx+1} - a_{nx} x_{nx} - a_{nx} \delta x_{nx} \over dt_{nx+1/2}} + H_0 (x_{nx+1} + \delta x_{nx+1} - x_{mw,nx+1}) - v_{mw,nx+1}. \end{eqnarray} \tag{ A44 }$$

Recall that the orbit of the reference galaxy, and thus its contribution to the redshift, is held constant during the variation.

Setting cz' = cz gives

$$\begin{eqnarray} \delta x_{nx+1} &=& {a_{nx} \delta x_{nx} \over 1 + H_0 dt_{nx+1/2}}. \end{eqnarray} \tag{ A45 }$$

In three dimensions, the angular coordinates of the particle must also remain unchanged after the variation, and so now there are three equations constraining the coordinate shifts at the final two timesteps:

$$\begin{equation} cz^{\prime } = cz,\qquad \theta ^{\prime } = \theta ,\qquad \phi ^{\prime } = \phi , \end{equation} \tag{ A46 }$$

where θ and ϕ are defined relative to the reference galaxy at z = 0:

$$\begin{eqnarray} \tan \theta _{nx+1} &=& {x_{2,nx+1} - x_{mw,2,nx+1} \over x_{1,nx+1} - x_{mw,1,nx+1}} \end{eqnarray} \tag{ A47 }$$

$$\begin{eqnarray} \tan \phi _{nx+1} &=& {x_{3,nx+1} - x_{mw,3,nx+1} \over \sqrt{(x_{1,nx+1} - x_{mw,1,nx+1})^2 + (x_{2,nx+1} - x_{mw,2,nx+1})^2 }}. \end{eqnarray} \tag{ A48 }$$

The θ and ϕ constraints impose a relationship between the three δx_{k, nx + 1}. The θ constraint gives

$$\begin{eqnarray} \theta ^{\prime }_{nx+1} &=& \theta _{nx+1} \nonumber \\ {x^{\prime }_2 \over x^{\prime }_1} &=& {x_2 \over x_1} \nonumber \\ {x_{2,nx+1} - x_{mw,2,nx+1} + \delta x_{2,nx+1} \over x_{1,nx+1} - x_{mw,1,nx+1} + \delta x_{1,nx+1}} &=& {x_{2,nx+1} - x_{mw,2,nx+1} \over x_{1,nx+1} - x_{mw,1,nx+1}} \nonumber \\ \delta x_{2,nx+1} &=& {x_{2,nx+1} \over x_{1,nx+1}} \delta x_{1,nx+1}. \end{eqnarray} \tag{ A49 }$$

The combined θ and ϕ constraints further give

$$\begin{eqnarray} \phi ^{\prime }_{nx+1} &=& \phi _{nx+1} \nonumber \\ {x_1^{\prime 2} + x_2^{\prime 2} \over x_3^{\prime 2}} &=& {x_1^2 + x_2^2 \over x_3^2} \nonumber \\ \delta x_{3,nx+1} &=& {x_{3,nx+1} \over x_{1,nx+1}} \delta x_{1,nx+1}. \end{eqnarray} \tag{ A50 }$$

As for the third constraint, the redshift in three dimensions is

$$\begin{eqnarray} cz &=& (\boldsymbol {v}_{nx+1} - \boldsymbol {v}_{mw,nx+1}) \cdot \hat{r}_{nx+1} + H_0 r_{nx+1} \nonumber \\ &=& \sum _k \left ({x_{k,nx+1} - a_{nx} x_{k,nx} \over dt_{nx+1/2}} - \boldsymbol {v}_{mw,nx+1} \right) \cdot \hat{x}_{k,nx+1} + H_0 \sqrt{\sum _k (x_{k,nx+1} - x_{mw,k,nx+1})^2}, \end{eqnarray} \tag{ A51 }$$

where

$$\begin{equation} \hat{x}_{k,nx+1} = {x_{k,nx+1} \over r_{nx+1}} = {x_{k,nx+1} \over \sqrt{\sum _{k^{\prime }} (x_{k^{\prime },nx+1} -x_{mw,k^{\prime },nx+1})^2}}. \end{equation} \tag{ A52 }$$

After the position shifts,

$$\begin{eqnarray} &&\scriptsize{cz^{\prime } = \sum _k \left ({x_{k,nx+1} + \delta x_{k,nx+1} - a_{nx} x_{k,nx} - a_{nx} \delta _{k,nx} \over dt_{nx+1/2}} - \boldsymbol {v}_{mw,k,nx+1} \right) \cdot \hat{x}^{\prime }_{k,nx+1} + H_0 \sqrt{\sum _k (x_{k,nx+1} - x_{mw,k,nx+1} + \delta x_{k,nx+1}) ^2},}\nonumber\\ \end{eqnarray} \tag{ A53 }$$

where

$$\begin{eqnarray} \hat{x}^{\prime }_{k,nx+1} &=& {x_{k,nx+1} + \delta x_{k,nx+1} \over \sqrt{\sum _{k^{\prime }} (x_{k^{\prime },nx+1} - x_{mw,k^{\prime },nx+1} + \delta x_{k^{\prime },nx+1})^2}}. \end{eqnarray} \tag{ A54 }$$

Setting cz' = cz and using the approximation δx_k/x_k ≪ 1, we find, after some algebra,

$$\begin{equation} \delta x_{k,nx+1} = {x_{k,nx+1} \over \alpha } \sum _{k^{\prime }} x_{k^{\prime },nx+1} a_{nx} \delta x_{k^{\prime },nx}, \end{equation} \tag{ A55 }$$

where

$$\begin{eqnarray} \alpha &=& \sum _{k^{\prime \prime }} x_{k^{\prime \prime },nx+1}(x_{k^{\prime \prime },nx+1} (2 -\beta) - a_{nx} x_{k^{\prime \prime },nx} - v_{mw,k^{\prime \prime },nx+1} dt_{nx+1/2}) \end{eqnarray} \tag{ A56 }$$

$$\begin{eqnarray} \beta &=& 1 - H_0 dt_{nx+1/2} - \sum _{k^{\prime \prime \prime }} {x_{k^{\prime \prime \prime },nx+1} a_{nx} x_{k^{\prime \prime \prime },nx} + x_{k^{\prime \prime \prime },nx+1} v_{mw,^{\prime \prime \prime },nx+1} dt_{nx+1/2} \over r_{nx+1}^2}. \end{eqnarray} \tag{ A57 }$$

It can be shown that this reduces to the one-dimensional case when we set two of the three position coordinates to zero.

With the above, the position shifts δx_{k, nx + 1} can be replaced in Equations (A41)–(A42) by their equivalent expressions in terms of δx_{k, nx}:

$$\begin{eqnarray} 0 &=& S_{k,n_x} + \sum _{k^{\prime }}S_{k,n_x;k^{\prime },n_x} \left [ A_{k^{\prime },n_x} + \sum _{k^{\prime \prime }} B_{k^{\prime },n_x;k^{\prime \prime }}\delta x_{k^{\prime \prime },1}\right ] + S_{k,n_x;k,n_x-1}\left [ A_{k,n_x-1} + \sum _{k^{\prime \prime }} B_{k,n_x-1;k^{\prime \prime }}\delta x_{k^{\prime \prime },1}\right ] \nonumber\\ && + S_{k,n_x;k,n_x+1} x_{k,nx+1} \alpha ^{-1} \sum _{k^{\prime }} x_{k^{\prime },nx+1} \left [ A_{k^{\prime },n_x} + \sum _{k^{\prime \prime }} B_{k^{\prime },n_x;k^{\prime \prime }}\delta x_{k^{\prime \prime },1}\right ]. \end{eqnarray} \tag{ A58 }$$

$$\begin{eqnarray} 0 &=& T_k + \sum _{k^{\prime }}T_{k,k^{\prime \prime }}\delta x_{k^{\prime \prime },1}, \nonumber \\ T_k &=& S_{k,n_x} + \sum _{k^{\prime }}S_{k,n_x;k^{\prime },n_x}A_{k^{\prime },n_x} + S_{k,n_x;k,n_x-1}A_{k,n_x-1}, + S_{k,n_x;k,n_x+1} x_{k,nx+1} \alpha ^{-1} \sum _{k^{\prime }} x_{k^{\prime },nx+1} A_{k^{\prime },n_x} \nonumber\\ T_{k,k^{\prime \prime }} &=& \sum _{k^{\prime }}S_{k,n_x;k^{\prime },n_x}B_{k^{\prime },n_x;k^{\prime \prime }} + S_{k,n_x;k,n_x-1} B_{k,n_x-1;k^{\prime \prime }}, + S_{k,n_x;k,n_x+1} x_{k,nx+1} \alpha ^{-1} \sum _{k^{\prime }} x_{k^{\prime },nx+1} B_{k^{\prime },n_x;k^{\prime \prime }}. \end{eqnarray} \tag{ A59 }$$

Note that the actual value of the redshift has been nowhere specified and must be arranged by hand, as by choosing the initial trial orbits to have the input redshifts. However, the approximation δx_k/x_k ≪ 1 used to arrive at Equation (A55) can fail, particularly in the first few iterations of the action, and so in practice the redshifts may drift away from their input values as the stationary point in the action is reached. This drift is corrected by periodically adjusting the particle positions at timestep nx + 1 according to Equation (A51).