Constraint of Void Bias on Primordial non-Gaussianity

TL;DR

This paper investigates how cosmic voids' large-scale bias, influenced by primordial non-Gaussianity, can improve constraints on PNG parameters, showing voids provide complementary and enhanced information beyond halos.

Contribution

It demonstrates that void bias exhibits scale-dependent corrections similar to halos and that including voids improves PNG constraints by a factor of two.

Findings

Void bias shows scale dependence similar to halos.

Large voids anti-correlate with dark matter density.

Including voids tightens PNG parameter constraints by a factor of two.

Abstract

We study the large-scale bias parameter of cosmic voids with primordial non-Gaussian (PNG) initial conditions of the local type. In this scenario, the dark matter halo bias exhibits a characteristic scale dependence on large scales, which has been recognized as one of the most promising probes of the local PNG. Using a suite of $N$-body simulations with Gaussian and non-Gaussian initial conditions, we find that the void bias features scale-dependent corrections on large scales, similar to its halo counterpart. We find excellent agreement between the numerical measurement of the PNG void bias and the general peak-background split prediction. Contrary to halos, large voids anti-correlate with the dark matter density field, and the large-scale Gaussian void bias ranges from positive to negative values depending on void size and redshift. Thus, the information in the clustering of voids can be complementary to that of the halos. Using the Fisher matrix formalism for multiple tracers, we demonstrate that including the scale-dependent bias information from voids, constraints on the PNG parameter $f_{\rm NL}$ can be tightened by a factor of two compared to the accessible information from halos alone, when the sampling density of tracers reaches $4 \times 10^{-3} \, h^3 \mathrm{Mpc}^{-3} $.

Equations

$$\Phi =\varphi +f_{\mathrm{NL}}({\varphi }^2-⟨{\varphi }^2⟩),$$

$$\mathcal{F}(\nu,\delta_{\rm v},\delta_{\rm c})=\sqrt{\frac{2}{\pi}}\exp\Big{(}-\frac{\nu^{2}}{2}\Big{)}\exp\Big{(}-\frac{|\delta_{\rm v}|}{\delta_{\rm c}}\frac{\mathcal{D}^{2}}{4\nu^{2}}-2\frac{\mathcal{D}^{4}}{\nu^{4}}\Big{)}.$$

$$n_{\mathrm{v}}=\frac{1}{V_{\mathrm{L}}}\frac{d\ln \nu }{d\ln R_{\mathrm{L}}}\nu \mathcal{F}(\nu ),$$

$${\delta }_{\mathrm{s}}(\mathbf{k})=\mathcal{M}(k){\varphi }_{\mathrm{s}}(\mathbf{k})(1+2f_{\mathrm{NL}}{\varphi }_{\mathrm{l}}).$$

$$\mathcal{M}(k)=\frac{2}{3}\frac{k^{2}T(k)D(z)}{{\Omega }_m{H}_0^2},$$

$$b_{\mathrm{v}}^{\mathrm{G}}=1+\frac{{\nu }^2-1}{{\delta }_{\mathrm{v}}}+\frac{{\delta }_{\mathrm{v}}\mathcal{D}}{4{\delta }_{\mathrm{c}}^2{\nu }^2}.$$

$$b_{\mathrm{v}}^{\mathrm{NG}}$$

$$b_{\rm v}^{\rm NG}(k)=\frac{3f_{\rm NL}\Omega_{\rm m}H_{0}^{2}}{k^{2}T(k)D(z)}\bigg{(}\nu^{2}-1-\frac{|\delta_{\rm v}|\mathcal{D}^{2}}{2\delta_{\rm c}\nu^{2}}-\frac{8\mathcal{D}^{4}}{\nu^{4}}\bigg{)}.$$

$$b_{\mathrm{v}}^{\mathrm{NG}}=\frac{3f_{\mathrm{NL}}{\Omega }_{\mathrm{m}}{H}_0^2}{k^{2}T(k)D(z)}{\delta }_{\mathrm{v}}(b_{\mathrm{v}}^{\mathrm{G}}-1).$$

$$⟨{\delta }_i(\mathbf{k}){\delta }_j({\mathbf{k}}^{\prime })⟩=P_{ij}(k){\delta }_{\mathrm{D}}(\mathbf{k}+{\mathbf{k}}^{\prime }),$$

$$\displaystyle F_{f_{\rm NL}f_{\rm NL}}=V\int\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{2}\mathrm{Tr}\Big{(}\Sigma^{-1}\frac{\partial\Sigma}{\partial f_{\rm NL}}\Sigma^{-1}\frac{\partial\Sigma}{\partial f_{\rm NL}}\Big{)},$$

Full text

I Introduction

Primordial non-Gaussianity (PNG) offers an important probe into the physics of inflation Bartolo et al. (2004); Liguori et al. (2010); Chen (2010), as it allows to constrain the production mechanism of primordial perturbations which seed the structures we observe in the Universe today. Furthermore, PNG can be directly related to primordial interactions taking place at energies as high as $10^{14}$ GeV, providing a unique window into the particle content of inflation Chen (2010); Arkani-Hamed and Maldacena (2015). In the local PNG model, the Bardeen potential $\Phi$ is given by Gangui et al. (1994); Verde et al. (2000); Komatsu and Spergel (2001)

[TABLE]

where $\phi$ is a Gaussian potential field and $f_{\rm NL}$ parametrizes the strength of non-Gaussianity. Using maps of the CMB, the Planck collaboration derived a stringent constraint of $f_{\rm NL}=0.8\pm 5.0$ Ade et al. (2016).

Besides the CMB, the large-scale structure (LSS) is another frontier in constraining PNG. It can be detected using the galaxy bispectrum Verde et al. (2000); Scoccimarro et al. (2004); Sefusatti and Komatsu (2007). However, to extract the feeble PNG signal from it, dominant contributions from late-time non-Gaussianities due to dark matter nonlinearities, galaxy biasing, and redshift-space distortions must be modeled well. Recovering information from the bispectrum is also hampered by its large covariance Chan and Blot (2017). On the other hand, it has been discovered that the halo bias exhibits a strong scale dependence on large scales in the local PNG scenario Dalal et al. (2008). In contrast, for Gaussian initial conditions the large-scale halo bias remains scale-independent. Scale-dependent halo bias from PNG has been extensively investigated with numerical simulations Dalal et al. (2008); Desjacques et al. (2009); Pillepich et al. (2010); Grossi et al. (2009); Giannantonio and Porciani (2010); Hamaus et al. (2011); Scoccimarro et al. (2012); Biagetti et al. (2017). Since the inflationary consistency relation implies that no scale-dependent halo bias is generated in single field inflation Maldacena (2003); Tanaka and Urakawa (2011); Creminelli et al. (2011); Baldauf et al. (2011); Dai et al. (2015); de Putter et al. (2015), its detection in the bias of tracers on large scales (modulo projection effects Yoo (2010); Bonvin and Durrer (2011); Challinor and Lewis (2011); Jeong et al. (2012)) offers a means to rule out single-field inflation. This unique feature has been applied to constrain local PNG using galaxy survey data Slosar et al. (2008); Afshordi and Tolley (2008); Xia et al. (2011); Ross et al. (2013); Giannantonio et al. (2014); Leistedt et al. (2014); Ho et al. (2015). Although contaminations from late-time non-Gaussianity are relatively mild in the linear regime, it has been realized that the low-$k$ part of the power spectrum is susceptible to observational systematics, such as stellar contamination Ross et al. (2013); Pullen and Hirata (2013); Huterer et al. (2013). After carefully eliminating systematics, the current bound on $f_{\rm NL}$ is $5\pm 21$ (cross correlation between various data sets Giannantonio et al. (2014)) and $-39<f_{\rm NL}<23$ (quasars Leistedt et al. (2014)). The constraints from future surveys are expected to tighten by one to two orders of magnitude Hamaus et al. (2011); Dore et al. (2014); Alonso and Ferreira (2015); de Putter and Dore (2017); Raccanelli et al. (2015); Yamauchi et al. (2014); Ferraro and Smith (2015); Ferramacho et al. (2014). This can be achieved by combining multiple tracers of the LSS, which allows one to cancel out the dominant cosmic variance contribution on large scales McDonald and Seljak (2009); Seljak (2009); Hamaus et al. (2011, 2012); Abramo and Leonard (2013).

Almost all studies on PNG focus on tracers with positive bias parameters (except Sekiguchi and Yokoyama (2012); Uhlemann et al. (2018); Castorina et al. (2018)), with galaxies as the prime example. Voids are distinct from halos, because their large-scale bias ranges from positive to negative values as the void size increases Hamaus et al. (2014); Chan et al. (2014). In recent years, various clustering properties of voids have been measured using galaxy samples: redshift-space distortions around voids Paz et al. (2013); Hamaus et al. (2015); Cai et al. (2016); Hawken et al. (2017); Hamaus et al. (2016); Achitouv et al. (2017); Hamaus et al. (2017), the configuration-space void bias Clampitt et al. (2016), the tracer bias around voids Pollina et al. (2018), and the Baryon Acoustic Oscillations from voids Kitaura et al. (2016). Because the clustering of voids enables us to probe a range of bias that is not accessible to halos, in this work we investigate void bias in the presence of local PNG and its potential constraining power on $f_{\rm NL}$. We consider voids defined both in dark matter, as well as halo density fields in real space.

II Theory

We derive the void bias in PNG using the peak-background split formalism Kaiser (1984); Mo and White (1996); Sheth and Tormen (1999); Desjacques et al. (2018). To do so we consider the response of the void size distribution to long wavelength perturbations. The void size distribution can be modeled as a first-crossing distribution problem in the excursion set formalism Bond et al. (1991). To be concrete, let us consider a simple first crossing distribution $\mathcal{F}$ for voids Sheth and van de Weygaert (2004):

[TABLE]

The peak significance of a void is defined as $\nu\equiv|\delta_{\rm v}|/\sigma_{R_{\rm L}}$, where $\delta_{\rm v}$ is void formation threshold (in the spherical collapse model the linearly extrapolated value is -2.72 Blumenthal et al. (1992)) and $\sigma_{R_{\rm L}}$ is the RMS of the density field evaluated using a top-hat filter of size $R_{\rm L}$, the Lagrangian size of the void. $\mathcal{D}\equiv|\delta_{\rm v}|/(\delta_{\rm c}+|\delta_{\rm v}|)$ denotes the so-called void-and-cloud parameter with $\delta_{\rm c}$ being the halo collapse threshold (1.68 in the spherical collapse model Gunn and Gott (1972)). $\mathcal{F}$ (modulo a Jacobian) is the fraction of Lagrangian space volume characterized by $R_{\rm L}$ that will turn into voids of size $R_{\rm v}$ in Eulerian space. We can map the Eulerian void size to the Lagrangian one using the spherical collapse relation in Bernardeau (1994). The void size distribution can then be written as

[TABLE]

where the Lagrangian void volume is $V_{\rm L}=(4\pi/3)R_{\rm L}^{3}$.

To derive the effect of a long density fluctuation mode on small-scale ones in the local PNG case, we first split the Gaussian potential $\phi$ into long- and short-wavelength perturbations $\phi=\phi_{\rm l}+\phi_{\rm s}$. The small-scale Bardeen potential $\Phi_{\rm s}$ becomes $\Phi_{\rm s}\approx\phi_{\rm s}+2f_{\rm NL}\phi_{\rm l}\phi_{\rm s}$. The small-scale overdensity is obtained via the Poisson equation

[TABLE]

The factor $\mathcal{M}$ reads

[TABLE]

where $\Omega_{\rm m}$ and $H_{0}$ are the matter density and Hubble parameter at present time, $T(k)$ is the transfer function, and $D(z)$ the growth factor normalized to the scale factor in the matter-dominated era. By considering the response of $n_{\rm v}$ to a long-wavelength perturbation, we obtain the linear Gaussian bias from Eq. (2) as Sheth and van de Weygaert (2004); Chan et al. (2014)

[TABLE]

In addition to $b_{\rm v}^{\rm G}$, there is an extra contribution $b_{\rm v}^{\rm NG}$ in the presence of local PNG. Following Slosar et al. (2008), in the local PNG model a long-wavelength perturbation effectively rescales the amplitude of the small-scale fluctuations [Eq. (4)]. The amplitude of linear fluctuations is parametrized by $\sigma_{8}$, defined as the RMS of the density field at scales of $8\,\mathrm{Mpc}\,h^{-1}$. We can write the PNG bias as a response of $n_{\rm v}$ to the local value of $\sigma_{8}$ as Slosar et al. (2008)

[TABLE]

For the case of halos one typically assumes the universality of the halo mass function and replaces the local $\sigma_{8}$ amplitude by $\sigma_{R_{\rm L}}$. Following the same approach for voids, with Eq. (2) we obtain

[TABLE]

This is similar to the PNG void bias derived in the high-peak limit derived in Sekiguchi and Yokoyama (2012):

[TABLE]

It is the analog of the well-known PNG halo bias Dalal et al. (2008); Matarrese and Verde (2008); Slosar et al. (2008) (with $\delta_{\rm v}$ replaced by $\delta_{\rm c}$) and agrees with Eq. (8) in the high-peak limit, which is valid for large void sizes.

III Numerical results

The $N$-body simulation used in this work contains $1536^{3}$ particles in a box of 2000 $\,\mathrm{Mpc}\,h^{-1}$ side length. There are three sets of different initial conditions: $f_{\rm NL}=0$, $250$, and $-250$, with each eight realizations. The cosmological parameters are $\Omega_{\rm m}=0.3$, $\Omega_{\Lambda}=0.7$, $n_{s}=0.967$, and $\sigma_{8}=0.85$. To compute $\partial\ln n_{\rm v}/\partial\ln\sigma_{8}$ numerically, we use two sets of Gaussian simulations with $\sigma_{8}=0.83$ and $0.87$. They share the same cosmological parameters as the other Gaussian simulations except $\sigma_{8}$ and there are two realizations each. The initial particle displacements are implemented using 2LPTic Crocce et al. (2006); Scoccimarro et al. (2012) at $z=99$, and then evolved with Gadget2 Springel (2005). Halos are identified using the halo finder Rockstar Behroozi et al. (2013). For more details, we refer the readers to Biagetti et al. (2017). Void catalogs are extracted using the void finder VIDE Sutter et al. (2015), which is based on ZOBOV Neyrinck (2008) using a watershed algorithm Platen et al. (2007).

As a consequence of the typically large extent of voids, their number density is generally low and shot noise can be substantial. Furthermore, exclusion effects are significant and cause strong scale dependence on relatively large scales in the void auto-power spectrum Hamaus et al. (2014); Chan et al. (2014). Thus, the scale-dependent bias from PNG is most apparent when voids are cross-correlated with other tracer species. The cross-power spectrum between the species $i$ and $j$ is given by

[TABLE]

where $\delta_{i}$ and $\delta_{j}$ are their overdensities, and $\delta_{\rm D}$ is the Dirac delta function. On large scales it can be expressed as $P_{ij}=b_{i}b_{j}P_{\rm mm}+\mathcal{E}_{ij}$, where $b_{i}$ are the linear bias parameters (Gaussian or non-Gaussian), $P_{\rm mm}$ is the matter power spectrum, and $\mathcal{E}_{ij}$ the shot noise matrix Hamaus et al. (2010).

Fig. 1 shows the cross bias between matter and voids, $b_{\rm mv}\equiv P_{\rm mv}/P_{\rm mm}$, both in the Gaussian and the PNG simulations ($f_{\rm NL}=250$) for a range of void sizes at $z=0$. These voids are identified in the matter density field in real space with a sampling density of $0.005\,(\,\mathrm{Mpc}^{-1}\,h)^{3}$ and binned into different effective radii $R_{\rm v}$ of bin-width $10\,\mathrm{Mpc}\,h^{-1}$. The Gaussian void bias becomes independent of scale at low $k$, so we fit it with a constant up to $k<0.03\,\mathrm{Mpc}^{-1}\,h$. It decreases from $b_{\rm v}^{\rm G}\sim 1$ to 0 when the void size reaches $R_{\rm v}\sim 25\,\mathrm{Mpc}\,h^{-1}$, and goes negative when the void size further increases. The precise values of the bias parameters depend on the sampling density (c.f. Chan et al. (2014)). As shown in Chan et al. (2014), the Gaussian bias from Eq. (6) can qualitatively describe the trend in the simulations provided that $\delta_{\rm v}$ is chosen less negative than its spherical collapse value of $-2.72$., e.g. $\delta_{\rm v}\simeq-1.0$ in Chan et al. (2014) (see also Pisani et al. (2015)).

Furthermore, the PNG cross bias exhibits strong scale-dependence in the low-$k$ regime. We over-plot the prediction from Eq. (7) and find an excellent agreement with the numerical measurement. The derivative $\partial\ln n_{\rm v}/\partial\ln\sigma_{8}$ is computed by finite differencing $n_{\rm v}$ measured in the Gaussian simulations with $\sigma_{8}=0.83$ and $0.87$. From Fig. 1 it is evident that the sign of the additional PNG bias changes from negative, when $b_{\rm v}^{\rm G}\gtrsim 0$, to positive, when $b_{\rm v}^{\rm G}\lesssim 0$. Contrary to Castorina et al. (2018), where zero-bias tracers are suggested to be very sensitive to PNG, our numerical results suggest the PNG signal to almost vanish for this type of tracer.

We go on to test Eq. (9), which is expected to be a good approximation for large voids. In Fig. 2 we present the results for $\delta_{\rm v}$ = $-2.7$ and $-1.0$, which both fail to accurately describe the numerical results. Poor agreement is also found using Eq. (8). Eq. (9) implies that the PNG signal should flip sign at $b_{\rm v}^{\rm G}=1$ instead of $b_{\rm v}^{\rm G}\sim 0$. If voids tend to remain at their Lagrangian position, then $b_{\rm v}^{\rm G}$ is simply equal to the Lagrangian bias and the sign indeed flips at $b_{\rm v}^{\rm G}=0$, but recent work Massara and Sheth (2018) shows that voids do move along with the LSS. However, to arrive at Eq. (8) [or (9)], we implicitly replace $\sigma_{\rm 8}$ by $\sigma_{R_{\rm L}}$ in the definition of $\nu$. These $\sigma$’s are evaluated at different scales. While interchanging them appears to be fine for halos, it leads to erroneous predictions for voids. Hence, for voids one should apply the more general peak-background split result of Eq. (7). A further advantage of Eq. (7) is that it is independent of the definition criterion for voids, because it only depends on the measured void size distribution.

In order to circumvent the difficulty in observing $P_{\rm mv}$ in galaxy surveys, we consider the cross-power spectrum between halos and voids, $P_{\rm hv}$. We define $b_{\rm hv}\equiv P_{\rm hv}/P_{\rm mm}$. On large scales, where linear bias is valid, it is natural to expect $b_{\rm hv}\approx P_{\rm mh}P_{\rm mv}/P^{2}_{\rm mm}$. In Fig. 3 we plot the numerical measurements of $b_{\rm hv}$ between a halo bin of mean mass $M_{\rm h}=1.1\times 10^{13}\,M_{\odot}{h}^{-1}$ and two different void-size bins, $R_{\rm v}=$ 15 and 35 $\,\mathrm{Mpc}\,h^{-1}$. These void samples are the same as the ones used in Fig. 1. There is significant scale dependence in the PNG case relative to the Gaussian one on large scales. We also show $P_{\rm mh}P_{\rm mv}/P^{2}_{\rm mm}$ using the numerical power spectra. Indeed $b_{\rm hv}$ agrees well with $b_{\rm mh}b_{\rm mv}$, but on smaller scales where nonlinearity and nonlinear biasing kick in, we expect deviations from this simple relation (although it is not apparent for this halo mass bin).

IV Fisher forecast

In this section we study how much additional information on $f_{\rm NL}$ can be gained by including voids in a multitracer analysis using the Fisher matrix formalism (see Heavens (2009); Tegmark et al. (1997); Albrecht et al. (2006) for a review). For simplicity, we consider $f_{\rm NL}$ as the only free parameter, as degeneracies with other cosmological parameters are expected to be of minor importance when multiple tracers of the same underlying density field are considered McDonald and Seljak (2009); Seljak (2009); Hamaus et al. (2011). Assuming the Fourier modes of the overdensities of halos and voids are Gaussian distributed, the Fisher matrix of the fields reads

[TABLE]

where $V$ is the survey volume and $\Sigma(k)$ the covariance matrix of the multitracer density field Hamaus et al. (2011). The elements $\Sigma_{ij}(k)$ are simply given by all possible auto- and cross-power spectra $P_{ij}(k)$ from Eq. (10). We consider a survey volume of 8$\,\mathrm{Gpc}^{3}h^{-3}$ at $z=1$, and compare constraints on $f_{\rm NL}$ from a multitracer survey of halos, voids, and a combination of the two. To be more realistic, voids are constructed from the distribution of halos in real space rather than the dark matter. In each case, the multitracer technique takes advantage of sampling variance cancelation in the primordial modes, which are fairly uncorrelated on large scales Seljak (2009); Hamaus et al. (2011).

In calculating $F_{f_{\rm NL}f_{\rm NL}}$, we directly use the numerical power spectra, including their shot noise contribution. The $f_{\rm NL}$-response derivative $\partial\Sigma/\partial f_{\rm NL}$ is obtained via finite differencing of the $f_{\rm NL}=250$ and $f_{\rm NL}=-250$ power spectrum measurements. We remove the scale-independent contribution due to differences in shot noise between the catalogs. To increase the accessible range in number densities, we use halos containing at least five particles. As we are primarily interested in the large-scale modes of the density field, the fact that these halos are not well resolved is of minor importance. Using higher-resolution simulations we have verified that our forecast is robust up to an uncertainty of at most $10\%$111We used the Carmen simulation from the LasDamas project (http://lss.phy.vanderbilt.edu/lasdamas/overview.html) to perform a partial resolution study. The LasDamas simulations are Gaussian, so we rely on the $f_{\rm NL}$-response measured at lower resolution, but obtain more accurate power spectra.. We divide the full halo sample (minimum 5 particles) into five mass bins, with mean halo masses $M_{\rm h}=1.33\times 10^{12}$, $4.86\times 10^{12}$, $1.10\times 10^{13}$, $3.34\times 10^{13}$, and $1.44\times 10^{14}\,\,M_{\odot}{h}^{-1}$. Voids are extracted from the halo samples with various mass thresholds (minimum 5, 10, 20, or 50 particles). The resulting void samples are further divided into three bins of void size $R_{\rm v}=$ [0, 20], [20, 40], and [40, 80] $\,\mathrm{Mpc}\,h^{-1}$. A detailed modeling of the PNG signals from these void samples is beyond the scope of this paper and will be investigated elsewhere.

Fig. 4 presents the constraints on $f_{\rm NL}$ via $\sigma_{f_{\rm NL}}\equiv 1/\sqrt{F_{f_{\rm NL}f_{\rm NL}}}$, using all Fourier modes up to $k_{\rm max}=0.08\,\mathrm{Mpc}^{-1}\,h$. The increasing number densities correspond to halo samples with at minimum 5, 10, 20, or 50 particles, respectively. The halo constraint is obtained using the available mass bins depending on the minimum particle threshold222Because the lowest halo mass bin corresponds to 5-20 particles, the halo constraint for minimum 10-particles is not shown. For the joint constraint in this case, we use all the halo mass bins. . The constraints from voids alone are weak compared to those from halos due to their high shot noise, saturating for $n_{\rm h}\gtrsim 7\times 10^{-4}\,\,\mathrm{Mpc}^{-3}\,h^{3}$. When the number density of halos increases, the distribution of void sizes shifts to smaller voids. As we are keeping the bins of void sizes fixed, this means we are losing the largest voids that contribute the strongest signal on $f_{\rm NL}$, due to their very negative bias amplitudes. However, when voids and halos are combined in a multitracer analysis, the joint constraints improve appreciably thanks to the large cross-correlations and the low shot noise between them. In particular, when the halo number density reaches $4\times 10^{-3}\,\,\mathrm{Mpc}^{-3}\,h^{3}$, the error on $f_{\rm NL}$ is reduced by almost a factor of 2 compared to the halo case. This result is in good agreement with the analytical Fisher forecast based on Eq. (29) in Hamaus et al. (2011), shown as the green dotted line in Fig. 4. In that calculation, we used the $f_{\rm NL}$ response from the cross-power spectra of each tracer with the matter field and assumed Poisson shot noise.

V Conclusion

Using a suite of $N$-body simulations, we have demonstrated that in the local PNG model, voids exhibit a scale-dependent bias on large scales, just like halos. Although the standard calculation that is analogous to the halo case [Eq. (8) or (9)] fails to describe the simulation results for voids, the general peak-background split prediction [Eq. (7)] yields an excellent agreement. Furthermore, based on the Fisher matrix formalism for multiple tracers we have shown that by combining the clustering information from voids and halos, constraints on $f_{\rm NL}$ can substantially be tightened, as long as the number density of tracers is sufficiently high. Our simplistic analysis using a volume of $8\,\mathrm{Gpc}^{3}h^{-3}$ and tracer densities up to $4\times 10^{-3}\,\,\mathrm{Mpc}^{-3}\,h^{3}$ already renders $f_{\rm NL}$ constraints of $\mathcal{O}$(a few) achievable. Although we only show the results in real space, the effects from redshift-space distortions on the clustering of voids on large, linear scales is well studied Paz et al. (2013); Hamaus et al. (2015); Cai et al. (2016); Hawken et al. (2017); Hamaus et al. (2016); Achitouv et al. (2017); Hamaus et al. (2017); Chuang et al. (2017). Optimizing the binning strategy in constructing multiple tracers from halos and voids will most likely yield further gains Hamaus et al. (2010, 2011). Future surveys, such as Euclid, will have access to even larger volumes and higher densities of tracers, opening up the possibility to significantly improve upon current CMB constraints on PNG with the help of cosmic voids. The latter are contained in the survey data anyway and hence provide additional information at no cost.

Acknowledgements

We thank Vincent Desjacques for sharing his numerical simulations with us for the initial phase of the study, and Neal Dalal for fruitful discussions. KCC acknowledges the support from the National Science Foundation of China under the grant 11873102. NH acknowledges support from the DFG cluster of excellence “Origins” and the Trans-Regional Collaborative Research Center TRR 33 “The Dark Universe” of the DFG. NH is grateful to the Institut Henry Poincaré for hospitality during the “CosmoMath” trimester in the fall 2018. MB acknowledges support from the Nederlands Organisation for Scientific Research (NWO) under the VENI grant 016.Veni.192.210.