Synchronizing the Consistency Relation

We study the $N$-point function of the density contrast to quadratic order in the squeezed limit during the matter-dominated (MD) and radiation-dominated (RD) eras in synchronous gauge. Since synchronous gauge follows the free-fall frame of observers, the equivalence principle dictates that in the gradient approximation for the long-wavelength mode there is only a single, manifestly time-independent consistency relation for the $N$-point function. This simple form is dictated by the initial mapping between synchronous and local coordinates, unlike Newtonian gauge and its correspondingly separate dilation and Newtonian consistency relations. Dynamical effects only appear at quadratic order in the squeezed limit and are again characterized by a change in the local background, also known as the separate universe approach. We show that for the 3-point function the compatibility between these squeezed-limit relations and second-order perturbation theory requires both the initial and dynamical contributions to match, as they do in single-field inflation. This clarifies the role of evolution or late-time projection effects in establishing the consistency relation for observable bispectra, which is especially important for radiation acoustic oscillations and for establishing consistency below the matter-radiation equality scale in the MD era. Defining an appropriate angle and time average of these oscillations is also important for making separate universe predictions of spatially varying local observables during the RD era, which can be useful for a wider range of cosmological predictions beyond $N$-point functions.


Introduction
The inflation era in the early universe sets the initial condition of the cosmological perturbations, which seed the cosmic microwave background (CMB) anisotropies and the large scale structure (LSS) [1][2][3][4][5][6]. The cosmological perturbations are often characterized by the initial comoving curvature perturbations 1 ζ I . Their power spectrum P ζ I has now been precisely measured on the CMB scales [7], revealing a nearly scale-invariant shape. However, this information about the initial condition is essentially kinematic, following from the approximate de Sitter symmetry of inflation, and the underlying dynamics of inflation remains largely untested.
To gain a deeper understanding of the microphysics of inflation, it is essential to measure higher-order correlations or non-Gaussianity of the curvature perturbations. Such measurements provide a more direct window into the dynamical content of inflation, and can help distinguish between different inflationary models [8]. Primordial non-Gaussianity can be observed in the CMB bispectrum [9][10][11][12] and the statistics of large-scale structure including scale-dependent halo bias [13,14], and upper limits on its presence have already ruled out some multi-field and low sound speed inflationary parameter space.
In single-field inflation models, the three-point correlation function of the curvature perturbations in the squeezed limit satisfies the following consistency relation [15] (see also [16][17][18][19][20]): where δ D is the Dirac delta function and n s (k S ) − 1 = d ln k 3 S P ζ I (k S ) d ln k S , (1.2) characterizes the tilt, with k S = (k 1 − k 2 )/2 and q ≪ k 1 , k 2 . This relation stems from the fact that the long-wavelength mode, ζ I (q), acts just as the rescaling of the spatial coordinates up to O(q). This consistency relation can be extended to observable correlation functions after inflation, e.g., the squeezed bispectrum ⟨ζ I δδ⟩, which describes the correlation of the late-time density perturbations δ in the presence of the long-wavelength mode ζ I . If future measurements of the squeezed three-point function in the CMB anisotropies and the LSS discover a deviation from this evolved (late-time) consistency relation, it would rule out single-field inflation models. Previous studies of the late-time consistency relation have so far focused on Newtonian gauge [21,22] and often in the context of a separate "Newtonian consistency relation" at O(q) [23,24]. In this paper, we present a new analysis of the consistency relation in synchronous gauge, which we show has the same form as the rescaling of coordinates in the inflationary consistency relation. Since the synchronous gauge is defined as the free-fall frame of initially comoving observers, the equivalence principle requires that the consistency relations are trivial/kinematic spatial coordinate transformations up to O(q). In addition, synchronous time slicing clarifies the separate roles of initial non-Gaussianity versus the evolution of short-wavelength fluctuations in establishing the observable late-time consistency relation. 2 Moreover, the synchronous calculations at O(q 2 ) reveal the dynamical effect of the evolution of short-wavelength density perturbations in the presence of a long-wavelength mode, often approximated in the so-called separate universe approach where it is considered as a local background cosmology [25,[28][29][30][31]. This is because the separate universe construction is conceptually related to a free-fall frame and contains the dynamics of a long-wavelength adiabatic mode. Consequently, synchronous gauge has the computational advantages in the derivation of the consistency relation. Unlike in Newtonian gauge, the consistency relation in synchronous gauge can be derived with time-independent coordinate transformations, which are manifestly non-dynamical and highlight the equivalence principle by maintaining the same free-fall frame throughout the superhorizon evolution. 3 As a result, the adiabatic or gauge conditions for the long mode associated with the coordinate transformation become simple in synchronous gauge, and do not require any conditions on the matter content of the universe, such as the neglect of the anisotropic stress of the neutrinos, as is common in other techniques.
To verify that the synchronous gauge consistency relation is satisfied for single-field inflation and show how it can be violated under other types of initial non-Gaussianity, we compare it to second-order density perturbation calculations assuming a perfect fluid during matterdominated (MD) and radiation-dominated (RD) eras. Second-order perturbations have been extensively studied in the context of both the CMB [34][35][36][37][38][39][40][41][42][43][44][45][46][47] and the LSS [48][49][50][51][52][53][54][55]. In general, the second-order density perturbations can be divided into the homogeneous part, whose evolution is the same as the first-order density perturbations, and the inhomogeneous part, whose evolution is sourced by the convolution integrals of the first-order perturbations. Importantly, the homogeneous part includes the primordial non-Gaussianity, whereas the inhomogeneous part leads to a purely non-primordial, late-time non-Gaussianity due to the non-linear evolution of the perturbations. We show both the homogeneous and inhomogeneous parts are needed for the consistency relation to hold and that the inhomogeneous evolution in RD is responsible for establishing the consistency relation in MD below the horizon scale at matterradiation equality. Finally we show that separate universe approach gives the same result as the second-order calculation at O(q 2 ), if we consider an angle and time average of the latter. This paper is organized as follows. In §2, we derive the late-time consistency relation in synchronous gauge and show how its extension to O(q 2 ) can be interpreted using the separate universe approach. Then, we discuss the compatibility of these relations with the evolution of second order perturbations, first from arbitrary initial non-Gaussianity conditions, and then from single-field initial non-Gaussianity in §3. We devote §4 to the conclusion. A number of appendices contain additional technical details and derivations. In Appendix A, we review the removal of the delta function from the consistency relation. In Appendix B, we solve for second-order perturbations in synchronous gauge. In Appendix C, we present and review the analogous results in Newtonian gauge, and in Appendix D we perform gauge transformations to relate the results in Newtonian and synchronous gauges as well as the initial conditions from comoving and uniform density gauge. Appendix E collects the main variables used in this paper.

Consistency Relations and Separate Universe in Synchronous Gauge
The physical essence of the cosmological consistency relations is the equivalence principle. These relations express the fact that a local observer's experience is insensitive to particular gravitational perturbations-those that can be generated by coordinate transformations. As such, it is natural to expect that the kinematic constraints that the consistency relations capture will be most simply stated in coordinates adapted to locally freely falling observers, which are precisely those of synchronous gauge.
In this section, we derive the kinematic consequences of these spatial coordinate transformations in synchronous gauge that mimic the effect of a long-wavelength curvature perturbation of wavenumber q to O(q), and translate them into soft theorems satisfied by correlation functions of local observables. These same techniques of absorbing the long-wavelength mode into the background, now thought of as defining a separate universe or local Friedmann-Lemaître-Robertson-Walker (FLRW) background, also mimic the truly dynamical effects of the long-wavelength curvature perturbation on local observables at O(q 2 ).

Synchronous Gauge and Free-Fall Frame
We are interested in the behavior of fluctuations around a spatially flat FLRW background metric δg µν = g µν −ḡ µν , whereḡ µν is the conformally flat background metric with the scale factor a. Synchronous gauge is defined by the requirement that the perturbations of the metric satisfy δg 0µ = 0. The line element can therefore be parameterized as where the metric perturbations are solely in the spatial components. Note that the metric perturbations here denote (long-wavelength) first-order perturbations, while later in §3 we will discuss the higher-order impact on short-wavelength observables in the presence of longwavelength mode. In general, E ij contains components that transform as a scalar, vector, and tensor under the spatial rotation symmetries of the background. While the long-wavelength tensor mode in E ij induces locally observable tidal effects as soon as it enters the horizon [33,56], we focus on the dynamics of scalar perturbations (in particular the second-order scalar perturbations induced by the first-order scalar perturbations) in this paper. Then, we restrict to the scalar component of E ij , which we write as whereÊ is a 3-scalar function. 4 Throughout this work, indices of 3-vectors and 3-tensors, such as spatial vectors, momentum vectors, and their spatial derivatives, are raised or lowered by the comoving background spatial metric δ ij , which is the same convention as in [57]. Repeated spatial indices of 3-vectors and 3-tensors are likewise summed over even when they both are lowered or raised. Synchronous gauge is so named because it is the coordinate system that is established by a set of free-fall observers who initially synchronize their clocks from their initial spatial coordinate positions. As such, synchronous gauge retains the freedom to specify the time slice of the synchronization and those initial spatial coordinates, leaving the residual coordinate freedom x µ → x µ + ξ µ with Notice that the residual gauge freedom represented by the spatial functions γ(x) and h i (x) is non-dynamical, reflecting its origin in the choice of initial observers. It is advantageous to choose for the time slicing synchronous observers that are initially at rest with respect to the background expansion as a → 0, which fixes the residual temporal freedom γ(x) = 0. Since these test observers are in free fall, as long as gravity is the only relevant force this means that their frame coincides with the one that moves with the matter, namely comoving gauge where the matter stress-energy tensor is set to T 0 i = 0 (see Eq. (D. 22)). In general, this approximation holds above the Jeans scale of the matter r s (see [58] for the generalization beyond general relativity). For example, in the RD era where this scale is comparable to the horizon, r s ∼ η, corrections for a Fourier wavenumber q due to the change in slicing are suppressed by (qη) 2 . In the MD era and linear theory, the Jeans scale of the pressureless matter vanishes, and so the correction is entirely negligible. Specifically, the 3-curvatures on comoving and synchronous slicings physically coincide. In linear theory, since the background is homogeneous, they coincide at each coordinate point as well, regardless of the choice of synchronous spatial coordinates h i (x). In Fourier space, this implies that the comoving curvature perturbation ζ satisfies 5
Instead, we have split the spatial metric perturbation into the trace and traceless pieces, and absorbed all of the trace into Ψ. 5 The comoving curvature perturbation is denoted by R in some literature (e.g. Ref. [59]).
where the notation O(q n ) throughout denotes corrections of order (qr s ) n when qr s ≪ 1. Notice that in the formal limit η → 0, we have r s → 0 and the quadratic correction vanishes.
In the inflationary context, the minimum value of r s corresponds to the horizon at the end of inflation. Since this is negligible compared to wavelengths considered here, we will denote this initial value as ζ = ζ I and drop any O(q 2 ) corrections initially. In Fourier space, the anisotropic metric perturbation is expressed as While a common convention is to choose h i such that lim q→0Ê = 0 at finite ζ, for pedagogy we retain the residual gauge freedom. Without loss of generality, we can parameterize it by the fraction of the invariant curvature perturbation that is carried individually by Ψ orÊ as See Eq. (D.43) for the specific relationship between the spatial coordinate threading for different f .

Consistency Relation up to O(q)
To introduce the consistency relation, we adopt the approach used in [18,21]. For the adiabatic mode, the long-wavelength perturbation with the wavenumber q above the Jeans scale should act on short-wavelength perturbations as the coordinate transformation with deviations only at O(q 2 ). If the perturbations originate from single-field inflation, then this ansatz applies to the initial conditions for the short-wavelength perturbations as well.
In this case, the whole correlation function of the short-mode density perturbations in the presence of the long mode, represented by the initial comoving curvature perturbation ζ I , is the same as the one without the long-wavelength curvature perturbation after an appropriate coordinate transformation x →x: where the notation O(q n ) in correlators implicitly means that the Fourier modes that are included in this long-wavelength ζ I satisfy |q · (x i − x j )| ≪ 1 in addition to qr s ≪ 1, i.e., that we can Taylor expand ζ I (x i − x j ) and treat its effects as a background modulation for the short-wavelength modes order by order. Supposing that the coordinate transformation induced by ζ I is perturbatively small, we obtain to O(q) where we have used ⟨ζ I ⟩ = 0. Here ∆δ(x a , η a ) ≡ δ(x a , η a ) − δ(x a , η a ), which correlates with ζ I due to the coordinate change it represents. In Fourier space, the consistency relation can be expressed as where ∆δ(k a , η a ) ≡ d 3 x a e −ika·xa ∆δ(x a , η a ). In synchronous gauge, the task therefore is to characterize the coordinate change that models the "background wave" influence of ζ I . In the passive coordinate transform approach, the coordinate shift corresponds to a change in the small-scale density perturbation as Note that under this transformation and keeping only first-order terms in ξ i , the spatial metric transforms as In general, we seek to model this change at the level of a gradient approximation for ζ I as so that the Fourier-space consistency relation is maintained to O(q). Unlike Newtonian gauge (see Appendix C) and comoving gauge (where the translation of spatial coordinates is time and anisotropic stress dependent [60]), this coordinate change is purely spatial and time independent, but also carries the residual gauge freedom parameterized by f in Eq. (2.6). We account for the most general case by separating this coordinate change as where T denotes Ψ (the trace component) andT denotesÊ (the trace-free component). We can consider each piece independently, and then compose the most general consistency relation from this linear combination, as we shall show explicitly for the three-point function in §2.2.3. We loosely call the spatial threading with f = 0 and f = 1 as the (initially) "isotropic" and "anisotropic" cases, respectively, even though evolution makesÊ ̸ = 0 and Ψ ̸ = 0 for any f as the perturbations cross the horizon scale.

Isotropic-Synchronous Gauge (f = 0)
In this case, the initial comoving curvature perturbation is represented by the synchronous gauge metric Ψ as ζ I = −Ψ + O(q 2 ). Given the isotropy of the spatial metric, we can only allow a local rescaling of the coordinates, consisting of spatial dilation and special conformal transformations (SCT), with , it follows that ζ I transforms as from which we can associate the dilation λ and SCT b parameters that would be required to absorb ζ I into the coordinate change: λ = ζ I (0) and 2b i = ∂ i ζ I (0). Since this is purely a timeindependent spatial coordinate redefinition, this transformation behaves exactly as ζ would for the growing mode of adiabatic perturbations 6 to O(q) without further restrictions due to the gauge conditions or the matter content, including any anisotropic stress it carries (see Appendix C for the Newtonian case, where further restrictions must be taken into account). The above coordinate change then transforms the small-scale density perturbation as In Fourier space, this transformation can be expressed as Substituting this into Eq. (2.9), we obtain the consistency relation where we have neglected the O(q 2 ) contributions and used Note that the correlators on both sides of Eq. (2.19) contain momentum-conserving delta functions. Using standard operations to remove the momentum-conserving delta functions which we review in Appendix A, we obtain 6 A coordinate transformation cannot change the local number density ratios of different particle species, which would be required to mimic an isocurvature mode.
where ⟨· · ·⟩ = (2π) 3 δ D (k tot ) ⟨· · ·⟩ ′ and the differential operators correspond to the dilation and SCT generators (for the scaling weight of zero). As we shall discuss in §2.2.3, this removal of the momentum-conserving delta function does not mean that we remove momentum conservation. The consistency relation still enforces conservation, but there are various ways to send q → 0 that can be used to simplify its construction and interpretation, especially for the O(q) term.

Anisotropic-Synchronous Gauge (f = 1)
We now discuss the case of f = 1, for which ζ I = −Ê/3 + O(q 2 ). SinceÊ is an anisotropic metric perturbation, this case requires an anisotropic rescaling of coordinates to model. This leads us to consider where λ and b i are a constant scalar and vector, respectively. The derivative of ξ iT is given by This implies that ∂ i ξ iT = 0, which keepsΨ (1) = Ψ (1) = 0 up to O(q). With this transformation, the trace-free part of the spatial metric transforms as To see the transformation ofÊ, let us focus on the perturbation at finite momentum q in the gradient approximation around x i = 0, where we usedÊ(k) = A(2π) 3 δ D (k − q). We note that, for the single Fourier long mode, we can rewrite the consistency relation, Eq. (2.8), as The reason for the minus sign of −q in the right-hand-side (RHS) is that it gives a nonzero value after the ensemble average with E ij (x; q); for a spectrum of modes this is equivalent to conjugation and the condition that E ij (x) is real. To erase E ij (x; −q) with the coordinate transformation, we then equate λ = −A/3 and 2b i = iq i A/3. SinceÊ corresponds to the comoving curvature perturbation ζ in the case of f = 1 using Eq. (2.4), we find Substituting λ = 3ζ(0), 2b i = −3iq i ζ(0) and using Eq. (2.9), we obtain We then remove the delta function in both sides using a generalization of the standard procedure, which we review in Appendix A. We finally obtain 7 In the case of a general mixture of trace and trace-free pieces, the consistency relation is simply a linear combination of Eqs. (2.21) and (2.33), which we will illustrate next with the consistency relation for the three-point function.

Three-Point Consistency Relation
Let us illustrate the advantages of synchronous gauge by explicitly considering the simplest consistency relation, namely that of the 3-point function. The 3-point function in the squeezed 7 Note that the consistency relation in the anisotropic-synchronous gauge (2.33) takes the same form as the tensor consistency relation in comoving gauge [18,20], except that the traceless projector q i q j q 2 − 1 3 δij is replaced with a traceless-transverse polarization tensor. In contrast, the anisotropic scalar perturbations in synchronous gauge are traceless but not transverse. limit involves derivative operations on the 2-point function in momentum space, which require us to specify how momentum is conserved in taking the squeezed limit.
The suppression of the momentum-conserving delta functions does not mean that we break momentum conservation as q → 0 but rather that we have a choice as to what momenta to leave fixed as we take this limit. For example, consider the momentum-conserving triangles k 1 + k 2 + q = 0 in Fig. 1. We can either take q → 0 at fixed k 1 or take the midpoint approach and fix k S = (k 1 − k 2 )/2. While Eqs. (2.21) and (2.33) hold for any such choice and the resulting expressions are the same once converted to the same associations for k, the separate contributions of the O(q 0 ) and O(q) terms differ. We can use this fact to further simplify the consistency relation and highlight its physical content in synchronous gauge.
As pointed out by Refs. [18,61], originally in the context correlators in comoving gauge evaluated on the same time surface, the principal advantage of holding k S fixed is that the O(q) terms vanish, leaving the relations for absorbing ζ I in the gradient approximation exactly the same as for the homogeneous approximation. In the equal-time correlator case where η = η 1 = η 2 , the two-point function on the RHS becomes P δ (k S , η) and has symmetric k 1 and k 2 derivatives (since k S = (k 1 +k 2 )/2+O(q 2 )), so that momentum conservation guarantees the vanishing of terms with a q i k i a = 0 + O(q 2 ). For unequal-time correlators the derivation of the consistency relation in comoving gauge from coordinate transformations is more involved due to the time-dependent translation required to maintain the gauge condition of an isotropic spatial metric [60], which we will refer to as "isotropic threading", especially in the presence of neutrino anisotropic stress. 8 In synchronous gauge, this notational advantage and conceptual simplicity applies to unequal-time (η 1 ̸ = η 2 ) correlators as well, since the underlying coordinate transformation has no time dependence. The explicit, most general, three-point consistency relation in synchronous gauge to O(q) takes the form 9 lim q→0 ⟨ζ I (q)δ(k 1 , η 1 )δ(k 2 , η 2 )⟩ ′ P ζ I (q)P δ (k S , η 1 , η 2 ) = −n δ (k S , η 1 , η 2 ) + 3f µ 2 n δ (k S , η 1 , η 2 ) + (1 − 3µ 2 ) , (2.34) where µ ≡k S ·q is the angle between the short and long modes, and the unequal-time power spectrum is defined as with its (time-dependent) spectral tilt Notice that for f = 0, which is the usual choice for a fully gauge-fixed synchronous gauge, the entire consistency relation to O(q) is a dilation. Hence in this synchronous gauge, the socalled Newtonian consistency relation is also entirely subsumed into a dilation (cf. Eq. (C.15)). Specifically, the consistency relation for f = 0 can be rewritten compactly as The synchronous-gauge interpretation is that local observers cannot tell if they are moving in the gradient of the long-wavelength curvature perturbation because they are freely falling, even when comparing at different times. The local coordinates that they would establish would remove the dilation that Eq. (2.37) represents, leaving no means of locally measuring any impact of the long-wavelength mode to O(q).
In fact, if we choose f = 1 for which there is no isotropic rescaling of local spatial scales from the initial metric, the consistency relation becomes 10 Since the angle average of µ 2 is 1/3, we have where the overline indicates the angle average. The angle-averaged consistency relation then becomes trivial, and this f = 1 choice makes it clear that a local observer cannot detect the presence of a long-wavelength mode in the gradient approximation. 9 See [62] for the analogous result derived for the inflationary bispectra for comoving slicing and anisotropic threading. 10 The same anisotropic angular dependence is found in the squeezed three-point function of ζ in solid inflation [63,64].

Separate Universe at O(q 2 )
The advantage of employing synchronous gauge for adiabatic perturbations becomes even more apparent at O(q 2 ). Here the impact of the long-wavelength ζ I cannot be merely absorbed into free-fall coordinates, since it changes the local density and 3-curvature that determines short-wavelength observables (see Appendix D.3). On the other hand, this background wave still appears as a local change in the FLRW background under which short-wavelength fluctuations evolve. The isotropic impact of this background can be described by a change in cosmological parameters that is sometimes called the separate universe approach.
Here the advantage of synchronous gauge is that its coordinates are closely related to the description of separate-universe observers comoving with the local expansion, as discussed in [31,58]. In this case synchronous gauge remains the local frame established by free-fall observers, whereas comoving gauge can differ even in the choice of time-slicing if the matter fields experience non-gravitational forces.
In the separate universe construction, synchronous perturbations are reabsorbed into the local background by a change in the background density ρ SU = ρ(1 + δ L ) from the longwavelength δ L , the local scale factor a SU = a(1 − Ψ), the local FLRW curvature where we have assumed that K = 0 in the global background, and the local Hubble parameter Note that the subscript "SU" indicates a local quantity in the separate universe. The smallscale density fluctuations relative to the local background density then evolve in this separate universe in a manner given by this change in these cosmological parameters. This construction applies even if the small-scale observables are in the nonlinear regime, and has been successfully implemented in cosmological simulations [28,30,65,66]. The only difference between the separate universe and synchronous descriptions is the absorption of the synchronous metric fluctuation into the local scale factor. In the separate universe construction this would introduce an additional evolving dilation of short-wavelength scales since Ψ evolves with time at O(q 2 ) (see [30], Eq. (44)), which is absent in synchronous coordinates and moreover cannot be absorbed into the initial definition of synchronous spatial coordinates with the gauge parameter f . Effectively in the fluid mechanics analogy, synchronous spatial coordinates act as Lagrangian coordinates in contrast with the separate universe Eulerian coordinates.
In addition, the absorption of the metric fluctuation into the scale factor only applies to its spatial trace Ψ and so the impact of the anisotropic pieceÊ is eliminated by angle averaging, in the same way as we treated the f = 1 case in Eq. (2.39). In fact, the separate universe description can be generalized to encapsulate the full angular dependence by using an anisotropic background cosmology [67]. This has also been successfully implemented in simulations [68][69][70][71], but for simplicity we employ the FLRW separate universe and angleaveraged correlation functions here.

Second-Order Perturbations in Synchronous Gauge
In this section, we assess the compatibility between the consistency relation and secondorder perturbation theory results in synchronous gauge and its dependence on initial non-Gaussianity. Before going to the concrete calculations, we summarize the methodology of the second-order perturbations. See Appendix B for a detailed derivation.
For simplicity we assume here a perfect fluid where the energy-momentum tensor is expressed as where ρ is the energy density, P is the pressure, and u µ is the four-velocity. The general expression of the energy-momentum tensor with nonzero anisotropic stress is given in Appendix B. However, an advantage of the synchronous gauge consistency relation is that anisotropic stress does not introduce any conceptual differences in its construction unlike Newtonian (see Appendix C) or comoving gauges [60] (see also Eq. (D.22)). Anisotropic stress from neutrinos or radiative viscosity does of course inhibit its analytic verification using second-order perturbation theory and allows only numerical checks in general, which is why we omit it here. We expand the density perturbation up to second order as where the overbar indicates the background value and the superscripts denotes the order in perturbations. We then define the n-th order density perturbation δ (n) as Similarly, we expand the metric perturbations up to second order as, We also define the transfer functions of perturbations relative to the initial comoving curvature ζ I into the initial k-dependence and growth pieces, T m/r and D m/r , respectively: where y ≡ kη with k ≡ |k|, and the subscripts, "m" and "r", represent a quantity during the MD and RD eras, respectively. In the RD era, T r (k) = 1. In the MD era, T m (k) represents the change in the spectrum from the end of inflation to the epoch of MD (η > η m ) due to the sub-Jeans-scale evolution during the prior RD epoch. In the literature, this is usually referred to as the transfer function for the large-scale structure and has the limiting forms The detailed form of T m is not relevant for this work, but it can be computed using the standard Einstein-Boltzmann codes or be approximated with a fitting formula [72]. The density growth function D m/r then transfers these "initial" perturbations to the final perturbations at η, and its explicit form is given below. Unlike [21,22], we explicitly account for subhorizon evolution in the RD regime and its impact on separating the pieces of the consistency relation that come from initial non-Gaussianity and evolution. Consequently, the power spectrum of the density perturbation P δ is related to P ζ I by At second order, the density perturbation can be expressed as m/r is the part of the second-order kernel that evolves the second-order density field from some initial surface η m/r to η under the first-order sources. Note that we have used p 1 ·p 2 = k 2 (1−u 2 −v 2 )/2 with p 2 = k−p 1 to express the momentum dependence of the kernel I in terms of u and v. The second-order contribution at the initial surface evolves under the homogeneous equations, which are the same as in linear theory. In general, a fraction of the initial value at η m/r goes into each of the homogeneous modes as determined by matching at the boundary. Here we have assumed all of the initial contributions are in the growing mode as appropriate for taking η m/r → 0. More generally, when there are first-order-squared sources before η m/r , this expression holds for η ≫ η m/r when only the growing mode component survives, with δ (2) (k, η m/r ) as the initial value of that component at η m/r (see Appendix B for details).
In general, we can characterize the initial or homogeneous piece by specifying an initial non-Gaussianity. For the three-point study, this involves the initial three-point function. For the RD era, we can generically characterize the homogeneous term through the initial three-point function at the end of inflation or the start of radiation domination, η = η r , by assuming that all k-modes are outside the horizon at η r . For an initial non-Gaussianity obeying the single-field inflationary consistency relation (1.1), the initial three-point function also satisfies the general consistency relation (2.34), where n δ = n s + 3. Here O(q 2 ) generally contains (qη r ) 2 and for single-field inflation (q/k S ) 2 terms. It may also contain terms associated with the horizon scale at which slow-roll conditions are violated during inflation [73]. We shall see that all of these O(q 2 ) terms are parametrically smaller than those generated from the inhomogeneous sourcing at late times. For a contrasting case, where the consistency relation can be violated initially, we consider local non-Gaussianity from multi-field inflation in Appendix C.2. Likewise, we characterize the effective initial condition for the calculation in the MD era as Here η m ≳ η eq is the beginning of the MD era whereafter we omit the radiation component.
We shall see that dynamical evolution in the RD universe through the inhomogeneous term preserves the consistency relation, so that single-field inflation models satisfy where likewise we shall see that the O(q 2 ) term contains (qη m ) 2 terms from the RD evolution as well as (q/k S ) 2 terms from single-field inflation. The preceding RD epoch does change well after the equality but leaves n δ independent of η 1 and η 2 during that era (η 1 >η m and η 2 >η m ). Note that even for local non-Gaussianity where the inflationary consistency relation is violated, the dilation piece that arises from n δ (k S ) − (n s + 3) or T m will remain (see Appendix C.2.1, Eq. (C.26)).

Consistency Relation and Initial Non-Gaussianity
In the following, we will compare the second-order perturbations and the consistency relation by computing the three-point function up to O(q). Since we focus on the case where the long mode is always outside the Jeans scale, i.e., q ≪ η −1 in RD and q ≪ η −1 eq in MD, we set T m/r (q) = 1 for the long mode. We can then calculate the left-hand-side (LHS) of the consistency relation by expressing the second-order perturbations as Next, we re-express the wavenumbers k 1 , k 2 in term of k S and µ ≡q ·k S using which follows from q + k 1 + k 2 = 0. We then expand the inhomogeneous term in small q to obtain where y 1 = k S η 1 , y 2 = k S η 2 here and throughout. In the following, we will give the concrete expressions of B m/r as appropriate. Note that the O(q) consistency relation involves s ≤ 1.

Matter-Dominated Era
We first consider the MD era. The growth function of the first-order density perturbations is given by 15) and the inhomogeneous kernel of the second-order density perturbations is given by [48,54] where we have taken η m → 0 as described in Appendix B (see Eq. (B.51)). Plugging this into Eq. (3.12) and expanding in the squeezed limit, we find that the inhomogeneous contribution to the three-point function vanishes up to O(q), i.e., B [s,t] m = 0 for s = 0, 1 in Eq. (3.14). Since it immediately follows that for an initial non-Gaussianity or homogeneous term that satisfies the consistency relation Eq. (3.10), the late-time consistency relation of Eq. (2.34) holds: to O(q). In particular at this order in q, the consistency relation for f = 0 becomes a pure dilation of short-wavelength fluctuations in the long-wavelength mode, whereas for f = 1 its angle average vanishes entirely, as expected.

Radiation-Dominated Era
Next, we summarize the expressions in the RD era. The growth function during the RD era is given by where we have definedỹ ≡ y/ √ 3. We will refer to the oscillatory features in this growth function as acoustic oscillations and the change from the initial power law behavior as the mode crossing the Jeans scale. The expressions of the inhomogeneous kernel, again taking the initial surface η r → 0, can be found in Appendix B. For brevity, we here directly go to the expressions of B [s,t] r whose nonzero values for s ≤ 1 are Unlike in the MD era, the inhomogeneous part contributes to the three-point function up to O(q). This is because the RD growth function D r depends on scale so that For the consistency relation to hold at late times, it is not sufficient that the initial non-Gaussianity obeys the consistency relation Eq. (3.9): the second-order perturbations must evolve in the correct way such that the sum of the two pieces obeys dilation symmetry. Let us check this preservation of the consistency relation explicitly. At O(q), if we assume the initial piece obeys the consistency relation (3.8), then the total s = 0 contribution goes as whereas the total s = 1 contribution vanishes: The sum then implies that the three-point function satisfies the consistency relation for all times η 1 , η 2 ≥ η r . Physically, the inhomogeneous piece dilates the acoustic scale of the oscillatory small-scale perturbations, while the homogeneous or initial piece dilates the initial scale dependence. Therefore in the subsequent MD era, an inflationary non-Gaussianity that obeys the consistency relation once evolved through RD preserves the consistency relation for the "initial conditions" of the MD regime (3.11), despite the change in n δ (k S , η m , η m ) ̸ = n s +3 induced by the Jeans scale, or horizon at matter-radiation equality that is reflected in T m (k). Moreover, even if the inflationary initial condition does not obey the consistency relation, e.g., in the local non-Gaussianity case, the dilation symmetry acting on T m (k) is enforced by the RD inhomogeneous evolution, which is independent of the initial non-Gaussianity. We elaborate on these issues in Appendix C.2.2.

Separate Universe and Averaging
In this section, we verify the consistency of the O(q 2 ) three-point correlation and the separate universe approach described in §2.3. In the separate universe approach, we absorb the long-wavelength perturbation into the background, and follow the evolution of the shortwavelength perturbations in the modified background including the O(q 2 ) effects of ζ I on the spatial curvature of the separate universe. Since the separate universe is taken to be locally isotropic, this describes the angle-averaged three-point correlation [25,[28][29][30][31]74].

Matter-Dominated Era
First, let us consider the MD era case. As discussed in §2.3 (see [25,75] for details), there are in general two effects for the perturbation evolution in the separate universe approach. The first is that the changes in cosmological parameters of the local background, Eqs. (2.40) and (2.41), modify the growth of the short-wavelength density fluctuations. Converting ζ I to δ L the long-wavelength mode in real space using Eq. (3.5), we find The second is that the local density fluctuation is measured with respect to the local mean and that is perturbed as 1 + δ L , which adds 1 to the result, giving 34/21. We can then write the short-wavelength mode in the presence of the long-wavelength mode, δ L , as where the δ(k) is the density perturbation with respect to the global background and δ(k)| 0 is the density perturbation without the long-wavelength mode. Using this expression and scaling the evaluation time using D m from Eq. (3.15), we obtain Next, let us calculate the same quantity using the second-order perturbations. The inhomogeneous term contributes at O(q 2 ) from the B (3.33) The homogeneous term carries the initial non-Gaussianity and the non-Gaussianity generated during the RD era. For single-field inflation, the former is of O(y 2 1 y 2 2 q 2 /k 2 S ) due the O(q 2 /k 2 S ) corrections to H m in Eq. (3.11), and the latter is of O(y 2 1 y 2 2 q 2 η 2 m ) due to the O(q 2 η 2 m ) corrections. Both of these are negligible in comparison to the inhomogeneous contributions so long as the correlation is evaluated at late times η 1,2 ≫ η m . Converting the initial curvature and power spectra to density fluctuations using Eqs. (3.5) and (3.6) we obtain This is consistent with the separate universe approach, Eq. (3.32).

Radiation-Dominated Era
During the RD era, the separate universe approach again in principle gives a change in the growth function of small-scale perturbations D r through a change in cosmological parameters [73]. 11 However, gravitational growth of small-scale perturbations stabilizes at the Jeans scale and thereafter the growth of the long-wavelength density perturbation δ L ∝ η 2 has no further impact on the amplitude of small-scale perturbations. In this case, the shortwavelength acoustic oscillations to leading order in k S η ≫ 1 is given by whereỹ S ≡ k S η/ √ 3 and we have used D r → −4 cosỹ S from Eq. (3.19) in Eq. (3.5). Note that the local cosmological parameters do change the local conformal time and hence the frequency of the oscillation by O(δ L ) [76]. If we were to expand this relation for a small δ L , it would appear that the phase shift provides the larger effect sin(ỹ S ) × O(ỹ S δ L ) since it continues to accumulate over time. On the other hand, this shift is only observable by comparing modes at very different times. In addition, as discussed in §2.2.3 the geometric effects induced by momentum conservation q + k 1 + k 2 = 0 also alter the frequencies for the three-point correlation and accumulate over time.
We can cleanly extract the more relevant amplitude change (1 + δ L ) due to the rescaling of the local background by averaging over a cycle of the oscillation and comparing the shortwavelength modes at equal time in the three-point correlation. We therefore get where the overline denotes the cycle average of the acoustic oscillation. 11 Although the derivation of the expression in [73] was done in comoving gauge for the short-wavelength fluctuations, this becomes the same as the expression in synchronous gauge in the subhorizon limit. On the other hand the long-wavelength comoving density perturbation δL|com = 4δL| synch /3 so that the shortwavelength response in Eq. (3.36) to δL|com would no longer take the simple separate universe form.
Next, let us see the expressions of the three-point function from the second-order perturbations. In the subhorizon limit of the short-wavelength modes at equal times (η 1 = η 2 ), the O(q 2 ) contribution from the B Subleading terms in the subhorizon limit are given in Appendix B. Similarly to the MD case the contribution from the homogeneous term in Eq. (3.9) for single-field inflation is generically suppressed as 1/(k S η) 2 and (η r /η) 2 . The oscillatory terms in brackets vanish under time averaging, leaving where the double overline represents both angle and time averaging. Finally we convert the initial curvature power spectrum to the time-averaged density power spectrum using Eq. (3.6) to obtain which is consistent with the leading-order separate universe expectation, Eq. (3.36).

Conclusion
In this paper, we have thoroughly analyzed the squeezed N -point function to quadratic order in the wavenumber of the soft mode in synchronous gauge. The equivalence principle implies that synchronous gauge is the natural coordinate system to analyze the impact of longwavelength fluctuations on short-wavelength observables because it follows the free-fall frame of observers that are initially at rest with respect to the expansion. It is then manifest that there are no dynamical effects of the gradient of the long-wavelength curvature perturbation on short-wavelength modes in synchronous gauge. The entire impact of the long-wavelength mode at O(q) is a time-independent change of the local spatial coordinates, making the form of the consistency relation simple, even for late-and unequal-time correlation functions.
With the freedom to define initial spatial coordinates in synchronous gauge, for the threepoint function we can make this effect as a pure dilation or remove it entirely through angle averaging. This single relation automatically includes the so-called Newtonian consistency relation that expresses the equivalence principle in Newtonian gauge.
At quadratic order O(q 2 ), this background wave method also provides the perturbation framework for the separate universe approach. In this approach, while synchronous time slicing still defines the constant-time surfaces of the local universe, the long-wavelength mode modifies the local 3-curvature and expansion rate. We have extended these separate universe methods for the RD universe, for which special care must be taken in the averaging of acoustic oscillations with respect to angle and time.
We have also analyzed the relationship between these kinematic and dynamical constraints on the squeezed three-point function and those predicted by second-order density perturbations. In doing so, we paid special attention to the role of the dynamical evolution sourced from first-order perturbations versus initial non-Gaussianity, which we refer to as inhomogeneous vs. homogeneous contributions. In particular, we have clarified the role of the evolution during the RD era in providing the initial conditions for the MD era for scales below the horizon during radiation domination. During the RD era the initial second-order contributions are responsible for the dilation of the initial scale dependence of perturbations, whereas sourced evolution causes a dilation of the acoustic scale in the radiation and matter transfer functions. Once this combination of effects is established in the prior RD epoch, it is maintained simply by the homogeneous evolution during the MD era. We have clarified how this combination enters by explicitly keeping the matter transfer function T m throughout the calculations and the change in the tilt of the density power spectrum that it implies.
This decomposition of contributions to the consistency relation also has implications for the squeezed limit of the observable angular bispectrum of the CMB anisotropy [32,77]. In this case, the long-wavelength mode which traces the large angle temperature fluctuation dilates the angular coordinates of the small-scale temperature power spectrum, shifting their acoustic peaks from a projection of these k-space relations from the perspective of the observer today, which is not local in the plasma at recombination. One can likewise decompose its pieces into an initial piece that carries the initial tilt n s − 1 of the spectrum and a first-order sourced piece that carries the acoustic scale. Beyond single-field inflation where the former piece can break, the latter piece still contributes as we have shown explicitly for local non-Gaussianity in Appendix C and D. A breaking of this joint dilation of all scales would indicate a falsification of single-field inflation.
Beyond the consistency relation, the RD background wave and separate universe techniques introduced here are useful for the prediction of the dependence of any local quantity on the long-wavelength curvature perturbation, e.g., the formation of CMB spectral distortions [73] and the spatial dependence of the dark matter abundance relative to the radiation. We leave these and other applications to future work. and its founder Fred Kavli. W.H. was supported by U.S. Dept. of Energy contract DE-FG02-13ER41958 and the Simons Foundation. We acknowledge the use of the Mathematica package xAct [78].
We begin with the isotropic case. On the RHS of Eq. (2.19), the derivatives act on the delta function. First, let us treat the O(q 0 ) term in the RHS of Eq. (2.19). Using the relation where note again the prime means the removal of the momentum-conserving delta function, Next, let us examine the O(q) terms in Eq. (2.19). The term a 3q · ∂ ka δ D (k tot ) can be rewritten as 3N q · ∂ ktot δ D (k tot ), while the second derivative of the delta function can be expressed as In addition, Eq. (2.19) also has the first derivative of delta function multiplied by the first derivative of ⟨· · ·⟩ ′ :  From this, we can see that the term with the derivatives on the delta function in Eq. (2.32) becomes zero due to the contraction with the traceless projector ( . We can then remove the delta function without changing the O(q 0 ) expression, giving At O(q), the RHS of Eq. (2.32) has the second derivative of the delta function, which can be expressed as Also, the RHS of Eq. (2.32) has the first derivative of delta function multiplied by the first derivative of ⟨· · ·⟩ ′ : where we have used the rotational invariance of the N -point function. Using these expressions, we can rewrite the RHS of Eq. (2.32) as As

B Second-Order Relations in Synchronous Gauge
In this appendix, we provide details of the second-order calculation in synchronous gauge. We first present relevant second-order quantities in an arbitrary gauge in §B.1, and then solve the Einstein equations in the MD and RD eras for synchronous gauge in §B.2. We also provide explicit expressions for the three-point function in the squeezed limit at O(q 2 ).

B.1 Second-Order Perturbation Theory
The line element perturbed around the flat FLRW metric can be written as A perturbed quantity can be expanded as ij··· , where the overbar indicates its background value. The inverse metric components up to second order are We remind the reader of our convention that the indices of 4-tensors are raised and lowered with the metric g µν , while those of 3-tensors such as B i and C ij above are raised and lowered with δ ij . To minimize confusion, in this appendix we will lower all the 3-indices whenever possible and sum over the repeated lowered indices. It is convenient to decompose the metric components into their 3-scalar, vector, and tensor parts as Up to second order, the Einstein tensor components take the form with H ≡ a ′ /a, and the term multiplying the Kronecker delta in the spatial components G i j is given bŷ where the comma X ,i is the shorthand expression of ∂ i X. The energy-momentum tensor of a fluid with density ρ and pressure P is given by where u µ = dx µ /dτ is the 4-velocity given by up to second order, v i = ∂ iv is the comoving 3-velocity withv the velocity potential, and Σ µν denotes the anisotropic stress, which can be non-vanishing in the presence of e.g. freestreaming neutrinos. It is subject to the conditions Σ µ µ = 0, Σ µν u ν = 0, and its spatial components can be decomposed as Up to second order, we have whereρ andP denote the background quantities obeying the Friedmann equations

B.2 Synchronous Gauge
We now specialize to synchronous gauge by setting Φ = B i = 0. We focus on scalar perturbations, and therefore set all vector and tensor perturbations to zero, and relabel the spatial trace of the metric as Ψ = − 1 3 C to make contact with the notation used in the rest of the paper. This leaves Ψ andÊ as the metric perturbations with the residual gauge freedom of choosing the time slicing and spatial threading of synchronous coordinates reflecting the choice of initial free-fall observers. In synchronous gauge, the second-order Einstein equations take the form whereN ij ≡ 3 2 ∇ −2 D ij and we have used the background equations and the second-order source terms S µν consist of first-order quantities for the variables O ∈ {Ê, Ψ, δ}. The differential operators L O take the form in terms of the dimensionless variable y ≡ kη. The order of each differential operator indicates the number of pure gauge modes on top of the two physical (growing and decaying) solutions; there are two gauge modes forÊ and C, and a single gauge mode for δ. The source terms are given by Since the source terms start at second order, the first-order perturbations obey the homogeneous equations (B.30), which admit analytic solutions for a constant w in terms of the generalized hypergeometric functions. At early times, the independent solutions have the following power-law behavior: For the two metric variables, the constant mode in the curvature perturbation Ψ (1) + 1

3Ê
(1) corresponds to the growing mode in the density solution, which scales as y 2 . As discussed in §2.1 (see Eq. (2.6)), there is a residual spatial gauge freedom that is fixed by the choice of initial observers as to how to assign the curvature to the two metric perturbations individually.
In the cases where the constant curvature is assigned to one or the other entirely, f = 0 or f = 1 in the main text, then the other will scale as y 2 . The other two power-law solutions correspond to the decaying mode and the residual gauge freedom in setting the initial time surface. We shall make these assignments more concrete in the MD and RD cases below, but in general we set them to zero. We will be particularly interested in solving the equation for the second-order density perturbation in the MD and RD eras with w = 0 and w = 1/3. We will express the first-and second-order solutions in terms of the initial first-order curvature perturbation as The inhomogeneous kernel I inhom.
m/r captures the part that is generated by the dynamical evolution from the initial time surface to late times, and the homogeneous part I hom.
m/r is fixed by the initial condition at early times. In the following, we will solve the inhomogeneous solution using the Green's function, and leave the homogeneous solution to be specified by the type of initial non-Gaussianity. (1) (k, η) = cÊy 2 + dÊy −3 + eÊy −1 + fÊ ,

B.2.1 Matter-Dominated Era
We set the decaying mode of adiabatic perturbations to zero by setting d O = 0. In terms of the choice of synchronous observers, e O corresponds to the initial time slicing of synchronization, and our choice of e O = 0 specifies that the observers are at rest with respect to the background expansion. The time-independent gauge freedom f O then corresponds to their initial spatial coordinates, which we parameterize according to (2.6) as Focusing on the second-order density perturbation, we integrate the Green's function G δ over the source from η m to η to obtain where y m ≡ kη m . Notice that the lower boundary at η m (or y m ) contributes to terms that correspond to the growing mode c δ , decaying mode d δ , and temporal gauge mode e δ of linear theory (B.45) at η > η m . These homogeneously evolving contributions represent the evolution of the impact of the sources at the boundary. Conversely, with a model where there are sources at η < η m there will also be homogeneous contributions to each of these modes at the boundary. At η ≫ η m only the growing mode component of contributions from η m or before survive. Moreover if we take η m → 0 these terms formally disappear, the inhomogeneous term includes all sourced contributions, and Eq. (B.51) reduces to Eq. (3.16). The homogeneous term then represents the "true" unsourced initial conditions, which we assume are initially in the growing mode. Hereafter we implicitly assume this approach of sending the initial surface to zero time. Interestingly, the inhomogeneous solution is independent of the purely spatial gauge freedom parameterized by f , whereas the consistency condition depends on f . As discussed in the main text and shown explicitly below, the MD consistency relation comes directly from non-Gaussianity at the initial surface and not from the sourced second-order evolution.
With η m → 0, the homogeneous part I hom.
δ,m provides a growing mode solution c δ , a time-independent function of u and v, that depends on the initial second-order perturbation. While this can be fully fixed with a full specification of this solution, its implications for the three-point function at late times can be fixed by specifying the initial three-point function as given by Eq. (3.10) through H m .
We here summarize the expressions of the O(q 2 ) inhomogeneous terms of the three-point function, parameterized as where B m denotes the contribution from the inhomogeneous solution. The nonzero coefficients for s = 2 are Notice that since the inhomogeneous contributions have no s = 0, 1 components, the validity of the consistency relation relies entirely on the initial condition. If H m satisfies the consistency condition initially, as it does when this initial non-Gaussianity arises from single-field inflation, then it is trivially preserved by matter-dominated evolution. This holds because matterdominated evolution has no intrinsic scale and the consistency relation in synchronous gauge is then a pure remapping of initial scales.

B.2.2 Radiation-Dominated Era
The first-order solutions in the RD era take the form Integrating the Green's function G δ over the source from η r to η and then sending η r → 0 gives the inhomogeneous solution for the second-order density perturbation, which can be expressed in the squeezed limit as Since D r (kη) ∼ η 2 at early times, this way of organizing various terms manifests that the second-order density kernel vanishes as we approach the initial time surface at η r ≈ 0. More precisely, the above solution scales as η 4 as η → 0.
Unlike the MD case, the inhomogeneous solution in RD has a nonzero contribution at O(q 0 ) at late times. The parameterization of the three-point function is given by , where D i ≡ D r (y i ), S i ≡ỹ i sinỹ i , n s is evaluated at k S , and n ′ s = dn s (k S )/dk S . Note that the terms proportional to µ 4 vanish for the f = 0 gauge choice.

C Relations in Newtonian Gauge
In this appendix, we focus on the three-point function in Newtonian gauge in parallel to that in synchronous gauge, discussed in the main text and Appendix B. We discuss the consistency relation up to O(q) and the second-order density perturbations up to O(q 2 ) in Newtonian gauge. For the second-order perturbations, we use the results in [80], which takes the non-Gaussianity to be spatially local in uniform density gauge. Unlike in the main text and Appendix B, we do not explicitly split the contributions into the homogeneous and the inhomogeneous ones in this appendix. Instead, we use the total kernel for the second-order density perturbation. Note that, even if we separate the homogeneous and inhomogeneous contributions from the total kernel in Newtonian gauge, this does not correspond to the homogeneous and the inhomogeneous contributions in synchronous gauge. This is because there is no time-slicing invariant way to split the homogeneous and the inhomogeneous contributions. Since our synchronous condition sets its time slicing to be the same as comoving and uniform density gauge as η → 0, that split has a well-defined meaning for local non-Gaussianity (see Appendix D.1).
We take the following metric notation in Newtonian gauge: where the dots mean higher order in perturbation theory.

C.1 Dilation and Newtonian Consistency Relation
Unlike in synchronous gauge, the consistency relation in Newtonian gauge requires not only a time-independent or initial spatial coordinate transformation but also time-dependent space and time transformation, complicating both its form and an interpretation [21,22,81].
Requiring the coordinate transformation to maintain the gauge condition g 0i = 0 and g ij | i̸ =j = 0, we obtain:η Under this coordinate transformation, Φ and Ψ in Eq. (C.1) transforms as By choosingΦ =Ψ = 0, the new coordinate absorbs the long-wavelength modes of Φ and Ψ up to O(q) if However, for the new coordinates to absorb the physical mode, the parameters λ, ϵ and ξ must satisfy the adiabatic condition [82], which requires the long-wavelength mode to satisfy the Einstein equations with a finite wavenumber q. These conditions can be expressed as where recall thatv is the velocity potential and where σ is the anisotropic stress. Substituting Eqs. (C.6) and (C.7) into Eq. (C.8), we obtain v = −(ϵ + x · ξ) . (C.10) Also, the long-wavelength comoving curvature (ζ = −Ψ + Hv) can be expressed as From this, we can see λ = ζ I (x = 0) and 2b i = ∂ i ζ I (x = 0). Using Eq. (C.9), we obtain In the absence of anisotropic stress (σ = 0), we can solve for ϵ and ξ i in terms of the growth function ofv, which satisfies (C.14) Namely, Eqs. (C.12) and (C.13) lead to ϵ = −V ′ (η)ζ I and ξ i = −V ′ (η)∂ i ζ I . Then following [21,22], the consistency relation in Newtonian gauge under the perfect fluid condition is given by where 13 a , (C. 17) and recall that D and K i were defined in (2.22) and (2.23). Throughout, we put the subscript "N" when we would like to emphasize that the Newtonian gauge expressions are different from the synchronous gauge expressions. Although V has an arbitrary additive, temporallyconstant term, it does not appear in the final results because of the momentum conservation a k i a = 0. We here use K instead of D for the growth function to clarify that it is for Newtonian gauge, which are defined by 13 Note that the last term in Eq. (C.17) is omitted in [22] (their Eq. (92)).
Under the perfect fluid condition Φ (1) = Ψ (1) , the growth functions during the MD and RD eras become [80] K δ,m (y) = 6 5 + y 2 10 , Given V ′ (η) = η/5 during the MD era, D N and K N become while during the RD era we have V ′ (η) = η/3 and For the MD era, the squeezed three-point function then takes the form Here we follow the convention that n s denotes the tilt of the initial curvature power spectrum P ζ I andñ s (k) denotes its tilt in the MD era after evolution through the RD era. The correspondence to the tilt of the synchronous gauge density power spectrum in the main text isñ s (k) = n δ (k) − 3 for η > η m .
The O(q) piece of Eq. (C.26) is the so-called Newtonian consistency relation during MD. We note that this equation is exact rather than leading order in y 1,2 ≫ 1 as in Ref. [22] due to our inclusion of the last term in Eq. (C.17), and corrects typographical errors in Ref. [21] (their Eq. (62)). It also generalizes the results in both [21,22] to scales below the horizon at matter-radiation equality by including T m and the modified tiltñ s . We remark that O(q) terms are nonzero when η 1 ̸ = η 2 , while those in synchronous gauge are always zero. This highlights one of the advantages of synchronous gauge. In synchronous gauge, the time coordinate is always the proper time of free-falling observers at each spatial point. This leaves only the initial time-independent spatial coordinate choice to distinguish local and synchronous coordinates, even when comparing perturbations at different times. On the other hand, in Newtonian gauge the consistency relation is associated with the time-dependent coordinate transformation, which leads to time-dependent operators in the consistency relation. These time-dependent operators give nonzero O(q) terms for η 1 ̸ = η 2 even in the k S reference convention.
Similarly, the O(q 0 ) and O(q) terms during the RD era are given by where recallỹ = y/ √ 3. Here we have again expressed the results in terms of the tilt of the initial curvature spectrum, which is related to the tilt of the synchronous gauge density power spectrum at η r in the main text as n s = n δ (k, η r , η r ) − 3. Using the above expressions, we can express the three-point function in the largeỹ 1 andỹ 2 limit up to O(q) as lim q→0, y 1,2 ≫1 As with the MD epoch, this Newtonian consistency relation in the RD era has an unequal-time contribution at O(q).

C.2 Three-Point Function from Local Non-Gaussianity
Next, we show the expressions of the second-order density perturbations in Newtonian gauge. Unlike in Appendix B, we take a specific choice of initial non-Gaussianity in this appendix, which amounts to specifying the homogeneous term of second-order density perturbations (see Appendix D). In particular, we follow [80] in choosing the local-type non-Gaussianity where the second-order curvature in uniform density gauge ζ (2) δ at η → 0 is parameterized as [45] ζ (2) Note that b NL = 0 corresponds to the single-field, slow-roll, inflationary prediction. Strictly speaking, for this to apply in the MD epoch given an initial local non-Gaussianity from inflation, all modes should be well outside the horizon at matter-radiation equality. We shall see that applying this ansatz directly in the MD era on smaller scales, without evolving second-order perturbations through the RD era, will miss the full impact of dilation on scales where T m (k) ̸ = 1 and technically violate the consistency relation for any b NL . We express the second-order perturbation as where we have used the symbol J for the kernel, as opposed to I in synchronous gauge, to clarify the difference of the gauge. This leads to the three-point function with kernels given below (see Eqs. (C.34) and (C.43)). In general we can again expand this expression in the squeezed limit in powers of q as lim q→0 ⟨ζ I (q)δ(k 1 , η 1 )δ(k 2 , η 2 )⟩ ′ P ζ I (q)T 2 m/r (k S )P ζ I (k S ) = s,tC where we used the notationC to parameterize the total squeezed bispectrum for local non-Gaussianity in Newtonian gauge, as opposed to B used in synchronous gauge, which just includes the inhomogeneous part but applies to any type of initial non-Gaussianity. In Appendix D, we will also see that the results in this appendix once gauge transformed to synchronous gauge are consistent with those in Appendix B for the homogeneous contributions of local non-Gaussianity.

C.2.1 Matter Era
In the MD era, the second-order Newtonian-gauge kernel J is given by [54,80] 14 The second-order kernel J in this paper corresponds to A 2 uvI in [80], where A = 3/5 for the MD era and A = 2/3 for the RD era.
Taking the angle average in the limit of y 1 ≫ 1 and y 2 ≫ 1, we obtain lim q→0, y 1,2 ≫1 where we have used n δ (k, η 1 , η 2 ) =ñ s (k) + 3 + O(1/y 2 1,2 ) in Newtonian gauge. Unlike in synchronous gauge, the growth function in Newtonian gauge depends on the horizon scale and therefore the tilt of the density power spectrum depends on η 1 and η 2 , though its timedependent parts are subdominant in y 1,2 ≫ 1. In the case of a subhorizon long-wavelength mode qη ≫ 1, we can rewrite this expression as lim q→0, y 1,2 ≫1, qη≫1 Note that, except for the term proportional to n δ , this expression is the same as Eq. (3.34) for synchronous gauge. The n δ term denotes the effect of the dilation, which comes from the difference between the spatial coordinates in Newtonian and synchronous gauges, whereas the time slicing shift causes only a small effect in the MD era due to the subhorizon growth of density fluctuations. The n δ term comes from the combination of O(y 4 (v − 1)) term in Eq. (C.35) and the expansion of K δ,m (kη)T 2 m (k)P ζ I (k) with respect to k S . We will concretely see that the O(y 4 (v − 1)) term is cancelled by the spatial gauge shift from Newtonian to synchronous gauge in Appendix D.

C.2.2 Radiation Era
The homogeneous and inhomogeneous kernels during the RD era in the squeezed limit (u, |v − 1| ≪ 1) are given by where note againỹ = y/ √ 3. See [80] for the full expression of J r . The nonzero coefficients in (C.33) up to O(q) are then given bỹ K δ,r (y 1 )K δ,r (y 2 ) , K δ,r (y 1 )K δ,r (y 2 ) Notice that unlike specifying local non-Gaussianity directly in MD, here there is a selfconsistent choice that verifies the consistency relation: b NL = (1 − n s )/4. Substituting this into Eqs. (C.44) and (C.45) reproduces Eq. (C.28). This initial consistency condition is then dynamically preserved for any later η 1 , η 2 in the RD era. For this value of the initial local non-Gaussianity, evolution through the RD era then restores consistency in the MD era as well by supplying the missing conversion of b NL → (1 −ñ s (k S ))/4 that includes the dilation of T m (k S ) discussed below Eq. (C.38).
Using the above expressions, we can express the three-point function in the large y 1 and y 2 limit up to O(q) as lim q→0, y 1,2 ≫1 ⟨ζ I (q)δ(k 1 , η 1 )δ(k 2 , η 1 )⟩ ′ P ζ I (q)P ζ I (k S ) = 64(ỹ 2 cosỹ 1 sinỹ 2 +ỹ 1 sinỹ 1 cosỹ 2 ) 3 which is the same as Eq. (C.29). Notice that the initial consistency relation plays no role in this y 1,2 ≫ 1 limit. As discussed at the end of §3.1, in the RD era the effect of dilation on the scale of the sinusoidal acoustic oscillations dominates over that of the initial power spectrum, though both pieces are explicitly verified if the consistency condition holds initially using the exact coefficients above.
Expanding the coefficients at O(q 2 ) as the nonzero components for ℓ = 0 arẽ C [2,0] r,0 (y 1 , y 2 ) = and the nonzero components for ℓ > 0 arẽ C [2,0] r,1 (y 1 , y 2 ) =C [2,2] r,2 (y 1 , y 2 ) where K i ≡ K δ,r (y i ) andS i ≡ sinỹ i . Then, the time-angle averaged correlation function with Similar to the MD era case, the origin of the (n s − 1)-dependent term is the O((v − 1)ỹ 2 cosỹ) term in Eq. (C.43). Apart from that, the (n s − 1)-independent term is also different from that in synchronous gauge, Eq. (3.38). In Appendix D, we will also see that the (n s − 1) dependent or independent terms in the RD era are modified with the spatial and time shift, respectively, in the transformation from Newtonian to synchronous gauge.

D Relations Between Gauges
In this appendix, we show how the various consistency relations and separate universe expressions are related in different gauges via a second-order gauge transformation. We begin in §D.1 with the active view of gauge transformations at second order and discuss their interpretation as a change in time slicing and spatial threading associated with the chosen gauge constraints. In §D.2, we then transform the second-order Newtonian gauge treatment for local non-Gaussianity in Appendix C to synchronous gauge. By doing so, we gain geometric intuition as to why the Newtonian gauge consistency relations are simplified in synchronous gauge. In addition, this provides a cross-check for the general synchronous gauge derivation in Appendix B, which is valid for any type of initial non-Gaussianity. Finally in §D.3, we relate the curvature perturbation to the 3-geometry of the time slicing and show how initial non-Gaussianity in uniform density gauge and the various spatial coordinatizations of synchronous gauge are related.

D.1 Second-Order Gauge Transformations
Let us first summarize the gauge transformation properties of the quantities that we use (see Ref. [57] for a detailed derivation). Unlike in the main text, we take the active approach for interpreting ξ µ as the generator of the gauge transformation. A (tensorial) quantity T transforms as where the tilde denotes a gauge-transformed quantity and £ ξ is the Lie derivative with respect to ξ µ . Hereafter, we neglect third-order perturbations. Then, the perturbation of the quantity is given by Similarly, the generating vector ξ µ = (α, β ,i ) can be expanded as ,i + Up to second order, perturbations transform under (D.1) as The Lie derivatives with respect to ξ µ of a scalar φ, a vector v µ , and a tensor t µν are given by Since the energy density ρ is a scalar and the backgroundρ does not transform, the density contrast δ = ρ/ρ obeys the following scalar transformation rule: The 4-velocity vector u µ obeys the vector transformation rule. From the i-component, we can obtain the transformation of the velocity potential, v (1) The metric obeys the tensor transformation rule. The metric perturbations are given by g µν =ḡ µν + δg µν . Then, we obtain δg µν = δg µν +ḡ µν,λ ξ λ +ḡ µλ ξ λ ,ν +ḡ λν ξ λ ,µ + δg µν,λ ξ λ + δg µλ ξ λ ,ν + δg λν ξ λ ,µ + 1 2ḡ µν,λα ξ λ ξ α + 1 2ḡ µν,λ ξ λ ,α ξ α +ḡ µλ,α ξ α ξ λ ,ν +ḡ λν,α ξ α ξ λ ,µ +ḡ λα ξ λ ,µ ξ α ,ν From the first-order part of the 0i and ij components, we get 15) and the first and the second order parts of the 00 components give To go from Newtonian to synchronous gauge, we must erase Φ while keeping B ,i = 0. From Eq. (D.16), we can then relate Φ in Newtonian gauge to the gauge transformation parameter α as By imposing B ,i = 0 in both side of Eq. (D.13), we can also relate α to β at first order as where C e is an integration constant. Before going to the concrete gauge transformation from Newtonian to synchronous, let us briefly comment on uniform density and comoving gauges here. The uniform density gauge condition at first order is given by Meanwhile, the comoving gauge condition at first order is given bŷ To go from synchronous to uniform density gauge, the first-order gauge transformation parameters need to satisfy where the subscript "S" denotes quantities in synchronous gauge. Because δ S ∝ (qη) 2 for the growing mode outside the horizon, we have α (1) = 0 up to O(q). Using (D.10) we can then see that α (2) = 0 to O(q) as well, for the synchronous f = 0 gauge whereÊ (1) S = O((qη) 2 ). Since the time slicing is independent of the induced density perturbation from the β (1) change in the spatial threading between surfaces, the time slicing of synchronous gauge for any f and uniform density gauge coincide in second-order perturbation theory in this limit as well. We shall see in Eqs. (D.51) and (D.96) that this also has the consequence that ζ δ = −Ψ for f = 0, since the intrinsic 3-geometries of the same time slice must agree. The equivalence of uniform density and comoving gauge slicing and ζ δ = ζ follows the same logic because the comoving gauge density perturbation scales as (qη) 2 as well.
Synchronous gauge differs from uniform density and comoving gauges in the spatial threading of the time slices. To move from synchronous gauge to comoving gauge, the firstorder gauge transformation parameter should satisfy We can again see that the time slicing remains the same to leading order in the superhorizon limit ifv (1) S = 0 initially. For both uniform density and comoving gauges, the spatial coordinates differ due to the isotropic conditionÊ (1) = 0. SinceÊ S depends on ∇ 2 σ from Eq. (B.22), the above gauge parameter β depends on σ itself. This indicates that, in order to obtain the isotropic threading of spatial coordinates in uniform density or comoving gauges in the presence of the time-dependent anisotropic stress by a coordinate change, we need to perform the time-dependent spatial coordinate change β (1) that depends on σ at leading O(q 0 ). This complicates the derivation and interpretation of the consistency relation in uniform density or comoving gauges for unequal-time correlators (see also the discussion in the beginning of §2.2.3).

D.2 Second-Order Kernel from Local Non-Gaussianity
We now derive second-order results for local non-Gaussianity in synchronous gauge by using the results of Newtonian gauge and performing the gauge transformation between the two. By doing so we also extract the inhomogeneous kernel in synchronous gauge which is independent of the initial non-Gaussianity.

D.2.1 Matter-Dominated Era
We first discuss the gauge transformation in the MD era. The synchronous gauge conditions, Eqs. (D.17) and (D.18), leave two free temporal integration constants (free spatial functions) that correspond to the initial choice of synchronous time slicing and spatial coordinates. Throughout, we fix the initial time slicing by imposing that lim η→0 α = 0, which corresponds to e O = 0 in Appendix B. For the spatial coordinate gauge fixing, we first take lim η→0 β = 0 which corresponds to taking fÊ = 0 in Appendix B or equivalently f = 0. 15 For this case, we obtain from Eqs. (D.17) and (D.18), Also, we are neglecting the decaying mode in Newtonian gauge, which sets d O = 0 in synchronous gauge and leaves only the growing mode. During the MD era, we can rewrite Eq. (D.10) as At first order, this gives where we here denote the quantity in synchronous gauge by the subscript "S" and the above expression is consistent with Eq. (B.48). At second order, we can express Eq. (D. 25) in Fourier space as , (D.27) with u = p 1 /k, v = p 2 /k. Note that we have omitted the time arguments, and we will also sometimes omit the wavenumber arguments when they are not important. Here and in the rest of this appendix, we use the shorthand convention to denote implicit convolution. The U and V terms correspond to the terms proportional to α (1) and β (1) in Eq. (D.25), respectively, and are explicitly given by Let us here determine α (2) . Using Eq. (D.17), we can obtain In Newtonian gauge, Φ (2) is related to Ψ (2) as [80] ,j ) . (D. 33) In Fourier space, we have and Ψ (2) during the MD era is given by [80] Ψ (2) Combining these expressions, we get Plugging these into Eq. (D.31), we get the following equation for α (2) : whose solution is given by where we have set the integration constant to zero, consistent with our choice of initial time slicing for synchronous observers. From this expression, we can see that J α,m is subdominant in the limit y ≫ 1 compared to V m , which has O(y 4 ) terms (see Eq. (D.30)). Combining these expressions, we finally obtain the kernel of δ The O(y 4 ) terms, which determine the O(q 2 ) expressions of the three-point function in the late time limit (y ≫ 1), appear only in J m and V m (see Eq. (C.41) for J m ). In particular, the O(y 4 (v − 1)) term in V m cancels the O(y 4 (v − 1)) in J m . More specifically, the n δ terms in the three-point function in the Newtonian gauge at O(q 2 ) in Eq. (C.42) are cancelled by the gauge transformation term coming from β ,k δ ,k in Eq. (D.25). In other words, the appearance of the O(q 2 ) term as the dilation effect in Newtonian gauge originates from the difference between the spatial coordinates in Newtonian and synchronous gauges induced by the first-order long-wavelength mode, as anticipated below Eq. (C.42).
We have so far derived the kernel in synchronous gauge with f = 0 by the gauge transformation from Newtonian gauge. We here perform the gauge transformation from f = 0 to a general f in synchronous gauge, which leads toÊ = 0 → −3f ζ I in the superhorizon limit.
From Eq. (D.14), this gauge transformation is just the spatial shift given by From Eq. (D.10), we see that the second-order density is transformed as ,i . (D. 44) In Fourier space, we can approximate the gauge shift between synchronous gauges as 2β (1) ,i δ where we have focused on the kernel for ζ I (q)δ (2) (k 1 )δ (1) (k 2 ) with β(q), u = q/k 1 , and v = k 2 /k 1 . The kernel from the gauge term then becomes up to O(u 2 ), where we have symmetrized the kernel with respect to u and v in the first line. The second-order kernel for general f is then given by Using this kernel, we finally obtain the following expression for the local non-Gaussianity ansatz in synchronous gauge: The nonzero coefficients ofQ [s,t] in s ≤ 2 arẽ f 800 (−105 − 15n 2 δ + n 3 δ − 14k S n ′ δ + n δ (71 + 3k S n ′ δ + k 2 S n ′′ δ ))y 2 1 y 2 2 .
(D. 49) Notice again that if we mimic the single-field initial condition by replacing b NL → (4 − n δ (k S ))/4, thenQ [0,t] terms become consistent with Eq. (2.34), though b NL is constant in k S in the local non-Gaussianity model (C.30). For a more general value of b NL , there remains a O(q 0 ) piece that does not simply reflect the relationship between local and synchronous initial spatial coordinates. Instead, the squeezed bispectrum represents a modulation of the amplitude of the local short-wavelength power spectrum in the long-wavelength mode, and produces locally observable effects such as the scale dependence of halo bias [13]. Moreover, there are no inhomogeneous contributions up to O(q) in the MD era, so all effects are determined by the choice of an initial non-Gaussianity. Finally, with these full expressions for local non-Gaussianity in synchronous gauge via gauge transformation from Newtonian gauge, we can compare the results to Appendix B, for the same initial non-Gaussianity. To make contact with those results, we first need to extract the initial non-Gaussianity and define the homogeneous term at some initial surface η m . Using Eq. (D.47) and imposing I inhom. m = 0 at η m → 0, we can find the homogeneous part as Notice that, in η m → 0, the homogeneous piece is purely in the growing mode D m since all of the sourced terms are grouped into the inhomogeneous term. We have also confirmed that this expression of the homogeneous kernel can be obtained directly by using the expressions of the uniform density second order curvature in synchronous gauge. Specifically, c where the first line comes from the equivalence of the time slicing between these gauges and that the spatial coordinates are chosen in the same way to make the spatial metric isotropic, i.e., isotropic threading. Subtracting the homogeneous piece from Eq. (D.47) returns the inhomogenous result from Appendix B, and demonstrates the consistency of these second-order relations between synchronous gauge and Newtonian gauge via gauge transformation. Finally, the explicit homogeneous bispectrum terms, or equivalently H m , can be extracted fromQ − B for local non-Gaussianity

D.2.2 Radiation-Dominated Era
Next, we discuss the gauge transformation in the RD era. Similar to the MD case, we take this in two steps. We first perform the gauge transformation to the synchronous gauge with f = 0, and then to the one with a general f . From Eqs. (D.17) and (D.18), we can obtain We have imposed α, β → 0 in the limit of η → 0 to determine the integration constants, similar to the MD era case. Note again these conditions correspond to taking d O = e O = f = 0 in Appendix B. During the RD era, we can rewrite Eq. (D.10) as From the first-order equation, we obtain where U and V correspond to the terms proportional to α (1) and β (1) in Eq. (D.56), respectively, and are given by In the RD era, Φ (2) is related to Ψ (2) as [80] Φ (2) = Ψ (2) + 4(Φ (1) ) 2 − F (2) r , (D.62) with the corresponding Fourier space expression For Ψ (2) , we have where J Ψ,r (u, v, y) in the squeezed limit (u ≪ 1, |v − 1| ≃ O(u)) is given by [80] J Ψ,r (u, v, y) ≃ 4 −8 − 3b NL + 2ỹ 2 sinỹ + (8 + 3b NL )ỹ cosỹ 3ỹ 3 Combining these expressions, we get This allows us to rewrite Eq. (D.60) as Solving this equation, we obtain where we have again set the integration constant to zero to focus on the growing mode.
Combining these expressions, we finally obtain the kernel of δ S for f = 0: For the time-and angle-averaged three-point function ⟨ζδδ⟩ at O(q 2 ) in the limit ofỹ ≫ 1, only the − 4 η αδ in Eq. (D.56) modifies the term independent of (n s − 1) in the Newtonian result Eq. (C.52), and the δ ,k β ,k in Eq. (D.56) cancels the (n s − 1) term in Eq. (C.52). This can be easily checked by calculating the time-and angle-averaged three-point function only with the two terms, − 4 η αδ and δ ,k β ,k . Similar to Eq. (D.46), the kernel for the gauge transformation to a general f in the RD era is given by The kernel for a general f is then given by With this kernel, the local non-Gaussianity ansatz in synchronous gauge becomes ⟨ζ I (q)δ(k 1 , η 1 )δ(k 2 , η 2 )⟩ ′ P ζ I (q)P ζ I (k S ) = ∞ s,t=0Q , where D i ≡ D r (y i ), S i ≡ỹ i sinỹ i , n ′′ s = d 2 n s (k S )/dk 2 S , and we have considered the limit of η r /η 1,2 ≪ 1 for simplicity, as we did in §3.
To see the correspondence to Appendix B, similar to the MD era case, we here split the kernel as I r = I hom. r + I inhom.
r . Imposing I inhom. = 0 at η r → 0 in Eqs. (D.73)-(D.75), we obtain the following expression of the homogeneous kernel: where we have also confirmed that this expression of the homogeneous kernel can be obtained directly by using the expressions of the uniform density second order curvature in synchronous gauge, similar to the MD era case. Also, after some calculations, one can find that the inhomogeneous part (I inhom.

D.3 Second-Order Curvature and 3-Geometry
In this section, we clarify the relation between the intrinsic spatial curvature of the constant time slice R 3 and the curvature perturbation ζ δ at second order, and examine the spatial gauge dependence of these quantities. Note that the intrinsic 3-geometry depends only on the choice of time slicing, which we have shown is initially the same for synchronous, comoving and uniform density gauge as η → 0, but its spatial coordinatization and representation in terms of the trace and tracefree part of the spatial metric can differ.
For the separate universe approach, the intrinsic spatial curvature from long-wavelength mode gives the local FRW curvature in Eq. (2.40) As usual, the first-order intrinsic curvature is given by the Laplacian of the curvature perturbation −Ψ− 1 3Ê . In synchronous gauge, the second-order intrinsic curvature can be expressed as  where we have also shown the full y-dependent part. For f ̸ = 0, we thus see that R 3 only differs due to the change in spatial coordinates assigned to the same invariant curvature, in the same way as for the density perturbation itself, hence all of the dependence on f can be absorbed into c (2) δ . We can also compare the above to the second-order curvature perturbation in uniformdensity gauge. Note the quantities that are associated with the expansion of 3-curvature related metric quantities at first order are no longer directly related to the 3-curvature at second order. In this case, the "curvature perturbation" ζ δ is defined to be the trace of the spatial metric in uniform density gauge. From perturbations defined with alternate gauge conditions of B =Ê = 0, it reads in the superhorizon limit [57,83] lim y→0 ζ (2) In synchronous gauge where B = 0, the density perturbation scales as y 2 , and thus does not survive in the y → 0 limit. In the same limit, the choice of f = 0 for the threading setsÊ = 0.