{'fraction_non_alphanumeric': 0.04659508026077905, 'fraction_numerical': 0.016463047170328166, 'mean_word_length': 4.159270774448841, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Local stochastic non-Gaussianity and N-body simulations

Large-scale clustering of highly biased tracers of large-scale structure has emerged as one of the best observational probes of primordial non-Gaussianity of the local type (i.e. f_{NL}^{local}). This type of non-Gaussianity can be generated in multifield models of inflation such as the curvaton model. Recently, Tseliakhovich, Hirata, and Slosar showed that the clustering statistics depend qualitatively on the ratio of inflaton to curvaton power \xi after reheating, a free parameter of the model. If \xi is significantly different from zero, so that the inflaton makes a non-negligible contribution to the primordial adiabatic curvature, then the peak-background split ansatz predicts that the halo bias will be stochastic on large scales. In this paper, we test this prediction in N-body simulations. We find that large-scale stochasticity is generated, in qualitative agreement with the prediction, but that the level of stochasticity is overpredicted by ~30%. Other predictions, such as \xi independence of the halo bias, are confirmed by the simulations. Surprisingly, even in the Gaussian case we do not find that halo model predictions for stochasticity agree consistently with simulations, suggesting that semi-analytic modeling of stochasticity is generally more difficult than modeling halo bias.


Introduction
One of the most exciting prospects for cosmology in the near future is the ability to constrain the physics of inflation [1,2,3,4,5,6,7], thus probing energy scales which are far beyond the reach of accelerator experiments. The simplest choice of initial conditions, namely stochastic initial fluctuations which are adiabatic, scalar, Gaussian, and scale-invariant, has been ruled out at the ≈3σ level. Current observations are consistent with either a power-law initial power spectrum which is redder than scale invariant (n s − 1 ≈ −0.04), and marginally consistent with initial conditions which are scale-invariant but contain contributions from tensor modes (r ≈ 0.2) [8]. The next few years will bring a wealth of new data which will sharpen this picture considerably.
Primordial non-Gaussianity has emerged as a particularly powerful probe of inflation due to the ability to rule out large qualitative classes of models. For example, there is a theorem [9,10,11] which states that in all models of single field inflation whose power spectrum is nearly scale invariant, the 3-point function ζ(k 1 )ζ(k 2 )ζ(k 3 ) is observationally indistinguishable from zero in "squeezed" triangles (i.e. min(k i ) max(k i )). However, a detectably large squeezed 3-point function is naturally generated in other models, such as the ekpyrotic scenario [12,13,14,15]. Observational constraints on Gaussianity to date have mainly focused on the 3-point function and have parameterized the deviation from zero by three parameters f local N L , f equil N L , f orthog N L [16,17,18,19,20,21,22]. The current WMAP constraints from [8]  In this paper, we will focus on local non-Gaussianity and use the notation f N L = f local N L throughout. In this case, the initial curvature fluctuation is of the form ζ(x) = ζ G (x) + 3 5 f N L (ζ G (x) 2 − ζ 2 G (x) ), where ζ G is a Gaussian field [23,24,25]. Recently, this type of non-Gaussianity has been studied extensively in the context of large-scale structure, beginning with a pioneering paper by Dalal et al. [26], which showed that scale-dependent halo bias is generated on large scales. Subsequently, this prediction has been confirmed and extended, in both analytical and simulation-based studies [27,28,29,30,31,32,33,34,35,36,37,38]. In particular, in [33] an improved expression for the f N L dependence of the bias was introduced (the improved expression agrees with the original expression from [26] in the limit k → 0). In [28], the constraint f local N L = 20 ± 25 was obtained from observations of large-scale halo clustering in SDSS (the measurement is obtained from a variety of tracer objects, but the statistical weight is dominated by the high-z photometric quasar sample). One qualitative finding which will be particularly relevant for this paper is that the scale-dependent bias is non-stochastic, in the sense that the correlation coefficient between halos in different mass bins is equal to one, after shot noise has been subtracted.
The curvaton model is a two-field model of inflation in which the source of initial curvature fluctuations is not the inflaton, but a second field σ whose contribution to the energy density during inflation is subdominant [39,40,41,42,43,44]. Most studies of the curvaton model have only considered the case where the curvaton contribution to the primordial curvature fluctuation ζ is much larger than the inflaton contribution. In this case, the curvaton field can give rise to non-Gaussianity of the local type at a detectable level (f N L is essentially a free parameter of the model).
Recently, Tseliakhovich, Hirata, and Slosar considered the more general case in which the ratio ξ of inflaton and curvaton contributions to the primordial curvature fluctuation ζ can be significantly different from zero [45] (see also [46]). Applying the same theoretical arguments which predict scaledependent, non-stochastic halo bias for ξ = 0, the authors argue that scale-dependent stochastic halo bias should be present in the more general case where f N L and ξ are both nonzero. A two-component hybrid model with similar observational signatures was studied in [47]; in this case ξ is typically of order one if the non-Gaussianity is large enough to be detectable.
The purpose of this paper is to analyze N -body simulations whose initial conditions contain curvaton and inflaton contributions, and study the dependence of halo clustering and halo stochasticity on the parameters {f N L , ξ} of the model. Although our main interest is the case ξ = 0, we also present results for the curvaton model with ξ = 0 as a baseline for comparison.

Initial conditions
In this subsection we review the curvaton model, in the same generality as [45].
The curvaton is assumed to decay before dark matter freezout, so that no dark matter isocurvature mode is generated, and the adiabatic curvature fluctuation ζ is a sum of inflaton and curvaton contributions: We assume that the fields ζ i and ζ c are uncorrelated and that their power spectra P ζ i , P ζc , are proportional, so that we can define a parameter ξ = (P ζ i /P ζc ) 1/2 which is independent of scale. The power spectra of ζ i , ζ c are thus related to the power spectrum P ζ of the total curvature fluctuation by: The power spectrum P ζ is taken to be of power-law form (k 3 /2π 2 )P ζ (k) = ∆ 2 ζ (k/k piv ) ns−1 with parameters ∆ ζ , k piv given in §1.
We assume that ζ i is a Gaussian field, but ζ c is a non-Gaussian field of "local type", i.e.
where ζ c,G is a Gaussian field. Non-Gaussianity of local type is generated if the curvaton potential V (σ) is assumed quadratic in σ. Throughout this paper, we will take {f N L , ξ} to be the parameters of the curvaton model.
To get some intuition for this parameterization, it is useful to note that the power spectrum, bispectrum, and connected higher-point functions of the initial curvature fluctuation ζ depend on f N L and ξ as follows: To lowest order in f N L , the 3-point function is proportional to f N L , with no ξ dependence 1 . This makes it easy to interpret the CMB bispectrum constraint from [8] as a constraint f N L = 32 ± 21 (1σ error), with no constraint on ξ. The CMB trispectrum constraint from [50] can similarly be interpreted as a constraint τ N L = (1.35 ± 0.98) × 10 4 on the combination of parameters τ N L = ( 6 5 f N L ) 2 (1+ξ 2 ), with the caveat that the bispectrum and trispectrum estimators are not statistically independent and there are subtleties in combining the assoicated parameter constraints [51].
The effect of nonzero ξ is to boost the amplitude of the N -point correlation functions (where N ≥ 4) relative to the amplitude of the 3-point function. One interesting consequence is that the 4-point function can be made large while keeping the 3-point function within observational limits on f N L . (In fact, halo stochasticity, or "boosting" the amplitude of the halo-halo power spectrum P hh relative to the amplitude of the matter-halo power spectrum P mh , can be viewed as a formal consequence of boosting the primordial 4-point function relative to the 3-point function.)

Halo clustering and the peak-background split
The peak-background split formalism is a heuristic argument for predicting correlation functions in which one or more scales is large compared to the scales which are relevant for spherical collapse [52,53,54,55,56]. This formalism can be applied to study non-Gaussian halo clustering in the curvaton model [28,45]. We will review this calculation here, and extend it by including 1-halo terms which will be relevant for the stochasticity results to be presented later.
Let us write the inflaton contribution ζ i to the initial curvature fluctuation as a sum of longwavelength and short-wavelength pieces: ζ i = (ζ i,l + ζ i,s ). We analogously write the Gaussian field ζ c,G as a sum (ζ c,l + ζ c,s ). The total initial curvature fluctuation is then given by 1 We have defined fNL in Eq. (4) with the extra factor of (1+ ξ 2 ) so that the 3-point function will have this property.
and we have assumed the long and short wavelength parts of ζ c are uncorrelated. Let us make the approximation that the terms in the third line of (9) are negligible. (The ζ 2 c,l term will not be important for our purposes since we will only use the small-scale component of Eq. (9); the main effect of the ζ 2 c,s term is to change the constant part of the halo bias, which is a free parameter anyway.) The first line is the long-wavelength part of ζ, which is unchanged from the Gaussian case (i.e. it does not depend on f N L ). The second line is the short-wavelength part; we find that the effect of the non-Gaussianity is to modulate the small-scale curvaton mode ζ c,s by a factor (1 + 6 5 f N L (1 + ξ 2 ) 2 ζ c,l ) which depends on the long-wavelength curvaton mode ζ c,l . In the peak-background split picture, we interpret Eq. (9) as saying that the small-scale matter power spectrum is no longer spatially constant in a non-Gaussian cosmology, but rather a local quantity which varies with position. If we consider a large box at position x, then the average small-scale power spectrum in the box is given by Following notation from [28], we will parameterize the amplitude of the small-scale power spectrum by σ 8 , the RMS of the linear density field at z = 0 with 8h −1 Mpc tophat smoothing, and rephrase Eq. (10) by writing σ 8 as a function of position x: Now let us ask how the number density n h of halos varies on large scales in the peak-background split picture. The density of halos n h (x) in a large box at position x will differ from the mean densitȳ n h for two reasons: first, because the local matter density ρ m (1 + δ l (x)) in the box differs from the mean ρ m , and second because the local value of σ 8 differs from the mean via Eq. (11). Combining these effects we can write: The (1 + δ l ) prefactor comes from converting Lagrangian to Eulerian space. The second term is proportional to the (scale-independent) derivative (∂ logn h /∂δ l ) of the mass function with respect to the background density, i.e. the Lagrangian halo bias. These two terms are present in a Gaussian cosmology and represent the usual halo bias which is constant on large scales. The third term represents the effects of primordial non-Gaussianity, which gives an extra scale-dependent contribution.
Taking the Fourier transform of Eq. (12) and dropping second-order terms, we get where δ h (k) = n h (k)/n h is the fractional halo overdensity, and in the last line we have defined where The quantity α(k, z) relates the matter overdensity δ m (k, z) to the initial curvature ζ(k) in linear perturbation theory: δ m (k, z) = 3 5 α(k, z)ζ(k). Note that we have defined inflaton and curvaton contributions to the matter density by δ i (k, z) = 3 5 α(k, z)ζ i (k) and δ c (k, z) = 3 5 α(k, z)ζ c (k), even though strictly speaking, Poisson's equation applies only to the sum of the two fields.
The peak-background split expression (13) for δ h (k) applies on scales which are large compared to scales relevant for spherical collapse. It is also incomplete, in the sense that it treats the halo overdensity as a continuous field, and ignores stochastic variations due to random halo locations. If we compute matter-halo and halo-halo power spectra using this expression, then we get: where the primes refer to different halo masses and P lin (k) = 9 25 α 2 (k)P ζ (k) is the linear theory matter power spectrum. We have included the superscript "2H" because these expressions omit 1-halo terms. Using standard machinery from the halo model [57,58,59,60], it is straightforward to calculate 1-halo contributions to these power spectra. In non-overlapping mass bins, let n i be the number density of halos in the i-th bin, and let f i denote the total fraction (by mass) of dark matter in halos in mass bin i. Then we get: 2 where we have defined Although the 1-halo terms are generally smaller than the 2-halo terms on the angular scales we will study in this paper (k ≤ 0.04 h Mpc −1 ), they are the leading souce of stochasticity (aside from shot noise) predicted by the halo model in the Gaussian case. In the non-Gaussian case, the term N G ) represents extra stochasticity on large scales for nonzero f N L and ξ > 0. This source of 2 These expressions neglect convolution by the halo density profiles, but this can be neglected for purposes of this paper, where we only study clustering on large scales (k ∼ < 0.04 h Mpc −1 ). stochasticity can be understood intuitively: the non-Gaussian part of the bias traces the curvaton field δ c on large scales, and if ξ > 0, this field is not 100% correlated to the matter overdensity δ.
There is a useful simplification to the above expression if we assume a "universal" halo mass function of the form dn/dM = (ρ m /M )f (ν)(dν/dM ). Here, ν(M, z) = (δ 2 c /σ 2 (M, z)), where δ c = 3(12π) 2/3 /20 ≈ 1.69 is the threshhold for spherical collapse, and σ 2 (M, z) is the variance of the linear matter overdensity after tophat smoothing on the scale corresponding to halo mass M . For a universal mass function, the relation holds [28], so that b N G and b G are related by Mass functions obtained from simulations are roughly universal [61,55,62,63], and we will generally assume that Eq. (24) holds throughout the paper.

N -body simulations
To study halo clustering in the curvaton model, we performed collisionless N -body simulations using the GADGET-2 TreePM code [64]. Simulations were done using periodic box size R box = 1600 h −1 Mpc, particle count N p = 1024 3 , and force softening length R s = 0.05(R box /N 1/3 p ). With these parameters and the fiducial cosmology from §1, the particle mass is m p = 2.92 × 10 11 h −1 M . Results in this paper were obtained from two simulations with f N L = 0, and two simulations for each choice of f N L ∈ {±250, ±500} and ξ ∈ {0, 1} (for a total of 18 simulations). For given curvaton model parameters f N L , ξ, we simulate initial conditions as follows. First, we simulate Gaussian fields ζ i and ζ c,G in Fourier space with power spectra given by Eqs. (2), (3). We then compute the non-Gaussian curvaton field in real space by ζ c = ζ c,G + 3 . (When generating the initial conditions, all Fourier transforms are computed on a grid with N 3 p elements.) We apply the transfer function T (k) to the total curvature fluctuation ζ = ζ i + ζ c to obtain the Newtonian potential Φ(k) at the initial redshift z ini = 100 of the simulations. Finally, we obtain initial particle positions using the Zeldovich approximation [65]. (At z ini = 100, transient effects due to use of this approximation should be negligible [66]. ) We group particles into halos using an MPI parallelized implementation of the friends-of-friends (FOF) algorithm [67] with link length L FOF = 0.2R box N −1/3 p . For a halo containing N FOF particles, we assign a halo position given by the mean of the individual particle positions, and a halo mass given by: The second term is recommended in [68] to minimize particle resolution artifacts when estimating the mass function using a FOF halo finder.
In this paper we will analyze matter and halo power spectra, using redshifts and halo mass bins defined in Tab. 1. We estimate power spectra by assigning particle positions (or halo positions) to a real-space grid with 1536 3 points using the Cloud-in-Cell algorithm, taking the Fourier transform, and averaging the power over all Fourier modes in a k-bin. This is described in more detail in Appendix A.  Table 1: Mass bins used throughout this paper, with halo number density n, total fraction (by mass) f of dark matter in halos in each bin, and Gaussian bias b G estimated from simulation. The fitting procedure used to estimate b G and assign statistical errors is described in §4 and uses only wavenumbers k ≤ 0.04 h Mpc −1 .

Halo bias
For a halo mass bin i, we define the bias parameter In this section, we will compare values of b mi (k) estimated from simulation with the predicted form: In writing down this predicted form, we have omitted 1-halo terms derived previously in Eqs. (19)- (21). We find that including 1-halo terms does not qualitatively affect any conclusions from this section, but the simplified form in Eq. (27) is convenient for comparison with the rest of the literature. In Fig. 1 we show the dependence of b mi (k) on the non-Gaussianity parameters f N L and ξ, for several choices of redshift and mass bin, on large angular scales (k ≤ 0.04 h Mpc −1 ).
Assigning error bars in Fig. 1 is nontrivial. On large angular scales, variations in P mi (k) and P mm (k) are highly correlated (since variations in both power spectra are mainly due to sample variance) and therefore mostly cancel when we take the ratio to estimate b mi (k). In Appendix A we show in detail how to assign error bars in a way which accounts for this correlation. The resulting where (∆b I ) is derived by replacingn h in Eq. (14) with a non-Gaussian mass function and b G is the Gaussian (i.e. f N L = 0) bias. This functional form is not precisely equivalent to Eq. (27): the two differ at O(f 2 N L ) by the term (2δ c f N L /α(k, z))(∆b I ). We have not investigated whether one of the functional forms is a better fit than the other, but we anticipate that an O(f 2 N L ) difference will be negligible for values of f N L which are observationally relevant. However, we do find that including an O(f N L ) term which is constant in k, either via the last term in Eq. (28) or by allowing b 0 to differ from b G in Eq. (27), is needed to obtain a good fit to simulation. The magnitude of this scale-independent f N L correction we find is in qualitative agreement with what one would get using the non-Gaussian mass-function of [69] in Eq. (14).
We find that for some choices of redshift and halo mass bin, the fits return bad χ 2 values (Tab. 2), i.e. we find statistically significant disagreement between the simulations and the predicted form of the bias in Eq. (27). Detailed inspection of the bad fits shows that, in all cases, the prediction tends to overestimate the magnitude of the non-Gaussian bias on large scales, but the discrepancy between simulation and prediction is only ≈10% of the total non-Gaussian bias in the worst case. The theoretical assumptions made in deriving the non-Gaussian bias, i.e. a universal mass function and the peak-background split relation between halo bias and the mass function, also fail at roughly this level [70,71], so a ≈10% discrepancy is not at all surprising.
Third, we can ask whether the bias b mi (k) is independent of ξ for fixed f N L , as predicted by Eq. (27). To test this, we let b(k) and b (k) denote estimates of the bias from two N -body simulations  with ξ = 0 and ξ = 1 (and the same value of f N L ). We define a χ 2 statistic by summing over k-bins: For almost all redshifts and mass bins, we find acceptable χ 2 values, indicating that the bias estimates from the two simulations are consistent within statistical errors. (There is one exception: we find an anomalous χ 2 for f N L = −500 and z = 2, but inspection of the bad fit shows that the bias only differs by ≈ 10% between ξ = 0 and ξ = 1. The bias is larger in the ξ = 1 case, at low k.) Our conclusion in this section is that the prediction for the non-Gaussian halo bias in Eq. (27) is an impressive fit to the simulations across a wide range of curvaton model parameters. Although we do detect statistically significant deviations from the prediction at the ≈ 10% level, this is typical for results based on general arguments such as the peak-background split.

Halo stochasticity
For halo mass bins i, j, we define the stochasticity parameter If the halos are perfectly non-stochastic tracers of the dark matter (or more precisely, if the only source of stochasticity is shot noise) then both the diagonal (i.e. i = j) and non-diagonal (i.e. i = j) components of r ij (k) will be zero. As discussed in §2.2, the halo model predicts that the leading contribution to r ij (k) on large scales arises from 1-halo terms in the matter-halo power spectra. If ξ > 0, then we expect large r ij since the non-Gaussian part of the bias will be stochastic. In Fig. 2, we show estimates of the diagonal components r ii (k) from simulation, for several choices of redshift, mass bin, and model parameters (f N L , ξ). Similarly to the case of halo bias from the previous section, assigning error bars is nontrivial, since estimates of P mm , P mi , and P ij are all highly correlated on large scales. Our procedure for estimating the stochasticity and assigning error bars is given in Appendix A and results in error bars which are much smaller than sample variance would suggest.
In this section, we will compare the stochasticity estimated from simulations to the predicted form which follows from the halo model calculations in §2.2. As in the previous section, we will separate the issue of whether this prediction agrees with simulation into three separate questions. First, does the prediction in Eq. (31) agree with simulations in the Gaussian case (f N L = 0)? Surprisingly, we do not even find agreement at a qualitative level. Although there are a few choices of redshift and mass bin where the halo model prediction roughly fits the data, there are more cases where there is no resemblance (example good and bad fits are shown in Fig. 3). 4 In this paper, we have not attempted to propose a general model for r ii in the Gaussian case. Therefore, when we compare the prediction for r ij to simulation in the non-Gaussian case, our approach is to estimate the change in stochasticity ∆r ij between non-Gaussian and Gaussian initial conditions, as a function of (f N L , ξ), and compare with the prediction for ∆r ij obtained from Eq. (31).
The second question we can ask is, how does the stochasticity r ii depend on f N L , in the case ξ = 0? In this case, we find that Eq. (31) predicits an f N L dependence which is small compared to the statistical errors of our simulations. To test this prediction, we define a χ 2 statistic by: wherer ii (k) andr ii (k) denote stochasticity estimates from N -body simulations with f N L = 0 and f N L = 0 respectively (taking ξ = 0 in the non-Gaussian simulation). For f N L = ±500 and all redshifts and mass bins considered in this paper, we find good χ 2 values, i.e. no statistically significant dependence of the stochasticity r ii on f N L (provided ξ = 0). This agrees with the halo model       prediction (31), even though the halo model does not correctly predict the actual value of r ii as previously remarked. Third, we can ask, how does the stochasticity depend on ξ? Since we do not have a model for the Gaussian stochasticity, we ask whether the quantity i.e. the excess stochasticity over Gaussian, is correctly modeled by the peak-background split prediction: (This is actually an approximation to the prediction obtained by differencing Eq. (31) between values of (f N L , ξ), but we find that this approximation is within statistical errors of the simulation, so it is a convenient simplification.) We find that the prediction in Eq. (34) systematically overestimates (∆r ii ). Empirically, we find that if we scale the prediction for (∆r ii ) by a multiplicative constant q, then the modified prediction is an excellent fit to the simulations (i.e. when q is treated as a free parameter, all fits have good χ 2 values), for all choices of redshfit, halo mass bin, and f N L . In Fig. 4, we show an example fit; it is seen that the peak-background split overpredicts the level of stochasiticity, but an excellent fit is obtained by simply scaling the peak-background split prediction. In Tab. 3, we tabulate values of q obtained by fitting Eq. (35) to the simulations, together with statistical errors from the fits. One general trend evident in this table is that q is an increasing function of f N L . This simply means that we are considering large enough values of f N L that we are sensitive to O(f 3 N L ) terms, if we think of (∆r ii ) as a power series in f N L . However, it is clear from Tab. 3 that if we interpolate to f N L → 0, most values of q are still < 1. This indicates that, even to leading order O(f 2 N L ), we are seeing ≈ 30% disagreement between the prediction (34) and simulation. We have not attempted to study O(f 3 N L ) terms in detail, or compare them between theory and simulation, since we are already seeing disagreement at leading order O(f 2 N L ). Since the typical value is q ≈ 0.7, our interpretation is that the peak-background split generally overpredicts non-Gaussian stochasticity by ≈30%. Although some discrepancy is expected (e.g. in the previous section we found ∼ < 10% discrepancies in predictions for halo bias), this level of disagreement is somewhat uncomfortable and should probably be incorporated when constraining ξ from observations. It is not clear how to interpret this discrepancy theoretically, or whether it is related to the discrepancy that we found previously in the Gaussian case (Fig. 3). Modeling stochasticity using a semianalytic framework such as the peak-background split or halo model appears to be more difficult than modeling halo bias.

Discussion
In this paper, we have compared semianalytic predictions for halo clustering to N -body simulations, in the two-field inflationary model from [45], in which the initial curvature fluctuation is a sum of  Note that curvaton models also generally predict a cubic contribution of the form (g N L ζ 3 G ) to the initial curvature. We have not considered such a term here since the peakbackground split analysis differs significantly from the f N L case, and defer study of the g N L term to future work [74].
The halo bias b(k) = P mh (k)/P mm (k) in simulation is found to agree very well with the peakbackground split prediction (27) on scales k ≤ 0.04 h Mpc −1 , for a range of redshifts and halo masses. In the Gaussian case, the bias is constant in k at the percent level. In the non-Gaussian case, we find deviations from the functional form for b(k) predicted by the peak-background split which are small ( ∼ < 10%) but statistically significant for our simulation volume. We interpret this as agreement with the prediction, since the peak-background split is expected to break down at the ∼10% level.
We also compare the shot noise subtracted halo stochasticity parameter measured in simulation to semianalytic predictions. In this case the results are more puzzling; some of the semianalytic predictions are confirmed and others are not. In the Gaussian case (f N L = 0), the halo model makes a prediction for r(k) (the leading contribution is from the 1-halo term in P mh ), but this prediction does not consistently agree with simulation. The excess stochasticity (relative to the halo model prediction (31)) observed in simulations is consistent with an additive Poisson-like contribution to the halo-halo power spectrum (i.e. ∆P hh (k) is independent of k on large scales, but can depend on redshift and halo mass).
In the non-Gaussian case, the halo model also predicts that the stochasticity does not depend on f N L (if ξ = 0); we find that this is true in the simulations. The last prediction is a specific functional form (34) for the difference stochasticity (∆r ii ) in stochasticity between a model with ξ > 0 and a Gaussian cosmology. We find that the simulations deviate from this prediction by an overall multiplicative factor q ≈ 0.7, a significant enough disagreement that it should probably be incorporated when constraining the two-parameter curvaton model from observations. It is not clear whether the disagreements between theory and simulation in the Gaussian and non-Gaussian cases are related; it is not straightforward to compare the two since the excess stochasticity appears to be "additive" in the Gaussian case and "multiplicative" in the non-Gaussian case.
Large-scale halo clustering has emerged in the last few years as one of the most powerful probes of primordial non-Gaussianity of one of the local types (i.e. either f local N L , g local N L , or the two-field local type considered here). In this paper, we have confirmed qualitative predictions from the peak-background split ansatz, showing that the peak-background split picture is very useful for relating local-type primordial non-Gaussianity to observations. However, detailed comparison reveals differences which are large enough (≈ 30% in this case) to be important for data analysis, highlighting the need for simulations. The next few years should bring a mixture of theoretical, simulation-based, and observational work which will greatly sharpen our observational constraints on the physics of the early universe.

A Estimators for b and r
Throughout this paper, we have given estimates of the halo bias b i (k) = P mi (k)/P mm (k) and stochasticity r ii (k) = (P ii (k) − 1/n i )/P mm (k) − (P mi (k)/P mm (k)) 2 , with error bars that are used for parameter fitting and computing χ 2 statistics. In this appendix we describe our estimator methodology, in particular the calculation of error bars.
In principle, error bars could be assigned by running multiple N -body simulations and using the Monte Carlo scatter in estimates of b or r, but this is impractical since the "error on the error" would be 2/N mc where N mc is the number of N -body simulations, so computing error bars with 10% accuracy would require running N mc ≈ 200 simulations. Therefore, an analytic prescription for the error bars is necessary.
One can see intuitively that in the limit of zero stochasiticity (r → 0) and zero shot noise (n i → ∞), the statistical errors on b and r should go to zero, since there will be no scatter around the mean relation P ii = b i P mi = b 2 i P mm . Put another way, the statistical errors should not receive contributions from sample variance, since sample variance cancels when we take ratios of power spectra. This will be reflected in our final expressions for the error bars (Eqs. (50) and (51) below), in which all terms contain prefactors of r or (1/n).
In a finite volume V we use the Fourier conventions With these conventions, the infinite-volume two-point function δ(k)δ(k ) * = P (k)(2π In a k-bin b, we estimate power spectra using the estimator where N k = k∈b 1 is the number of Fourier modes in the k-bin, and indices α, β can denote either the matter overdensity field or a halo mass bin i. For the rest of the appendix, we fix the redshift, the k-bin b, and the halo mass bin i. We will make the approximation that variation in power spectra across a single k-bin is not important, and use the compressed notation P αβ = P αβ (k). We will also use compressed notations n = n i for the halo number density, b = P mi /P mm for the bias, and r = (P ii − 1/n i )/P mm − (P mi /P mm ) 2 for the stochasticity.
Define estimators for halo bias and stochasticity by: The factors of (N k − 2)/(N k − 1) and (N k − 2)/N k in the second equation are ad hoc for now, but we will show (Eqs. (48), (49) below) that they ensure thatr is an unbiased estimator of r in the case where N k is not 1. To calculate Var(b) and Var(r), we will make the approximation that the matter overdensity δ m and halo overdensity δ i are Gaussian fields. 5 We would first like to characterize the joint PDF of the power spectrum estimatorsP mm ,P mi andP ii . Let us first define normalized fields whose power spectra are normalized to P 11 = P 12 = 1 and P 12 = 0. The two sets of power spectrum estimators are related by: The joint PDF of the random variables (N kP11 ), (N kP12 ), and (N kP22 ) is well-studied; it is known as the Wishart distribution. A convenient way to characterize this distribution is via the Bartlett-Cholesky decomposition, which states that: where the random variable χ 0 is χ 2 -distributed with N k degrees of freedom, χ 1 is χ 2 -distributed with (N k − 1) degrees of freedom, α is distributed as a unit Gaussian, and the three variables χ 0 , χ 1 , and α are statistically independent.
Combining Eqs. (44) and (45), we find that These equations, combined with the statement at the end of the previous paragraph which gives the joint PDF of χ 0 , χ 1 , and α, completely characterize the sampling PDF of the estimatorsP mm ,P mi andP ii . Armed with this characterization, it is easy to calculate Var(b) and Var(r). We write the estimatorsb andr (defined in Eqs. (41), (42)) in terms of the variables χ 0 , χ 1 , and α (using Eq. (46)), It is then straightforward to calculate the mean and variance of the quantities on the right-hand side, obtaining: Var(r) = 2 N k − 4 r 2 + 2(N k − 2) (N k − 1)(N k − 4) r + 1 nP mm 2 (51) For these calculations the expectation value χ m = 2 m Γ(m+N/2)/Γ(N/2), where χ is χ 2 -distributed with N degrees of freedom, is useful.
This calculation is exact even in the case where N k is not 1, and this level of precision appears to be necessary, e.g. we find that if terms of order (1/N k ) in the variance are neglected, then a few jackknife tests in Tab. 4 below fail. We also note that both estimators are unbiased (i.e. b = b and r = r, justifying the factors of (N k − 1) and (N k − 2) in the definition (42) ofr, which were ad hoc until now. (The estimators would be biased if these factors were omitted.) There is one final wrinkle: in order to apply these expressions for Var(b) and Var(r), we need to know the stochasticity r. In this paper we do not propose a general model for stochasticity, finding for example that the nonzero stochasticity seen in simulations for Gaussian initial conditions is not fit well by the halo model prediction ( §5), and makes a non-negligible contribution to the estimator variance. Therefore, we infer the stochasticity directly from the simulation itself, by making the approximation r ≈P ii − 1/n P mm − P mî P mm 2 (52) on the right-hand sides of Eqs. (50), (51).
Putting the results of this appendix together, our final estimates for Var(b) and Var(r) are given by: We have used these expressions to assign error bars throughout this paper, e.g. in Figs. 1 and 2. We conclude this appendix with an end-to-end test of our estimates for Var(b) and Var(r). If we run two independent N -body simulations with the same values of f N L and ξ, then the differences (b −b ) between bias estimates, and differences (r −r ) between stochasticity estimates, should be consistent with zero. This jackknife test is intended to check correctness of the error bars without requiring a model for the expected values of the bias and stochasticity.
More precisely, for each halo mass bin i, we define a χ 2 statistic by summing over k-bins b: and analogously with the stochasticity estimatorr replaced by the bias estimatorb. Results from the jackknife tests are shown in Tab. 4. The results appear consistent with χ 2 statistics, indicating that our error estimates are accurate. (The most anomalous χ 2 value in the table is 27.6 with 14 degrees of freedom, corresponding to a p-value of 1.6%. Since there are 100 entries in the table, an entry which is anomalous at this level is expected.)