Good and Proper: Self-similarity of N-body Simulations with Proper Force Softening
DDraft version February 19, 2021
Typeset using L A TEX twocolumn style in AASTeX63
Good and Proper: Self-similarity of N -body Simulations with Proper Force Softening Lehman H. Garrison , Michael Joyce, and Daniel J. Eisenstein Center for Computational Astrophysics, Flatiron InstituteSimons Foundation, 162 Fifth Ave.New York, NY 10010, USA Laboratoire de Physique Nucl´eaire et de Hautes ´EnergiesUPMC IN2P3 CNRS UMR 7585Sorbonne Universit´e, 4, place Jussieu, 75252 Paris Cedex 05, France Center for Astrophysics | Harvard & Smithsonian60 Garden St, Cambridge, MA 02138, USA
ABSTRACTAnalysis of self-similarity in scale-free N -body simulations reveals the spatial and temporal scales forwhich statistics measured in cosmological simulations are converged to the physical continuum limit.We examine how the range of scales in which the two-point correlation function is converged dependson the force softening length and whether it is held constant in comoving or proper coordinates. Wefind that a proper softening that reaches roughly 1/30th of the inter-particle spacing by the end of thesimulation resolves the same spatial and temporal scales as a comoving softening of the same lengthwhile using a third fewer time steps, for a range of scale factors typical to ΛCDM simulations. Weadditionally infer an inherent resolution limit, set by the particle mass and scaling as a − / , beyondwhich reducing the softening does not improve the resolution. We postulate a mapping of these resultswith spectral index n = − INTRODUCTIONDetailed comparison of large-scale structure surveydata to numerical simulations is a cornerstone of preci-sion cosmology. In particular, cosmological N-body sim-ulations provide a robust accounting of the gravitationaldynamics that govern the clustering of galaxies in a gen-eral relativistic framework. Although non-gravitationalphysics is ignored in such simulations, they provide ascaffolding for the placement of galaxies with techniqueslike the halo occupation distribution (HOD) that canmarginalize over a plausible range of un-modeled physics(e.g. Zheng et al. 2005; Wechsler & Tinker 2018 for a re-cent review).As the resolution of simulations increase, techniquesthat make use of small-scale structure have ballooned,including sub-halo abundance matching (SHAM, Con-roy et al. 2006; Vale & Ostriker 2006) and tech-niques that use the matter field of the simulations(e.g. GRAND-HOD; Yuan et al. 2018) rather than as-suming dynamical equilibrium and a profile like NFW
Corresponding author: Lehman Garrisonlgarrison@flatironinstitute.org (Navarro et al. 1997) or Einasto (Einasto 1965). Mod-ern galaxy surveys resolve increasingly faint, less mas-sive objects, residing in smaller halos. In such a context,the convergence of the N -body simulation must be es-tablished with respect to discretization parameters likethe number of particles N , the associated mean inter-particle spacing L/N / for a given box size L , and soft-ening length (cid:15) that regularizes the small-scale force.The precision of simulations matters, too: emulatorsand other techniques that rely on the derivatives of ob-servables with respect to cosmology require relative ac-curacy. If the numerical errors in N -body simulationsare cosmology-dependent, then such derivatives will notbe faithfully reproduced. There are hints of such cos-mology dependence, as seen in the redshift dependenceof errors in the Euclid code comparison project, for ex-ample (Schneider et al. 2016).In principle, the convergence of N -body simula-tions to the relevant physical limit—the Vlasov-Poissonsolution—can be assessed by extrapolating the dis-cretization parameters appropriately, such as increasingparticle count N . In practice, an obstacle to establish-ing convergence in this way is the computational costof increasing particle density while also decreasing theforce softening length (cid:15) (e.g. Power et al. 2003). Indeed, a r X i v : . [ a s t r o - ph . C O ] F e b Garrison et al. the the appropriate extrapolation for softening does notsend (cid:15) to zero while holding the particle spacing fixed, asthis amplifies two-body scattering. Instead, (cid:15) must betaken small compared to the physical scales of interest,but still large compared to the particle spacing. Such aregime is, in practice, numerically inaccessible.Comparison of the N -body simulation to a referencesolution is therefore desirable but elusive in practice.One class of test that circumvents some of these diffi-culties is analysis of self-similarity of scale-free simula-tions (e.g. Efstathiou et al. 1988; Colombi et al. 1996;Jain & Bertschinger 1998; Smith et al. 2003; Widrowet al. 2009). Scale-free simulations are initialized witha power-law power spectrum in an Ω M = 1 (EdS) cos-mology, such that there is only one physical scale in thebox—the scale of the onset of non-linearity. The clus-tering of matter on small scales at early times should bea rescaling of the clustering on large scales at late times,else the simulation is said to break self-similarity.All simulations will break self-similarity at some level,due to UV cutoffs (finite N or softening), IR cutoffs(finite L ), or other preferred scales imprinted by themethod (e.g. cell size in the force solver). But trackingthese deviations from self-similarity, and validating thata simulation method has not imprinted any additionalscales on the clustering, provides bounds on the range ofmasses, times, and scales that may possibly be trustedfrom a simulation.With this in mind, we have begun a program of ex-ploring the resolution and convergence of scale-free sim-ulations and their data products. In Joyce et al. (2020),we initiated this program with an analysis of the two-point correlation function (2PCF) of matter in a singlesimulation, employing an n = − a − / at late times. Further, it con-cluded that these resolution limits at small scales areset essentially by the mass resolution and should be in-sensitive to softening provided it is sufficiently small.In this work, we vary the softening to identify directlyits role in setting the small-scale resolution cut-off, andin particular disentangle its contribution from that offinite mass resolution. We test several different softeninglengths, fixed in comoving coordinates, and several fixedin proper coordinates.Of course, ΛCDM simulations are not power-law, EdScosmologies, and contain many scales besides the onsetof non-linearity. But the non-linear scale is present inboth, and does provide a route to creating a mapping from scale-free simulations to ΛCDM simulations, as inSmith et al. (2003) or Joyce et al. (2020). We discussthis below briefly and will explore it more fully in futurework.Self-similarity in a scale-free test does not necessarily“prove” that a method is correct or free of errors. Thereexist classes of systematic errors that are themselvesscale-free—for example, a global time step of any fixedlogarithmic size, even implausibly large, will exhibit self-similarity. However, the primary method we are testingin this work—force softening—manifestly breaks self-similarity. We thus have a strong theoretical expecta-tion that we should be able to identify scales where self-similarity is broken, and that scales where it is preservedare converged with respect to force softening.The plan of this paper is as follows. In Section 2, welay out some theoretical groundwork on N -body forcesoftening. In Section 3, we introduce the scale-free unitsin which we will analyze the simulations and define thesimulation parameters—a set of identical simulations ex-cept for the force softening length and technique. In Sec-tion 4, we identify the length scales and epochs for whichthe simulations exhibit self-similarity, and consider howthese scales change with softening. We summarize inSection 5. SOFTENINGForce softening in N -body simulations regularizes thesmall-scale gravitational force to mitigate the effects oftwo-body scattering. To the extent that N -body sim-ulations attempt to model collisionless Vlasov-Poissondynamics, collisional encounters are non-physical andshould be suppressed. However, softening also modifiesthe small-scale growth of structure, so it should be keptmuch smaller than the physical scales of interest. Thechoice of softening length is thus a balance between am-plifying unphysical effects at smaller scales and “wash-ing out” clustering on large scales. A smaller softeninglength additionally imposes the computational cost ofintegrating tighter orbits, requiring more time steps.Choice of softening length has been studied exten-sively (e.g. Knebe et al. 2000; Power et al. 2003; Die-mand et al. 2004; Power et al. 2016; Mansfield &Avestruz 2021). In principle, the proper procedure toselect a softening length is to choose the halo mass ofinterest, select a softening a few times smaller than theexpected core radius of such objects, and then choose N to densely sample the softening length; i.e. (cid:15) (cid:29) L/N / .Such a procedure results in smooth, well-sampled dy-namics, with limited two-body relaxation and the UVcutoff set only by the softening length. Of course,increasing N to such extremes is computationally in- elf-Similarity of Proper Softening /r at afinite radius), correspondence to some analytic densityprofile, and avoidance of dynamical instabilities. Forcosmological N -body simulations, the softening lengthseems to matter more than the functional form, exceptto note that non-compact softening laws can extend thesuppression of growth to scales many times larger than (cid:15) (Garrison et al. 2019).In Plummer softening (Plummer 1911), the F ( r ) = r /r force law is modified as F ( r ) = r ( r + (cid:15) p ) / , (1)where (cid:15) p is the softening length. This softening is veryfast to compute but is not compact, meaning it neverexplicitly switches to the exact r − form at any radius.Therefore, many N -body codes, such as gadget (Springel et al. 2020) and Abacus (Garrison et al.,in prep; Garrison et al. 2019 for a recent description),adopt a spline softening.
Abacus uses the follow-ing spline softening from Garrison et al. (2016), whichswitches to the exact form past radius (cid:15) s : F ( r ) = (cid:2) − r/(cid:15) s ) + 6( r/(cid:15) s ) (cid:3) r /(cid:15) s , r < (cid:15) s ; r /r , r ≥ (cid:15) s . (2)The softening scales (cid:15) s and (cid:15) p imply different mini-mum dynamical times (an important property, as thissets the time step necessary to resolve orbits). We al-ways choose the softening length as if it were a Plum-mer softening and then internally convert to a softeninglength that gives the same minimum pairwise dynami-cal time for the chosen softening method. For our spline,the conversion is (cid:15) s = 2 . (cid:15) p .While the softening length is usually taken to be fixedin comoving coordinates, it is sometimes also fixed inproper coordinates, especially in hydrodynamical simu-lations where galactic structures are effectively “frozenout” of cosmic expansion. Formally, N -body leapfrogintegration is usually defined in terms of comoving co-ordinates and canonical momentum, such that any time-varying softening formally breaks the symplectic natureof the integration (Quinn et al. 1997). However, the pro-cedure of choosing a time step is often non-commutativewith the kick and drift operators and therefore breaks Scale factor a / a C o m o v i n g s o f t e n i n g ProperComoving
Figure 1.
Softening length versus scale factor for the fidu-cial simulations with comoving and proper softenings. Theproper softening length has a greater comoving value for allepochs considered (Eq. 3 with (cid:15) prop = 0 . the symplectic property anyway. Time stepping schemesare an active area of work in N -body (e.g., recently,Springel et al. 2020; Zhu 2021). For practical purposes,in cosmological simulations, the computational economyof adaptive time stepping is often preferred over the for-mal advantage of symplecticity.We test both comoving and proper softening in thiswork. Comoving softening is constant in time, since thesimulation uses comoving coordinates. Proper softeningincreases in comoving coordinates towards earlier times.Since cosmological simulations typically run for a factorof 10 to 100 in scale factor, proper softening can becomeinappropriately large at early times if unchecked (e.g. asoftening of 1 / z = 0will span three particles in the initial configuration at z = 99). This will suppress linear growth of structurebeyond what already happens due to the discretizationof the mass field (Garrison et al. 2016). Furthermore, itwill introduce transients, because the Lagrangian per-turbation theory assumes a continuum of mass elementsthat evolve without softening; the introduction of soft-ening at the start of the simulation changes the growthrate discontinuously (unless one solves the Lagrangianperturbation theory with direct force evaluations, as inGarrison et al. 2016).Therefore, we cap the proper softening at high redshiftas (cid:15) ( a ) = min (cid:16) a (cid:15) prop a , . (cid:17) . (3) (cid:15) is the comoving epsilon employed by the simulation atscale factor a , for a given proper softening (cid:15) prop definedat a . This is shown in Figure 1 for (cid:15) prop = 0 .
3. Ef-fectively, at early times the softening is comoving andfixed to 0.3 in units where the particle spacing is unity,so the worst of the growth rate modifications should be
Garrison et al. suppressed as it is smaller than the mean interparticlespacing. Then, at some intermediate redshift, when theproper softening is smaller than 0.3, the softening beginsshrinking in inverse proportion to the scale factor. Notethat with this definition, larger (cid:15) prop does not changecomoving plateau value, but only delays the transitionfrom comoving to proper evolution.Fixing the softening in proper coordinates, (cid:15) prop = a(cid:15) com , has the effect of increasing the comoving soft-ening length at earlier times. One should expect theusual effects of increasing softening: suppressed growthof structure on small scales, but decreased relaxationfrom impulsive encounters.In a simulation with only one particle species, masssegregation from two-body relaxation is not a concern.However, impulsive encounters may still inject kineticenergy into a fraction of particles, stealing potential en-ergy from the remaining particles, leading to core col-lapse and/or evaporation. In practice, the time scale ofcore collapse is tens to hundreds of relaxation times andis unlikely to occur on cosmological time scales; evapo-ration takes even longer (Knebe et al. 2000; Binney &Knebe 2002; Binney & Tremaine 2008). But the popu-lation of particles with high kinetic energy is a genericfeature that may still be present, perhaps manifestingas “puffier” halos. SCALE-FREE SIMULATIONS3.1.
Definitions
In scale-free N -body simulations, the particles are ini-tialized with a power-law power spectrum of index n inan Ω M = 1 background cosmology. The linear powerspectrum grows as: P L ( k, a ) = Aa k n (4)for scale factor a and amplitude A .The only physical scale in the problem is the non-linear scale, defined by the onset of large-amplitude mat-ter clustering. In configuration space, we may define thisusing σ ( R, a ), the variance of density in spheres of ra-dius R : σ ( R, a ) = (cid:90) d k (2 π ) ˜ W T ( kR ) P L ( k, a )= 9 Aa R − (3+ n ) π / Γ[(3 + n ) / − n ) / n − n − , (5)where ˜ W T is the Fourier transform of the sphericaltophat window function. Eq. 5 is written in a numeri-cally stable form that avoids divergence of the Γ func-tions for integer n ∈ ( − ,
1) (Garrison 2019). We may define R nl as the scale at which the dimen-sionless variance reaches unity, using Eq. 5 to find aself-similarity relation: σ ( R nl , a ) = 1 a R − (3+ n )nl ∝ R nl ∝ a / (3+ n ) . (6)Eq. 6 defines a map of redshift to scale. In other words,any function of scale, such as the two-point correlationfunction, must be a rescaling of the same function at adifferent epoch. Indeed, any dimensionless function ofscale r/R nl must be constant so long as self-similarityholds. 3.2. Abacus N -body Code We employ the
Abacus N -body code (described inGarrison et al. 2019; code paper in prep.) for the sim-ulations in this work. Abacus is a high-performance,GPU-based code that uses an analytic separation of thenear- and far-field forces to solve the former with exactsummation and the latter with a high-order multipolemethod ( p = 8 in this work). The result is highly accu-rate forces, such that the remaining error is primarily thefinite time step size of the leapfrog integration (Garri-son et al. 2019). In Joyce et al. (2020), we demonstratedthat a time step parameter of η acc = 0 .
15 is sufficientfor the 2PCF of matter to be converged to 0.1% for an n = − Abacus employs a global leapfrog time step, sharedby all particles. The size of the time step, ∆ a , is cho-sen at the beginning of each time step based on the cellwith the smallest v rms /a max , scaled by the time stepparameter η . A global v rms /a max is also computed andused instead of the cell-based value if it is larger—thisguards against abnormally cold cells causing catastroph-ically small time steps.Additionally, the time step is capped to 0.03 in unitsof d ln a , ensuring at least 33 steps per e -fold of the scalefactor. This criterion dominates at early times, beforeclustering has brought about small dynamical times inthe centers of halos. This ensures that the accuracy ofearly linear-theory growth.3.3. Simulation Parameters
The fiducial simulation is a scale-free realization of a n = −
2, Ω M = 1 cosmology with N = 1024 particles.We work in units of the mean inter-particle spacing, or L/N / = 1.The first output epoch, a , is chosen based on the top- elf-Similarity of Proper Softening σ ( R = 1 , a ) = 0 . . (7)The value of 0 .
56 = 1 . / √ M nl = 0 . . (8) M nl may be defined proportionally to the non-linearlength scale R nl : M nl ∝ R . (9)From this and Eq. 6, we may determine the spacing ofour outputs in scale-factor:∆ log a nl = d ln ad ln M ∆ log M nl , = 3 + n M nl , = 1 / . (10)In the last line, we used the index n = −
2. We produce38 outputs, or about a factor of 5 × in M nl and 80in R nl .The simulations use 2LPT initial conditions with theconfiguration-space method of Garrison et al. (2016).The initial scale factor, a i , is chosen so that the densityfluctuations are small compared to unity: σ ( R = 1 , a i ) = 0 . . (11)We employ particle linear theory (PLT) corrections inthe initial conditions to null out transients that arisefrom the particle discretization of the matter densityfield (Joyce & Marcos 2007; Garrison et al. 2016). Onsmall scales, we apply a wavevector-dependent rescalingof the amplitude of the power spectrum such that thesimulations arrive at the correct linear theory value at a P LT , even though the growth rate remains unavoidablymodified. We choose a P LT = a ; we have tested varia-tions in this choice, including disabling rescaling entirely,which generally worsen self-similarity but do not affectour qualitative conclusions.This work focuses on variations in softening and timestep, which are summarized in Table 1.3.4. Softening Prescriptions
All softening in this work uses
Abacus spline soften-ing, defined in Equation 2. As discussed in Section 2,proper softening can become inappropriately large atearly times, spanning multiple particle spacings and
Table 1.
Simulations used in this work. For proper softening, (cid:15) ( a )corresponds to (cid:15) prop in Eq. 3 (with a defined by Eq. 7). N Softening (cid:15) ( a ) η Comment1024 Comoving 1 /
30 0.15 Fiducial comoving0.075, 0.3 Time step variations1 /
60, 1 /
15 Softening variations1024 Proper 0.3 0.15 Fiducial proper0.075, 0.3 Time step variations0.1, 0.5, 0.9 Softening variations thereby modifying the linear growth of structure. There-fore, the softening is capped at early times to 0.3, asdetailed in Equation 3.In our fiducial simulations the proper softening isgreater than the comoving softening for all epochs con-sidered, starting 9 × greater for a < a and ending atan almost identical value. This increases the smallestdynamical times in the simulation, requiring the simu-lation to take fewer time steps—an advantage in com-putational economy. In Section 4.5, we will discuss thatthe fiducial proper simulation takes 65% as many timesteps as the fiducial comoving.Because we are employing compact spline softening,there is an exact switchover radius from softened forceto 1 /r (see the discussion in Sec. 2). For a comovingsoftening of 0.3 (the maximum value allowed by Eqn. 3),that switchover radius is about 0.65. It is not 0.3 be-cause all our softening lengths are quoted as effectivePlummer values to facilitate comparison between dif-ferent softening schemes. Since this value is less than1, we expect that the most serious effects of softeningmodifying the initial growth rate discontinuously fromthat assumed by the initial conditions generator will besmall. 3.5. Analysis
As in Joyce et al. (2020), our analysis methodologyis to evaluate the matter 2PCF of the simulations atall epochs, apply the self-similar rescaling to cast thecomparison in units of R nl , and observe how quicklythe simulations converge to the self-similar solution fordifferent length scales. This is presented in Section 4.We use the Corrfunc code of Sinha & Garrison(2019, 2020) to compute the correlation functions. Theradial bins are scaled self-similarly for each epoch, sono interpolation or centering corrections are requiredto compare measurements across epochs. Particles aredown-sampled by a factor of 2 before pair counting, andlate-time, large-scale bins are elided. These bins are not
Garrison et al. r ( r ) R nl Figure 2.
Correlation functions from 38 epochs of the fidu-cial simulation with comoving softening. Line color corre-sponds to non-linear scale R nl , a label of the epoch, withlighter color indicating later epoch. Oscillations from thenear-lattice structure of the initial conditions are evident atearly times near r = 1 (the particle spacing scale). Thecomoving softening length, (cid:15) = 1 /
30, is marked. needed to assess small-scale convergence and indeed be-gin to exhibit finite box-size effects, even at 1/100th ofthe box scale, due to the very red n = − RESULTS4.1.
Correlation Functions
Figures 2 & 4 show the raw 2PCF measurements forthe fiducial simulations with comoving and proper soft-ening, respectively. Figures 3 & 5 show their self-similarrescalings, using Eq. 6.First, we see that the rescaled correlation functionslargely exhibit self-similarity—that is, they stack. Butthere are prominent deviations demonstrating that self-similarity is broken. On small scales, we see “fanning”,or flattening, due to a combination of softening and finiteparticle mass. At somewhat larger scales, we see dipsand oscillations, more clearly seen in the raw correlationfunctions. There, we observe that these early-time oscil-lations fall around r = 1, and are thus the “memory” ofthe near-lattice structure of the initial conditions (theparticle spacing is unity). This feature is washed outafter a factor of a few in R nl as the simulation Poisson-izes.To more precisely assess how quickly different scalesconverge, we take vertical slices through these rescaledfigures and plot each slice as a panel of, e.g., Figure 6and subsequent figures. Flat lines—constant correlationamplitude as a function of epoch—indicate convergenceto self-similarity. The range of epochs and r/R nl scalesthat exhibit this convergence are what we refer to as r / R nl ( r ) Figure 3.
The self-similar rescaling of Figure 2. The lineslie atop one another, except where self-similarity is broken.On the left, the “fanning” is due to softening and finite massresolution, and the oscillations from the initial particle latticeare still seen at early times. By taking vertical slices throughthis plot, we may assess the rate of convergence as a functionof epoch (Figure 6 and onwards). r ( r ) R nl Figure 4.
Same as Figure 2, but for the fiducial simulationwith softening held constant in proper coordinates. The co-moving softening therefore shrinks towards later epochs, asshown by the dashed line. “resolved scales”. The impact of variations in softeningon resolved scales will be discussed in Section 4.A few small, sharp wiggles are visible on the smallestscales of the 2PCF measurements, below the softeningscale. This is due to the lossy compression applied to theparticle data that truncates the precision of the positionsto 14 bits (as cell offsets, so 22.7 bits of global position).4.2.
Comoving vs. Proper Softening
In Figure 6, we compare the self-similarity of the fidu-cial simulation with comoving softening ( (cid:15) = 1 /
30) tothat with proper softening ( (cid:15) ( a ) = 0 . r/R nl (different panels), elf-Similarity of Proper Softening r / R nl ( r ) Figure 5.
The self-similar rescaling of Figure 4. which, when there is convergence to self-similarity, cor-respond to different ξ amplitudes.First, we see that the overall self-similarity—the rangeof plateaued ξ values—is strikingly similar between co-moving and proper. Indeed, there is even mild evidencethat the proper softening increases the range of resolvedepochs (e.g. second panel). However, at early epochs(small R nl ), the clustering with proper softening falls offmore steeply. In other words, the comoving and propersimulations converge to the same limit over the samescales, but outside of those scales, the proper softeningfalls off more steeply.Additionally, despite the proper softening being largerthan the comoving softening for the entire durationof the simulation (Figure 1)—starting from 9 × largerat initial output a —the proper softening arrives at amore clustered state on small scales at late times (firstpanel). This is somewhat surprising, as considering thegrowth of modes, a larger softening decreases the growthrate. But considering halo dynamics, it is plausible thatproper softening decreases two-body relaxation, which isdependent on the amplitude of the impulsive scattering,leading to more concentrated halos.We can quantify the range of resolved scales by iden-tifying the epoch R nl at which the correlation functionfirst reaches 5% of its plateau value. Repeating this foreach slice of r/R nl (each panel), we can map the smallestresolved scales versus a/a . This is shown in Figure 7.The axes are no longer in rescaled units; this analysisconnects scale-free units back to comoving units of theinter-particle spacing.We see that the comoving and proper softeningslargely agree which scales are resolved, with some evi-dence of better convergence (smaller resolved r ) with theproper softening at early and late times, as observed inFig. 6. But the measurements are noisy due to binningand finite number of output epochs, so we do not fo- cus on detailed differences—as small as one bin in mostcases—but rather broad agreement.The resolved r propagates to progressively smallerscales as a − / to a good approximation. This supportsthe hypothesis set forth in Joyce et al. (2020) that theresolution is set by the two-body relaxation of clusteringthat is approximately constant in proper coordinates—so-called “stable clustering” (Efstathiou et al. 1988;Colombi et al. 1996; Jain & Bertschinger 1998; Widrowet al. 2009).The comoving and proper softening lengths aremarked by dashed lines in the figure. The proper soft-ening evolves with time in comoving units, but at allresolved epochs, the resolution limit is above the propersoftening. So although we achieve similar or betterresults than with comoving softening by employing aproper softening that is larger for all epochs, it is likelythat the proper softening must still be smaller than theresolution limit set by the particle mass, lest the sim-ulation be softening-limited. This is supported by ourfindings in the next section.While we have focused on the range of resolved scales,we should also revisit the behavior outside those scales.As mentioned above, proper softening suppresses the un-resolved clustering at early times (although sometimesenhances it at late times). The suppression effect canbe large—a factor of several, as one can see comparingFigures 2 & 4 by eye. Where using such scales can-not be avoided, it is already common practice to applyflexible analysis to marginalize over the details of suchclustering. These results confirm the importance of suchrobust analysis methods and understanding where suchunresolved scales begin.In Fig. 6, for large r/R nl , the comoving and propersoftenings agree perfectly, as we would expect for scaleswell above the UV cutoff. Deviations from self-similaritystart to appear on the largest scales at the latest times,which are most likely finite box size effects. The lengthscales at which such effects start to appear, around1/300th of the box size, are smaller than normal large-box, ΛCDM intuition would indicate, but are not sur-prising given the power-law power spectrum with a red, n = − Softening Length
We vary the softening length of the comoving simula-tion in Figure 8 and the proper simulation in Figure 10.In the comoving case, we consider (cid:15) = 1 /
60 and 1 / /
30. We find that the rangeof resolved scales substantially improves from (cid:15) = 1 / /
30, showing nearly a factor-of-two improvementfor large correlation amplitudes (small r/R nl ). Halv- Garrison et al. r / R nl = 0.00715 1000400600800 r / R nl = 0.0227200300400 r / R nl = 0.0718 100406080 r / R nl = 0.21510820 r / R nl = 0.681 11.251.522.5 r / R nl = 2.160.1 10.2 0.3 0.6 2 3 6 R nl r / R nl = 6.46 0.1 10.2 0.3 0.6 2 3 6 R nl r / R nl = 20.5 Figure 6.
The self-similar convergence of the fiducial simulation for both comoving and proper softening. Each panel representsa vertical slice in r/R nl (see Fig. 3), with flat lines indicating self similarity. The range of resolved scales is remarkably similarin both, but with proper softening exhibiting a sharper rolloff at early times. Each panel shows 0.5 dex in ξ . ing the softening again from 1 /
30 to 1 /
60 shows onlymarginal improvement; the range of resolved epochs in-creases only about 10%. In unresolved epochs, the clus-tering amplitude does increase, but again, not as muchas the first halving from 1 /
15 to 1 / a − / at late times. This result is consis-tent with a hypothesis that the mass resolution limitsthe resolved scales through two-body scattering. There-fore, if we were to run our simulations for longer—to higher clustering amplitudes—the convergence wouldpropagate to small enough length scales that the comov-ing softening would once again be the limiting factor.Therefore, the choice of softening length for comovingsoftening must be made not only with respect to massresolution but also with respect to the range of scalefactors over which the simulation will be run.Proper softening, which scales as a − , shrinks fasterthan the resolution, so it does not share this consid-eration (Fig. 7 shows these scalings). Instead, it mustbe kept smaller than the resolved scales at early times;e.g. near a , or whatever the first epoch of analysis is. elf-Similarity of Proper Softening a / a R e s o l v e d r comprop ( a ) ComovingProper a Figure 7.
The minimum resolved scale r versus epoch forthe fiducial comoving and proper simulations. r is given inunits of the comoving inter-particle spacing, and the epochis relative to the first output epoch a (Eq. 7). The comov-ing simulation has (cid:15) = 1 /
30 (blue dashed line shows 3 × thisvalue) and the proper (cid:15) ( a ) = 0 . a − / . As long as this condition is satisfied, the proper soften-ing will always be smaller than the resolved scales foranalysis epochs.Figure 9 quantifies the resolved scales versus epochfor different values of the comoving softening length. Aswith Fig. 7 for the comparison of comoving and propersoftening, we compare the epochs at which the simula-tion reaches 5% of the plateau value in Fig. 8 and infer aminimum resolved scale in units of the comoving inter-particle spacing. We focus on overall behavior ratherthan detailed differences, because our precision is lim-ited by binning and the spacing of our output epochs.We see that (cid:15) = 1 /
30 and (cid:15) = 1 /
60 agree rather well,with the minimum resolved scale improving by muchless than a factor of two—indeed, only by one bin inepoch. Doubling the softening length from (cid:15) = 1 /
30 to (cid:15) = 1 /
15 degrades the resolution by a similarly smallamount, except at late epochs where the degradationbecomes large. Indeed, we see nearly the full factor of 2one would expect if the resolution were determined fullyby the softening at these late epochs. As anticipatedin our discussion above, the effect is likely only seen atlate epochs because the resolution set by the particlesmass propagates to smaller scales as the simulation pro-gresses (hypothesized as a − / ). When this resolutionapproaches a few times the softening length, the soften-ing length becomes the limiting factor. In this case, weobserve degradation at about 3 × softening length. We vary the softening of the proper softening simula-tion in Figure 10, using values (cid:15) ( a ) = 0 . (cid:15) ( a ) = 0 . (cid:15) ( a ) = 0 .
9, in addition to the fiducial value of (cid:15) ( a ) = 0 .
3. These epsilon values are given at a , theepoch of first output, but recall that the resulting co-moving value is capped at early times (Eq. 3).Figure 11 shows the inferred resolution for each soft-ening. The first factor-of-3 decrease from (cid:15) ( a ) = 0 . (cid:15) ( a ) = 0 . . .
1, makes almost nodifference—at most one bin of epoch. The clustering inthe unresolved regions does benefit somewhat, however(Fig. 10). We infer diminishing returns past (cid:15) ( a ) = 0 . (cid:15) ( a ) = 0 . Time Step
Time step accuracy is difficult to assess with self-similarity analysis (see fig. 2 and discussion in Joyceet al. 2020). In particular, if a simulation uses a con-stant, global leap frog time step in log( a ), then the sim-ulation would be expected to be self-similar for any stepsize, due to the scale-free nature of logarithmic timestepping, even though the clustering is far from the phys-ical solution. The Abacus time step scheme is not quitethis simple—using a constant log time step for the ear-liest times and switching to a dynamical criterion at theonset of clustering ( § η = 0 .
3, 0 .
15 (the fiducial), and0 . . ξ/ξ . ratio as we halve from η = 0 . η = 0 .
15 and then again to η = 0 . Garrison et al. r / R nl = 0.0071520003000 = 1/60= 1/30= 1/15 r / R nl = 0.0227 1000600800 r / R nl = 0.0718200300400 r / R nl = 0.215 100406080 r / R nl = 0.68110820 r / R nl = 2.16 11.251.522.5 r / R nl = 6.460.1 10.2 0.3 0.6 2 3 6 R nl r / R nl = 20.50.1 10.2 0.3 0.6 2 3 6 R nl Figure 8.
The self-similar convergence with respect to comoving softening length. constant as a function of r/R nl . Or, perhaps more in-tuitively, the time step error is only dependent on thevalue of ξ , rather than the physical scale. Therefore,Figure 12 shows one 2PCF ratio for each simulation forseveral r/R nl values (the same ones considered in theconvergence of the softening length). In detail, there issome small evolution of the 2PCF ratio even in scale-freecoordinates, so the mean ratio across epochs is plotted,along with the standard deviation as an error bar.We see that halving the time step from our fiducialchoice in this work, η = 0 .
15, to 0.075 makes almost nodifference, with the largest difference of 0.2% seen in thesmallest scales. Therefore, we consider our fiducial timestep converged for the purposes of this work.The proper softening simulation is seen to exhibitslightly less dependence on time step than the comov-ing simulation. At early times, the proper simulationhas likely taken more time steps relative to the dynam-ical times of collapsing objects because
Abacus ’s early log-constant time step. The η = 0 . η = 0 .
15 simulation takes 108 time steps—clearlynot double. Therefore, one might posit that the propersoftening’s relative insensitivity to η is because the first100 steps make up a larger fraction of the time steps.By the time the last output is reached, the propersimulation has taken 667 steps with η = 0 . η = 0 .
15, and the comoving simulation 997 and1964 steps—ratios of 1.88 and 1.97, respectively. Whilethe direction of this effect accords with our hypothesis, itis unclear if this small difference can explain the propersimulation’s error at η = 0 . e -fold of the simulation, where the correlation amplitudesreach several thousand. The fiducial comoving simula-tion with η = 0 .
15 takes about 1000 time steps to com- elf-Similarity of Proper Softening a / a R e s o l v e d r = 1/15= 1/15= 1/30= 1/60 Figure 9.
Same as Fig. 7 but for comoving softening only,showing variations in softening length. As we observe inFig. 8, halving the softening length from 1/15 to 1/30 pro-duces a noticeable gain in resolved scales (smaller r ) at lateepochs where the resolution, nominally improving as a − / (black dashed line) approaches the softening. Halving againto 1/60 produces a much smaller gain; indeed, this gain isonly one bin in epoch. The dark purple dashed line showsthe largest of the three softening lengths, 1 / plete the last e -fold, and the fiducial proper about 850.An AbacusSummit ΛCDM simulation would take about1200 steps for the last e -fold if it were to use that sametime step parameter, with a comoving softening of 1 / × h − M (cid:12) .4.5. Run Time
To reach the final output, the fiducial simulation took1964 steps with comoving softening and 1252 steps withproper softening. This is approximately 65% as manysteps, coming at apparently no cost to the range of re-solved scales, but suppressed clustering in unresolvedscales, as discussed in Section 4.2. SUMMARYWe have evaluated the self-similarity of scale-free N -body simulations with an n = − (cid:15) = 1 /
15 to (cid:15) = 1 / (cid:15) = 1 /
30 to (cid:15) = 1 /
60, showsmuch less improvement.We have found, as reported by Joyce et al. (2020), thatresolution propagates to smaller scales at later times as a − / . This explains the importance of a smaller soft-ening at later times: the softening is not the limitingfactor until this resolution, set by the particle mass, ap-proaches the softening scale. This is consistent with themass-dependent resolution being mediated by two-bodyrelaxation of halos whose clustering is constant in propercoordinates (stable clustering). Once this mass resolu-tion limit is reached, we infer there is no gain in resolvedscales to be made by taking the resolution smaller.Considering only proper softening and varying thesoftening length, we similarly have found diminishingreturns beyond our fiducial choice. We note that ourchoice runs close to the early-time resolution limit setby the particle mass, such that increasing further de-grades resolution in this regime. However, since propersoftening shrinks in comoving coordinates faster thanthe resolution ( a − versus a − / ), the late-time cluster-ing should be less susceptible.With the mass resolution limit scaling as a − / , itis natural to consider whether a softening length thatscales in the same way would be advantageous. Whilethis may match the two UV scales well, we do considerproper softening to have a better physical motivationthan such a softening. With proper softening, we sup-press secular evolution of halo profiles in proper coordi-nates due only to the change in softening length.With respect to time step, we have shown that an Abacus time step size of η = 0 .
15 is converged for theaccuracy required in this work, and that the time steperror is approximately self-similar—in other words, de-pendent only on the amplitude of the 2PCF, not theepoch. The proper softening shows relatively less sen-sitivity to the time step parameter than the comovingsoftening.With respect to ΛCDM, one might construct a map-ping to these n = − a for theΛCDM simulation using the epoch of collapse of the firsthalos (Eq. 7), as elaborated in Joyce et al. (2020). Wewill explore the accuracy of such a mapping in futurework, including the sensitivity to spectral index. How-ever, because the small-scale ΛCDM index is more neg-2 Garrison et al. r / R nl = 0.0071520003000 ( a ) = 0.1( a ) = 0.3( a ) = 0.5( a ) = 0.9 r / R nl = 0.0227 1000600800 r / R nl = 0.0718200300400 r / R nl = 0.215 100406080 r / R nl = 0.68110820 r / R nl = 2.16 11.251.522.5 r / R nl = 6.460.1 10.2 0.3 0.6 2 3 6 R nl r / R nl = 20.50.1 10.2 0.3 0.6 2 3 6 R nl Figure 10.
The self-similar convergence with respect to proper softening length. The softening values (cid:15) ( a ) are defined at theepoch a of initial output. ative than n = −
2, a ΛCDM simulation will have moreclustering at scales above the radius whose mass vari-ance defines a . Since structure is collapsing on largerscales for a given a —the mass variance is flatter withrespect to scale—we expect the sensitivity to softeninglength to be smaller. Therefore, we expect that ourfindings with respect to required softening length areconservative.Aside from spectral index, our analysis was also re-stricted to the 2PCF. We will explore other summarystatistics and spectral indices in future work. Such ex-tensions will allow robust extrapolation of these resultsto the ΛCDM simulations that underlie much of the in-terpretation of large-scale structure data. ACKNOWLEDGMENTSDJE is supported by U.S. Department of Energy grantDE-SC0013718 and as a Simons Foundation Investiga-tor. Abacus development has been supported by thosefunds and by NSF AST-1313285 and Harvard Universitystartup funds. We would like to thank the co-authors of Abacus : Nina Maksimova, Doug Ferrer, Marc Metch-nik, and Philip Pinto.
Software:
NumPy (Harris et al. 2020), astropy (As-tropy Collaboration et al. 2018), Matplotlib (Hunter2007) DATA AVAILABILITYThe correlation function measurements used in thiswork are available upon request. The underlying simu-lation data is substantially larger but may also be madeavailable upon reasonable request. elf-Similarity of Proper Softening a / a R e s o l v e d r ( a ) = . ( a ) = 0.9( a ) = 0.5 ( a ) = 0.3( a ) = 0.1 Figure 11.
Same as Figs. 7 & 9 but for proper soften-ing only, showing variations in softening length. As seen inFig. 10, taking (cid:15) ( a ) from 0 . . . . r . Going from 0 . . . a − / for thelower resolution limit as set by the particle mass. r / R nl = 0.007150.970.980.991.00 / . r / R nl = 0.0227 0.970.980.991.00 r / R nl = 0.0718 / . ProperComoving r / R nl = 0.215 Figure 12.
Time step convergence for the fiducial comoving and proper simulations at four different r/R nl values (panels).Ratios are plotted relative to the fiducial simulation with η = 0 .
15 (i.e. the middle point). The behavior is nearly independentof epoch in r/R nl units, but we show the mean and standard deviation over epochs for completeness. REFERENCES