Large-Scale Tides in General Relativity

Density perturbations in cosmology, i.e. spherically symmetric adiabatic perturbations of a Friedmann-Lema\^itre-Robertson-Walker (FLRW) spacetime, are locally exactly equivalent to a different FLRW solution, as long as their wavelength is much larger than the sound horizon of all fluid components. This fact is known as the"separate universe"paradigm. However, no such relation is known for anisotropic adiabatic perturbations, which correspond to an FLRW spacetime with large-scale tidal fields. Here, we provide a closed, fully relativistic set of evolutionary equations for the nonlinear evolution of such modes, based on the conformal Fermi (CFC) frame. We show explicitly that the tidal effects are encoded by the Weyl tensor, and are hence entirely different from an anisotropic Bianchi I spacetime, where the anisotropy is sourced by the Ricci tensor. In order to close the system, certain higher derivative terms have to be dropped. We show that this approximation is equivalent to the local tidal approximation of Hui and Bertschinger (1995). We also show that this very simple set of equations matches the exact evolution of the density field at second order, but fails at third and higher order. This provides a useful, easy-to-use framework for computing the fully relativistic growth of structure at second order.


Introduction
Cosmological observations ranging from the cosmic microwave background (CMB) to the clustering of galaxies strongly suggest that our universe is well described by a Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime on large scales, with metric ds 2 = a 2 (τ )[−dτ 2 + γ ij dx i dx j ] . (1.1) Here a(τ ) is the scale factor, τ is conformal time (with dτ = dt/a), and γ ij is the metric of a maximally symmetric three-dimensional space with either positive, negative, or vanishing curvature. However, observations also prove that there are perturbations to the spacetime Eq. (1.1) whose treatment is essential in order to correctly interpret cosmological observations. Interestingly, a special class of perturbations of Eq. (1.1), namely spherically symmetric adiabatic perturbations, are locally exactly equivalent to a different FLRW solution a → a F , where a F obeys the Friedmann equations, as long as their wavelength is much larger than the sound horizon of all fluid components. This is commonly known as the "separate universe" picture [2][3][4][5][6][7][8][9][10][11][12], and holds at fully nonlinear order [13]. That is, an observer within such a "perturbed" spacetime cannot distinguish the spacetime from exact FLRW by any local measurements, where local means on spatial scales much smaller than the wavelength of the perturbation. Note that this holds for all time, i.e. the observer could keep making spatially local measurements for several Hubble times and would continue to find agreement with FLRW.
However, no such relation to an exact solution of Einstein's equations is known for anisotropic adiabatic perturbations around FLRW. While the spherically symmetric perturbations mentioned above can be equivalently seen as density or curvature perturbations, the anisotropic case is equivalent to large-scale tidal fields. It is sometimes argued that the exact solution corresponding to this case is a Bianchi I spacetime, where δ ij + H ij (τ ) is a trace-free symmetric positive-definite matrix. Indeed, the motion of a non-relativistic test particle in such a spacetime is equivalent to that in an FLRW metric perturbed, in conformal-Newtonian gauge, by a tidal potential perturbation that can locally be written as (e.g., Sec. 7 of [14]) where dots denote derivatives with respect to time t. However, as we will show below, this is not the correct physical picture. In a dust-filled universe (and indeed any universe dominated by an ideal fluid), the Bianchi I solution Eq. (1.2) leads to a rapidly decaying anisotropy H ij ∝ aH (see Sec. 2). Any other behavior would require significant anisotropic stress. The actual tidal fields in our universe on the other hand grow in conjunction with the perturbations in the matter density. It follows that the tidal fields in our universe are not sourced by any component of the Ricci tensor, but are instead encoded by the Weyl tensor. We show this by providing a closed, fully relativistic set of evolutionary equations involving the matter density, velocity shear, local scale factor a F , and local tidal field E F ij . E F ij is a specific component of the Weyl tensor. It can be interpreted as the electric Weyl tensor evaluated on the geodesic or simply the "local tidal field" (motivated by the expression for the 00−component of the metric perturbation in our framework, shown later). Specifically, these quantities are defined in the conformal Fermi (CFC) frame [14,15], which ensures that each one of them is a local observable from the point of view of an observer comoving with the matter fluid. Beyond elucidating the physical interpretation of tidal fields in the relativistic context, the result is also useful for estimating post-Newtonian corrections appearing in nonlinear cosmological perturbation theory. The separate universe picture proves that all local gravitational effects of isotropic adiabatic metric perturbations (at all orders in perturbation theory) are captured by the spatial curvature in the comoving frame [13]. Our results prove that the corresponding quantity for anisotropic perturbations is the electric part of the Weyl tensor. Moreover, written in terms of locally observable quantities, the evolutionary equations only contain terms that are familiar from the subhorizon, Newtonian calculation (although, of course, they do contain post-Newtonian terms once expressed in a specific gauge such as Poisson gauge). However, in order to a obtain closed set of equations, certain terms need to be dropped. This local tidal approximation (LTA) [1] only recovers the correct physical evolution of the tidal field at linear order, and the density field at second order. Hence, anisotropic (tidal) adiabatic perturbations are significantly more complex than isotropic (density) perturbations.
The outline of this paper is as follows. In Section 3, we briefly recap the idea of Conformal Fermi Coordinates (CFC). In Section 4, we derive a closed system of equations for the nonlinear evolution of the density and tidal fields in the CFC frame from the corresponding covariant system of equations. In Section 5, we present the perturbative solutions for the quantities in the system up to second order and compare them to the known solutions in standard perturbation theory. We also relate the system to the collapse of a homogeneous irrotational ellipsoid and point out that, when restricting to the leading order perturbation in the CFC metric, our system matches the Local Tidal Approximation introduced by Hui and Bertschinger [1] for ellipsoidal collapse. In Section 6, we consider the rest-frame matter three-point function as an application of our results. In Section 7, we summarize our findings and discuss their applications.
2 Why a Bianchi I spacetime does not describe large-scale perturbations in our Universe We begin from the Bianchi I metric Eq. (1.2), and rotate the spatial coordinates to the frame where H ij is diagonal. Then, we can write for all t. In other words, we write Eq. (1.2) as Jacobs [16] derived the solution of the Einstein equations for the ansatz Eq. (2.3), assuming a perfect-fluid stress tensor (he also considered the case of a uniform magnetic field, which we do not discuss here), and obtained [Eq. (10) there] where H =ȧ/a, while Π, Σ are the anisotropy parameters. The isotropic scale factor a(t) satisfies the standard Friedmann equation. Thus, apart from an unobservable constant mode, which can be removed by a trivial redefinition of spatial coordinates, the anisotropy always decays asΠ,Σ ∝ a −3 . The tidal field experienced by a local comoving observer in this spacetime [Eq. (1. 3)] then scales as and analogously for Σ. That is, in a universe whose stress energy content is given by a perfect fluid (more specifically, in the absence of significant anisotropic stress), any initial tidal field described by a Bianchi I spacetime decays rapidly, ∝ aH. In case of initial conditions from inflation, which are set when a given mode exits the horizon, this contribution disappears exponentially fast after horizon exit, as expected for a decaying mode. This is clearly very different from actual tidal fields, which do not decay outside the horizon, and whose importance in the growth of structure remains relevant up to late times.

Recap of conformal Fermi coordinates
We begin by recapitulating the gist of the Conformal Fermi Coordinates (CFC), t F , x i F , rigorously defined in [14]. Quantities defined with respect to CFC shall be denoted by a subscript F . Take an observer free-falling in some general spacetime. His worldline is then a timelike geodesic in said spacetime. For most applications of these coordinates, we take the spacetime to be perturbed Friedmann-Lemaître-Robertson-Walker (FLRW), although this is not a necessary assumption. We construct a coordinate system centered on the observer, such that he always sees an unperturbed FLRW spacetime on the geodesic, with corrections going as the spatial distance from him squared, O x i F 2 . Note that a power expansion in x i F requires that |x i F | be less than the typical scale of variation of h F µν , which we shall call Λ −1 . Since in a realistic cosmological setting, there exist metric perturbations on very small scales, we need to coarse-grain the metric (and the stress energy tensor) on a spatial scale Λ −1 . Then, the only contributing modes in the resulting metric perturbations have wavenumbers k Λ. CFC is then valid over a finite region [14]. Since the absolute scale of the coarse-graining is not of relevance to the results of this paper, we will explicitly indicate the scale Λ.
We can now proceed to construct the CFC frame for the coarse-grained metric. We start by taking the tangent U µ = dx µ /dt F to his worldine, i.e. the fluid 4-velocity, to be the time direction and the hypersurface composed of all vectors orthogonal to his worldline to be the constant-time hypersurface, with the observer being at the spatial origin x i F = 0. The orthonormal tetrad thus chosen at a point on his worldline is then parallel transported along the latter, such that these properties are preserved. Then, given a scalar a F (x) that is positive in a neighborhood of the geodesics, the spatial coordinate for a given t F is defined as follows: t F , x i F is the point at which we arrive when starting from P = (t F , 0) (on the observer's worldline), we move along the spatial geodesic ofg F µν for a proper distance of a F (P ) δ ij x i F x j F in the direction defined w.r.t. the spatial components of the observer's tetrad.
For any given spacetime scalar a F (x) 1 , the metric can thus be made to take the form Here,R F µναβ (0) is the Riemann tensor of the conformally-related metric, evaluated on the central geodesic.
So far, we have not specified the local scale factor a F . As shown in [13,14], the natural, physically motivated choice is to define the local Hubble rate along the geodesic through a F is then uniquely defined, up to an arbitrary overall multiplicative constant, by integrating the Hubble rate along the observer's geodesic. Apart from reducing to a F = a for an unperturbed FLRW spacetime, this choice ensures that H F as well as h F µν are strictly local observables from the point of view of the observer. 2

Nonlinear evolution of density and tidal fields
In this section, we derive a closed system of evolutionary equations for the density ρ, velocity divergence Θ, velocity shear σ µν , and tidal field E µν . Importantly, our relations will be derived at fully nonlinear order and fully relativistically, without assuming v 2 ≪ c 2 or subhorizon scales k ≫ H := aH as usually done in large-scale structure studies. Moreover, by giving expressions in the CFC frame, our results correspond directly to local observables from the point of view of a comoving observer.
We will assume that the cosmological fluid is pressureless (CDM) and perfect. It follows that there is no anisotropic stress and the stress-energy tensor takes the form T µν = ρU µ U ν , where ρ is the proper energy density in the fluid rest frame. However, our results also hold in the presence of pressure, as long as pressure perturbations can be neglected. This holds trivially for a cosmological constant. Moreover, it is valid as long as the long-wavelength perturbations considered are outside the sound horizon of all fluid components.
Let us introduce the projection tensor P µν (U ) := g µν + U µ U ν , which projects on the subspace orthogonal to the fluid 4-velocity (in CFC, P µν will become trivial). The floworthogonal part of the velocity gradient can be decomposed as [17][18][19] where Θ := P µν ∇ µ U ν is the expansion scalar. It describes the change in the volume of a sphere of test particles centered on the geodesic. The shear tensorσ µν is the traceless symmetric part of the velocity gradient tensor. It describes the rate of distortion of a sphere of test particles into an ellipsoid. Here, we have neglected the vorticity ω µν , which is the antisymmetric part of the (flow-orthogonal) velocity gradient. Since, for a single barotropic fluid, vorticity is not generated, and any initial vorticity decays, we set ω µν to zero throughout. Including this decaying mode is a trivial extension.

Covariant equations
Our goal is to generalize the separate universe picture, by finding a closed system of equations for the evolution of a homogeneous ellipsoid. The Friedmann equations governing the background evolution of an FLRW spacetime are a special case of the Raychaudhuri equation, While intrinsically a purely geometric relation derived from the Ricci identity, we have used the Einstein equations to replace the Ricci scalar with the trace of the stress-energy tensor on the r.h.s.. The Ricci identity holds for any Levi-Civita connection Γ α µν and is given by Here, [∇ µ , ∇ ν ] is the commutator of the covariant derivatives and R α βµν is the Riemann tensor of the metric g µν . It describes the difference between parallel-transporting U α in the direction ∇ µ then ∇ ν and vice versa. Contracting Eq. (4.3) with P ν α U µ then yields Eq. (4.2). We complement this with energy-momentum conservation, projected with the 4-velocity which for our CDM fluid reduces to the evolutionary equation for the rest-frame energy density, ρ.
In addition, we need an equation governing the evolution of the velocity shear, which can again be derived from the Ricci idenity Eq. (4.3), namely by contracting with U µ P σν P ρα and then P µσ . This yields and C ανµβ is the Weyl tensor [20], i.e. the traceless part of the Riemann tensor which describes the contributions of nonlocal sources to spacetime curvature. As we will see, its electric part E µν can be understood as the invariant definition of the local tidal field. Further, is the traceless part of the spatially-projected Ricci tensor. Eq. (4.5) shows that there are two sources of velocity shear: the electric part of the Weyl tensor, andR µν , which by virtue of the Einstein equations is proportional to the trace-free part of the velocity-orthogonal stress tensorT µν . In the absence of anisotropic stress in the fluid rest frame, which is the case for baryons and cold dark matter,R µν vanishes. Clearly, in the actual universe velocity shear is sourced by the electric part of the Weyl tensor, which is the relativistic generalization of the Newtonian tidal tensor (∂ i ∂ j − δ ij /3∇ 2 )Ψ. Now, in order to obtain a closed system of evolutionary equations, we need an equation governing E µν . This can be obtained by combining the Bianchi identity and the Einstein equation, This equation is the only one in the set that involves spatial-more precisely, fluid-orthogonalderivatives acting on the Riemann tensor (via the Weyl tensor) and the density ρ. Note that fluid-orthogonal derivatives are simply spatial derivatives in the frame comoving with the fluid. We are interested in describing the evolution of a long-wavelength perturbation. In order to obtain a closed system of evolutionary equations, we thus discard fluid-orthogonal derivatives acting on the density and the Riemann tensor and neglect them in the following (of course, throughout we retain terms involving any derivatives acting on the velocity). This is the key approximation made in our derivation, and we will discuss its implications in detail below.

Closed system in CFC frame
At this point, it is not obvious that the four relations Eqs. (4.2), (4.4), (4.5), and (4.8) can be rewritten to form a closed local system. However, this becomes clear once expressing all relations in the CFC frame. We denote˙= d dt F and H F :=ȧ F a F . We shall also adopt the convention where 3-indices are raised and lowered by the Kronecker delta, δ ij , whereas 4-indices are raised and lowered by P µν , which is the effective spatial metric for the observer in covariant form.
Further, it is sufficient to evaluate all quantities on the geodesic x i F = 0. This is because we are free to choose a geodesic through any given point. Thus, all terms that scale, after taking all spatial derivatives, as x i F n , n ≥ 1, vanish. This applies in particular to the Then, without making any further approximations, the nonlinear tensor equations simplify greatly. On the geodesic, we have U µ = a −1 F (1, 0, 0, 0), and the projection tensor simply becomes Firstly, Eq. (4.2) becomesḢ and ∂ F i ≡ ∂/∂x i F . Note that we have pulled out a factor of a −1 F in the definition of the velocity shear. Further, recall that in CFC, we define Θ := 3H F on the geodesic, such that the spatial velocity divergence, ∇ k,F v k F , is absorbed into the definition of a F . Thus, σ F ij is trace-free.
Secondly, Eq. (4.4) becomesρ which is unsurprisingly the familiar continuity equation. Thirdly, Eq. (4.5) becomeṡ Here, we have introduced E F ij as the relevant component of the Weyl tensor in CFC frame, which is directly related to the perturbation h F 00 of the 00-component of the CFC metric (see below). Further, we have used the fact that This vanishes on the geodesic when we evaluate R F µν using the trace-reversed Einstein equation because the CFC frame is defined to be the fluid rest frame, and we assume no anisotropic stress as discussed above.
Finally, in order to evaluate Eq. (4.8), note that the Weyl tensor is trace-free over any two indices and has the same symmetry properties as the Riemann tensor. Moreover, is a fluid-orthogonal derivative of the Riemann tensor and thus neglected as explained after Eq. (4.8). Moreover, in the CFC construction this term is naturally higher order, as for the leading expression given in Eq.
vanishes when evaluated on the observer's worldline. It follows then that the 0ij0−component of the LHS of Eq. (4.8) becomes , as the spatial derivatives acting on U µ yield nonnegligible terms. Since our goal is to derive a closed system in terms of local gravitational observables in the framework of the CFC, we choose Eq. (4.16) as the local approximation to Eq. (4.8). In order to emphasize the subtle distinction between E µν (U ) and C F i0j0 in the context of our local approximation, we have introduced the new symbol E F ij in Eq. (4.14). The physical interpretation of E F ij becomes clear when deriving its relation to the perturbation h F 00 of the 00-component of the metric in CFC Eq. (3.2). Using the transformation law of the Riemann tensor under a conformal rescaling of the metric, one obtains where we have used Eq. (4.10) in the second line. Using the definition of the Weyl tensor as trace-subtracted version of the Riemann tensorR F , we then obtain [see also Eq. (3.31) in [13]] where we have used Eq. (4.14). Note that since the Weyl tensor is invariant under conformal rescaling, we haveC F = C F . We see that E F lm = C F 0l0m /a 2 F is exactly the trace-free part of −∂ l ∂ m h F 00 , which is the local tidal field acting on non-relativistic bodies in the CFC frame. Any other Ricci-contribution to the local tidal field would have to be due to anisotropic stress.
The four equations Eqs. (4.10)-(4.18) now clearly form a closed, local system of ordinary differential equations governing the evolution of the four unknowns ρ F , H F , σ F , and E F along the fluid geodesic:Ḣ Note that E F has dimensions 1/length 2 , while σ F has dimensions 1/length. Given initial conditions for ρ F , H F (or equivalently, curvature K F ), σ F ij , and E F ij , Eq. (4.21) can be integrated numerically without any further approximations. Since we have used the Einstein equations through Eq. (4.2) and Eq. (4.8) (with all components of the latter derived in App. A), all constraints are self-consistently fulfilled at leading order in derivatives. In the next section, we will consider the general perturbative solution around a flat matterdominated (Einstein-de Sitter) universe. The fully relativistic system Eq. (4.21) clarifies the physical meaning of the locally observable impact of tidal fields. That is, tides are a manifestation of the Weyl tensor, which encodes the nonlocal part of gravitational interactions, and thus of the large-scale inhomogeneities in the matter distribution. This is to be contrasted with homogeneous and anisotropic Bianchi I spacetimes, where the anisotropy is sourced by the trace-free part of the Ricci tensorR ij . As we argued above, this in fact vanishes everywhere for a pressureless (and indeed any ideal) fluid. Thus, the construction leading to Eq. (4.21), rather than a Bianchi I solution, are the proper relativistic model of long-wavelength density and tidal perturbations.
Note that the system simplifies further when setting σ F = 0 = E F initially, corresponding to a spherically symmetric system. The symmetry is preserved so that σ F and E F remain zero. We then obtainḢ These are just the second Friedmann equation and continuity equation of the FLRW spacetime. One can then show (again, assuming scales much larger than the sound horizon of the fluid) that the first Friedmann equation is satisfied as well [13]. This proceeds in the usual way by multiplying Eq. (4.22) withȧ F , integrating once, and using the continuity equation to yield where K F , the local spatial curvature, is an integration constant, hence is conserved. This is the well-known separate universe picture: long-wavelength isotropic adiabatic perturbations are indistinguishable from an FLRW spacetime (with different cosmological parameters) by local observations. Eq. (4.22) is also equivalent to the spherical collapse equation [21]. Note that Eq. (4.22) is exact given the said assumptions, while in the anisotropic case Eq. (4.21) is an approximation whose accuracy we will discuss further below. First, however, we will proceed to solve the general, anisotropic case.

Perturbative solution in Einstein-de Sitter
We now consider a perturbed Einstein-de Sitter universe, where in the background the Hubble rate satisfiesH 2 = 8πGρ/3 ∝ā −3 . We expand all quantities into orders of perturbation, denoted by [n] , n = 1, 2, · · · . Quantities with an overbar are evaluated in the background Einstein-de Sitter universe. Correspondingly, we write the density as Throughout, we assume that initial conditions are set at sufficiently early times (see Sec. 5.1 below), so that we can restrict to the fastest growing modes throughout. This is merely for calculational convenience; it is straightforward to keep the subleading modes when solving Eq. (4.21).

Initial conditions
We briefly describe how the growing-mode initial conditions for Eq. (4.21) can be determined. For this, we assume they are set at sufficiently early times, so that linear perturbation theory is accurate. We can then relate the CFC-frame quantities to those in a given global coordinate system at linear order in perturbations. Specifically, we consider two frequently used gauge choices: conformal-Newtonian gauge, defined through and comoving gauge, which is frequently used for calculations during inflation: Here, spatial slices are chosen such that T 0 i = 0, hence the designation "comoving." In both cases, we have only included scalar perturbations. The reason is that vector perturbations are pure decaying modes in both cases, so that they are irrelevant for the fastest growing modes. Tensor modes on the other hand are propagating modes whose nonlinear evolution we do not expect to be described correctly by the local approximation employed in Eq. (4.21). The linear evolution of a tensor-mode-induced E F ij and its effect on small-scale perturbations was derived in [14,22]. We stress again that this simplification is merely for convenience and not required within the CFC formalism; one can straightforwardly include vector and tensor modes and their associated decaying modes.
Let us first consider the simpler, isotropic case, where the initial conditions can be either specified through δ F or K F . As shown in [13], one has Note that K F = const at all times (not only during matter domination), as long as it is outside the sound horizon of all fluid components. During matter domination and at linear order, one further has where Ω m0 is the matter density parameter today.
Next, consider the anisotropic case. Restricting to the growing mode, it is sufficient to provide initial conditions for E F ij . In conformal-Newtonian gauge, we have at linear order in perturbations where we have used the CFC metric constructed at linear order by [13], and Eq. (4.20). 3 Moreover, using that the two spacetime potentials in conformal-Newtonian gauge are equal in the absence of anisotropic stress, and setting vorticity to zero, it is easy to show that Eq. (4.18) yields the correct evolution of E F ij at linear order [cf. Eq. (4.50) in [20]]. Note that unlike the curvature K F , E F ij is in general not conserved during cosmic evolution even outside the sound horizon of all fluid components. This is a qualitative difference to the isotropic case, where the effect of the long-wavelength perturbation is described by K F which is constant at all orders on large scales. However, in the particular case of matter domination, E F ij is conserved at linear order (see below). In this case, initial conditions can be simply specified by using the well-known relation for growing-mode adiabatic perturbations during matter domination, R = (5/3)Φ, so that This relation can be used immediately to initialize a calculation for the nonlinear evolution of the tidal field E F ij for modes that enter the horizon during matter domination. In general, one should follow the linear evolution of E F ij via Eq. (5.6) until matter domination (for example, using a Boltzmann code), at which point E F ij approaches a constant and can be matched to the nonlinear calculation in matter domination which we describe next.

Perturbative solution up to second order
We begin with the linear evolution in CFC, which matches that of standard (subhorizon) perturbation theory [23], and was derived in [13]. In terms of the matter density perturbation and electric Weyl tensor, we obtain Here, we have introduced the linearly-extrapolated initial overdensity δ F (x F , t 0 ) and the local tidal field, E [1] F ij (x F ), for later convenience. Note that δ [1] F ∝ā, while E [1] F ij = const at linear order. One could thus think of E F ij as the analog for anisotropic perturbations of the spatial curvature K. Unlike the latter however, we will see that E F ij is not conserved at nonlinear order. Moreover, even at linear order the conservation of E F ij only holds in a flat matter-dominated universe.
The local scale factor and velocity shear are at linear order given by Note thatā(t F ) = (t F /t 0 ) 2/3 -the familiar scale factor in the unperturbed FLRW spacetime. Continuing to solve the equations Eq. (4.21) perturbatively, and keeping only the fastest growing mode (i.e. the term with the highest power of t F ), gives, at second order in perturbations, F lm E lm F [1] + · · · (5.14) where the left-hand side is evaluated at t F and x F ?, while on the right-hand side the time dependence is completely encoded in the factors ofā =ā(t F ). Instead of writing the solutions in terms of the initial conditions at some early time, we have followed the customary choice in cosmological perturbation theory of phrasing the solution in terms of the linearlyextrapolated initial overdensity δ [1] L (x F ) and initial electric Weyl contribution E [1] F introduced above. The results Eqs. (5.12)-(5.15) can be formally continued straightforwardly to higher order. Further, in order to describe the solution in a ΛCDM background, one can perform the standard very accurate approximation of replacingā n with the linear growth factor [D(t F )] n .

Connection to standard perturbation theory
It is instructive to compare our result Eq. (5.13) to the second order density perturbation in standard (subhorizon) perturbation theory [9,23], which yields δ(x, t) = δ [1] where q[x, t] denotes the Lagrangian position corresponding to the given Eulerian coordinate, and where we have used the Einstein-de Sitter background. The second term in Eq. (5.16) describes the second order growth of density perturbations in the absence of tidal fields; its coefficient is exactly what is obtained from the second order expansion of the spherical collapse solution. The third term encodes the tidal effects on the density perturbations. In standard Eulerian perturbation theory, the first, linear term is expanded around the Eulerian position yielding a displacement term δ [1] (q[x, t], t) = δ [1] (x, t) − s [1] (x, t) · ∇δ [1] (x, t) , (5.18) where s [1] is the linear displacement from the Lagrangian to the Eulerian position. Since the CFC calculation corresponds to working in Lagrangian coordinates (as the origin of the coordinate system comoves with the fluid, evidenced by the fact that v F = 0 on the geodesic), this displacement does not appear in Eq. (5.13); that is, the CFC formalism automatically resums all the displacement terms appearing in Eulerian perturbation theory. Using that we immediately see that Eq. (5.13) matches the standard perturbative calculation Eq. (5.16). However, in the derivation of Eq. (5.13) we have not assumed that the scale of the perturbation is much smaller than the Hubble horizon. This shows that, when measured in the proper rest frame, the second order evolution of the matter density in the presence of tidal fields is exactly as given by the standard result which is seemingly only valid on subhorizon scales. This fact was already known for isotropic perturbations, in which case the evolution is determined by a local Friedmann equation [13]. Eq. (5.13) generalizes this result to the anisotropic case, i.e. the inclusion of tidal fields. However, the agreement between the evolution predicted by the closed system Eq. (4.21) and standard perturbation theory no longer holds at higher order. This can already be seen in the result for σ F ij at second order, Eq. (5.15). The SPT prediction for σ ij , which corresponds to the derivative with respect to Lagrangian coordinates of the fluid velocity v, can be derived by using that v is equal to the convective (or Lagrangian) time derivative w.r.t. τ of the Lagrangian displacement s. This yields (e.g. App. B of [24]) F lm E lm F [1] . (5.20) Clearly, this differs from Eq. (5.15). In particular, the SPT result is spatially nonlocal (the same holds when deriving the SPT result for the nonlinear evolution of E F ij ). The differences go back to the terms neglected when evaluating the evolution equation Eq. (4.8) for the electric Weyl tensor component. When neglecting these terms, we were able to obtain a closed system that is local in space around the fluid geodesic. However, the gravitational evolution of density and tidal fields within a region, when followed over a finite duration of time, is not completely described by the local tidal and density field. This is encoded in the nonlocal terms appearing at second order in σ ij in standard perturbation theory, which, apart from dropping post-Newtonian terms, does not make approximations in Eq. (4.8). Note that the nonlocal term appears in the density perturbation δ F only at third order (see also [25,26]); in fact the nonlocal contributions to δ [3] F are proportional to σ [1]lm σ [2] lm [24]. Finally, we see that the nonlocal terms disappear in the case of spherical symmetry, in which case Eq. (4.21) recovers the exact separate universe result which matches perturbation theory to all orders.
In the next section, we will connect our results to those from other, previously considered local approaches to nonlinear gravitational evolution of non-spherically symmetric systems.

Connection to the local tidal approximation and ellipsoidal collapse
Interestingly, the system Eq. (4.21) has been derived before, in the subhorizon limit, by Ref. [1] who referred to this as the Local Tidal Approximation (LTA). This approximation originated in a series of attempts at improving the local approximation to study large-scale structure and at clarifying the relation between general relativity and Newtonian dynamics.
The first was the Zeldovich approximation (ZA) [27]. Bertschinger and Jain [28] proposed the non-magnetic approximation (NMA). In the NMA, the covariant magnetic part of the Weyl tensor is set to zero, H ij = 0. However, Refs. [29,30] showed that H ij cannot be consistently neglected in the Newtonian limit. That is, H ij has a Newtonian counterpart after all. Later, Hui and Bertschinger proposed the LTA. Here they defined a new quantity M ij , composed of a combination of terms in the tidal evolution equation, and set it to 0. They found that whereas NMA is exact for spherical perturbations but not cylindrical ones; LTA is exact for both and, more generally, for any growing-mode perturbations whose gravitational equipotential surfaces have constant shape with time.
In the LTA approximation, the authors imposed the condition where conformal-Newtonian gauge is adopted and E ij ≡ (∂ i ∂ j − (δ ij /3)∇ 2 )Φ is the tidal field in the subhorizon approximation. Note that in the subhorizon limit, the velocity-orthogonal projection is equivalent to simply restricting to spatial coordinates. H ij here is the magnetic component of the Weyl tensor. The tensor M ij can also be written as The resultant E ij in LTA is in fact our local tidal field E F ij . Ref. [1] discussed the LTA in the context of the collapse of a homogeneous isolated ellipsoid [31,32]. In this model, one neglects the gravitational effect of the ellipsoid on the surrounding matter. Ref. [33] studied the validity of this approximation by performing an N-body simulation which allows for the backreaction of the ellipsoid on its environment. Specifically, the environment consisted of a small homogeneous negative density perturbation to compensate for the mass contained within the overdense ellipsoid. He found that the abovementioned approximation is quite accurate until the late stages of collapse. The advantage of neglecting gravitational backreaction on the environment is that a closed, semi-analytical solution can be obtained.
In particular, working in the conformal-Newtonian gauge and taking the subhorizon limit, the evolution of an irrotational, isolated homogeneous ellipsoid embedded in an FLRW background is governed by where R i are the proper axis lengths of the ellipsoid;ρ is the local homogeneous (or, mean) density around the ellipsoid, such that ρ :=ρ (1 + δ) ; t is the proper time; and α i are defined by The peculiar velocity field inside the ellipsoid is where x i are comoving spatial coordinates. These relations can be used to derive the local Hubble rate, velocity shear, and electric Weyl tensor component [1], which in CFC become Given this one-to-one mapping to quantities in the system Eq. (4.21), the isolated ellipsoid can be considered as another local approximation to the evolution of non-spherically symmetric perturbations, even though it is not truly local as the α i , and thus E F ij , obey an integral equation Eq. (5.24). As shown in [1] however, the LTA, and hence the closed system in CFC derived in Sec. 4.2, recovers the evolution of Eqs. (5.23)-(5.24) only up to first order. Thus, given the N-body results of [33], the LTA describes the actual evolution of an isolated ellipsoidal perturbation much less accurately than Eqs. (5.23)-(5.24). This difference can be attributed to the perfectly local approximation made in LTA, whereas the isolated ellipsoid does take into account the finite extent of the perturbation. However, in practice, perturbations in the universe cannot be considered isolated, at least until they have reached sufficient density to dominate the local gravitational potential. As discussed in Sec. 5.3, the correct evolution of tidal fields is fundamentally nonlocal starting at second order. This is qualitatively different from isotropic perturbations, where the LTA reduces to the "separate universe" which provides an exact solution for large-scale adiabatic perturbations.
6 Application: the rest-frame matter three-point function As derived in Sec. 5, the matter density perturbation in CFC frame up to second order in perturbations is given by δ F (q, t) = D(t)δ [1] L (q) + 17 21 [D(t)δ [1] L (q)] 2 + 2 7 [D(t)K [1] ij (q)] 2 , (6.1) where K [1] ij is defined through Eq. (5.17). Here, we have only assumed adiabatic Gaussian initial conditions. q is a spatial coordinate which labels the geodesic around which the CFC frame is constructed, while t denotes proper time along the geodesic. We will choose q to denote the spatial position on a constant-proper-time slice at t = 0, that is, at the end of inflation.
We now would like to relate this result to the observed, late-time statistics of the matter density field, without resorting to the sub-horizon approximation made in most large-scale structure studies; specifically, we are interested in the leading signature of nonlinear evolution, the three-point function. Unfortunately, the precise prediction for these statistics depends on how the matter density field is measured, for example through gravitational lensing (which strictly measures the deflection of photon geodesics), or tracers such as galaxies (whose relation to matter is complicated by bias, redshift-space distortions, and other effects; see [34] for a review). Since these issues go beyond the scope of this paper, we here assume the idealized case of observers on different geodesics communicating their local rest-frame matter density as well as proper time to a distant observer on their future light cone.
Let us define the displacement s as parametrizing the difference between the spatial position x, as observed by the distant observer, of the fluid geodesic relative to the initial position q as a function of proper time: x = q + s(q, t) . (6.2) Note that s depends on how the distant observer measures the spacetime position of the geodesic. In general, s describes the effect of large-scale gravitational perturbations on the geodesic, which are given by an integral over the gradient of the gravitational potential, as well as the details of the distant observer's measurement. For example, if they use apparent photon arrival directions and redshifts, then s includes Doppler shift, gravitational redshift, and deflection by structures along the line of sight. The explicit expression of s further depends on the coordinate (gauge) choice that is used to calculate the displacement, although the end result is independent of the gauge; see [35] for a brief overview. For this reason, we will not derive s explicitly here. However, independently of the gauge choice, s is first order in cosmological perturbations. Hence, we can expand the CFC-frame density field up to second order as Moreover, since only the cross-correlation of s with the density field δ F enters in the leading three-point function, it is sufficient to consider the longitudinal contribution to s, which we write as where S is a scalar function which, for a fixed Lagrangian position q, only depends on proper time t. Note that S has dimension of length 2 .
With a slight generalization of the results in App. A of [36], we can then write the three-point function of the CFC-frame matter density as − µ 12 ∂ξ Sδ (r 1 , t) ∂r 1 ∂ξ 0 (r 2 , t) ∂r 2 + (1 ↔ 2) + 2 cycl. perm. . (6.5) , while µ ij ≡ r i · r j /r i r j , and we have defined Finally, P L (k, t) is the power spectrum of the linearly extrapolated rest-frame matter density perturbations (this is equivalent to the density perturbations in synchronous-comoving gauge), while P Sδ (k, t) is the cross-power spectrum between S and δ, again at linear order. While Eq. (6.5) only represents an idealized case, which does not include the mapping from the local rest-frame matter density to a realistic observable, it clearly illustrates the simplicity of the result obtained, without restriction to subhorizon scales, when using CFC to describe nonlinear gravitational evolution. Given the assumption stated at the beginning of the section, Eq. (6.5) is fully accurate at second order, and not restricted to the squeezed limit.

Summary and discussion
The conformal Fermi coordinates (CFC) are a convenient construction designed to explicitly isolate the leading locally observable gravitational effects of long-wavelength perturbations in the cosmological context. We show that using this construction, there is a natural way to derive a closed system of ordinary differential equations [Eq. (4.21)] describing the evolution of a long-wavelength perturbation. Here, long-wavelength means that the scale of the perturbation is much larger than the sound horizon of all fluid components (note that the sound horizon for pressureless matter is the nonlinear scale, r s ∼ 1/k NL [37]). This system is exactly equivalent to the local tidal approximation (LTA). Moreover, the CFC frame provides a fully relativistic realization of this closed system; that is, the results are valid on horizon or super-horizon scales without any post-Newtonian corrections. As shown in [13,15], the CFC also allows for a direct connection to the initial conditions from inflation at nonlinear order.
This construction clarifies the anisotropic generalization of the separate universe result for spherically symmetric long-wavelength perturbations. That is, while the latter are entirely described locally by a curved FLRW spacetime, anisotropic long-wavelength perturbations contain a nonzero electric component E µν of the Weyl tensor as well. Physically, this is entirely different from a Bianchi I spacetime, where the anisotropy is encoded in the tracefree part of the spatial Ricci tensor. The latter is negligible for ideal fluids. We thus argue that the Bianchi I picture is not the proper local physical representation of long-wavelength perturbations in our Universe.
Solving the system of evolutionary equations up to and including second order, we have found that the solution for the density field is exactly equivalent to the result of standard, sub-horizon perturbation theory. Since our closed system in the CFC frame recovers the exact nonlinear evolution of isotropic perturbations, and recovers the correct linear evolution of anisotropic perturbations, our result for δ F represents the proper rest-frame matter density at second order in perturbations, including all relativistic corrections (assuming matter domination holds, and that primordial decaying modes can be neglected). This follows from the fact that anisotropic perturbations only contribute to the density at quadratic or higher order. Note that this statement includes any vector and tensor metric perturbations that are generated by anisotropic perturbations at second order [38][39][40]; these can only contribute to the rest-frame matter density at third order (we did not however include primordial vector and tensor modes, which would contribute at second order). A corollary is that, if we interpret the results of standard N-body simulations in this frame, any post-Newtonian corrections to the measured density field have to be third order in perturbation theory. This provides a strong constraint on their numerical relevance. On the other hand, our closed system fails for local tidal fields and velocity shear already at second order. Further, post-Newtonian corrections to non-ideal fluids, such as neutrinos, and on gravitational lensing are in general less suppressed than those for the matter density [41][42][43].
In this context, one might wonder whether this approach could be used to study the generation of gravitational waves from large-scale structure. However, for this one needs to define what gravitational wave means. A natural, physical definition is to derive the transverse-traceless component of metric perturbations in the far-field limit [44,45]. Unfortunately, the local nature of the CFC construction implies that we cannot derive the far-field limit of these metric perturbations in this approach.
As a first, simple application of these results, we have derived the leading three-point function of the CFC-frame matter density field in Sec. 6 [Eq. (6.5)]. While this expression is idealized in the sense that it assumes that local observers directly communicate their local density to a distant observer, it is fully valid on arbitrarily large scales, and not restricted to specific configurations such as the squeezed limit. A further interesting possible application of these results is the implementation of N-body simulations with a long-wavelength, non-spherically symmetric perturbation imposed; that is, the anisotropic generalization of the "separate universe simulations" presented in [5][6][7][10][11][12]. In principle, the effect of the electric Weyl tensor component E F ij can be simply included by adding an external force ∝ E F ij x j to any particle with position x. This however breaks the periodic boundary conditions inherent in conventional N-body simulations, and can thus cause numerical problems. We leave this issue for future work. Note also that, unlike the isotropic case, where the superimposed long mode can be made nonlinear [11], this is nontrivial in the anisotropic case, as the correct nonlinear evolution of E F ij is nonlocal [see discussion after Eq. (5.20)]. This is indeed the main caveat to all local approximations to the nonlinear evolution of (non-spherically symmetric) perturbations. At second order, the gravitational evolution of tidal field and velocity shear is spatially nonlocal. This comes as no surprise, given that we have shown that the tidal forces are due to the Weyl tensor, which captures those parts of the full Riemann tensor that are not locally related to the matter distribution. In case of the density field, this nonlocality only appears at third order. Thus, in all applications restricted to sufficiently low order, local approximations are exact; however, starting at third order in the matter (and galaxy) density field, spatial nonlocality is an unavoidable feature of nonlinear gravitational evolution. further like to thank Simon Nakach for helpful discussions. FS acknowledges support from the Marie Curie Career Integration Grant (FP7-PEOPLE-2013-CIG) "FundPhysicsAndLSS." A Completeness of Eq. (4.21) We now show that the, perhaps surprisingly, simple form of the evolution equation for E F ij is the single nontrivial component of Eq. (4.8) in CFC at leading order in derivatives. This implies that Eq. (4.21) consistently contains all constraints from the full covariant Einstein system, with no additional constraints at leading order in derivatives. We begin with the full covariant equation for the Weyl tensor, Eq. (4.8) : evaluated in CFC, and consider all the distinct combinations of space-time indices, up to symmetries of the Weyl tensor (which are the same as those of the Riemann tensor). Henceforth, we shall drop the "F" subscript. CFC shall be assumed throughout. The line element of the CFC metric is given by Eq. (3.1) We shall first consider the LHS of the system, The relevant Christoffel symbols of Eq. (A.2) evaluated on the geodesic are [14], Γ µ 0ν | geo = H| geo δ µ ν , Γ 0 ij | geo = H| geo δ ij (A.5) For {µ, ν, λ} = {0, 0, 0}, C 00κ0 is automatically zero given the symmetries of the Weyl tensor. For {µ, ν, λ} = {0, 0, i}, the only non-vanishing Weyl component is C 0ik0 . Since h µν = O(x 2 ) and in leading order CFC we restrict to 2 spatial derivatives, regarding spatial derivatives acting on the metric, only terms of the form ∂ 2 k h µν survive. Also, any spatial derivative of the Weyl tensor is neglected in leading order CFC. We further note that the Weyl tensor is tracefree w.r.t. any 2 indices. It follows that We shall now consider the RHS of Eq. (4.8) for the components with non-trivial LHS, adopting an ideal pressureless fluid stress-energy tensor (again, we really only require the absence of pressure perturbations This is a higher-derivative contribution, as it involves a spatial (fluid-orthogonal) derivative on the stress-energy tensor and is thus equivalent to a third spatial derivative acting on the metric. Consequently, the only apparently non-trivial component of the system, besides our equation for E ij , isĊ ij 0 At this order in derivatives, this component is entirely decoupled from the other quantities in Eq. (4.21) and corresponds to a constant that can be set to zero.