Extending the domain of validity of the Lagrangian approximation
MMon. Not. R. Astron. Soc. , 1–38 () Printed 13 November 2018 (MN L A TEX style file v2.2)
Extending the domain of validity of the Lagrangianapproximation
Sharvari Nadkarni-Ghosh (cid:63) and David F. Chernoff † Department of Physics, Cornell University, Ithaca, NY 14853 USA Department of Astronomy, Cornell University, Ithaca, NY 14853 USA
ABSTRACT
We investigate convergence of Lagrangian Perturbation Theory (LPT) by analysingthe model problem of a spherical homogeneous top-hat in an Einstein-deSitter back-ground cosmology. We derive the formal structure of the LPT series expansion, workingto arbitrary order in the initial perturbation amplitude. The factors that regulate LPTconvergence are identified by studying the exact, analytic solution expanded accord-ing to this formal structure. The key methodology is to complexify the exact solution,demonstrate that it is analytic and apply well-known convergence criteria for powerseries expansions of analytic functions. The “radius of convergence” and the “time ofvalidity” for the LPT expansion are of great practical interest. The former describesthe range of initial perturbation amplitudes which converge over some fixed, futuretime interval. The latter describes the extent in time for convergence of a given ini-tial amplitude. We determine the radius of convergence and time of validity for a fullsampling of initial density and velocity perturbations.This analysis fully explains the previously reported observation that LPT fails topredict the evolution of an underdense, open region beyond a certain time. It alsoimplies the existence of other examples, including overdense, closed regions, for whichLPT predictions should also fail. We show that this is indeed the case by numericallycomputing the LPT expansion in these problematic cases.The formal limitations to the validity of LPT expansion are considerably morecomplicated than simply the first occurrence of orbit crossings as is often assumed.Evolution to a future time generically requires re-expanding the solution in overlappingdomains that ultimately link the initial and final times, each domain subject to itsown time of validity criterion. We demonstrate that it is possible to handle all theproblematic cases by taking multiple steps (LPT re-expansion).A relatively small number ( ∼
10) of re-expansion steps suffices to satisfy the timeof validity constraints for calculating the evolution of a non-collapsed, recombination-era perturbation up to the current epoch. If it were possible to work to infinite La-grangian order then the result would be exact. Instead, a finite expansion has finiteerrors. We characterise how the leading order numerical error for a solution generatedby LPT re-expansion varies with the choice of Lagrangian order and of time step size.Convergence occurs when the Lagrangian order increases and/or the time step sizedecreases in a simple, well-defined manner. We develop a recipe for time step controlfor LPT re-expansion based on these results.
Key words: cosmology: theory – large-scale structure of Universe. (cid:63)
E-mail: [email protected] † E-mail: chernoff@astro.cornell.educ (cid:13)
RAS a r X i v : . [ a s t r o - ph . C O ] S e p Sharvari Nadkarni-Ghosh and David F. Chernoff
Understanding the non-linear growth of structure in an expanding universe has been an active area of research for nearlyfour decades. Simulations have been instrumental in illustrating exactly what happens to an initial power spectrum of smallfluctuations but analytic methods remain essential for elucidating the physical basis of the numerical results. Perturbationtheory, in particular, is an invaluable tool for achieving a sophisticated understanding.The Eulerian and Lagrangian frameworks are the two principal modes of description of a fluid. The fundamental dependentvariables in the Eulerian treatment are the density ρ ( x , t ) and velocity v ( x , t ) expressed as functions of the grid coordinates x and time t , the independent variables. In perturbation theory the dependent functions are expanded in powers of a smallparameter. For cosmology that parameter typically encodes a characteristic small spatial variation of density and/or velocitywith respect to a homogeneous cosmology at the initial time. As a practical matter, the first-order perturbation theorybecomes inaccurate when the perturbation grows to order unity. Subsequently one must work to higher order to handle thedevelopment of non-linearity (see Bernardeau et al. 2002 for a review) or adopt an alternative method of expansion.In the Lagrangian framework, the fundamental dependent variable is the physical position of a fluid element or particle(terms used interchangeably here). The independent variables are a set of labels X , each of which follows a fluid element, andthe time. Usually X is taken as the position of the element at some initial time but other choices are possible. In any case,the physical position and velocity of a fluid element are r = r ( X , t ) and ˙ r ( X , t ), respectively. Knowledge of the motion ofeach fluid element permits the full reconstruction of the Eulerian density and velocity fields. In cosmological applications ofLagrangian perturbation theory (LPT), just like Eulerian perturbation theory, the dependent variables are expanded in termsof initial deviations with respect to a homogeneous background. The crucial difference is that the basis for the expansionis the variation in the initial position and position-derivative not the variation in the initial fluid density and velocity. TheEulerian density and velocity may be reconstructed from knowledge of the Lagrangian position using exact non-perturbativedefinitions. A linear approximation to the displacement field results in a non-linear expression for the density contrast. TheLagrangian description is well-suited to smooth, well-ordered initial conditions; a single fluid treatment breaks down onceparticle crossings begin, caustics form and the density formally diverges.First-order LPT was originally introduced by Zel’Dovich (1970) to study the formation of non-linear structure in cos-mology. In his treatment the initial density field was taken to be linearly proportional to the initial displacement field (the“Zeldovich approximation”). These results were extended by many authors (Moutarde et al. 1991; Buchert 1992; Bouchet et al.1992; Buchert & Ehlers 1993; Buchert 1994; Munshi et al. 1994; Catelan 1995; Buchert 1995; Bouchet et al. 1995; Bouchet1996; Ehlers & Buchert 1997). The work pioneered by Bouchet focused on Zeldovich initial conditions and established thelink between LPT variables and statistical observables. The work by Buchert as well as the paper by Ehlers & Buchert (1997)formalised the structure of the Newtonian perturbative series for arbitrary initial conditions. A general relativistic version ofthe Zeldovich approximation was developed by Kasai (1995) and other relativistic descriptions of the fluid in its rest framewere investigated by Matarrese & Terranova (1996) and Matarrese, Pantano & Saez (1993, 1994). LPT has been used formany applications including, recently, the construction of non-linear halo mass functions by Monaco (1997) and Scoccimarro& Sheth (2002).Not much has been written about the convergence of LPT although LPT expansions are routinely employed. Sahni &Shandarin (1996) pointed out that the formal series solution for the simplest problem, the spherical top-hat, did not convergefor the evolution of homogeneous voids. Figure 1 illustrates the conundrum that the LPT approximations diverge from theexact solution in a manner that worsens as the order of the approximation increases. The details will be described in the nextsection.This paper explores LPT convergence for the spherical top-hat and identifies the root cause for the lack of convergence.The analysis naturally suggests a means of extending the range of validity of LPT. This generalisation of LPT guaranteesconvergence to the exact solution of the model problem at all times prior to the occurrence of the first caustic.Tatekawa (2007) attempted to treat the divergence by applying the Shanks transformation to the LPT series. Althoughnon-linear transformations can sum a divergent series, the correct answer is not guaranteed; comparison of several differentmethods is usually necessary to yield trustworthy results. Other approaches include the Shifted-Time-Approximation (STA)and Frozen-Time-Approximation (FTA) which have been investigated by Karakatsanis, Buchert & Melott (1997). Theseschemes modify lower order terms to mimic the behavior of higher order terms and/or extend the range of applicability intime. None of these techniques are considered here.The organisation follows: § § § § § c (cid:13) RAS, MNRAS000
Understanding the non-linear growth of structure in an expanding universe has been an active area of research for nearlyfour decades. Simulations have been instrumental in illustrating exactly what happens to an initial power spectrum of smallfluctuations but analytic methods remain essential for elucidating the physical basis of the numerical results. Perturbationtheory, in particular, is an invaluable tool for achieving a sophisticated understanding.The Eulerian and Lagrangian frameworks are the two principal modes of description of a fluid. The fundamental dependentvariables in the Eulerian treatment are the density ρ ( x , t ) and velocity v ( x , t ) expressed as functions of the grid coordinates x and time t , the independent variables. In perturbation theory the dependent functions are expanded in powers of a smallparameter. For cosmology that parameter typically encodes a characteristic small spatial variation of density and/or velocitywith respect to a homogeneous cosmology at the initial time. As a practical matter, the first-order perturbation theorybecomes inaccurate when the perturbation grows to order unity. Subsequently one must work to higher order to handle thedevelopment of non-linearity (see Bernardeau et al. 2002 for a review) or adopt an alternative method of expansion.In the Lagrangian framework, the fundamental dependent variable is the physical position of a fluid element or particle(terms used interchangeably here). The independent variables are a set of labels X , each of which follows a fluid element, andthe time. Usually X is taken as the position of the element at some initial time but other choices are possible. In any case,the physical position and velocity of a fluid element are r = r ( X , t ) and ˙ r ( X , t ), respectively. Knowledge of the motion ofeach fluid element permits the full reconstruction of the Eulerian density and velocity fields. In cosmological applications ofLagrangian perturbation theory (LPT), just like Eulerian perturbation theory, the dependent variables are expanded in termsof initial deviations with respect to a homogeneous background. The crucial difference is that the basis for the expansionis the variation in the initial position and position-derivative not the variation in the initial fluid density and velocity. TheEulerian density and velocity may be reconstructed from knowledge of the Lagrangian position using exact non-perturbativedefinitions. A linear approximation to the displacement field results in a non-linear expression for the density contrast. TheLagrangian description is well-suited to smooth, well-ordered initial conditions; a single fluid treatment breaks down onceparticle crossings begin, caustics form and the density formally diverges.First-order LPT was originally introduced by Zel’Dovich (1970) to study the formation of non-linear structure in cos-mology. In his treatment the initial density field was taken to be linearly proportional to the initial displacement field (the“Zeldovich approximation”). These results were extended by many authors (Moutarde et al. 1991; Buchert 1992; Bouchet et al.1992; Buchert & Ehlers 1993; Buchert 1994; Munshi et al. 1994; Catelan 1995; Buchert 1995; Bouchet et al. 1995; Bouchet1996; Ehlers & Buchert 1997). The work pioneered by Bouchet focused on Zeldovich initial conditions and established thelink between LPT variables and statistical observables. The work by Buchert as well as the paper by Ehlers & Buchert (1997)formalised the structure of the Newtonian perturbative series for arbitrary initial conditions. A general relativistic version ofthe Zeldovich approximation was developed by Kasai (1995) and other relativistic descriptions of the fluid in its rest framewere investigated by Matarrese & Terranova (1996) and Matarrese, Pantano & Saez (1993, 1994). LPT has been used formany applications including, recently, the construction of non-linear halo mass functions by Monaco (1997) and Scoccimarro& Sheth (2002).Not much has been written about the convergence of LPT although LPT expansions are routinely employed. Sahni &Shandarin (1996) pointed out that the formal series solution for the simplest problem, the spherical top-hat, did not convergefor the evolution of homogeneous voids. Figure 1 illustrates the conundrum that the LPT approximations diverge from theexact solution in a manner that worsens as the order of the approximation increases. The details will be described in the nextsection.This paper explores LPT convergence for the spherical top-hat and identifies the root cause for the lack of convergence.The analysis naturally suggests a means of extending the range of validity of LPT. This generalisation of LPT guaranteesconvergence to the exact solution of the model problem at all times prior to the occurrence of the first caustic.Tatekawa (2007) attempted to treat the divergence by applying the Shanks transformation to the LPT series. Althoughnon-linear transformations can sum a divergent series, the correct answer is not guaranteed; comparison of several differentmethods is usually necessary to yield trustworthy results. Other approaches include the Shifted-Time-Approximation (STA)and Frozen-Time-Approximation (FTA) which have been investigated by Karakatsanis, Buchert & Melott (1997). Theseschemes modify lower order terms to mimic the behavior of higher order terms and/or extend the range of applicability intime. None of these techniques are considered here.The organisation follows: § § § § § c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation nd th th st rd th th th (cid:72) t (cid:76) b (cid:72) t (cid:76) Figure 1.
The time-dependent scale factor b of an initial spherical top-hat perturbation is plotted as a function of the background scalefactor a . The perturbation is a pure growing mode, i.e. the density and velocity perturbations vanish at t = 0. The black dotted line isthe exact solution. The smooth blue lines are the LPT results obtained by working successively to higher and higher order. Series witheven (odd) final order lie below (above) the exact solution. Roughly speaking, LPT converges only for a < ∼ .
2. Beyond that point thehigher order approximations deviate from the exact solution more than lower order ones . the time of validity may be extended by re-expanding the solution in overlapping domains that ultimately link the initial andfinal times, each domain subject to an individual time of validity criterion. The feasibility of this method is demonstrated insome examples. § This section describes the governing equations, the initial physical conditions, the formal structure of the LPT series solutionand the order-by-order solution.
Consider evolution on sub-horizon scales after recombination in a matter-dominated universe. A Newtonian treatment ofgravity based on solving Poisson’s equation for the scalar potential and on evaluating the force in terms of the gradient ofthe potential gives an excellent approximation for non-relativistic dynamics. When there are no significant additional forceson the fluid element (e.g. pressure forces) then it is straightforward to eliminate the gradient of the potential in favour of ¨ r ,the acceleration. The governing equations are ∇ x · ¨ r = − πGρ ( x , t ) (1) ∇ x × ¨ r = 0 (2)where ρ ( x , t ) is the background plus perturbation density, G is Newton’s gravitational constant and ∇ x is the Eulerian gradientoperator. In the Lagrangian treatment, the independent variables are transformed ( x , t ) → ( X , t ) and the particle position r = r ( X , t ) adopted as the fundamental dependent quantity. For clarity note that x refers to a fixed Eulerian grid not acomoving coordinate. The starting physical configuration is a compensated spherical perturbation in a homogeneous background cosmology. Theperturbation encompasses a constant density sphere about the centre of symmetry and a compensating spherical shell. Theshell that surrounds the sphere may include vacuum regions plus regions of varying density. Unperturbed background extendsbeyond the outer edge of the shell. Physical distances are measured with respect to the centre of symmetry. At initial time t the background and the innermost perturbed spherical region (hereafter, “the sphere”) have Hubble constants H and H p ,and densities ρ and ρ p , respectively. Let r b, ( r p, ) be the physical distance from the centre of symmetry to the inner edge ofthe background (to the outer edge of the sphere) at the initial time. Let a , b be the initial scale factors for the background c (cid:13) RAS, MNRAS , 1–38
Sharvari Nadkarni-Ghosh and David F. Chernoff and the sphere respectively. Two sets of Lagrangian coordinates Y = r b, /a and X = r p, /b are defined. A gauge choice sets a = b . Appendix A provides a figure and gives a somewhat more detailed chain of reasoning that clarifies the constructionof the physical and Lagrangian coordinate systems. The initial perturbation is characterised by the independent parameters δ = ρ p ρ − δ v = H p H − . (3)Finally, assume that the background cosmology is critical Ω = 1. The perturbed sphere hasΩ p = 1 + δ (1 + δ v ) . (4)The physical problem of interest here is the future evolution of an arbitrary initial state unconstrained by the past history. Ingeneral, the background and the perturbation can have different big bang times. Initial conditions with equal big bang timeswill be analysed as a special case of interest and imply an additional relationship between δ and δ v .While the previous paragraphs summarise the set up, they eschew the complications in modelling an inhomogeneoussystem in terms of separate inner and outer homogeneous universes. For example, matter motions within the perturbed innerregion may overtake the outer homogeneous region so that there are problem-specific limits on how long solutions for thescale factors a ( t ) and b ( t ) remain valid. The appendix shows that there exist inhomogeneous initial configurations for whichthe limitations arising from the convergence of the LPT series are completely independent of the limitations associated withcollisions or crossings of inner and outer matter-filled regions. A basic premise of this paper is that it is useful to explore thelimitations of the LPT series independent of the additional complications that inhomogeneity entails. During the time that the spherical perturbation evolves as an independent homogeneous universe it may be fully describedin terms of the motion of its outer edge r p . Write r p ( t ) = b ( t ) X (5)where b ( t ) is the scale factor and X is the Lagrangian coordinate of the edge. The initial matter density of the homogeneoussphere ρ ( X, t ) = ρ p = ρ (1 + δ ). The physical density of the perturbation at time t is ρ ( X, t ) = ρ ( X, t ) J ( X, t ) J ( X, t ) (6)where the Jacobian of the transformation relating the Lagrangian and physical spaces is J ( X, t ) = det (cid:18) ∂(cid:126)r∂ (cid:126)X (cid:19) . (7)Since eq. (5) implies J ( X, t ) = b ( t ) and the choice a = b implies J ( X, t ) = a the perturbation matter density at latertimes is ρ p ( t ) = ρ (1 + δ ) a b ( t ) . (8)Substituting for ρ p and r p in eq. (1) gives¨ bb = − H a (1 + δ ) b (9)with initial conditions b ( t ) = a and ˙ b ( t ) = ˙ a (1 + δ v ). The curl of the acceleration (i.e. eq. (2)) vanishes by sphericalsymmetry. The corresponding equation for the background scale factor is¨ aa = − H a a (10)with initial conditions a ( t ) = a and ˙ a ( t ) = ˙ a = a H . The solution for b ( t ) will be expressed in terms of its deviationsfrom a ( t ).In summary, the physical setup is an Ω = 1 background model and a compensated spherical top-hat (over- or underdense).The properties of interest are the relative scale factors a ( t ) /a and b ( t ) /a (the choice of a is arbitrary and b = a ). Theevolution of the relative scale factors is fully specified by H , H p and Ω p at time t . The perturbed physical quantities, H p and Ω p , may be equivalently specified by a choice of δ and δ v . Appendix A contains a systematic description and enumeratesdegrees of freedom, parameters, constraints, etc. c (cid:13) RAS, MNRAS000
Sharvari Nadkarni-Ghosh and David F. Chernoff and the sphere respectively. Two sets of Lagrangian coordinates Y = r b, /a and X = r p, /b are defined. A gauge choice sets a = b . Appendix A provides a figure and gives a somewhat more detailed chain of reasoning that clarifies the constructionof the physical and Lagrangian coordinate systems. The initial perturbation is characterised by the independent parameters δ = ρ p ρ − δ v = H p H − . (3)Finally, assume that the background cosmology is critical Ω = 1. The perturbed sphere hasΩ p = 1 + δ (1 + δ v ) . (4)The physical problem of interest here is the future evolution of an arbitrary initial state unconstrained by the past history. Ingeneral, the background and the perturbation can have different big bang times. Initial conditions with equal big bang timeswill be analysed as a special case of interest and imply an additional relationship between δ and δ v .While the previous paragraphs summarise the set up, they eschew the complications in modelling an inhomogeneoussystem in terms of separate inner and outer homogeneous universes. For example, matter motions within the perturbed innerregion may overtake the outer homogeneous region so that there are problem-specific limits on how long solutions for thescale factors a ( t ) and b ( t ) remain valid. The appendix shows that there exist inhomogeneous initial configurations for whichthe limitations arising from the convergence of the LPT series are completely independent of the limitations associated withcollisions or crossings of inner and outer matter-filled regions. A basic premise of this paper is that it is useful to explore thelimitations of the LPT series independent of the additional complications that inhomogeneity entails. During the time that the spherical perturbation evolves as an independent homogeneous universe it may be fully describedin terms of the motion of its outer edge r p . Write r p ( t ) = b ( t ) X (5)where b ( t ) is the scale factor and X is the Lagrangian coordinate of the edge. The initial matter density of the homogeneoussphere ρ ( X, t ) = ρ p = ρ (1 + δ ). The physical density of the perturbation at time t is ρ ( X, t ) = ρ ( X, t ) J ( X, t ) J ( X, t ) (6)where the Jacobian of the transformation relating the Lagrangian and physical spaces is J ( X, t ) = det (cid:18) ∂(cid:126)r∂ (cid:126)X (cid:19) . (7)Since eq. (5) implies J ( X, t ) = b ( t ) and the choice a = b implies J ( X, t ) = a the perturbation matter density at latertimes is ρ p ( t ) = ρ (1 + δ ) a b ( t ) . (8)Substituting for ρ p and r p in eq. (1) gives¨ bb = − H a (1 + δ ) b (9)with initial conditions b ( t ) = a and ˙ b ( t ) = ˙ a (1 + δ v ). The curl of the acceleration (i.e. eq. (2)) vanishes by sphericalsymmetry. The corresponding equation for the background scale factor is¨ aa = − H a a (10)with initial conditions a ( t ) = a and ˙ a ( t ) = ˙ a = a H . The solution for b ( t ) will be expressed in terms of its deviationsfrom a ( t ).In summary, the physical setup is an Ω = 1 background model and a compensated spherical top-hat (over- or underdense).The properties of interest are the relative scale factors a ( t ) /a and b ( t ) /a (the choice of a is arbitrary and b = a ). Theevolution of the relative scale factors is fully specified by H , H p and Ω p at time t . The perturbed physical quantities, H p and Ω p , may be equivalently specified by a choice of δ and δ v . Appendix A contains a systematic description and enumeratesdegrees of freedom, parameters, constraints, etc. c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation (cid:72) ∆ , ∆ v (cid:76) (cid:205) (cid:72) (cid:68) , Θ (cid:76) (cid:68) Θ c (cid:45) Θ c (cid:43) initially expandinginitially contracting (cid:45) (cid:45) (cid:45) ∆ ∆ v (cid:45) (cid:45) (cid:45) ∆ ∆ v Figure 2.
Phase diagram of density and velocity perturbations ( δ, δ v ). Physical initial conditions require − < δ < ∞ and −∞ < δ v < ∞ .The left panel highlights the qualitatively different initial conditions. The shaded (unshaded) region corresponds to closed (open) modelwith negative (positive) total energy. For small ∆, models with θ − c < θ < θ + c are closed. Initially expanding and contracting models areseparated by the dashed horizontal line ( δ v = − δ and δ v . The solid blue line corresponds tothe “Zeldovich” condition i.e no perturbation at t = 0. The points ( − , − , ( − , .
5) and (0 ,
0) are unstable, stable and saddle fixedpoints of the phase space flow. The flow lines (indicated by the blue vectors) converge along the Zeldovich curve either to the stable fixedpoint at ( − , .
5) or move parallel to the Zeldovich curve to a future density singularity. Further discussion follows in § § The initial density and velocity perturbations are taken to be of the same order in the formalism developed by Buchert(1992, 1994), Buchert & Ehlers (1993) and Ehlers & Buchert (1997). We assume the same ordering here. Write the initialperturbation ( δ , δ v ) in terms of magnitude ∆ and angle θ ∆ = (cid:112) δ + δ v (11)so that δ = ∆ cos θ (12) δ v = ∆ sin θ. (13)To map physical perturbations ( δ, δ v ) in a unique manner to (∆ , θ ) adopt the ranges ∆ (cid:62) − π < θ (cid:54) π . Figure 2 (leftpanel) shows the phase space of initial perturbations. Since density is non-negative the regime of physical interest is δ (cid:62) − δ v >
0, lie above the horizontal dashed line. The right panel of figure 2 summarises the overall evolution of thesystem. The initial choice of δ and δ v dictates the trajectory in the plane. Cosmologically relevant initial conditions generallyassume there to be no perturbation at t = 0. We adopt the name “Zeldovich” initial conditions for models that satisfy thiscondition. This establishes a specific relation between δ and δ v which is indicated by the sold blue line. The exact mathematicalrelationship is given in § δ, δ v ), the system as it evolves traces out a curve in phasespace indicated by the blue arrows. There are three fixed points visible. The origin ( δ, δ v ) ≡ (0 , − , −
1) is a unstable node and thevacuum, expanding model at ( − , .
5) is a degenerate attracting node. Far to the right and below the dashed line the modelscollapse to a future singularity. The phase portrait illustrates that the trajectories either converge to the vacuum, expandingmodel or to the singular, collapsing model. The equations that govern the flow and further relevance of the Zeldovich solutionis discussed in § § c (cid:13) RAS, MNRAS , 1–38
Sharvari Nadkarni-Ghosh and David F. Chernoff
The scale factor is formally expanded b ( t ) = ∞ (cid:88) n =0 b ( n ) ( t )∆ n (14)where b ( n ) denotes an n -th order term. The initial conditions are b ( t ) = a ( t ) (15)˙ b ( t ) = ˙ a (1 + δ v ) = ˙ a (1 + ∆ sin θ ) . (16)Substitute the expansion for b ( t ) into eq. (9), equate orders of ∆ to give at zeroth order ¨ b (0) + 12 H a b (0)2 = 0 (17)which is identical in form to eq. (10) for the unperturbed background scale factor. The initial conditions at zeroth order: b (0) ( t ) = a (18)˙ b (0) ( t ) = ˙ a . (19)The equation and initial conditions for b (0) ( t ) simply reproduce the background scale factor evolution b (0) ( t ) = a ( t ). Withoutloss of generality assume that the background model has big bang time t = 0 so that a ( t ) = a (cid:16) tt (cid:17) / = a (cid:16) H t (cid:17) / . (20)At first order ¨ b (1) − H a b (1) a = − H a cos θa (21)and, in general,¨ b ( n ) − H a b ( n ) a = S ( n ) (22)where S ( n ) depends upon lower order approximations ( b (0) , b (1) . . . b ( n − ) as well as θ . The first few are: S (2) = − H a a (cid:2) b (1) (cid:8) b (1) − a cos θ (cid:9)(cid:3) (23) S (3) = − H a a (cid:104) b (1) (cid:110) − (cid:0) b (1) (cid:1) + 6 ab (2) + 3 ab (1) cos θ (cid:111) − a b (2) cos θ (cid:105) (24) S (4) = − H a a (cid:104)(cid:0) b (1) (cid:1) (cid:110) (cid:0) b (1) (cid:1) − ab (2) − ab (1) cos θ (cid:111) +6 a b (1) (cid:8) b (3) + b (2) cos θ (cid:9) + 3 a (cid:0) b (2) (cid:1) − a b (3) cos θ (cid:105) . (25)These terms can be easily generated by symbolic manipulation software. The initial conditions are b (1) ( t ) = 0 (26)˙ b (1) ( t ) = ˙ a sin θ (27)and for n > b ( n ) ( t ) = 0 (28)˙ b ( n ) ( t ) = 0 . (29)The ordinary differential equations for b ( n ) may be solved order-by-order.To summarise, the structure of the hierarchy and the simplicity of the initial conditions allows the evaluation of thesolution at any given order in terms of the solutions with lower order. This yields a formal expansion for the scale factor ofthe sphere b = ∞ (cid:88) n =0 b ( n ) ( t )∆ n (30) c (cid:13) RAS, MNRAS000
The scale factor is formally expanded b ( t ) = ∞ (cid:88) n =0 b ( n ) ( t )∆ n (14)where b ( n ) denotes an n -th order term. The initial conditions are b ( t ) = a ( t ) (15)˙ b ( t ) = ˙ a (1 + δ v ) = ˙ a (1 + ∆ sin θ ) . (16)Substitute the expansion for b ( t ) into eq. (9), equate orders of ∆ to give at zeroth order ¨ b (0) + 12 H a b (0)2 = 0 (17)which is identical in form to eq. (10) for the unperturbed background scale factor. The initial conditions at zeroth order: b (0) ( t ) = a (18)˙ b (0) ( t ) = ˙ a . (19)The equation and initial conditions for b (0) ( t ) simply reproduce the background scale factor evolution b (0) ( t ) = a ( t ). Withoutloss of generality assume that the background model has big bang time t = 0 so that a ( t ) = a (cid:16) tt (cid:17) / = a (cid:16) H t (cid:17) / . (20)At first order ¨ b (1) − H a b (1) a = − H a cos θa (21)and, in general,¨ b ( n ) − H a b ( n ) a = S ( n ) (22)where S ( n ) depends upon lower order approximations ( b (0) , b (1) . . . b ( n − ) as well as θ . The first few are: S (2) = − H a a (cid:2) b (1) (cid:8) b (1) − a cos θ (cid:9)(cid:3) (23) S (3) = − H a a (cid:104) b (1) (cid:110) − (cid:0) b (1) (cid:1) + 6 ab (2) + 3 ab (1) cos θ (cid:111) − a b (2) cos θ (cid:105) (24) S (4) = − H a a (cid:104)(cid:0) b (1) (cid:1) (cid:110) (cid:0) b (1) (cid:1) − ab (2) − ab (1) cos θ (cid:111) +6 a b (1) (cid:8) b (3) + b (2) cos θ (cid:9) + 3 a (cid:0) b (2) (cid:1) − a b (3) cos θ (cid:105) . (25)These terms can be easily generated by symbolic manipulation software. The initial conditions are b (1) ( t ) = 0 (26)˙ b (1) ( t ) = ˙ a sin θ (27)and for n > b ( n ) ( t ) = 0 (28)˙ b ( n ) ( t ) = 0 . (29)The ordinary differential equations for b ( n ) may be solved order-by-order.To summarise, the structure of the hierarchy and the simplicity of the initial conditions allows the evaluation of thesolution at any given order in terms of the solutions with lower order. This yields a formal expansion for the scale factor ofthe sphere b = ∞ (cid:88) n =0 b ( n ) ( t )∆ n (30) c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation which encapsulates the Lagrangian perturbation treatment. The right hand size explicitly depends upon the size of theperturbation and time and implicitly upon a , H , and θ . This hierarchy of equations is identical to that generated by thefull formalism developed by Buchert and collaborators when it is applied to the top-hat problem. The convergence propertiesin time and in ∆ are distinct; a simple illustrative example of this phenomenon is presented in Appendix B. The series solution outlined in the previous section does not converge at all times. Figure 1 is a practical demonstration ofthis non-convergence for the case of an expanding void. An understanding of the convergence of the LPT series is achievedby extending the domain of the expansion variable ∆ from the real positive axis to the complex plane.
The differential eq. (9) and initial conditions for the physical system are¨ b ( t ) = − H a (1 + ∆ cos θ ) b ( t ) b ( t ) = a ˙ b ( t ) = ˙ a (1 + ∆ sin θ ) (31)where t , b ( t ), ∆ and all zero-subscripted quantities are real. This set may be extended by allowing ∆ and b to become complexquantities, denoted hereafter, ∆ and b , while the rest of the variables remain real. The complex set is¨ b ( t ) = − H a (1 + ∆ cos θ ) b ( t ) b ( t ) = a ˙ b ( t ) = ˙ a (1 + ∆ sin θ ) . (32)The theory of differential equations (for example, Chicone 2006) guarantees that the solution to a real initial value problemis unique and smooth in the initial conditions and parameters of the equation and can be extended in time as long as there areno singularities in the differential equation (hereafter, the maximum extension of the solution). First, note that each complexquantity in eq. (32) may be represented by a real pair, i.e. b = u + iv by pair { u, v } = {(cid:60) b , (cid:61) b } and ∆ = x + iy by pair { x, y } = {(cid:60) ∆ , (cid:61) ∆ } . The basic theory implies continuity and smoothness of solution u and v with respect to initial conditionsand parameters x and y . Second, observe that the Cauchy-Riemann conditions u x = v y and u y = − v x are preserved by theform of the ordinary differential equation. Since the initial conditions and parameter dependence are holomorphic functionsof ∆ it follows that b ( t, ∆ ) is a holomorphic function of ∆ at times t within the maximum extension of the solution.Inspection shows that the differential equation is singular only at b = 0. For a particular value of ∆ = ∆ (cid:48) , the solutionto the initial value problem can be extended to a maximum time t mx such that b ( ∆ (cid:48) , t mx ) = 0 or to infinity. The existence ofa finite t mx signals that a pole in the complex analytic function b ( ∆ , t ) forms at ∆ = ∆ (cid:48) and t = t mx . For times t such that t (cid:54) t < t mx , the solution b ( ∆ , t ) is analytic in a small neighbourhood around the point ∆ (cid:48) . Of course, there may be poleselsewhere in the complex ∆ plane.The relationship between the original, real-valued physical problem and the complexified system is the following. In theoriginal problem ∆ is a real, positive quantity at t . LPT is a power series expansion in ∆ about the origin (the point ∆ = 0).LPT’s convergence at any time t can be understood by study of the complexified system. Consider the complex disk D centredon the origin and defined by | ∆ | < ∆. At t each point in D determines a trajectory b ( ∆ , t ) for the complexified systemextending to infinity or limited to finite time t = t mx ( ∆ ) because of the occurrence of a pole. The time of validity is definedas T (∆) = min D t mx , i.e. the minimum t mx over the disk. Since there are no poles in D at t the time of validity is the spanof time when D remains clear of any singularities. If a function of a complex variable is analytic throughout an open diskcentred around a given point in the complex plane then the series expansion of the function around that point is convergent(Brown & Churchill 1996). The LPT expansion for the original problem converges for times less than the time of validitybecause the complex extension b ( ∆ , t ) is analytic throughout D for t < T (∆). If ∆ < ∆ then, in an obvious notation, thedisks are nested D (∆ ) ⊂ D (∆ ) and the times of validity are ordered T (∆ ) (cid:62) T (∆ ).This idea is shown in figure 3. No singularities are present for the initial conditions at t ; at t a singularity is present outside the disk but it does not prevent the convergence of the LPT expansion with ∆ equal to the disk radius shown; at t a singularity is present in the disk or on its boundary and it may interfere with convergence.A distinct but related concept is the maximum amplitude perturbation for which the LPT expansion converges at theinitial time and at all intermediate times up to a given time. The radius of convergence R ∆ ( t ) is the maximum disk radius ∆for which t > T (∆). Because the disks are nested if t < t then R ∆ ( t ) (cid:62) R ∆ ( t ). c (cid:13) RAS, MNRAS , 1–38
Sharvari Nadkarni-Ghosh and David F. Chernoff
Im[ Δ ] time t = t Re[ Δ ] t = t Δ =0 t = t Figure 3.
This figure is a schematic illustration of how the time of validity is determined. The initial conditions imply a specific, real∆ at time t . The LPT series is an expansion about ∆ = 0, convergent until a pole appears at some later time within the disk of radius∆ (shown in cyan) in the complex ∆ plane. Typically, the pole’s position forms a curve (blue dashed) in the three dimensional space( (cid:60) [ ∆ ] , (cid:61) [ ∆ ] , t ). The black dots mark the pole at times t and t . At t the pole does not interfere with the convergence of the LPT series;at t it does. The time of validity may be determined by a pole that appears within the disk without moving through the boundary (notillustrated). The time of validity and the radius of convergence are inverse functions of each other. If the initial perturbation isspecified, i.e. ∆ is fixed, and the question to be answered is “how far into the future does LPT work?” then the time ofvalidity gives the answer. However, if the question is “how big an initial perturbation will be properly approximated by LPTover a given time interval?” then the radius of convergence provides the answer.Finally, note that one can trivially extend this formalism to deal with time intervals in the past.
The following recipe shows how to calculate the radius of convergence R ∆ ( t ) and the time of validity T (∆) efficiently. Fix a , H , t and θ ; these are all real constants set by the initial conditions. Assume that it is possible to find b ( ∆ , t ) for complex ∆ and real t by solving eq. (32). There exist explicit expressions for b as will be shown later.Start with t = t and R ∆ ( t ) = ∞ . The iteration below maps out R ∆ ( t ) by making small increments in time δt . • Store old time t previous = t , choose increment δt and form new time of interest t = t previous + δt . • Locate all the ∆ which solve b ( ∆ , t ) = 0. The roots correspond to poles in the complex function. Find the root closestto the origin and denote its distance as | ∆ near | . • The radius of convergence is R ∆ ( t ) = min( | ∆ near | , R ∆ ( t previous )). • Continue.Since R ∆ is decreasing, the inversion to form T (∆) is straightforward. Figure 4 shows a schematic cartoon of the con-struction process. c (cid:13) RAS, MNRAS000
The following recipe shows how to calculate the radius of convergence R ∆ ( t ) and the time of validity T (∆) efficiently. Fix a , H , t and θ ; these are all real constants set by the initial conditions. Assume that it is possible to find b ( ∆ , t ) for complex ∆ and real t by solving eq. (32). There exist explicit expressions for b as will be shown later.Start with t = t and R ∆ ( t ) = ∞ . The iteration below maps out R ∆ ( t ) by making small increments in time δt . • Store old time t previous = t , choose increment δt and form new time of interest t = t previous + δt . • Locate all the ∆ which solve b ( ∆ , t ) = 0. The roots correspond to poles in the complex function. Find the root closestto the origin and denote its distance as | ∆ near | . • The radius of convergence is R ∆ ( t ) = min( | ∆ near | , R ∆ ( t previous )). • Continue.Since R ∆ is decreasing, the inversion to form T (∆) is straightforward. Figure 4 shows a schematic cartoon of the con-struction process. c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation Complex Δ plane Root Plot Im[ Δ ] Re[ Δ ]Roots at t Roots at t Δ T( Δ ) t R Δ (t) Figure 4.
A schematic illustration of the radius of convergence and the time of validity. The left panel shows the location of poles inthe complex ∆ plane at times t and t , denoted by orange squares and green dots, respectively. At a fixed time, the pole nearest theorigin determines the disk (black circle) within which a series expansion about the origin converges. The right panel shows | ∆ | for t and t . The black line is R ∆ ( t ), the minimum | ∆ | calculated for a continuous range of times (where t , the initial time, lies far to theleft). The arrows show how the time of validity is inferred for a given ∆. The usual parametric representation provides an efficient method to construct an explicit complex representation for b ( ∆ , t ). The original system eq. (31) depends upon a , H , θ and ∆. The assumed Einstein-deSitter background has a > a >
0; as defined, the perturbation amplitude ∆ (cid:62) θ with − π < θ (cid:54) π . The quantity (1 + ∆ cos θ ) is proportional to total density and must be non-negative. Thesign of ˙ b is the sign of 1 + ∆ sin θ and encodes expanding and contracting initial conditions.Briefly reviewing the usual physical solution, the integrated form is˙ b = H a (cid:20) (1 + ∆ cos θ ) b + (1 + ∆ sin θ ) − (1 + ∆ cos θ ) a (cid:21) . (33)The combination E (∆ , θ ) = (1 + ∆ sin θ ) − (1 + ∆ cos θ ) (34)is proportional to the total energy of the system. If E >
E < E = 0 which separates open and closed regions. For infinitesimal ∆ the line of divisionhas slope tan θ = 1 /
2. Models with θ ∈ [ θ − c , θ + c ] = [ − π + tan − (1 / , tan − (1 / − . , .
46] are closed while those outsidethis range are open.There are four types of initial conditions (positive and negative E , positive and negative ˙ b ) and four types of solutions,shown schematically in figure 5. The solutions have well-known parametric forms involving trigonometric functions of angle η or iη (see Appendix C). The convention adopted here is that the singularity nearest the initial time t coincides with η = 0 andis denoted t + bang ( t − bang ) for initially expanding (contracting) solutions (see figure 5). The time interval between the singularityand t is t age = | t − t ± bang | (cid:62) c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff open (cid:72) E (cid:62) (cid:76) closed (cid:72) E (cid:60) (cid:76) dbdt (cid:72) t (cid:61) t (cid:76) (cid:62) t (cid:43) bang t (cid:43) bang t t coll b (cid:72) t (cid:76) open (cid:72) E (cid:62) (cid:76) closed (cid:72) E (cid:60) (cid:76) dbdt (cid:72) t (cid:61) t (cid:76) (cid:60) t coll t t (cid:45) bang t (cid:45) bang b (cid:72) t (cid:76) Figure 5.
Scale factor as a function of time. The initial conditions ( b = a = 1 and varying ˙ b ) are given at time t (dashed blue line).The left (right) panel illustrates initially expanding (contracting) models. t ± bang corresponds to η = 0; t coll to η = 2 π . For expandingsolutions t age = t − t + bang is the time interval since the initial singularity and t coll is the future singularity for closed models. Forcontracting solutions t age = t − bang − t is the time until the final singularity and t coll is the past singularity for closed models. The parametric solution for the models can be written as b ( η, ∆ , θ ) = a θ )[ − E (∆ , θ )] (1 − cos η ) t ( η, ∆ , θ ) = t ± (cid:18) H (1 + ∆ cos θ )[ − E (∆ , θ )] / ( η − sin η ) − t age (∆ , θ ) (cid:19) . (35)The plus and minus signs give the solution for initially expanding and initially contracting models respectively. Parameter η is purely real for closed solutions and purely imaginary for open solutions. The distance to the nearest singularity is t age = (cid:90) b = a b =0 db [˙ b ] / = 1 H (cid:90) y =1 y =0 dy [(1 + ∆ cos θ ) y − + E (∆ , θ )] / . (36)The second equality uses eq. (33) and the substitution y = b/a . To extend the above parametric solution to the complex plane, one might guess the substitution ∆ → ∆ e iφ where − π < φ (cid:54) π in eq. (35) and eq. (36). The physical limit is φ = 0. However, this leads to two problems. First, the integral for t age can havemultiple extensions that agree for physical φ = 0 but differ elsewhere including the negative real axis. This is tied to the factthat the operations of integration and substitution ∆ → ∆ e iφ do not commute because of the presence of the square root inthe expression for t age . A second related problem is the presence of multiple square roots in the parametric form for t . Thesegive rise to discontinuities along branch cuts such that one parametric form need not be valid for the entire range of φ , butinstead the solution may switch between different forms. Directly extending the parametric solution is cumbersome.However, the original differential eq. (32) is manifestly single-valued. The equation can be integrated forward or backwardnumerically to obtain the correct solution for complex ∆ . One can then match the numerical solution to the above parametricforms to select the correct branch cuts. This procedure was implemented to obtain the form for all ∆ and θ . The main resultis that the solution space for all θ and ∆ is completely spanned by complex extensions of the two real parametric forms whichdescribe initially expanding and contracting solutions. The expressions for t age and details are given in Appendix C3.The traditional textbook treatment relating physical cosmological models with real Ω > < η → iη in the parametric forms and one verifies that this exchanges closed and open solutions.However, starting from the second order differential equation it is straightforward to use the same type of reasoning as aboveto construct an explicit analytic continuation from one physical regime to the other. c (cid:13) RAS, MNRAS000
Scale factor as a function of time. The initial conditions ( b = a = 1 and varying ˙ b ) are given at time t (dashed blue line).The left (right) panel illustrates initially expanding (contracting) models. t ± bang corresponds to η = 0; t coll to η = 2 π . For expandingsolutions t age = t − t + bang is the time interval since the initial singularity and t coll is the future singularity for closed models. Forcontracting solutions t age = t − bang − t is the time until the final singularity and t coll is the past singularity for closed models. The parametric solution for the models can be written as b ( η, ∆ , θ ) = a θ )[ − E (∆ , θ )] (1 − cos η ) t ( η, ∆ , θ ) = t ± (cid:18) H (1 + ∆ cos θ )[ − E (∆ , θ )] / ( η − sin η ) − t age (∆ , θ ) (cid:19) . (35)The plus and minus signs give the solution for initially expanding and initially contracting models respectively. Parameter η is purely real for closed solutions and purely imaginary for open solutions. The distance to the nearest singularity is t age = (cid:90) b = a b =0 db [˙ b ] / = 1 H (cid:90) y =1 y =0 dy [(1 + ∆ cos θ ) y − + E (∆ , θ )] / . (36)The second equality uses eq. (33) and the substitution y = b/a . To extend the above parametric solution to the complex plane, one might guess the substitution ∆ → ∆ e iφ where − π < φ (cid:54) π in eq. (35) and eq. (36). The physical limit is φ = 0. However, this leads to two problems. First, the integral for t age can havemultiple extensions that agree for physical φ = 0 but differ elsewhere including the negative real axis. This is tied to the factthat the operations of integration and substitution ∆ → ∆ e iφ do not commute because of the presence of the square root inthe expression for t age . A second related problem is the presence of multiple square roots in the parametric form for t . Thesegive rise to discontinuities along branch cuts such that one parametric form need not be valid for the entire range of φ , butinstead the solution may switch between different forms. Directly extending the parametric solution is cumbersome.However, the original differential eq. (32) is manifestly single-valued. The equation can be integrated forward or backwardnumerically to obtain the correct solution for complex ∆ . One can then match the numerical solution to the above parametricforms to select the correct branch cuts. This procedure was implemented to obtain the form for all ∆ and θ . The main resultis that the solution space for all θ and ∆ is completely spanned by complex extensions of the two real parametric forms whichdescribe initially expanding and contracting solutions. The expressions for t age and details are given in Appendix C3.The traditional textbook treatment relating physical cosmological models with real Ω > < η → iη in the parametric forms and one verifies that this exchanges closed and open solutions.However, starting from the second order differential equation it is straightforward to use the same type of reasoning as aboveto construct an explicit analytic continuation from one physical regime to the other. c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation In addition, note that the differential equation and its solution remain unchanged under the simultaneous transformations ∆ → − ∆ and θ → θ + π . Every complex solution with − π < θ (cid:54) < θ (cid:54) π and vice-versa. For determining the radius of convergence and the time of validity the whole disk of radius | ∆ | is searched forpoles so it suffices to consider a restricted range of θ to handle all physical initial conditions. The condition b = 0 signals the presence of a pole. Inspection of the parametric form shows that this condition can occuronly when η = 0 or η = 2 π . The corresponding time t ( ∆ , θ ) = (cid:40) t ± (cid:16) πH (1+ ∆ cos θ )[ − E ( ∆ )] / − t age ( ∆ ) (cid:17) ( η = 2 π ) t ∓ t age ( ∆ ) ( η = 0) (37)is immediately inferred. Since the independent variable t is real the transcendental equation (cid:61) t ( ∆ , θ ) = 0 (38)must be solved. It is straightforward to scan the complex ∆ plane and calculate t to locate solutions. Each solution gives aroot of b = 0 and also implies the existence of a pole at the corresponding ∆ . Note that relying upon the parametric solutionsis a far more efficient method for finding the poles than integrating the complex differential equations numerically. We haveverified that both methods produce the same results.In practice, we fix θ , scan a large area of the complex ∆ plane, locate all purely real t and save the { ∆ , t } pairs. Theseare used to create a scatter plot of | ∆ | as a function of time (hereafter the “root plot”). Generally, the location of the polesvaries smoothly with t and continuous loci of roots are readily apparent. Finding R ∆ and T (∆) follows as indicated in figure4. Root plots were calculated for a range of angles 0 (cid:54) θ (cid:54) π . Since the root plots depend upon | ∆ | they are invariant under θ → θ − π and this coverage suffices for all possible top-hat models. For the results of the full survey in θ see Appendix D.The theoretical radius of convergence R ∆ ( t ) and time of validity T (∆) follow directly.This section analyses the theoretical convergence for specific open and closed models derived from the root plots. Theseestimates are compared to the time of validity inferred by numerical evaluation of the LPT series. The range of models withlimited LPT convergence is characterised. The concept of mirror models is introduced to elucidate a number of interconnectionsbetween open and closed convergence. The physical interpretation of roots introduced by the complexification of the equationsbut lying outside the physical range are discussed. Finally, the special case where the background and the perturbation havethe same big bang time is analysed. Figure 6 shows R ∆ ( t ) for θ = 2 .
82 and initial scale factor a = 10 − . All ∆ yield expanding open models for this θ ; one choicecorresponds to the model whose LPT series appeared in figure 1 (∆ = 0 . θ = 2 . a = 10 − ). The x-axis is log a and isequivalent to a measure of time. The y-axis is log | ∆ | , i.e. the distance from the origin to poles in the complex ∆ plane. Inprinciple, future evolution may be limited by real or complex roots. The blue solid line and the red dotted line indicate realand complex roots of η = 2 π respectively. The cyan dashed and pink dot-dashed lines indicate the real and complex roots of η = 0 respectively. Future evolution is constrained by real roots (blue and cyan) in this example.The time of validity is the first instance when a singularity appears within the disk of radius ∆ in the complex ∆ plane.For the specific case, starting at ordinate ∆ = 10 − , one moves horizontally to the right to intersect the blue line and thenvertically down to read off the scale factor a v = a [ T (∆)] = 0 . δ v = 1 (or ∆ = 1 / sin θ ) at which point the root switches from η = 2 π below to 0 above.Such a switch might occur if varying the initial velocity transposes an expanding closed model into a contracting closed model.But it is expected to occur at δ v = − § c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff (cid:45) (cid:45) (cid:45) (cid:45) a v (cid:45) (cid:45) (cid:45)
101 Log (cid:72) a (cid:72) t (cid:76)(cid:76) L og (cid:200) (cid:68) (cid:200) Figure 6. R ∆ for θ = 2 .
82 and a = 10 − (vertical dashed line). To determine the time of validity for LPT expansion with a given ∆,move horizontally to the right of a = a following the dashed line with arrow and locate the first coloured line with ordinate equal to ∆and then move vertically down to read off the scale factor at the time of validity a v . The specific case illustrated (∆ = 10 − ) matchesthat of the model with problematic convergence in figure 1. The time of validity is correctly predicted. The meaning of the colours isdiscussed in the text. Coloured version of the figure is available online. openclosed (cid:45) (cid:45) (cid:45) (cid:45) a v a c a v (cid:61) a c (cid:45) (cid:45) (cid:68) rc (cid:68) E (cid:61)
01 Log (cid:72) a (cid:72) t (cid:76)(cid:76) L og (cid:200) (cid:68) (cid:200) Figure 7. R ∆ for θ = 0 .
44 and a = 10 − . The line ∆ E = 0 separates open and closed models. The scale factor at the time of validityis a v . For closed models the scale factor at time of collapse is a c . Blue solid line and red small dashed line denote real and complex rootsof η = 2 π , respectively. The cyan dashed lines denotes the real roots of η = 0. When the first singularity encountered is real, a v = a c ,the time of validity is the future time of collapse. However, when the singularity is complex the time of validity is less than the actualcollapse time. In the range ∆ rc < ∆ < ∆ E =0 , there are closed models with a v < a c . Figure 7 presents R ∆ ( t ) for models with θ = 0 .
44 and a = 10 − . There are several new features. Over the angular range θ − c < θ < θ + c the cosmology is closed for small ∆ (see shaded region in figure 2 near ∆ = 0). Conversely, a straight line drawnfrom ∆ = 0 within this angular range must eventually cross the parabola E = 0 except for the special case θ = 0. Since thevelocity contribution to energy E ∝ ∆ while the density contribution ∝ − ∆ it is clear that eventually E > E =0 , is a function of θ . Below the brown horizontal dot-dashed line in figure 7 the models are closed,above they are open (line labelled ∆ = ∆ E =0 ).The root plot has, as before, blue solid and red dotted lines denoting the distance to real and complex ∆ poles, respectively,for η = 2 π . The cyan dashed line denotes real roots for η = 0 and does not restrict future evolution.For small ∆ real roots determine the time of validity. These roots correspond exactly to the model’s collapse time. Inother words, the time of validity is determined by the future singularity. For example, for ∆ = 0 .
01, the root plot predictsthat a series expansion should be valid until the collapse at a = 5 . a v = a c ” on the x-axis . This predictionis confirmed in the left hand panel of figure 8. The root diagram is consistent with the qualitative expectation that smalloverdensities should have long times of validity because collapse times are long: lim ∆ → T (∆) → ∞ .As ∆ increases from very small values, i.e. successively larger initial density perturbations, the collapse time decreases.Eventually the velocity perturbation becomes important so that at ∆ = ∆ rc a minimum in the collapse time is reached. For∆ E =0 > ∆ > ∆ rc the collapse time increases while the model remains closed. As ∆ → ∆ E =0 the collapse time becomesinfinite and the model becomes critical. All models with ∆ > ∆ E =0 are open. c (cid:13) RAS, MNRAS000
01, the root plot predictsthat a series expansion should be valid until the collapse at a = 5 . a v = a c ” on the x-axis . This predictionis confirmed in the left hand panel of figure 8. The root diagram is consistent with the qualitative expectation that smalloverdensities should have long times of validity because collapse times are long: lim ∆ → T (∆) → ∞ .As ∆ increases from very small values, i.e. successively larger initial density perturbations, the collapse time decreases.Eventually the velocity perturbation becomes important so that at ∆ = ∆ rc a minimum in the collapse time is reached. For∆ E =0 > ∆ > ∆ rc the collapse time increases while the model remains closed. As ∆ → ∆ E =0 the collapse time becomesinfinite and the model becomes critical. All models with ∆ > ∆ E =0 are open. c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation (cid:68) (cid:61) Θ(cid:61) (cid:72) t (cid:76) b (cid:72) t (cid:76) (cid:68) (cid:61) Θ(cid:61) (cid:72) t (cid:76) b (cid:72) t (cid:76) Figure 8.
The exact solution (black, dashed) and LPT expansions of successively higher order (blue) for two expanding, closed modelswith θ = 0 .
44. The left hand panel has ∆ = 0 .
01. LPT converges to the exact solution at all times up to the singularity at a = 5 .
5. Theright hand panel has ∆ = 0 .
2. LPT does not converge beyond a = 0 . The root diagram shows that for ∆ > ∆ rc , the time of validity is determined by complex not real ∆ for η = 2 π . Closedmodels with ∆ rc < ∆ < ∆ E =0 have a time of validity less than the model collapse time. For example, for ∆ = 0 .
2, the collapseoccurs at a = 0 .
94 but convergence is limited to a (cid:54) .
38. This prediction is verified in the right panel of figure 8.The convergence of LPT expansions for some closed models is limited to times well before the future singularity. Thisgeneral behaviour is observed for θ − c < θ < θ + c and ∆ rc < ∆ < ∆ E =0 where both ∆ rc and ∆ E =0 are functions of θ . AppendixD provides additional details. The parametrization of the perturbation in terms of ∆ > − π < θ (cid:54) π and the complexification of ∆ → ∆ can give riseto poles anywhere in the complex ∆ space. When R ∆ is determined by a pole along the real positive axis, a clear interpretationis possible: the future singularity of the real physical model exerts a dominant influence on convergence. LPT expansions forclosed models with ∆ < ∆ rc are limited by the future collapse of the model and are straightforward to interpret.The meaning of real roots for open models is less clear cut. The roots determining R ∆ at large t are negative real andsmall in magnitude. Negative ∆ lies outside the parameter range for physical perturbations taken to be ∆ >
0. Nonethelessthe mapping (∆ , θ ) → ( − ∆ , θ ± π ) preserves ( δ , δ v ) and the original equations of motion. The poles of the models withparameters (∆ , θ ) and (∆ , θ ± π ) are negatives of each other. Let us call these “mirror models” of each other.For infinitesimal ∆ if the original model is open then the mirror model is closed. Figure 2 shows that the ∆ E =0 line hassome curvature (in fact, it is a parabola) whereas the mirror mapping is an exact inversion through ∆ = 0. Small ∆ pointsare mapped between open and closed; large ∆ points may connect open models to other open models.If the original model is open with limiting pole which is negative real of small magnitude then it corresponds to a futuresingularity of the closed mirror model. For example, the closed model with parameters (∆ = 0 . , θ = 0 .
44) in the left panelof figure 8 and the open model with parameters (∆ = 0 . , θ = 0 . − π ) shown in the left panel of figure 9 are mirrors. Thetime of the validity of the open model equals the time to collapse of its closed mirror.The notion of mirror models explains other features of the root diagrams. The time of validity of open models waspreviously discussed using figure 6 ( θ = 2 . θ = 2 . − π = − .
32. As ∆ increases the sequence of mirrormodels crosses the δ v = − η = 2 π to η = 0. This explains the switch in root label from blue solid to cyan dashedseen in figure 6, which occurs at δ v = 1 in the original model.The symmetry of the mirroring is not limited to cases when ∆ is real. It applies for complex ∆ , too. For example, themodels in the right panels of figures 8 and 9 are mirrors of each other. Their time of validity is the same and determined by c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff (cid:68) (cid:61)
Θ (cid:61) (cid:45)Π (cid:72) t (cid:76) b (cid:72) t (cid:76) (cid:68) (cid:61) Θ (cid:61) (cid:45)Π (cid:72) t (cid:76) b (cid:72) t (cid:76) Figure 9.
Mirror models of the closed models of figure 8. Each graph shows the exact solution (black, dashed) and the LPT expansionto successively higher orders (blue) of one mirror model. The original model and the mirror have the same time of validity for the LPTexpansion. complex roots which are negatives of each other. These singularities are non-physical and have no interpretation in terms ofthe collapse of any model yet they limit the LPT convergence in the same way.Figure 10 shows the areas of phase space where complex roots determine the time of validity in light red. The area withinthe parabola (light blue) contains closed models. Most of the light blue region has a time of validity determined by real roots,i.e. the time to the future singularity. The area with both light blue and red shading encompasses closed models with theunexpected feature that the time of validity is less than the time to collapse.The area outside the parabola contains open models. The time of validity of the unshaded region is determined byreal roots. The original observation of LPT’s non-convergence for an underdensity (Sahni & Shandarin 1996) is an examplethat falls in this region. For small amplitude perturbations the time of validity is simply related by mirror symmetry to theoccurrence of future singularities of closed models. The right hand plot in figure 9 is an example of an open model with timeof validity controlled by complex roots (red shading outside the parabola).Finally, some open models (especially those with large ∆) have mirrors that are open models. Figure 11 shows mirrormodels (∆ = 2, θ = 17 π/
36) and (∆ = 2, θ = 17 π/ − π ). These are initially expanding and contracting solutions respectively.The root plot in figure 12 predicts that the series is valid until a v = 0 . η = 0 (cyan line) sets the timeof validity and corresponds to the bang time (the future singularity) of the initially contracting model.In all cases, the analysis correctly predicts the convergence of the LPT series. The large expanse of phase space shaded light red in figure 10 suggests that complex roots should play a ubiquitous rolein LPT applications but the situation is somewhat more subtle. For good physical reasons purely gravitational cosmologicalcalculations often start with expanding, small amplitude, growing modes at a finite time after the big bang. The absenceof decaying modes implies that the linearized perturbations decrease in the past . A non-linear version of this condition isthat the perturbation amplitude is exactly zero at t = 0. The same condition can be formulated as “the background and theperturbation have the same big bang time” or “the ages of the perturbation and the background are identical.” The conditionis1 H (cid:90) y =1 y =0 dy [(1 + ∆ cos θ ) y − + E (∆ , θ )] / = 23 H . (39) Our analysis is restricted to the case of initially expanding models, i.e. near ∆ = 0. For initially contracting closed models, similarphysical arguments motivate a consideration of the behaviour near the initial singularity (not the future bang time). For initiallycontracting open models the epoch of interest is t → −∞ . These models have large ∆ and are not described by the linear limit discussedin the text. c (cid:13) RAS, MNRAS000
36) and (∆ = 2, θ = 17 π/ − π ). These are initially expanding and contracting solutions respectively.The root plot in figure 12 predicts that the series is valid until a v = 0 . η = 0 (cyan line) sets the timeof validity and corresponds to the bang time (the future singularity) of the initially contracting model.In all cases, the analysis correctly predicts the convergence of the LPT series. The large expanse of phase space shaded light red in figure 10 suggests that complex roots should play a ubiquitous rolein LPT applications but the situation is somewhat more subtle. For good physical reasons purely gravitational cosmologicalcalculations often start with expanding, small amplitude, growing modes at a finite time after the big bang. The absenceof decaying modes implies that the linearized perturbations decrease in the past . A non-linear version of this condition isthat the perturbation amplitude is exactly zero at t = 0. The same condition can be formulated as “the background and theperturbation have the same big bang time” or “the ages of the perturbation and the background are identical.” The conditionis1 H (cid:90) y =1 y =0 dy [(1 + ∆ cos θ ) y − + E (∆ , θ )] / = 23 H . (39) Our analysis is restricted to the case of initially expanding models, i.e. near ∆ = 0. For initially contracting closed models, similarphysical arguments motivate a consideration of the behaviour near the initial singularity (not the future bang time). For initiallycontracting open models the epoch of interest is t → −∞ . These models have large ∆ and are not described by the linear limit discussedin the text. c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation (cid:45) (cid:45) (cid:45) ∆ ∆ v Figure 10.
The red shaded region denotes part of phase space where complex roots play a role. The solid blue line represents the initialconditions which correspond to the background and perturbation having the same big bang time. The black solid parabola separates theclosed and open models. Coloured version online. (cid:68) (cid:61) Θ(cid:61) Π (cid:144) (cid:72) t (cid:76) b (cid:72) t (cid:76) (cid:68) (cid:61) Θ(cid:61) Π (cid:144) (cid:45)Π (cid:72) t (cid:76) b (cid:72) t (cid:76) Figure 11.
Two open models which are mirrors of each other. Each plot shows the exact solution (black, dashed) and the LPT seriesexpansion to successively higher orders (blue). The left panel is an initially expanding, open model whose convergence is limited to scalefactors less than a v = 0 . a v = 0 . This is a nonlinear relationship between the two initial parameters ∆ and θ which is shown by a thick blue line on the phasespace diagram in figure 10. We have adopted the name “Zeldovich” initial conditions for the top-hat models that satisfy theequal bang time relation. There are a variety of definitions for Zeldovich initial conditions given in the literature. Generally,these agree at linear order. This one has the virtue that it is simple and easy to interpret. Note that the blue curve does notintersect the region of phase space where complex roots occur except, possibly, near ∆ = 0.In the limit of small ∆ eq. (39) becomes∆(3 sin θ − cos θ ) = 0 . (40) c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff a v (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) a (cid:72) t (cid:76)(cid:76) L og (cid:200) (cid:68) (cid:200) Figure 12. R ∆ for θ = 17 π/
36 and a = 10 − . The blue solid and cyan dashed lines denoted real roots with η = 2 π and η = 0 respectively.For ∆ = 2, the time of validity is set by the root with η = 0, which is the bang time of the mirror model with θ = 17 π/ − π . See figure11 for the evolution of both models. The solutions are θ = θ Z ± where θ Z + = 2 .
82 and θ Z − = π − θ Z + = − .
32. The second quadrant solution θ Z + correspondsto open models while its mirror in the fourth quadrant θ Z − to closed models. Only when ∆ → − , .
5) and overdense points move towards collapse. The distribution of initial density and velocityperturbations yields a cloud of points in phase space but complex roots never play a role because nothing displaces individualpoints from the Zeldovich curve. Each moves at its own pace but stays near the curve. Alternatively, it is well known that tidalforces couple the collapse of nearby points. These interactions amplify the initial inhomogeneities leading to the formationof pancakes and filaments. As time progresses motions transverse to the Zeldovich curve will grow. If these deviations aresufficient they may push some points into areas with complex roots. In a subsequent paper, we will explore these issues forgeneral inhomogenous initial conditions.
To overcome the constraints above, an iterative stepping scheme that respects the time of validity is developed for LPT. Theinitial parameters at the first step determine the solution for some finite step size. The output at the end of the first stepdetermines the input parameter values for the next step and so on.
Choose the background ( a , H , Ω = 1, Y ) and the perturbation ( b = a , H p , Ω p , X ) at initial time t . The perturbedmodel is fully characterised by H p and Ω p or by δ = ρ p /ρ − δ v, = H p /H − and θ . Extra subscriptshave been added to label steps.LPT converges for times t < T (∆ , θ ). Use LPT to move forward to time t ∗ satisfying t < t ∗ < T (∆ , θ ). At t ∗ ,the background and perturbed scale factors and time derivatives are a ∗ , b ∗ , ˙ a ∗ , and ˙ b ∗ . The fractional density and velocityperturbations with respect to the background are δ ∗ = (1 + δ ) (cid:16) a ∗ b ∗ (cid:17) − c (cid:13) RAS, MNRAS000
Choose the background ( a , H , Ω = 1, Y ) and the perturbation ( b = a , H p , Ω p , X ) at initial time t . The perturbedmodel is fully characterised by H p and Ω p or by δ = ρ p /ρ − δ v, = H p /H − and θ . Extra subscriptshave been added to label steps.LPT converges for times t < T (∆ , θ ). Use LPT to move forward to time t ∗ satisfying t < t ∗ < T (∆ , θ ). At t ∗ ,the background and perturbed scale factors and time derivatives are a ∗ , b ∗ , ˙ a ∗ , and ˙ b ∗ . The fractional density and velocityperturbations with respect to the background are δ ∗ = (1 + δ ) (cid:16) a ∗ b ∗ (cid:17) − c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation δ v, ∗ = ˙ b ∗ /b ∗ ˙ a ∗ /a ∗ − . (42)Re-expand the perturbation around the background model as follows. First, let the time and Lagrangian coordinate for thebackground (inner edge of the unperturbed sphere) be continuous: t = t ∗ and Y = Y . These imply a = a ∗ and ˙ a = ˙ a ∗ ,i.e. the scale factor and Hubble constant for the background are continuous.At the beginning of the first step we assumed a = b . This is no longer true at the end of the first step. Define a newLagrangian coordinate X = X b ∗ /a ∗ , new scale factor b = a ∗ , and new scale factor derivative ˙ b = ˙ b ∗ a ∗ /b ∗ . These definitionsleave the physical edge of the sphere and its velocity unaltered r physical, ∗ = b ∗ X = b X (43)˙ r physical, ∗ = ˙ b ∗ X = ˙ b X . (44)The re-definitions relabel the fluid elements with a new set of Lagrangian coordinates and re-scale the scale factor. Theperturbation parameters are unchanged δ = δ ∗ and δ v, = δ v, ∗ because physical quantities are unmodified. Consequently,∆ = ∆ ∗ and θ = θ ∗ . To examine how Lagrangian re-expansion works consider how the Lagrangian parameters ∆ and θ would vary if they wereevaluated at successive times over the course of a specific cosmological history. Let δ ( t ) and δ v ( t ) be defined via eq. (3) andapply the second-order equations of motion eqs. (9) and (10) to derive the coupled first-order system dδdt = − t δ v (1 + δ ) (45) dδ v dt = 13 t { (1 + δ v )(1 − δ v ) − (1 + δ ) } (46)where all occurrences of δ and δ v are functions of time. From δ ( t ) and δ v ( t ) one infers the parameters, ∆( t ) and θ ( t ). Thesehave the following simple interpretation: a Lagrangian treatment starting at time t (cid:48) has ∆ = ∆( t (cid:48) ) and θ = θ ( t (cid:48) ) in the LPTseries.Since the system is autonomous it reduces to a simple flow in phase space. The flow has three fixed points at ( δ, δ v ) = (0 , − , − − , . ,
0) shows it is a saddle fixed point. The tangent to the E = 0 curve at the origin is the attracting direction andthe tangent to the equal big bang curve is the repelling direction. The fixed point at ( − , .
5) is a degenerate attracting nodeand that at ( − , −
1) is an unstable node. The flow vectors are plotted in the left panel of figure 13. The blue shaded regionindicates closed models and red shaded region indicates models where complex roots limit the time of validity for LPT.Note that the flow lines smoothly cover the whole phase space. The interpretation is that the continuous relabelling ofLagrangian coordinates and re-scaling of the scale factor has the potential to overcome the convergence limitations discussedthus far. Otherwise one might have seen ill-defined or incomplete flows or flows that were confined to a given region.
The right panel of figure 13 zooms in on the area near the origin. Initial points that correspond to open models starting nearthe origin approach the Zeldovich curve and asymptotically converge to the strong attractor at ( δ, δ v ) = ( − , . δ → ∞ . In the asymptotic limit, the solution to (46) is given by δ ∼ δ v + K withintegration constant K . From figure 13, the flow lines of closed models that start in the vicinity of the origin trace a parabolicpath that is parallel and essentially equivalent to the Zeldovich curve.The flow shows where re-expansion is needed. Closed model flow lines that start near the origin never pass through the redshaded region where complex roots play a role; the time of validity equals the time to collapse and no re-expansion is needed.However, closed models that originate in the red region must be re-expanded. The flow suggests that they eventually moveinto the blue region. So even though a closed model may initially have an LPT series with limited convergence, re-expansionmakes it possible to move into the part of phase space where a single step suffices to reach collapse. This section and the next examine the feasibility of extending a solution from recombination to today. The results will beapplied to fully inhomogeneous evolution in future paper.Let the asymptotic time of validity for an open model be expressed in dimensionless form χ = lim t →∞ H ( t ) T (∆( t ) , θ ( t )).Here, ∆ → (cid:112) / θ → tan − ( − /
2) = 2 .
677 and T (∆ , θ ) is determined by the future time to collapse of the closed mirror c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff (cid:45) (cid:45) (cid:45) (cid:45) ∆ ∆ v (cid:45) (cid:45) (cid:45) (cid:45) ∆ ∆ v Figure 13.
The left panel shows streamlines of the flow described by eq. (46). The colour coding of the plot is same as figure 10. Theright panel zooms in on the area near the origin which is where all models are located at sufficiently early times. At late times, openmodels move away from the origin towards the attracting fixed point at ( δ, δ v ) = ( − , . (cid:68) (cid:61) (cid:68) (cid:61) Θ (cid:61) Θ (cid:61) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45)
101 Log (cid:72) a (cid:72) t (cid:76)(cid:76) L og (cid:200) (cid:68) (cid:200) Figure 14.
Extending the time of validity of LPT. The first step has ∆ = 10 − and θ = 2 .
82 and implies scale factor at the time ofvalidity a v = 0 . , θ ) = (0 . , . model. The result is χ = 2 .
62 (numerical results in Appendix D), the time of validity is proportional to the characteristic ageof the background and individual steps grow larger and larger.An example shows that the basic effect can be seen even before the asymptotic regime is achieved. Figure 14 sketchesthe first two steps where the assumed model parameters at the first step are (∆ , θ ) = (0 . , . a = 0 . , θ ) = (0 . , .
68) or ( δ , δ v, ) = ( − . , . T ) andafter ( T (cid:48) ) a step α = T (cid:48) T . (47)Figure 15 shows α evaluated along the continuous flow as a function of scale factor for three different starting initial conditions.Since α > t i the time after N steps is roughly t ∼ α N t i > N t i .Consider, for example, the number of steps needed to extend an open solution from recombination to today. Let t f ( t i ) c (cid:13) RAS, MNRAS000
68) or ( δ , δ v, ) = ( − . , . T ) andafter ( T (cid:48) ) a step α = T (cid:48) T . (47)Figure 15 shows α evaluated along the continuous flow as a function of scale factor for three different starting initial conditions.Since α > t i the time after N steps is roughly t ∼ α N t i > N t i .Consider, for example, the number of steps needed to extend an open solution from recombination to today. Let t f ( t i ) c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation (cid:72) t (cid:76) Α Figure 15.
The ratio of successive times of validity ( α ) vs. a ( t ). The dashed, dot-dashed and dotted lines indicate three initial startingpoints (0 . , . , (0 , , ( − . , .
2) respectively. The ratio converges to about 3.6 and the time of validity increases geometrically with N . be the final (initial) time of interest where t f /t i ∼ a f /a i ∼ . . Estimating α = 3 implies N ∼ log . ∼
10 steps areneeded. This numerical result for N is an overestimate and one can do better. It is important to recall that it based on anarbitrarily high order expansion which achieves an exact solution. If one is limited to calculations of finite Lagrangian orderand imposes a maximum numerical error at the end of the calculation then more than N steps may be required. At least N steps are needed for series convergence and more than N steps may be needed for error control.One can extend any open model to an arbitrary future time while respecting the time of validity of the LPT series. Thenumber of steps is governed by a geometric progression.One can also extend any closed model to the future singularity while respecting the time of validity of the LPT series.Only a single step is needed for a closed model when the root is real (blue shaded region of figure 13). When it is complex(the region shaded both blue and red) the model flows first toward the node at (0 ,
0) (∆ decreases) and ultimately reaches theregion of real roots. Multiple steps will generally be necessary to escape the region of complex roots. An approximate fit (eq.(D2)) shows that χ ∼ T (∆ , θ ) H ( t ) ∝ ∆ β for small ∆ where β < − .
5. Both χ and the time of validity increase as the nodeis approached. Time advances at least as quickly as a geometric progression and this is analogous to the manner in whichthe open model steps towards its limit point. However, unlike the open case, once the trajectory crosses into the blue region(assuming it does not lie exactly on the unstable attracting trajectory) a single final step is needed. The specific number ofsteps will depend upon the starting initial conditions but will be small because of the property of geometric progression. LPT re-expansion can solve the problematic convergence in previously analysed open and closed models.Open models have asymptotic values of ∆ and θ and simple evolution. The first section below includes numerical resultsthat provide a practical demonstration of the success of LPT re-expansion in this case. Convergence as Lagrangian orderincreases and/or time step size decreases is observed qualitatively.Closed models have a somewhat more complex behaviour (before and after turnaround). The second section providesboth a qualitative and quantitative discussion of convergence. The scaling of the leading order error and the time step controlwhich are derived are of general applicability. Figure 16 investigates the effect of time step and order on the evolution of the open model introduced in figure 1 (∆ = 0 . θ = 2 . a = 10 − ). The series convergence breaks down at a = 0 . a = 1 using successively higher LPT series orders. As expected, higher order terms do not improve the accuracy ofthe description because the time of validity is violated. The middle panel employs three steps to reach a = 1, each respectingthe time of validity. Now the LPT series with higher order improves the accuracy just as one desires. The right panel employssix steps to reach a = 1, each respecting the time of validity. Again, higher order improves the description. Note that morefrequent re-expansion, i.e. smaller steps in time, improve the errors at fixed LPT order. c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff nd th th st rd th th th (cid:72) t (cid:76) b (cid:72) t (cid:76) (cid:72) t (cid:76) b (cid:72) t (cid:76) (cid:72) t (cid:76) b (cid:72) t (cid:76) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45)
101 a (cid:72) t (cid:76) L og (cid:200) b a pp r ox (cid:144) b e x ac t (cid:45) (cid:200) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45)
10 a (cid:72) t (cid:76) L og (cid:200) b a pp r ox (cid:144) b e x ac t (cid:45) (cid:200) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45)
10 a (cid:72) t (cid:76) L og (cid:200) b a pp r ox (cid:144) b e x ac t (cid:45) (cid:200) Figure 16.
LPT re-expansion of an open model with ∆ = 0 .
01 and θ = 2 .
82. The top three figures show the scale factor for the sameinitial conditions calculated with one step (left), three steps (middle) and five steps (right). The black dots indicate the position of thetime steps. In the middle and right panels, the solution was advanced 9 /
10 and 1 / b ( t ) including the unphysical negative ones. The order of the LPT expansion arecolour-coded according the top left figure. The single step expansion does not respect the time of validity whereas both the three andsix step examples do. The original expansion does not converge over the full time range whereas the re-expansions do. Coloured versiononline. Figure 17 investigates the closed model introduced in figure 8 (∆ = 0 . , θ = 0 . , a = 10 − ). The time of validity is determinedby a complex root. The first panel shows that the series begins to diverge at a = 0 .
38 well before the collapse singularity isreached at a = 0 . t = t and ending at t f has leading order error for the m -th order Lagrangian approxi-mation ∝ ( t f /t − m +2 ∆ m +1 , where ∆ is the value at the initial time. If the same small interval is covered in N smallersteps, the error after N steps scales as N − m ( t f /t − m +2 ∆ m +1 (see Appendix E for details). If the step size increases ina geometric sequence such that δt/t is a constant for each intermediate step, then t f = t (1 + δt/t ) N and the error after N steps scales as N ( t f /t − δt/t ) m +1 ∆ m +1 . This leads to the interpretation that the error per intermediate step scales as( δt/t ) m +1 ∆ m +1 . Define (cid:15) = ( δt/t )∆. The leading order error scales as (cid:15) m +1 which is numerically small if (cid:15) <
1. The sum ofall the missing higher order terms is finite if δt < T , i.e respects the time of validity.In a practical application, the initial and final times are not close. A reasonable time step criterion is to choose (cid:15) < δt for a given ∆. Other choices are possible but δt must always be less than thetime of validity. If (cid:15) is held fixed throughout the evolution, then the net error after N steps for the m -th order approximation ∝ (cid:15) m +1 N .The number of steps required to go from the initial to the final time can be estimated. As a special case assume that ∆ Typically, the numerical coefficient is of order unity and varies with m as well as the particular value of θ . For the purposes of adiscussion of the scaling of the error term, we assume the numerical coefficients to be constant as m and θ vary.c (cid:13) RAS, MNRAS000
1. The sum ofall the missing higher order terms is finite if δt < T , i.e respects the time of validity.In a practical application, the initial and final times are not close. A reasonable time step criterion is to choose (cid:15) < δt for a given ∆. Other choices are possible but δt must always be less than thetime of validity. If (cid:15) is held fixed throughout the evolution, then the net error after N steps for the m -th order approximation ∝ (cid:15) m +1 N .The number of steps required to go from the initial to the final time can be estimated. As a special case assume that ∆ Typically, the numerical coefficient is of order unity and varies with m as well as the particular value of θ . For the purposes of adiscussion of the scaling of the error term, we assume the numerical coefficients to be constant as m and θ vary.c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation is constant. The time step criterion implies that the number of steps to move from the initial time t = t to the final time t f for given (cid:15) is N = log( t f /t ) / log(1 + ( (cid:15)/ ∆)). For limited total intervals ( t f − t << t ) and small steps ( (cid:15)/ ∆ <<
1) theexact answer reduces to N ∼ ( t f − t )∆ /(cid:15) = ( t f − t ) /δt . Here δt = (cid:15)t ∆ does not grow appreciably over the interval so theestimate for N is a maximum. In this limit, the net error ∝ (cid:15) m ∆. The leading order error for the m-th order Lagrangianscheme decreases at least as quickly as (cid:15) m .In more general situations the value of ∆ varies. Once the closed model turns around ∆ increases without bound. Forfixed (cid:15) the step size δt decreases monotonically to zero as t → t coll where t coll is the time of the future singularity. At anyorder it would take infinitely many steps to follow the solution up until collapse. Consider the problem of tracking the solutionup to a large, finite value of ∆ = ∆ f . This moment corresponds to a fixed time t f < ∼ t coll in the exact solution. The numberof steps N < N max ∼ t f /δt f where δt f is the step size for the system near ∆ f ; δt f ∝ (cid:15)/ ∆ f . The leading order error after N steps at the m -th Lagrangian order ∝ (cid:15) m +1 N < (cid:15) m +1 N max ∼ (cid:15) m ∆ f . This method of step control forces the leading ordererror at fixed time t f < t coll to decrease as the Lagrangian order m increases and/or the control parameter (cid:15) decreases.The second and third panels in figure 17 show the runs with (cid:15) = 0 . (cid:15) = 0 . >
100 so as to avoid the infinite step regime. This required 32 steps for the (cid:15) = 0 . (cid:15) = 0 . m increases. This is true in a quantitative as well as qualitative sense. For example, in the second panel at a = 0 .
64 a plot ofthe log of the absolute error is approximately linear in m , as expected.(ii) A comparison of the same coloured lines in the middle and right panels shows that smaller (cid:15) implies better accuracy.Again, this is true in a quantitative as well as qualitative sense. For example, the observed ratio of errors at a = 0 .
64 for the9-th order calculations is 5 × − . To evolve up to this time with (cid:15) = 0 . (cid:15) = 0 . . / . (22 / ∼ × − , the same order of magnitude as theobserved ratio.These comparisons lead to the important conclusion that the leading order error for LPT re-expansion varies with Lagrangianorder and time step as theoretically expected.It is clear that considerable benefit accrues not only from implementing higher order Lagrangian schemes but also bylimiting time step size (which must always be less than the time of validity). For simple examples like the top-hat it is feasibleto work to very high Lagrangian order but this is not likely to be true in the context of more complicated, inhomogeneousproblems. On the other hand, marching forward by many small time steps using LPT re-expansion is generally feasible. Inthe example above the initial perturbation is ∆ = 0 . ∼ − . For the same (cid:15) the practical application requires more steps for the phase before turnaround but the netincrease is only a modest logarithmic factor. In fact, most of the steps in the example were taken after turnaround and thetotal number varies with the depth of the collapse. This will continue to be true for the practical calculation. The choice ofstep size and order for such applications will be the subject of a forthcoming paper. We have investigated the time of validity of Lagrangian perturbation theory for spherical top-hat cosmologies with generalinitial conditions. Using techniques from complex analysis we showed that the time of validity is always limited for openmodels. We also discovered a class of closed models whose time of validity is less than their time to collapse. We introducedthe concept of the mirror model and derived a symmetry principle for the time of validity of mirror models. For small initialperturbations the time of validity of LPT series expansion of an open model corresponds to the collapse time of a closedmirror model.A qualitative analogy is useful. A single LPT series expansion is similar to a single step in a finite difference approximationfor advancing a hyperbolic partial differential equation like the wave equation. The time of validity of the LPT expansion isanalogous to the Courant condition which guarantees stability. In LPT the constraint is an acceleration-related time-scale; inthe wave equation it is a sound-crossing time-scale.We developed the method of LPT re-expansion which overcomes the limitations intrinsic to a single expansion. Wedemonstrated how to iteratively re-expand the solution so as to link convergent series expressions that extend from initialto final times. The time of validity of the expansions set the minimum number of re-expansion steps ( ∼
10) necessary forcosmological simulations starting at recombination and proceeding to the present epoch. Finite as opposed to infinite order c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff nd th th st rd th th th th nd th th st rd th th th th nd th th st rd th th th th nd th th st rd th th th th nd th th st rd th th th th nd th th st rd th th th th nd th th st rd th th th th nd th th st rd th th th th nd th th st rd th th th th (cid:72) t (cid:76) b (cid:72) t (cid:76) (cid:72) t (cid:76) b (cid:72) t (cid:76) (cid:72) t (cid:76) b (cid:72) t (cid:76) (cid:45) (cid:45) (cid:45)
50 a (cid:72) t (cid:76) L og (cid:200) b a pp r ox (cid:144) b e x ac t (cid:45) (cid:200) (cid:45) (cid:45) (cid:45)
50 a (cid:72) t (cid:76) L og (cid:200) b a pp r ox (cid:144) b e x ac t (cid:45) (cid:200) (cid:45) (cid:45) (cid:45)
50 a (cid:72) t (cid:76) L og (cid:200) b a pp r ox (cid:144) b e x ac t (cid:45) (cid:200) Figure 17.
LPT re-expansion of a closed solution with ∆ = 0 . θ = 0 .
44. The top three figure show the scale factor calculated witha single step (left) and multiple steps with (cid:15) = 0 . (cid:15) = 0 . (cid:15) ). The bottom figuresshow the errors for all LPT approximations to b ( t ) including the unphysical negative ones. The order of the expansion is colour-codedas in the top left figure. The single step expansion does not respect the time of validity whereas both the other cases do. The black dotsindicate the position of the time steps. The original expansion does not converge over the full time range whereas the re-expansions do.Coloured version online. Lagrangian expansions required extra steps to achieve given error bounds. We characterised how the leading order numericalerror for a solution generated by LPT re-expansion varied with the choice of Lagrangian order and of time step size. Weprovided a recipe for time step control for LPT re-expansion based on these results.Our long-term goal and motivation for this study is to develop a numerical implementation of LPT re-expansion for fullyinhomogeneous cosmological simulation. Top-hats with Zeldovich initial conditions have special properties with respect toLPT convergence. We found that all underdense models must be treated by re-expansion while none of the overdense onesneed be. However, during the course of an inhomogeneous simulation the density and irrotational velocity perturbations (withrespect to a homogeneous background cosmology) at an arbitrary point will generally not fall on the top-hat’s Zeldovich curve.Hence, the convergence of LPT in inhomogeneous applications must be guided by the analysis of more general models. Top-hats with arbitrary initial conditions are the simplest possibility and constitute the main focus in this paper. The limitationson LPT convergence which we have elucidated in this generic case are considerably more complicated than in the top-hatwith Zeldovich initial conditions. Our plan is to use the generic time of validity criterion to determine the time-stepping forinhomogeneous evolution. This should allow us to develop high-precision simulations with well-defined control of errors. Thepractical impact of a refined treatment of LPT convergence is not yet clear.The convergence issues we have dealt with should not be confused with the breakdown when orbit crossing takes place andthe Jacobian of the transformation from Lagrangian to physical coordinates becomes singular. At that time the flow becomesmulti-streamed and much of the simplicity and advantage of the Lagrangian approach vanishes. The aim of the current workis to make sure it is possible to reach the epoch of multi-streamed flow but offers nothing new on how to proceed beyond it. Infact, it may be necessary to include an effective pressure term in the equations to account for the velocity dispersion induced byorbit crossing (Adler & Buchert 1999; Buchert, Dominguez & J.Perez-Mercader 1999) or to adopt alternative approximationsfor the basic dynamics (such as the adhesion approximation; see Sahni & Coles 1995 for a review and references therein) tomake progress. c (cid:13) RAS, MNRAS000
44. The top three figure show the scale factor calculated witha single step (left) and multiple steps with (cid:15) = 0 . (cid:15) = 0 . (cid:15) ). The bottom figuresshow the errors for all LPT approximations to b ( t ) including the unphysical negative ones. The order of the expansion is colour-codedas in the top left figure. The single step expansion does not respect the time of validity whereas both the other cases do. The black dotsindicate the position of the time steps. The original expansion does not converge over the full time range whereas the re-expansions do.Coloured version online. Lagrangian expansions required extra steps to achieve given error bounds. We characterised how the leading order numericalerror for a solution generated by LPT re-expansion varied with the choice of Lagrangian order and of time step size. Weprovided a recipe for time step control for LPT re-expansion based on these results.Our long-term goal and motivation for this study is to develop a numerical implementation of LPT re-expansion for fullyinhomogeneous cosmological simulation. Top-hats with Zeldovich initial conditions have special properties with respect toLPT convergence. We found that all underdense models must be treated by re-expansion while none of the overdense onesneed be. However, during the course of an inhomogeneous simulation the density and irrotational velocity perturbations (withrespect to a homogeneous background cosmology) at an arbitrary point will generally not fall on the top-hat’s Zeldovich curve.Hence, the convergence of LPT in inhomogeneous applications must be guided by the analysis of more general models. Top-hats with arbitrary initial conditions are the simplest possibility and constitute the main focus in this paper. The limitationson LPT convergence which we have elucidated in this generic case are considerably more complicated than in the top-hatwith Zeldovich initial conditions. Our plan is to use the generic time of validity criterion to determine the time-stepping forinhomogeneous evolution. This should allow us to develop high-precision simulations with well-defined control of errors. Thepractical impact of a refined treatment of LPT convergence is not yet clear.The convergence issues we have dealt with should not be confused with the breakdown when orbit crossing takes place andthe Jacobian of the transformation from Lagrangian to physical coordinates becomes singular. At that time the flow becomesmulti-streamed and much of the simplicity and advantage of the Lagrangian approach vanishes. The aim of the current workis to make sure it is possible to reach the epoch of multi-streamed flow but offers nothing new on how to proceed beyond it. Infact, it may be necessary to include an effective pressure term in the equations to account for the velocity dispersion induced byorbit crossing (Adler & Buchert 1999; Buchert, Dominguez & J.Perez-Mercader 1999) or to adopt alternative approximationsfor the basic dynamics (such as the adhesion approximation; see Sahni & Coles 1995 for a review and references therein) tomake progress. c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation ACKNOWLEDGMENTS
S. N. thanks Varun Sahni for discussions on convergence of Lagrangian theory at the IUCAA CMB-LSS summer school, PaulGrabowski, Sergei Dyda and Justin Vines for useful conversations and Saul Teukolsky for feedback on the manuscript. Theauthors would like to thank Thomas Buchert for useful comments on the paper. This material is based upon work supportedby the NSF under Grant No. AST-0406635 and by NASA under Grant No. NNG-05GF79G. c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff
REFERENCES
Adler S., Buchert T., 1999, A&A, 343, 317Bernardeau F., Colombi S., Gazta˜naga E., Scoccimarro R., 2002, Phys. Rep., 367, 1Bouchet F. R., 1996, astro-ph/9603013Bouchet F. R., Colombi S., Hivon E., Juszkiewicz R., 1995, A&A, 296, 575Bouchet F. R., Juszkiewicz R., Colombi S., Pellat R., 1992, ApJ, 394, L5Brown J. W., Churchill R., 1996, Complex Variables and Applications. WeilyBuchert T., 1992, MNRAS, 254, 729Buchert T., 1994, MNRAS, 267, 811Buchert T., 1995, astro-ph/9509005Buchert T., Dominguez A., J.Perez-Mercader 1999, A&A, 349, 343Buchert T., Ehlers J., 1993, MNRAS, 264, 375Catelan P., 1995, MNRAS, 276, 115Chicone C., 2006, Ordinary Differential Equations with Applications, (section 1.1, 1.12), second edn. Texts in AppliedMathematics, New York: Springer-VerlagEhlers J., Buchert T., 1997, General Relativity and Gravitation, 29, 733Karakatsanis G., Buchert T., Melott A., 1997, A&A, 326, 873Kasai M., 1995, Phys. Rev. D, 52, 5605Landau L. D., Lifshitz E. M., 1975, The Classical Theory of Fields, fourth edn. Butterworth HeinemannMatarrese S., Pantano O., Saez D., 1993, Phys. Rev. D, 47, 1311Matarrese S., Pantano O., Saez D., 1994, MNRAS, 271, 513Matarrese S., Terranova D., 1996, MNRAS, 283, 400Monaco P., 1997, MNRAS, 287, 753Moutarde F., Alimi J., Bouchet F. R., Pellat R., Ramani A., 1991, ApJ, 382, 377Munshi D., Sahni V., Starobinsky A. A., 1994, ApJ, 436, 517Sahni V., Coles P., 1995, Phys. Rep., 262, 1Sahni V., Shandarin S., 1996, MNRAS, 282, 641Scoccimarro R., Sheth R. K., 2002, MNRAS, 329, 629Tatekawa T., 2007, Phys. Rev. D, 75, 44028Tolman R. C., 1934, Proceedings of the National Academy of Science, 20, 169Zel’Dovich Y. B., 1970, A&A, 5, 84
APPENDIX A: FORMAL SET-UP OF THE SPHERICAL TOP-HAT
We intend to study an inhomogeneous universe. It contains a single, compensated spherical perturbation evolving in a back-ground cosmology. To describe two spatially distinct pieces of the inhomogeneous universe (the background and the centralperturbation) we invoke the language of homogeneous cosmology.
A1 Description of the background
The origin of the coordinate system is the centre of the sphere. The background system at the initial time t is set by thephysical size of the inner edge r b, , the velocity ˙ r b, and density parameter Ω . The Lagrangian coordinate system is extendedlinearly throughout space once the Lagrangian coordinate of the inner edge is fixed. Let the Lagrangian coordinate of theinner edge be Y = r b, a . (A1)Either choose the initial background scale factor a and determine the coordinate system or, alternatively, fix Y and inferthe background scale factor. In either case, the scale factor embodies the gauge freedom associated with the radial coordinatesystem.The future evolution of the inner edge of the background is given by r b ( t ) = a ( t ) Y . The velocity at the initial timesatisfies ˙ r b, = ˙ a Y . The density at any later time is ρ b ( t ) = ρ b a a , (A2)and the Hubble parameter for the background is c (cid:13) RAS, MNRAS000
The origin of the coordinate system is the centre of the sphere. The background system at the initial time t is set by thephysical size of the inner edge r b, , the velocity ˙ r b, and density parameter Ω . The Lagrangian coordinate system is extendedlinearly throughout space once the Lagrangian coordinate of the inner edge is fixed. Let the Lagrangian coordinate of theinner edge be Y = r b, a . (A1)Either choose the initial background scale factor a and determine the coordinate system or, alternatively, fix Y and inferthe background scale factor. In either case, the scale factor embodies the gauge freedom associated with the radial coordinatesystem.The future evolution of the inner edge of the background is given by r b ( t ) = a ( t ) Y . The velocity at the initial timesatisfies ˙ r b, = ˙ a Y . The density at any later time is ρ b ( t ) = ρ b a a , (A2)and the Hubble parameter for the background is c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation H = ˙ r b, r b, = ˙ a a . (A3)The evolution of the scale factor is¨ aa = − πGρ b a a = − H a Ω a (A4)The quantities, r b, , ˙ r b, , Ω and t along with the choice of the coordinate system, completely specify the background universe. A2 Description of the innermost perturbation
The perturbation can be described by four physical quantities: the physical position r p, and velocity ˙ r p, of the edge (or theratio H p = ˙ r p, /r p, ), the density parameter Ω p at the initial time t . The Lagrangian coordinate system for the perturbationis X = r p, b ( t ) . (A5)It can be linearly extended throughout space.Like a , b ( t ) embodies the gauge freedom associated with the choice of the coordinate system. Without loss of generality,one can pick this gauge to satisfy b ( t ) = a . (A6)Note that the Lagrangian coordinate systems for the background and perturbation are different.Let ρ and ρ p, denote the densities of the background and perturbation respectively. Define the perturbation parameters δ = ρ p ρ b − δ v = H p H − p = (1 + δ )(1 + δ v ) . (A9) A3 Inhomogeneous model
Figure A1 shows how an overdense and underdense innermost sphere may be embedded with compensation in a homogeneousbackground universe. The assumption that the background cosmology evolves like a homogeneous model, fully describedin terms of its Hubble constant and density, imposes consistency conditions. At the initial instant the “inner edge” of theunperturbed background distribution is at physical distance r b, from the centre of the sphere. The region with r > r b, willevolve like an unperturbed homogeneous cosmology as long as(i) the mass within equals the mass that an unperturbed sphere would contain;(ii) matter motions within the perturbed region do not overtake the inner edge of the homogeneous region.These conditions which are obvious in the Newtonian context have general relativistic analogues (Landau & Lifshitz 1975).Next, consider the innermost perturbed spherical region. At the initial time let r p, be the “outer edge” of this region.The physical properties and evolution of the innermost region are fully described in terms of its Hubble constant and densityas long as its outer edge does not overtake matter in surrounding shells. While this is obvious in a Newtonian context thereexists a relativistic analogue (Tolman 1934; Landau & Lifshitz 1975).The inhomogeneous model is incomplete without specification of the transition region between the innermost sphere andthe background. For the background to evolve in an unperturbed fashion the mass within r b, must be exactly 4 πρ r b, / δ > r p, < r < r b, so that ( ρ p /ρ ) = ( r b, /r p, ) = ( Y /X ) . The evolution of each matter-filled region proceeds independentlyas long as the trajectories of the inner and outer edges do not cross. When δ <
0, a more complicated transition is required.For example, one choice is to nest sphere, empty shell and dense shell (see figure A1) so that the mass within r b, matchesthat of the unperturbed background. In this case ρ p r p, = fρ r b, for some f < − f is placedin the dense shell). Varying the specifics of the compensation region while keeping the properties of the sphere fixed leaves δ and δ v , as defined above, invariant.For fixed δ and δ v the solution b ( t ) is independent of the details of the transition. Nonetheless, variation in f , r b, /r p, and Y /X all go hand-in-hand. Hence, the extent of time that the sphere’s evolution may be treated as independent of thematter-filled outer regions also varies. A basic premise of this paper is that it is meaningful to determine the limitations c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff r p,0 = b X = a X r b,0 = a Y Background VacuumOverdensePerturbation ! p,0 ! ! r p,0 = b X = a X r b,0 = a Y VacuumUnderdensePerturbation ! p,0 Figure A1.
A cartoon showing the physical set-up of the problem. arising from the convergence of the LPT series independently of limitations associated with crossing of separate matter-filledregions. For a given a δ and δ v this separation can be achieved for specific constructions by choosing the radius and (hencevelocity) of the inner sphere and the energy of the compensating region appropriately. A4 Number of degrees of freedom for the innermost sphere
If the innermost sphere corresponds to an overdensity then the compensating region can be a vacuum as shown in figure A1.Having picked the co-ordinate system, having selected equal initial times for the background and perturbation (not equalbang times but equal times at which we give the background and perturbation values), and required the correct amount ofmass, only two degrees of freedom remain: δ and δ v .To reiterate, the background and the perturbation can have different big bang times. Setting them equal would imply arelationship between δ and δ v and leave a single free parameter.If the innermost sphere corresponds to an underdensity then the compensating region is not vacuum but a spherical shell.In this case, in addition to δ and δ v , one must specify f or, equivalently, r p, . But the solution for b ( t ) is independent of thesize of the innermost sphere so, again, only two degrees of freedom remain. A5 Preventing shell crossing
There are two sorts of limitations for the solution of b ( t ). One is the calculation-dependent limitation arising from theconvergence properties of the Lagrangian series expansion. It involves the scale factors only. The other is a physical limitationarising from collisions of the innermost region with surrounding non-vacuum regions (either the background or a compensatingshell). We show that it is possible to delay the epoch of collisions indefinitely without altering the evolution of the innermostregion.Fix H , H p, , ρ and ρ p, . This implies that the expansion parameters in LPT, δ and δ v , and the time of validity of theLPT solution are all fixed. Consider the case of an overdensity surrounded by vacuum. To stave off the collision of the outeredge of the innermost region with inner edge of the homogeneous background hold r b, fixed and reduce r p, . The velocity˙ r p, = H p r p, becomes arbitrarily small. The time for the edge to reach any fixed physical distance increases without bound.Shell crossings may be put off indefinitely. However, we have altered the mass within the innermost edge of the backgroundso we add back a thin, dense shell just inside r b, and set it on a critical trajectory outward. This accomplishes our goal. c (cid:13) RAS, MNRAS000
There are two sorts of limitations for the solution of b ( t ). One is the calculation-dependent limitation arising from theconvergence properties of the Lagrangian series expansion. It involves the scale factors only. The other is a physical limitationarising from collisions of the innermost region with surrounding non-vacuum regions (either the background or a compensatingshell). We show that it is possible to delay the epoch of collisions indefinitely without altering the evolution of the innermostregion.Fix H , H p, , ρ and ρ p, . This implies that the expansion parameters in LPT, δ and δ v , and the time of validity of theLPT solution are all fixed. Consider the case of an overdensity surrounded by vacuum. To stave off the collision of the outeredge of the innermost region with inner edge of the homogeneous background hold r b, fixed and reduce r p, . The velocity˙ r p, = H p r p, becomes arbitrarily small. The time for the edge to reach any fixed physical distance increases without bound.Shell crossings may be put off indefinitely. However, we have altered the mass within the innermost edge of the backgroundso we add back a thin, dense shell just inside r b, and set it on a critical trajectory outward. This accomplishes our goal. c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation The case of the underdensity surrounded by a compensating shell is identical. First, we must make sure that the com-pensating shell does not overrun the homogeneous model. Choose the shell to be thin, fix its initial physical distance fromthe centre and adjust is velocity (based on how the interior mass changes) to give a critical solution. The two power laws,one for the compensating shell and one for the innermost boundary of the homogeneous model, cannot cross in the future.Second, as above, note that reducing r p, reduces the outward velocity of the edge so that it takes more time to reach theinitial position of the compensating shell. The time can be made arbitrarily long.The limitations in LPT convergence are completely distinct from those associated with physical collisions in inhomoge-neous model. APPENDIX B: SERIES EXPANSIONS FOR A FUNCTION OF TWO VARIABLES
In this section we elucidate by example some qualitative features of the expansion of b ( t, ∆), the central quantity in theLagrangian treatment of the top-hat. We assume a very simple form denoted f ( t, ∆) and look at convergence with respect toexpansions in t and ∆. Let f ( t, ∆) = t / (cid:16) t + ∆ (cid:17) / . (B1)The series expansion of this function around ∆ = 0 at fixed t is f ∼ t / + t / ∆ − t / ∆ + 581 t / ∆ − t / ∆ + 22729 t / ∆ + O (∆ ) (B2)which is supposed to mimic the Lagrangian expansion in ∆. One can also expand the function as a series in t around t = t i f ∼ (cid:114) ∆ + 1 t i t / i + (2∆ t i + 1)( t − t i )3 (cid:0) ∆ + t i (cid:1) / t / i + (cid:0) − ∆ t i − ∆ t i − (cid:1) ( t − t i ) (cid:0) ∆ + t i (cid:1) / t / i (∆ t i + 1)+ (cid:0) t i + 6∆ t i + 12∆ t i + 5 (cid:1) ( t − t i ) (cid:0) ∆ + t i (cid:1) / t / i (∆ t i + 1) + O (cid:0) ( t − t i ) (cid:1) . (B3)Both expansions involve the complex power z / . There are two branch cuts which extend to z = 0 so at ∆ = − /t thefunction is not analytic. Additionally, the expansion in t is not analytic at t = 0.The efficacy of various expansions are illustrated in figure B1. In all the plots the black dotted line indicates the exactfunction. The top left panel shows successively higher order series approximations in ∆ as a function of t for the specific case∆ = 1 /
10. The question here is whether the pole at a given time lies with a disk of radius 1 /
10? The location of the pole is∆ = − /t so the answer is “yes” when t >
10. This pole interferes with the convergence of the series expansion for ∆ = 1 / t < t = 1 /
10. The question here is how big a perturbation will convergeat t = 1 /
10? Since the location of the pole is ∆ = − /t the radius of convergence at the indicated time is 10. Perturbationswith | ∆ | >
10 are not expected to converge and the figure shows that this is indeed the case.The bottom left panel shows the series in t expanded around t i = 2 for fixed ∆ = 1 /
10. The poles are at t = −
10 and t = 0 in the complex t plane. The expected radius of convergence is min( | − | , | − ( − | ) = 2 or t i − < t < t i + 2. Asseen in the plot, the series converges only in the expected range (0 , t expanded around t i = 2 for ∆ = − /
3. The poles are at t = 3 and t = 0 inthe complex t plane. The expected radius of convergence is min( | − | , | − | ) = 1 or t i − < t < t i + 1. As seen in the plot,the series converges only in the expected range (1 , APPENDIX C: PARAMETRIC SOLUTION
The background model has scale factor a and Hubble constant H = ˙ a /a . The model, perturbed in density and velocity,is parameterized by ∆ and θ and has scale factor b ( t ). For the choice of coordinate system given in the text the second orderequation for b is¨ bb = − H a (1 + ∆ cos θ ) b (C1)with the initial conditions that at t = t , b ( t ) = a , ˙ b ( t ) = ˙ a (1 + ∆ sin θ ). The scale factor a and the velocity of thebackground ˙ a at the initial time t are positive. The parametrization of ˙ b ( t ) allows either positive or negative values where∆ is non-negative and − π < θ (cid:54) π . The quantity (1 + ∆ cos θ ), proportional to total density, is non negative.This equation once integrated is c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff f (cid:72) t , (cid:144) (cid:76) (cid:45) (cid:45) (cid:45)
10 0 10 20 300.00.20.40.60.8 (cid:68) f (cid:72) (cid:144) , (cid:68) (cid:76) (cid:45) (cid:45) (cid:45) (cid:45) f (cid:72) t , (cid:45) (cid:144) (cid:76) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) f (cid:72) t , (cid:45) (cid:144) (cid:76) Figure B1.
Series expansions in t and ∆ for an illustrative function f ( t, ∆) (see text). The black dotted line indicates the exact function f and the blue solid lines indicate successive approximations. The top left and right panels are series expansions in ∆ around ∆ = 0plotted as a function of t (for ∆ = 1 /
10) and function of ∆ (for t = 1 /
10) respectively. The bottom left and right panels are seriesexpansions in the t around t = 2 plotted as functions of t for ∆ = − /
10 and ∆ = − / ˙ b = H a (cid:20) (1 + ∆ cos θ ) b + (1 + ∆ sin θ ) − (1 + ∆ cos θ ) a (cid:21) . (C2)The combination E (∆ , θ ) = (1 + ∆ sin θ ) − (1 + ∆ cos θ ) (C3)is proportional to the total energy and determines the fate of the system. If E (∆ , θ ) >
0, the model is open and if E (∆ , θ ) < E , positive and negative ˙ b ) are shownin figure 5. C1 Initially Expanding Solutions
The expanding case with ˙ b > E >
0) has solution b ( η, ∆ , θ ) = a θ ) E (∆ , θ ) (cosh η −
1) (C4) t ( η, ∆ , θ ) = 12 H (1 + ∆ cos θ ) E (∆ , θ ) / (sinh η − η ) + t + bang (∆ , θ ) (C5)and the singularity b = 0 occurs at η = 0. For closed models ( E <
0) the solution is b ( η, ∆ , θ ) = a θ ) | E (∆ , θ ) | (1 − cos η ) (C6) t ( η, ∆ , θ ) = 12 H (1 + ∆ cos θ ) | E (∆ , θ ) | / ( η − sin η ) + t + bang (∆ , θ ) . (C7)For closed models, the convention adopted sets η = 0 at the singularity nearest in time to t . For both models, the time at η = 0 is denoted t + bang . For closed models the time at η = 2 π is denoted t + coll .At the initial time the solutions (both open and closed) satisfy b ( t ) = a , ˙ b ( t ) = ˙ a (1+∆ sin θ ) and t = t . The condition b ( t ) = a sets the value of the parameter at the initial time η . The velocity condition is then manifestly satisfied from theform of eq. (C2). The condition t = t at η = η sets the value of the bang time c (cid:13) RAS, MNRAS000
0) the solution is b ( η, ∆ , θ ) = a θ ) | E (∆ , θ ) | (1 − cos η ) (C6) t ( η, ∆ , θ ) = 12 H (1 + ∆ cos θ ) | E (∆ , θ ) | / ( η − sin η ) + t + bang (∆ , θ ) . (C7)For closed models, the convention adopted sets η = 0 at the singularity nearest in time to t . For both models, the time at η = 0 is denoted t + bang . For closed models the time at η = 2 π is denoted t + coll .At the initial time the solutions (both open and closed) satisfy b ( t ) = a , ˙ b ( t ) = ˙ a (1+∆ sin θ ) and t = t . The condition b ( t ) = a sets the value of the parameter at the initial time η . The velocity condition is then manifestly satisfied from theform of eq. (C2). The condition t = t at η = η sets the value of the bang time c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation t + bang = t − (cid:40) H (1+∆ cos θ ) | E (∆ ,θ ) | / ( η − sin η ) E < H (1+∆ cos θ ) E (∆ ,θ ) / (sinh η − η ) E > . (C8)The bang time for the model can also be written as t + bang = t − (cid:90) b = a b =0 db (˙ b ) (1 / , (C9)where ˙ b is given by eq. (C2) with the sign for the square root positive. The age of the model since its birth is t age (∆ , θ ) = (cid:90) b = a b =0 db (˙ b ) (1 / = (cid:90) η = η η =0 db/dη · dη (˙ b ( η )) (1 / . (C10)Inserting the appropriate parametric solution, one can verify that the bang times obtained from (C8) and (C9) are identical.Generally t + bang (cid:54) = 0.The velocity at the initial time is˙ b = ˙ a | E | / (cid:26) sin η − cos η E < sinh η cosh η − E > . (C11)First, ˙ b > η >
0. Second, if the age of the model increases, η increases. For the open solution if η varies from 0 to ∞ time increases from t + bang to ∞ . For a single cycle of the closed solutions, η increases from 0 to 2 π and time increases from t + bang to t + coll .In summary, the parametric solutions solve eq. (C1) and eq. (C2) for the specified initial conditions. As a final usefulstep, rewrite eq. (C9) by defining y = b/a t + bang = t − H (cid:90) y =1 y =0 dy [(1 + ∆ cos θ ) y − + E (∆ , θ )] / (C12)which follows from eq. (C2) and uses the same positive square root convention. C2 Initially Contracting Solutions
Next, consider the case ˙ b <
0. The parametric solution for
E > b ( η, ∆ , θ ) = a θ ) E (∆ , θ ) (cosh η −
1) (C13) t ( η, ∆ , θ ) = 12 H (1 + ∆ cos θ ) E (∆ , θ ) / ( − sinh η + η ) + t − bang (∆ , θ ) (C14)and for E < b ( η, ∆ , θ ) = a θ ) | E (∆ , θ ) | (1 − cos η ) (C15) t ( η, ∆ , θ ) = 12 H (1 + ∆ cos θ ) | E (∆ , θ ) | / ( − η + sin η ) + t − bang (∆ , θ ) . (C16)Again, for closed models, the convention adopted is that the singularity nearest to t corresponds to η = 0. The time at η = 0is t − bang and the collapse time for closed models is t − coll .The parametric form of the solutions satisfies eq. (C1) and eq. (C2). Just as in the previous case, the initial conditionsset η and t − bang . Since the singularity at η = 0 lies to the future of t , t − bang = t + (cid:90) b = a b =0 db (˙ b ) (1 / . (C17)where, ˙ b is given by eq. (C2). The sign of the square root is chosen to be positive and the integral is a positive quantitywhich is added to t . For closed models the singularity at η = 2 π lies to the past of t at t − coll . In this case (see figure 5) thelabelling implies t − coll < t < t − bang . Although this might seem backwards, it facilitates combining the open and closed modelsinto one complex function as was done in the positive ˙ b case. The initial velocity is˙ b = ˙ a | E | / (cid:26) sinh η − cosh η E > sin η cos η − E < b < η >
0. For the age of the model to increase, η must decrease. Conversely, if η increases,the time in the open model decreases from t to −∞ and the time in the closed model decreases from t to t − coll .A table summarising the properties of the physical solutions with η = | η | ζ = ηζ follows. c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff
Closed Open If η increases If t increases t bang − t ˙ b > ζ = 1 ζ = i t increases from t + bang η increases to ∞ or 2 π < b < ζ = 1 ζ = i t decreases from t − bang η decreases to 0 > C3 Analytic Extension of the exact solution in parametric form
The differential eq. (32) was solved numerically over the range 0 (cid:54) θ (cid:54) π , 0 < δ <
100 and − π < φ (cid:54) π where ∆ = ∆ e iφ . Foreach value of (∆ , φ, θ ), the numerical solution matched one of the two possible parametric forms.Omitting the explicit functional dependence on ∆ and θ the following abbreviations are useful j = (1 + ∆ cos θ ) (C19) h = (1 + ∆ sin θ ) j (C20) E = ( h − j . (C21)The two possible parametric forms that agree with the numerical solution are b ( η ) = a j [ − E ] (1 − cos η ) (C22) t ( η ) = t ± (cid:18) H j [ − E ] / ( η − sin η ) − t age (cid:19) (C23)where t age = 1 H (cid:18)(cid:112) j √ h − j [ E ] / sinh − (cid:114) Ej (cid:19) . (C24)The branch cut lies along the negative real axis for all fractional powers and from − i ∞ to − i and + i to i ∞ for the inversesinh function.The prescription for the correct form is for the choice of the ± sign in t eq. (C23) and denoted t + and t − . The correctform depends upon θ , φ , arg[ h ] (the arg is defined to be between − π and π ) and the (real) value j = j when φ = 0 or π . Thefigure C1 shows the upper half plane for the perturbation partitioned into areas where the complex extension of the solutionhas one of two forms. The lower half plane has the same structure inverted through the origin. The horizontal red dashed linedenotes ∆ sin θ = 1 and the vertical red dashed lines denote ∆ cos θ = ±
1. In some areas a single form applies as marked butin the central area both occur. The detailed prescription is t = (cid:54) θ (cid:54) π/ φ = π, | ∆ | sin θ < j < t − otherwise t + π/ < θ (cid:54) π < φ < π and arg h > t + − π < φ < h < t + φ = 0 and cos θ > θ < j > t + φ = π and cos θ < θ > j < t − otherwise t − (C25) c (cid:13) RAS, MNRAS000
1. In some areas a single form applies as marked butin the central area both occur. The detailed prescription is t = (cid:54) θ (cid:54) π/ φ = π, | ∆ | sin θ < j < t − otherwise t + π/ < θ (cid:54) π < φ < π and arg h > t + − π < φ < h < t + φ = 0 and cos θ > θ < j > t + φ = π and cos θ < θ > j < t − otherwise t − (C25) c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation t (cid:43) t (cid:43) t (cid:43) except at Φ(cid:61)Π t (cid:43) except at Φ(cid:61) t (cid:43) t (cid:45) Both t (cid:43) and t (cid:45) (cid:45) (cid:45) ∆ (cid:61) (cid:68) Cos Θ ∆ v (cid:61) (cid:68) S i n Θ Figure C1.
This figure describes one aspect of the analytic extension of the exact solution. For a given real ∆, the complex extension∆ → ∆ e iφ obeys eq. (C25) with two possible forms t + and t − . The choice depends on φ , ∆, θ . For some (∆, θ ) a single form is sufficientfor all φ ; for other values both forms are needed. This figure illustrates how the upper half plane is partitioned based on this property.Coloured version online. APPENDIX D: NUMERICAL SOLUTIONSD1 Algorithm
The initial conditions are parameterized by ∆ > − π < θ (cid:54) π . The transformation θ → π ± θ and ∆ → − ∆ leaves thesolution unchanged. At any time the roots for θ and π ± θ are negatives of each other. The root plot only depends upon theabsolute value of the root so the plots for θ and π ± θ are identical. It is sufficient to consider the upper half plane.For a given θ the algorithm to map out R ∆ ( t ) is the following: Vary ∆ from 0 to an arbitrarily large value ( ∼ ∆ = ∆ e iφ by varying the phase angles φ over the range 0 to 2 π . For each ∆ evaluate t ( η ) at η = 0 and η = 2 π calculated according to eq. (C25). Finally, hunt for solutions that set the imaginary part of t to 0.This last step involves one-dimensional root-finding in φ at fixed ∆. A solution leads to a specific pair ( t, ∆ ) that is a pole inthe function b ( ∆ , t ).Roots with t > t limit future evolution; those with t < t limit backwards evolution. Both sets are shown in the results.Roots are classified based on whether they are real or complex. For closed models the real roots can represent a singularitythat is nearby ( η = 0) or far away ( η = 2 π ) from t . This classification at the initial time is independent of whether thesingularity is in the past or future and is independent of whether the model is expanding or contracting. For open models thereal roots are always considered nearby ( η = 0).In what follows the numerical answers are first described in qualitative terms. In the next section simple analytic estimatesfor the time of validity are developed.Figure D1 shows the root plots on a log-log scale. Sixteen panels, each with a particular value of θ listed at the top, aredisplayed. The x-axis is log H t and the y-axis is the log of the distance of the singularity from the origin in the complex ∆ plane. The initial time, H t = 2 /
3, is marked by the vertical black dashed line.For each θ , the shaded region indicates the range of ∆ that gives rise to closed models. Figure 2 shows that closed modelsoccur only for θ < θ + c = 0 .
463 in the upper half plane so only some of the root plots have shading and then only at smaller ∆.The colour coding of the dots indicates four types of roots: real and complex roots where η = 2 π are in blue and red,respectively; real and complex roots with η = 0 are in cyan and pink, respectively. The radius of convergence at the initialtime t is infinite, i.e. the Lagrangian series is exact at the initial time by construction. At times very close to the initialtime the root loci lie off the plot. Only the roots to the right of H t are relevant for forward evolution and, conversely, onlythose to the left are relevant for backwards evolution. The discussion is focused on the case of forward evolution but it isstraightforward to consider the restrictions on marching backwards in time.The phase of the root (of smallest magnitude) appears in figure D3. When closed models have real roots they are positive;when open models have real roots they are negative. However, some open and closed models also possess complex roots. Theset of models with complex roots (of smallest magnitude) is evident from the shading in figure 10. The phase of each root ofsmallest magnitude in figure D1 is indicated by the colour shading in figure D3.There are horizontal dashed lines with colours green, blue and purple in figures D1 and D3 indicating | δ v | = 1, | δ | = 1 c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff
Θ (cid:61) (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Figure D1.
Root plots for θ in the range 0 (cid:54) θ (cid:54) π . In each plot the abscissa is log H t and the ordinate is the logarithm of themagnitude of the root. The vertical black dashed line marks the initial time. The shaded area corresponds to closed models. The blueand red points show real and complex roots with η = 2 π , respectively. The cyan and pink show real and complex roots with η = 0,respectively. The green and purple dashed lines are | δ v | = 1 (∆ = | sin θ | − ) and | δ | = 1 (∆ = | cos θ | − ), respectively. The blue dashedline indicates the switch between two forms and a single form of the parametric solution at ∆ = | θ − csc θ | . Coloured version online. and the transition between one and two complex forms, respectively. For each θ the lines mark the implied, special value of∆. These dashed lines also appear with the same colour coding in figure D2.The roots in figure D1 will be analysed in the range 0 < θ (cid:54) π/ π/ < θ (cid:54) π/ π/ < θ (cid:54) π . D1.1 (cid:54) θ (cid:54) π/ θ = 0; the blue dots indicate real roots with η = 2 π ; the blue shading indicates a closedmodel; the phase is positive (top left panel in figure D3). Only a single branch is evident. In sum, each root is the collapsetime of a closed, pure density perturbation. For an expanding model, η = 2 π implies that the root is the future singularity.That θ = 0 is a special case can be seen by consulting figure D2: a ray starting at ∆ = 0 with θ = 0 never intersects any ofthe other lines of the diagram. In general, each time a ray crosses one of the lines there is qualitative change in the propertiesof the roots. c (cid:13) RAS, MNRAS000
Root plots for θ in the range 0 (cid:54) θ (cid:54) π . In each plot the abscissa is log H t and the ordinate is the logarithm of themagnitude of the root. The vertical black dashed line marks the initial time. The shaded area corresponds to closed models. The blueand red points show real and complex roots with η = 2 π , respectively. The cyan and pink show real and complex roots with η = 0,respectively. The green and purple dashed lines are | δ v | = 1 (∆ = | sin θ | − ) and | δ | = 1 (∆ = | cos θ | − ), respectively. The blue dashedline indicates the switch between two forms and a single form of the parametric solution at ∆ = | θ − csc θ | . Coloured version online. and the transition between one and two complex forms, respectively. For each θ the lines mark the implied, special value of∆. These dashed lines also appear with the same colour coding in figure D2.The roots in figure D1 will be analysed in the range 0 < θ (cid:54) π/ π/ < θ (cid:54) π/ π/ < θ (cid:54) π . D1.1 (cid:54) θ (cid:54) π/ θ = 0; the blue dots indicate real roots with η = 2 π ; the blue shading indicates a closedmodel; the phase is positive (top left panel in figure D3). Only a single branch is evident. In sum, each root is the collapsetime of a closed, pure density perturbation. For an expanding model, η = 2 π implies that the root is the future singularity.That θ = 0 is a special case can be seen by consulting figure D2: a ray starting at ∆ = 0 with θ = 0 never intersects any ofthe other lines of the diagram. In general, each time a ray crosses one of the lines there is qualitative change in the propertiesof the roots. c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation (cid:68)(cid:61) (cid:200) Θ(cid:45) cos Θ (cid:200) (cid:68)(cid:61) (cid:200) Θ(cid:45) cos Θ (cid:200) (cid:68) sin Θ(cid:61) (cid:68) sin Θ(cid:61)(cid:45) (cid:68) cos (cid:61) (cid:68) cos (cid:61)(cid:45) (cid:68) E (cid:61) (cid:68) rc (cid:68) rc P (cid:45) (cid:45) (cid:45) (cid:45) ∆ ∆ v Figure D2.
Several conditions determine the nature of the roots in phase space. The most significant are schematically illustrated here.The green horizontal lines are | δ v | = ∆ | sin θ | = 1; purple vertical lines are | δ | = ∆ | cos θ | = 1; the black curved line is the E = 0 criticalsolution. The red lines ∆ rc mark where real roots associated with closed models (or closed mirror models) transform to complex roots.The blue dashed lines mark the division between one and two complex forms (see also figure C1). Physical models lie to the right of δ = −
1. Expanding models lie above δ v = −
1. The intersection δ = δ v = 1 occurs at θ = π/
4. The point P near θ = 0 .
84 is the meetingof δ v = 1 and ∆ rc . Coloured version online. For 0 < θ (cid:54) π/ < θ < θ + c in figure D2 (tan θ + c = 1 / E = 0 line at the origin). Eventually such a line will cross the blackline which is the E = 0 critical solution labelled ∆ E =0 . For small ∆ the models are closed; for larger ∆ they are open. Infigure D1 this distinction corresponds to the the blue shading (closed models) at small ∆ versus the unshaded (open) modelsat large ∆.Within the shaded region note that two branches of real roots are present beyond a given time; at large t (asymptotically)the lower branch is ∆ → → ∆ E =0 . The lower branch sets the time of validity for small ∆. Eachroot is the collapse time of a closed model which has both density and velocity perturbations at the initial time.As ∆ increases the time of validity inferred from the lower branch decreases. At the critical point ∆ = ∆ rc , the two realbranches merge and connect to a branch of complex roots (intersection of red and blue points). For ∆ > ∆ rc , the complexroots determine the time of validity even though the upper branch provides a real root. The complex roots do not have a c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff direct physical interpretation in terms of future singularities of physical models. On figure D2 the ray emanating from theorigin at shallow angle crosses the red dashed line labelled ∆ rc at this critical point.Physically, when ∆ exceeds ∆ rc , the velocity perturbation dominates the density perturbation in the sense that thecollapse time begins to increase. The real root corresponds to the future singularity of the model. As ∆ increases further, thesolution eventually becomes critical (infinite collapse time). The particular value where this occurs is ∆ E =0 and it correspondson figure D2 to the ray crossing the labelled black line. Within the entire range ∆ rc < ∆ < ∆ E =0 the complex root determinesthe time of validity. So, even though any model in this range is closed and possesses a real future singularity, the time ofvalidity is determined by the complex root. This gives the sliver on figure 10 which is the overlap of light red and blue shadings.Both ∆ rc and ∆ E =0 decrease as θ → θ + c as is evident from figure D2 and both vanish at θ + c . On figure D1 the real rootscompletely disappear and only the complex roots are present, i.e. the two real branches have been pushed out to infinite times.The panel with θ = 5 π/
36 = 0 .
436 is numerically closest to the critical case θ + c = 0 .
464 and the real branches are just barelyvisible at the right hand edge.For the rest of the upper half plane θ + c < θ (cid:54) π the ray no longer intersects any closed models.For θ + c < θ < π/ θ = π/ t the two branches have ∆ → → ∆ E =0 andare completely analogous to the real branches just discussed for closed models. The separation between the two real branchesincreases as θ increases and the solution loci shifts upwards in ∆. And just as before the two branches join and meet a complexbranch. The second red dashed line ∆ rc in figure D2 shows the real to complex transition for the roots for the open models.This behaviour might be expect to continue for π/ < θ < π but there is an additional complication: the analyticextension involves two forms. As the ray sweeps counterclockwise in figure C1 it crosses δ v = 1 (horizontal dashed line andthe curved blue line. These are also schematically illustrated in figure D2. D1.2 π/ < θ < π/ θ < π/ θ increases, the point where ∆ rc meets δ v = 1 islabelled P.For a fixed θ consider increasing ∆ from small values near the origin to ∞ . The order in which this ray intersects thegreen ( δ v = 1), purple ( δ = 1), red (∆ rc ) and blue (one or two complex forms) curves will correlate with the change in roots.The roots are negative real for small ∆. They correspond to the collapse time of a closed mirror model. Increase ∆and ignore ∆ rc . When the δ v = 1 line is crossed, the sign of the closed mirror model’s velocity switches from expanding tocontracting. This just means that the labelling of the future singularity switches from further away ( η = 2 π ) to nearer ( η = 0).Now recall ∆ < ∆ rc implies real roots and, by definition, δ v = ∆ sin θ . Hence, ∆ rc ( θ ) > / sin θ implies that the label switchoccurs just as outlined. On figure D2 rays counterclockwise of point P belong to this case. This is responsible for the switchfrom blue (real η = 2 π ) to cyan (real η = 0) roots at the green line in figure D1 for θ = π/ π/ rc ( θ ) < / sin θ the roots are already complex and the label switch occurs between the correspondingcomplex roots. There are no pictured examples in figure D1.In the previous section with the 0 < θ (cid:54) π/
4, the physical interpretation of ∆ rc (as ∆ increases) was that the velocitycontribution to the perturbation became dominant in the original model if the model was closed or in the mirror closed modelif the original model was open. In latter case the mirror models were initially expanding. Now, the same idea continues toapply in the regime π/ < θ < .
84. Here the transition from real to complex roots occurs before the δ v = 1 line is crossed.The significance of ∆ rc is that it marks the increasing importance of velocity perturbations in the closed expanding mirrormodels.However, for 0 . < θ < π/ δ v = 1, the mirror model swaps from η = 2 π to 0and the roots (real) corresponds to the real future singularity of a closed, contracting model. As ∆ increases further, first themirror model becomes critical and then an open model contracting to a future singularity. While the magnitude of ∆ growslarger than a critical value the velocity perturbation dominates the mirror model dynamics. When ∆ > ∆ rc the roots switchfrom real to complex. At this point the contracting mirror model can be open, closed or critical.Note in figure D2 that ∆ rc asymptotes to the vertical purple line δ = 1. The corresponding mirror model hits the line δ = − θ < ∼ π/ η = 0 to η = 2 π . The roots remain complex and since there is no physical interpretationand it is irrelevant whether they belong to η = 0 or η = 2 π . c (cid:13) RAS, MNRAS000
84. Here the transition from real to complex roots occurs before the δ v = 1 line is crossed.The significance of ∆ rc is that it marks the increasing importance of velocity perturbations in the closed expanding mirrormodels.However, for 0 . < θ < π/ δ v = 1, the mirror model swaps from η = 2 π to 0and the roots (real) corresponds to the real future singularity of a closed, contracting model. As ∆ increases further, first themirror model becomes critical and then an open model contracting to a future singularity. While the magnitude of ∆ growslarger than a critical value the velocity perturbation dominates the mirror model dynamics. When ∆ > ∆ rc the roots switchfrom real to complex. At this point the contracting mirror model can be open, closed or critical.Note in figure D2 that ∆ rc asymptotes to the vertical purple line δ = 1. The corresponding mirror model hits the line δ = − θ < ∼ π/ η = 0 to η = 2 π . The roots remain complex and since there is no physical interpretationand it is irrelevant whether they belong to η = 0 or η = 2 π . c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation In figure D1 the panel with θ = π/ π/
36 show these transitions: the blue to cyan transition at the green dashedline is the mirror model switch from expanding to contracting; the cyan to pink transition at the purple dashed line is themirror model moving through δ = −
1; the pink to red transition is the switch from two to one complex roots and η = 0 to η = 2 π . D1.3 θ = π/ θ = π/
2, only real roots of η = 0 are present for large ∆. This is a special case in that a ray only intersects one specialline δ v = 1 in the upper half plane. D1.4 π/ < θ (cid:54) π All models in this range also correspond to open models. Like the previous cases, small ∆ have real, negative roots with η = 2 π . The mirror models in this case lie in the fourth quadrant. The crossover of real roots from η = 0 to η = 2 π occurs at δ v = 1, however, unlike in the earlier case, the line δ = − η = 2 π roots for small∆ are collapse times of initially expanding closed mirror models and the η = 0 are future singularities of initially contractingclosed and open mirror models for intermediate and large values of ∆ respectively. D2 Numerical Results
Here we present numerical formulas that give the time of validity for any initial ∆ and θ . Real roots occur for small ∆ when0 < θ < π/
2; and they occur for all ∆ when π/ (cid:54) θ (cid:54) π or θ = 0. Real roots correspond to past or future singularities ofphysical models and are known exactly.Figure D1 shows that complex roots occur 0 < θ < π/
2. In the range π/ < θ < π/ π . We can approximate these roots as real, negative roots. Conversely, figure D3 also showsthat in the range 0 < θ (cid:54) π/ π . These roots are complex only when ∆ > ∆ rc . First, we fit ∆ rc by∆ rc,app ( θ ) = (cid:12)(cid:12) .
41 csc θ (cos θ − θ ) + 3 . θ − θ )(sin θ ) . (cid:12)(cid:12) . (D1)We cannot approximate the time of validity with the results for physical cases but it turns out that the numerically derivedtime of validity is insensitive to θ in the range 0 < θ < π/ H t app (∆) = 23 + . . < ∆ (cid:54) . . < ∆ (cid:54) . . < ∆ (cid:54) . < ∆ (cid:54) × ∆ . ∆ > . (D2)Using these quantities, the table below gives an approximation to the time of validity, T app , for all values of θ and ∆. Thetimes for collapse and the bang times are equivalent to eq. (C23) and reproduced here for convenience: t coll (∆ , θ ) = t + 12 H (1 + ∆ cos θ )[ − E (∆ , θ )] / (2 π ) − t age (∆ , θ ) (D3) t − bang (∆ , θ ) = t + t age (∆ , θ ) (D4)where, t age (∆ , θ ) = 1 H (cid:112) (1 + ∆ cos θ ) (cid:114) (1 + ∆ sin θ ) (1 + ∆ cos θ ) − H (1 + ∆ cos θ )[ E (∆ , θ )] / sinh − (cid:114) E (∆ , θ )(1 + ∆ cos θ ) , (D5) E (∆ , φ, θ ) = (1 + ∆ sin θ ) − (1 + ∆ cos θ ) (D6)The error in the fit is estimated as E = T − T app T . (D7)If E > E <
0, then the approximation overestimates the time of validity. Using the above fits the worst case is
E (cid:39) − . δt which satisfies δt < . T app so that the inaccuracy in the approximation is irrelevant. c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff
Table D1.
Approximation to time of validity, T app (∆ , θ ), for 0 (cid:54) θ (cid:54) π . Note that ∆ rc,app is an approximation to ∆ rc in eq. (D1).Parameter range T app < ∆ < ∆ rc,app E (∆ , θ ) < t coll (∆ , θ )0 < θ < π/ < ∆ < ∆ rc,app E (∆ , θ ) > t coll ( − ∆ , θ )∆ > ∆ rc,app t app (∆ , θ )0 < ∆ < | sin θ | t coll ( − ∆ , θ ) π/ < θ (cid:54) π/ | sin θ | (cid:54) ∆ (cid:54) (cid:12)(cid:12) θ − cos θ sin θ cos θ (cid:12)(cid:12) (cid:60) (cid:2) t − bang ( − ∆ , θ ) (cid:3) ∆ > (cid:12)(cid:12) θ − cos θ sin θ cos θ (cid:12)(cid:12) (cid:60) [ t coll ( − ∆ , θ )]0 < ∆ (cid:54) | sin θ | t coll ( − ∆ , θ ) π/ (cid:54) θ (cid:54) π ∆ > | sin θ | t − bang ( − ∆ , θ ) APPENDIX E: ERROR CHARACTERISATION OF THE LAGRANGIAN SERIES
We want to characterise the errors associated with calculating the solution at time t f given some fixed initial conditions attime t . Errors arise because any real calculation involves truncating the Lagrangian expansion. We want to compare the errorsthat result from different choices of truncation order and of the number of re-expansion steps assuming all series expansionsare convergent (i.e. all respect the time of validity). Let N m represent the final physical coordinate generated with a m -thorder calculation using N steps. Ultimately, we seek to characterise differences like N m − N (cid:48) m (cid:48) . The quantity 1 ∞ is the exactanswer. E0.1 Single step
The Lagrangian series solution for a single step has the form1 ∞ = a ( t ) (cid:32) ∞ (cid:88) i =1 b ( i ) ( t ) a ( t ) ∆ i (cid:33) X , (E1)where each b ( i ) satisfies¨ b ( i ) − H a b ( i ) a = S ( i ) (E2)and initial conditions are specified at t = t . The initial conditions at each order and the forms for the first few S ( i ) are givenin the text.If t f − t = δt << t , then the solutions can be expanded in the small parameter δt/t . The solutions are b (1) ( t ) /a ( t ) ∼ c (1) δtt , (E3) b ( i ) ( t ) /a ( t ) ∼ c ( i ) (cid:16) δtt (cid:17) i +1 for i (cid:62) . (E4)The coefficients c ( i ) depend on the angle θ and have a weak dependence on the Lagrangian order. The difference between theexact answer and the m -th order approximation for a single step is1 ∞ − m = (cid:32) ∞ (cid:88) i = m +1 c ( i ) (cid:16) δtt (cid:17) i +1 ∆ i (cid:33) X . (E5)As long as t f is within the time of validity of LPT, by definition, the LPT series converges and from the equation above, theleading order error scales as ∼ ( δt/t ) m +2 ∆ m +1 (order terms first by powers of ∆ and then by powers of δt/t ). E0.2 Multiple steps
In general for a practical application one is limited to working at a finite Lagrangian order. In such cases, it becomes necessaryto ask if convergence can be achieved by working at a finite Lagrangian order with increasing number of steps.First, we outline the calculation. The initial data is subscripted by “0”. For example, let the initial perturbed scale factorbe b = b ( t ), the initial background scale factor a , the initial density contrast δ and the initial velocity perturbation δ v .The Lagrangian expansion parameter ∆ and angle θ follow from the relations δ = ∆ cos θ and δ v = ∆ sin θ . Thephysical coordinate is r = b X ; for given r the initial Lagrangian coordinate X is fixed by choosing b to be equal to a . c (cid:13) RAS, MNRAS000
In general for a practical application one is limited to working at a finite Lagrangian order. In such cases, it becomes necessaryto ask if convergence can be achieved by working at a finite Lagrangian order with increasing number of steps.First, we outline the calculation. The initial data is subscripted by “0”. For example, let the initial perturbed scale factorbe b = b ( t ), the initial background scale factor a , the initial density contrast δ and the initial velocity perturbation δ v .The Lagrangian expansion parameter ∆ and angle θ follow from the relations δ = ∆ cos θ and δ v = ∆ sin θ . Thephysical coordinate is r = b X ; for given r the initial Lagrangian coordinate X is fixed by choosing b to be equal to a . c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation Consider taking N steps from initial to final time with an m -th order Lagrangian expansion. Assume that the final timeis within the time of validity of the Lagrangian expansion. For definiteness, let the j -th time be t j = t β j/N where β = t f /t (so t N is just the final time t f ). This geometric sequence of increasing steps is well-suited for an expanding background witha small growing perturbation. The scaling of differences like N m − ∞ , ( N + 1) m − N m and N m +1 − N m with N and m areall of interest. We expect the same scaling of these differences with N and m for any uniformly refined set of time steps.A finite order Lagrangian expansion accurate to order m is a truncated representation of the full Lagrangian solution b ( t ) = m (cid:88) i =0 b ( i ) ( t )∆ i . (E6)At the beginning of the first step the scale factor at t is advanced to t and written as b ( t → t ; ∆ , θ ). Note the explicitdependence on the perturbation parameters at t . Abbreviate the scale factor and its derivative for the truncated expressionas b and ˙ b . The background scale factor at time t is a . At the end of the first step the Lagrangian coordinate X and thenew b are inferred as described in the main body of the text by re-scaling quantities calculated at t . The new b is not b .The net result for the full step t → t is X = X ba (E7) b = a (E8)˙ b = ˙ bb a (E9) δ = (1 + δ ) (cid:16) a b (cid:17) − δ v = a ˙ b ˙ a b − . (E11)The newly defined quantities subscripted by “1” will be used to initiate the next step. The updated perturbations imply newLagrangian expansion parameter and angle according to∆ = (cid:112) δ + δ v (E12)cos θ = δ ∆ (E13)sin θ = δ v ∆ . (E14)The new physical position is r = b X = bX . In a numerical calculation the truncated b is exact to floating point precisionbut contains an error because of the omitted orders; in a symbolic calculation b is known to order ∆ m .The next step from t → t involves a similar update with b = b ( t → t ; ∆ , θ ) X = X ba (E15) b = a (E16)˙ b = ˙ bb a (E17) δ = (1 + δ ) (cid:16) a b (cid:17) − δ v = a ˙ b ˙ a b − . (E19)The new physical position is r = b X . This iterative scheme repeats for a total of N steps. It ultimately yields an approxi-mation to the position at the final time denoted N m = b N X N .A difference like ( N + 1) m − N m may be calculated numerically for various N and m and the scaling fitted and inferred.In addition, one can approach the problem symbolically. To write N m we need to expand the final result in powers of ∆ .Note, for example, that ∆ and θ are known as expansions in ∆ with coefficients that depend upon θ . Perturbation-relatedquantities are re-written systematically in terms of initial quantities. For example, b ( t → t ; ∆ , θ ) may be expanded inpowers of ∆ with coefficients depending upon θ . Next, all occurrences of ∆ and θ are replaced by expansions in powersof ∆ and coefficients depending upon θ . All terms up to and including ∆ m are retained in the final result. This procedureis systematically repeated until all quantities are expressed in terms of initial data. Finally, the difference ( N + 1) m − N m iscalculated symbolically. Similar strategies allow construction of all the differences of interest.To make analytic progress assume that f t = β − t f − t ) /t << N + 1) m − N m many “lower order” terms will coincide. Consider an ordering of terms by the powers of ∆ (first) and bypowers of f t (second). Define the leading-order difference to be the first non-vanishing term proportional to ∆ p f qt for smallest c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff p and then smallest q . It is straightforward to apply this ordering to simplify the differences like ( N + 1) m − N m . The leadingorder differences satisfy the following simple equalities | ∞ − N m | ∼ | g N,m || N m +1 − N m | ∼ | g N,m | (E20) | ( N + 1) m − N m | ∼ | g N +1 ,m − g N,m | where g N,m = K N,m cos θ sin m θ ∆ m +1 f m +2 t K N,m = 19 (cid:16) − (cid:17) m (cid:18) N − m m N m +1 (cid:19) . These differences can be compared with the numerical differences for which no expansion in f t is carried out.We verified the analytical scaling by the following numerical experiment. The parameters of the problem at the startingtime t are ∆ = 1 / θ = − π/
4. The final time of interest t f is close to the initial time so that ( t f − t ) /t = 1 /
4. The m -thorder Lagrangian approximation is evaluated at this fixed final time with successively increasing number of steps. Values of m from 1 to 4 and values of N from 10 to 50 were considered. For geometric time steps ( t i +1 − t i ) /t i = β /N − i and denoted δt/t below.The results are plotted in figure E1. The points indicate the numerical data points and the solid lines indicate theanalytical functions defined in eq. (E20). The numerical calculation was done with a high enough precision that even smallerrors of the order of 10 − are not contaminated by floating point errors. The agreement between the numerical experimentand the symbolic differences is very good.Thus, the scaling of the errors implies that for a small total time step, any finite order Lagrangian scheme will yieldconvergent results upon taking multiple steps. Conversely, for a fixed number of steps, a higher order Lagrangian calculationwill give better results.It is useful to express the scaling in terms of the individual small step size δt/t . Under the assumptions that ( t f − t ) /t is small, ( t f − t ) /t ∼ Nδt/t . The scaling | ∞ − N m | ∼ N − m ∆ m +1 (( t f − t ) /t ) m +2 can be re-written as | ∞ − N m | ∼ N (( t f − t ) /t ) · ∆ m +1 ( δt/t ) m +1 , which can be interpreted as an error of (( t f − t ) /t )∆ m +1 ( δt/t ) m +1 per step. In the text,the quantity (cid:15) = ∆ δt/t is kept constant. For fixed initial and final times, the error scales ∝ N(cid:15) m +1 . If ∆ does not changeappreciably then the error is ∝ ∆ (cid:15) m . Convergence is attained when (cid:15) → c (cid:13) RAS, MNRAS000
4. The m -thorder Lagrangian approximation is evaluated at this fixed final time with successively increasing number of steps. Values of m from 1 to 4 and values of N from 10 to 50 were considered. For geometric time steps ( t i +1 − t i ) /t i = β /N − i and denoted δt/t below.The results are plotted in figure E1. The points indicate the numerical data points and the solid lines indicate theanalytical functions defined in eq. (E20). The numerical calculation was done with a high enough precision that even smallerrors of the order of 10 − are not contaminated by floating point errors. The agreement between the numerical experimentand the symbolic differences is very good.Thus, the scaling of the errors implies that for a small total time step, any finite order Lagrangian scheme will yieldconvergent results upon taking multiple steps. Conversely, for a fixed number of steps, a higher order Lagrangian calculationwill give better results.It is useful to express the scaling in terms of the individual small step size δt/t . Under the assumptions that ( t f − t ) /t is small, ( t f − t ) /t ∼ Nδt/t . The scaling | ∞ − N m | ∼ N − m ∆ m +1 (( t f − t ) /t ) m +2 can be re-written as | ∞ − N m | ∼ N (( t f − t ) /t ) · ∆ m +1 ( δt/t ) m +1 , which can be interpreted as an error of (( t f − t ) /t )∆ m +1 ( δt/t ) m +1 per step. In the text,the quantity (cid:15) = ∆ δt/t is kept constant. For fixed initial and final times, the error scales ∝ N(cid:15) m +1 . If ∆ does not changeappreciably then the error is ∝ ∆ (cid:15) m . Convergence is attained when (cid:15) → c (cid:13) RAS, MNRAS000 , 1–38 xtending the domain of validity of the Lagrangian approximation Θ (cid:61) (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) Θ (cid:61) Π (cid:45) (cid:45) (cid:45) (cid:72) H t (cid:76) L og (cid:200) (cid:68) (cid:200) (cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230) (cid:230)(cid:230) (cid:230)(cid:230) (cid:230)(cid:230) (cid:230)(cid:230) (cid:230)(cid:230) (cid:230)(cid:230) (cid:230)(cid:230) (cid:230)(cid:230)(cid:230)(cid:230) (cid:45) (cid:45) (cid:45) (cid:45) (cid:64) (cid:68) (cid:68) Im (cid:64) (cid:68) (cid:68) Figure D3.
Roots with η = 2 π plotted in the complex ∆ plane for 0 < θ (cid:54) π . These values of θ correspond to those in figure D1. Thecolour codes the complex phase of the roots ( ∆ = ∆ e iφ ). The real positive ( φ = 0) and negative ( φ = π ) roots are shown in red and cyanrespectively. The complex roots can have any colour other than these two and the bottom figure provides the coding. By comparison withfigure D1 one sees that all open models with real roots are cyan (negative); likewise all closed models with real roots are red (positive).Note, however, that there exist complex roots for both open and closed models. Coloured version online.c (cid:13) RAS, MNRAS , 1–38 Sharvari Nadkarni-Ghosh and David F. Chernoff (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) N L og (cid:200) (cid:165) (cid:45) N m (cid:200) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) N L og (cid:200) N m (cid:43) (cid:45) N m (cid:200) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) N L og (cid:200) N (cid:43) m (cid:45) N m (cid:200) Figure E1.
The three panels show the log of the errors | ∞ − N m | , | N m +1 − N m | and | ( N + 1) m − N m | vs. N . The final time t f isthe same for all these comparisons. The dots correspond to the data generated by the numerical experiment and lines correspond to theanalytical formulas given in eq. (E20). The lines from top to bottom correspond to m = 1 , , , m = 1 , , m , increasing the number of steps improves convergence. Conversely,for a fixed N , increasing the Lagrangian order m improves convergence. c (cid:13) RAS, MNRAS000