The Astrophysical Journal, 858:87 (15pp), 2018 May 10 https://doi.org/10.3847/1538-4357/aabb53
© 2018. The American Astronomical Society. All rights reserved.
Minkowski Tensors in Two Dimensions: Probing the Morphology and Isotropy of the
Matter and Galaxy Density Fields
Stephen Appleby1, Pravabati Chingangbam2, Changbom Park1 , Sungwook E. Hong3 , Juhan Kim4 , and Vidhya Ganesan2,5
1 School of Physics, Korea Institute for Advanced Study, 85 Hoegiro, Dongdaemun-gu, Seoul 02455, Republic of Korea; stephen@kias.re.kr
2 Indian Institute of Astrophysics, Koramangala II Block, Bangalore 560 034, India
3 Korea Astronomy and Space Science Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea
4 Center for Advanced Computation, Korea Institute for Advanced Study, 85 Hoegiro, Dongdaemun-gu, Seoul 02455, Republic of Korea
5 Indian Institute of Science, Bangalore 560 034, India
Received 2017 September 15; revised 2018 March 5; accepted 2018 March 31; published 2018 May 10
Abstract
We apply the Minkowski tensor statistics to two-dimensional slices of the three-dimensional matter density field.
The Minkowski tensors are a set of functions that are sensitive to directionally dependent signals in the data and,
furthermore, can be used to quantify the mean shape of density fields. We begin by reviewing the definition of
Minkowski tensors and introducing a method of calculating them from a discretely sampled field. Focusing on the
statistic W1,12 —a 2×2 matrix—we calculate its value for both the entire excursion set and individual connected
regions and holes within the set. To study the morphology of structures within the excursion set, we calculate the
eigenvalues λ1, λ2 for the matrix W1,12 of each distinct connected region and hole and measure their mean shape
using the ratio b º ál l ñ. We compare both W1,12 1 2 and β for a Gaussian field and a smoothed density field
generated from the latest Horizon Run 4 cosmological simulation to study the effect of gravitational collapse on
these functions. The global statistic W1,12 is essentially independent of gravitational collapse, as the process
maintains statistical isotropy. However, β is modified significantly, with overdensities becoming relatively more
circular compared to underdensities at low redshifts. When applying the statistics to a redshift-space distorted
density field, the matrix W1,12 is no longer proportional to the identity matrix, and measurements of its diagonal
elements can be used to probe the large-scale velocity field.
Key words: dark matter – galaxies: evolution – methods: statistical
1. Introduction N dimensional field (Doroshkevich 1970; Adler 1981; Gott
One of the fundamental axioms implicit within the standard et al. 1986, 1989; Hamilton et al. 1986; Melott et al. 1989;
cosmological model is that the distribution of matter in the Ryden et al. 1989; Park & Gott 1991; Park et al. 1992;
universe is statistically isotropic and homogeneous when Matsubara 1994, 1996; Schmalzing et al. 1996; Kerscher et al.
smoothed over suitably large scales. This condition is very 2001; Park et al. 2005). Since they are scalar quantities, they
well observed in the early epoch of radiation and matter cannot inform us of any directionally dependent information
domination, where fluctuations in the dark matter density eld contained within the data. The concept of Minkowskifi
are small. However, the scale at which this remains true at low functionals can be generalized to vector and tensor counterparts
redshifts is less clear, as nonlinear gravitational evolution (McMullen 1997; Alesker 1999; Beisbart et al. 2002; Hug et al.
generates a complex web of structures. We expect alignment of 2008; Schroder-Turk et al. 2010, 2013); these quantities are
structures due to their position within the cosmic web (Lee & typically defined as integrals of some higher-rank quantity over
Pen 2002; Aubert et al. 2004; Patiri et al. 2006; Hahn et al. the boundary of an excursion set. As such, they contain
2007; Lee et al. 2008; Paz et al. 2008; Codis et al. 2014, 2015a, information not present in the standard Minkowski functionals.
2015b) and a bias in the clustering properties of galaxies. The In particular, they can be used to identify globally anisotropic
dark matter field exhibits coherent structures even at very large signals in the data, as well as provide a measure of the shape of
scales (~100 h-1 Mpc), so the scale at which alignments cease the peaks/troughs of a density field when applied to individual
to become significant remains an open question. connected regions and holes in an excursion set. In both
Furthermore, when introducing an observer, one can state instances, the Minkowski tensors measure directions associated
that the observed distribution of dark matter tracers (typically with a boundary.
galaxies) are neither homogeneous nor isotropic; selection The application of Minkowski tensors to cosmology is a
effects generate a nontrivial radial profile in the observed relatively new phenomenon. Ganesan & Chingangbam (2017)
number density, and redshift-space distortion effects will applied a Minkowski tensor that encodes shape and alignment
modify the apparent positions of galaxies along the line of information of structures to the two-dimensional cosmic
sight. Line-of-sight effects will generate a bias in the detection microwave background (CMB) fields. The authors showed
of structures perpendicular to the line of sight. Isotropy of the that the 2015 E-mode Planck data (Adam et al. 2016) exhibit a
galaxy sample is lost via masks and boundaries. If we can higher than 3σlevel of anisotropy or alignment of hot spots
measure the degree of anisotropy in data sets, then we can and cold spots. Analytic expressions of translation invariant
potentially minimize observational systematics and, in the case Minkowski tensors for Gaussian random fields in two
of redshift-space distortion, constrain the growth rate. dimensions have been derived in Chingangbam et al. (2017).
The N+1 Minkowski functionals are a set of scalar In this work, we apply the Minkowski tensors to two-
quantities that characterize the morphology and topology of an dimensional slices of the dark matter density field. We first
1
The Astrophysical Journal, 858:87 (15pp), 2018 May 10 Appleby et al.
Figure 1. Marching-squares algorithm. Each set of four adjacent pixels for the discretized field δij can take one of 16 possible combinations of “in” or “out” states.
Black circles denote points in our grid for which the density lies inside the excursion set δij>σ0ν, and white circles are “out” points with δij<σ0ν. We label each
case with an integer 1„Nc„16. Between each “in” and “out” state, we linearly interpolate between corners of the box to find the point along the edge of the square
that satisfies δ=σ0ν. We then connect these vertices, shown as solid black arrows. The arrow defines the boundary of the excursion region and is directed such that it
always flows anticlockwise around the “in” states. The cases Nc=7 and Nc=10 are ambiguous, as discussed in the text.
review the generalization of the Minkowski functionals to their Our starting point is a discrete two-dimensional density field
tensor equivalents. We then ask how these quantities can be δij on a regular lattice spanned by i, j subscripts, 1  i  Npix,
used to test the isotropy of the field. Throughout this paper, we 1„j„Npix, where Npix is the number of grid points in one
focus on two-dimensional slices of a three-dimensional dimension. The domain is chosen to be a square with periodic
volume; in a companion paper, we consider the three- boundary conditions, but this condition is not necessary. We
dimensional generalization of these statistics. define a dimensionless density threshold n = dc s0, where δc is
In the following sections, we provide a thorough explanation a constant density threshold and σ0 is the rms fluctuation of δij.
of our construction of the Minkowski functionals and tensors A perimeter of constant density δc=σ0ν defines an excursion
by generating the boundary of an excursion set in two set of the field. We can label each (i, j) pixel as either “in” or
dimensions. We then define the Minkowski tensors and show “out” of the excursion set according to δij>νσ0 or δij<νσ0,
how they can be calculated for a discrete field and bounding respectively. Our intention is to construct a closed boundary
perimeter. We apply our algorithms to a Gaussian random field, perimeter that separates in/out pixels.
connecting our numerical results to analytic predictions We adopt the method of marching squares (Mantz
wherever possible. We close by applying the statistics to the et al. 2008). The method performs a single sweep through
late-time gravitationally evolved dark matter field using the the entire grid systematically from one corner to the opposite.
latest Horizon Run simulation. At each grid point (i, j), we generate a square from its adjacent
pixels—they are (i, j), (i + 1, j), (i, j + 1), (i + 1, j + 1).
2. Generating the Boundary of an Excursion Set: Two Each of these pixels can be either “in” or “out” of the
Dimensions excursion set, so there are 24=16 possible unique states of the
We begin with a discussion of our construction of a square. In Figure 1, we exhibit the standard 16 states, where we
bounding perimeter enclosing a subset of a density field in two use the integer 1„Nc„16 to define each case as labeled.
dimensions. Our analysis in this section will closely follow that Each point denotes a δij vertex; black are “in” states δij>σ0ν,
of Schroder-Turk et al. (2010), but we detail the method for and white are “out” states δij<σ0ν. A closed bounding
completeness. perimeter is then constructed based on the 16 cases by linearly
2
The Astrophysical Journal, 858:87 (15pp), 2018 May 10 Appleby et al.
Figure 3. We exhibit two squares, (i, j) and (i, j + 1), which represent cases
Figure 2. Ambiguous cases 7 and 10. The “in” states can either be connected Nc=2 and Nc=6 in Figure 1. By interpolation, we generate the black arrows,
(7 and 10) or disconnected (7b and 10b). Using either the original case or its which are the vectors eij and eij+1. The brown arrows denote the unit normals to
corresponding “b” state will yield a closed curve. To determine which to use, eij and eij+1; the angle between them is bi. We note that βi is directed
we calculate the density at the center of the square. If it is “in,” then we assume following the anticlockwise conventions of our method.
that the “in” states are linked; otherwise, we assume that the “out” states are
linked (the “b” states).
we define the three Minkowski functionals W0,1,2 as
1 1
W0 =
A ò da = å ∣DAn ∣, (1)interpolating along the edges of the squares. Specifically, one Q A n
can note that whenever an “in” and “out” state are joined along 1 1
an edge of a square, we linearly interpolate between the values W1 = ò dℓ = å ∣e∣, (2)4A ¶Q 4A e
of δ at these vertices along the edge to the point at which
δ=σ0ν is reached. This defines a vertex in the bounding 1 1W2 = kdℓ = å bi, (3)
perimeter. Vertices are then joined according to Figure 1 this 2pA ò— ¶Q 2pA i
defines the edge components of the boundary, which
where ... da and ... dℓ are integrals over the area and
correspond to the solid arrows in the figure. Finally, we use òQ ò¶Q
trigonometry to calculate the area enclosed by the bounding perimeter of an excursion set, respectively, and κ is the local
perimeter in each square (the shaded region in each case in curvature. Here ∣DAn∣ is the area of the shaded region in each
Figure 1). The perimeter of the boundary, exhibited as black pixel square, ån is the sum over all pixel squares, and ∣e∣ is the
arrows in the figure, is directed such that the arrow always length of the boundary in each square (the length of the solid
flows anticlockwise around the “in” vertices δij>σ0ν. black arrows in Figure 1). The sum åe indicates the sum over
There is a caveat to the method: there is an ambiguity all edge segments in the discrete boundary, and A is the total
regarding the cases Nc=7 and Nc=10. In these cases, each area of the two-dimensional plane. The angle βi between
edge of the square will have a vertex belonging to the excursion normals of adjacent perimeter segments is exhibited in Figure 3
boundary, and we can join these vertices in two different ways. ; the genus is simply the sum of all such angles. The sum åi
In Figure 2, we exhibit the ambiguity, labeling the squares indicates the sum over all vertices in the bounding perimeter.
Nc=7, 7b and Nc = 10, 10b. Regardless of which of these The genus is a topological quantity that measures the number
cases we choose to adopt, the method will always yield a of connected regions minus the number of holes.
closed bounding perimeter. Furthermore, Nc=7, 10 are rare The above algorithm will allow us to calculate the Minkowski
configurations when calculating the bounding surface of fields functionals and their generalizations, the Minkowski tensors,
that are smoothed over more than a few pixel lengths. which will be defined in Section 3. These quantities describe the
Nevertheless, one must still account for the ambiguity. We global properties of the excursion set. However, the total
select either 7 or 7b by estimating the value of δ at the center of excursion set will be composed of a set of disconnected “in” and
the square simply as the mean of the four vertices. If this value “out” subregions; see Figure 1. To calculate the properties of the
is “in” (d > s0n), then we assume that the two “in” vertices of subregions, we apply a simple friends-of-friends algorithm to the
the square belong to the same excursion region (that is, we density field. For each δij that is inside the excursion set, we
select Nc=7). Otherwise, we select 7b. We perform a similar assign all points di1,j1 as belonging to the same subregion if
operation for case 10. they are also within the excursion set, and we repeat the
Once we have generated the vertices that define the bounding operation iteratively on these points. The only caveat is again the
perimeter of the excursion set, we can calculate its total length cases Nc=7 and Nc=10 in Figure 1—if two “in” vertices are
and enclosed area. These two quantities are proportional to the linked diagonally in the box they share (for example, δij and
Minkowski functionals. The final Minkowski functional in two di+1,j+1 are inside the excursion set and di,j+1 and di+1,j are out),
dimensions is the genus; this can be generated by first then we test whether the box is 7, 10 or 7b, 10b by calculating
calculating the normals to the bounding perimeter. Then, the the central value of the density in the box. If the square is 7, 10,
genus is linearly related to the sum of the angles between the then the diagonal “in” vertices are assumed to belong to the same
normals of the adjacent perimeter edge sections. Specifically, excursion region, otherwise they are not assigned as friends.
3
The Astrophysical Journal, 858:87 (15pp), 2018 May 10 Appleby et al.
Note that they may still ultimately be linked via our algorithm,
just not through a 7b or 10b box.
Once we have assigned each “in” grid point to a particular
excursion subregion (there are cid distinct subregions), we can
calculate the area, perimeter, and genus of each one,
constructing a distribution of Minkowski functionals for each
density threshold ν. Furthermore, we can perform the same
algorithm but instead tracking the “out” states—this will yield
the properties of the individual holes in the field.
The ability of this algorithm to accurately reproduce the
bounding perimeter of an excursion set decreases for
structures that are poorly resolved, specifically for peaks
that have sizes roughly equal to our pixel resolution. As
cosmological density fields exhibit structure on all scales, we
must be careful to check that numerical artifacts do not impact
our results. In Appendix A, we highlight two principal sources Figure 4. We exhibit a schematic diagram of the area Q of an excursion set, its
of numerical error and attempt to quantify the size of these boundary ∂Q as a green solid line, and the vectors r, nˆ and ê used in the
effects. We find that smoothing the field over more than five construction of the Minkowski tensors.
pixel lengths is sufficient to ensure that marching squares
reconstruct the excursion-set boundary of the dark matter field Further information may be extracted from higher-rank
to better than 1% accuracy for thresholds −4„ν„4. generalizations m + n > 2, but we do not consider these
quantities in this work.
3. Minkowski Tensors: Definition There exist relations between Equations (7) and (10) and
The Minkowski functionals are scalar quantities. In Wj , where  is the identity matrix and Wj, ( j=0, 1, 2) are
McMullen (1997), Alesker (1999), Beisbart et al. (2002), the scalar Minkowski functionals (McMullen 1997):
Hug et al. (2008), and Schroder-Turk et al. (2010, 2013), the W  = 2W1,1, (11)
vector and tensor generalizations were constructed; we direct 0 1
the reader to these works for the details of their definition. pW  = W 0,21 1 + W
1,1
2 , (12)
These statistics were applied to cosmology in Beisbart et al. 2
(2001, 2001), Mecke et al. (1994), Schmalzing & Buchert W  = 2W 0,2. (13)
(1997), and Ganesan & Chingangbam (2017). 2 2
The Minkowski tensors of rank (m, n) of a field in a flat two- These relations imply that W 0,2 1,12 and W1 carry no additional
dimensional space are given by information relative to the scalar Minkowski functionals. The
W1,12 andW
0,2
1 tensors carry new information, with the sum being
Wm,0
1
= ò rmda, (4) related to W1 according to Equation (12). It is sufficient to0 A Q measure one of these two tensors, with the other containing no
new information. Here W1,12 is related to W
0,2
1 via a rotation,
m,n 1W1 = ò rm Ä nˆndℓ, (5) W 0,21 = pTW1,1T t2 2, where T is the p 2 rotation matrix and T
t is
4A ¶Q its transpose.
1 The tensor W
1,1
2 can be reexpressed as (Chingangbam
Wm,n2 =
2pA ò r
m Ä nˆnkdℓ, (6) et al. 2017)
¶Q
1,1 1
where r is the two-dimensional position vector and n̂ is the unit W2 = ò eˆ Ä eˆ dℓ, (14)2pA ¶Q
normal to the tangent vector of the bounding perimeter. We
schematically present the vectors r, nˆ and ê—the unit tangent where ê is the unit tangent to the curve. For a discretized field,
vector to the boundary—in Figure 4. this formula can be expressed in component form as
The rank-zero Minkowski tensors are the standard Minkowski
1,1 1 1 -1
functionals. Our focus is on rank-two Minkowski tensors (W2 )ij = (rinj + rjni)kdℓ = ∣e∣ eiej,
m + n = 2, speci cally the subset that is translationally 4p
å
A ò¶Q 2pA
fi e
invariant: (15)
where the i, j subscripts run over the standard two-dimensional
1
W1,11 = ò r Ä nˆdℓ, (7) x1, x2 orthogonal coordinates. The sum is over all edge4A ¶Q
1 segments of the excursion-set perimeter, ei is the length of
W1,12 = ò r Ä nˆkdℓ, (8) the boundary segment in the ith direction, and ∣e∣ is the length2pA ¶Q of the two-dimensional vector e. The diagonal components of
1 (W1,10,2 2 )ij are proportional to the (squared) total length of theW1 = ò nˆ Ä nˆdℓ, (9)4A ¶Q excursion-set bounding perimeter in the ith direction, and the
off-diagonal component is the cross term. The existence of a
0,2 1W = ò nˆ Ä nˆkdℓ. (10) preferred direction in the excursion boundary will manifest as2 2pA ¶Q an inequality between the diagonal components of (W1,12 )ij.
4
The Astrophysical Journal, 858:87 (15pp), 2018 May 10 Appleby et al.
Figure 5. We exhibit a 500 ´ 500(h-1 Mpc)2 subset of a Gaussian random
field as a heat map. The field has a flat power spectrum with a Gaussian
smoothing kernel. Our algorithm generates a bounding perimeter of constant ν;
we exhibit two examples, ν=σ0 and ν=1.4σ0, as black/red contours.
The function (W1,12 )ij can be defined not only over the entire
excursion-set perimeter but also over each distinct subregion
(both connected region and hole). The principal axes of each
subregion will be aligned in different directions; so, to measure
the shapes of these structures, we extract the eigenvalues of
(W1,12 )ij. The result is a pair of (λ1, λ2) values for each
connected region and hole in the set.
We define the mean ratio of eigenvalues of all individual
excursion subregions as
l2 l lbc º b º
2
h btot º
2 , (16)
l1 c l1 h l1 tot
where áñc,h,tot denotes the sample average over all individual
connected regions, holes, and combined connected regions and
holes, respectively. Therefore, bc,h,tot  1 provides information
regarding the mean shape of the excursion regions. Here
bc,h,tot = 1 corresponds to a perfectly isotropic average shape,
and any value less than unity indicates some level of
anisotropy, either ellipticity or a more general departure from Figure 6. Minkowski functionals of a two-dimensional Gaussian field. The
blue points and error bars are the mean and error on the mean of N =100
isotropy. realrealizations of a Gaussian field with a flat power spectrum. The solid black line
Additional information is contained within W1,12 relative to is the analytic prediction.
the scalar Minkowski functionals. The statistic is invariant
under translations, and a perfectly isotropic field would
correspond to a diagonal matrix with equal components. Any We apply our two-dimensional marching-squares algorithm
departure from this equivalence will signify a preferred to the resulting δij. In Figure 5, we exhibit a small
direction in the bounding perimeter of the excursion set. 500 ´ 500(h
-1 Mpc)2 subset of the density field. We also
exhibit an example of our algorithm: we apply a density
4. Applications: Two-dimensional Gaussian Random Field threshold ν=σ0, ν=1.4σ0 and find the boundaries of the
excursion set. They are exhibited as black/red lines in Figure 5.
We test our algorithm by applying it to a Gaussian random From these boundaries, we construct the Minkowski func-
field. For such a field, the Minkowski functionals can be tionals and tensors.
calculated analytically, and we can also compare the shape of The scalar Minkowski functionals are exhibited as a function
the field in the vicinity of peaks to known analytic results. of normalized density threshold ν in Figure 6. We generate
We generate a two-dimensional Gaussian random field δk in Nreal=100 realizations of a Gaussian random field; the blue
Fourier space with a constant power spectrum (Gaussian white points are the mean of these realizations, obtained using our
noise). This field is then converted to its real-space counterpart algorithm. The error on the mean is smaller than the points. The
via a fast Fourier transform algorithm. We generate the field solid black line is the theoretical expectation value. The
over a 3150 ´ 3150(h-1 Mpc)2 area, adopting a 2048×2048 accurate reproduction of the theoretical curves serves as a
equispaced grid over this range, yielding a resolution consistency check of our method.
 = 1.54 h-1 Mpc. We smooth the field in the plane with a In Figure 7, we exhibit the matrix components of (W1,12 )ij. In
Gaussian kernel, using a smoothing scale RG = 15 h-1 Mpc. the top panel, we exhibit the mean and error on the mean of
5
The Astrophysical Journal, 858:87 (15pp), 2018 May 10 Appleby et al.
Figure 7. Top panel: matrix components of the Minkowski tensor (W1,12 )ij;
(i, j) = (1, 1) (blue), (1, 2) (yellow), and (2, 2) (red). The solid black curve is
the theoretical prediction for W1/π, which should match the diagonal
components. Bottom panel: (i, j) = (1, 1) (blue) and (2, 2) (red) components Figure 8. Statistic β for holes and connected regions (top and middle panels)
of the fractional residuals defined in Equation (17). These quantities are and the combined mean of both holes and connected excursion regions (bottom
consistent with zero for all thresholds probed. panel). The black circles and green squares indicate the mean value of this
statistic after we have made different cuts to our sample of connected regions
diagonal components (i, j) = (1, 1), (2, 2). We also show the and holes (according to the size of the region) for a fixed smoothing scale
R = 15 h-1 Mpc. We note that the statistic is insensitive to these cuts. The red
theoretical prediction for W1/π, which should match the
G
squares are the same statistic with a smaller smoothing scale R = 10 h-1G Mpc
diagonal components for an isotropic Gaussian random field. applied. In this case, the percolation of the field occurs at a slightly different
We find close agreement between the isotropic expectation value of ν, but the statistics behave similarly otherwise. We also exhibit the
value and our numerical reconstruction. We exhibit the off- theoretical prediction (Equation (18)) obtained using peak statistics as a solid
diagonal component of the matrix, (i, j) = (1, 2), nding black curve in the bottom panel. The theoretical curve approaches ourfi numerical result for high threshold values.
consistency with zero.
In the bottom panel, we exhibit the fractional difference
In what follows, we use the quantities D(W1,12 ) ,D(W1,111 2 )22,
( 1,1) and (W1,1 1,1p W - W  2 )12 áW2 ñ to study the sensitivity of the statistic W1,12
D( 1,1) º 2 ij 1 ijW2 ij , (17) to galaxy bias, gravitational evolution, and redshift-space
W1 distortion, where áW1,1ñ = [(W1,1) + (W1,12 2 11 2 )22] 2 is the aver-
age of the diagonal components of the matrix. We note that the
which should be consistent with zero for an isotropic Gaussian functionsD(W1,1) ,D(W1,1) , and (W1,1) áW1,1ñ will not be
field. The error bar increases with ∣n∣ 11 22 12due to the smaller 2 2 2 2Gaussian-distributed but will be symmetric with respect to the
perimeter of the excursion set, leading to larger statistical peaks of their probability distributions.
fluctuations. Note that from the definitions in Equations (2) We exhibit the mean and error on the mean of β in Figure 8,
and (15), the sum D(W1,1) 1,12 11 + D(W2 )22 must be zero. generated from the Nreal=100 realizations. The top and
6
The Astrophysical Journal, 858:87 (15pp), 2018 May 10 Appleby et al.
middle panels show the statistic for holes βh and connected Table 1
regions βc, respectively, and the bottom panel shows the Fiducial Parameters Used in the Horizon Run 4 Simulation
average for the combined holes and connected regions βtot. The Parameter Fiducial Value
number of distinct connected regions and holes used to
calculate the averages varies greatly as a function of ν and is Ωmat 0.26
related to the Betti numbers (Chingangbam et al. 2012; Park ΩΛ 0.74
n 0.96
et al. 2013). For the smoothing scales and volumes probed in s
σ8 0.794
this work, we have ~(103) distinct connected regions/holes
at ν±1. The function β is only defined in the domain
0<β„1; therefore, it is not a Gaussian-distributed variable. been constructed by analytically calculating the mean shape of
However, we have checked that the probability distribution the peaks of a two-dimensional Gaussian field. This calculation
function of β is not significantly skewed and the mean is an has been performed previously in Bond & Efstathiou (1987;
appropriate proxy for its peak. see also Aurich et al. 2011 for later work and Bardeen et al.
In Figure 8, we exhibit the statistic β after making various 1986 for the three-dimensional case), and we quote the results
cuts to the excursion-set sample. To minimize spurious relevant to our analysis in Appendix B. In the large ν threshold
numerical artifacts, we have adopted a highly resolved plane
( -1 )2 limit, ál2 l1ñc is expected to beof size 3150 ´ 3150 h Mpc with a total number of pixels
Npix = 2048 ´ 2048 and a smoothing scale RG = 15 h-1 Mpc. 2pò cos2 f (k2 - (k2 - 1)cos2m m f)-3 2 dfWith this choice, we smooth over nearly 10 pixels. To further ál l ñ = 02 1 c ,
test for numerical artifacts, we make cuts to our sample. 2pò cos2 f (1 - (1 - k2m)cos2 f)-3 2 df
Specifically, the black circles and green squares in Figure 8 0
represent the mean ál2 l1ñ from a sample with Acut = 0, 4 2, (18)
where  = 1.54 h-1 Mpc is the pixel size and Acut is the area where κm is defined in Equation (34) and is related to the
cut that we apply to the excursion regions. So, for the black expectation value of the ellipticity of a peak e , defined in
points, we use the entire sample to calculate ál2 l1ñ, and for the
m
green squares, we remove all excursion regions (holes and Equation (33). In deriving Equation (18), one assumes that the
connected regions) that have an area A < 4 2 before contours of constant density in the vicinity of a peak are
calculating ál l ñ. As discussed in Appendix A, we apply elliptical. For a Gaussian field, Equation (18) is valid for βtot at2 1
area cuts to test that no spurious anisotropic signals are large ν, as in this regime βtot;βc. Similarly, due to the
generated as a result of including poorly resolved excursion ν→−ν symmetry of a Gaussian field, Equation (18) is also
subsets in the average quantities bc,h,tot. We find that the valid for βh at extreme negative ν values. We stress that the
statistics are practically independent of any area cut that we excursion-set boundary will only trace the peaks and troughs of
impose, indicating that the well-resolved objects dominate our the field in the high ∣n∣ threshold limit. There is no general
sample for the thresholds probed. If poorly resolved regions correspondence between peaks and connected regions or
become dominant, then at high ∣n∣ one would observe a troughs and holes.
spurious decrease in bc,h,tot. We also exhibit the same statistics
for a field smoothed on a smaller scale RG = 10 h-1 Mpc: we 5. Minkowski Tensors Applied to Simulated Galaxy
observe that βtot is insensitive to RG. As discussed in Catalogs
Appendix B, this result is expected for a Gaussian white-noise
field. We now consider the Minkowski tensors of the low-redshift
The top and middle panels of Figure 8 present very different dark matter density field and study the effect of galaxy bias,
1,1
behavior on either side of ν=0; this is due to the fact that we gravitational evolution, and redshift-space distortion on W2
initially have a single hole (or connected region) with structures and βtot. We apply our statistics to the Horizon Run 4
embedded. In this regime—the right- and left-hand sides of the simulation data. Before continuing, we briefly describe the
top and middle panels, respectively—the mean value of β is simulation.
dominated by a single region that is roughly the size of the Horizon Run 4 is the latest data release from the Horizon
7
entire plane.6 This single excursion region undergoes rapid Run project. It is a dense, cosmological-scale N-body
percolation into many structures, which is exhibited by the simulation that gravitationally evolved N = 6300
3 particles in
-1 3
rapid change at ν=±1 in the figures. Following this, the a V = (3150 h Mpc) volume box. The cosmological para-
statistic β is dominated by distinct holes for ν <−1 and β meters used can be found in Table 1, and the details of theh A c
by distinct connected components for ν >1. simulation are discussed in Kim et al. (2009, 2015). We useA
Regardless of the A that we use, the mean shape of each two-dimensional slices of snapshot data at z=0.2 of thicknesscut
individual connected region and hole has a value β ∼0.6 Δ. The field is smoothed in the plane of the data using atot
(bottom panel; Ganesan & Chingangbam 2017), which Gaussian kernel of width RG. We vary both Δ and RG in what
increases with ∣n∣, as we expect. This suggests that the mean follows.
shape is becoming increasingly circular with increasing density Rather than use the dark matter particle data, we adopt the
threshold, but β remains significantly smaller than unity even mock-galaxy catalog constructed in Hong et al. (2016). Mocktot
at large ∣n∣. galaxies are assigned by the most bound halo particle–galaxy
In the bottom panel, we exhibit the theoretical prediction correspondence scheme. The survival time of satellite galaxies
(Equation (18)) as a solid black line. The theoretical curve has after merger is calculated by adopting the merger timescale
model described in Jiang et al. (2008). We take a fiducial
6 When calculating βtot, we do not include connected regions or holes that
have sizes of the same order of magnitude as the total area of the plane. 7 http://sdss.kias.re.kr/astro/Horizon-Runs
7
The Astrophysical Journal, 858:87 (15pp), 2018 May 10 Appleby et al.
galaxy number density of n¯ = 10-3(h-1 Mpc)-3 by applying a using galaxies as tracers of the underlying field, this effect will
lower-mass cut. be intertwined with galaxy bias. Fixing a constant galaxy
From the galaxy distribution, we generate a density field by number density at each redshift generates a galaxy distribution
first generating a regular grid of size 2048×2048 in the x1, x2 with a bias that is roughly constant with redshift. On the other
plane and slices of width Δ in the x3 direction (taken as the line hand, a more highly biased point distribution will better trace
of sight). We bin the mock galaxies according to the x3 slice to the high peaks of the underlying density field, which will be
which they belong and then in the two-dimensional pixelated more spherical. Hence, we can expect the tilt in βtot as a
grid according to a cloud-in-cell scheme. Taking each slice in function of νA to increase with increasing bias.
turn, we use the average number of galaxies n̄ per pixel to In the left panels of Figure 10, we exhibit D(W1,12 )11,
define the two-dimensional density field dij = (nij - n¯) n¯, D(W1,12 ) , (W1,1 1,122 2 )12 áW2 ñ, and βtot at three epochs. We take a
where i, j indices run over the 2048×2048 lattice. Next, the Gaussian random field with a linear ΛCDM dark matter power
slice is smoothed over the plane using a two-dimensional spectrum as the initial condition of the simulation and calculate
Gaussian of width RG. For each slice, we calculate the β andW1,1tot 2 for this field and for the Horizon Run 4 snapshot
Minkowski tensor W1,12 and βtot. Rather than use the boxes at z=1 and 0.2, fixing (D, RG) = (30, 15)h-1 Mpc.
conventional overdensity threshold ν to define the excursion We use all galaxies in the simulation, fixing the number
set, we instead adopt the area threshold νA parameter, which is density n¯ ~ 1.5 ´ 10-2(h-1 Mpc)-3.
defined as We exhibit D(W1,12 )11 and D(W1,12 )22 at different redshifts,
¥ finding no evidence of evolution. Similarly, the off-diagonal1
f = ò exp[-t2 2]dt, (19) component (W1,1A 2 )12 remains consistent with zero. This implies2p nA that the relationshipW1,12 µ W1 is not affected by gravitational
evolution. However, both (W1,12 )ij and W do evolve withwhere fA is the fractional area of the field above νA. The 1νA redshift: the scalar Minkowski functional W1 is skewed due to
parameterization eliminates the non-Gaussianity in the one- the non-Gaussianity generated by the effect of gravity
point function (Gott et al. 1987; Weinberg et al. 1987; Melott (Matsubara 1994, 2003). They evolve in such a way that the
et al. 1988). This choice allows us to compare excursion sets in relationship W1,12 µ W1 is preserved.
the Gaussian and non-Gaussian fields that occupy the Here βtot becomes increasingly tilted relative to its Gaussian
same area. form with decreasing redshift. In Figure 10, we exhibit both βtot
and the residual Δβtot, which is the fractional difference
5.1. Varying Smoothing Scales RG, Δ between βtot as measured from the galaxy catalogs at
z=1,0.2 and the initial condition. Here Δβ varies
In the top panels of Figure 9, we exhibit the quantities tot
D(W1,1) , D(W1,1) approximately linearly with νde ned in Equation (17) for the z=0.2 A and is Δβtot;0.1 for high2 11 2 22 fi thresholds ∣nA∣  4. The increasing signal with time indicatesHorizon Run 4 mock-galaxy density field. In the left panels, we
D = -1 that overdense patches of fixed area become increasinglyfix the slice thickness 30 h Mpc and vary the Gaussian
-1 circular as collapse occurs. The underdense regions ν <0 ofsmoothing scale in the plane RG = 20, 15, 10 h Mpc (green,
A
the same area become less spherical relative to the initial
yellow, and blue points). In the right panels, we fix
R = 15 h-1 −1
condition.
G Mpc and vary Δ=40, 30, 20 h Mpc (cyan, In the right panels, we plot D(W1,1) , D(W1,1) ,
yellow, and red points). The error bars are constructed as the 2 11 2 221,1 1,1
error on the mean calculated using N =75 slices of the eld. (W2 )12 áW2 ñ, and βtot for the z=0.2 snapshot data, takingslice fi
BothW1,1 and W are reconstructed from the data. In the middle different mass cuts to the galaxy distribution to generate galaxy2 1
( 1,1) 1,1 ( 1,1) ( 1,1) catalogs with number density n̄ = 5.0 ´ 10
-3, 1.0×10−3, and
panels, we exhibit W2 12 áW2 ñ. Here D W2 11, D W2 22, 5.0×10−41,1 1,1 (h
−1Mpc)−3. One can observe no significant depend-
and (W2 )12 áW2 ñ are all consistent with zero, which means ence of mass cut and number density onD(W1,12 )11,D(W1,12 )22, or
that the relationW1,12 µ W1 remains true for the gravitationally (W1,1) áW1,12 12 2 ñ; however, as we decrease the galaxy numberevolved nonlinear density field. density, the value of βtot decreases for νA<0. We also exhibit
We exhibit βtot in the bottom panels of Figure 9. This Δβtot, which is the fractional difference between βtot as measured
quantity is sensitive to both Δ and RG. The most significant with the n¯ = 1.0 ´ 10-3, 5.0 ´ 10-4(h-1 Mpc)-3 samples and
effect of gravitational evolution on βtot is in the large νA n¯ = 5 ´ 10-3(h-1 Mpc)-3. The change in number density n̄
regime, where overdensities become increasingly spherical due affects the shape of βtot(νA) predominantly in the νA<0 regime.
to gravitational collapse. In contrast, underdense regions One can observe that both gravitational collapse and galaxy bias
characterized by νA<0 become less spherical for excursion affect the shape of the βtot(νA) curve similarly.
sets of fixed νA. The tilt in βtot(νA) indicates that the holes are
less circular than those in a Gaussian field occupying the same
area, and the overdensities are more circular. The βtot decreases 5.3. Redshift-space Distortion
for negative thresholds νA<0 as RG is lowered but is only
> Our results indicate that the statisticweakly sensitive to for 0. D(W
1,1
2 ) is insensitive toRG νA gravitational collapse, and this remains true regardless of our
choice of Δ and RG smoothing scales. The matrix retains the5.2. Redshift Evolution and Galaxy Bias
relation W1,12 µ W1 , as the effect of gravity introduces no
One can study the redshift evolution of W1,12 and βtot by preferred direction. However, as stated in the Introduction, the
calculating these statistics for slices of snapshot data at dark matter field that we observe via galaxy tracers is not
different redshifts. One should observe an initially symmetric isotropic: a preferred direction exists due to the redshift-space
βtot at high redshift, which becomes increasingly tilted due to distortion effect along the line of sight. We close this section by
gravitational collapse with decreasing z. However, as we are considering how the global properties of the field are modified
8
The Astrophysical Journal, 858:87 (15pp), 2018 May 10 Appleby et al.
Figure 9. Top panels: fractional differencesD(W1,12 )11,D(W1,12 )22 defined in Equation (17). We exhibit these quantities as a function of nA, which is related to the area
fraction of the field. The dependence ofD(W1,12 )11,D(W1,12 )22 onΔ and RG is negligible. Middle panels: off-diagonal component (W1,12 )12 divided by the average of the
two diagonal components áW1,12 ñ. This quantity is consistent with zero for all RG and Δ values. Bottom panels: btot calculated using all connected regions and holes at
each n threshold. Left panels: D = 30 h-1 Mpc, R = 20, 15, 10 h-1 Mpc (green, yellow, and blue). Right panels: R = 15 h-1 Mpc, D = 40, 30, 20 h-1A G G Mpc
(cyan, yellow, and red). Error bars denote the error on the mean from Nslice = 75 fields.
when we introduce a preferred direction to the data. For this the plane is perpendicular to the line of sight, and θlos=0
purpose, we take the z=0.2 snapshot data and apply a corresponds to a density plane aligned exactly with the line of
redshift-space distortion to the position of each galaxy by sight. The introduction of the velocity correction to the galaxy
adjusting their position in the x2, x3 directions via the relation positions generates a global anisotropy in the field, which the
(1 + z) Minkowski tensor W1,12 is sensitive to.x2¢ = x2 + v2 cos qlos, (20)
H (z) In Figure 11, we exhibit D(W
1,1
2 )11, D(W1,12 )22 (top panel),
(W1,1) 1,12 12 áW2 ñ (middle panel), and βtot (bottom panel) for the
¢ (1 + z)x = x + v sin q , (21) real-space field (yellow squares) and slices of the redshift-space3 3 3 ( ) losH z distorted field aligned perpendicular (red circles) and parallel
(green circles) to the line of sight. We fix Δ=30 h−1 Mpc and
where v2,3 are the velocities in the x2,3 directions and θlos is the RG=15 h
−1 Mpc. These values were chosen to ensure that the
angle of the data plane relative to the line of sight. We always field is in the mildly nonlinear regime in both smoothing
generate data slices along the x3 axis, so varying 0„θlos„ planes. In the bottom panel, we exhibit both βtot for the three
π/2 is equivalent to varying the observer line of sight with cases and the fractional residuals Δβtot between the real-space
respect to the plane. Here θlos=π/2 is the standard case where value of βtot and the redshift-space distorted values (so, for
9
The Astrophysical Journal, 858:87 (15pp), 2018 May 10 Appleby et al.
Figure 10. Left panels: redshift evolution of the statistics D(W1,12 )11, D(W1,12 ) 1,122, (W2 )12 áW1,12 ñ (top and middle), and btot (bottom) for the Horizon Run 4 snapshot
data at z=0.2 and 1 (yellow and gray) and the Gaussian initial condition (white). We have used fiducial smoothing parametersD = 30 h-1 Mpc, R = 15 h-1G Mpc.
All galaxies in the simulation are used as density tracers at z=1 and 0.2, with a total number density n¯  1.5 ´ 10-2(h-1 Mpc)-3. Here Dbtot is the fractional
difference between btot as measured at z = 1, 0.2 and the initial condition. Right panels: D(W1,12 )11, D(W1,1) 1,1 1,12 22, (W2 )12 áW2 ñ, and btot for z=0.2 snapshot data,
taking different mass cuts to the galaxy sample to yield number density n̄ = 5.0 ´ 10-3, 1.0×10−3, and 5.0×10−4 (h−1 Mpc)−3 (blue, yellow, and green). Here
Dbtot is the fractional difference between b -3tot measured using n¯ = 1.0 ´ 10 , 5.0 ´ 10-4(h-1 Mpc)-3 galaxy catalogs and the most dense
sample n¯ = 5.0 ´ 10-3(h-1 Mpc)-3.
example, the green circles in the lower panel represent the ofD(W1,12 ) 1,111,D(W2 )22 from zero when measured in either real
fractional residual Db = (btot,rsd (qlos = 0) - btot,real) btot,real). or redshift space with qlos = p 2. The shape of individual
The effect of redshift-space distortion is markedly different excursion regions as described by βtot is modified by ∼1% but
for the two planes. If we align the data plane perpendicular to is not systematically higher or lower than its real-space value.
the line of sight (red circles), then the effect of linear redshift- When we align the data slice parallel to the line of sight, the
space distortion is to increase the density contrast, as galaxies effect of peculiar velocities will be to increase the ellipticity of
in the vicinity of the slice boundary will be scattered into/out overdensities. In contrast to the θlos=π/2 case, the effect of
of the slice for overdense/underdense regions. The shapes of redshift-space distortion will now generate a globally preferred
connected regions and holes will change, but their orientations direction in the excursion-set boundary. The effect on
will remain random. As such, we can expect W1,12 to remain individual excursion regions is small. In the bottom panel of
insensitive to redshift-space distortion when we take Figure 11, we observe the fractional change Δβtot in βtot as
θlos=π/2. This agrees with our numerical result. In the top measured in real space and the θlos=0 plane in redshift space
panel of Figure 11, we find no statistically significant departure (green circles). Here Δβtot is negative in the range
10
The Astrophysical Journal, 858:87 (15pp), 2018 May 10 Appleby et al.
nonlinear fingers-of-God effects, which modify the shape of
overdensities parallel to the line of sight, are also contributing
to the signal. The asymmetry about νA=0 also implies that
the W1,1 matrix no longer satis es the relation W1,12 fi 2 µ W1—
non-Gaussianity of the velocity field affects W1 and W1,12
differently, and additional information can be extracted by
measuring both.
Here W1,12 is sensitive to redshift-space distortion and not
gravitational collapse because the latter effect is statistically
isotropic, in principle, at any scale. The density field will
undergo collapse, but no preferred direction will be generated
in the excursion-set boundary in real space. This makes the
Minkowski tensor an ideal candidate to measure the large-scale
properties of the velocity field.
We expect that the linear Kaiser effect will generate a
constant shift in DW1,12 and the fingers of God a tilt as a
function of νA. It follows that measurements of W1,12 can be
used to simultaneously constrain the redshift-space distortion
parameter b = f b and the velocity dispersion of gravitation-
ally bound galaxies. The next stage of this analysis requires a
theoretical prediction of the Minkowski tensors in redshift
space. A real-space analysis was conducted in Chingangbam
et al. (2017); the generalization to redshift space will be
considered elsewhere.
6. Summary
In this work, we have studied the morphological properties
of two-dimensional density fields. For this purpose, we have
adopted the Minkowski tensor W1,12 . To use this statistic, we
must first generate a bounding perimeter of constant density,
which defines an excursion set. We adopted the method of
marching squares, the details of which are described in the text.
We studied the diagonal and off-diagonal elements ofW1,12 for a
Gaussian random field, finding that this matrix is proportional
to the identity matrix and the scalar Minkowski functional W1.
We then considered the W1,12 statistic applied to individual
subsets of the excursion set. For every threshold ν, we
constructed the matrix W1,12 for each distinct connected region
and hole, and from them we extracted the eigenvalues λ1,2.
Figure 11. Statistics D(W1,12 )11, D(W1,1) , (W1,1 1,12 22 2 )12 áW2 ñ (top and middle These quantities inform us of the shape of individual excursion-
panels), and btot (bottom panel) as a function of nA for the Horizon Run 4
= set regions. We calculated the mean eigenvalue ratio ál l ñ asz 0.2 snapshot box, where we have introduced a redshift-space distortion 2 1
along the line of sight. Here q a function of ν and, in the large ∣n∣ limit, related this quantity tolos is the angle of the line of sight relative to the
plane of the data, with qlos = p 2 being the usual case where the plane is the mean ellipticity of the field in the vicinity of a peak. We
perpendicular to the line of sight. The yellow squares represent the statistics in found reasonable agreement between theory and numerical
real space, and the red/green circles represent the redshift-space distorted fields application of our algorithm in the large threshold limit.
perpendicular/parallel to the line of sight. The shape of individual objects btot
also exhibits a dependence on q , with excursion subsets becoming less However, the statistic βtot is a more general measure of shapelos
circular in redshift space relative to the real-space field for qlos = 0. In the than the ellipticity; it is a property of the excursion-set
lower right panel, we exhibit the fractional residual Dbtot between btot as boundary and makes no assumption regarding its shape.
measured in real and redshift space. Finally, we applied the Minkowski tensor to mock-galaxy
data and considered how it is modified by gravitational
−3<νA<3, which indicates that structures in real space are collapse. We found that the mean eigenvalue ratio βtot is
more spherical; however, the difference is a roughly ∼1% particularly sensitive to the effect of gravity, the dominant
effect that slowly increases with increasing νA. effect being a tilt that indicates that connected regions become
Although the effect on each individual excursion-set region increasingly circular relative to holes occupying the same area.
is small, it is coherent in the sense that, statistically, all In contrast, the matrixW1,12 defined over the entire excursion set
overdensities/underdensities will be distorted in the same is essentially insensitive to gravitational collapse, as the process
direction. AsW1,12 is a measurement of the preferred directions introduces no preferred direction.
in the global excursion-set bounding perimeter, the distortion However, when the data contain a large-scale anisotropic
generates a cumulative signal in this statistic. In the top panel signal, W1,12 will exhibit strong sensitivity. When we corrected
of Figure 11, we observe this effect: D(W1,12 )22 exhibits an mock-galaxy positions to account for redshift-space distortion
∼8% departure from the isotropic limit, tilting as a function of and repeated our analysis using slices of the density field
νA. The asymmetry of D(W1,12 )22 about νA=0 indicates that oriented by angle θlos relative to the line of sight, we found that
11
The Astrophysical Journal, 858:87 (15pp), 2018 May 10 Appleby et al.
the diagonal components of W1,12 are significantly modified: a
distinctive functional dependence on nA develops with over-
densities preferentially aligning along the line of sight. The
anisotropy manifests as both a change in amplitude ofD(W1,12 )ii
and a roughly linear dependence on νA. The fact that the
statistic is sensitive to anisotropy in the data, and only very
weakly to the non-Gaussianity of the late-time field, makes it a
promising candidate to study the velocity perturbations in
redshift space. We consider the redshift-space theoretical
expectation of the Minkowski tensors in a forthcoming Figure 12. Example of the ambiguity implicit within our algorithm. Any
publication. density peak or nonmonotonic behavior of the field over scales smaller than our
To observe this signal using real data, we require a resolution ò will not be detected. Here we exhibit two distinct cases of a density
measurement of the density field in planes parallel to the line field: in both, the four vertices at which we measure dij are “out” of the
of sight. Upcoming galaxy surveys such as DESI (Aghamousa excursion region, but in the right panel, there is a maxima internal to the square.
Our algorithm can never detect such small-scale features and will always select
et al. 2016) and LSST (Abell et al. 2009) will provide volume- the left case in this example.
limited galaxy samples over gigaparsec volumes from which
we can take subsets of the field perpendicular and parallel to
the line of sight. Existing surveys such as HectoMAP (Geller once. This is tantamount to the statement that the density field
et al. 2011; Geller & Hwang 2015) provide spectroscopic is monotonically increasing or decreasing on the scale of our
galaxy catalogs over cosmological scales in slices parallel to spatial resolution ò. As a result, our method will not be able to
the line of sight and can be used to extract the redshift-space distinguish certain examples of the density field; for example,
distortion signal predicted in this work. in Figure 12, we exhibit two squares that cannot be
Photometric redshift uncertainties will be the dominant distinguished, and our algorithm will always adopt the left
source of contamination to the signal, as they will also scatter panel. The right panel shows a density peak internal to the
galaxy positions along the line of sight. As the redshift-space square as a gray shaded area; we cannot reconstruct such a peak
distortion effect is present when smoothing over large scales
~ -1 using marching squares. Critical points are manifestly higher-RG 15 h Mpc, spectroscopic catalogs will be better suited
1,1 order phenomena, which cannot be modeled via linearto measuring W2 . However, upcoming photometric catalogs interpolation.
can still potentially be used to extract information from W1,12 . Although our method will miss this small-scale behavior, we
The fingers of God introduce a tilt inW1,12 as a function of νA, always smooth the field over at least three pixel lengths. The
as galaxies in overdense regions will predominantly experience smooth field will generically be monotonic over scales~( ).
the effect. It follows that redshift-space distortion and However, the extremes of the distribution (the high threshold
photometric redshift uncertainty could potentially be disen- peaks, for example) are likely to occupy a small area, and we
tangled, as the latter will not possess the same sensitivity to will inevitably fail to detect some of these objects. With the
density fluctuations and hence will not generate the same tilt in Horizon Run 4 mock-galaxy data, we can test the significance
DW1,12 (nA). Photometric redshift catalogs are characterized by of this issue using the following method.
large number densities and cosmological-scale volumes, and The mock-galaxy data are a point distribution. We take the
using them will provide better statistics relative to spectro- three-dimensional Horizon Run 4 mock galaxies and bin them
scopic samples. A detailed study of photometric redshift into two-dimensional slices of thickness Δ, as before. We then
contamination will be presented elsewhere. generate two grids xi, y ¢ ¢j and xi , yj = xi +  2, yi +  2 in
the two-dimensional plane and bin the galaxies according to a
The authors thank the Korea Institute for Advanced Study cloud-in-cell scheme for each grid. The resulting density fields
for providing computing resources (KIAS Center for Advanced are denoted δij and d¢ij, respectively. For the density field δ , weComputation Linux Cluster System) for this work. This work ijperform the marching-squares algorithm as described in the
was supported by the Supercomputing Center/Korea Institute main body of the text, but we now perform an additional check
of Science and Technology Information, with supercomputing whenever cases Nc=1 or Nc=16 (displayed in Figure 1) areresources including technical support (KSC-2013-G2-003), and encountered.
the simulation data were transferred through a high-speed Our algorithm will always predict δ<ν and δ>ν for the
network provided by KREONET/GLORIAD. central values of Nc=1 and Nc=16, respectively—consistent
with no small-scale structures on scales~( ). This is because
Appendix A we use a simple linear interpolation scheme to predict the
Sources of Numerical Error density between the (i, j) pixels. We now test the center of these
There are two issues with the marching-squares algorithm boxes by using di¢j as the value of the density field at xi+ò/2,
that are capable of generating spurious numerical artifacts, yj+ò/2. Whenever we encounter the case Nc=1, we check if
which we briefly discuss. di¢j > n . If this inequality is satisfied, then we can say that there
is some structure on the scale of our resolution ò that the
algorithm has failed to detect. Similarly, for the case N =16,
A.1. Topological Ambiguity Associated with Marching-squares cwe test di¢j < n , in which case there is a hole in the excursionAlgorithm set that has not been detected. We count the total number of
The first problem is our choice of interpolation scheme. We holes and connected regions that the code fails to detect in the
are assuming that in/out states joined along the edges of entire plane at each density threshold ν. We then divide this
squares in the δij grid will always cross the threshold δ=νσ0 number by the total number of holes and connected regions that
12
The Astrophysical Journal, 858:87 (15pp), 2018 May 10 Appleby et al.
Figure 13. Fraction of missed connected regions and holes fmissed as a function
of nA for three different box resolutions, ò=1.54, 3.00, and 6.15 h−1 Mpc
(yellow squares, blue circles, and green triangles). This statistic informs us of Figure 14. Example of the boundary that the marching-squares algorithm will
spurious numerical errors in our reconstruction of the genus due to the generate for a small excursion-set region of the order of the size of a single
marching-squares algorithm. The number of missing structures is negligible for pixel ~( 2). The discrete nature of our algorithm generates a polygon that
our choice of residual resolution ò=1.54 h−1 Mpc but increases sharply with will not accurately represent the true, smooth boundary.
increasing ò.
pixel size, as the reconstructed bounding perimeter for these
the code successfully finds during the course of the algorithm; objects will not accurately represent the smooth contour of the
we denote this fraction fmissed. Efficacy of the method requires continuous field d (x1, x2). As an example, in Figure 14, we
fmissed = 1. exhibit an excursion set represented by a single point in our
We repeat this calculation for Nslice=75 slices of the three- grid. The boundary is exhibited as a polygon, but the difference
dimensional density field and calculate the average fraction between the discretized perimeter and the smooth underlying
fmissed as a function of νA. We exhibit this quantity in field is likely to be large in this instance. This difference will
Figure 13. We repeat the calculation for three different spatial lead to large errors in our numerical reconstruction of the
resolutions, ò=1.54, 3.00, and 6.15 h−1 Mpc (yellow squares, morphological properties of the field whenever small excur-
blue circles, and green triangles), fixing the smoothing scale sion-set regions dominate the total excursion region. This is
RG=15 h
−1 Mpc in the plane. We observe that the fraction of likely to occur at large threshold values ∣n∣.
missed structures is negligible for our fiducial resolution We can estimate the magnitude of this discretization effect
ò=1.54 h−1 Mpc, but the effect increases with ò. Additionally, directly. To do so, we take a regular grid and generate a smooth
the fraction of missed structures increases with increasing ∣nA∣. circular density field. Defining the center of a circle
This is to be expected, as the typical area occupied by peaks rcen = (x1,cen, x2,cen), we define d (x1, x2) as
and minima decreases at large thresholds. The total number of
missed holes and connected regions is between 5% and 10% d (x1, x2)
d
= cen , (22)
when using ò=6.15 h−1 Mpc, corresponding to smoothing 1 + (x 2 21 - x1,cen) + (x2 - x2,cen)
over 2.4 pixels. We must smooth over at least five pixel lengths
to ensure that the number of structures missed by the algorithm where δcen is an arbitrary constant. For this field, surfaces of
is ∼1%. constant d (x1, x2) = dc will generate an excursion set of
It is difficult to provide a physical interpretation of fmissed, as constant radial distance from rcen. Increasing the threshold δc
our test does not reveal all cases in which the algorithm can will generate smaller circular excursion sets. For a circle, we
fail. For example, we have only calculated the density field in can trivially calculate all properties of the excursion set, and
the exact center of the squares; peaks of the field may occur at specifically, we have βc=1. Therefore, as we decrease theany point. Furthermore, all 16 cases in Figure 1 can exhibit radius of the circular density field, we can ascertain the extent
nonlinear behavior of the field on scales of order ò that can
modify the genus, and we have only considered failures to which our numerical reconstruction of β deviates from unity.
associated with N =1, 16. However, we can argue that the In Figure 15, we exhibit the residual Δβ=β−1 as a functionc
center of the boxes N =1, 16 are most likely to exhibit of the radius of the circular region in units of pixel size: r/ò.c
irregularities (being maximally distant from our interpolation The blue circles and error bars denote the average and rms
points); hence, fmissed provides a conservative indicator of the fluctuations of Nreal=100 realizations, randomly placing the
failure rates in all boxes. We conclude that the fiducial center of the circular density rcen within the box.
smoothing and resolution scales adopted in this work are One can observe close agreement between theory and
sufficient to minimize this particular spurious numerical numerical approximation when the circle is well resolved,
artifact. r>4ò, but order ~(10%) discrepancies are apparent for
r∼ò. This numerical artifact will artificially decrease the
isotropy of the density peaks. However, when we calculate
A.2. Finite Resolution Effect bc,h,tot for a stochastic density field, we should not simply
A second source of numerical contamination arises for remove poorly resolved excursion-set regions from the sample
excursion-set regions occupying an area of the order of the and calculate the mean properties of the remaining set, because
13
The Astrophysical Journal, 858:87 (15pp), 2018 May 10 Appleby et al.
Equation (23) in terms of its eigenvalues ω1,2,
d ( 1ri) = d (0) - å w 2i ri , (24)
2 i
where ri is a coordinate basis in this rotated frame. We fix
ω1>ω2 in what follows. A surface of constant δ(r)=δc
corresponds to an ellipse with axes
⎡ 2(d (0) - d ) ⎤1 2
ai = ⎢ c⎣ ⎦⎥ . (25)wi
The ratio of the axes of the ellipse a2/a1 is then simply given by
a2 w= 1 , (26)
a1 w2
and we define the ellipticity as
Figure 15. Fractional departure of the statistic β from its theoretical w1 - we = 2 . (27)
expectation value β=1 for a single circular density region as a function of 2(w1 + w2)
the ratio of the radius of the circle to the resolution of our two-dimensional
grid. For circular regions that are not well resolved, r<3ò, we find errors For a Gaussian random field d (x1, x2) with power spectrum>10% in our numerical reconstruction of the statistics.
P2D(k), one can calculate the conditional probability of the
ellipticity e of a peak of height n ; it is given by (Bond &
doing so could introduce a selection bias. The shape of an p
excursion set is correlated with its area: large regions are more Efstathiou 1987)
likely to be less circular. Making size cuts to our sample will (1 - 4e2)8ede e-4x 2* e2
bias the resulting β statistic. So we have two competing effects: P (e∣np)de =
if we simply calculate the mean statistic ál2 l1ñ for all regions 1 + 8(1 - g2)e2 G (g, x*)
in our sample, then we will observe a spurious anisotropic ⎡ ⎛ ⎞⎤
signal at high ∣n∣ threshold values where the peaks will typically ⎢ 1 x´ ⎢1 - erfc ⎜ * ⎟⎥,be small. However, if we make an area cut and eliminate small ⎣ 2 ⎝⎜ 2(1 - g2)( ⎟1 + 8(1 - g2)e2) ⎠⎦⎥
regions from the average, then we will be preferentially
selecting large excursion regions in our sample. In the main (28)
body of the text, we vary the cut and examine its effect on our ⎡ ⎛ ⎞⎤
statistics. As we smooth the field over an increasingly large G (g, x 2 2 ⎢ 1 x*) = (x* - g )⎢1 - erfc ⎜ * ⎟⎥number of pixels, this numerical artifact will become less ⎣ 2 ⎝⎜ 2(1 - g2) ⎠⎟⎦⎥
significant, and we find that the effect is negligible if we
2 2
smooth over several pixel lengths. e-x* [2(1-g )]2
We have repeated the above test using elliptical density + x*(1 - g ) 2p (1 - g2)
fields randomly located within a two-dimensional plane. We
again find percent-level agreement between our numerical e-x
2
* (3-2g
2) ⎡ ⎤1 ⎛ x ⎞
algorithm and analytic predictions, subject only to the + ⎢⎢1 - erfc2 ⎜⎜ * ⎟⎥2 2 ⎟⎥,
condition that the ellipse is well resolved (with minor axis 3 - 2g ⎣ 2 ⎝ 2(1 - g )(3 - 2g ) ⎠⎦
e>4ò). This indicates that our numerical error will not be (29)
sensitive to the shape of the excursion-set regions.
where g = s21 s2s0, x* = gnp, and s0,1,2 are cumulants of the
density field
Appendix B kdk
Mean Shape of a Peak in a Two-dimensional s2 º ò P (k)k2jj 2D . (30)
Gaussian Field 2p
A peak in a two-dimensional eld d (x To predict the mean ellipticity of all peaks above the densityfi 1, x2) can be
characterized by its height, which we define as ν , and its threshold ν, we also require the number density of peaks  ;p p
ellipticity e. In the vicinity of a peak at xi=0, we can expand this is given by (Bond & Efstathiou 1987)
the density field as
( ) 2
dn
p np dn = Ae-n p 2
p
p G (g, gnp), (31)
2p
d (x) = d (0) 1+ å zij xi xj, (23)2 ij where A is a normalizing factor.
For all peaks above a particular threshold ν, we can therefore
where i, j subscripts run over the two-dimensional x , x estimate the probability distribution of e as1 2
coordinate system. Here ζij=∇i∇jδ is a matrix composed of ¥
second derivatives of δ. We can diagonalize ζij and rewrite
P (e) = A0 ò p(np)P (e∣np)dnp, (32)n
14
The Astrophysical Journal, 858:87 (15pp), 2018 May 10 Appleby et al.
where A0 is a normalization factor. For a given ν threshold, we Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376
find the most likely ellipticity e as the expectation value Aurich, R., Janzer, H. S., Lustig, S., & Steiner, F. 2011, IJMPD, 20, 2253m Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, ApJ, 304, 15
1 2 Beisbart, C., Buchert, T., & Wagner, H. 2001, PhyA, 293, 592
em = ò eP (e)de. (33) Beisbart, C., Dahlke, R., Mecke, K., & Wagner, H. 2002, in Morphology of0 Condensed Matter, ed. K. Mecke & D. Stoyan (Berlin: Springer), 238
It remains to relate e to the axis ratio of an ellipse and then Beisbart, C., Valdarnini, R., & Buchert, T. 2001, A&A, 379, 412m Bond, J. R., & Efstathiou, G. 1987, MNRAS, 226, 655
calculate the Minkowski tensor for such a shape. Chingangbam, P., Park, C., Yogendran, K. P., & van de Weygaert, R. 2012,
The most probable ellipticity em is related to the most ApJ, 755, 122
probable axis ratio k º (a a ) as Chingangbam, P., Yogendran, K. P., Joby, P. K., et al. 2017, JCAP, 1712, 023m 1 2 m Codis, S., Dubois, Y., Pichon, C., Devriendt, J., & Slyz, A. 2014, in Proc. IAU
⎛ a ⎞ 1 - 2e 11, The Zeldovich Univere: Genesis and Growth of the Cosmic Web
km º
1 m
⎝⎜ ⎟⎠ = . (34) (Cambridge: Cambridge Univ. Press), 437a2 m 1 + 2em Codis, S., Gavazzi, R., Dubois, Y., et al. 2015b, MNRAS, 448, 3391
Codis, S., Pichon, C., & Pogosyan, D. 2015a, MNRAS, 452, 3369
For an ellipse with axes a1,2, we can calculate the Minkowski Doroshkevich, A. G. 1970, Ap, 6, 320
tensor W1,1 analytically (Schroder-Turk et al. 2010), Ganesan, V., & Chingangbam, P. 2017, JCAP, 6, 0232 Geller, M. J., Diaferio, A., & Kurtz, M. J. 2011, AJ, 142, 133
W1,1 = diag( f 1,1 (a , a ), f 1,1 (a , a )), (35) Geller, M. J., & Hwang, H. S. 2015, AN, 336, 4282 2 1 2 2 2 1 Gott, J. R., Dickinson, M., & Melott, A. L. 1986, ApJ, 306, 341
Gott, J. R., Weinberg, D. H., & Melott, A. L. 1987, ApJ, 319, 1
where Gott, J. R., III, Miller, J., Thuan, T. X., et al. 1989, ApJ, 340, 625
Hahn, O., Carollo, C. M., Porciani, C., & Dekel, A. 2007, MNRAS, 381, 41
2p 2
f 1,1 (a , a ) 1= a 2 2 cos fdfa ò . Hamilton, J. S. A., Gott, J. R., & Weinberg, D. 1986, ApJ, 309, 12 1 2 1 22 0 (a 2 2 2 2 3 2 Hong, S. E., Park, C., & Kim, J. 2016, ApJ, 823, 1031 - (a1 - a2 )cos f) Hug, D., Schneider, R., & Schuster, R. 2008, St. Petersburg Math. J, 19, 137
(36) Jiang, C. Y., Jing, Y. P., Faltenbacher, A., Lin, W. P., & Li, C. 2008, ApJ,
675, 1095
Note that the matrix is diagonal only in the coordinate basis Kerscher, M., Mecke, K., Schmalzing, J., et al. 2001, A&A, 373, 1
aligned with the axes of the ellipse. In this case, f 1,1 (a , a ) Kim, J., Park, C., Gott, R., & Dubinski, J. 2009, ApJ, 701, 15472 1 2 Kim, J., Park, C., L’Huillier, B., & Hong, S. E. 2015, JKAS, 48, 213
and f 1,12 (a2, a1) correspond to the eigenvalues of this matrix. Lee, J., & Pen, U.-L. 2002, ApJL, 567, L111
The ratio of these eigenvalues can be compared to b in the Lee, J., Springel, V., Pen, U.-L., & Lemson, G. 2008, MNRAS, 389, 1266tot
∣n∣ Mantz, H., Jacobs, K., & Mecke, K. 2008, JSMTE, 2008, P12015large limit. Matsubara, T. 1994, arXiv:astro-ph/9501076
For a Gaussian white-noise density eld, γ and hence W1,1fi 2 Matsubara, T. 1996, ApJ, 457, 13
are independent of the smoothing scale RG and power spectrum Matsubara, T. 2003, ApJ, 584, 1
amplitude, which are the only parameters in the analysis. When McMullen, P. 1997, Rend. Circ. Palermo, 50, 259
applying these statistics to a cosmological dark matter eld, Mecke, K. R., Buchert, T., & Wagner, H. 1994, A&A, 288, 697fi
1,1 Melott, A. L., Cohen, A. P., Hamilton, A. J. S., Gott, J. R., & Weinberg, D. H.W2 will depend on both RG and the cosmological parameters 1989, ApJ, 345, 618
Wmat, ns via the γ dependence of Equations (28) and (31). Melott, A. L., Weinberg, D. H., & Gott, J. R. 1988, ApJ, 328, 50
Park, C., & Gott, J. R. 1991, ApJ, 378, 457
Park, C., Gott, J. R., Melott, A. L., & Karachentsev, I. D. 1992, ApJ, 387, 1
ORCID iDs Park, C., Kim, J., & Gott, J. R. 2005, ApJ, 633, 1
Park, C., Pranav, P., Chingangbam, P., et al. 2013, JKAS, 46, 125
Changbom Park https://orcid.org/0000-0001-9521-6397 Patiri, S. G., Cuesta, A. J., Prada, F., Betancort-Rijo, J., & Klypin, A. 2006,
Sungwook E. Hong https://orcid.org/0000-0003-4923-8485 ApJL, 652, L75
Juhan Kim https://orcid.org/0000-0002-4391-2275 Paz, D., Stasyszyn, F., & Padilla, N. 2008, MNRAS, 389, 1127
Ryden, B. S., Melott, A. L., Craig, D. A., et al. 1989, ApJ, 340, 647
Schmalzing, J., & Buchert, T. 1997, ApJL, 482, L1
References Schmalzing, J., Kerscher, M., & Buchert, T. 1996, in Proc. of the International
School of Physics Course CXXXII, ed. S. Bonometto, J. R. Primack, &
Abell, P. A., Julius, A., Anderson, S. F., et al. 2009, arXiv:0912.0201 A. Provenzale (Oxford: IOS Press), 281
Adam, R., Ade, P. R., Aghanim, A., et al. 2016, A&A, 594, A9 Schroder-Turk, G., Kapfer, S., Breidenback, B., Beisbart, C., & Mecke, K.
Adler, R. 1981, The Geometry of Random Fields (New York: Wiley) 2010, JMic, 238, 57
Aghamousa, A., Aguilar, J., Ahlen, S., et al. 2016, arXiv:1611.00036 Schroder-Turk, G. E., Mickel, W., Kapfer, S. C., et al. 2013, NJPh, 15, 083028
Alesker, S. 1999, Geometriae Dedicata, 74, 241 Weinberg, D. H., Gott, J. R., & Melott, A. L. 1987, ApJ, 321, 2
15