Theory

The Delaunay Tessellation Field Estimator (DTFE) and its extension to phase space (PS-DTFE) are mathematical tools for the reconstruction of the density and velocity field of a discrete point set. We review here the derivations underlying these reconstruction methods.

The Delaunay Tessellation Field Estimator

Delaunay Tessellation Field Estimator (DTFE) was developed by Willem Schaap and Rien van de Weijgaert (Schaap and Weijgaert 2004, Willem Schaap 2007). The DTFE method dynamically adapts to the variation of the density and geometry of the point set. The method is used in cosmological datasets (see for example the astronomy picture of the day, 7 November 2007) as it can simultaneously capture the geometry of the voids, walls, filaments and clusters and preserve the total mass of the pointset. The total mass is not preserved in the closely related natural neighbor interpolation based on the Voronoi tessellation. For an efficient implementation of the DTFE method in C++ see the code by Marius Cautun.

Consider the pointset $\mathcal{P}$ consisting of $N$ labelled points $p_i \in \mathbb{R}^d$, the velocity set $\mathcal{V}$ consisting of the $N$ velocities $v_i \in \mathbb{R}^d$, and the mass set $\mathcal{M}$ consisting of the masses $m_i \in \mathbb{R}$ corresponding to the points in $\mathcal{P}$, with $i=1,\dots,N$. Using the points and masses, we reconstruct the density field $\rho:\mathbb{R}^d \to \mathbb{R}$. Using both the points and the velocities, we construct the velocity field $v:\mathbb{R}^d\to \mathbb{R}$.

Density reconstruction

Given the pointset $\mathcal{P}$, consider the Delaunay tesselation $\mathcal{D}$ consisting of $N_T$ labelled simplices $D_i$, with $i=1,\dots, N_T$. In two dimensions, a simplex is a triangle spanning three points in $\mathcal{P}$. In three dimensions, a simplex is a tetrahedron spanning four points in $\mathcal{P}$.

Let's assume we can associate a (to be determined) density estimate $\rho_i$ to each point in $\mathcal{P}$. Given a simplex $D \in \mathcal{D}$ spanned by the vertices $p_{l_0},\dots, p_{l_d} \in \mathcal{P}$, with the corresponding densities $\rho_{l_0}, \dots, \rho_{l_d}$, we construct a linear interpolation of the density field in the simplex as

\[\rho(x) = \rho_{l_0} + [\nabla \rho] (x-p_{l_0})\,,\]

with $x \in D$ and the gradient vector $\nabla \rho\in \mathbb{R}^d$ associated to the simplex determined by the $d$ linear consistency relations $\rho(p_{l_i}) = \rho_{l_i}$ for $i=1,\dots,d$. In matrix notation,

\[\nabla \rho = \begin{pmatrix} p_{l_1}-p_{l_0}\\ \vdots\\ p_{l_d}-p_{l_0}\\ \end{pmatrix}^{-1} \begin{pmatrix} \rho_1-\rho_0\\ \vdots\\ \rho_d-\rho_0\\ \end{pmatrix}\,.\]

This system of equations is solvable when none of the points $p_{l_0}, \dots, p_{l_d}$ are collinear. The integral over the linear interpolation yields

\[\int_D \rho(x)\mathrm{d}x = \frac{V(D)}{d+1} \sum_{i \in D} \rho_i \,,\]

with $V(D)$ the volume of the simplex $D$ and $i$ the labels of the vertices of $D$ in $\mathcal{P}$. Note that the volume of a simplex $D$ spanned by $p_{l_0}, \dots, p_{l_d}$ can be expressed as the determinant

\[V(D) = \frac{1}{d!} \begin{vmatrix} p_{l_1}-p_{l_0}\\ \vdots\\ p_{l_d}-p_{l_0}\\ \end{vmatrix}\,.\]

The integral over density – assuming the density vanishes outside the convex hull of the Delaunay tesselation – takes the form

\[\begin{align} \int \rho(x) \mathrm{d}x &= \sum_{i=1}^{N_T} \int_{D_i} \rho(x)\mathrm{d}x\\ &= \frac{1}{d+1} \sum_{i=1}^{N_T} V(D_i) \sum_{j \in D_i} \rho_j\,, \end{align}\]

where the first sum runs over the simplices of the tessellation and the second sum runs over the vertices of a given simplex. Note that $\rho_i$ enters the sum for each simplex for which it is a vertex. These simplices form the star $W_i$ of the point $p_i$

Using this observation, we reorder the double sum, by collecting the terms involving $\rho_i$ leading to the terms $\rho_i(V(D_{l_0}) + \dots + V(D_{l_n})) = \rho_i V(W_i)$, with the $D_{l_i}$'s forming the star of $p_i$. The integral over the density is now a single sum over the points in $\mathcal{P}$, i.e.,

\[\int \rho(x) \mathrm{d}x =\frac{1}{d+1} \sum_{i=1}^N \rho_i V(W_i)\,.\]

The key observation in DTFE is that when we chose the natural estimate of the density,

\[\rho_i = \frac{(d+1) m_i}{V(W_i)},\]

by which the density at $p_i$ only depends on the mass at the point and the local geometry, the integral over the density reduces to the total mass of the point set

\[\int \rho(x)\mathrm{d}x = \sum_{i=1}^N m_i\]

as one might expect of a density reconstruction method.

Now, given a point $x \in \mathbb{R}^d$, we can reconstruct the density field at this point by finding the simplex in which the point lies and evaluating the linear interpolation of the density in the simplex.

Velocity reconstruction

The velocity field of the dataset is analogously reconstructed with a linear interpolation in the Delaunay cells. Given a simplex $D$ spanned by the vertices $p_{l_0}, \dots, p_{l_d}$ and the associated velocities $v_{l_0},\dots,v_{l_d}$, we write the velocity field in the simplex $D$ as

\[v(x) = v_{l_0} + \nabla v (x - p_{l_0})\]

with the velocity gradient $\nabla v$ associated to the simplex determined by the linear consistency relations $v(p_{l_i}) = v_{l_i}$ for $i=1,\dots,d$. In matrix notation,

\[\nabla v = \begin{pmatrix} p_{l_1}-p_{l_0}\\ \vdots\\ p_{l_d}-p_{l_0}\\ \end{pmatrix}^{-1} \begin{pmatrix} v_1-v_0\\ \vdots\\ v_d-v_0\\ \end{pmatrix}\,.\]

Now, given a point $x \in \mathbb{R}^d$, we can reconstruct the velocity field at this point by finding the simplex in which the point lies and evaluating the corresponding linear interpolation in the simplex.

Note that the gradient $\nabla v$ is a piecewise constant function. Given the gradient $\nabla v$ for each simplex, it is natural to evaluate velocity deformation modes. In three dimensions, we evaluate the divergence $\theta$, the shear $\sigma_{ij}$ and the vorticity $\omega = \epsilon^{ijk} \omega_{ij}$ defined by

\[\begin{align} \theta &= \nabla \cdot v = \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} + \frac{\partial v_z}{\partial z},\\ \sigma_{ij} &= \frac{1}{2} \left[\frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i}\right] - \frac{1}{3} (\nabla \cdot v) \delta_{ij},\\ \omega_{ij} &=\frac{1}{2} \left[\frac{\partial v_i}{\partial x_j} - \frac{\partial v_j}{\partial x_i}\right] \end{align}\]

with the Kronecker delta $\delta_{ij}$ and the Levi-Civita symbol $\epsilon^{ijk}$.

For a point $x \in \mathbb{R}^d$, the reconstruction amounts to looking up the corresponding simplex and velocity deformation modes.

The Phase-Space Delaunay Tessellation Field Estimator

The Phase-Space Delaunay Tessellation Field Estimator (PS-DTFE) extends the DTFE method to phase-space, inspired by work by [...]. For more details we refer to the publication Phase-Space Delaunay Tesselation Field Estimator.

The evolution of our universe can be described in Lagrangian fluid dynamics in terms of the Lagrangian map $x_t(q) = q + s_t(q)$ mapping a point in the space of initial conditions (Lagrangian space) to a point in the current universe (Eulerian space). The displacement field $s_t(q)$ captures the displacement of a particle starting at $q$ in time $t$. Large-scale structure formation can be visualized with the folding of a three-dimensional dark matter sheet in six-dimensional phase-space. The projection of the dark matter sheet into Eulerian space leads to the formation of multistream regions bounded by caustics in which the density spikes.

Given the Lagrangian map $x_t$, the density in a point $x'$ in Eulerian space is given by

\[\rho(x') = \sum_{q \in x_t^{-1}(x')} \frac{\rho_u}{|\det( \nabla x_t(q))|} = \sum_{q \in x_t^{-1}(x')} \frac{\rho_u}{|\det( I + \nabla s_t(q))|} \,,\]

where the sum ranges over all points in Lagrangian space that stream into $x'$ at the given time $t$. The caustics form at the boundaries of the multistream regions where the determinant of the deformation tensor vanishes, i.e. $\det(\nabla x_t)=0$.

The DTFE method successfully estimates the density in the single-stream regions. However, in multistream regions, the Delaunay tessellation may identify particles as neighbours in Eulerian space that where far separated in Lagrangian space (see the central panel of the figure above). To circumvent this problem, we evaluate the Delaunay tessellation of the particle coordinates in the early (high-redshift) configuration (that is, in Lagrangian space), where the universe was still in a single-stream state. We then evolve this tessellation with the displacement field and use the evolved tessellation to estimate the density and velocity fields in Eulerian space (see the right panel of the figure above). To this end, we identify the simplices (tetrahedra in three dimensions) in which a point $x'$ resides. This is efficiently implemented using a Bounding Volume Hierarchy (BVH), a tree structure on a set of tetrahedra spanning the dark matter sheet. Having identified the incoming simplices, the reconstruction of the density and velocities is analogous to the DTFE method, with the estimates in the individual streams summed to according to the discussion above.