The phasespace structure of nearby dark matter as constrained by the SDSS
Abstract
Previous studies using numerical simulations have demonstrated that the shape of the cosmic web can be described by studying the Lagrangian displacement field. We extend these analyses, showing that it is now possible to perform a Lagrangian description of cosmic structure in the nearby Universe based on largescale structure observations. Building upon recent Bayesian largescale inference of initial conditions, we present a cosmographic analysis of the dark matter distribution and its evolution, referred to as the dark matter phasespace sheet, in the nearby universe as probed by the Sloan Digital Sky Survey main galaxy sample. We consider its stretchings and foldings using a tetrahedral tessellation of the Lagrangian lattice. The method provides extremely accurate estimates of nearby density and velocity fields, even in regions of low galaxy density. It also measures the number of matter streams, and the deformation and parity reversals of fluid elements, which were previously thought inaccessible using observations. We illustrate the approach by showing the phasespace structure of known objects of the nearby Universe such as the Sloan Great Wall, the Coma cluster and the Boötes void. We dissect cosmic structures into four distinct components (voids, sheets, filaments, and clusters), using the Lagrangian classifiers diva, origami, and a new scheme which we introduce and call lich. Because these classifiers use information other than the sheer local density, identified structures explicitly carry physical information about their formation history. Accessing the phasespace structure of dark matter in galaxy surveys opens the way for new confrontations of observational data and theoretical models. We have made our data products publicly available.
I Introduction
The accurate description of largescale structure formation, from which the cosmic web originates (Bond, Kofman & Pogosyan, 1996), is one of the major goals of modern cosmology. Traditionally, the Lagrangian and Eulerian descriptions of a discretized fluid of collisionless dark matter are comparatively discussed. In the Eulerian picture, quantities are tracked at fixed positions as particles move. In the Lagrangian framework, the motion and evolution of individual elements is followed. Therefore, particles can be thought of not only as mass tracers, but also as elements of a threedimensional dark matter sheet, which stretches, folds and distorts in a sixdimensional velocityposition phase space.
Using cosmological simulations, several recent papers have demonstrated the power of working within the dark matter sheet. Tessellations of the initial Lagrangian space allow the identification of stream crossings and improve density measurements (Abel, Hahn & Kaehler, 2012; Shandarin, Habib & Heitmann, 2012), provide accurate estimates of cosmic velocity fields (Hahn, Angulo & Abel, 2015), and yield a new approach for simulating dark matter (Hahn, Abel & Kaehler, 2013). In addition, the study of caustics (Arnold, Shandarin & Zeldovich, 1982; White & Vogelsberger, 2009; Vogelsberger & White, 2011; Hidding, Shandarin & van de Weygaert, 2014; Feldbrugge, Hidding & van de Weygaert, 2014), the twodimensional edges along which the dark matter sheet folds, brings valuable insights into structure formation (Neyrinck, 2012). Strong observational motivations exist for a Lagrangian understanding of the cosmic web. These include studying the dependence of galaxy properties on their evolving environment (e.g. Blanton et al., 2005; Libeskind et al., 2015), probing the effect of the dynamic largescale structure on the cosmic microwave background (e.g. Lavaux, Afshordi & Hudson, 2013; Planck Collaboration, 2014), testing the standard generalrelativistic picture of gravitational instability (Falck, Koyama & Zhao, 2015), and identifying the most promising regions to observe a possible dark matter annihilation signal (e.g. Gao et al., 2012).
Unfortunately, the Lagrangian picture cannot be readily applied to observational data. Difficulties are of two sorts: (i) the need for reconstructions of the dark matter density from galaxies, accounting for typical observational effects such as biases, survey mask and selection functions; and (ii) the need for corresponding initial conditions, which requires the inclusion of a physical picture of structure formation into the data model. For these reasons, most cosmic web analyses have been limited so far to simulations (see e.g. Cautun et al., 2014, for a recent study). However, recently, Jasche et al. (2010) and Leclercq, Jasche & Wandelt (2015b, hereafter borg sdss–Tweb) addressed the first issue. Using real data, they presented Eulerian classifications of cosmic environments in constrained realizations of the largescale structure. Importantly, they demonstrated capability of propagating observational uncertainties to cosmic web classification. Taking advantage of the inclusion of a physical model within the inference process, borg sdss–Tweb also presented the first probabilistic analysis of protostructures present in the initial density field. However, a complete investigation of the problem, including a probabilistic description of phasespace properties of dark matter based on real data, has not yet been performed in the literature.
This paper describes the first analysis of the Lagrangian dark matter sheet in the nearby Universe, as probed by galaxies of the northern galactic cap of the Sloan Digital Sky Survey (SDSS) main sample. Capitalizing on the full phasespace information, we produce highlydetailed and accurate maps of cosmic density and velocity fields, even in regions only poorly sampled by galaxies. We also analyze additional information specific to collisionless fluids, which was so far inaccessible in observations: the number of matter streams, the deformation and parity of fluid elements, and the presence of caustics. We present the first application of the Lagrangian cosmic web classifiers diva and origami to the real Universe, as well as a new algorithm, lich, which we design to take into account the heterogeneous (potential and vortical) nature of flows.
This chronocosmography project builds upon the inference of the past and present cosmic structure in the SDSS (Jasche, Leclercq & Wandelt, 2015, hereafter borg sdss–inference), performed using the borg algorithm (Jasche & Wandelt, 2013). It also relies on a large set of constrained realizations of the Universe, produced using the cola method (Tassev, Zaldarriaga & Eisenstein, 2013). These necessarily introduce some degree of extrapolation with respect to the original inference results, since they contain aspects of the physical model that have not been constrained by the data (in particular, in this work, vortical matter flows). However, these extrapolations are not arbitrary, but are methodical predictions that need to be consistent with the observational constraints. Indeed, a complex forward model naturally introduces correlations between its constrained and unconstrained aspects, which may yield reliable predictions for the physics that has not been captured by the simpler data model. Our methodology therefore predicts the properties of the Lagrangian dark matter sheet based on fusing prior information from simulations and data constraints. Such an approach has been previously demonstrated to accurately recover the properties of voids at the level of the dark matter distribution (Leclercq et al., 2015). Importantly, all of our predictions are fully probabilistic, which means that their degree of speculativeness is controlled and can be checked against unconstrained simulations. From a Bayesian perspective, all the results presented in this work can then be used as prior information for updating our knowledge in light of additional observations, or for making optimal decisions in the presence of uncertainty (see Leclercq, Jasche & Wandelt, 2015a).
The structure of this paper is as follows. We start by describing our tools for analyzing the dark matter sheet in section II. This includes a review of the Lagrangian displacement field, the description of the existing Lagrangian classifiers diva and origami and the introduction of a new one that we call lich, and the definition of estimators based on tessellating the Lagrangian particle grid. In section III, we then infer, characterize and classify the SDSS volume in Lagrangian space. Finally, in section IV, we show the consequences of these results in Eulerian coordinates and demonstrate that resulting maps of the redshiftzero cosmic web have an interesting level of predictive power. We summarize our study and conclude in section V.
Ii Description of the Lagrangian sheet
ii.1 The Lagrangian displacement field
The fundamental object studied in this paper is the displacement field , which maps the initial (Lagrangian) position of particles q to their final (Eulerian) position (see e.g. Bernardeau et al., 2002, for a classic review):
(1) 
As any smooth vector field, it is possible to split into a scalar and a rotational part (the Helmholtz decomposition, Chan, 2014),
(2) 
where is the scalar potential and is the vector potential, which obey the Poisson equations:
(3)  
(4) 
From analytic approaches and numerical simulations, it has been known for some time that the displacement field is almost fully potential; in particular, contains all the information up to second order in Lagrangian perturbation theory (2LPT). In this paper, we estimate and by computing respectively and , where is a fast Fourier transform (FFT) on the Lagrangian mesh on which the displacement field is defined. For later use, we define .
We denote by the Jacobian determinant of the transformation between Lagrangian and Eulerian coordinates,
(5) 
where the deformation tensor can be written as the identity tensor plus the shear of the displacement field . The characteristic polynomial of is written as
(6) 
where the are the three eigenvalues (real or complex) of (i.e. the solutions of ). The Jacobian is then .
The characteristic equation can be written
(7) 
where we call the the Lagrangian invariants (in analogy to the rotational invariants introduced by Wang et al., 2014, who use the shear of the Eulerian velocity field instead of the Lagrangian displacement field). They are given in terms of the eigenvalues by (Chong, Perry & Cantwell, 1990)
(8)  
(9)  
(10) 
from which one can derive the explicit formulation in terms of the coefficients of (summation over repeated indices is implied):
(11)  
In the previous equations, we have used the decomposition of into a symmetric part (the rate of deformation tensor) and an antisymmetric part (the spin tensor),
(12) 
with
(13)  
(14) 
Physically, denotes the rate of deformation of a Lagrangian fluid element and the rate of rotation.
The surface in Lagrangian invariants space that divides real and complex solutions for is given by the equation
(15) 
Using the discriminant of this quadratic equation, one finds two real solutions and () if and only if . Therefore, at fixed , the regions where all the eigenvalues of are real correspond to
(16) 
In this case, the displacement field is purely potential and we denote the eigenvalues by . We also define the ellipticity and prolateness (e.g. Peacock & Heavens, 1985; Bardeen et al., 1986) as
(17)  
(18) 
Ellipsoids are prolatelike if or oblatelike if . The limiting cases are for prolate spheroids and for oblate spheroids. In the following, we will also consider and , for which the boundary conditions are still valid.
After orbitcrossing, an antisymmetric part of is generated and the eigenvalues of are no longer all real. We parametrize them as and where are real. The invariants then simply read (see equations (8)–(10))
(19)  
(20)  
(21) 
As we will show in section II.2, the Lagrangian invariants allow a complete classification of the cosmic web, including both potential and vortical flow patterns (as with the rotational invariants, see Chong, Perry & Cantwell, 1990; Wang et al., 2014).
Finally, we note that it is possible to translate the Lagrangian invariants to the eigenvalues , using Vieta’s formulas for cubic equations:
(22)  
(23)  
(24) 
and for ,
(25) 
ii.2 Lagrangian webtype classifications
In our earlier analysis of the volume covered by the SDSS main sample galaxies (borg sdss–Tweb), we adopted the Eulerian cosmic web classifier known as the Tweb (Hahn et al., 2007; ForeroRomero et al., 2009). In this framework, structures are classified according to the sign of the eigenvalues of the tidal field tensor , defined as the Hessian of the rescaled gravitational potential :
(26) 
where obeys the reduced Poisson equation
(27) 
being the local density contrast. A void point corresponds to no positive eigenvalue, a sheet to one, a filament to two, and a cluster to three positive eigenvalues.
Related Hessianbased techniques are the Vweb (Hoffman et al., 2012) and classic (Kitaura & Angulo, 2012). More broadly, Eulerian cosmic web algorithms include topological (AragónCalvo et al., 2010; Sousbie, 2011), stochastic (González & Padilla, 2010; Tempel et al., 2014), and multiscale field (AragónCalvo et al., 2007; Cautun, van de Weygaert & Jones, 2013; AragonCalvo & Yang, 2014) methods, all of which yield different physical insights from the phasespace analysis subject of this paper.
In this section, we discuss three alternative methods that can be used for producing Lagrangian classifications of the cosmic web: diva, lich and origami.
ii.2.1 Diva
In analogy with the Zel’dovich (1970) “pancake” theory, Lavaux & Wandelt (2010) propose to define cosmic structures using , the eigenvalues of the shear of the displacement field, . In the diva (DynamIcal Void Analysis) scheme, a particle is defined as belonging to a Lagrangian void, sheet, filament or cluster, if, respectively, zero, one, two or three of the are negative. In order to deal only with real eigenvalues, diva assumes that the flow is potential, i.e. that is equal to its symmetric part .
The diva and Tweb methods share strong similarity in their formalism. However, a major difference is that the Tweb operates on voxels of the discretized domain, whereas diva operates on Lagrangian patches (numerically approximated by particles). As shown in borg sdss–Tweb (sections III and IV), the Tweb can be used to classify both earlytime and latetime structures; however the Tweb classification of earlytime structures is performed at the level of voxels of the discretized initial density field, whereas the diva classification of Lagrangian structures is performed at the level of particles of the dark matter sheet.
ii.2.2 Lich
This section introduces a new cosmicweb classification procedure based on the Lagrangian invariants, which we call lich (Lagrangian Invariants Classification of Heterogeneous flows). As we will show, lich is a generalization of diva to the case of heteregenous (i.e. potential and vortical) flows. lich uses for the Lagrangian invariants the same formalism as introduced by Wang et al. (2014) for the rotational invariants; the reader is therefore referred to this paper for more details.
Let us first consider the case of a potential flow. The real solution condition, equation (16), must be fulfilled. The case where all the are negative (a diva cluster) further corresponds to for all (see equations (8)–(10)). A lich potential cluster is therefore defined by
(28) 
Similarly, a lich potential void, which generalizes a diva void (all positive), corresponds to
(29) 
Two regions in Lagrangian invariants space correspond to a diva sheet (, ):
(30)  
(31) 
They define the lich potential sheets. Two regions also correspond to a diva filament (, ):
(32)  
(33) 
They define the lich potential filaments.
The case of a vortical flow corresponds to either (in which case there is no solution to equation (15)), either . In these regions in Lagrangian invariants space, an additional frontier is important: (i.e. ), which corresponds to the plane
(34) 
The remaining regions can then be labeled as vortical voids, vortical sheets, vortical filaments or vortical clusters. Figure 1 summarizes the lich classification procedure in Lagrangian invariants space. For the detailed criteria defining all the regions there, the reader is referred to Appendix A in Wang et al. (2014).
In section III.1.2, we will also use a simplified version of the full lich procedure, where particles fall into one of the usual four categories (void, sheet, filament or cluster) by simply ignoring the potential/vortical subdivision. Note that this simplified lich procedure is different from diva, since it does not use the initial assumption that the flow is potential.
ii.2.3 Origami
Another way of classifying particles is to consider them as vertices of an initially regular grid, distorted by gravity as the largescale structure forms. Shellcrossing happens when Lagrangian cells collapse or invert. A particle’s origami (OrderReversIng Gravity, Apprehended Mangling Indices) morphology is defined by the number of orthogonal axes along which shellcrossing has occurred (Knebe et al., 2011; Falck, Neyrinck & Szalay, 2012). More precisely, void, sheet, filament, and cluster particles are defined as particles that have been crossed along zero, one, two, and three orthogonal axes, respectively.
In practice, the computer implementation of origami works as follows: a particle has been crossed along an axis if there exists a particle in the same row of the Lagrangian grid, such that and have opposite sign. Such particlecrossing are checked along the three orthogonal axes of the original Cartesian grid, as well as three triplets of orthogonal axes, obtained by rotating the original grid by 45 along one of its axes. For more details, the reader is referred to section 2.1 in Falck, Neyrinck & Szalay (2012).
ii.3 Tessellation of the Lagrangian grid
Abel, Hahn & Kaehler (2012) and Shandarin, Habib & Heitmann (2012) proposed new techniques to trace the dark matter phasespace sheet and demonstrated that they give a wealth of information previously unavailable in simulations. In this paper, we will demonstrate that these techniques can now be used in dataconstrained realizations of the nearby Universe. The core idea is to use a tessellation of the uniform Lagrangian lattice of particles. Each fundamental cube is decomposed into tetrahedral cells that use the particles as their vertices. These cells are tracked in the final conditions, using the particle IDs which are setup so that they encode the initial position of the particle.
Several possible decompositions of the unit cube exist: for example, Shandarin, Habib & Heitmann (2012) use five tetrahedra of two different volumes. In this work, however, we will follow Abel, Hahn & Kaehler (2012); Hahn, Angulo & Abel (2015) and use the Delaunay decomposition into six identical tetrahedra (see figure 2). Specifically, the connectivity list is , , , , and where the unit cube vertices are labeled as in figure 2.
Initially, the volume of each tetrahedron is where is the side length of the fundamental cube. In the final state, the same tetrahedron has a volume given by a determinant (or crossproduct)
(35) 
where A, B, C, D label the final positions of the tetrahedron vertices. Note that the volume can be negative if the tetrahedron has been turned insideout. In fact, the number of sign reversals of an elementary Lagrangian volume defines the flipflop field (e.g. Neyrinck, 2012; Shandarin & Medvedev, 2016). Each tetrahedron carries 1/6 of the mass of a particle , the mass density associated is therefore . At any particle position, there are 24 tetrahedra that share this vertex (, see figure 2).
Tessellating the dark matter phasespace sheet allows one to define a new quantity on the Lagrangian grid, the primordial stream density,
(36) 
where the summation runs on the 24 tetrahedra that share the particle’s initial position q as one of their vertices. Altogether, these tetrahedra have an initial volume of .
In addition to the primordial stream density, tracking the final configuration of the tessellated dark matter phasespace sheet allows one to perform the Lagrangian transport of any quantity defined on a particle basis, which can be used to define a variety of estimators on the Eulerian grid. This is achieved by finding all the tetrahedra that contain the point x at which properties are to be determined. Each such tetrahedron is given four coefficients , , , for the four particles that constitute its vertices (A, B, C, D). Fields are interpolated inside each tetrahedron using Shepard’s method (inverse distancesquared weighting); namely, each tetrahedron deposits to the point x the quantity
(37) 
where and .
In this work, we consider the following coefficients and associated estimators:

simply yields the number of tetrahedra enclosing point x. Since we start from a complete tessellation of the considered volume, this is also the number of matter streams contributing at point x, which we call the Eulerian secondary stream density, following Abel, Hahn & Kaehler (2012).

yields an estimator of the final Eulerian density field (see Abel, Hahn & Kaehler, 2012), which we note .

if particle P is of type , 0 otherwise yields an estimator of the density contained in structure of type , i.e. . The field represents the probability for a final Eulerian point to belong to structure type . This method can therefore be used to translate Lagrangian classifications into Eulerian coordinates. As it is computationally expensive, we will introduce a simplified alternative in section IV.2.

where is the velocity of particle P along the direction yields an estimator of the component of the massweighted velocity field. The resulting field divided by yields the velocity field estimated from the dark matter phasespace sheet (see Hahn, Angulo & Abel, 2015).
Instead of explicitly projecting all the tetrahedra to the Eulerian grid, it is possible to approximate each of them by four pseudoparticles (in red in figure 2; see the discussion in Hahn, Abel & Kaehler, 2013, section 2.4, for where to place these pseudoparticles). In this context, the original particles are called the flow tracers and the pseudoparticles are called the mass tracers. Flow tracers uniquely determine the position of mass tracers, which can be deposited to the Eulerian grid by standard interpolation techniques such as the cloudincell (CiC, Hockney & Eastwood, 1981) algorithm.
Iii Lagrangian analysis of the SDSS volume
As mentioned in the introduction, this paper is part of a series following the recent application of the Bayesian inference framework borg to the SDSS main sample galaxies (borg sdss–inference). We further rely on the set of 1,097 constrained simulations presented in borg sdss–Tweb. These realizations are obtained, starting from borg initial conditions, by nonlinear filtering with cola (Tassev, Zaldarriaga & Eisenstein, 2013). The initial density field, defined on a voxel grid with side length of Mpc/, is populated by dark matter particles placed on a regular Lagrangian lattice. Particles are evolved with 2LPT to the redshift of then with 30 cola timesteps from to . The displacement field is obtained on the Lagrangian lattice by subtracting the final position of particles from their initial position, and any necessary quantity involving derivatives is computed using FFTs on the same lattice.
In this section, we first analyze the Lagrangian properties of individual samples (section III.1), then show how the full set of constrained samples can be used to propagate uncertainties to all inferred quantities (section III.2).
iii.1 Analysis of individual samples
iii.1.1 Stream properties
The first row of figure 3 shows the stream properties in one of the samples of our set. The slices show a region of the celestial equatorial plane encompassing the Sloan Great Wall (Gott et al., 2005; Einasto et al., 2011).^{1}^{1}1In the plots of this paper, we kept the coordinate system of borg sdss–inference. The flow is characterized in terms of the divergence of the displacement field , the norm of its vector part , and the primordial stream density (see equation (36)). quantifies the angleaverage stretching and distortion of the dark matter sheet (Neyrinck, 2013) and contains all the information on as long as the flow remains potential. can be used to quantify how much vorticity is generated when structure formation enters the multistream regime, and the primordial stream density counts these streams.
It is important to note that, since the displacement field is curlfree up to second order in Lagrangian perturbation theory (see e.g. Bernardeau et al., 2002), and since borg relies on 2LPT as a proxy for gravitational dynamics (Jasche & Wandelt, 2013), any vorticity present in our largescale structure realizations is not explicitly constrained by the data. The field therefore corresponds one plausible realization, compatible both with the constrained initial conditions and with the nonlinear structure formation model (cola).
Our constrained samples exhibit the same behavior as previously observed in simulations. Specifically, regions of high primordial stream density nicely correspond to Lagrangian clusters, the “lakes” where . This is also where a rotational part of the displacement field is generated. Conversely, singlestream regions correspond to the “mountains” where . From the behavior of the field (limited at most to Mpc/ in the regions of highest stream density), we also note that we can take advantage of the absence of substantial vorticity on large scales. This is consistent with the results of Chan (2014): at and up to Mpc/, the contribution to the power spectrum of coming from its curl component remains negligible compared to the nonlinear evolution of its scalar part. Therefore, for the sake of our analysis (except when explicitly measuring vorticity), it will be safe to assume that all the information on is contained in its divergence , or equivalently, that .
iii.1.2 Lagrangian classifications
We then assign a Lagrangian structure type to each particle, according to the diva, lich, and origami procedures. Lavaux & Wandelt (2010) discuss the possibility of smoothing the displacement field at different Lagrangian scales, in order to probe the hierarchical cosmic web. This allows the description of structure formation at different dynamical epochs, each Lagrangian scale corresponding to a different collapse time. In contrast, explicit Eulerian multiscale formalisms (AragónCalvo et al., 2007; Cautun, van de Weygaert & Jones, 2013; AragónCalvo et al., 2010; AragonCalvo & Szalay, 2013) probe the hierarchy of structures at a given time. The adhesion approximation (Gurbatov, Saichev & Shandarin, 1989; Weinberg & Gunn, 1990; Kofman, Pogosian & Shandarin, 1990; Kofman et al., 1992; Gurbatov, Saichev & Shandarin, 2012; Hidding, van de Weygaert & Shandarin, 2016) provides a Lagrangian description of the hierarchical evolution of structures, but it is only correct in the linear and weakly nonlinear regime, contrary to the numerical cola model used in this work. The main effect of smoothing the Lagrangian grid is to change the unconditional number of particles in different structures; data constraints do not dramatically alter this result (Lavaux & Wandelt, 2010). Therefore, while this would be perfectly feasible in our framework, we choose to focus our analysis on the scale defined by our setup, without an additional smoothing step. This choice corresponds to a resolution of Mpc/ for the Lagrangian grid ( particles in a Mpc/ cubic box) and Mpc/ for discretized Eulerian maps (on a voxel grid).
The second row of figure 3 represents Lagrangian structures on the grid of particles, in one of our constrained samples. From left to right, the panels show the diva, lich, and origami morphologies of the dark matter sheet. As expected, all the classifications correlate well with the stream properties shown in the first panel (in particular with ). Clusters (in red) correspond to the “lakes” of high primordial stream density where , while less complex structures are regions with fewer streams and higher values of . Singlestream regions, which correspond to origami void particles, form a percolating system, as noted by Falck & Neyrinck (2015).
The difference between diva and origami classifications (notable primarily for voids, then for sheets and filaments) is caused by the different phasespace criterion used for defining structures: particlecrossings that have already happened for origami, versus those nonlocally predicted from the deformation tensor for diva. A large singlestream void percolates with origami, while filament and sheet regions are limited: this reflects the current configuration of the cosmic web. Conversely, with diva, particles around filaments and sheets typically get the signature of these structures, on the basis of the dynamical trend.
As shown in section II.2.2, lich is a generalization of diva to the case of heterogeneous (potential and vortical) flows — applying lich under the assumption of a potential flow yields diva. Thus, the comparison between diva and lich allows to check this assumption. Comparing the diva, lich, and panels shows that the two classifications differ in the regions where vorticity is generated in . Specifically, with lich, most cluster regions are vortical ( of cluster particles), except the smallest ones, and surrounded by a filamentary shell. of filament particles are vortical, in the surroundings clusters; the rest are potential. Sheet regions are potential ( of particles). Void regions are mostly potential ( of particles), but with some spinning regions. When ignoring the subclassification into potential and vortical structures, but without assuming a potential flow, the simplified lich procedure agrees with diva for of particles, which is an additional argument for the accuracy of the potential flow approximation. The difference is mostly driven by clusters and voids, with of lich cluster and of lich void particles misidentified by diva (almost all of which as filaments and sheets, respectively).
iii.2 Uncertainty quantification
As in our previous work (Jasche, Leclercq & Wandelt, 2015; Leclercq, Jasche & Wandelt, 2015b; Lavaux & Jasche, 2016), the variation among constrained samples constitutes a Bayesian quantification of uncertainties (coming in particular from the survey mask and selection effects). This is illustrated in figure 4, where we show the particlewise posterior mean of different quantities on the Lagrangian grid. The first row shows the eigenvalues of the rate of deformation tensor (the symmetric part of the shear of the displacement field ). These are used by diva to classify structures. The second row shows derived quantities: the (reduced) ellipticity and prolateness , and the Jacobian of the transformation from Lagrangian to Eulerian coordinates, . For each particle, it is the volume of the associated fluid element times its parity (Neyrinck, 2012; Shandarin & Medvedev, 2016). White/blue particles have the same parity as in the initial conditions, while black/orange particles have swapped parity. Caustics are clearly visible wherever blue and orange regions are juxtaposed, for example in the patches that collapse to form halos and in the shellcrossed walls surrounding voids. In some clusters, several parity inversions can be observed. The last row shows the posterior mean for the Lagrangian invariants , , , which are used by lich to classify structures. These cosmographic results constitute the first characterization of the nearby Lagrangian dark matter sheet, and contain a wealth of information on the dynamic largescale environment of SDSS main sample galaxies. They come as fully probabilistic, which means that uncertainties are quantified at the level of the map.
Uncertainty quantification can also be selfconsistently propagated to Lagrangian structure type classification, as we now show. Let us denote by one of the Lagrangian classifiers (diva, lich, or origami). In a specific realization, uniquely characterizes the Lagrangian grid, meaning that it provides particlewise scalar fields that obey the following conditions for each particle :
(38) 
where for diva and origami ( void, sheet, filament, cluster); for lich ( potential void, potential sheet, potential filament, potential cluster, vortical void, vortical sheet, vortical filament, vortical cluster). By applying to the complete set of borgcola realizations, we can quantify the degree of belief in structure type classification in terms of particlewise scalar fields that obey the following conditions for each particle :
(39) 
Here, are the posterior probability distribution functions (pdfs) for particle to belong to structure , as defined by classifier , conditional on the data . These are estimated by counting the relative frequencies at each individual particle position within the set of realizations (see also Jasche et al., 2010; borg sdss–Tweb). Therefore, the cosmic webtype pdf mean on the Lagrangian grid is given by
(40) 
where labels one of the samples, is the result of classifier on the th sample at particle position (i.e. a unit vector containing zeros except for one component, which indicates the structure type), and is a Kronecker symbol.
Using the above definitions, we obtain probabilistic maps of structures on the Lagrangian grid describing the early Universe in the volume probed by SDSS main sample galaxies. More precisely, for each classifier , we obtain a probability mass function (pmf) at each particle position on the Lagrangian grid, indicating the possibility to encounter a specific structure type at that position. These pmfs consist of numbers in the range , , that sum up to one for each particle. They are approximated by the mean of each structure pdf (see equation (40)).
In figure 5, we show slices through these webtype posterior probabilities. The rows represent the structures inferred using diva, lich, and origami, respectively. For comparison, the reader is referred to figure 6 in borg sdss–Tweb, which shows structures inferred in the initial density field using the Tweb. Note that the present results are defined on the Lagrangian grid of particles, whereas the Tweb result is defined on the voxel grid on which initial density fields are discretized. The plot shows the anticipated behavior, with values close to certainty (i.e. zero or one) in regions covered by data, while the unobserved regions (at high redshift or out of the survey boundary) approach a uniform value corresponding to the prior. It is worth noting that different classifiers can have very different prior probabilities for the same structure types (these are given in table II in Leclercq et al., 2016). The good agreement between diva and lich at the level of individual samples, coming from the absence of a substantial vector part in , is also correctly propagated to the posterior mean.
A voxelwise quantification of the information gain due to SDSS galaxies is provided by the relative entropy of corresponding posteriors and priors (borg sdss–Tweb; Leclercq et al., 2016). We applied this technique and found that the Lagrangian webtype posteriors presented in figure 5 yield much higher information gain than the application of the Tweb to initial density fields. Averaged over the constrained volume, this number is Sh for diva, Sh for lich, and Sh for origami, as compared to Sh for the Tweb. This consideration substantiates the nonlocal transport of information along Lagrangian trajectories operated by the structure formation model, and demonstrates that a probabilistic description of the Lagrangian dark matter sheet yields a much higher information content for structure type classification than previous methods.
Iv Consequences in Eulerian coordinates
In this section, we describe potentially observable consequences of our inference of the nearby Lagrangian dark matter phasespace sheet. To do so, we analyze the properties of our constrained simulations in Eulerian coordinates in light of previous Lagrangian results. In other words, we demonstrate that our analysis of dark matter phasespace properties has an interesting degree of predictive power.
Formally, the problem can be expressed as follows. Before the data are considered, the distribution of an unknown observable is its prior predictive distribution,
(41) 
where are the threedimensional initial conditions. After the data have been observed, the distribution for is the posterior predictive distribution,
(42)  
In the above equations, we have used the borg inference of initial conditions, i.e. where is one of the samples. We have also assumed that does not further depend on the data once the initial conditions are known, and that is a deterministic function of the initial conditions, i.e. . Therefore, under the assumption that inferred initial conditions contain all the information on , our procedure generates samples of .
In the following, we start by looking at phasespace properties traced at Eulerian coordinates, which yields improved estimators for density fields (section IV.1), then we show how to translate Lagrangian classifications to Eulerian coordinates (section IV.2), discuss the properties of cosmic velocity fields (section IV.3), and finally present some cosmographic results (section IV.4). For all results, we also show that observational uncertainties are propagated in a nonlinear and nonGaussian fashion by the cola forward model, from which we get a feeling for how precise the predictions are.
iv.1 Phasespace properties traced at Eulerian coordinates
iv.1.1 Eulerian secondary stream density
The simplest phasespace property that can be translated into Eulerian coordinates is the number of matter streams. As discussed in section II.3, this is simply achieved by counting, in the final configuration of the tessellated Lagrangian dark matter sheet, the number of tetrahedra that contain any given position. The resulting field is called the secondary stream density. Since the dark matter phasespace sheet is discretized using a finite number of particles (and hence of tetrahedra), the secondary stream density is only a lower bound on the actual number of dark matter streams at any single location.
We applied this estimator to our full set of constrained realizations. Figure 6 shows the resulting posterior mean for the Eulerian secondary stream density. A slice of 3 Mpc/ (one voxel) thickness around the celestial equator is shown, and the galaxies of the northern cap of the SDSS main galaxy sample that are contained in the same volume are overplotted.^{2}^{2}2The transformation between our set of Cartesian coordinates and equatorial J2000.0 spherical coordinates can be found in appendix B of borg sdss–inference.
As expected from simulations (e.g. Ramachandra & Shandarin, 2015), the number of streams correlates well with the galaxy density field: regions of high galaxy density are the ones where several matter flows contribute. For example, in the inner regions of the Sloan Great Wall, which is clearly visible in the middle of the slice, at least 6 matter streams are converging. At large distances ( Mpc/), comparing the posterior mean for the Eulerian secondary stream density and galaxy number counts reveals that the predictions remain accurate even in regions poorly sampled by galaxies. This comes from the capacity of the borg algorithm to accurately recover structures even in regimes of low signaltonoise (see the discussion in borg sdss–inference, ).
Our result readily provides a probability distribution for the number of largescale matter streams at any galaxy position. This information, which is not available from raw survey data, is of potential interest for various cosmological and astrophysical analyses. It is expected to be crucially linked to the formation history of galaxies. If dark matter annihilates, then clear signals are predicted to come from nearby rich galaxy clusters (Gao et al., 2012). In fact, if dark matter is absolutely cold, its mass density (and therefore the annihilation luminosity) should be infinite at caustics (White & Vogelsberger, 2009). Galaxies living in streamcrossed regions are also the ones for which the “reconstruction” procedure used to correct for largescale flows in baryon acoustic oscillations analyses (Eisenstein et al., 2007; Padmanabhan et al., 2012; Burden, Percival & Howlett, 2015) will fail, as it relies on the singlestream assumption and on the Zel’dovich approximation.
iv.1.2 Density estimators
As shown by Abel, Hahn & Kaehler (2012) and Hahn, Abel & Kaehler (2013) in simulations, instead of simply assigning dark matter particles to the grid via CiC deposit, knowledge of their initial positions can be exploited to give better density estimators. Specifically, the tessellated dark matter phasespace sheet yields a density estimator when the four vertices of each tetrahedron are given a weight of (see section II.3). As a simpler alternative, tetrahedral mass elements can be approximated by pseudoparticles (the mass tracers) which are assigned to the grid by standard CiC.
We applied these two new density estimators to our set of constrained simulations. The results are shown in figure 7 in comparison to earlier results. There, the first column represents the raw borg output, for which the structure formation model is 2LPT (borg sdss–inference). The second column represents the standard CiC assignment of usual dark matter particles (flow tracers) in our borgcola realizations, as obtained in borg sdss–Tweb. The third and fourth columns are the new density estimators: the third raw corresponds to the CiC assignment of mass tracers, and the fourth row to the explicit projection of all tetrahedra. For each density estimator, we show one sample (first row), the posterior mean (second row) and the posterior standard deviation (third row). Tetrahedra can be projected to an arbitrarily fine grid. Similarly, there are 24 times more mass tracers than flow tracers, but as argued by Hahn, Abel & Kaehler (2013), the pseudoparticle approximation technique yields more than a factor 24 increase in effective mass resolution, when comparing CiC deposit of flow tracers versus mass tracers. For these reasons, we have been able to double the resolution of the Eulerian grid: the third and fourth columns show a grid of voxels versus voxels for the first and second columns.
In individual samples, the plot shows the anticipated behavior, with visually indistinguishable observed and unobserved regions. This indicates the joint performance of borg, cola, and phasespace estimators to augment unobserved regions with physically and statistically correct information. The ensemble of dataconstrained realizations also permits to correctly propagate observational uncertainties. As expected, all mean density fields exhibit highlydetailed structures where data constraints are available, and approach the cosmic mean value in unobserved regions. The standard deviation reveals a high degree of correlation between signal and noise, which is typical of nonlinear inference problems. In particular, in cola realizations, the densest clusters which cannot be reproduced by 2LPT correspond to high variance regions. In some small patches away from the wellconstrained regions, the density estimated from explicit projection of tetrahedra has high posterior mean and variance. These are numerical artifacts driven by one or a few samples, and due to a bias of the estimator in the densest regions (the linearlyinterpolated phasespace sheet no longer tracks the true distribution function there, see Abel, Hahn & Kaehler, 2012).
The two new density estimators extend the nonlinear filtering procedure originally introduced by Leclercq et al. (2015) and reviewed in Leclercq (2015, chapter 7). Starting from inferred density fields, we go forward in time using an arbitrary structure formation model, , combined with a density estimator, . The result is a nonlinear function, noted , which maps the initial density field into a final density field at a scale factor :
(43) 
The borg data model uses 2LPT and {CiC of flow tracers} (Jasche & Wandelt, 2013). In Leclercq et al. (2015) and borg sdss–Tweb, we improved the physical model over 2LPT by using respectively Gadget2 (Springel, 2005) and cola (Tassev, Zaldarriaga & Eisenstein, 2013); but we kept {CiC of flow tracers}. Here, we keep the physical model fixed, cola, but we improve instead the density estimator, using {CiC of mass tracers} or {projection of tetrahedra}. As shown, this alternative approach also yields improvements of the composite function .
iv.2 Translation of Lagrangian classifications to Eulerian coordinates
In section III.1.2, we have considered classifications of the dark matter sheet, discretized by the simulation particles. In order to classify structures at redshift zero, it is possible to translate the results of Lagrangian classifiers into Eulerian coordinates. As mentioned in section II.3, a natural way to do this is to use the tessellation of the Lagrangian lattice into tetrahedra, which allows the translation of any particlewise quantity to the Eulerian grid. However, since this is rather computationally expensive, we introduce below a simplified procedure to translate Lagrangian classifications into Eulerian coordinates. The latter is less precise at the level of individual samples, but we have checked that it yields negligible differences with the respect to the tessellation approach, as far as the posterior mean (figure 8) is concerned.
The simplified approach to translate Lagrangian classifications into Eulerian coordinates works as follows. We consider particles at their final positions and assign them to the Eulerian grid (of voxels in our case) using the CiC scheme. For structure type , each particle gets a weight of if its Lagrangian classification is , and 0 otherwise. In each voxel, we then count the relative frequencies for each massweighted structure type. In this fashion, in each realization , we obtain probabilities, in the range , for a final Eulerian voxel to belong to each of the structure types. We denote them by . In general, we have ; this sum is if and only if, in realization , no particle was assigned to the cell by the CiC scheme (i.e. it is empty as well as its eight neighbors). These voxels are flagged as noisy and discarded.
It is then possible to go from these persample pdfs to a global pdf for the entire analysis, in a similar way as before: we count the relative frequencies of the probabilities at each Eulerian spatial position within the set of samples. With this definition, the webtype posterior mean is given by
(44) 
where labels one of the samples, and the normalization in voxel is
(45) 
In figure 8, we show slices through the redshiftzero webtype posterior mean for the different structures as defined by diva, lich, and origami. As before, the plot demonstrates that we are able to propagate typical observational uncertainty to structure type classification. Consistently with the results of section III.1.1, where we argued for the accuracy of the approximation of to a potential field, we find only small differences between diva voids, sheets, and filaments, and the corresponding lich potential structures. However, unsurprisingly, lich allows a more physical detection and characterization of clusters, where the potential flow assumption breaks down. At this point, the reader is referred to figure 3 in borg sdss–Tweb for our map of the same volume, inferred using the Tweb definition. The Tweb and diva maps are visually similar, with an overall smoother structure for the voids defined by diva, which are sharply separated by walls. diva and lich also visually excel in highlighting the filamentary structure of the cosmic web. In contrast, with origami, most of the volume is filled by voids and more complex, shellcrossed structures are rarely identified.
iv.3 Cosmic velocity fields in the SDSS volume
Understanding largescale velocity fields is of central importance for modern cosmology: it is linked to galaxy and largescale structure formation (Pichon & Bernardeau, 1999), is crucial in modeling redshiftspace distortions (Percival & White, 2009) and in measuring the kinetic SunyaevZel’dovich effect (Shao et al., 2011; Lavaux, Afshordi & Hudson, 2013). Unfortunately, determining the continuous velocity field is an intrinsically complicated problem, even in simulations, because (i) it is only sampled at discrete locations, with poor resolution in lowdensity regions; and (ii) it is multivalued in regions where multiple streams contribute. Typical observational uncertainties (mask, selection effects, galaxy biases) add an additional layer of complexity. Recently, Hahn, Angulo & Abel (2015) showed that the first issues can be overcome in simulations by using the Lagrangian phasespace sheet. We apply this technique to our set of constrained simulations and demonstrate that we can also circumvent the observational complications. In doing so, we provide unprecedentedly accurate reconstruction of cosmic velocity fields in the volume of the nearby Universe covered by the northern galactic cap of the SDSS main galaxy sample and its surroundings.
We applied the estimators described in section II.3 to our set of constrained simulations, and obtained the components of the velocity field , , , on Eulerian grids of and voxels. As shown in the beginning of section IV, under the assumption that velocities are conditionally independent of the data once initial conditions are known, these samples are drawn from . As before, each of them gives one plausible realization of the velocity field, and the variation between samples quantifies the remaining uncertainty. In figure 9, we show the posterior mean for the three spherical components of the velocity field, in the celestial equator. The interested reader may want to compare with figure 7 in borg sdss–inference, which showed the velocity field as probed by 2LPT particles in one particular sample. The present work constitutes a substantial improvement over these results, since the velocity field is now predicted by a fully nonlinear structure formation model instead of 2LPT, wellsampled even in lowdensity regions, and with more straightforward uncertainty quantification. The recovered velocity field is fully compatible with the known structure of the nearby Universe; showing, for example, the infall of matter onto the Sloan Great Wall.
As noted in borg sdss–inference, velocity fields are derived a posteriori with borg, i.e. they are a prediction of physical models given the inferred initial density field.^{3}^{3}3In contrast, algorithms like virbius (Lavaux, 2016) explicitly exploit velocity information contained in the data to constrain the result, and reconstructions based on surveys of peculiar velocities have produced highly appealing results (e.g. Courtois et al., 2013; Tully et al., 2014). However, since velocity fields are predominantly governed by largescale modes, which are accurately recovered by borg (borg sdss–inference), any improvement in the inference scheme is not expected to crucially change our results.
To characterize the accuracy at which our velocity fields are predicted, we analyzed their spectral properties against that of reference simulations produced with the same setup (presented in borg sdss–Tweb, ). From the Cartesian velocity components,^{4}^{4}4We do not use the divergence and curl directly computed using Lagrangian transport, but rather from the velocity field itself, for the reasons discussed in Hahn, Angulo & Abel (2015, section 5.1). we computed the divergence and the curl,
(46)  
(47) 
using FFTs on the Eulerian grid of side length 750 Mpc/ and voxels. We computed their respective power spectra, given by
(48)  
(49)  
To test the proportionality between and the density contrast , valid at first order in Lagrangian perturbation theory, we also computed the cross and autopower spectra
(50)  
(51) 
In this case, is estimated by standard CiC deposit of dark matter particles to the grid, and we apply a correction factor for aliasing effects and suppression of smallscale power close to the Nyquist wavenumber (Jing, 2005).^{5}^{5}5We here rely on the CiC scheme and not on the phasespace density fields estimators, following the arguments of Hahn, Angulo & Abel (2015, section 5.2).
In figure 10, we show , , and as a function of the wavenumber . The measurements from our simulations agree with the results of Hahn, Angulo & Abel (2015). The plots confirm the agreement between unconstrained and constrained density fields () at all scales (already discussed in borg sdss–Tweb).^{6}^{6}6We note that a lack of small scale power with respect to theoretical predictions, for , is observed in both unconstrained and constrained realizations. This is a gridding artifact due to the finite mesh size used to compute power spectra. This scale corresponds to around one quarter of the Nyquist wavenumber. Remarkable agreement with simulations can also be observed for and correlations. Deviations can only be observed for . Below this scale, structure formation becomes notably nonlinear, meaning that the 2LPT model used in the borg inference process breaks down, which explains the lack of smallscale power when the constrained initial conditions are run forward with a fully nonlinear code. The lack of and correlations does not exceed 30% for . Since the velocity field is curlfree in 2LPT, the vorticity is a unique prediction of the physical model accounting for full gravity. Strikingly, the prediction for still agrees well with the result of unconstrained simulations, meaning that cola filtering reinstates most of the final vorticity, and does so consistently with the inferred phases (see the panel in figure 3). Discrepancies are limited to an almost scaleindependent negative offset of 35%, which can be explained by known limitations of the borg data model: use of 2LPT instead of fully nonlinear gravity, only partial fitting of galaxy bias, absence of redshiftspace distortions and lightcone effects.
Summing up this section, at the demonstrated degree of accuracy, our results go beyond the state of the art for the characterization of cosmic velocity fields in the SDSS volume, including its low galaxydensity regions. They include the first quantitative prediction of the vorticity field, and come with the possibility of uncertainty quantification at the level of fields and summary statistics.
iv.4 Cosmography
For illustration purposes, as well as to demonstrate the power of inferring the phasespace structure of dark matter, we discuss in this section some cosmographic results concerning wellknown objects of the Local Universe. We have already discussed the Sloan Great Wall (Gott et al., 2005; Einasto et al., 2011), at celestial equatorial coordinates , , Mpc/, visible in figures 3, 6 and 9.
In figure 11, we show slices through the supergalactic plane (SGP). The supergalactic coordinate system is defined by the position of its North Pole () in the direction of Galactic coordinates (; , Lahav et al., 2000). The portion of the SGP () which is constrained in our analysis corresponds to , the rest is not covered by the SDSS main galaxy sample. The different panels show the posterior mean for: the density estimated from mass tracers (see section IV.1.2), the secondary stream density (see section IV.1.1), and the radial velocity field (see section IV.3). Galaxies in a Mpc/thick slice are overplotted. Some major structures of the Local Universe are clearly visible: the Coma cluster (, Mpc/), the UrsaMajoris supercluster (, Mpc/), and two farther clusters in the direction of UrsaMajor (, Mpc/) and UrsaMajor/ComaBerenices (, Mpc/) (e.g. Abell, Corwin & Olowin, 1989). One can also notice part of Leo supercluster (, Mpc/), which is mostly below the SGP (at ). Three voids intersect the SGP: the wellknown Boötes void (Kirshner et al., 1981), whose center is above the SGP (, , Mpc/) and which is around Mpc/ in radius (Batuski & Burns, 1985), a void in the direction of UrsaMajor (, Mpc/), and one in the direction of ComaBerenices (, Mpc/) (e.g. Fairall, 1998). Part of the VirgoBoötesHercules (VBH) filament, which extends from to at Mpc/ (Lavaux & Jasche, 2016) is also visible. For comparison, see also figure 4 in Carrick et al. (2015) and figure 7 in Lavaux & Jasche (2016).
As expected, structures observed in galaxies correlate well with the reconstructed dark matter density. At this mass resolution, up to 5 matter streams converge in the inner regions of Coma, the UrsaMajoris supercluster, and the VBH filament. The voids in ComaBerenices and in UrsaMajor are singlestream regions, as well as most voids of the Local Universe. However, our analysis reveals some substructure of the Boötes void. Reconstructed velocity fields gives further dynamical insight into these structures. Superclusters are still accreting matter, as seen in radial velocity “dipoles”. The ongoing emptying of voids is also clearly observable; for example, dark matter inside the Boötes void is streaming out towards its surrounding overdensities, the Coma cluster and the VBH filament.
V Summary and conclusions
This work discusses a fully probabilistic, chronocosmographic analysis of the Lagrangian dark matter sheet in a volume covered by the northern cap of the SDSS main galaxy sample. It is a data application of advanced tools, the application of which was so far limited to postprocessing simulations and deemed unfeasible in observations. The result is a wealth of new information about the dark matter underlying the observed galaxy distribution. Importantly, all these results are probabilistic, i.e. come with a quantification of uncertainty at the level of maps.
Specifically, we characterized the number of streams and the flow properties in Lagrangian and Eulerian space; we explained how to detect caustics; we improved estimates of the density field and showed that this can be thought of as an extended nonlinear filtering; and we obtained extremely accurate reconstructions of the velocity field and its derivatives (velocity divergence and vorticity). We introduced the Lagrangian classifier lich, a generalization of diva to the case of heterogeneous flows. lich is a more precise assignment of webtypes — distinguishing potential and vortical structures — that we have used to check the accuracy of the potential flow approximation. Using diva, lich, and origami, we presented the first maps of nearby Lagrangian cosmic web elements, and showed how to translate the result into redshiftzero Eulerian maps. With these definitions, identified structures have a direct physical interpretation in terms of the entire history of matter flows toward and inside them.
All the results shown in this paper are predictions based on fusing data constraints imprinted in the inference products (the initial conditions) and physical information expressed in the numerical structure formation model. At various places, we have checked the accuracy of these predictions. All posterior results are consistent with the respective priors, and deviate only in expected ways where the data model used for inference of the initial conditions has known shortcomings. This analysis therefore supports the standard paradigm of largescale structure formation. If our tests suggest that data models are not yet accurate enough to allow inference of cosmological model parameters, they reinforce the concept underlying borg, i.e. the general picture of including physical structure formation within survey data analysis. In particular, our first prediction of vorticity matches expectations at the % level, which is remarkable given that the inference model only relies on the assumption of a potential flow. A measure of the information content of Lagrangian cosmic web maps shows that it is vastly enhanced with respect to Eulerian ones, due to information transport along Lagrangian trajectories. The natural subsequent question of choosing the optimal cosmic web classifier for certain analyses is not addressed in this work; the required informationtheoretic framework is introduced in Leclercq et al. (2016).
This work constitutes the first threedimensional mapping of the nearby Lagrangian dark matter sheet, as constrained by the SDSS main galaxy sample. The nature of these maps is unprecedented, as they are dataconstrained with a high degree of control on uncertainties. The understanding of the phasespace properties of dark matter is a pathway to treating it in a more physical fashion at cosmological scales, beyond the phenomenological description of the CDM model. In particular, density fields and locations of caustics can be used for dark matter indirect searches in ray data from the Fermi satellite (see Gao et al., 2012); velocity fields can be used to account for peculiar velocity effects in supernova cosmology (see Davis et al., 2011) and to characterize the kinetic SunyaevZel’dovich effect (see Shao et al., 2011); tidal and vorticity fields can be used to gain better control on intrinsic alignments in the context of weak lensing (Codis et al., 2015) and CMB lensing (Larsen & Challinor, 2016). In addition to these specific examples, we believe that vast scientific opportunities to use our maps may emerge from the community. These data products, complementing to our earlier release of Eulerian density fields (borg sdss–inference) and Tweb maps (borg sdss–Tweb), are provided as online material along with this paper. They are available from the first author’s website, currently hosted at http://icg.port.ac.uk/leclercq/data.php.
Beyond the specific region probed in this work, our study demonstrates the possibility of mapping the dark matter sheet in the real Universe. This possibility, so far limited to simulations, is uniquely offered by Bayesian physical inference of the initial conditions from largescale structure surveys, which allows a selfconsistent treatment of the necessary aspects: handling of observational effects and description of structure formation. Thus, our result may have important implications for largescale structure data analyses, as they allow new confrontations of observational data and theoretical models, and therefore novel tests of the standard paradigm of cosmic web formation and evolution.
Note added: After submission of this paper, several related pieces of work appeared in the literature. These include a new cosmic web classifier using the eigenvalues of the Hessian of the multistream field (Ramachandra & Shandarin, 2017), a new derivation of the theory of caustics (Feldbrugge et al., 2017) and a comparison project of different cosmic web algorithms (Libeskind et al., 2017).
Statement of contribution
FL implemented cola and the required cosmic web analysis tools, performed the study, produced the maps and wrote the paper. JJ developed borg and lead the SDSS analysis; GL developed diva; JJ and GL contributed to the development of computational tools. BW was involved in the conception and design of Bayesian largescale structure inference. BW and WP contributed to the interpretation of results. All the authors read and approved the final manuscript.
Acknowledgements.
FL would like to thank Oliver Hahn, Mark Neyrinck and Rien van de Weygaert for discussions that stimulated this analysis. The manuscript benefited from constructive comments and useful suggestions from the reviewer. We are grateful to Bridget Falck and Mark Neyrinck for making publicly available their origami code, which was used in this work. Numerical computations were done on the Sciama High Performance Compute (HPC) cluster which is supported by the ICG, SEPNet and the University of Portsmouth. This work has also made use of the Horizon Cluster hosted by the Institut d’Astrophysique de Paris; we thank Stéphane Rouberol for running smoothly this cluster for us. FL and WP acknowledge funding from the European Research Council through grant 614030, Darksurvey. JJ is partially supported by a Feodor Lynen Fellowship by the Alexander von Humboldt foundation. BW acknowledges funding from an ANR Chaire d’Excellence (ANR10CEXC00401) and the UPMC Chaire Internationale in Theoretical Cosmology. This work has been done within the Labex Institut Lagrange de Paris (reference ANR10LABX63) part of the Idex SUPER, and received financial state aid managed by the Agence Nationale de la Recherche, as part of the programme Investissements d’avenir under the reference ANR11IDEX000402. This research was supported by the DFG cluster of excellence “Origin and Structure of the Universe”.References
References
 Abel, Hahn & Kaehler (2012) (Abel, Hahn & Kaehler, 2012) T. Abel, O. Hahn, R. Kaehler, Tracing the dark matter sheet in phase space, Mon. Not. R. Astron. Soc. 427, 61 (2012), arXiv:1111.3944 [astroph.CO].
 Abell, Corwin & Olowin (1989) (Abell, Corwin & Olowin, 1989) G. O. Abell, H. G. Corwin, Jr., R. P. Olowin, A catalog of rich clusters of galaxies, Astrophys. J. Supp. 70, 1 (1989).
 AragonCalvo & Szalay (2013) (AragonCalvo & Szalay, 2013) M. A. AragonCalvo, A. S. Szalay, The hierarchical structure and dynamics of voids, Mon. Not. R. Astron. Soc. 428, 3409 (2013), arXiv:1203.0248 [astroph.CO].
 AragonCalvo & Yang (2014) (AragonCalvo & Yang, 2014) M. A. AragonCalvo, L. F. Yang, The hierarchical nature of the spin alignment of dark matter haloes in filaments, Mon. Not. R. Astron. Soc. 440, L46 (2014), arXiv:1303.1590 [astroph.CO].
 AragónCalvo et al. (2007) (AragónCalvo et al., 2007) M. A. AragónCalvo, B. J. T. Jones, R. van de Weygaert, J. M. van der Hulst, The multiscale morphology filter: identifying and extracting spatial patterns in the galaxy distribution, Astron. & Astrophys. 474, 315 (2007), arXiv:0705.2072.
 AragónCalvo et al. (2010) (AragónCalvo et al., 2010) M. A. AragónCalvo, E. Platen, R. van de Weygaert, A. S. Szalay, The Spine of the Cosmic Web, Astrophys. J. 723, 364 (2010), arXiv:0809.5104.
 Arnold, Shandarin & Zeldovich (1982) (Arnold, Shandarin & Zeldovich, 1982) V. I. Arnold, S. F. Shandarin, I. B. Zeldovich, The large scale structure of the universe. I  General properties One and twodimensional models, Geophysical and Astrophysical Fluid Dynamics 20, 111 (1982).
 Bardeen et al. (1986) (Bardeen et al., 1986) J. M. Bardeen, J. R. Bond, N. Kaiser, A. S. Szalay, The statistics of peaks of Gaussian random fields, Astrophys. J. 304, 15 (1986).
 Batuski & Burns (1985) (Batuski & Burns, 1985) D. J. Batuski, J. O. Burns, Finding lists of candidate superclusters and voids of Abell clusters, Astronom. J. 90, 1413 (1985).
 Bernardeau et al. (2002) (Bernardeau et al., 2002) F. Bernardeau, S. Colombi, E. Gaztañaga, R. Scoccimarro, Largescale structure of the Universe and cosmological perturbation theory, Physics Reports 367, 1 (2002), astroph/0112551.
 Blanton et al. (2005) (Blanton et al., 2005) M. R. Blanton, D. Eisenstein, D. W. Hogg, D. J. Schlegel, J. Brinkmann, Relationship between Environment and the Broadband Optical Properties of Galaxies in the Sloan Digital Sky Survey, Astrophys. J. 629, 143 (2005), astroph/0310453.
 Bond, Kofman & Pogosyan (1996) (Bond, Kofman & Pogosyan, 1996) J. R. Bond, L. Kofman, D. Pogosyan, How filaments of galaxies are woven into the cosmic web, Nature 380, 603 (1996), astroph/9512141.
 Burden, Percival & Howlett (2015) (Burden, Percival & Howlett, 2015) A. Burden, W. J. Percival, C. Howlett, Reconstruction in Fourier space, Mon. Not. R. Astron. Soc. 453, 456 (2015), arXiv:1504.02591.
 Carrick et al. (2015) (Carrick et al., 2015) J. Carrick, S. J. Turnbull, G. Lavaux, M. J. Hudson, Cosmological parameters from the comparison of peculiar velocities with predictions from the 2M++ density field, Mon. Not. R. Astron. Soc. 450, 317 (2015), arXiv:1504.04627.
 Cautun, van de Weygaert & Jones (2013) (Cautun, van de Weygaert & Jones, 2013) M. Cautun, R. van de Weygaert, B. J. T. Jones, NEXUS: tracing the cosmic web connection, Mon. Not. R. Astron. Soc. 429, 1286 (2013), arXiv:1209.2043 [astroph.CO].
 Cautun et al. (2014) (Cautun et al., 2014) M. Cautun, R. van de Weygaert, B. J. T. Jones, C. S. Frenk, Evolution of the cosmic web, Mon. Not. R. Astron. Soc. 441, 2923 (2014), arXiv:1401.7866.
 Chan (2014) (Chan, 2014) K. C. Chan, Helmholtz decomposition of the Lagrangian displacement, Phys. Rev. D 89, 083515 (2014), arXiv:1309.2243.
 Chong, Perry & Cantwell (1990) (Chong, Perry & Cantwell, 1990) M. S. Chong, A. E. Perry, B. J. Cantwell, A general classification of threedimensional flow fields, Physics of Fluids 2, 765 (1990).
 Codis et al. (2015) (Codis et al., 2015) S. Codis, R. Gavazzi, Y. Dubois, C. Pichon, K. Benabed, V. Desjacques, D. Pogosyan, J. Devriendt, A. Slyz, Intrinsic alignment of simulated galaxies in the cosmic web: implications for weak lensing surveys, Mon. Not. R. Astron. Soc. 448, 3391 (2015), arXiv:1406.4668.
 Courtois et al. (2013) (Courtois et al., 2013) H. M. Courtois, D. Pomarède, R. B. Tully, Y. Hoffman, D. Courtois, Cosmography of the Local Universe, Astronom. J. 146, 69 (2013)