Cosmicflows Constrained Local UniversE Simulations
Abstract
This paper combines observational datasets and cosmological simulations to generate realistic numerical replicas of the nearby Universe. These latter are excellent laboratories for studies of the nonlinear process of structure formation in our neighborhood. With measurements of radial peculiar velocities in the Local Universe (cosmicflows2) and a newly developed technique, we produce Constrained Local UniversE Simulations (CLUES). To assess the quality of these constrained simulations, we compare them with random simulations as well as with local observations. The cosmic variance, defined as the mean onesigma scatter of celltocell comparison between two fields, is significantly smaller for the constrained simulations than for the random simulations. Within the inner part of the box where most of the constraints are, the scatter is smaller by a factor 2 to 3 on a 5 Mpc scale with respect to that found for random simulations. This onesigma scatter obtained when comparing the simulated and the observationreconstructed velocity fields is only 104 4 km s i.e. the linear theory threshold. These two results demonstrate that these simulations are in agreement with each other and with the observations of our neighborhood. For the first time, simulations constrained with observational radial peculiar velocities resemble the Local Universe up to a distance of 150 Mpc on a scale of a few tens of megaparsecs. When focusing on the inner part of the box, the resemblance with our cosmic neighborhood extends to a few megaparsecs ( 5 Mpc). The simulations provide a proper Large Scale environment for studies of the formation of nearby objects.
keywords:
Techniques: radial velocities, Cosmology: largescale structure of universe, Methods: numerical1 Introduction
The formation of structures in the Universe, from tiny fluctuations at the era of recombination to the large diversity observed today, is a highly nonlinear process. Its multiscale nature is best studied by numerical experimentation. Cosmological simulations of structure formation rely on the cosmological principle which assumes the homogeneity of the Universe on large enough scales. The random nature of the primordial Gaussian perturbation field however implies that separate patches of the Universe are not identical. In that context, properties of patches of a few megaparsecs on aside vary widely. To overcome this cosmic variance, statistical comparisons to study the formation of structures are based on large observational datasets (e.g. Stoughton et al., 2002; Abazajian et al., 2003, 2009) and large cosmological simulations (e.g., Klypin et al., 2011; Alimi et al., 2012; Prada et al., 2012; Angulo et al., 2012; Watson et al., 2014; Klypin et al., 2014; Skillman et al., 2014). In order to also resolve small scale structures in the entire required large computational boxes, the mass resolution must be sufficiently high, resulting in simulations that can be time consuming and expensive.
An alternative approach is to reduce the cosmic variance by focusing on the numerical study of structure formation in the nearby Universe and a direct comparison of theoretical models with local observations. There is a double advantage of producing simulations resembling the Local Universe. First, our neighborhood is without any doubt the bestobserved volume of the Universe and as such hides its own cosmological treasures, often as fossils from its early epochs. Second, because a very large box size is not required for such simulations, the desired high resolution can be more easily achieved without being overly time consuming. However, this double advantage requires a correct modeling of the initial conditions based on the observed structures in the Local Universe. Standard cosmological simulations are obtained with initial conditions drawn from a random realization of the primordial perturbation field for a given cosmological model. Observational data of the Local Universe are used as additional constraints on these initial conditions so that resulting simulations resemble the Local Universe. The resulting constrained simulations attempt to describe the evolution of the observed structures in the nearby Universe.
The billions of data points that characterize the initial conditions of a cosmological simulation cannot be constrained by only thousands of observational data. The main aim of the different existing reconstruction techniques is to reduce the cosmic variance of the resulting constrained simulations. Local object candidates can then be identified to study their formation and evolution in the proper environment. Since the introduction of the POTENT method and the first attempt to reconstruct initial conditions from sparse observational data of the velocity field traced by galaxies (Dekel et al., 1990; Bertschinger et al., 1990; Nusser & Dekel, 1992) a lot of progress has been made over the last two and half decades. The constrained realization technique proposed by Hoffman & Ribak (1991) constitutes such a significant step forward in developing further the reconstruction method. With this technique, Ganon & Hoffman (1993) have constructed constrained initial conditions using the POTENT data that led to the first body simulation that mimicked the matter distribution around the Local Group in a 256 Mpc box (Kolatt et al., 1996).
The first step in generating initial conditions for constrained simulations consists in reconstructing today’s threedimensional density field from sparse and noisy observational data of the local galaxies. These observational data can be either the positions and radial peculiar velocities of galaxies (Kravtsov et al., 2002; Klypin et al., 2003; Sorce et al., 2014a) or redshift catalogs (Heß et al., 2013). In a second step the initial linear density field must be retrieved to derive the initial conditions for the simulations. There are two possible ways: backwards as described is section 2 of this paper or forwards as recently proposed by Kitaura (2013); Heß et al. (2013); Jasche & Wandelt (2013); Wang et al. (2013). In the latter case, the initial density field is sampled from a probability distribution function consisting of a Gaussian prior and a likelihood. A complete overview about the different methods is given in Wang et al. (2014).
Coming back to the observational dataset used here to constrained the initial conditions, this paper is part of the CLUES project^{1}^{1}1http://www.cluesproject.org/ (Constrained Local UniversE Simulations, Gottlöber et al., 2010; Yepes et al., 2014) which focuses on using peculiar velocities as constraints. Within this project, indeed, a number of constrained simulations of the Local Universe, using peculiar velocities from Karachentsev et al. (2004); Willick et al. (1997); Tonry et al. (2001) as observational constraints, have been run. Although estimating peculiar velocities constitutes an observational challenge, there is a double advantage in using them: first they are highly linear and correlated on large scales ; second, they are excellent tracers of the underlying gravitational field as they account for both the baryonic and the dark matter.
However, the first generation of these constrained simulations was affected by a substantial shift in the positions of objects recovered at redshift zero. To reduce this shift a new technique has been developed to account for the cosmic displacement field (Doumler et al., 2013a, b, c) and the noisy velocity information (Sorce et al., 2014a). The method has been applied to the first catalog of the Cosmicflows project^{2}^{2}2http://www.ipnl.in2p3.fr/projet/cosmicflows/ (Tully et al., 2008) producing simulations resembling the Local Universe down to a few megaparsecs within 30 Mpc (Sorce et al., 2014a).
At the end of 2013 (Tully et al., 2013), a second catalog of the Cosmicflows project was released. Superior in size (number of constraints, 8000 against 2000) and extent ( 150 Mpc against 30 Mpc) to the first catalog, it constitutes an ideal supplier of observational data to construct initial conditions constrained by the observed peculiar velocities of galaxies. Uniting the second release of the Cosmicflows data and the newly developed technique to build more accurate constrained initial conditions, this paper aims at demonstrating by how much the cosmic variance is decreased with respect to random simulations and how far the resemblance with the Local Universe can be extended.
The structure of the paper is as follow. In the second section, we briefly review the methodology described in Sorce et al. (2014a) and Sorce (2015) which combines several techniques to reduce biases in the observational catalog and to build constrained initial conditions. The third section presents the resulting constrained simulations of the Local Universe: the cosmic variance is estimated, the Large Scale Structure is compared with galaxies from the 2MASS redshift catalog (Huchra et al., 2012) and with the threedimensional reconstruction of the overdensity and velocity fields of the Local Universe at redshift zero. In addition, the simulated Laniakea Supercluster of galaxies, a basin of attraction of local velocity flows, is compared with the observed Laniakea Supercluster (Tully et al., 2014). Finally, before concluding, dark matter halo candidates for well known nearby clusters, such as Virgo, are identified in the simulations.
2 Methodology
In this section we discuss how the initial conditions have been constructed from the observational data: a set of positions and radial velocities of galaxies. Important steps in this procedure are the grouping of galaxies and the minimization of biases.
2.1 Observational data: Cosmicflows2
Cosmicflows2 is the second generation observational catalog of galaxy distances built by the Cosmicflows collaboration. Published by Tully et al. (2013), it contains more than 8,000 accurate galaxy distances. Distance measurements come mostly from the TullyFisher relation (Tully & Fisher, 1977) and the Fundamental Plane methods (Colless et al., 2001). Cepheids (Freedman et al., 2001), Tip of the Red Giant Branch (Lee et al., 1993), Surface Brightness Fluctuation (Tonry et al., 2001), supernovae of type Ia (Jha et al., 2007) and other miscellaneous methods also contribute to this large dataset but to a minor extent ( 12%) although they have individually higher weights because of smaller errors.
The final goal of the paper is to build initial conditions constrained by this catalog in order to run cosmological simulations, working above the scale of galaxy virial motions in clusters (nonlinear displacements) is required. Therefore, in this paper, the grouped version of cosmicflows2 is used. With a method similar to that described in Tully (2015a, b), 552 groups and 4303 single galaxies can be identified in the dataset shrinking the number of constraints to 4855. The resulting number density is large enough to construct constrained initial conditions as demonstrated by Doumler et al. (2013a, b, c).
Figure 1 shows the normalized cumulative distribution of distances in cosmicflows2 as well as the mean and median distances. The catalog extends up to 230 Mpc but on average constraints are within 66 Mpc ( 61 Mpc for the median). In fact, about 60 % of the constraints are within a sphere of 70 Mpc radius, approximately 80 % are within 100 Mpc, and 98 % are within 160 Mpc. From this distribution, the inner box is expected to be the most constrained part of the simulation and beyond 160 Mpc the random component is anticipated to overcome the constraints. Intermediate results should be found in between.
2.2 Minimizing Biases in the Observational Catalog
Tully et al. (2013) warned us that cosmicflows2 is affected by biases with effects that cannot be ignored anymore as the effects are stronger with increasing distance. There are four biases known:

In the literature the first bias has been given a number of terms that are used interchangeably. These are Problem I, Selection Effect/Bias, r against V, Distancedependent, Frequentist, Calibration problem, Mbias of the second kind (Kaptney, 1914; Malmquist, 1922; Han, 1992; Sandage, 1994; Teerikorpi, 1997, 1993, 1990; Hendry & Simmons, 1994; Willick, 1994; Teerikorpi, 1995). This bias is analogous to a selection effect in magnitude (dim galaxies are selectively excluded from the observational sample) resulting in underestimated distances. This bias has been nulled in the observational cosmicflows2 catalog (e.g. Tully & Courtois, 2012; Sorce et al., 2013, 2014b).

The second bias, referred to as Homogeneous Malmquist Bias or Problem II, General Malmquist Bias, Geometry Bias, V against r, Classical, Bayesian, Inferreddistance problem, Mbias of the first kind in the literature (Kaptney, 1914; Malmquist, 1920; LyndenBell et al., 1988; Han, 1992; Teerikorpi, 1997; Sandage, 1994; Teerikorpi, 1993, 1990; Hendry & Simmons, 1994; Teerikorpi, 1995; Strauss & Willick, 1995), is due to the fact that the number of observable galaxies from our perspective increases with the distance. There are more galaxies available to scatter inward due to errors than outward, creating the tendency to locate galaxies closer than they should be, namely to underestimate distances.

On top of this last bias, there is a lognormal error distribution or an asymmetry in the distribution of fractional errors on distances because distances are derived from distance moduli via a logarithmic function. Thus if a same galaxy is located farther than it should be rather than closer, the error on the distance is larger although this is not reflected by the assigned fractional errors.
These biases lead mainly to a major infall onto the Local Volume. An iterative method to minimize the infall and reduce spurious nongaussianities in the radial peculiar velocity distribution was applied to obtain a new distribution of radial peculiar velocities and corresponding distances (Sorce, 2015).
2.3 Building Initial Conditions
To build initial conditions for dark matter only numerical simulations (a set of particles with velocities and positions at a starting redshift) constrained by peculiar velocities, we rely on four techniques assuming a prior cosmological model:

The WienerFilter method (WF, Zaroubi et al., 1995) to reconstruct the cosmic displacement field required to account for the displacement of constraints from their precursors’ locations,

the Constrained Realization (CR, Hoffman & Ribak, 1991) of Gaussian field technique to produce overdensity fields constrained by observational data, adding a random realization to compensate for the missing power spectrum. These latter are then converted into white noise that can be used to increase the resolution,

the resolution is increased by adding some random small scale features in the white noise, then the white noise is converted back to build initial conditions for cosmological simulations.
These four steps can be summarized in a set of equations:
(1) 
where = x, y, z and C are the constraints plus their uncertainties and is the growth rate (D the growth factor and a the scale factor). Brackets denote correlation functions depending solely on the assumed prior cosmological model (here the power spectrum) ; is the velocity field and is the displacement field. The ’WF’ exponent means that a field is obtained with the WienerFilter method,
(2) 
where is the approximate location of the constraints’ precusors (linear theory at first order valid down to 2 Mpc), is the measured position of the constraints. Note that the left of Eq. 2 is an extension of equation 14 in Nusser et al. (1991).
(3) 
where stands for overdensity fields. stands for the random realization field and each represents a random constraint drawn from . represents the constrained realization field.
(4) 
where stands for the white noise in Fourier space, represents the power spectrum, n is the number of particles and V the volume of the simulated box. This last equation is first used in reverse to convert the overdensity fields in white noises, then random small scale features are added to increase the resolution and the equation is used again to obtain the higher resolution density fields to prepare the initial conditions.
We apply the whole scheme to the observational catalog cosmicflows2 within the framework of Planck cosmology (=0.307, =0.693, H=67.77, , Planck Collaboration et al., 2014).
3 Constrained Local UniversE Simulations
In the second section, we noted that the largest distance of the cosmicflows2 catalog is 230 Mpc, thus the size of the computational box should be sufficiently large to avoid spurious effects due to periodic boundary conditions where the observational data still have a constraining power. Tests have shown that a 500 Mpc box meets this requirement. Such a computational volume extends the study up to the Shapley supercluster (the location of the farthest constraint) and with particles, the mass resolution is sufficient to resolve large groups and clusters of galaxies. Two of the simulations were also run at higher resolutions (1024 particles) to check that none of the conclusion drawn in this paper are affected by the 512 choice. As none of them are, we settle with the reasonable choice of 512 particles. A total of 25 constrained simulations and 15 random simulations have been performed in order to study the residual cosmic variance.
The initial conditions of these different simulations have been constructed in various ways:

The first fifteen constrained initial conditions are built out of different random realization fields plus the observational dataset (see equations 1 to 3). This first step uses the Wiener Filter, the reverse Zel’dovich approximation and the constrained realization technique (Zaroubi et al., 1999; Doumler et al., 2013c; Sorce et al., 2014a; Hoffman & Ribak, 1992) to generate density grids in accordance with the expected minimum scale on which the constraints are effective. The corresponding white noise field is used as an input to increase the resolution to particles with the Ginnungagap code^{3}^{3}3https://github.com/ginnungagapgroup/ginnungagap . This full ‘MPI+OpenMP’ parallel code adds random small scale fluctuations in real space to increase the resolution to any level within a given cosmology. The resolution limit is dictated only by the total memory of the supercomputer. A final simulation is then characterized by two random seeds: the random realization field and the added small scale features ;

The additional ten constrained initial conditions share the same random realization field , but different seeds are used to increase the resolution to 512. They are thus expected to differ only on scales smaller than that of the input white noise field ;

Finally a set of fifteen random, i.e. not constrained, initial conditions has been constructed. They share the same seeds (same ) as the fifteen constrained initial conditions.
The simulations based on all these initial conditions have been performed with Gadget3 (Springel, 2005) from redshift 60 to redshift 0 with a 25 kpc force resolution.
In Figure 2, we first compare the resulting power spectra at to the linear power spectrum of the chosen cosmology. As expected the 15 random simulations (grey area) scatter at large scales around the linear input power spectrum (blue solid line). The 15 constrained simulations (dotted black line) tend to have less power on large scales, an effect which decreases with the box size as shown in the Appendix, because a smaller and smaller fraction of the box is constrained. The bottom panel of Figure 2 represents the power spectra divided by the Planck power spectrum as well as the mean values. Although on the low side, the power spectra of constrained simulations are within the scatter obtained with those of random simulations. Their mean is on average smaller by a factor 1.3 (factor that decreases with the scale) than that of random simulations on the large scales. As for the ten constrained simulations built out of the same random realization field, they share the same Large Scale power spectrum (red dashed lines) in agreement with the fact that the random features added to increase the resolution affect only the small scales. In the middle panel, the ratio of the power spectra to their mean for the three samples (15 constrained simulations, 10 constrained simulations with the same seed in and 15 random simulations) is displayed. By construction this is indistiguishable from a straight line for the 10 simulations where only small scale structures have been added.
Next, the Amiga’s Halo Finder (Knollmann & Knebe, 2009) is applied to each simulation to compile a list of dark matter halos. In Figure 3, using M defined with respect to the critical density, the cumulative mass functions of the different simulations are compared with the same color code as that of Figure 2. The blue color now stands for the Tinker cumulative mass function as defined by Tinker et al. (2008) using the online mass calculator of Murray et al. (2013) and the Planck cosmological parameters as defined in section 2 of this paper. In the left panel of the figure, the cumulative mass functions for the entire 500 Mpc boxes are shown (top), as well as their cumulative mass functions divided by the Tinker cumulative mass function (middle) and their scatters around their respective mean (bottom). The cumulative mass functions of the different simulations overlap and, as expected a smaller scatter is observed for the constrained simulations sharing the same random realization field (in red). In order to evaluate the effect of constraints on the cumulative mass functions, these latter are derived in a sphere of radius 160 Mpc centered on the original box, namely on the center of the box. A radius of 160 Mpc is a reasonable choice as it encompasses 98% of the constraints (see section 2). The corresponding cumulative mass functions are shown in the right panel of Figure 3. One can clearly see that at high mass ranges the cumulative mass functions of the constrained simulations are on the low side compared to the cumulative mass functions of the random realizations computed in the same sphere. These latter, contrary to the former, tend to scatter symmetrically around the Tinker cumulative mass function. As expected the scatter of the cumulative mass functions in the sphere is smaller for constrained simulations compared to random simulations (bottom of the right panel of Figure 3) and is the smallest for constrained simulations sharing the same random realization field. In this respect the cosmic variance is reduced in constrained simulations with respect to random simulations. We investigate more thoroughly the residual of the cosmic variance in the following subsection.
3.1 Reduced Cosmic Variance
In this section we compare and quantify the cosmic variance within the different sets of simulations, i.e. constrained and random. To this end, a cloudincell scheme on a grid is applied to the particle distributions of the initial conditions (z=60) and of the simulations at z=0 with a subsequent Gaussian smoothing on a scale of 5 Mpc. Normalized by the mean density, the resulting smoothed density fields of any pair of constrained and random simulations are compared cell by cell. For each pair we build a densitydensity plot (the density field of a first simulation versus the density field of a second simulation). If the two simulations were identical all points would follow the 1:1 relation. We define the cosmic variance between two simulations as the onesigma, hereafter 1, scatter (or standard deviation) around this 1:1 relation. We repeat this procedure for the 105 pairs of the 15 random simulations and those of the 15 constrained simulations as well as the 45 pairs of the 10 constrained simulations. Then we calculate the mean and the variance of the 1 scatters. The result is three points (filled dark grey, black and light grey circles for each one of the simulation types: random, constrained, constrained sharing the same random realization) in each of Figures 4 and 5, with error bars at the xaxis value of 500 Mpc at z=60 and z=0 respectively. We repeat this procedure in smaller subboxes of size 400, 300, 250 Mpc, etc, centered on the original box in order to measure the cosmic variance in the smaller central volumes where most of the observational data are, i.e. where they are the most effective.
In Figure 4, at the starting redshift, the mean scatter of the random initial conditions is independent of the size of the subbox with increasing variance in smaller subboxes. On the contrary the mean scatter of the 15 constrained initial conditions starts to decrease substantially below 300 Mpc. As the median distance in cosmicflows2 is only 61 Mpc and 98% of the measurements are within 160 Mpc, the standard deviation between density fields in the initial conditions is reduced in the expected large volume. The 10 constrained initial conditions with the same random realization field and added small scale features are on a line because adding modes at the scale of 0.98 (500/512) Mpc does not affect the smoothed density distribution shown here. It only slightly increases the variance in the smallest subbox.
The top of Figure 5 shows the mean 1 scatters at redshift as a function of the subbox size for the different types of simulations (random and constrained). The effect of the nonlinear clustering is visible on small scales. As the system evolves and becomes more nonlinear, a larger fraction of the subbox is covered by low density regions (voids) than high density regions (clusters). Thus, voids become dominant in the volume. The densities of these nearly empty regions tend asymptotically to zero regardless of the initial field while that of high density regions, rising from small but positive differences in the initial field, are magnified. Consequently, the probability to compare a cell in a low density region with one in another low density region, with similar values, increases, reducing the scatter for each pair of random simulations. When considering subbox smaller than 100 Mpc, the 1 scatters decrease on average for the random simulations although the variance of these scatters increases because there is still the probability to find a high density region in one of the pair simulation versus a low density region in the other simulation of the pair. This is a limitation of the comparison method because the smaller the subbox considered the higher the probability to find the same kind of field even if the probability to hit different types of fields is non null. The effect is not that pronounced for the constrained simulations because of the existence of similar (by construction) structures close to the center of the box (such as the Great Attractor and the Virgo cluster) and disappears for the constrained simulations with the same random realization field. We repeat the same procedure for the three components of the velocity field in the bottom panel of Figure 5 and make the same observation. In addition, as a reference to assess the low values of the 1 scatters and thus the small discrepancies between the constrained simulated velocity fields, one can consider the validity of the WienerFilter reconstructed velocity field which is [100150] km s (Sorce, 2015).
To summarize, Figures 4 and 5 reveal that the cosmic variance, measured with the 1 scatter of celltocell comparisons, is considerably reduced for constrained simulations, by a factor 2 to 3 on a scale of 5 Mpc for both density and velocity fields in the inner part of the box, when compared to those obtained for random simulations. In addition, the error bars of these 1 scatters are smaller by a factor at least 2 when considering constrained simulations with respect to random simulations. As an example, taking a subbox of 150 Mpc, the cosmic variance is decreased from 0.5 to 0.3 for the density fields normalized by the mean density at z=0 and from 150 to 50 km s for the velocity fields.
In Figure 5, we compare density and velocity fields, in (sub)boxes, of constrained and random simulations. The observations, however, approximate a spherical distribution with a large radial dependence as demonstrated in Figure 1. Thus, corners of the (sub)boxes appear at first to be less constrained than other regions at the edges of the (sub)boxes. One might argue that spherical subvolumes would be better suited for comparisons between random and constrained simulations. Using spherical regions to derive 1 scatters, we find no difference.
Next, we apply different smoothing to the density fields normalized by the mean density to determine the constraining power, of the observational dataset combined with the method to build initial conditions, not only as a function of the box size but also as a function of the resolution. We define the resolution as where is the value of the Gaussian smoothing in Mpc, and the constraining power as the ratio of the mean 1 scatter obtained for the random simulations to that derived for the constrained simulations. We give to values ranging from 1 to 8 Mpc. Figure LABEL:fig:constrpower display the constraining power of the observational peculiar velocities combined with the method to build initial conditions in a resolution versus subbox size plane: the darker the grey, the higher the ratio of the 1 scatters, the higher the constraining power is (the more the cosmic variance is decreased in the constrained simulations with respect to the random simulations). As expected, the larger the smoothing (the smaller the resolution), the higher the constraining power is. On the other hand, with smaller smoothing (higher resolution), the cosmic variance is larger. In addition, the larger the subbox size, the smaller the constraining power is because of the decreasing number of constraints. One can also notice that the effect of nonlinearities becomes prominent as the smoothing decreases: small scales are essentially not constrained by the observational data. The non linear theory threshold is reached, on smaller scales shell crossing has wiped the initial correlations.
In the next section we compare the constrained simulations to the observed Local Universe.
3.2 Simulation of the Large Scale Structure
In the previous subsection a study of the cosmic variance has demonstrated that the constrained simulations are remarkably similar to each other in the constrained part of the simulation box. To compare these simulations with the observed Local Universe, an observer is assumed to be at the center of the box and the three supergalactic coordinates are defined similarly to observational supergalactic coordinates. Comparing observations with simulations is not an easy task. Two possibilities are available although they both involved their own limitations: 1) comparing the Large Scale Structure in simulations to the distribution of observed galaxy surveys which are however magnitude limited, 2) comparing the fields of the simulations with the reconstructed ones obtained with the WienerFilter technique from the observational data. They constitute however only the linear fields and tend to the null fields in absence of data or in presence of noisy data. These limitations highlight the importance of the simulations which give access not only to the formation history but also to the full fields including nonlinearities of the Local Universe.
Before comparing reconstructions, simulations and redshift surveys, we begin with a description of the observed structures in both the reconstruction and in one of the chosen randomly simulation as shown in Figure 7. Note that the choice of the simulation has no impact on the following discussion as simulations present the same Large Scale structure in the constrained 200 Mpc radius area, namely high density regions and voids in this area are in every simulation. In this figure both the reconstructed and simulated (over)density (contour) and velocity (arrows) fields are displayed in a 5 Mpc thick slice of the XY supergalactic plane. On top of the fields, red dots represent galaxies from the 2MASS redshift catalog in a 10 Mpc thick slice. These galaxies are superimposed for comparisons purposes and one can notice fingersofgod in the galaxy distribution. Several wellknown structures can be identified in both the reconstructed and the simulated Local Universe: PerseusPisces (PP), Shapley as well as Coma superclusters but also voids such as the Sculptor void. In addition to the major structures and voids, the Zone of Avoidance (ZOA) due to our MilkyWay dust is marked. Note, that no structure have been reconstructed in that zone beyond 50 Mpc from the center of the box due to a lack of information in the observed data. However, the simulation shows structures in this region, in particular connections between objects above and below the ZOA.
Little is known about structures in the Zone of Avoidance but the simulation reproduces quite well the observations also in that zone:

A potential supercluster, at a distance about 60 Mpc in the SGX direction, situated in the zone of obscuration (KraanKorteweg et al., 1994) is on an extension of the filament departing from Hydra and Antlia clusters going across the Zone Of Avoidance to reach the region of the Great Attractor (Centaurus Supercluster, KraanKorteweg et al., 1994).

KraanKorteweg et al. (1994) noted a clustering at a distance greater than 100 Mpc in the SGX direction, in the zone hidden by our galaxy dust, a potential connection between the Horologium and Shapley Superclusters. The simulation contains a high density zone beyond 100 Mpc in that direction.
Regarding comparison with redshift surveys, there is qualitatively a good agreement between the 2MASS redshift catalog and the simulation on Figure 7. To assess this agreement we use the cosmic web based on the velocity shear tensor (Hoffman et al., 2012) rather than on the gravitational tidal tensor (e.g Hahn et al., 2007) or on the displacement tensor (Lavaux & Wandelt, 2010). Eigenvalues of the velocity shear tensor permits to determine whether a region of the Universe with such a velocity field constitutes a knots, a filament, a sheet or a void. It is thus straightforward to determine the position of a galaxy in the cosmic web. With a null threshold and the definition used in Hoffman et al. (2012) for the velocity tensor, three negative eigenvalues correspond to a void while three positive value correspond to a knot. Two negative and one positive value constitute a sheet while the opposite stands for a filament. If simulations are in good agreement with observations, there should be approximately the same number of galaxies in filament and sheets ( 3545 %) then less in knots and in voids ( 10 %) (e.g. ForeroRomero & González, 2015; Libeskind et al., 2012, although they choose a slightly higher than zero threshold, while we choose zero, they checked that general results are quasi independent of the threshold choice as long as the chosen value is reasonable). On average of the fifteen cosmic webs computed from the fifteen different constrained simulations, 61 % of the galaxies are in knots, 352 % are in filaments, 482 % are in sheets and 101 % are in voids. The galaxies are distributed as expected giving agreement with observations.
As for comparisons with the reconstruction, coming back to Figure 7, the WienerFilter reconstructs fairly well the Local Universe in the center of the box although it shows only the linear fields and tends to the null field in absence of data or in presence of noisy data, thus reconstruction and simulation agree qualitatively very well: the reconstruction presents more feature in the center relatively to its edges but the loss of precision with the distance from the center of the box is the cause. The simulation allows one to go deeper into the Zone of Avoidance and to extend further the study of the Large Scale Structure and, more importantly, it supplies the whole density field, including nonlinearities. An estimation of the agreement between reconstruction and simulations can be made with celltocell comparisons between the velocity fields which, unlike densities, are highly linear. The 1 scatter is on average of the order of 100150 km s (i.e. 23 Mpc, the linear theory threshold, in terms of displacement). From this comparison between a simulation and the reconstruction, it can be concluded that the major attractors and voids of the Local Universe are properly simulated.
The variance between the different constrained simulations is relatively low (see previous subsection). If a Large Scale Structure feature is present in one of the simulations, there is a high probability to recover it in all the other realizations. The Large Scale Environment is robustly simulated. Since this Large Scale environment has been suggested to play an essential role in the formation and evolution of local objects (e.g. GarrisonKimmel et al., 2014), these constrained simulations are ideal to study local objects. First, we turn towards a recently discovered Large Scale feature in our neighborhood, the Laniakea supercluster of galaxies (Tully et al., 2014).
3.3 Example of the Laniakea Supercluster
In this subsection, we focus on a particular structure of the Local Universe, the Laniakea Supercluster of galaxies discovered and defined in Tully et al. (2014). This supercluster is constituted of a local basin of attraction and every objects with a perturbative motion toward it. Local flows encompassed in that region converge onto the local attractor. To compare the simulated superclusters with the observedreconstructed one, the divergent field is evaluated. For a chosen volume, the divergent field corresponds to velocities due solely to densities in this volume. In order to remove most of the non linear components, the density and velocity grids are smoothed with a 5 Mpc Gaussian. With grids obtained via a cloudincell scheme and a 60 Mpc sphere radius centered at [47,13,5] Mpc similarly to the definition given by Tully et al. (2014), we are able to find the local basin of attraction in the different simulations. Enlarging the radius, the contours of the simulated Laniakea superclusters of galaxies are recovered. The result obtained with one realization is given on the right panel of Figure 8. The left panel shows the WienerFilter reconstructed supercluster. In this figure, the gradient of colors correspond to densities while streamlines stand for the velocity fields. From blue to red the density increases, namely voids are in blueblack and high density regions are in yellowred. The simulation and the reconstruction look alike. The Laniakea supercluster streamlines are very well simulated and converge in a similar location about [,11,0] Mpc. The Laniakea supercluster is surrounded by cosmic gravitational streams flowing towards the PerseusPisces superclusters on the positive SGX side, towards Shapley on the negative SGX side and towards Coma on the positive SGY side in both the simulation and the reconstruction. Tully et al. (2014) estimate the mass of the supercluster at around 6.5 10 on the scale of this paper, in relatively good agreement with simulations where the number of dark matter particles contained in the Laniakea supercluster, simplified as a 60 Mpc sphere radius centered at [47,13,5] Mpc, gives a mass of approximately 2 0.3 10 over the different realizations.
To evaluate the agreement between the WienerFilter reconstruction and the simulations, we proceed as in the previous section with celltocell comparisons within the Laniakea supercluster region. The 1 scatter is 104 km s on average with a standard deviation of 4 km s, the median is identical to the mean. Celltocell comparisons in this region between constrained simulations give on average 1 scatters of 45 6 km s for the velocity fields and of 0.29 0.02 for the density fields at z=0. By comparison, the average 1 scatters obtained when comparing random simulations in that region are about 145 35 km s for the velocity fields and of 0.43 0.09 for the density fields.
3.4 Observed Clusters & Simulated Dark Matter Halos
Finally, we turn our attention to the study of halos at redshift zero. Lists of dark matter halos obtained with the Amiga Halo Finder (Knollmann & Knebe, 2009) are used to match halos in simulations with clusters in the Local Universe. Virgo is one of the largest cluster in our neighborhood and it is the closest. It is natural to look for candidates of that cluster as it should be the most efficiently constrained. In addition, attempts to find candidates for Centaurus (restricted to one component of a complex in the region), Coma and Perseus are made.
We begin with the Virgo cluster. Candidates are found at less than 34 Mpc from its observational position in all the simulations. The replicas have masses (M with respect to the critical density) for the simulations between 2.7 and 4.3 10 h which is in good agreement with current values given for this cluster (e.g. 2  6 10 in Planck cosmology, Karachentsev & Nasonova, 2010), especially considering that within simulations, masses of cluster candidates can vary between half and twice the cluster mass (Ludlow & Porciani, 2011) and that the estimated mass of a cluster depends on the method/definition used, as does the mass of a simulated halo. Several parameters can be compared between cluster candidates and the observed cluster. To estimate the agreement between candidates and observed Virgo cluster, we define the inaccuracy as the difference between the simulated and the observed (mass, distance, coordinate) or WienerFilter reconstructed (velocity) component divided by: 1) the mean distance in cosmicflows2 for the coordinates and distances, 2) the mean velocity of halos of approximately the same mass as Virgo candidates for the velocities and their components, 3) the mean estimated mass of the observed clusters for the masses. Note that the definition we choose for the inaccuracy differs from that of the relative change or difference as we normalize neither by the reference (i.e. observed or reconstructed) parameter nor by a function of the simulated and observed or reconstructed parameters. We justify this choice by the fact that furtherfromtheboxcenter halos would artificially appear on Figure 9 in better agreement with the observed clusters than their counterparts using the latter definitions (i.e. division by the distance). The same artificial observation would happen if more massive (faster) halos were normalized by the mass (velocity) of the halos. For this reason, the normalization is made with reference values independent of the cluster characteristics. The compared parameters are the three supergalactic coordinates and the square root of the sum of their squares, namely the distance, the three components of the velocity and the 3D velocity and the mass (M with respect to the critical density for the simulations).
Figure 9 gathers the mean and scatter of the inaccuracies of the parameters as defined above for the Virgo candidates found in the different constrained simulations. The top panel shows the mean inaccuracies (filled black circles) and their standard deviations (error bars) for the Virgo candidates in the 15 simulations with a different random realization field while the bottom panel gives inaccuracies and standard deviations of Virgo candidates found in the ten simulations sharing the same random realization field. There are two observations: 1) there is a very good agreement between the observedreconstructed and simulated parameters as inaccuracies are close to zero. 2) the standard deviation between the different simulations are quite small (less than 1015% of the considered parameter). As expected, the ten simulations which differ only on the random small scale features give Virgo candidates in slightly better agreement with each other (scatter less than 510% of the considered parameter).
We proceed similarly for Centaurus, Coma and Perseus although we expect candidates to be somewhat more scattered because the simulations, although constrained on the Large Scale, are not as well constrained as they are in the very center of the box, near Virgo candidates. Rather than using the estimated masses, we settle for looking for halos more massive than 5 10 close to the positions of the observed clusters. If present these halos confirm a high density regions where expected. In only one of the fifteen constrained simulations, we cannot locate a halo massive enough in a reasonable radius around the observed position of Coma and Centaurus clusters. As for the Centaurus location, in two simulations we have two candidates of approximately the same mass nearby (the distance between the two is less than 3 Mpc), we select the closest to the estimated observed position. Note that finding two halos close to each other is not surprising as Centaurus is constituted of several massive components. We then proceed as for Virgo candidates, i.e. we derive the inaccuracies for Centaurus, Coma and Perseus candidates although we fix the mass estimate to 5 10 as the reference for these three other clusters. We justify our choice by the fact that 1) simulations are not as well constrained at the positions of these other clusters as they are at Virgo’s position, 2) the estimated mass of these clusters varies with the measurement method and the defined boundaries of the cluster (as an example, the estimated mass of the Coma cluster is 5.1 10 according to weak lensing measurements in Gavazzi et al. (2009) while it is 1.4 10 when using the radius of second turn around as define by Tully (2010)) and 3) the definition of the observed mass is different from the mass of simulated halos. Results are given in Figure 9 with three different colors for the three different clusters. They are slightly shifted to the right on the xaxis with respect to datapoints obtained for Virgo candidates for a clearer visibility. As for Virgo candidates, inaccuracies are small in terms of positions in agreement with the fact that the Large Scale Structure is well constrained. Scatters are larger than for the Virgo candidates but that was expected as we are looking at less constrained regions. This is not surprising that the higher scatter in velocities is observed for the Centauruslike halos as Centaurus is an ensemble of objects rather than a compact object: nonlinear motions are involved in that zone and because these regions are very dense, there are different massive halos with a high probability to pick a different one in each one of the realization. Looking at the bottom of Figure 9 which shows candidates found in the simulations sharing the same random realization field, the scatters are decreased by at least a factor 2.
Only objects with low mass are heavily modified by the random small scale features added to increase the resolution while the Large Scale Structure is quite unaffected by them. Considering that a goal of the CLUES project is to build a large statistical sample of Local Grouplike entities, selecting simulations containing all the appropriate clusterslike objects, we will be able to build a factory of lookalikes of the Local Group in the proper environment to study their statistical properties due to both the Large Scale Structure (constrained with different random realizations) and to the small scale structures (constrained but sharing the same random seed) in a decoupled way.
4 Conclusion
The first generation of simulations constrained by observational peculiar velocities produced by the CLUES project was affected by a substantial shift in the positions of objects recovered at redshift zero. In this paper, we present a double improvement with respect to this first generation. First, we use the second catalog of the observational Cosmicflows project which is superior in size (number of constraints), extent ( up to 150 Mpc) and accuracy compared to the previously used catalogs. Second, we add the newly developed techniques involving the grouping of galaxies, the minimization of biases, and the reverse Zel’dovich approximation based on the WienerFilter method, with the constrained realization technique to build more accurate constrained initial conditions. We are able to show that not only constrained simulations exhibit a lower cosmic variance than random simulations but also that they are in agreement with our cosmic neighborhood up to the nonlinear scale (23 Mpc).
To do so, we compare a set of 15 random and 25 constrained simulations (10 of the 25 simulations share the same Large Scale random phase and differ only by the small waves added to increase the resolution) of 512 particles within a 500 Mpc boxsize. A check with two 1024particles simulations showed that the results are not affected by the number of particles. We apply a cloudincell scheme to all the simulations and smooth the resulting velocity and density grids with a 5 Mpc Gaussian.
We define the cosmic variance as the onesigma, scatter (or standard deviation) in densitydensity plots (the field of a first simulation versus the field of a second simulation) obtained from celltocell comparisons between pairs of simulations of the same nature (random or constrained). We average the results over the different pairs and found that the 1 scatters obtained for constrained simulations are not only minimal when comparing the inner part of the boxes, where most of the constraints are, but they are also smaller by a factor 2 to 3 with respect to those found for random simulations. The best constrained part of the simulations is the inner box within approximately 100 Mpc for the smallest (clusters) scales, the resemblance extends to 300 Mpc on larger scales (5 to a few tens of megaparsecs). This agreement meets expectations as the cosmicflows2 catalog extends to 230 Mpc with 98% of the distance measurements within 160 Mpc and a median distance of 61 Mpc.
We found that on average the constrained simulations tend to have less power on large scales than the random simulations, although they are well within the expected scatter. This effect could be due to the observational data and/or to the way they are processed to build the initial conditions. We are working on an improved algorithm which reduces or removes this effect by taking into consideration the reduced accuracy of the reconstruction on large distances. This will be important when studying large scale velocity flows. However, results presented in this paper are in no way affected by this effect. In fact, adding artificially power on the largest modes would only slightly change the mass of the most massive objects. We will study and discuss this effect in more detail with a series of new simulations and discuss it in a forthcoming paper.
To compare the simulations with the observed Local Universe, we use celltocell comparisons between the reconstructed and the simulated velocity field. We find that simulations at redshift zero agree with the WienerFilter reconstruction obtained with the observations at 100150 km s or 23 Mpc, namely the linear theory threshold.
Taking as an example the Laniakea Supercluster of galaxies, defined as a local basin of attraction and all flows going towards it, we show that simulated and reconstructed Laniakea superclusters are in relatively good agreement. The mean 1 scatter obtained from celltocell comparisons between the reconstructed and simulated velocity fields is 104 4 km s. When comparing the simulations, the mean 1 scatter of the simulated Laniakea superclusters’ fields is 45 6 km s for their velocity fields and 0.29 0.02 for their density fields. By comparison, similar regions between random simulations differ by 145 35 km s for the velocity fields and 0.43 0.09 for the density fields.
Finally, we give an overview of Virgo candidates as well as other wellknown nearby clusters at redshift zero and show that again the scatter between simulated dark matter halo candidates within themselves and also with the observedreconstructed clusters in terms of position, velocity and mass is only of the order of 10%.
These comparisons show that simulations are in agreement between each other and above all with the reconstruction. Because the reconstruction recovers fairly well all the major attractors and voids of the Local Universe, they must also be present in all the simulations at redshift zero for these latter to be similar to the reconstruction, i.e. to the observations. The method to build more accurate constrained initial conditions is extremely efficient. We produced the first simulations constrained with observational radial peculiar velocities which resemble the Local Universe up to 150 Mpc, with increasing accuracy when reaching the inner part of the box.
Acknowledgements
The simulations have been performed at the Leibniz Rechenzentrum (LRZ) in Munich and at Jülich Supercomputing Centre. We thank the anonymous referee whose comments helped improving the manuscript. JS acknowledges support from the Alexander von Humboldt Foundation. SG and YH acknowledge support from DFG under the grant GO563/211. GY acknowledges support from the Spanish MINECO under research grants AYA201231101 and FPA201234694. YH has been partially supported by the Israel Science Foundation (1013/12). HC acknowledges support from the Lyon Institute of Origins under grant ANR10LABX66 and from CNRS under PICS06233. RBT received support from the NASA Astrophysics Data Analysis Program.
References
 Abazajian et al. (2003) Abazajian K. et al., 2003, AJ, 126, 2081
 Abazajian et al. (2009) Abazajian K. N. et al., 2009, ApJS, 182, 543
 Alimi et al. (2012) Alimi J.M. et al., 2012, ArXiv eprints: 1206.2838
 Angulo et al. (2012) Angulo R. E., Springel V., White S. D. M., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 426, 2046
 Bertschinger et al. (1990) Bertschinger E., Dekel A., Faber S. M., Dressler A., Burstein D., 1990, ApJ, 364, 370
 Colless et al. (2001) Colless M., Saglia R. P., Burstein D., Davies R. L., McMahan R. K., Wegner G., 2001, MNRAS, 321, 277
 Dekel (1994) Dekel A., 1994, ARA&A, 32, 371
 Dekel et al. (1990) Dekel A., Bertschinger E., Faber S. M., 1990, ApJ, 364, 349
 Doumler et al. (2013a) Doumler T., Courtois H., Gottlöber S., Hoffman Y., 2013a, MNRAS, 430, 902
 Doumler et al. (2013b) Doumler T., Gottlöber S., Hoffman Y., Courtois H., 2013b, MNRAS, 430, 912
 Doumler et al. (2013c) Doumler T., Hoffman Y., Courtois H., Gottlöber S., 2013c, MNRAS, 430, 888
 ForeroRomero & González (2015) ForeroRomero J. E., González R., 2015, ApJ, 799, 45
 Freedman et al. (2001) Freedman W. L. et al., 2001, ApJ, 553, 47
 Ganon & Hoffman (1993) Ganon G., Hoffman Y., 1993, ApJ, 415, L5
 GarrisonKimmel et al. (2014) GarrisonKimmel S., BoylanKolchin M., Bullock J. S., Lee K., 2014, MNRAS, 438, 2578
 Gavazzi et al. (2009) Gavazzi R., Adami C., Durret F., Cuillandre J.C., Ilbert O., Mazure A., Pelló R., Ulmer M. P., 2009, A&A, 498, L33
 Gottlöber et al. (2010) Gottlöber S., Hoffman Y., Yepes G., 2010, ArXiv eprints: 1005.2687
 Hahn et al. (2007) Hahn O., Carollo C. M., Porciani C., Dekel A., 2007, MNRAS, 381, 41
 Han (1992) Han M., 1992, ApJ, 395, 75
 Hendry & Simmons (1994) Hendry M. A., Simmons J. F. L., 1994, ApJ, 435, 515
 Heß et al. (2013) Heß S., Kitaura F.S., Gottlöber S., 2013, MNRAS, 435, 2065
 Hoffman et al. (2012) Hoffman Y., Metuki O., Yepes G., Gottlöber S., ForeroRomero J. E., Libeskind N. I., Knebe A., 2012, MNRAS, 425, 2049
 Hoffman & Ribak (1991) Hoffman Y., Ribak E., 1991, ApJ, 380, L5
 Hoffman & Ribak (1992) Hoffman Y., Ribak E., 1992, ApJ, 384, 448
 Huchra et al. (2012) Huchra J. P. et al., 2012, ApJS, 199, 26
 Hudson (1994) Hudson M. J., 1994, MNRAS, 266, 468
 Jasche & Wandelt (2013) Jasche J., Wandelt B. D., 2013, MNRAS, 432, 894
 Jha et al. (2007) Jha S., Riess A. G., Kirshner R. P., 2007, ApJ, 659, 122
 Karachentsev et al. (2004) Karachentsev I. D., Karachentseva V. E., Huchtmeier W. K., Makarov D. I., 2004, AJ, 127, 2031
 Karachentsev & Nasonova (2010) Karachentsev I. D., Nasonova O. G., 2010, MNRAS, 405, 1075
 Kitaura (2013) Kitaura F.S., 2013, MNRAS, 429, L84
 Klypin et al. (2003) Klypin A., Hoffman Y., Kravtsov A. V., Gottlöber S., 2003, ApJ, 596, 19
 Klypin et al. (2014) Klypin A., Yepes G., Gottlöber S., Prada F., Hess S., 2014, ArXiv eprints: 1411.4001
 Klypin et al. (2011) Klypin A. A., TrujilloGomez S., Primack J., 2011, ApJ, 740, 102
 Knollmann & Knebe (2009) Knollmann S. R., Knebe A., 2009, ApJS, 182, 608
 Kolatt et al. (1996) Kolatt T., Dekel A., Ganon G., Willick J. A., 1996, ApJ, 458, 419
 KraanKorteweg et al. (1994) KraanKorteweg R. C., Cayette V., Balkowski C., Fairall A. P., Henning P. A., 1994, 67, 99
 Kravtsov et al. (2002) Kravtsov A. V., Klypin A., Hoffman Y., 2002, ApJ, 571, 563
 Landy & Szalay (1992) Landy S. D., Szalay A. S., 1992, ApJ, 391, 494
 Lavaux & Wandelt (2010) Lavaux G., Wandelt B. D., 2010, MNRAS, 403, 1392
 Lee et al. (1993) Lee M. G., Freedman W. L., Madore B. F., 1993, ApJ, 417, 553
 Libeskind et al. (2012) Libeskind N. I., Hoffman Y., Knebe A., Steinmetz M., Gottlöber S., Metuki O., Yepes G., 2012, MNRAS, 421, L137
 Ludlow & Porciani (2011) Ludlow A. D., Porciani C., 2011, MNRAS, 413, 1961
 LyndenBell et al. (1988) LyndenBell D., Faber S. M., Burstein D., Davies R. L., Dressler A., Terlevich R. J., Wegner G., 1988, ApJ, 326, 19
 Murray et al. (2013) Murray S. G., Power C., Robotham A. S. G., 2013, Astronomy and Computing, 3, 23
 Nusser & Dekel (1992) Nusser A., Dekel A., 1992, ApJ, 391, 443
 Nusser et al. (1991) Nusser A., Dekel A., Bertschinger E., Blumenthal G. R., 1991, ApJ, 379, 6
 Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A16
 Prada et al. (2012) Prada F., Klypin A. A., Cuesta A. J., BetancortRijo J. E., Primack J., 2012, MNRAS, 423, 3018
 Sandage (1994) Sandage A., 1994, ApJ, 430, 1
 Skillman et al. (2014) Skillman S. W., Warren M. S., Turk M. J., Wechsler R. H., Holz D. E., Sutter P. M., 2014, ArXiv eprints: 1407.2600
 Sorce (2015) Sorce J. G., 2015, MNRAS, 450, 2644
 Sorce et al. (2014a) Sorce J. G., Courtois H. M., Gottlöber S., Hoffman Y., Tully R. B., 2014a, MNRAS, 437, 3586
 Sorce et al. (2013) Sorce J. G. et al., 2013, ApJ, 765, 94
 Sorce et al. (2014b) Sorce J. G., Tully R. B., Courtois H. M., Jarrett T. H., Neill J. D., Shaya E. J., 2014b, MNRAS, 444, 527
 Springel (2005) Springel V., 2005, MNRAS, 364, 1105
 Stoughton et al. (2002) Stoughton C. et al., 2002, AJ, 123, 485
 Strauss & Willick (1995) Strauss M. A., Willick J. A., 1995, Physics Reports, 261, 271
 Teerikorpi (1990) Teerikorpi P., 1990, A&A, 234, 1
 Teerikorpi (1993) Teerikorpi P., 1993, A&A, 280, 443
 Teerikorpi (1995) Teerikorpi P., 1995, Astrophysical Letters and Communications, 31, 263
 Teerikorpi (1997) Teerikorpi P., 1997, ARA&A, 35, 101
 Tinker et al. (2008) Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709
 Tonry et al. (2001) Tonry J. L., Dressler A., Blakeslee J. P., Ajhar E. A., Fletcher A. B., Luppino G. A., Metzger M. R., Moore C. B., 2001, ApJ, 546, 681
 Tully (2010) Tully R., 2010, ArXiv eprints: 1010.3787
 Tully (2015a) Tully R. B., 2015a, AJ, 149, 54
 Tully (2015b) Tully R. B., 2015b, AJ, 149, 171
 Tully et al. (2014) Tully R. B., Courtois H., Hoffman Y., Pomarède D., 2014, Nature, 513, 71
 Tully & Courtois (2012) Tully R. B., Courtois H. M., 2012, ApJ, 749, 78
 Tully et al. (2013) Tully R. B. et al., 2013, AJ, 146, 86
 Tully & Fisher (1977) Tully R. B., Fisher J. R., 1977, A&A, 54, 661
 Tully et al. (2008) Tully R. B., Shaya E. J., Karachentsev I. D., Courtois H. M., Kocevski D. D., Rizzi L., Peel A., 2008, ApJ, 676, 184
 Wang et al. (2014) Wang H., Mo H. J., Yang X., Jing Y. P., Lin W. P., 2014, ApJ, 794, 94
 Wang et al. (2013) Wang H., Mo H. J., Yang X., van den Bosch F. C., 2013, ApJ, 772, 63
 Watson et al. (2014) Watson W. A., Iliev I. T., Diego J. M., Gottlöber S., Knebe A., MartínezGonzález E., Yepes G., 2014, MNRAS, 437, 3776
 Willick (1994) Willick J. A., 1994, ApJS, 92, 1
 Willick et al. (1997) Willick J. A., Courteau S., Faber S. M., Burstein D., Dekel A., Strauss M. A., 1997, ApJS, 109, 333
 Yepes et al. (2014) Yepes G., Gottlöber S., Hoffman Y., 2014, New Astr. Rev., 58, 1
 Zaroubi et al. (1999) Zaroubi S., Hoffman Y., Dekel A., 1999, ApJ, 520, 413
 Zaroubi et al. (1995) Zaroubi S., Hoffman Y., Fisher K. B., Lahav O., 1995, ApJ, 449, 446
Appendix
In this appendix, we investigate the effect of the constrained boxsize on the power spectrum at large scales (small wavelengths). To do so, we produce 600 random and 600 constrained initial conditions at low resolution (grid size 128). One third of these initial conditions is built within a 250 Mpc boxsize, one third within 500 Mpc and the last third is built within a 1 Gpc boxsize. Values of the power spectra of all these initial conditions are then derived at different wavelengths (0.025, 0.05 and 0.075 Mpc) and their distributions are plotted as histograms on Figure 10. From left to right, the wavelength increases (or the scale decreases from 250 to 80 Mpc) and the dotted line gives the value of the Planck power spectrum at the different wavelengths respectively. We observe that constraining, indeed, tend to decrease the values of the power spectrum (solid lines versus dashed lines) at large scales. As expected, the difference between the power spectrum values of the random and the constrained initial conditions tend to be erased with the increase in boxsize (from black to red), i.e. when smaller and smaller regions of the box are constrained.