{'fraction_non_alphanumeric': 0.030534160248074423, 'fraction_numerical': 0.006977093127938381, 'mean_word_length': 4.825903263403263, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Relativistic Interpretation of Bias in Newtonian Simulations

Observables of cosmic structures are usually not the underlying matter field but biased tracers of matter, such as galaxies or halos. We show how the bias found in Newtonian N-body simulations can be interpreted in terms of the weak-field limit of General Relativity (GR). For this we employ standard Newtonian simulations of cold dark matter and incorporate GR/radiation via a weak-field dictionary that we have recently developed. We find that even when a simple local biasing scheme is employed in the Newtonian simulation, the relativistic bias becomes inherently scale-dependent due to the presence of radiation and GR corrections. This scale-dependence could be in principle observed on large scales in upcoming surveys. As a working example, we apply our methodology to Newtonian simulations for the spherical collapse and recover permille-level agreement between the approaches for extracting the relativistic bias on all considered scales.


Introduction
The bias of structures (halos, galaxies) compared to the underlying dark matter distribution is a fundamental ingredient for the analysis of large-scale structure surveys, and has been extensively studied in both Newtonian [1][2][3][4][5][6][7] and relativistic [8][9][10][11][12] frameworks. However, it is not only a requisite to analyse surveys, but measuring the luminous galaxies w.r.t. the underlying dark matter density offers a second observable containing precious information from the large-scale structure. This is especially true when observing two (or more) populations that are biased independently [13], and because the bias is sensitive to the primordial non-Gaussianity [14], thereby opening a window into the early Universe.
On sufficiently small cosmological scales a Newtonian analysis is expected to be accurate, at least in the most simple model of a universe filled only with cold dark matter, however the next generation of galaxy surveys will probe large volumes, comparable to the cosmological horizon, where a simplified Newtonian analysis is no longer accurate. In addition, the early Universe was dominated by radiation, including photons and massive neutrinos. The latter are also expected to play an important role in structure formation. All these contributions are inherently relativistic and thus beyond the scope of a Newtonian framework.
During the past years, several approaches to simulate the relativistic evolution of the Universe have been developed. In [15,16], employing the relativistic code gevolution, Nbody simulations in the weak-field limit of General Relativity (GR) were performed (see [17] for the inclusion of massive neutrinos). Subsequently, photons and massless neutrinos, i.e., purely relativistic species that are of particular importance on large cosmological scales (and rather at early times), have been incorporated essentially by three different approaches. One of these approaches was established by a modified Newtonian N-body code [18], dubbed cosira, that makes use of a linear GR dictionary, based on the so-called N-body gauge [19]. Another approach was established by a modification of the aforementioned relativistic code gevolution [20]. Lastly, the authors of [21,22] introduced the Newtonian motion (Nm) gauge approach, originally based on a linear GR dictionary, which allows the interpretation of an unmodified Newtonian simulation on a perturbed relativistic space-time. Shortly later, the Nm gauge has been generalised in [23], employing a GR dictionary that is valid within the weak-field limit of GR, and furthermore in [24] including massive neutrinos. In the present paper we employ the non-linear Nm approach and study how the bias in GR can be connected to the bias observed in a Newtonian N-body simulation.
So far, there exists no weak-field description for the bias nor a bias description that incorporates the effect of radiation -both issues that we aim to address in the present paper. In [8][9][10][11][12], in the absence of radiation, the relativistic bias has been investigated to linear order in cosmological perturbation theory (CPT) [25][26][27]. It has been found that the linear bias receives a GR correction stemming from a long-wave perturbation (itself directly related to the comoving curvature perturbation), a result that we also recover in the radiation-free and linear limit. In [28] the relativistic bias has been investigated in a gauge-invariant formulism at second order in CPT, accounting for both volume and source effects. Recently, some confusion exists in the literature about appropriate gauge choices to deal with bias beyond the linear order. In particular, the synchronous-comoving gauge choice (see e.g. [29]) was discussed in [30] who concluded that this gauge should be avoided for bias definitions, although, to our opinion, for the incorrect reasons. While [30] claimed that mass conservation is violated in the synchronous-comoving gauge at second order, it was soon shown in [31,32] that the term that seemingly violates mass conservation is the consequence of employing a Lagrangian-coordinates approach for the relativistic bias. Mass conservation is a fundamental principle in relativistic and Newtonian fluid mechanics and holds irrespective of the chosen fluid frame, and therefore we disagree with the reasoning of [30].
Although a second-order bias model in the synchronous-comoving gauge may lead to a consistent description on the largest scales, however, for the present purpose of incorporating effects from radiation and the non-linear collapse, this gauge choice should be avoided, essentially for two reasons. The first reason for that is that the gauge conditions of the synchronous-comoving gauge become inconsistent when fluids with pressure (e.g., [27]; or non-zero vorticity [33]) are included. Secondly, bias relates to the event of matter collapse where non-linearities in the density (and velocity) inevitably become large. In a weak-field description this process can be modelled accurately as long as the underlying metric perturbations remain small. The synchronous-comoving gauge however is not a suitable candidate to perform such a weak-field limit as in that gauge the non-linear collapse is imprinted in generally large metric perturbations (related to the relativistic Lagrangian displacement field, cf. [31]). For similar reasons, also comoving-orthogonal gauges (e.g., the total matter gauge) should be avoided for non-linear bias descriptions, since in such gauges the (generally large) velocity is proportional to the space-time component of the metric. Therefore, to include the leading non-linearities in a relativistic bias description, the question of gauge choice becomes important. In the present paper, we discuss certain gauge choices that in particular do not violate the smallness assumptions of the metric perturbations, and additionally allow for large non-linearities in the density and velocity variables. This comes with the advantage that our approach remains accurate until much smaller scales, but sacrifices some precision on the largest scales (where however a linear treatment may suf-fice). Furthermore our approach is not limited to a simple dark matter cosmology and allows us to include the impact of photons or neutrinos on the relativistic bias.
The paper is organised as follows. In the following we introduce our notations and conventions. In section 3 we provide a summary of the used gauge choices, related to the recently developed Nm framework, as well as briefly discuss the resulting relativistic equations (section 3.4) and outline our used Eulerian bias description (section 3.5). In section 4 we discuss the computation of the relativistic bias from Newtonian simulations. Further results on the bias, in particular the study of the radiation-free (late-time) limit and the case of massive neutrinos, are discussed in section 5. We confront our theoretical predictions with a relativistic numerical simulation in section 6, and conclude in 7.
Notation: Newtonian quantities, such as the density and velocity (obtained from the Newtonian simulation), are labelled with the index "N", whereas relativistic quantities come with no special index. Biased quantities have the index "X". We set the speed of light to unity. Most of our equations are formulated in real space, except some of the results in section 4 where we occasionally outline results in Fourier space -noted in the text and through the explicit dependence on the wave vector |k| = k.

Metric convention in our weak-field description
The metric line element in a yet unspecified gauge has the following metric coefficients of the scalar type Here, a is the cosmic scale factor, δ ij the Kronecker delta, we make use of the conformal time τ , and we have defined the normalised gradient operator∇ i ≡ −(−∇ 2 ) −1/2 ∇ i which translates into −ik i in Fourier space, wherek i ≡ k i /|k|. By employing the normalised gradient operator, we guarantee that all scalar perturbations in the above metric are of the same order in our expansion scheme. We work in the weak-field limit of general relativity, amounting to a double expansion scheme in metric perturbations and spatial derivatives. The weak-field approach encompasses the largest cosmological scales where linear perturbation theory delivers an accurate description, but also incorporates the leading-order contributions from smaller scales where non-linear effects are important. We only demand the smallness of the metric perturbations while other perturbations, especially the matter density, may become arbitrarily large. See [23] for a more thorough discussion about our used weak-field definition which is in agreement with the employed weak-field scheme of [34][35][36][37].

Recap of Newtonian motion (Nm) approach
Let us briefly outline the Nm framework [21][22][23][24] which we use as the weak-field dictionary for the Newtonian simulation. The basic idea of the Nm approach is to enable conventional Newtonian (N-body) simulations to obtain a relativistic description that includes the effect of radiation on cosmic structure formation. Massless radiation can be described by linear cosmological perturbation theory accurately at all times and scales relevant for structure formation, and consequently the first Nm dictionary employed only a linear dictionary [21]. In this approach the linear radiation source provides feedback to the fully non-linear matter field, that itself is evolved within a Newtonian simulation. Subsequently this dictionary was generalised to the weak-field limit. The weak-field Nm framework essentially provides a unified method for resolving the largest scales where linear perturbation theory provides an accurate description as well as the non-linear scales where it includes the leading non-linear contributions.
In the present paper we employ the weak-field Nm formulation, relevant for studying the small scales that lead to the formation of structures. A consistent weak-field description requires that the gauge conditions do not violate the underlying smallness assumptions of the metric potentials. For such gauges the relativistic evolution equations for matter include several relativistic corrections that are not present in a Newtonian description. For example, the relativistic momentum conservation is coupled to the radiation energy density and to the metric perturbations -relativistic effects that are not present in the Newtonian momentum conservation (the Euler equation).
We use the spatial gauge degree of freedom in GR to select so-called Nm gauges, that, by definition, imply an equation for the relativistic momentum conservation that precisely matches its Newtonian counterpart. This way, the GR corrections are embedded in the underlying space-time that can be accounted for as a post-processing of the Newtonian simulation simply by using a weak-field Nm dictionary. For the relativistic equations for matter in the presently employed Nm gauge, see section 3.4.

Used temporal gauge choice
The Nm idea is related to the spatial gauge fixing, described above. For any temporal gauge choice that does not violate the weak-field smallness conditions, there exists a unique Nm gauge. One possible temporal gauge choice that satisfies this criteria is the Poisson temporal gauge condition (see e.g., [38,39]) which we will employ in the present paper, with the operator K = (−∇ 2 ) 1/2 which corresponds to the magnitude k = |k| in Fourier space. As a consequence of this temporal gauge fixing, the density in our Nm gauge is identical to the density in the Poisson gauge.

The N-boisson gauge
Besides of the Nm gauge we make use of another gauge in the present paper, the N-boisson gauge. We define that gauge through the temporal gauge condition of the Poisson gauge, is the gauge-invariant comoving curvature perturbation. This is a special weakfield gauge in which conventional simulations with CDM and Λ can be brought in agreement with GR -provided that radiation plays no significant role for the gravitational dynamics anymore. For a standard ΛCDM cosmology, the latter is the case for sufficiently late times (z < 50). In a sense, that gauge behaves as a so-called N-body gauge [19] which implies that for vanishing radiation, the relativistic momentum conservation of matter agrees with the Newtonian Euler equation. In the spirit of the temporal gauge condition, the N-boisson gauge has also advantageous properties in the weak-field sense, i.e., all metric perturbations remain naturally small.
The N-boisson gauge has a different spatial gauge condition as the Poisson gauge, and because of that the analysis in the N-boisson gauge differs slightly from the one of [40] in the Poisson gauge. Nonetheless, their analysis can be connected to the N-boisson gauge, since in their linear GR dictionary, after the Newtonian simulation has been terminated, N-body particles are displaced to the spatial position of the N-boisson gauge (see also the related discussions in [20]).

Relativistic equations of motion in the Nm and N-boisson gauge
As outlined above, we employ the Nm and N-boisson gauges that make use of an identical temporal coordinate, whilst the spatial coordinates are related via a simple spatial gauge transformation ( ). To outline the difference between these two gauges it is instructive to fix for the moment only the temporal gauge condition. In such a not yet completely fixed gauge the relativistic momentum conservation of CDM reads (vector perturbations are suppressed for simplicity), where v i = ∇ i v cdm , and Σ cdm is the anisotropic stress resulting from overlapping matter streams. The above equation is simply relating the gravitational acceleration of fluid trajectories (the l.h.s.) to the Hubble drag term, the anisotropic stress and certain metric perturbations (all r.h.s.). This relativistic result (derived in [23]) should be contrasted with the Newtonian momentum conservation which reads where the index "N" denotes Newtonian quantities (as determined through Newtonian simulations), and Φ N is the cosmological potential subject to the usual Newtonian Poisson equation. At this stage, having eq. (3.2) as a starting point gives us two options that allows us to find relativistic trajectories that are formally identical with the Newtonian trajectories: 1. We fix the spatial gauge by demanding that the relativistic corrections on the r.h.s. of (3.2) are exactly identical with the Newtonian potential, i.e., leading to the relativistic momentum conservation in the Nm gauge which is formally identical with the Newtonian momentum conservation, even at times when radiation can not be neglected.

We fix the spatial gauge according to
completing the gauge fixing for the N-boisson gauge, which also results in having a Poisson equation in Newtonian form (see below). The resulting relativistic momentum conservation then still contains a relativistic term ∼ (∂ τ + H) K −1∇ iζ , however, crucially, that relativistic term vanishes when the effect of radiation can be neglected. This is precisely the case at sufficiently late times (z < 50), which thus makes the N-boisson gauge a suitable gauge for analyses at late times.
As regards to the mass conservation, since the temporal gauge conditions in the discussed two gauges are identical, mass conservation for both the Nm and N-boisson gauge reads where we have defined a new density, the counting density, with δ = (ρ −ρ)/ρ being the actual relativistic density contrast. The counting density will be useful later when relating relativistic densities to the density as counted in Newtonian simulations. Finally, note that the Nm spatial gauge condition features a residual gauge degree of freedom, which suitably can be linked to the initial conditions (ICs) that are used to initialise the Newtonian simulation; see section 4.1. in [21] for a thorough discussion on ICs. In the present paper, we fix the remaining gauge freedom by identifying the spatial gauge at the initial time with the N-boisson gauge.

Eulerian bias description
In the present paper we investigate a dictionary that links results from Newtonian matter simulations to a relativistic bias description. That dictionary also requires the specification of a Newtonian bias scheme that relates the simulated Newtonian density to the tracer density obtained in that simulation. We employ a local Eulerian bias scheme, defined in real space via [7] δ N X = b 1 δ N + 8) where δ N X is the density of the tracer density X (e.g., halos, galaxies), b 1 , b 2 and b s are constant Newtonian bias coefficients, δ N the fully non-linear density as measured from Newtonian simulations, and s N = 2∂ i ∂ j Φ N /(3H 2 ) − δ ij δ N /3 a function that extracts the tidal effects from the non-linear density. Here, H is the conformal Hubble parameter, Φ N the Newtonian cosmological potential satisfying the usual Poisson equation, and ∂ i the Newtonian spatial derivative w.r.t. Eulerian component x i .
we seek to find the relativistic counterpart of δ N X , the relativistic tracer density that we call δ X in the following.
First, we can embed a Newtonian simulation into General Relativity by using the Nm gauge. This embedding is uniquely defined by the ICs of the corresponding Newtonian simulation. This means that the measured bias parameters from a Newtonian simulation are linked to exactly one Newtonian motion gauge. Employing the weak-field dictionary the relation between the Newtonian density, δ N , and the relativistic density, δ, is (Nm and N-boisson gauge, with radiation) (4.1) which evidently agrees with the aforementioned counting density defined in eq. (3.7). While the pure matter densities are related to their GR counterparts by this very simple dictionary, a similar relation for the tracer density δ X does not exist in general. This question is linked to understanding the process of tracer formation from the underlying dark matter distribution including the impact of relativistic corrections and may depend significantly on the chosen gauge.
Using a Newtonian motion gauge leads to small deformations of the coordinate system beyond the homogenous expanding space-time assumed in the Newtonian simulations. While this is an important correction on large scales close to the horizon, we find that in all Nm gauges that are weak-field compatible (cf. the discussion in [23]), the physics on non-linear scales is not affected by relativistic corrections. For example while the distance between two distant objects in GR may be significantly different from the distance found in a Newtonian simulation, local distance measurements remain almost unaffected.
From this it follows that a halo finder that is taking into account all the information present in GR will find exactly the same halos as the Newtonian one. The reason being that the process of defining a halo is completely local and the relativistic interpretation may change the distances between halos, but not the local physics. This statement remains true as long as H T is small compared to the non-linear matter densities on scales that are relevant for the formation of the tracers.
We thus can employ the simple geometrical relation, taking into account the volume deformation (i.e., the deformation of the spatial 3-volume due to relativistic diagonal corrections) that alters the mass density according to δ N X = δ X + 3H L .
(Nm and N-boisson gauge, with radiation) (4. 2) The last equation is one of our main results in the present paper; it provides the link between the Newtonian and relativistic tracer density. This equation states that the tracer density is only changed by the large-scale GR corrections that move around individual halos compared to their Newtonian positions, but the GR dictionary does not change whether a local structure is identified as a halo/galaxy. Since we consider a local Newtonian bias scheme, which amounts at the leading order to δ N X = b 1 δ N , we find the following simple relation for the relativistic bias , (Nm and N-boisson gauge, with radiation) (4.3) where we have used our GR-dictionary (4.1) and the above Eq. (4.2). Observe that an additional bias proportional to the local volume perturbation H L appears on top of the Newtonian bias b 1 that is now multiplied with the relativistic density. This additional term describes that in GR, not only high-density regions can collapse, but also regions that start on a large volume perturbation but with a somewhat smaller initial matter overdensity. This is directly opposed to the Newtonian case where a collapse is uniquely related to the initial overdensity alone. This additional correction from GR to the relativistic bias is multiplied by (b 1 − 1) stating that it results from a source that does not follow the dark matter particles or halos; it naturally has to vanish for tracers with bias b 1 = 1.
A key signature of the relativistic bias is its generic scale dependence. Transforming for the next few equations into Fourier space, we can recast the GR bias into a bias based on the relativistic density δ alone: where T H L and T δ are respectively the transfer functions of H L and δ. The relativistic bias is necessarily scale dependent due to the different evolution of the density and volume perturbation.
Including the higher-order terms for the relativistic bias is straightforward. Using the GR dictionary and the same steps as outlined above, we find in real space with γ ≡ −(∂ τ + H)K −2Ḣ T + 8πGa 2 K −2 Σ and the Bardeen potential Φ, where it should be noted that, due to the used spatial gauge fixing, in the last term Φ+γ is nothing but Φ N . Our analysis is consistent to leading order with weak-field relativity, however some of the terms in (4.6) are subdominant and may be neglected. We find that H L is weak-field suppressed compared to δ while γ is suppressed relative to Φ. This implies that we may neglect all relativistic corrections in the higher-order bias. We thus have to a good approximation with s ij = (2∂ i ∂ j Φ/(3H 2 ) − δ ij δ/3). The same necessarily applies not only to second order, but to all higher-order corrections in the weak-field limit. It should be noted that the bias in general relativity does depend on the chosen gauge and this is the bias in the Newtonian motion gauge corresponding to the simulation that was used to determine b 1 , b 2 and b s . However, the Newtonian motion gauge shares the temporal gauge choice of the Poisson gauge and therefore this relation may be connected to the Poisson gauge bias with very little effort. The densities and the potential Φ remain identical, while only H L has to be replaced with Φ − H T /3.

Further results on the relativistic bias
Most studies of bias in general relativity employ a simplified Universe, dominated by matter plus a cosmological constant. For a standard ΛCDM Universe, this simplification is valid at sufficiently late times. This makes for a particularly interesting and simple special case that we discuss in the following section. Then, in section 5.2 we switch on radiation and provide the full analysis.

Bias in the radiation-free (late-time) limit
In the case of vanishing radiation there exists a remarkably simple Newtonian motion gauge, i.e., the N-boisson gauge (see section 3.3) with H T = 3ζ and ζ the comoving curvature perturbation, for which the densities are related via (N-boisson gauge, no radiation) (5.1) where Φ and ζ can be easily obtained by using standard linear Boltzmann codes. By contrast, to derive the density in the case where radiation can not be neglected, one is forced to calculate H T which is computationally more demanding (cf. [23]). Traditionally, the relativistic bias has been investigated in the relativistic counterpart of a Lagrangian-coordinates approach [8,11]. For vanishing vorticity (cf. [33]) that Lagrangian approach is established in the synchronous-comoving gauge. To linear order the relativistic Lagrangian bias is (synchronous-comoving gauge, no radiation, linear order) (5.2) where b L 1 is the relativistic bias in Lagrangian coordinates. Now we show how this standard linear equation can be related to our findings (which however hold beyond linear theory). For this we employ the Nm dictionary to connect the simulation density to the N-boisson gauge density, which, as discussed above, is identical with the density in the Poisson gauge (Poisson and N-boisson gauge, no radiation) (5.3) There are two simple but important observations to make the connection of the last result to eq. (5.2). First, observe that the 3(Φ − ζ) = HK −1 v is nothing but the temporal gauge generator that connects the Poisson (or N-boisson) gauge to the synchronous-comoving gauge. Secondly, from Newtonian bias arguments [6] it is well known that the linear Eulerian bias is related to the linear Lagrangian bias via a relation that makes much sense also in the present case. Indeed, the quantity 3(Φ − ζ) physically amounts to the local volume deformation, a Lagrangian quantity; explaining why the bias factor attached to 3(Φ − ζ) is a Lagrangian one. Furthermore, this directly shows that the two relations are compatible and that the simulation density δ N = δ + 3(Φ − ζ) is in fact identical to the synchronous-comoving gauge density (again, vanishing radiation and linear treatment is assumed). We have thus established a simple relation how the linear relativistic bias can be interpreted in other gauges. Before considering the case with non-vanishing radiation in the following section, let us briefly comment on the interpretation beyond the linear approximation. There has been significant work on bias to second order in perturbation theory (see [30,32]), but we work to leading order in weak-field gravity, where densities and velocities may become large while the metric potentials remain small. Our work thus is valid up to the very small scales, while our approach does not include second-order metric corrections on the large scales.
Within the weak-field limit, the synchronous-comoving gauge is no longer a viable choice, since the metric perturbations do not remain small and in fact are responsible for the nonlinear the growth of structures. The Nm gauge (or the Poisson gauge), on the other hand, remains well defined and can be used to describe the bias consistently to weak-field precision. For this reason our formulation describes the physics accurately on all relevant scales including non-linearities, while a description in synchronous-comoving gauge should be avoided on the non-linear scales due to its incompatibility with the weak-field approach.

The impact of relativistic species on bias
As evident by a host of astrophysical observations, the early Universe is radiation dominated, and, depending on the initial redshift considered in the simulation, radiation might introduce a non-negligible correction. In this case we can no longer use the N-boisson gauge, but have to use a dynamically evolving Nm gauge, leading to the bias (obtained from replacing ζ → H T in (5.3)): (Nm gauge, with radiation) (5.5) where H T is related to the spatial gauge transformation connecting the Nm gauge to the Poisson gauge. The simple relation to the synchronous gauge that we had found in the pure matter case is lost in the presence of radiation. The reason for this is because the synchronous-comoving gauge is no longer comoving to the matter particles in the presence of additional species that possess a pressure. Calculations of the bias based on the synchronouscomving gauge therefore can no longer be applied in more complex cosmologies, while our method (which is closely related to the Poisson gauge) can be adapted.
In the case of massive neutrinos, we found in a recent paper [24] that neutrinos have a non-trivial impact on structure formation on the smaller scales (also within the Nm gauge framework). The above assumption 4.2 for the bias can thus no longer be used since they neglect the non-trivial impact neutrinos have on the small-scale structure formation (especially for large neutrino masses). A solution to this is obtained by using the specifically designed backwards approach presented in [24]. In the following we briefly summarise the ingredients for that approach, and refer to the aforementioned papers for further details.
The so-called backscaling approach to N-body simulations makes use of ICs at time z = 0 which are then appropriately rescaled by the linear-growth function to the desired initialisation redshift. This way, not the actual ICs of the universe are prescribed, but instead of a fictitious universe that has the same radiation content as the Universe at redshift z = 0. This backscaling approach is actually the standard to initialise N-body simulations, since these simulations generally do not evolve for radiation perturbations. In [22] it has been shown, that N-body simulations from such backscaled ICs can be self-consistently embedded in a relativistic space-time.
For massive neutrinos the situation is more complex and in [24] a backwards approach for neutrinos was described, that coincides with backscaling only in the limit of massless neutrinos. It allows the absorption of small scale effects directly into the ICs and thus treat them at the level of the Newtonian simulation and not as a modification of the spacetime. This in turn can be used to solve the problem mentioned above. By employing a backwards method, the complex effects of neutrinos are automatically included in the Newtonian simulation and already present in the 'Newtonian' bias obtained from such a simulation. We are then able to relate the bias from such simulations to the relativistic Poisson-gauge bias using the simple relation (5.3).

Numerical Simulations
Our analysis connects the non-linear bias found in Newtonian simulations to the relativistic weak-field bias in Poisson gauge. However it also holds for other types of Newtonian simulations. In order to show the impact of relativistic corrections, we perform a suite of simulations for the spherical collapse of a top-hat overdensity, both using Newtonian equations and using general relativity in the Poisson gauge.
We consider an isolated spherical mass shell characterised by a small-scale overdensity δ S in a background Universe perturbed by a long mode δ L . By defining the mass with respect to the shell's physical volume M S = V R d 3 R δρ S we receive the shell mass, which sources our local gravitational potential. The geodesic equation for the physical radius of the shell R reads in the weak-field limitR where an overdot denotes a partial derivative w.r.t. conformal time. The last term describes the shell mass driving the collapse, while the Hubble dragä a describes the impact of the background evolution on the shell. The contributionsΦ + 2ȧ aΦ is the relativistic acceleration of the shell including the impact of the long mode in GR. For our simulations in Newtonian gravity we neglect these terms, while we keep them otherwise.
To compute the bias we fix the initial time z ini = 100 and the time of collapse z coll = 3, and solve the equation for a structure mass of M = 10 13 M sun for varying ICs δ ini S ∈ [0.14, 0.15] and δ ini L ∈ [0.01, 0.02]. From this we obtain the correlation between the small and large mode, which we can translate to the bias. For more details about the used method see [41].
Our results are shown in figure 1. We find that the Newtonian case leads to a nearly scale-independant bias as expected. Employing the full Poisson gauge dynamics however leads to a significant modification of the bias on scales larger than the matter-radiation equality scale. We then compare this with the relativistic bias obtained from interpreting the Newtonian one using Eq. (5.5). The Newtonian motion gauge interpretation of the Newtonian bias agrees with the relativistic Poisson gauge bias at the permille level on all analysed scales. The predicted scale dependance of the relativistic bias is clearly visible and is matched by the relativistic simulation perfectly.
We further may investigate the impact of relativistic species on the bias using our simulations. Instead of the N-boisson gauge that should be used in a pure matter Universe, we may employ the dynamical Newtonian motion gauge corresponding to the initialisation time of z = 100 in Eq. (5.5). We find that the effect is negligible on the relativistic bias; at least for initial times of redshift 100. But we may still confirm that our method does compute radiation effects consistently. Since the large-scale relativistic bias becomes very close to unity, small corrections can no longer be seen. However we may also reconstruct the Newtonian bias from the relativistic one using the inversion of Eq. (4.5) that returns an almost scale-invariant result that can be compared against the Newtonian simulation. The outcome is illustrated in figure 2. On the small scales both methods are consistent, but on the large scales we only find a perfect agreement between the reconstructed Newtonian bias and the numerical simulation when including also the impact of radiation on the Nm gauge, leading to a perfectly scale-independent Newtonian bias.

Conclusions
We have shown how the non-linear bias of Newtonian N-body simulations can be interpreted in terms of weak-field GR, being valid on all scales relevant for the formation of structures in the Universe. Our prescription allows the extraction of the bias in the Poisson gauge, which is a gauge that remains well defined in the weak-field limit of GR.
In contrast to previous studies of bias in GR, we do not employ perturbation theory, but a weak-field description that remains accurate on smaller scales, where non-linear effects become important. When restricted to linear theory and to a universe only filled with matter, our results agree with the standard ones from the literature. In particular, within that limit, we recover the well-known result that the simulation bias can be directly related to the relativistic bias in the synchronous-comoving gauge. By going beyond this simplified setup, by both including non-linear effects as well as the effects from relativistic species, we have shown that the GR dictionary has to be slightly updated, thus deviating from the standard results in the synchronous gauge. Apart from the employed Nm gauges, we also find that the Poisson gauge is a very suitable gauge for investigating the non-linear bias.
For the relativistic bias we detect a scale-dependency which kicks in on large scales, ranging from the scale of matter-radiation equality up to super-horizon scales. Essentially, this scale-dependency is the consequence of the matter collapse being not only governed by the matter overdensities but also affected by the metric perturbations. In particular, the Figure 2. Shown is the reconstruction of the Newtonian bias starting from the relativistic one. In the top plot the calculation using the full H L is displayed in green while the magenta dotted line depicts the pure matter computation H L = Φ − ζ. The bottom plot illustrates the ratio between the transfer functions of the two cases at z ini . Note that the axes of this plot span a much smaller range than in the previous one. The Newtonian bias is scale independent within the numerical error of our simulation. metric perturbations that are most relevant for the bias are the Bardeen potential Φ, which deforms the local spatial volume element (w.r.t. to a Euclidean volume element), as well as the comoving curvature perturbation ζ which acts as a long wave-length perturbation for the collapse.
With the present methodology, the relativistic bias can be extracted from all standard Newtonian N-body simulations. In particular, Newtonian simulations can be interpreted in terms of weak-field GR within two orthogonal approaches, called the forwards and backscaling approach. In the forward approach, the actual ICs for the Newtonian simulation are taken from the initialization redshift, whereas in the backscaling approach (which is very common in the literature) one initializes the Newtonian simulation by rescaling the linear final power spectrum back to the initialization redshift. Within our approach, we are able to provide a relativistic interpretation of the bias for both the forward and the backward approach. This is particularly simple in the backscaling approach as the complications are already present in the employed initial conditions and the simple dictionary may be used (see section 5.2).
Our weak-field approach allows not only the inclusion of non-linear effects but is also able to incorporate the effects of relativistic species. We have already seen that the relativistic large-scale feature affects scales beyond the matter-radiation equality scale, linking it to the presence of radiation. But this feature is mediated via the metric potentials 'remembering' that the Universe underwent a phase of radiation domination before. In addition there is the direct impact that the remaining photon content has on a matter collapse. This effect was in fact computed in [41] where a small but visible correction was found. Using the Nm gauge approach we confirm a feature on the largest scales, but we have shown that this particular feature is extremely small in the relativistic bias, and thus may safely be ignored. Massive neutrinos, by contrast, have very likely a non-negligible impact on the bias. We gave the first steps to incorporate the effects of massive neutrinos, and argued that a backwards Nm approach applied to an actual Newtonian simulation with neutrinos should deliver an accurate description for the bias.
As a working example, we have applied our weak-field approach to Newtonian simulations for the matter collapse of a spherical-top hat perturbation. As expected, a purely Newtonian analysis does not yield to any scale-dependent features for the bias, whereas a GR description indeed reveals a significant modification of the bias on large scales (see Figs. 1 and 2). By employing our weak-field description, we then recover permille-level agreement between the approaches on all considered scales.
Having obtained the framework to interpret the relativistic bias from Newtonian simulations, we plan to apply the techniques to actual N-body simulations (going beyond the spherical collapse). A particularly interesting topic would be to explore in depth the bias in the presence of massive neutrinos, and to investigate whether knowledge of the neutrino mass scale could be extracted from future large-scale structure surveys. Such avenues would require high-resolution studies, obtained either from (relativistic) N-body simulations that actively include the evolution of massive neutrinos [17], or from standard N-body simulations that incorporate the effects of massive neutrinos via the Nm-gauge framework [24].