# The SIBELIUS Project: E Pluribus Unum

Till Sawala, Stuart Mcalpine, Jens Jasche, Guilhem Lavaux, Adrian Jenkins, Peter H. Johansson, Carlos S. Frenk

► To cite this version:

Till Sawala, Stuart Mcalpine, Jens Jasche, Guilhem Lavaux, Adrian Jenkins, et al.. The SIBELIUS Project: E Pluribus Unum. Monthly Notices of the Royal Astronomical Society, 2021, 509 (1), pp.1432-1446. 10.1093/mnras/stab2684. hal-03191135

HAL Id: hal-03191135

https://hal.science/hal-03191135

Submitted on 4 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The SIBELIUS Project: E Pluribus Unum

Till Sawala<sup>1,1*</sup> Stuart McAlpine<sup>1</sup> Jens Jasche,<sup>2</sup> Guilhem Lavaux<sup>1,3</sup> Adrian Jenkins<sup>1,4</sup>
Peter H. Johansson<sup>1</sup> and Carlos S. Frenk<sup>4</sup>

Accepted 2021 September 8. Received 2021 September 8; in original form 2021 March 21

## ABSTRACT

We introduce 'Simulations Beyond The Local Universe' (SIBELIUS) that connect the Local Group (LG) to its cosmic environment. We show that introducing hierarchical small-scale perturbations to a density field constrained on large scales by observations provides an efficient way to explore the sample space of LG analogues. From more than 60 000 simulations, we identify a hierarchy of LG characteristics emanating from different scales: the total mass, orientation, orbital energy, and the angular momentum are largely determined by modes above $\lambda = 1.6$ comoving Mpc (cMpc) in the primordial density field. Smaller scale variations are mostly manifest as perturbations to the MW–M31 orbit, and we find that the observables commonly used to describe the LG – the MW–M31 separation and radial velocity – are transient and depend on specifying scales down to 0.2 cMpc in the primordial density field. We further find that the presence of M33/LMC analogues significantly affects the MW–M31 orbit and its sensitivity to small-scale perturbations. We construct initial conditions that lead to the formation of an LG whose primary observables precisely match the current observations.

**Key words:** methods: numerical – galaxies: formation – Local Group – dark matter – large-scale structure of Universe.

## 1 INTRODUCTION

The Local Group (LG), comprising our own Milky Way (MW), the neighbouring Andromeda galaxy (M31), the Large Magellanic Cloud (LMC) and M33, and at least a hundred lower mass galaxies (e.g. Karachentsev et al. 2004; McConnachie 2012; Dooley et al. 2017), has been the place of important discoveries throughout the history of modern cosmology. The correct understanding of galaxies as distinct stellar systems was first proposed in the early eighteenth century, by Immanuel Kant and Thomas Wright of Durham, in reference to M31 (e.g. Whitrow 1967). The question of whether galaxies are truly distinct systems culminated in the 'Great Debate of Curtis and Shapley in 1920, and was finally settled by measurements of the distances to cepheids in three LG galaxies, NGC 6822 (Hubble 1925), M33 (Hubble 1927), and M31 (Hubble 1929). In the 1970s, the measurement of the H I rotation curve of M31 suggested the existence of dark matter inside galaxies (Rubin & Ford 1970), and by the 1980s, the velocity dispersions of LG dwarf spheroidal galaxies indicated exceptionally high dark matter fractions (Faber & Lin 1983).

Today the LG continues to be one of the most sensitive laboratories to study astrophysical phenomena (Mateo 1998) and cosmological models. Its assembly history, dynamical evolution, and total mass are linked to important open questions in cosmology. In recent years, semi-analytical models (e.g. Bullock, Kravtsov & Weinberg 2001; Benson et al. 2002; Somerville 2002) and hydrodynamical

simulations that include baryonic physics (e.g. Okamoto, Gao & Theuns 2008; Sawala et al. 2013, 2016b; Garrison-Kimmel et al. 2019) have shown how the combined effects of reionization and supernova feedback suppress star formation in, and even the growth of, low-mass dark matter halos. At the same time, the rapid ejection of interstellar gas from the centre of low-mass halos following bursts of star formation may explain the transformation from a steep, NFW-like cusp (Navarro, Frenk & White 1996) into a shallower dark matter density profile (Eke, Navarro & Frenk 1998; Read & Gilmore 2005; Governato et al. 2010; Teyssier et al. 2013). While some of the astrophysical processes are still not completely understood, let alone numerically resolved, these baryonic effects have at least offered plausible explanations to some of the challenges raised to the ΛCDM model (Moore et al. 1999; Boylan-Kolchin, Bullock & Kaplinghat 2011). However, the uncertainty in the total mass of the LG of at least a factor of 2 (Li & White 2008; van der Marel et al. 2012) remains a major caveat for all predictions of the amount of substructure.

In the meantime, another question that has been asked since at least the seventeenth century, but which dates back to the Copernican revolution of the sixteenth century, has become more pertinent: How special is our own vantage point on the Universe, and what can we induce about cosmology by closely studying the finite patch of the Universe that surrounds us? Even more pressing, perhaps, is the inverse: How plausible is the formation of the LG in the assumed cosmological model? It has recently been pointed out that the LG has characteristics that make it exceedingly rare in ΛCDM: Only ~5–10 per cent of Milky-Way-sized halos in ΛCDM simulations

© The Author(s) 2021.

<footnote n="1"><sup>1</sup>Department of Physics, Gustaf Hällströmin katu 2, University of Helsinki, FI-00170 Helsinki, Finland</footnote>

<footnote n="2"><sup>2</sup>The Oskar Klein Centre, Department of Physics, Stockholm University, Albanova University Center, SE-106 91 Stockholm, Sweden</footnote>

<footnote n="3"><sup>3</sup>CNRS and Sorbonne Université, UMR7095, Institut d'Astrophysique de Paris, F-75014 Paris, France</footnote>

<footnote n="4"><sup>4</sup>Institute for Computational Cosmology, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote>*E-mail: till.sawala@helsinki.fi</footnote>
have satellites as massive as the LMC and SMC (Benson et al. 2002; Boylan-Kolchin et al. 2011; Tollerud et al. 2011; Robotham et al. 2012), the measured MW–M31 orbit is unusually radial (e.g. Fattahi et al. 2015), and only a very small fraction of the ΛCDM halos have satellite systems as anisotropic as the LG’s (e.g. Ibata et al. 2013; Pawlowski et al. 2015).

The above-mentioned baryonic physics seem to offer no possible path to answer these questions (Bullock & Boylan-Kolchin 2017; Pawlowski 2018). Furthermore, ‘simple’ modifications or extensions to the ΛCDM model, such as warm dark matter (e.g. Lovell et al. 2011), a modification of the primordial power spectrum (e.g. Enqvist, Sawala & Takahashi 2020), or self-interacting DM (e.g. Lovell et al. 2020) are no convincing alternatives. If the particular features of the Local Universe are as exceptional as they are claimed to be, they may hint at changes to our cosmological models that go beyond what has been so far envisioned. Alternatively, they should lead to an understanding of the unique but perhaps not implausible environment in which the Local Universe has evolved.

The aim of the SIBELIUS (Simulations Beyond the Local Universe) project is to elucidate the co-evolution of the LG, the Local Universe, and the large-scale structure (LSS). Similar constrained local universe simulations include the CLUES (Gottloeber, Hoffman & Yepes 2010; Libeskind et al. 2010; Carlesi et al. 2016; Sorce et al. 2016) and HESTIA (Libeskind et al. 2020) projects, which embed several LG analogues within the observed LSS. However, our method allows us to systematically examine a much larger set of initial conditions, to not only create a precise match to the LG’s orbit, but to also to reveal the relation between its parameters, and their dependence on the initial density field.

In this paper, we present our new method for creating initial conditions that reproduce the kinematics of the LG, and describe the sensitivity of the LG observables to different scales in the initial conditions. We make extensive use of the methods described in Sawala et al. (2021) to create a series of hierarchical, random perturbation of the primordial density field to unveil the scales at which the observed properties of the LG are set, and identify the properties that characterize it.

This paper is organized as follows. In Section 2, we describe the methods: the parametrization of the phase information (2.1), the LSS constraints (2.2), their octree function representation (2.3), and the implementation of multilevel random variations (2.4). We describe the numerical simulations used (2.5) and the identification of LG analogues (2.6). Sections 3–5 present the results. In Section 3.1, we discuss the number of LG analogues identified when all unconstrained scales are varied, and in Section 3.2, we discuss multiple occurrences of the same object in our large random sample. In Section 3.3, we introduce a metric to compare LG analogues to observations. In Section 4.1, we identify the ‘characteristics’ of the LG as those quantities that are defined on the largest scales. In Section 4.2, we discuss the frequency of LMC/M33 analogues. In Section 5, we discuss smaller scale variation, and determine in Section 5.1 the scale that sets the LG’s orbital phase.

## 2 METHODS

All results in this paper are based on dark matter only (DMO) zoom-in simulations using initial conditions generated with the IC_GEN code (Jenkins 2010) for a ΛCDM Universe with the following cosmological parameters: $\Omega_0 = 0.307$, $\Omega_\Lambda = 0.693$, $\sigma_8 = 0.8288$, and $h = 0.6777$ (Planck Collaboration XVI 2014). The initial displacement field is calculated using second-order Lagrangian perturbation theory at $z = 169$, although the simulations are not sensitive to the starting

redshift. The unperturbed particle positions in the high-resolution region are realized as glasses, with identifiers that encode the initial (unperturbed) position of the particles. We use quadratic interpolation to compute the displacements.

### 2.1 Phase Information

In the ΛCDM model, the primordial density fluctuations after inflation are homogeneous and Gaussian-distributed. To specify the initial conditions for a cosmological ΛCDM simulation, in addition to the cosmological parameters, we define the dimensions of the periodic volume, the power spectrum, and choose the phases for the realization. The information describing the particular choices for the phases can be represented by a Gaussian white noise field. Convolving this white noise field in real space with a suitable, real non-negative filter, gives the ΛCDM density field (Salmon 1996).

In this work, we wish to construct a cosmological volume that contains a region resembling the observed LSS around us. We use the $S_8$ octree orthogonal basis function set introduced in Jenkins (2013) as the building blocks for representing the phase information. We choose this representation over the more obvious choice of Fourier modes because the octree representation is designed to be more efficient for building the kind of multi-scale initial conditions we require in the SIBELIUS project. The octree basis functions have compact support and their individual contributions to the ΛCDM density field are also well localized. This contrasts with Fourier modes that have equal weight everywhere. From a purely computational point of view, using Fourier modes is more expensive when we wish to specify initial conditions at high resolution only within a small region of a large periodic volume.

In Section 2.2, we will describe how we use LSS constraints to set the amplitudes for some of the physically largest octree functions. On the smallest scales, however, the LSS constraints are unimportant and the field is effectively unconstrained. By default, we set the unconstrained coefficients using the Panphasia Gaussian white noise field (Jenkins 2013). The choice of location within the Panphasia field is arbitrary but by making this choice we determine the physical sizes and locations of the octree cells relative to the simulation volume. We choose to embed the constraints that represent the local LSS out to a radius of ~300 cMpc within a cubic volume of side length 1000 comoving Mpc (hereafter cMpc). For consistency with Sawala et al. (2021), where we had identified a $100^3$ cMpc$^3$ volume with $12^3$ cells, we here identify the $1000^3$ cMpc$^3$ volume with a set of $15^3$ octree cells, yielding a power-of-2 ratio between the cell sizes in the two simulation volumes.

The number of root cells and the level at which they are defined within PANPHASIA provide an anchor point for the physical cell size: In all our simulations, the cubic cells at level 18 have a volume of ~66.67³ cMpc³; cells at Level L have side length $66.67 \times 2^{(18-L)}$ cMpc. In order to connect levels of the octree with a physical scale, in Sawala et al. (2021), we also computed $\lambda_{cut}(L)$, the wavelength at which the amplitude of fluctuations of the ΛCDM power spectrum is reduced by a factor of four when power at levels L and above are set to zero. Jenkins (2013) provides a more in-depth description of the public code PANPHASIA, while Jenkins & Booth (2013) contains a user guide.

In Sawala et al. (2021), we already explored the propagation of changes to the primordial density field to structures at $z = 0$. In this work, we apply this knowledge, first to explore the ensemble of LG analogues that can form within the constraints, and then to create variations of individual LG analogues by randomizing different levels of the white noise field.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Illustration of the phase information used in the simulations. For each constrained simulation (shown in Fig. 2), the set of constraints sets the coefficients up to level 21, with the remaining levels taken from PANPHASIA (yellow). For the 3.2-cMpc variants (top two rows in Fig. 3), we preserve the constraints, and vary the phase information in the unconstrained levels with different regions from PANPHASIA (red). For each set of 0.8-cMpc variations (bottom two rows in Fig. 3), we also preserve the information from a sample of the 3.2-cMpc variations, and independently randomize a small subregion at levels 24 and 25 (red). Note that, in the actual simulations, the $1000^3$ cMpc$^3$ volume is represented by $15^3$ cells at level 18, and the subregions are close to spherical rather than $n$-cubes.</figcaption></figure>

As illustrated in Fig. 1, the cosmological constraints described in Section 2.2 are located on levels 19–21 of the octree, constraining modes down to $\lambda_{\text{cut}}(L=21) = 6.48$ cMpc. As shown in Sawala et al. (2021), this is sufficient to determine the existence of individual $z=0$ halos above $\sim 10^{14}$M$_\odot$. For halos of $\sim 10^{15}$M$_\odot$, i.e. the size of the Virgo Cluster, there is an expected residual scatter in position of $\sim 1$ Mpc, and a scatter in mass of $\sim 5$ per cent, which are both below the observational uncertainties. However, the existence of unique $10^{12}$M$_\odot$ halos, i.e. in the mass range of both the MW and M31, requires defining the phase information to smaller scales, i.e. to higher levels of the octree. This offers an opportunity: Modifications to the phase information on the right scale allow us to create different ensembles of MW/M31 mass halos without affecting the observed LSS.

## 2.2 Constraints for the Local LSS

The phase information for the Sibelius project is created in multiple steps. Our starting point is the phase information obtained using the BORG algorithm (Bayesian Origin Reconstruction from Galaxies; Jasche & Wandelt 2013; Jasche & Lavaux 2019), which is represented on a regular $256^3$ mesh covering a volume of $1000^3$ cMpc$^3$, i.e. each volume element covers $\sim 3.9$ cMpc$^3$.

Like the more common Reverse Zel'dovich Algorithm (RZA; e.g. Hoffman & Ribak 1991; Sorce et al. 2014), the BORG algorithm is designed to infer the primordial density field from input data on the spatial galaxy distribution. However, in contrast to RZA, BORG does not rely on the reversibility of the linear regime of structure formation, but instead uses many forward simulations in a Markov chain to compute the relative likelihood for each initial density field to be the source of the observed galaxy distribution. Though, in principle, BORG could infer jointly the cosmological parameters, the specific initial conditions used here were derived for a fixed cosmology, which sets a prior on the power spectrum of

the primordial density fluctuations. The observational data used to perform the inference is the 2M++ compilation (Lavaux & Hudson 2011). The 2M++ galaxies are photometrically selected from the 2MASS sample (2 Micron All Sky Survey; Skrutskie et al. 2006). The redshift information has been assembled from the 2MASS Redshift Survey for the brightest galaxies (Huchra et al. 2012), the SDSS DR7 main galaxy sample (Sloan Digital Sky Survey Data Release 7; Abazajian & Survey 2009), and the 6dFGRS (6 degree Field Galaxy Redshift Survey; Jones et al. 2009). BORG is able to account for selection effects; thus, there is no sharp boundary to the constrained volume, which in practice reaches out to a distance of about 300 Mpc from the observer. BORG also infers self-consistently the bias relation between the galaxy population and the underlying dark matter density field through the use of empirical relations. The details of the inference mechanism are given in Jasche & Lavaux (2019).

We note that, by design, BORG does not yield a unique result, but instead employs a Markov chain to represent plausible samples of the probability that each initial condition yields the present universe. Each iteration of the Markov chain corresponds to a different density field and an associated likelihood value. Even when the overall likelihood converges, the local density field is still subject to significant variations arising from the uncertainties in the properties of the galaxy population in each volume element ($\sim 3.9$ cMpc$^3$).

A different realization would to differences in the local LSS, owing to selection effects, truncation of small scale density fluctuations and uncertainty on galaxy formation processes. Fig. 2 shows the projected dark matter density in spheres of 50 Mpc for eight different realizations of the constraints with comparable likelihood. Overlaid on each panel are the positions of galaxies in the Local Universe, from the 2MASS survey (Huchra et al. 2012). Table 1 lists the masses and positions of the closest analogues to the Virgo and Fornax clusters in each of the eight BORG realizations. It is worth noting that different observational techniques can give different mass estimates for the same object. For Virgo, Simionescu et al. (2017) find $M_{200, crit} = 1.05 \pm 0.02 \times 10^{14}$M$_\odot$ from X-ray measurements while Kashibadze, Karachentsev & Karachentseva (2020) find $6.3 \pm 0.9 \times 10^{14}$M$_\odot$ from cluster member kinematics, while the Fornax mass estimate of Nasonova et al. (2011) is derived from the Hubble flow around the cluster. All eight BORG realizations produce similar large-scale features, including the presence of both Virgo and Fornax at their approximate positions, but there is noticeable variation between individual realizations. We selected the last iteration (9350) of the chain, where Virgo is at a distance of 18.1 Mpc with $M_{200, crit} = 2.3 \times 10^{14}$M$_\odot$, and Fornax is at a distance of 21.0 Mpc with $M_{200, crit} = 2.0 \times 10^{14}$M$_\odot$. A more complete description of the cosmography and a detailed analysis of the local environment beyond the LG, including the properties of several nearby galaxy clusters, will be presented in a forthcoming paper (McAlpine et al. in preparation).

## 2.3 Converting the BORG constraints into an octree function representation

The LSS constraints on the phase information produced by BORG are expressed as a $256^3$ element Gaussian white noise field. Taking the Fourier transform of this field, we produce an equivalent octree representation by imposing the Fourier modes provided by BORG as a set of linear constraints on an unconstrained white noise field generated by octree functions.

We do this following the method outlined in Salmon (1996). This procedure applied to Gaussian white noise fields is particularly simple. In essence, we create an unconstrained Gaussian white noise field in the octree representation, following section 5.2 of Jenkins
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Projected dark matter density in spheres of $r = 50$ Mpc, centred on the fiducial observer at $z = 0$, in labelled iterations of the Markov chain of the BORG reconstruction. Overlaid in white are the positions of galaxies in the 2MASS redshift survey (Huchra et al. 2012) between $r = 10$ and $50$ Mpc. The observed clusters and voids are consistently reproduced, but there are noticeable differences between individual iterations.</figcaption></figure>

<table><caption>**Table 1.** Masses ($M_{200}$, crit) and positions relative to the fiducial observer of the two closest analogues to the Virgo and Fornax clusters for the eight BORG realizations with comparable global likelihood, whose projected density is shown in Fig. 2.</caption><thead><tr><th rowspan="2">Index</th><th colspan="4">Virgo</th><th colspan="4">Fornax</th></tr><tr><th><i>M</i><sub>200</sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th><i>D</i><br>(Mpc)</th><th>RA<br>(°)</th><th>Dec.<br>(°)</th><th><i>M</i><sub>200</sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th><i>D</i><br>(Mpc)</th><th>RA<br>(°)</th><th>Dec.<br>(°)</th></tr></thead><tbody><tr><td>7600</td><td>2.7</td><td>17.5</td><td>181</td><td>19</td><td>0.3</td><td>19.0</td><td>36</td><td>-15</td></tr><tr><td>7850</td><td>1.3</td><td>18.0</td><td>185</td><td>16</td><td>1.1</td><td>21.3</td><td>42</td><td>-20</td></tr><tr><td>8100</td><td>0.7</td><td>16.4</td><td>191</td><td>29</td><td>0.7</td><td>17.6</td><td>45</td><td>-27</td></tr><tr><td>8350</td><td>1.8</td><td>19.9</td><td>188</td><td>13</td><td>1.3</td><td>21.6</td><td>42</td><td>-23</td></tr><tr><td>8600</td><td>2.0</td><td>20.8</td><td>186</td><td>13</td><td>0.9</td><td>18.9</td><td>45</td><td>-28</td></tr><tr><td>8850</td><td>4.3</td><td>20.5</td><td>193</td><td>23</td><td>0.5</td><td>20.0</td><td>43</td><td>-27</td></tr><tr><td>9100</td><td>2.9</td><td>18.2</td><td>189</td><td>17</td><td>1.5</td><td>19.5</td><td>42</td><td>-31</td></tr><tr><td>9350</td><td>2.3</td><td>18.1</td><td>191</td><td>14</td><td>2.0</td><td>21.0</td><td>38</td><td>-22</td></tr><tr><td>Obs.</td><td>1.05</td><td>16.5</td><td>187</td><td>10</td><td>1.2</td><td>19.6</td><td>9</td><td>-35</td></tr></tbody></table>

*Note.* The final row shows the observed values for the Virgo cluster (Mei et al. 2007; Simionescu et al. 2017) and Fornax cluster (Nasonova, de Freitas Pacheco & Karachentsev 2011), respectively.

(2013) and compute the Fourier amplitudes of this unconstrained field. We then adjust each of the corresponding Fourier amplitudes to exactly match the corresponding Fourier mode predicted by BORG, and then convert back into the octree representation.

This procedure is not unique as it can be applied to different sets of octree functions depending on the choice of how deep to go within the tree. For this paper, we chose to use the octree information down to level 21. At this level the information content of the octree functions for the simulation volume consists of $240^3$ degrees of freedom, which is a close match to the BORG constraints. As the number of degrees of freedom in the octree representation are slightly less than in BORG, we miss some of the BORG shortest wave Fourier modes as constraints. The BORG constraints that we use are not unique at these scales in any case. This particular choice means that in practice all the information in the initial unconstrained octree Gaussian white noise

field is overwritten and the constraints are effectively confined to definite levels of the octree. This leaves all octree functions at deeper levels of the tree unconstrained, and thus free to be randomly varied.

## 2.4 Random variations

In a second step, as shown in Fig. 1, we take a sample from BORG and supplement it with random phase information at level 22 and above, i.e. randomizing phases with $\lambda_{cut} = 3.24$ cMpc and smaller. In total, we create 60 000 variations of the initial density field constituting an ensemble simulation (Aragon-Calvo 2016). We call this ensemble the ‘3.2-cMpc set’, and label individual variants with index *i*. The two top rows of Fig. 3 show the projected dark matter density in spheres of radius 5 Mpc centred on the fiducial observers in four ‘3.2-cMpc’ variations. As all variations share the same relevant phase information
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top two rows: projected dark matter density in spheres of $r = 5$ Mpc, centred on the fiducial observer, in four of the 60 000 3.2-cMpc random variations of constraint 9350. Due to the shared large-scale modes, all volumes show a similar filamentary structure. However, the identity of individual halos, determined on smaller scales, changes between volumes. Bottom two rows: spheres of $r = 2$ Mpc, centred on the fiducial observer, in four of the 0.8-cMpc variations, using the same constraints and the same 3.2-cMpc variation, $i = 24$ 493. Now, all volumes contain an LG analogue of similar mass and orientation, but varying separation.</figcaption></figure>

in the constrained region up to level 21 ($\lambda_{cut} = 6.48$ cMpc), all simulations in this sample have a similar intermediate scale density, including a planar structure and several visible filaments. However, randomizing the phase information on levels 22 and above yields different populations of halos. We describe the properties of the LG analogues found among the ‘3.2-cMpc’ set in Section 3.

As shown in Sawala et al. (2021), specifying the phase information up to and including octree level 22 ($\lambda_{cut} = 3.24$ cMpc) determines the existence of individual halos down to $M_{200} \sim 10^{12} M_{\odot}$, i.e. halos similar in mass to the MW and M31. After identifying pairs of such halos among the simulations of the 3.2-cMpc set,

we create smaller scale variations of promising candidates by independently randomizing the phase information at levels 24 ($\lambda_{cut} = 0.81$ cMpc) and 25 in a small region, as shown on the right-hand panel of Fig. 1. We refer to these sets as the 0.8-cMpc sets, and label individual variants with an additional index, $j$. The effect of variations at 0.8 cMpc can be seen in the bottom two rows of Fig. 3. Here, an LG analogue consisting of a pair of halos of $\sim 10^{12} M_{\odot}$ is nearly always present, but its properties vary. We describe the results of these simulation in Section 4. Still smaller scale variations, up to level 28 ($\lambda_{cut} = 0.05$ cMpc), are described in Section 5.

## 2.5 Numerical Setup

All simulations were performed with the GADGET-3 code, a variant of the publicly available code GADGET-2 (Springel 2005). To reduce the computational cost, before running the full set of 60 000 3.2-cMpc variants, we first ran a prototype set of 100, refining a cubic region of $20^3$ cMpc$^3$ around the origin of the fiducial observer within the constraints. Because of the partially identical phase information, the variation in the Lagrangian volumes is limited. In the next step, we computed the union of all particles within $r = 8$ Mpc of the fiducial observer at $z = 0$ in all 100 variants. This defined the refinement region (or ‘mask’), ensuring that the region of interest ($r = 5$ Mpc) is always included in the high-resolution region, while limiting computational cost. Extracting this region out of a $20^3$ cMpc$^3$ cube sampled with $80^3$ particles resulted in only $\sim 52^3$ particles required in the high resolution region. Each simulation required 4.4 core-hours on the COSMA-6 computer.

The mass resolution of the 3.2-cMpc variants is $2.4 \times 10^9 M_{\odot}$, i.e. a typical MW or M31 halo of $10^{12} M_{\odot}$, is resolved with $\sim 420$ particles. We used a $480^3$ FFT for computing the initial particle displacements inside the high-resolution region. The other variants were performed at twice higher spatial and eight times higher mass resolution, using masks based on the 3.2-cMpc simulations. In all simulations, the particles that represent the LSS and the remainder of the periodic volume outside the high-resolution region, progressively increase in mass with increasing distance. We tested that the results do not change qualitatively when using higher particle resolution or finer FFTs.

## 2.6 Identification and characterization of LG candidates

The distance between the MW and M31 has been measured using *Hipparcos* and *HST* observations of red clump stars as $d = 784_{-17}^{+13}$ kpc (Stanek & Garnavich 1998), using TRGB stars as $d = 785 \pm 25$ kpc (McConnachie et al. 2005), using Cepheids as $d = 765 \pm 28$ kpc (Riess, Fliri & Valls-Gabaud 2012). The galactocentric radial velocity has been accurately measured, $v_r = -109.3 \pm 4.4$ km s$^{-1}$ (van der Marel et al. 2012). The transverse velocity is more controversial. The only direct measurement of proper motions with *HST* (van der Marel et al. 2012) give $v_t = 17 \pm 17$ km s$^{-1}$, while van der Marel et al. (2019) give a combined estimate from Gaia DR2 and *HST* data of $v_t = 57_{-31}^{+35}$ km s$^{-1}$. Both of these imply a highly radial orbit. Conversely, indirect measurements using satellite kinematics give far higher values: $v_t = 164.4 \pm 61.8$ km s$^{-1}$ (Salomon et al. 2016), and a much lower eccentricity. In this work, we adopt a distance of $d = 770$ kpc, a radial velocity of $v_r = -109$ km s$^{-1}$, and a transverse velocity of $v_t < 40$ km s$^{-1}$. Of course, the methods and results described here remain valid for different parameters values.
In our simulations, halos and self-bound subhalos<fnref n="1" /> were identified in the $z=0$ snapshots using the FoF and SUBFIND algorithms. For identifying LG analogues, we consider all possible pairs of self-bound halos, each in the mass range $5 \times 10^{11}-5 \times 10^{12} M_\odot$, whose centres of potential are separated by 0.5–2 Mpc. We also require that there be no other halo more massive than the smaller of the pair within 2 Mpc.

In previous works, both those that studied the LG in arbitrary environments (e.g. Sawala et al. 2016a) as well as those including constraints (Carlesi et al. 2016), the orientation of the LG with regards to the LSS was considered arbitrary. However, in constrained simulations, we argue that reproducing the correct orientation of the LG with respect to the observed LSS is not a mere accessory. If the distinct features observed in the LG are linked to the local LSS, and in particular, if the anisotropy of the LG’s environment plays any role in its evolution, then the LG’s orientation with respect to this environment would appear to be crucial. And, as we will discuss in Section 4.1, the orientation is also set on fairly large scales.

Our simulation axes are aligned with the equatorial coordinate axes, placing the observed local structures at their observed locations. To parametrize the orientation of the LG analogues relative to the LSS, we calculate in turn the vectors pointing from the centre of each halo to the other, and quantify the angles, $\delta_{1,2}$, by which they differ from the true position of M31 relative to the MW. We then designate the halos as the ‘MW’ and ‘M31’ analogues, respectively, so that $\delta = \min(\delta_{1,2})$; $\delta$ thus ranges from 0° for an LG analogue with the correct orientation to 90° for an LG analogue whose orientation is perpendicular to the observed one. Having assigned the MW and M31 counterparts, we parametrize the mass ratio as $M31/MW = M_{M31}/M_{MW}$. Similarly to the total mass, the mass ratio is also uncertain, with most estimates within a factor of 2 of unity, and M31 likely to be the more massive of the pair (Li & White 2008; Phelps, Nusser & Desjacques 2013; Sofue 2015).

We also require that the position of the MW be within 5 Mpc of the fiducial observer defined by the cosmological constraints. Subsequently, as outlined in Section 3.1, we apply additional limits on the mass, mass ratio, separation, velocities, and orientation of the LG candidates.

## 3 VARIATIONS AT $\lambda_{\text{CUT}} = 3.2$ CMPC

A simulation where the density field is constrained down to only ~4 Mpc does not generally yield an LG at its centre. As explained in Section 2.4, to complement the phase information constrained by the local LSS, in a first step, we generated a set of 60 000 simulations, randomly varying the phase information at scales of $\lambda_{\text{cut}} \sim 3.2$ cMpc and smaller, which is levels of 22 and above.

### 3.1 Number of LG analogues

Searching the central $r=5$ Mpc sphere in all 60 000 simulations of the 3.2-cMpc set for LG analogues yielded a total of 48 323 halo pairs (hereafter labelled ‘all’) with an LG mass in the range $(1-8) \times 10^{12} M_\odot$, a separation of 0.5–2 Mpc, and fulfilling the condition that no other halo more massive than the lower-mass member be within 2 Mpc of the pair mid-point. As shown in Table 2, restricting

<table><caption>**Table 2.** Selection criteria and numbers of LG analogues among the total of 48 323 candidates from the 3.2-cMpc set.</caption><thead><tr><th></th><th>'Loose'</th><th>'Intermediate'</th><th>'Strict'</th></tr></thead><tbody><tr><td>M (10<sup>12</sup> M<sub>☉</sub>)</td><td>1.2, ..., 6.0<br>(46 054)</td><td>1.5, ..., 5.0<br>(39 090)</td><td>2.0, ..., 4.0<br>(24 108)</td></tr><tr><td>M31/MW</td><td>2/5, ..., 5<br>(40 490)</td><td>2/3, ..., 3<br>(28 630)</td><td>1, ..., 2<br>(14 162)</td></tr><tr><td>d (Mpc)</td><td>0.5, ..., 1.5<br>(34 693)</td><td>0.6, ..., 1.0<br>(15 356)</td><td>0.74, ..., 0.80<br>(1665)</td></tr><tr><td>v<sub>r</sub>(km s<sup>-1</sup>)</td><td>-200, ..., 0<br>(11 414)</td><td>-150, ..., -50<br>(4318)</td><td>-109, ..., -99<br>(597)</td></tr><tr><td>v<sub>t</sub>(km s<sup>-1</sup>)</td><td>&lt; 150<br>(38 885)</td><td>&lt; 100<br>(23 904)</td><td>&lt; 40<br>(10 792)</td></tr><tr><td>N(M,M31/MW, d, v<sub>r</sub>, v<sub>t</sub>)</td><td>6385</td><td>489</td><td>1</td></tr><tr><td>δ</td><td>&lt;45°</td><td>&lt; 30°</td><td>&lt; 15°</td></tr><tr><td>N(M,M31/MW, d, v<sub>r</sub>, v<sub>t</sub>, δ)</td><td>2309</td><td>82</td><td>0</td></tr></tbody></table>

*Notes.* M is the sum of the MW and M31 subhalo masses, d is the distance between the centres of potential, $v_r$ and $v_t$ are the relative velocities, and $\delta$ is the error in the position on the sky of M31. Numbers in brackets indicate the number of LG analogues from the 3.2-cMpc variations satisfying the corresponding criterion. N indicates the number that simultaneously satisfy all criteria in each column, either excluding or including the constraint on the orientation.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Fraction of simulation volumes from the 3.2-cMpc set with $N_{\text{LG}}$ LG analogues selected by different criteria. With the ‘loose’ criteria, 9.8 per cent of variants contain one LG analogue, 0.4 per cent contain two, and four volumes (0.0067 per cent) contain three LG pairs. Stricter limits on mass, mass ratio, separation, and velocity reduce the numbers, but even more so, the multiplicity: With the intermediate criteria, only two variants (0.0033 per cent) have two LGs and none have more than two.</figcaption></figure>

the selection criteria closer to the observed values has the effect of significantly reducing the number of LG analogues: our ‘loose’ criteria for mass, mass ratio, separation, radial velocity, and tangential velocity yield 6385 LG analogues, and only one LG analogue fulfils the strict ‘criteria’.

Once we also include constraints on the orientation, the numbers are reduced even further, and when we limit the orientation of the LG relative to the LSS to within 15°, none fulfils all strict criteria.

### 3.2 Multiplicity

With the broadest selection criteria, labelled ‘all’, a single simulation often contains more than one LG analogue within the central $r=5$ Mpc volume. Fig. 4 shows the number of variants from the 3.2-cMpc set containing N LG analogues for different selection criteria.

<footnote n="1">¹In the following, we use the term ‘halo’ when referring to a self-bound structure identified by the SUBFIND algorithm (Springel et al. 2001), and we consider the mass of the self-bound structure as the mass of the halo. The LG mass is defined as the sum of the two subhalo masses.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Duplicates of an LG analogue from the 3.2-cMpc variations. The top left-hand panel shows the projected DM density around the LG for $i = 622$, highlighting in red and yellow particles within 250 kpc of the centres of the MW and M31 analogue halos. The other three panels show the dark matter density in three other realizations of the 3.2-cMpc set, at the same $z = 0$ coordinates, and with the same particles highlighted. All four LG analogues share more than 50 per cent of particle IDs, i.e. they have formed from the same Lagrangian volume, and also have similar mass and orientation.</figcaption></figure>

As the selection criteria are tightened, the number of variants that contain multiple LG analogues decreases sharply, and more steeply than $N^{-2}$. With the ‘loose’ constraints, 9.8 per cent of variants have one LG, but only 0.4 per cent contain two and none contain more than three. With intermediate criteria, while there are 489 LGs, the chance of finding two in the same volume is ~1: 30 000.

More interestingly, we not only find multiple different LG analogues within the same $r = 5$ Mpc volume, we also find many instances of ostensibly the *same* LG forming in multiple volumes among the 3.2-cMpc variations. For the purpose of finding such duplicates, we consider two LGs in different simulations to be the same object if their centres of mass differ by less than 300 kpc, their total masses differ by less than 30 per cent, and their orientations are aligned to within 30°. We find that only ~10 per cent of LGs from either the loose or intermediate set have *no* duplicates, while almost ~30 per cent of LGs are found more than 10 times.

In Fig. 5, we show four ‘copies’ of an LG analogue from the 3.2-cMpc set. All four panels are centred on the same coordinates, and all four contain a pair of halos with similar orientation. As described in Section 2.5, all 3.2-cMpc variants use the same particle load, so that particles with the same IDs have identical (unperturbed) Lagrangian coordinates in all simulations from this set. Highlighted in red and yellow on all four panels are the positions of particles within 250 kpc of the centres of the MW and M31 analogues shown in the top left-hand panel. In each of the other three simulations, a large fraction of the same particles are also found in the two halos. In general, we find that copies of LGs share a large fraction of particles, indicating that they originate from largely overlapping Lagrangian volumes.

The fact that most of our LG analogues exist more than once suggest that, within the 60 000 random realizations of the 3.2-cMpc set, we are likely to have exhausted a large fraction of

the possible, truly distinct LG analogues that can form given this particular set of constraints. As we will discuss in the conclusions (Section 6), this has important implications for our future ability to test cosmological models. In the following sections, we make use of the fact that we can generate several variations of the same LG by deliberately introducing small-scale random variations to promising LG candidates identified among the 3.2-cMpc variants.

### 3.3 Quality of LG analogues

Fig. 6 illustrates the quality of individual LG analogues for all 8086 3.2-cMpc variations that satisfy the ‘loose’ criteria irrespective of orientation, as given in Table 2. On the x-axis, we show the index, $i$, while on the y-axis, we quantify the agreement with the observation in terms of the logarithm of the mass ratio mass (M31/MW, orange), the total mass (M, purple), orientation ($\delta$, brown), tangential velocity ($v_t$, green), radial velocity ($v_r$, red), and separation ($d$, blue). For each variable, the set of solid black lines correspond to our strict criteria (see Table 2), either the observational uncertainty, or in the case of the orientation, where the observational uncertainty is negligible, an acceptable error of $\delta = 15°$. The grey lines denote the intermediate and loose criteria.

The vertical positions of individual points of corresponding colour are normalized relative to the strict criteria, and show

$$ \Delta = \frac{x_{n,i} - x_n}{\sigma_n}, \qquad (1) $$

where $x_{n,i}$ is the value of variable $x_n$ for simulation $i$, $x_n$ is its observational value, and $\sigma_n$ is the observational uncertainty or acceptance region. For clarity, different variables are offset vertically. The grey points at the top of the figure show the sum of the individual errors that lie outside the observational bounds, i.e.

$$ \Sigma_i = \sum_n \left( \min \left( 1, \frac{|x_{n,i} - x_n|}{\sigma_n} \right) \right) - n. \qquad (2) $$

An LG candidate that lies within the bounds for all five variables would have a value of 0 and fall on the thick dashed line on the plot. The fact that all grey points lie above this line indicates that $\Sigma > 0$ for all candidates, i.e. all candidates from the 3.2-cMpc set are outside the bounds for at least one observable.

### 4 VARIATIONS AT $\lambda_{CUT} = 0.8$ CMPC

As discussed in the previous section, our set of 60 000 random variations produced tens of thousands of LG analogues within the constraints, and indeed, resulted in many variants of ostensibly the *same* objects, but it did not yield any that matches all of the strict criteria for the MW–M31 pair, and only one that matches the strict constraints neglecting the orientation. While it would be possible simply to create more tries, a back-of-the-envelope estimates of the number of tries required suggest that this would be prohibitively expensive. For example, dividing the number of tries that produced one LG analogue according to our ‘strict’ criteria irrespective of orientation by the fraction of the sky represented by an area of radius 15°, or extrapolating from the fraction of LGs that match the intermediate criteria with and without orientation, suggest that it would require over a million simulations to find a realization that satisfies all constraints simultaneously. Importantly, this does not imply that the LG, as defined here, is a particularly rare system, but simply reflects the fact that any system measured precisely enough is exceedingly unlikely to emerge from randomized initial conditions, making a single stage of randomization impractical.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Quality of LG analogues for 3.2-cMpc variations. Each vertical set of points represents an LG analogue from a simulation with index, *i*. Sets of horizontal stripes indicate the logarithm of the mass ratio (M31/MW, orange), the mass (M, purple), angular offset (δ, brown), tangential velocity (v_t, green), radial velocity (v_r, red), and separation (d, blue), respectively. For each observable, dashed lines indicate either the 1σ observational uncertainty, or the accepted error for our 'strict' criteria. Grey points on top show Σ, the total error outside the accepted range for all observables. While many variants fulfil the criteria for individual observables and some fulfil multiple criteria, no points fall on the thick dashed line, i.e. none of the 3.2-cMpc variants matches all six criteria simultaneously.</figcaption></figure>

Instead, we select the most promising candidates for which a change of smaller scales in the initial density field might produce variations that fulfil the constraints. For these, we keep the phase information fixed for two further levels above the constraints, 22 and 23, and independently vary levels 24 and above (corresponding to scales with $\lambda_{cut} < 0.8$ cMpc). As discussed in Section 2.4, we expect that a majority of these variants will contain the same individual ~$10^{12}$ M☉ halos, but with a scatter in their masses, positions, and velocities.

## 4.1 Characteristics of the LG

To select suitable candidates for refinement, it is important to determine first which properties of the LG are set on large scales in the primordial density field, and which ones vary when smaller scales are modified. As discussed in Section 2.6, observationally, the LG is characterized primarily and most precisely by the separation of the Milky Way and Andromeda, their radial velocity, their transverse velocity, and the LG's orientation relative to the LSS. It may seem natural to select those candidates among the 3.2-cMpc variations that come closest to matching all of these criteria. However, as shown in Fig. 7, we find that when we introduce variations on smaller scales, some of these observables still show a great amount of scatter, while others appear to be already set on the scales that determined the existence of the MW and M31 halos themselves. In total, we selected nine LG variants among the 3.2-cMpc sets, and created 1000 further variations of each, randomizing only phase information at scales of $\lambda_{cut} = 0.81$ cMpc and smaller.

The top panel of Fig. 8 shows the correlation of radial velocity and separation for the 3.2-cMpc variants (in grey), and for the nine sets of 0.8-cMpc variants, coloured according to index *i*. As expected, not all combinations are equally likely, even for the 3.2-cMpc set: More distant LG pairs are typically receding from one another, closer pairs

are typically approaching. However, the scatter in the individual 0.8-cMpc sets is remarkably reduced: while distances and radial velocities independently still show a large scatter among the variants in each set, their combinations only span a narrow region. As we discuss below, this is evidence that the 0.8-cMpc variations result in LGs of similar orbits, but differing orbital phases.

The centre panel in Fig. 8 shows the radial velocity and the magnitude of the transverse velocity. Compared to halo pairs with the same mass, separation, and radial velocity, the small transverse velocity of the LG is unusual (Li & White 2008; Fattahi et al. 2015). We find that for each set of 0.8-cMpc variations, the scatter in transverse velocity is narrower than the scatter in radial velocity, and this ratio is greater for the LG analogues on more radial orbits. Again, this is expected if each set of 0.8-cMpc variants produces LG analogues of similar, mostly radial orbits whose eccentricities do not vary much.

### 4.1.1 Integrals of motion

The distribution of distances and velocities along orbital trajectories suggests that we should consider as more fundamental the quantities that remain constant in orbital motion. If we consider the LG as a two-body system comprised of the MW and M31, these are the total mass, $M = M_{MW} + M_{M31}$, the specific orbital energy,

$$e_{\text{tot}} = \frac{1}{2} v^2 - \frac{GM}{r}, \qquad (3)$$

the specific orbital angular momentum,

$$l_{\text{orb}} = r \times v, \qquad (4)$$

and the reduced Laplace-Runge-Lenz vector, which defines the eccentricity and the orientation of the orbit,

$$a_{\text{LRL}} = \frac{v \times l_{\text{orb}}}{GM} - \hat{r}, \qquad (5)$$
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Probability density functions of mass, radial velocity, transverse velocity and separation (top row) and mass ratio, specific orbital angular momentum, specific orbital energy, and orbital spin parameter (bottom row), for the 3.2-cMpc variants (grey dashed lines), and for the nine sets of 0.8-cMpc variations, independently normalized and coloured according to the index, *i*. The quantity, *N*, gives the number of LG analogues matching the ‘loose’ constraints. For each set of 0.8-cMpc variants, the scatter is smaller than in the 3.2-cMpc variants in all observables, but while some observables (such as radial velocity and distance) still have comparable scatter, for other variables (e.g. mass, energy, and orbital angular momentum), the scatter is greatly reduced and depends on the set.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Radial velocity versus separation (top panel), transverse velocity (middle panel), and specific orbital angular momentum (bottom panel), for LG analogues in the 3.2-cMpc sample (grey dots), and in the nine sets of 0.8-cMpc variations, coloured according to *i*. Grey bands show our ‘strict’ criteria for *d*, $v_r$, $v_t$, and the corresponding specific orbital angular momentum, $I_{\text{orb}}$. There is still considerable scatter in both $v_r$ and *d* among the variants of each the 0.8-cMpc set, but they only occupy a narrow region of phase space, corresponding to LG orbits of similar energy but different orbital phase.</figcaption></figure>

where *v* is the velocity of M31 relative to the MW, *r* is the position of M31 relative to the MW, *r* = |*r*| is the distance between the two, and $\hat{r} = r/r$ is the unit vector from the MW to M31.

From Fig. 7 and the bottom panel of Fig. 8, it can be seen that the 0.8-cMpc variants show little scatter in specific angular momentum compared to the 3.2-cMpc variants. Fig. 9 shows the distributions of

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Specific potential and specific kinetic energy of all LG analogues in the 3.2-cMpc sample (grey dots), and of the LG analogues from the 0.8-cMpc variations, coloured according to =. Diagonal lines show constant total energy, and the black line marks zero energy orbits. Light, medium, and dark blue lines show the energy for LG analogues with *d* = 770 kpc, $v_\infty$ = -109 km s⁻¹, $v_t$ = 40 km s⁻¹, and with total masses of 2 × 10¹², 3 × 10¹², and 4 × 10¹² M☉, respectively. LGs in the same set appear to have similar total energy, although there is a slight trend for closer pairs to have more negative total energy.</figcaption></figure>

the 3.2- and 0.8-cMpc variants in the kinetic and potential energy plane. The LG analogues from a particular set, while greatly differing in both kinetic and potential energy, align along lines of similar total energy. This suggests that variances in radial velocity and separation reflect, to a large extent, a variance in orbital phase.

In Fig. 10, we show the positions on the sky of M31 in all LG analogues from the 3.2-cMpc sample (in grey), and in the nine sets of 0.8-cMpc variations, in a Mollweide projection. The assignment of the MW and M31 analogues by the LG's orientation limits the position of the M31 analogue to 90° around the true position of M31 (indicated by the black cross). Due to the common constraints,
<figure><img src="image_12.png" /></figure>

**Figure 10.** Mollweide projection of the position of M31 on the sky, for all LG analogues in the 3.2-cMpc sample (grey dots), and for the LG analogues from the 0.8-cMpc variations, coloured according to index = as in Figs 8 and 9. The black and grey dashed circles indicate 15° and 30° around the true position of M31 observed from the MW, marked by the black x. Because we assign the 'MW' and 'M31' analogues by orientation, the maximum difference from the true orientation is 90°. Points from the 0.8-cMpc variations cluster closely together and close to the true orientation.

the orientations among the 3.2-cMpc set is not uniform; some orientations are more likely than others. In each set of 0.8-cMpc variants, the scatter is significantly reduced. This indicates that the orientation of the orbital plane, and $a_{LRL}$, are also conserved along the orbit.

For all 0.8-cMpc variants, the orientation close to the observed value, and all orientations within a single set of 0.8-cMpc variants are similar to one another. This result also echoes the findings from the multiplicity of LG analogues among the 3.2-cMpc variations described in Section 3.2: When we find multiple pairs of halos of similar mass and in a similar position, their orientations also tend to be similar.

It seems clear that these constants of motion are set on larger scales than the canonical observables such as separation and velocities, and while they are more complicated to observe, it seems more natural to consider them as the defining properties of the LG. In Section 5, we will consider still smaller scale variations to determine the scales that set the phase of the orbit, but as we will see in the next section, the degree to which the LG is a true two-body system also affects the scatter in these integrals of motion.

## 4.2 Presence of massive third objects

In addition to M31 and the Milky Way, the LG contains two galaxies believed to be residing in halos of at least $10^{11}$ M$_{\odot}$: M33, at a distance of $203 \pm 37$ kpc from M31 (van der Marel et al. 2012), and the LMC at a distance of $\sim 50$ kpc (van der Marel et al. 2002) from the MW centre. With stellar masses of $(3-5) \times 10^9$ M$_{\odot}$, abundance matching places both in dark matter halos of $\sim 10^{11}$ M$_{\odot}$. H I kinematics further suggests that M33 could have a total mass as high as $5 \times 10^{11}$ M$_{\odot}$ if it resides in an NFW-like DM halo (Kam et al. 2017). For the LMC mass, timing constraints give a total mass of $2.5 \times 10^{11}$ M$_{\odot}$ (Peñarrubia et al. 2015); analogues by stellar mass in the EAGLE simulations indicate $M_{200} = 3.4_{-1.2}^{+1.8} \times 10^{11}$ M$_{\odot}$ (Shao et al. 2018), while its effect on the Orphan stream in the Milky Way suggests a mass of $1.38_{-0.24}^{+0.27} \times 10^{11}$ M$_{\odot}$ (Erkal et al. 2019).

In Fig. 11, we show histograms of the number of M33 or LMC analogues, which we define to be in the mass range (1–3) $\times 10^{11}$ M$_{\odot}$ within 1 and 2 Mpc of the LG centre, for the LGs from 3.2-cMpc variants, and for each of the nine sets of 0.8-cMpc variants.

While the nine sets of 0.8-cMpc variants collectively span a large range in the number of LMC/M33 analogues, we find that some sets

<figure><img src="image_13.png" /></figure>

**Figure 11.** Distribution of the number of objects in the mass range (1–3) $\times 10^{11}$ M$_{\odot}$ within 1 or 2 Mpc from the centre of the LG analogue, for the set of 3.2-cMpc variations (grey area), and for the nine sets of 0.8-cMpc variants (coloured lines). The average number of massive satellites in the 0.8-cMpc sets varies, but the scatter within the each set sets is smaller than that of the 3.2-cMpc set.

$(i = 12, 885, 21, 826, 36, 993, 43, 595)$ have no LMC/M33 counterparts within 1 Mpc in a large majority of cases, while other sets $(i = 24, 943, 37, 497, 41, 800)$ typically have one or more. To a large extent, the presence of such halos is already determined by the scales shared within each set.

As can be inferred from Figs 8 and 9, the number of massive additional halos appears to correlate with the scatter in the integrals of motion of a two-body orbit, as discussed in Section 4.1.1. For example, the sets of 0.8-cMpc variants with $i = 21, 826$ and $43, 595$ (shown by green and cyan lines or symbols), which typically have no LMC/M33 analogues within 1 Mpc, show very little scatter in total energy, mass, or angular momentum, while sets with a higher number of LMC/M33 analogues, such as $i = 37, 497$ and $41, 800$ (pink and olive lines and symbols), have a much greater scatter in the same quantities.

## 4.3 An LG from 0.8-cMpc variations

Fig. 12 shows the quality of the LG in the nine sets of 0.8-cMpc variants, analogous to Fig. 6 discussed in Section 3.3. We find that the typical mass error, and the error in transverse velocity, is much reduced compared to the 3.2-cMpc variations shown in Fig. 6. The error in orientation clusters close to 0°. Reflecting the discussion in Section 4.1, we find that the scatter in distance or radial velocity, independently, are not substantially smaller than in the 3.2-cMpc set. However, because of the total energy constraint, $d$ and $v_r$ are not independent for LG analogues in the 0.8-cMpc sets. For some sets, this means that no LG analogue matches the observations of both simultaneously, while for others, it significantly increases the fraction of LG analogues that match both observables simultaneously.

If it is possible for an LG to form inside the LSS constraints, independently randomizing all levels simultaneously would eventually yield an LG (the 0.8-cMpc variants are, strictly speaking, a proper subset of a (much larger) set of possible 3.2-cMpc variants). The 0.8-cMpc variants bifurcate: Because the scatter in the observables is much reduced, for almost all possible sets of 0.8-cMpc variations, it would be impossible to form an LG that matches the observations. Even among the set of nine promising 3.2-cMpc variants we selected, we find examples where 0.8-cMpc modifications will never lead to a
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Quality plot of 0.8-cMpc variations, *j*, of 3.2-cMpc variations, *i*, analogous to Fig. 6. There is very little scatter in total mass, angle or, in most cases, mass ratio, less scatter for $v_t$, but still a great amount of scatter in $v_r$ and *d*. The total error, $\Sigma$, is significantly reduced, and the remaining error is now dominated by $v_r$ and *d*. In some cases, the (largely conserved) energy of the orbit (e.g. *i* = 43 595) or angular momentum (e.g. *i* = 12 885) imply that no correct combination of *d* and $v_r$ or correct value of $v_t$ is possible. However, the set with *i* = 24 943 contains one LG analogue (*j* = 366) that matches all 'strict' constraints.</figcaption></figure>

formation of an LG, either because the transverse velocity is always
too great (i = 12 885), or because the total energy is never negative
enough (i = 21 826 and 43 595).

However, set *i* = 24 943 contains one LG analogue [*i* = 24 943, *j* = 366] that matches all observational constraints: *M*<sub>LG</sub> = 2.8 × 10<sup>12</sup> M<sub>☉</sub>, M31/MW = 1.07, *d* = 769 kpc, *v*<sub>r</sub> = −100 km s<sup>−1</sup>, *v*<sub>t</sub> = 37 km s<sup>−1</sup>, and *δ* = 9.4. As discussed in Section 4.2, similar to *i* = 28 471 and 37 497, variants from this set typically contain several LMC/M33 mass objects, and the scatter in the transverse velocity and the total energy of the two-body orbit is larger than for sets that contain none.

This does not necessarily imply that the presence of LMC/M33
analogues is required to form an LG with the observed parameters
(although it is necessary in our particular LSS constraints), but we
can conclude that the presence of the LMC and M33 significantly
affects the kinematics of the MW–M31 orbit.

5 SMALLER SCALE VARIATIONS

By varying higher levels of the octree, we can create even smaller
scale variations of existing LG analogues. From the 0.8-cMpc sets,
we selected four candidates, [*i* = 21 868, *j* = 1], [*i* = 41 800, *j* = 1], [*i* = 43 595, *j* = 1], and [*i* = 24 943, *j* = 366], the latter of which yielded
an LG that matches the observations. For each one, we created four
new sets of 50 variations each, by varying the phase information at
and above levels 25, 26, 27, and 28, corresponding to $\lambda_{cut} = 0.4, 0.2,
0.1$, and $0.05$-cMpc, respectively.

In Fig. 13, we illustrate the effect for [*i* = 24 943, *j* = 366]. The four columns show the effects of applying four different shifts in PANPHASIA at levels 25, 26, 27, and 28 of the octree (top to bottom). All variations shown contain an LG analogue, and as expected, the variation among them decreases from top to bottom with decreasing scale of the variation in the density field. At *L* ≥ 26, even some small individual substructures can be easily identified across the different variants.

<figure><img src="image_15.png" /><figcaption>**Figure 13.** The effect of varying the phase information on levels 25 ($\lambda_{\text{cut}}$ = 0.4 cMpc, top row), 26 ($\lambda_{\text{cut}}$ = 0.2 cMpc, second row), 27 ($\lambda_{\text{cut}}$ = 0.1 cMpc, third row), and 28 ($\lambda_{\text{cut}}$ = 0.05 cMpc, bottom row) on the projected density of LG candidate [*i* = 24 943, *j* = 366].</figcaption></figure>

**5.1 Constraining the MW–M31 orbital phase**

Fig. 14 shows the probability density functions of the mass, radial velocity, transverse velocity, and separation for 50 variations each at levels 25, 26, 27, and 28 ($\lambda_{\text{cut}}$ = 0.4, 0.2, 0.1, and 0.05, respectively) of the four LG variants from the 0.8-cMpc set listed above. The scatter in all observables decreases as the scale of the variations is decreased.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Probability density functions of LG observables; mass, radial velocity, transverse velocity, and separation for smaller scale variations of four different LG variants, [$i = 21$ 826, $j = 1$] (greens), [$i = 41$ 800, $j = 1$] (oranges), [$i = 43$ 595, $j = 1$] (blues), and [$i = 24$ 943, $j = 366$] (greys), analogous to Fig. 7. From light to dark, four different shades show variations at levels 25, 26, 27, and 28 ($\lambda_{cut} = 0.4, 0.2, 0.1$, and 0.05 cMpc, respectively). Reducing the scale of the variations in the initial conditions reduces the scatter in all observed variables. Note that, as the PDFs are independently normalized for each variable, the height of the curves can be compared for different variations on the same panel, but not across panels.</figcaption></figure>

It is worth noting that at $L = 25$, not all simulations have yielded an LG analogue that matches the ‘loose’ criteria. Variations at $L = 25$ of [$i = 21$ 826, $j = 1$] contain an LG analogue in 47/50 cases. For [$i = 41$ 800, $j = 1$], [$i = 24$ 943, $j = 366$], and [$i = 43$ 595, $j = 1$], the rate is 49/50, 44/50, and 41/50, respectively. Higher level variations of [$i = 21$ 826, $j = 1$], [$i = 41$ 800, $j = 1$], and [$i = 24$ 943, $j = 366$] always contain an LG that fits at least the loose criteria, but for [$i = 43$ 595, $j = 1$], the number remains around 40–45: The separation is very close to 500 kpc, and a fraction of LG analogues fall below the threshold.

The extent to which individual variables vary at each level depends on their values. For example, the scatter in separation, *d*, is largest for the sets that have the *lowest* separation. Conversely, the scatter in radial velocity is largest for the sets that have the most negative velocity and smallest for velocities close to zero. This can be readily understood if, as explained in Section 4.1, the small-scale variations are, to a large extent, perturbations to the phase of the LG’s orbit, whose integrals of motion have already been defined on larger scales. LG variants with large separation are close to apocentre, where the orbital velocity is low. In this case, a perturbation to the phase causes only a small change in position and velocity. On the other hand, LG variants with small separation are closer to pericentre, and the same perturbation to the phase causes a larger change in position and velocity.

To determine which scales in the primordial density field of the Local Universe were responsible for setting the phase of the LG’s orbit, we are thus primarily interested in perturbations close to the observed values, i.e. $d \sim 770$ kpc, $v_r \sim -110$ km s$^{-1}$. In Table 3, we list the sensitivity of different parameters to the variations of LG candidate [$i = 24$ 943, $j = 366$].

For all scales shown, the mass of the LG varies by $\sim 10$ per cent or less, well below the observational uncertainty. At $L = 25$, the scatter in separation and radial velocity is $\sim 200$ kpc and $\sim 40$ km s$^{-1}$,

<table><caption>**Table 3.** Sensitivity of the LG’s observables to small-scale variations of the initial density for candidate [$i = 24$ 943, $j = 366$].</caption><thead><tr><th>L</th><th>λ<sub>cut</sub><br>(cMpc)</th><th>M ± σ<sub>M</sub><br>(10<sup>12</sup>M<sub>☉</sub>)</th><th>d ± σ<sub>d</sub><br>(kpc)</th><th>v<sub>r</sub> ± σ<sub>v<sub>r</sub></sub><br>(km s<sup>-1</sup>)</th><th>v<sub>t</sub> ± σ<sub>v<sub>t</sub></sub><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>25</td><td>0.4</td><td>2.53 ± 0.28</td><td>1044 ± 205</td><td>-57 ± 36</td><td>30 ± 16</td></tr><tr><td>26</td><td>0.2</td><td>2.88 ± 0.12</td><td>857 ± 68</td><td>-67 ± 14</td><td>9 ± 11</td></tr><tr><td>27</td><td>0.1</td><td>2.88 ± 0.05</td><td>766 ± 29</td><td>-96 ± 8</td><td>37 ± 9</td></tr><tr><td>28</td><td>0.05</td><td>2.85 ± 0.04</td><td>774 ± 12</td><td>-101 ± 3</td><td>33 ± 4</td></tr></tbody></table>

*Notes.* *L* is the level of the octree on which coefficients are randomized, and $\lambda_{cut}$ is the corresponding cut-off scale for the $\Lambda$CDM power spectrum. *M*, *d*, *v*<sub>r</sub>, and *v*<sub>t</sub> are the median values of the LG mass, separation, radial velocity, and tangential velocity, respectively, while $\sigma_M$, $\sigma_d$, $\sigma_{v_r}$, and $\sigma_{v_t}$ are the 1σ scatter in the corresponding observables.

respectively. It is worth noting that the median separation and median radial velocity are those of an LG close to apocentre, resulting in smaller variance.

The variance in separation and radial velocity become comparable to the observational uncertainty only for levels $L = 27$. While the precise threshold appears somewhat arbitrary, we conclude that all scales down to $L = 26$, or $\lambda_{cut} = 0.2$ cMpc, have had a significant impact on setting the LG’s orbital phase, and its current relative velocity and separation. For comparison, in Sawala et al. (2021), we had computed the average scatter in the position and velocity of central halos of mass $\sim 10^{12}$ M<sub>☉</sub> to be $\sim 20$ kpc, and $\sim 3$ km s<sup>-1</sup>, respectively, from perturbations of the same scale. That these values are somewhat lower than those reported here shows that the relative positions and velocity of the MW and M31 are more sensitive to small-scale changes. The scale of $\lambda_{cut} = 0.2$ cMpc is also sufficient to determine the existence of individual halos down to $2 \times 10^9$ M<sub>☉</sub>. In the hierarchy of scales that determine the formation of the LG, the scales that set the phase of the LG’s orbit are thus significantly smaller than those that determine the presence of LMC/M33 analogues.

As is evident from Figs 14 and 15, the parameters of LG analogues initially chosen from the 0.8-cMpc set are not necessarily at the centres of the different variations. Every LG variant contains phase information at all levels that we can resolve in our simulations, and represents just one random sample of the possible variations. When we consider perturbations to the orbital phase of an LG analogue, for individual observables, especially for LGs that are close to pericentre, the expected ‘reversion to the mean’ becomes a ‘reversion to apocentre’, where a larger fraction of the orbital phase is spent.

# 6 CONCLUSION

Beginning with the BORG reconstruction of the local LSS, we have constructed hybrid constrained initial conditions that embed the LG in its cosmic environment. Using the techniques described in Sawala et al. (2021), our simulations supplement the large-scale constraints with independent randomized phase information below $\lambda_{cut} = 6.48$ cMpc within a 16-cMpc sphere around the observer. Depending on the selection criteria, we find many approximate LG analogues, and for most of them, we find several, and often many variants of the same object. This leads us to conclude that we have, most likely, explored a large fraction of the event space of LGs that can form with these particular constraints assuming $\Lambda$CDM.

We have identified the scales of the initial density field that determine the different LG properties. We find that the classical
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Quality plots, analogous to Figs 6 and 12, of higher order variations at levels 25, 26, 27, and 28 ($\lambda_{cut} = 0.4, 0.2, 0.1$, and $0.05$ cMpc, respectively) of four existing variants, $[i = 21\,826, j = 1]$ (top row), $[i = 41\,800, j = 1]$ (second row), $[i = 43\,595, j = 1]$ (third row), and $[i = 24\,943, j = 366]$ (bottom row). In all cases, the scatter in observables decreases further with decreasing scale of variation. By level 27, the variances in separation and radial velocity become comparable to the observational uncertainties. At this scale of variation, the phase of the orbit is set.</figcaption></figure>

observables, such as the separation and radial velocities are not only transient in time, but also subject to very small scale changes in the initial density field: we find that the phase information needs to be fixed down to $\lambda_{cut} = 0.2$ cMpc to constrain the phase of the orbit. On the other hand, as discussed in Section 4.1.1, the orientation of the LG with respect to the LSS, the total mass, and integrals of motion such as the orbital energy and angular momentum, are set on larger scales ($\lambda_{cut} = 0.8$ cMpc). We suggest to consider these as the defining characteristics of the LG.

In constructing and refining cosmological initial conditions, this hierarchy of observables also has practical consequences. We find that in our simulations, the scatter in LG observables introduced by small-scale variations of the initial density field depends strongly on the presence of additional halos. Furthermore, given our constraints on the larger scales, the scatter introduced by small scale variations is required in order to produce an LG with the observed properties. This in itself does not necessarily imply that the presence of the LMC/M33 is necessary for an MW-M31 pair to match the other observables. However, given how strongly the presence of such objects affects the orbit, and given that their formation is set on scales greater than those that set the orbital phase, it is clear that the actual LG, in its true cosmic environment, would not have the observed kinematics without the presence of the LMC and M33. Consequently, simulations that produce an LG with a precise match to the other observables in the absence of those additional massive objects, would miss important aspects of its environment.

This work uses dark matter only simulations. Given the level of detail we require to identify LG objects, this is justified, but using hydrodynamic simulations would allow us to select also by galaxy properties. In principle, there is nothing other than computational cost that would prevent us from employing hydrodynamic simulations in our hierarchical method for constructing initial conditions.

In this work, we have been interested in exploring the sensitivity of observables to small scale variations of the initial density field. Of course, in constructing higher level variations, it would be possible to select specific candidates at each level for refinement at smaller scales, and to create ever more precise matches to the observed LG. We will fully explore smaller scale variations, and their impact on substructures, in future works.

All constrained simulations of the Local Universe should eventually converge to the same result. One major caveat in our current work is the residual variance and uncertainty in the reconstruction of the local LSS. As discussed in Section 2.2, and in more detail in a forthcoming paper, the constrained Local Universe in the SIBELIUS simulation reproduces many features of the observed universe. However, the resemblance is not perfect. For example, for the Milky Way analogue in the system that matches the observed LG kinematics, we measure a galactocentric velocity in the ‘CMB frame’ of $380 \text{ km s}^{-1}$ towards equatorial coordinates $(118^{\circ}, -8^{\circ})$, somewhat lower than, and slightly misaligned with the measured velocity of $552 \text{ km s}^{-1}$ towards $(157^{\circ}, -24^{\circ})$; Kogut et al. 1993).

However, this may only be a temporary limitation. In principle, our method of iteratively randomizing smaller scales allows us to explore all possible LGs that can form in the Local Universe, and to do so in a very efficient way. If the $\Lambda$CDM model was manifestly incorrect, we might not have found a precise LG analogue within the constraints used. With improving constraints, we now have a clear pathway towards not only creating LG analogues, but eventually reproducing the LG. Alternatively, our method would also allow us to demonstrate that the LG cannot possibly form, given a set of
large-scale constraints and cosmological parameters. We envisage
this to become a powerful and comprehensive way of testing
cosmological and astrophysical assumptions.

ACKNOWLEDGEMENTS

We thank the referee, Noam Libeskind, for his careful reading of
our paper and for several suggestions that led to improvements in
our paper. TS is an Academy of Finland Research fellow. TS and
SM are supported by the Academy of Finland grant 314238. PHJ
acknowledges the support by the European Research Council via
ERC Consolidator Grant (no. 818930). CSF acknowledges support
from European Research Council (ERC) Advanced Investigator
grant DMIDAS (GA 786910). This work was also supported by
the Consolidated Grant for Astronomy at Durham (ST/L00075X/1).
GL acknowledges financial support from the ANR BIG4, under
reference ANR-16-CE23-0002. JJ acknowledges support by the
Swedish Research Council (VR) under the project 2020-05143 –
‘Deciphering the Dynamics of Cosmic Structure’.

This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, and ST/R002371/1, Durham University, and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure. This work also used facilities hosted by the CSC – IT Centre for Science, Finland. We gratefully acknowledge the use of open source software, including PY-SPHVIEWER (Benitez-Llambay 2015), MATPLOTLIB (Hunter 2007), and NUMPY (Harris et al. 2020).

DATA AVAILABILITY STATEMENT

The data underlying this paper will be shared on reasonable request
to the corresponding author.

REFERENCES

Abazajian K., Survey f. T. S. D. S., 2009, ApJS, 182, 543

Aragon-Calvo M. A., 2016, MNRAS, 455, 438

Benitez-Llambay A., 2015, *py-sphviewer*: Py-SPHViewer v1.0.0, doi: 10.5281/zenodo.21703

Benson A. J., Lacey C. G., Baugh C. M., Cole S., Frenk C. S., 2002, MNRAS, 333, 156

Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2011, MNRAS, 415, L40

Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343

Bullock J. S., Kravtsov A. V., Weinberg D. H., 2001, ApJ, 548, 33

Carlesi E. et al., 2016, MNRAS, 458, 900

Dooley G. A., Peter A. H. G., Yang T., Willman B., Griffen B. F., Frebel A., 2017, MNRAS, 471, 4894

Eke V. R., Navarro J. F., Frenk C. S., 1998, ApJ, 503, 569

Enqvist K., Sawala T., Takahashi T., 2020, J. Cosmol. Astropart. Phys., 2020, 053

Erkal D. et al., 2019, MNRAS, 487, 2685

Faber S. M., Lin D. N. C., 1983, ApJ, 266, L17

Fattahi A. et al., 2015, MNRAS, 457, 844

Garrison-Kimmel S. et al., 2019, MNRAS, 487, 1380

Gottloeber S., Hoffman Y., Yepes G., 2010, preprint (arXiv:1005.2687)

Governato F. et al., 2010, Nature, 463, 203

Harris C. R. et al., 2020, Nature, 585, 357

Hoffman Y., Ribak E., 1991, ApJ, 380, L5

Hubble E. P., 1925, ApJ, 62, 409

Hubble E. P., 1927, PAAS, 5, 261

Hubble E. P., 1929, ApJ, 69, 103

Huchra J. et al., 2012, ApJS, 199, 26

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Ibata R. A. et al., 2013, Nature, 493, 62

Jasche J., Lavaux G., 2019, A&A, 625, A64

Jasche J., Wandelt B. D., 2013, MNRAS, 432, 894

Jenkins A., 2010, MNRAS, 403, 1859

Jenkins A., 2013, MNRAS, 434, 2094

Jenkins A., Booth S., 2013, preprint (arXiv:1306.5771)

Jones D. H. et al., 2009, MNRAS, 399, 683

Kam S. Z., Carignan C., Chemin L., Foster T., Elson E., Jarrett T. H., 2017, AJ, 154, 41

Karachentsev I. D., Karachentseva V. E., Huchtmeier W. K., Makarov D. I., 2004, AJ, 127, 2031

Kashibadze O. G., Karachentsev I. D., Karachentseva V. E., 2020, A&A, 635, A135

Kogut A. et al., 1993, ApJ, 419, 1

Lavaux G., Hudson M. J., 2011, MNRAS, 416, 2840

Libeskind N. I. et al., 2020, MNRAS, 498, 2968

Libeskind N. I., Yepes G., Knebe A., Gottlöber S., Hoffman Y., Knollmann S. R., 2010, MNRAS, 401, 1889

Li Y.-S., White S. D. M., 2008, MNRAS, 384, 1459

Lovell M. R., Eke V. R., Frenk C. S., Jenkins A., 2011, MNRAS, 413, 3013

Lovell M. R., Hellwing W., Ludlow A., Zavala J., Robertson A., Fattahi A., Frenk C. S., Hardwick J., 2020, MNRAS, 498, 702

McConnachie A. W., 2012, AJ, 144, 4

McConnachie A. W., Irwin M. J., Ferguson A. M. N., Ibata R. A., Lewis G. F., Tanvir N., 2005, MNRAS, 356, 979

Mateo M. L., 1998, ARA&A, 36, 435

Mei S. et al., 2007, ApJ, 655, 144

Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, ApJ, 524, L19

Nasonova O. G., de Freitas Pacheco J. A., Karachentsev I. D., 2011, A&A, 532, A104

Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563

Okamoto T., Gao L., Theuns T., 2008, MNRAS, 390, 920

Pawlowski M. S., 2018, Mod. Phys. Lett. A, 33, 1830004

Pawlowski M. S., Famaey B., Merritt D., Kroupa P., 2015, ApJ, 815, 19

Peñarrubia J., Gómez F. A., Besla G., Erkal D., Ma Y.-Z., 2015, MNRAS, 456, L54

Phelps S., Nusser A., Desjacques V., 2013, ApJ, 775, 102

Planck Collaboration et al., 2014, A&A, 571, A16

Read J. I., Gilmore G., 2005, MNRAS, 356, 107

Riess A. G., Fliri J., Valls-Gabaud D., 2012, ApJ, 745, 156

Robotham A. S. G. et al., 2012, MNRAS, 424, 1448

Rubin V. C., Ford W. Kent J., 1970, ApJ, 159, 379

Salmon J., 1996, ApJ, 460, 59

Salomon J. B., Ibata R. A., Famaey B., Martin N. F., Lewis G. F., 2016, MNRAS, 456, 4432

Sawala T. et al., 2016a, MNRAS, 456, 85

Sawala T. et al., 2016b, MNRAS, 456, 85

Sawala T., Frenk C. S., Crain R. A., Jenkins A., Schaye J., Theuns T., Zavala J., 2013, MNRAS, 431, 1366

Sawala T., Jenkins A., McAlpine S., Jasche J., Lavaux G., Johansson P. H., Frenk C. S., 2021, MNRAS, 501, 4759

Shao S., Cautun M., Deason A. J., Frenk C. S., Theuns T., 2018, MNRAS, 479, 284

Simionescu A., Werner N., Mantz A., Allen S. W., Urban O., 2017, MNRAS, 469, 1476

Skrutskie M. F. et al., 2006, AJ, 131, 1163

Sofue Y., 2015, PASJ, 67, 75

Somerville R. S., 2002, ApJ, 572, L23

Sorce J. G. et al., 2016, MNRAS, 455, 2078

Sorce J. G., Courtois H. M., Gottlöber S., Hoffman Y., Tully R. B., 2014, MNRAS, 437, 3586

Springel V., 2005, MNRAS, 364, 1105

Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
Stanek K. Z., Garnavich P. M., 1998, ApJ, 503, L131

Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068

Tollerud E. J., Boylan-Kolchin M., Barton E. J., Bullock J. S., Trinh C. Q., 2011, ApJ, 738, 102

van der Marel R. P., Alves D. R., Hardy E., Suntzeff N. B., 2002, AJ, 124, 2639

van der Marel R. P., Fardal M., Besla G., Beaton R. L., Sohn S. T., Anderson J., Brown T., Guhathakurta P., 2012, ApJ, 753, 8

van der Marel R. P., Fardal M. A., Sohn S. T., Patel E., Besla G., del Pino A., Sahlmann J., Watkins L. L., 2019, ApJ, 872, 24

Whitrow G. J., 1967, QJRAS, 8, 48

This paper has been typeset from a *TEX/LATEX* file prepared by the author.