{'fraction_non_alphanumeric': 0.05378438543488631, 'fraction_numerical': 0.020273493401176658, 'mean_word_length': 4.179534692196829, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 2, '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}}
Perturbed Newtonian description of the Lema\^itre model with non-negligible pressure

We study the validity of the Newtonian description of cosmological perturbations using the Lemaitre model, an exact spherically symmetric solution of Einstein's equation. This problem has been investigated in the past for the case of a dust fluid. Here, we extend the previous analysis to the more general case of a fluid with non-negligible pressure, and, for the numerical examples, we consider the case of radiation (P=\rho/3). We find that, even when the density contrast has a nonlinear amplitude, the Newtonian description of the cosmological perturbations using the gravitational potential \psi and the curvature potential \phi is valid as long as we consider sub-horizon inhomogeneities. However, the relation \psi+\phi={\cal O}(\phi^2), which holds for the case of a dust fluid, is not valid for a relativistic fluid and effective anisotropic stress is generated. This demonstrates the usefulness of the Lemaitre model which allows us to study in an exact nonlinear fashion the onset of anisotropic stress in fluids with non-negligible pressure. We show that this happens when the characteristic scale of the inhomogeneity is smaller than the sound horizon and that the deviation is caused by the nonlinear effect of the fluid's fast motion. We also find that \psi+\phi= \max[{\cal O}(\phi^2),{\cal O}(c_s^2\phi \, \delta)] for an inhomogeneity with density contrast \delta whose characteristic scale is smaller than the sound horizon, unless w is close to -1, where w and c_s are the equation of state parameter and the sound speed of the fluid, respectively. On the other hand, we expect \psi+\phi={\cal O}(\phi^2) to hold for an inhomogeneity whose characteristic scale is larger than the sound horizon, unless the amplitude of the inhomogeneity is large and w is close to -1.


Introduction
The Newtonian description of cosmological perturbations plays a key role in studies of structure formation as it often allows an intuitive understanding of the dynamics. The availability of increasingly precise cosmological observations requires accurate theoretical tools in order to describe the evolution of cosmological perturbations. In this context, it might be useful to reconsider the properties of the perturbed Newtonian description of cosmological inhomogeneities. In general, the Newtonian description of general relativity is based on the slow-motion approximation. Therefore, a system with non-negligible pressure could be interesting in order to check the validity of the Newtonian description.
The Lemaître model [1] is a spherically symmetric solution of Einstein's equation which can be used to obtain the exact dynamical evolution of a perfect fluid. Therefore, the Lemaître model provides us with the possibility of testing how well the exact solution can be reproduced in terms of the perturbed Newtonian description. This problem has been considered in Refs. [2][3][4][5] for the case of a dust-dominated universe (see also [6][7][8]), which is usually named the Lemaître-Tolman-Bondi (LTB) model [1,9,10]. 1 Here, we aim at extending previous work to the more general case of a fluid with non-negligible pressure. Contrary to its pressureless LTB limit, the Lemaître model cannot be solved analytically. Therefore, we integrate the relevant equations numerically as done, for instance, in [12][13][14][15][16][17][18][19].
We particularize our analysis to the case of radiation (P = ρ/3) and study the corresponding evolution for a set of different initial conditions. We then study the validity of the description in terms of the perturbed Newtonian Friedmann-Lemaître-Robertson-Walker (FLRW) metric by performing the coordinate transformation numerically for the set of different initial conditions previously considered.
This paper is organized as follows: In Section 2 we derive the basic equations for the Lemaître model in a cosmological setup and in Section 3, for the case of a relativistic fluid, we solve numerically the equations and show the dynamics. In Section 4, we perform an exact coordinate transformation from the Lemaître model to the perturbed FLRW metric. Our results demonstrate the validity of the cosmological Newtonian description, which assumes |ψ|, |φ| ≪ 1, where ψ and φ are the gravitational potential and the curvature potential in the Newtonian gauge, respectively. Then we show that φ + ψ = O(φ 2 ) -which holds in the case of a dust fluid -is violated for a relativistic fluid with inhomogeneities at sub-sound horizon scales. In Section 5, we show that the nonlinear effect of the fluid's fast motion is responsible for this failure, which is clarified by considering the second-order perturbations in spatially conformally flat spacetime, as a generalization of the perturbed FLRW metric. This is an example of perturbations in a fluid with non-negligible pressure that effectively give rise to an anisotropic stress, which is an expected second-order effect [see, e.g., [20][21][22][23][24][25]. Therefore, by studying the Lemaître model we are able to follow the generation of anisotropic stress in an exact nonlinear way. Section 6 is devoted to summary and conclusions. In Appendix A we develop the perturbation theory for the Lemaître model. Throughout this paper we adopt c = 1, where c is the velocity of light.

The Lemaître model
In this Section we introduce the Lemaître model [see, e.g., 16] -a spherically symmetric solution of Einstein's equation -for the case of a perfect fluid with non-negligible pressure. 2 Using comoving coordinates, it is customary to write the line element as: where dΩ 2 = dθ 2 + sin 2 θdφ 2 and e ν is the lapse function, which cannot be gauged away when pressure gradients are present. Einstein's equation gives then: where a dot and a prime denote differentiation with respect to t and r, respectively, and ρ and P are the energy density and pressure of the perfect fluid. With the use of (2.3), equations (2.2) and (2.4) yield: respectively, where we defined the effective gravitating total mass M : As is clear from (2.6), M is related to the local density ρ through the Euclidean volume element. Consequently, M does not coincide with the invariant mass. We call M the effective gravitating mass because it enters the generalized Friedmann equation (2.22). See Refs. [16,18,27,28] for more details.
In the comoving coordinate we are using, the four-velocity of the perfect fluid is v µ = (e −ν , 0, 0, 0). The equations of motion of the perfect fluid are then: which are integrated as: , (2.12) where in the last step the equation of state P = wρ with a constant parameter w was assumed. We assume that the sound speed is c s = √ w, unless explicitly noted otherwise. In the above equations, ρ i (r) = ρ(t i , r), R i (r) = R(t i , r) and λ i (r) = λ(t i , r) are evaluated at the initial time t i , and ν 0 (t) = ν(t, r 0 ) is an arbitrary function which is determined once we fix the residual temporal gauge. By combining (2.9), (2.10) and (2.3) we obtaiṅ which with (2.10) yieldsλ (2.14) Then by integrating (2.3) one obtains the radial scale factor e λ as a function of the angular scale factor R and the lapse function ν: where E(r) is related to λ i (r) by e λ i (r) = R ′ i (r)/ 1 + 2E(r). Finally, the combination G θ θ − G r r from (2.4) and (2.5), which have not been used so far as a consequence of the Bianchi identity, gives: We will use this equation in order to check our numerical computations (see text below (3.10)).

Redefinition of scale-factor functions
We will find it useful to redefine the functions λ and R that describe the radial and angular scale factors by means of the functions a(t, r) and E(t, r) according to: R(t, r) = ra(t, r), e λ(t,r) = a(t, r)e E(t,r) . (2.17) Thanks to this redefinition the line element now reads: To our knowledge, this parametrization has not been used before. Equations (2.6) and (2.7) become: where the effective gravitating mass is now: From the previous equation one obtains: which can be used for the time evolution of a(t, r). Equation (2.12) is rephrased as where a i (r) = a(t i , r) and E i (r) = λ i (r) − ln a i (r). Equations (2.13)-(2.15) are rewritten as:

Initial and boundary conditions
We will now discuss the initial conditions necessary in order to integrate numerically [see, e.g., 16,19] the system of coupled differential equations presented in the previous section.
We adopt a gaussian initial density contrast: where the parameter A specifies the amplitude of the initial inhomogeneity, L is its comoving characteristic size,ρ i = 3H 2 i /8πG is the background (critical) density, H =ȧ/ā is the background Hubble parameter, andā(t) = a(t, r = ∞) is the background scale factor. 4 As a consequence of the chosen profile, the background FLRW model is recovered only at infinity -that is,ρ i = ρ i (t i , r = ∞) -rather than at a finite radius [see 29,30]. Also, as is clear from (2.26), if E(r = ∞) = 0, then the background model is spatially flat. We adopt t i = 0 and time is given in units of the Hubble time H −1 i , while the radial coordinate r is given in units of L. Furthermore, the dimensionless parameter ℓ = LH i specifies the characteristic scale of the spherical overdensity in units of the initial Hubble radius. In the present paper, we consider an initial overdensity, that is, A > 0.
To better understand the dynamical evolution and also in connection with the discussion of Section 5, it is useful to calculate the comoving sound horizon scale, which we estimate as: where the background solutionā(t) ∝ t 2/3(1+w) was used. Then we adopt ν 0 (t) = ν(t, r 0 ) = 0, where we set r 0 = ∞. Thanks to this choice, when w = 0 the metric reduces to the LTB metric in synchronous gauge. However, in the presence of the pressure the lapse function cannot be gauged away and is given by (2.11): We then fix the spatial gauge freedom left with a i (r) = 1. The initial condition for the gravitating mass is then from (2.19): , and, therefore, the last term in the previous equation is the average integrated density contrast.
Regarding E(t, r), from (2.26) we obtain: , (3.5) that is, the function E(r) gives the initial spatial curvature, as it is clear from (2.18). We will consider two different choices for the curvature function E(r) and therefore for E i (r).

Zero initial spatial curvature
The first case which we will consider is of zero initial spatial curvature, that is: An initial condition for E i (r) is equivalent to an initial condition forȧ i (r), that is, the initial velocity of the various shells. From (2.22) one finds: that is an initial velocity which is proportional to the initial integrated overdensity. This initial condition appears to be dominated by decaying modes as overdensities are expanding away faster than the background. In the previous equation e −ν i (r)ȧ i (r) = da/dτ is the derivative of a i (r) with respect to the proper time of the comoving observer/fluid as d/dτ = v µ ∂ µ = e −ν ∂/∂t.

Uniform initial expansion
The second case we will consider is of initial uniform expansion, that is: which implies: and E i (r) is given by (3.5). Clearly, it is E(r = ∞) = 0. This initial condition features a mixture of growing and decaying modes. The dynamical evolution relative to these two initial conditions was studied by Ref. [31] in the case of the pressureless LTB model. The initial spatial curvature is set to zero as explained in Section 3.1.1. In the top panels, the thick solid curve is the contrast δ(r), the dashed curve is the velocity V defined in (4.14), and the thin solid curve is the density contrast in the linear theory (see the Appendix A). See Figure 4 (left panels) for the time evolution of the density contrast. The second panels from the top show the corresponding profiles of the lapse function ν(r) and E(r). The details of the lower two panels are explained in Section 4. The third panels from the top show ξ(r) and H i ξ 0 (r), while the bottom panels show φ(r), ψ(r), and φ(r) + ψ(r). In these two panels, the curves are the results of the full nonlinear transformation defined by (4.8-4.11). See Section 3.2 for details.

Numerical Analysis
The validity of our numerical solutions is tested using the following two methods. First, we use the remaining extra equation (2.16), which can be rewritten as The terms in the first bracket of this equation do not include differentiation with respect to t, while the terms in the second bracket only include differentiation with respect to t. Therefore, we use (3.10) to check our numerical calculations. Figure 1 exemplifies a snapshot of our numerical results, demonstrating the balance between the first term (solid curve) and the second term (dashed curve). Second, in Appendix A we develop the perturbation theory  for the Lemaître model, whose prediction for the density perturbation is plotted (thin solid curves) against the fully nonlinear solutions in the top panels of Figures 2, 3, 5, 6 and in the panels of Figures 4 and 7. When the amplitude of the density contrast is not large (A < ∼ 0.1), the perturbation theory well reproduces the full numerical results based on the equations presented in Section 2. Figure 2 shows our numerical results for the case of a relativistic fluid (w = 1/3) with initial parameters A = 0.3 and ℓ = LH i = 1, that is, a characteristic scale of one initial Hubble radius. The initial spatial curvature is set to zero as explained in Section 3.1.1. The left and right panels of Figure 2 show snapshots at H i t = 5 (ā = 3.3) and at H i t = 20 (ā = 6.4), respectively. In the top panels, the thick solid curve is the density contrast δ(t, r) and the dashed curve is V (r) defined by equation (4.14), which is understood as the velocity field in the Newtonian description. One can see that the central overdensity propagates away, then a wave-like feature appears at a later time (see the left panels of Figure 4 for a more detailed time evolution). This is a consequence of the adopted initial conditions which are dominated by decaying modes and cause the density contrast to decrease. The linearized theory (thin solid curve, see Appendix A) is in agreement with the exact solution as the density contrast quickly becomes linear. In the second panels from the top, the solid curve is the lapse function ν(r) and the dashed curve is E(r). The details of the lower panels are explained in the next Section. In Figure 3 we repeat the calculation for the case of uniform initial expansion (see Section 3.1.2). All the other model parameters are unchanged. From the top panels of Figure 3, one sees that in this case the density contrast does not decrease as compared with Figure 2 (see also the right panels of Figure 4). This should be due to the fact that the initial conditions are dominated by growing modes and the characteristic scale of the inhomogeneity is initially larger than the sound horizon scale. Indeed, from (3.2) one obtains r s = L at H i t = 1, r s = 1.9L at H i t = 5 and r s = 3.7L at H i t = 20. The increase of r s at later times explains why the density contrast did not continue to grow but after reaching a maximum value of 0.7 at H i t = 3 it then decreases to about 0.3 at H i t = 20 (see also the right panels of Figure 4). The nonlinear central density contrast causes the deviation between the full numerical result and the linear perturbation theory as shown by the right panels of Figure 4. Figure 5 shows the evolution relative to the initial parameters A = 0.3 and ℓ = H i L = 0.01 at the times H i t = 0.005 (left panels) and H i t = 0.02 (right panels), respectively. Figure 6 shows the evolution relative to the same parameters of Figure 5 except for a larger initial contrast of A = 1. The sound horizon scale is in both cases r s = 58L and r s = 59L at H i t = 0.005 and H i t = 0.02, respectively, much larger than the inhomogeneity scale. We repeated the computation relative to Figures 5-6 for the case of the uniform initial expansion described in Section 3.1.2 and found that the evolution is basically unchanged. Therefore, we conclude that the dynamical evolution is insensitive to the initial conditions for the curvature function E(r) when the characteristic inhomogeneity scale is much smaller than the sound horizon scale. Equivalently, pressure gradients quickly overcome the initial velocity profile. By comparing the thick solid curve and the thin solid curve in the top panels of Figure 6, we  Figure 7. Time evolution of the density contrast relative to Figure 5 (left panels) and Figure 6 (right panels). In the left panels the initial amplitude is A = 0.3, while in the right panels the initial amplitude is A = 1.
see that the linear perturbation description is good at an earlier stage of the time evolution even for the case of a nonlinear initial overdensity contrast A = 1 (see also Figure 7). However, the discrepancy between linear and nonlinear dynamics becomes large at a later time, when wave-like features appear. The amplitude of the wave-like features disappear faster for the cases of A > ∼ 1 than for the cases of A < ∼ 1, which is understood as a nonlinear effect of strong gravity (compare the relative amplitudes in Figure 7). Summarizing, the linearized theory agrees better with the exact solution when the contrast is smaller -as is clear from

Newtonian potentials from coordinate transformation
We now consider whether the Lemaître model can be brought to the form of the perturbed FLRW universe in the Newtonian gauge. This problem was considered for the case of a dust fluid (w = 0) in [2][3][4]. To our knowledge, the case w = 0 has not been investigated so far.
From Figures 2 and 3 one can see that for a relativistic fluid (w = 1/3) with initial density contrast A = 0.3 and characteristic scale ℓ = 1, the relation ψ+φ = O(φ 2 ) is satisfied. In Figure 2 the potentials are well in the linear regime. This is explained by the fact that the adopted initial conditions are dominated by decaying modes and the density contrast quickly becomes linear. In Figure 3, instead, the initial conditions are dominated by growing modes and the density contrast grows. This causes larger values of the potentials which, at the initial times, are only marginally within the linear regime. This is essentially due to the fact that at initial time the inhomogeneity is horizon scale. Even though the amplitude of the potentials becomes mildly nonlinear reaching the numerical value of 0.1 ∼ 0.2, the Newtonian description works in the sense that ψ + φ = O(φ 2 ) is satisfied. Figures 5 and 6 show that in the case of sub-sound horizon perturbations the Newtonian potentials satisfy |φ|, |ψ| ≪ 1, which indicates that the perturbed Newtonian description is valid even when the density contrast is large (see Figure 6)). However, the relation ψ + φ = O(φ 2 ) is not satisfied, that is, there is an effective anisotropic stress, a well known prediction of second-order perturbation theory [see, e.g., [20][21][22][23][24][25]. Interestingly, the same conclusion is reached for a linear initial contrast (A = 0.3, see Figure 5) as well as for a nonlinear initial contrast (A = 1, see Figure 6), and, as said before, is independent of the adopted initial conditions for the curvature function E(r). In other words, the Lemaître model allows us to study in an exact nonlinear fashion the onset of anisotropic stress in fluids with non-negligible pressure.
We then Fourier expand the first-order quantities: where j ℓ (z) is the ℓ-th order spherical Bessel function of the first kind, so that equations (5.13-5.17) become: where , η denotes differentiation with respect to the conformal time defined by dη = dt/ā, H =ā ,η /ā and we used the background Friedmann equation H 2 = 8πGā 2ρ /3. In particular, asρ ∝ a −3(1+w) , the Friedmann equation yieldsā(t) ∝ t 2/3(1+w) ∝ η 2/(1+3w) . Note that it is Ψ + Φ = 0 at the first order of perturbative expansion. By solving equations (5.22) and (5.23) for ∆ and V , respectively, equation (5.26) gives: where we used Φ + Ψ = 0. The solution is given by: where we defined , (5.29) and J ν (z) and N ν (z) are the Bessel functions of the first and second kind, respectively, and α k and β k are the coefficients specified by the initial conditions. In the limit w = c 2 s , it is 30) which is consistent with the results of Ref. [32].

Second-order equations
Next, we consider the equations for the second-order perturbations. From the second-order Einstein's equation, we have δG r ), which gives: Since we aim at understanding the cause behind the violation of φ + ψ ≃ 0 (see the results of Section 4), we now focus on the relation between Φ 2 and Ψ 2 . Using the first-order equation which can be integrated in order to make explicit the source of Φ + Ψ: 5 Next we compare the two source terms inside the square brackets of (5.33). Using the firstorder solution, the first term is estimated as Regarding the second term, from (5.23) one has ( Φ ,η = c s k Φ): which allows us to conclude that: Therefore, the second term inside the square brackets of (5.33) is larger than the first term if: which simplifies to c s k/H ≫ O(1), as long as w is not close to −1. Thus the relation Φ 2 + Ψ 2 ∼ O(Φ 2 1 ) may become invalid when the inhomogeneity scale is smaller than the sound horizon and the contribution of the velocity field is significant.
We conclude the analysis by carefully estimating the value of Φ 2 + Ψ 2 . Introducing L as the typical scale of the inhomogeneity, from (5.32) one finds that the first term on the righthand side of (5.33) is of order Φ 2 1 . For scales smaller than the sound horizon (L ≪ c s H −1 ) from (5.35) one has V 1 ∼ c s Φ 1 /(1 + w)(LH) 2 . Consequently, in this limit, the second term on the right-hand side of (5.33) is of the order c 2 s Φ 1 2 /(1 + w)(LH) 2 . If we further assume that L ≪ H −1 , that is a sub-horizon inhomogeneity, from (5.22) one obtains for the density contrast ∆ 1 ∼ Φ 1 /(LH) 2 . Then, the second term on the right-hand side of (5.33) becomes c 2 s ∆ 1 Φ 1 /(1 + w). We can finally conclude that: for an inhomogeneity whose characteristic scale is smaller than both the Hubble horizon and the sound horizon, so long as w is not close to −1. The result of (5.38) is based on a perturbative expansion. Therefore, the second-order quantities must be smaller than the firstorder ones, that is, |Φ 2 |, |Ψ 2 | < |Φ 1 |, |Ψ 1 |. From (5.38) one then concludes that c 2 s Φ 1 ∆ 1 < Φ 1 , so that for consistency of the perturbative expansion we require ∆ 1 < 1.
On the other hand, if the characteristic scale of the inhomogeneity is larger than the sound horizon (L ≫ c s H −1 ), the second term on the right-hand side of (5.33) becomes Φ 1 2 /(1 + w) (see (5.36)), which is of the same order as the first term as long as w is not close to −1. The analytical results of this Section confirm what we have found numerically in the analysis relative to Figures 2-3 and 5-6.

Conclusions
We studied the Lemaître model with non-negligible pressure in a cosmological setup using a numerical method. We investigated the validity of the Newtonian description of the Lemaître model by transforming its metric into the perturbed FLRW metric. For the case of a relativistic fluid (w = 1/3), the Newtonian description is valid as the amplitudes of the potentials remain small |φ|, |ψ| ≪ 1, even when the density contrast is of the order O(1), as long as the inhomogeneity is sub horizon. However, the relation φ + ψ = O(φ 2 ), which holds in the linear cosmological perturbation theory, ceases to be valid. In the dust-dominated case, it is known that the deviation is expressed by the second-order potential, φ + ψ = O(φ 2 ). However, this is not the case for an inhomogeneous fluid with non-negligible pressure. This suggests that density inhomogeneities in a fluid with non-negligible pressure may effectively give rise to anisotropic stress. Therefore, the Lemaître model gives us the chance to study in an exact nonlinear fashion the generation of anisotropic stress in fluids with non-negligible pressure, thus extending previous work based on second-order perturbation theory [see, e.g., [20][21][22][23][24][25]. From the analysis of the second-order perturbations, we found that φ + ψ = O(φ 2 ) is not valid when the characteristic scale of the inhomogeneity is smaller than the sound horizon, in which case it is φ where δ is the density contrast, as long as w is not close to −1. As this estimation is based on a perturbative approach, δ < 1 is required for consistency. We also conclude that ψ + φ = O(φ 2 ) holds in general when the characteristic scale of the inhomogeneities is larger than the sound horizon scale, unless the amplitude of the inhomogeneities is nonlinear and w is close to −1.
As discussed in Section 3, at the initial time t = t i we choose a(t i , r) = 1, that is, ζ 1 (t i , r) = 0. Therefore, at the initial time t i from (A.8) we have: which, using (A.9) and (A.11), becomes: The Fourier coefficientsÃ k andB k are then given by: The following formula is useful to evaluate C k as we adopt δ i (r) = Ae −r 2 /L 2 (see Section 3):