---

NUMERICAL EVOLUTION OF GENERAL RELATIVISTIC VOIDSa

arXiv:hep-ph/9303295v1 23 Mar 1993FNAL–PUB–93/005-AMarch, 1993NUMERICAL EVOLUTION OF GENERAL RELATIVISTIC VOIDSaSharon L. VadasPhysics Dept. and Enrico Fermi Institute, The University of Chicago, Chicago, IL 60637and NASA/Fermilab Astrophysics CenterFermi National Accelerator Laboratory, Batavia, IL 60510AbstractIn this paper, we study the evolution of a relativistic, superhorizon-sizedvoid embedded in a Friedmann-Robertson-Walker universe.

We numericallysolve the spherically symmetric general relativistic equations in comoving,synchronous coordinates. Initially, the fluid inside the void is taken to behomogeneous and nonexpanding.

In a radiation-dominated universe, we findthat radiation diffuses into the void at approximately the speed of light as astrong shock—the void collapses. We also find the surprising result that thecosmic collapse time (the 1st-crossing time) is much smaller than previouslythought, because it depends not only on the radius of the void, but also onthe ratio of the temperature inside the void to that outside.If the ratioof the initial void radius to the outside Hubble radius is less than the ratioof the outside temperature to that inside, then the collapse occurs in lessthan the outside Hubble time.Thus, superhorizon-sized relativistic voidsmay thermalize and homogenize relatively quickly.

These new simulationsrevise the current picture of superhorizon-sized void evolution after first-orderinflation.Submitted to Phys. Rev.

DaPresented as a thesis to the Department of Physics, The University of Chicago, in partial fulfillmentof the requirements for the Ph.D. degree

I. IntroductionThe Big Bang model predicts the evolution of a homogeneous and isotropic universe.The confirmation of its predictions (e.g. the 3K Planck spectrum of the microwave back-ground, the primordial abundances of 4He, 3He, D, and 7Li from Big Bang Nucleosyn-thesis, and the number of light neutrino species) is stunning.

One is led to question whythe universe can be described so well by a homogeneous and isotropic model.Inflation can provide the answer to this question. It occurs when a scalar field σ hasnon-zero potential energy V (σ) which dominates the energy density of the universe.1,2 Thekey ingredients to all models are that the universe expands superluminally during inflationand that there is massive entropy generation afterward.

If the universe increases at least1027 times its original size, the flatness and horizon problems are solved; the universe isdynamically driven to homogeneity and isotropy. The most interesting class of models isfirst-order inflation.

Here, the scalar field is trapped in the false vacuum state of a stronglyfirst-order potential. Bubbles of true vacuum are nucleated at different spacetime points,and the end of inflation occurs when the universe is filled with true-vacuum bubbles ofvarying sizes.

The original model, Guth’s “old inflation”,1 does not work because it failsto percolate (fill) the universe with true-vacuum bubbles. More recent models of first-order inflation which modify the gravitational or particle sector (e.g.

“extended inflation”3) are promising as early-universe inflationary scenarios. Because the universe expandsas a power-law in time rather than exponentially, percolation is guaranteed to occur.The end of inflation occurs when true-vacuum bubbles of different sizes fill all ofspace.

A confusing mess of scalar field dynamics then occurs as the bubbles collide.4Because all the energy is contained in the bubble walls, reheat occurs when the σ-fieldgradient energy is converted into locally thermal radiation. After reheat, the standardhomogeneous and isotropic Big Bang model describes the evolution of that part of theuniverse that had contained horizon-sized bubbles, since these would have been ther-malized during reheat.

The very large superhorizon-sized bubbles (which were nucleatedearly on during inflation), however, have traditionally been a problem.5 After the σ-fieldin the bubble wall decays to relativistic particles, a nearly empty void is formed. Sincethese voids are much larger than the Hubble radius outside the void, it has been thoughtthat the inside of the largest superhorizon-sized voids would remain empty until afterrecombination, thus producing unobserved temperature fluctuations of order unity inthe microwave background.6,7,8 The longevity of these voids follows from assuming thata superhorizon-sized void expands conformally with spacetime.6,8 The earliest time atwhich thermalization could occur would then be when the Hubble radius outside thevoid is of order the size of the void, since this is the expected 1st-crossing time (i.e.

thetime for photons originally in the void wall to reach the origin). If the void has comovingsize r0, this occurs in time ∆t = t−ti ≃H−1out(ti)(c−1Rwall(ti)/H−1out(ti))2, where the initialvoid radius is Rwall(ti) ≡r0a(ti), cH−1out(ti) is the Hubble radius outside the void, a(t) isthe cosmic scale factor and ti is the cosmic time after reheating.

Other authors suggestthat the void would fill in with radiation, but this is estimated to occur on similar timescales.9 Because of this “big bubble problem”, first-order inflation models are fine-tunedin order to keep the production of large bubbles at a minimum. 6,8,9Motivated by the first-order (e.g.extended) inflation “big-bubble problem”, wepresent numerical studies of the evolution of a superhorizon-sized general relativistic voidembedded in a Friedmann-Robertson-Walker universe.

Although pressureless and thin-shell superhorizon-sized voids have been studied in the past,10,11 this is the first study ofgeneral relativistic voids with pressure and of arbitrary size and void wall structure. Weemphasize that our results are not dependent on any particular inflation model.

Thesenew simulations show that opposite sides of a superhorizon-sized relativistic void interact1

in a very short cosmic time, thereby suggesting that these voids can also thermalize andhomogenize on short time scales.In Section II, we discuss the general relativistic spherically symmetric metric in La-grangian gauge and synchronous coordinates, and present the equations to be solvednumerically. In Section III, the initial conditions, boundary conditions and numericaltechniques used are described.

In addition, remarks about the deceleration of a void wallare made. In Section IV, we derive the surprising result that the 1st-crossing time ofa superhorizon-sized general relativistic void can be vanishingly short.

Section V con-tains numerical solutions for pressureless dust and comparisons to exact Tolman-Bondisolutions. A test of the code for the relativistic Friedmann-Robertson-Walker solution ispresented in Section VI.

The numerical evolution of nonrelativistic, special relativistic,and general relativistic voids with pressure is examined in Section VII. Here it is foundthat a relativistic superhorizon-sized void collapses in the form of a strong shock movingat the speed of light.

Thus, the collapse time is approximately the 1st-crossing time.Finally, Section VIII contains a discussion of the results.II. Spherically Symmetric General Relativistic FluidsA.

The (Lagrangian) MetricThe most general spherically symmetric metric is12ds2 = c2l(t, r)dt2 + a(t, r)drdt + h(t, r)dr2 + k(t, r)dΩ2,(2.1)where dΩ2 = dθ2 + sin2 θdψ2 and c is the speed of light. To this metric we can applygeneral transformations of the type t = f1(t′, r′) and r = f2(t′, r′) without altering thespherical symmetry.

If we perform the necessary transformations to eliminate the drdtterm, then Eq. (2.1) can be written asds2 = −c2Φ2(t, r)dt2 + Λ2(t, r)dr2 + R2(t, r)dΩ2,(2.2)where we have dropped the primes, and where we require our coordinates to be comovingwith the fluid.

At time t, the metric function R(t, r) is the Eulerian distance that a fluidshell labeled by r is located from the center of coordinates. More precisely, 2πR(t, r) isthe spacelike circumference of a sphere centered on the origin which contains all particleswith comoving coordinate r. We have thus chosen the Lagrangian gauge with synchronouscoordinates (Gaussian normal coordinates).

Transformations of the form et = f(t) ander = g(r) can still be made. This metric was first used to study the general relativisticcollapse of stars to black holes or neutron stars during supernovae.13,14It is very important in numerical general relativity to carefully choose the appropri-ate gauge and coordinates to best match the physical problem to be studied, We havethus specifically chosen the Lagrangian gauge and synchronous coordinates to evolve ageneral relativistic void embedded in an expanding Friedmann-Robertson-Walker (FRW)universe.

In the Lagrangian gauge, we gain maximal coverage of the fluid in a numericalscheme. This is important since the void is embedded in an expanding universe, so thatwe would continuously lose mass shells in a Eulerian scheme.

In addition, the final re-sults are much more easily relatable to our own approximately FRW homogeneous andisotropic universe. We note that asynchronous coordinates could instead be used withthe Lagrangian gauge.

However, here time and space are mixed up so that the regionoutside a void is no longer spatially homogeneous—the necessary initial conditions arenot obvious and would need to be determined using comoving synchronous coordinatesto initially set up and evolve the void. Asynchronous slicing of space-time (e.g.

polarslicing15) is usually used numerically to study the collapse of a star to a black hole ina flat non-expanding universe. These coordinates are necessary to study the physically2

interesting mass zones outside the apparent horizon after a mass shell has crossed thishorizon. (In synchronous coordinates, once a mass shell crosses the apparent horizon,numerical integration stops so that the evolution of the mass shells outside this horizoncannot be studied).

Since we are evolving an underdense region here and do not expectapparent horizons to form, we do not need to resort to these coordinates.Because we wish to embed a void in a FRW universe, we now relate the metricfunctions from Eqn(2.2) to those from the familiar FRW metric:ds2 = −c2dt2 + a(t)2 dr21 −kr2/c2 + r2dΩ2!,(2.3)where a(t) is the cosmic scale factor, and k is −1, 0 or 1 for negative, zero and positivespatial curvature, respectively. The solutions describe a universe which is homogeneousand isotropic on each time slice.2 For the FRW metric then, Φ = 1, Λ = a(t)/√1 −kc−2r2and the Eulerian distance is R(t) = ra(t).In this paper, the particles are assumed to be everywhere in local thermal equilibrium,so that we can describe them as a fluid.

The stress-energy tensor for a viscous fluid withenergy density ρ and pressure p (measured in the frame of the fluid) is bT αβ = ρuαuβ + pP αβ −2ησαβ −3ζΘP αβ(2.4)where uα is the fluid 4−velocity, P αβ = uαuβ + gαβ is the projection operator, σαβ =1/2 (∇µuαP µβ+∇µuβP µα)−1/3 ΘP αβ is the shear viscosity tensor, Θ = ∇µuµ is the fluidexpansion coefficient, and η and ζ are arbitrary functions of r and t. If η = ζ = 0, the fluidis non-viscous. For Eqn(2.2) with comoving fluid 4-velocity uα = (−Φ−1, 0, 0, 0), σrr = βand σθθ = σψψ = −1/2 β, where β ≡2/3 ( ˙Λ/Λ −U/R), and Θ = ˙Λ/Λ + 2U/R.

(Forthe FRW metric, β = 0 and Θ = 3˙a/a). The stress tensor is diagonal, with componentsTtt = −ρ, Trr = (p−ζΘ)−2ηβ, and Tθθ = Tψψ = (p−ζΘ)+ηβ.

We will employ a scalarartificial viscosity, Q, so that η = 0 and ζΘ = −Q, where Q will be significantly non-zeroonly in areas of steep “velocity” gradients, as in a shock. This is essential for stabilizingnumerical shocks, as is well known.16 This viscosity will dissipate enough energy on smallscales so that the numerical solution approaches the exact solution in the limit that thegrid spacing approaches zero.B.

Fluids composed of massless particlesIf a fluid consists of (effectively) massless, locally thermalized particles, we can relate pto ρ through the equation of state p = p(ρ). For photons (or particles with mass µ havinglocal temperatures T ≫µ), p = ρ/3.

2,17 Setting Rαβ −gαβR = 8πGNTαβ, where Rαβis the Ricci tensor and GN is Newton’s constant, and using the conservation equations∇µT µα = 0, five independent equations are found: G00 = T00, G11 = T11, G01 =T01, T 0µ;µ = 0, and T 1µ;µ = 0. We defineΓ≡R′/Λ,(2.5)M′≡4πc−2ρR2R′(2.6)and U ≡Φ−1(∂R/∂t), where the prime denotes differentiation with respect to r. Thegeneral relativistic equations can be written˙R=ΦU(2.7)bWe set c = 1 for the rest of this subsection3

, ˙U=−Φ GNMR2+ 4πGN(p + Q)Rc2!−c2Γ2Φ(p + Q)′(ρ + p + Q)R′(2.8)˙M=−4π(p + Q)R2ΦU/c2(2.9)˙ρ=−Φ(ρ + p + Q)(R2U)′R2R′(2.10)Φ′=−Φ (p + Q)′ρ + p + Q(2.11)Γ2≡1 + (U/c)2 −2GNM/(Rc2),(2.12)where the dot denotes differentiation with respect to t, and where we have included analternate definition for Γ. There is also the auxiliary equation: ˙Λ = ΛΦU′/R.The quantity U ≡˙R/Φ describes a particle’s “velocity” as measured in the frame ofthe fluid, since for dr = dθ = dψ = 0, the infinitesimal proper time that each observermeasures is dτ =√−ds2 = cΦdt.For GN = 0, if a particle has velocity v, thenΓ = 1/q1 −(v/c)2 and U = Γv (see Eqn (C.9)); Γ and U represent the two non-trivialcomponents of the 4-velocity of the fluid.

If Γ ≫1, then U/c ≃Γ and the fluid is movingat relativistic velocities relative to a stationary observer.We can rewrite Eqn (2.6) in terms of the proper volume element d3V = 4πR2Λdr =4πR2R′dr/Γ. The “mass-energy” function M then becomesM(t, r) = c−2Zd3V Γρ.

(2.13)For GN = 0, because the volume along the radial direction is Lorentz-contracted by thefactor1/q1 −(v/c)2 = Γ and ρ is the energy density measured in the frame of thefluid, Γρ is just the energy density of a fluid parcel as measured by a stationary observer.Therefore, M(r, t) is just the total “mass-energy” contained within comoving coordinater at time t.To better understand these equations, we relate them to the FRW equations. Recallthe FRW results following Eqn (2.3): Φ = 1, R = ra and Λ = a/q1 −kr2/c2.

FromEqn(2.7), the “velocity” is U =˙R = r ˙a.Since R′ = a(t), Γ =q1 −kr2/c2 fromEqn(2.5). For a spatially flat (k = 0) FRW model then, Γ = 1 even though the fluid atR is moving away from the origin with velocity U.

Thus in general relativity, Γ in notthe relativistic gamma-factor of the fluid with respect to the origin. Using the fact thatM = 4πc−2ρR3/3, we see that Eqn (2.12) is just Friedmann’s equation, H2 = (˙a/a)2 =8πGNc−2ρ/3 −kc2/a2, with H ≡Φ−1 ˙a/a = U/R.In a spatially flat FRW universe then, U = c−1Rq8πGNρ/3 =q2GNM/R.

If weset U ≡UGRAV + UPEC, where UGRAV ≡2GNM/R is the gravitational velocity andUPEC is the peculiar velocity, Γ can be expressed in terms of these velocities: Γ2 =1 + (UPEC/c)2 + 2UPECUGRAV/c2.For a non-viscous fluid with p = (γ −1)ρ, Eqn(2.11) (the conservation of momen-tum equation, T µ1;µ = 0) can be integrated exactly to give Φ(t, r2) = Φ(t, r1)[ρ(t, r1)/ρ(t, r2)](γ−1)/γ. We are interested in evolving a void embedded in a FRW homogeneousand isotropic universe, so that p′ = ρ′ = 0 outside the void.

We define Φout and ρoutto be the spatially constant values of Φ and ρ outside the void. (In what follows, the4

subscripts “in” and “out” represent the spatially constant values inside and outside theevolving void region, respectively). Taking γ = 4/3, we findΦ = Φout ρoutρ!1/4.

(2.14)The fact that this equation can be solved exactly is important for calculating the 1st-crossing time for general relativistic voids, as will be discussed in Section IV.To gain some physical insight into Φ, we calculate the potential energy of a fluiddistribution. A comoving observer has dr = dθ = dψ = 0, so that the proper accelerationmeasured by this observer is ar = Φ−1 ˙U.

Since the force is radial, we can write ar =−∂φ/∂R, where φ is the potential. Using Eqn(2.8) and integrating, we can write thepotential as the following sum: φ(t, R) = φGRAV(t, R) + φFLUID(t, R), whereφGRAV=4πGNZ R0hM/R2 + c−2pRidR(2.15)φFLUID=Z R0 c2Γ2dp / (ρ + p)(2.16)for Q = 0.For a void with p = ρ/3, we obtain the usual gravitational potentialφGRAV ≃4πGNρinR2/3 = GNM/R inside the void, and φGRAV ≃4πGNρoutR2/3 ≃GNM/R outside.If in addition we choose Γ(ti, R) = 1 initially, then φFLUID(ti) =ln[ ρ(ti, R)/ρin(ti) ]1/4 = −ln[ Φ(ti, R)/Φin(ti) ].

(Note that the contribution to the fluidpotential is zero inside, but (potentially much) greater than zero outside the void. Thusit can substantially increase the already large potential outside the void).

Therefore, Φis proportional to the exponential of the fluid potential, φFLUID.C. Fluids composed of massive particlesSuppose instead we consider a fluid which consists of particles of mass µ and witharbitrary temperature T. Then, the total energy density of the fluid is the mass energydensity plus the internal energy density.

Denoting the (proper) mass density (µ timesthe number density) by n(t, r) and the internal energy per unit mass by ǫ(t, r),ρ = c2n(1 + ǫ/c2). (2.17)Here we have traded one unknown for two because in general the pressure depends notonly on ρ but also on n. For a fluid composed of relativistic (nonrelativistic) particles,ǫ/c2 > 1 (ǫ/c2 < 1).

If the fluid obeys the perfect gas law, then its pressure isp = (γ −1)nǫ. (2.18)For a highly relativistic species with γ = 4/3, the energy density is three times thepressure: ρ = nǫ = p/(γ −1) = 3p, whereas for a highly nonrelativistic fluid, the energydensity is much larger than the pressure: ρ ≃nc2 = p/[(γ −1)ǫ/c2] ≫p.

We assumethat the total number of particles per comoving volume is constant: ∇µ(nuµ) = 0.19Using Eqn(2.5), this can be integrated to givef(r)=4πnR2R′/Γ(2.19)≡r2(2.20)where f(r) is an arbitrary function depending only on the coordinate r. Specifying f(r)completely fixes the arbitrariness of the metric functions under transformations in r (as5

discussed after Eqn (2.2)). This particular definition for f is necessary in order to writethe difference schemes in a geometrical way that allows shocks and explosions to benumerically stable at the origin.16We can now rewrite the full set of general relativistic equations ((2.7)-(2.12)) as13,14˙R=ΦU(2.21)˙U=−Φ GNMR2+ 4πGN(p + Q)Rc2!−4πΓΦR2(p + Q)′wr2(2.22)˙M=−4π(p + Q)R2ΦU/c2(2.23)˙n=−nΦ(R2U)′R2R′(2.24)˙ǫ=−4πΦ(p + Q)(R2U)′Γ r2(2.25)Φ′=−Φ(p + Q)′nwc2,(2.26)where Γ is given by Eqn(2.12) and w ≡1 + (ǫ + p/n)/c2 is the relativistic enthalpy.Equations (2.21)-(2.26) (along with the definitions for Γ and w given in the previoussentence) are the set used in the numerical code.cWhen the kinetic energy of each particle is much less than its mass energy ǫ/c2 ≪1(or T/µ ≪1), we obtain the nonrelativistic Lagrangian fluid equations.21 (They can alsobe obtained by setting c2 →∞in Eqns (2.21)-(2.26)).

For future reference, in this limitΦ →1, Γ →1, and w →1 so that U = ˙R is the fluid velocity, M(t, r) = r3/3 is the totalmass within r and n = r2/(R2R′) is the mass density.The artificial viscosity used here is given by Equation (C.21), and is generalized fromthe expression used by previous workers:22,14Q=k2n(1 + ǫ/(Γc2))(U′)2dr2/ Γfor U′ < 0Q=0otherwise. (2.27)Consider the behavior of n and p in the limit that the entropy, S(t, r), within r isconserved.

Using the thermodynamic relation TdS = d(ρ/n) + pd(1/n), we find ˙ǫ =p ˙n/n2. Then, the pressure and internal energy for a shell labeled by r are related to theirinitial values by p ∝nγ and ǫ ∝p(γ−1)/γ.

In the ultrarelativistic limit with γ = 4/3,n ∝ρ3/4 and ǫ ∝ρ1/4. In addition, if the fluid is a ploytrope (i.e.

isentropic), S′ = 0.Using TdS = d(ρ/n) + pd(1/n) again, we find the familiar result20ǫ(t, r1)ǫ(t, r2) = n(t, r1)n(t, r2)!γ−1= p(t, r1)p(t, r2)!(γ−1)/γ. (2.28)For ultrarelativistic fluids, ρ = nǫ and therefore ǫ ∝ρ(γ−1)/γ and n ∝ρ1/γ.

If γ = 4/3,ǫ ∝ρ1/4 and n ∝ρ3/4. This is the property of a relativistic fluid, since then ρ ∝T 4 andn ∝T 3.cFor technical difficulties in the general relativistic case, we determine M via the˙M equation ratherthan the M ′ equation.

In addition, we determine n via the ˙n equation rather than through the analyticalsolution n = Γr2/(R2R′), because too much intrinsic viscosity is introduced for special relativistic voidsotherwise.6

D. Fluid Deceleration in the Void WallIn this section we show that if the outward peculiar velocity of the wall of a generalrelativistic void is very large, then the deceleration of this wall can be enormous. This“damping force” is responsible for slowing down and collapsing a superhorizon-sizedvoid.

Without this fluid force, the void would expand (not collapse) from gravitationalforces(see V).In section II-B, we saw that for special relativistic fluids (GN = 0), U is the netoutward momentum per particle mass µ.Similarly, if in the general relativistic caseU2 ≫U2GRAV = 2GNM/R, the gravitational attraction inward is negligible and thedynamics will be dominated by special relativistic effects; U can be loosely interpretedas the peculiar momentum per particle per mass. If Γ ≫1 and U > 0 in the wall ofa void, the wall moves outward with momentum much greater than the gravitationalattraction inward.

We now investigate what happens to this wall. From Eqn (2.22), the“conservation of momentum” equation isnwΦ−1 ˙U=−nw GNMR2+ 4πGN(p + Q)Rc2!−Γ2(p + Q)′R′.

(2.29)The functions Γ, Φ, n, p, Q, M, w and R′ are always positive. If U > 0, then at theinner edge of the void wall where (p + Q)′ > 0, ˙U will be negative—the fluid there isdecelerated.

This deceleration is due to both gravitational and fluid forces. We examinethe fluid force contribution only.

Using Eqns (2.12), (2.6) and (2.22) we find that˙Γ = −ΓUΦ(p + Q)′nwc2R′ . (2.30)Again, if (p + Q)′ > 0 and U > 0, ˙Γ will be negative and the fluid particles there losetheir energy per particle mass.Suppose the wall has a very large outward momentum so that U/c ≫q2GNM/(Rc2) >1.

Then Γ ≃U/c. Setting p = ρ/3 and Q = 0 as for a relativistic, non-viscous fluid,˙Γ = −Γ2 Φρ′4ρR′ ≃−Γ2ρ1/4out ρ′4 ρ5/4 R′ ,(2.31)where we have used the solution for Φ from Eqn (2.14).

We consider the deceleration offluid shells at the inner edge of a steep void wall. We can write ρ′/R′ ≃∆ρwall/∆Rwall,where ∆ρwall is the difference in the energy density over the width of the wall, and ∆Rwallis the thickness of the wall.

Then, since ∆ρwall ≃ρmax, where ρmax is the maximum wallenergy density, we can approximate ˙Γ by˙Γ ≃−Γ24 ρoutρ4maxρ5!1/4∆Rwall−1. (2.32)Initially, except for the first factor of Γ2, all factors on the right hand side of the previousequation are independent of Γ.

Thus, ˙Γ ∝−Γ2, which can be a very large dampingfactor! The second factor is proportional to ρoutρ4max/ρ5 ≤(ρout/ρ)5, so for the massshells on the innermost part of the wall (i.e.

those shells with the smallest values ofρ/ρout) this factor can be enormous. The third factor is ∆Rwall−1, so the thinner thewall, the faster it will slow down.

If the second and third factors change slowly enough7

with time, then the slow-down time for the wall is roughly independent of its initial valueof Γ0, since −R 1Γ0 dΓ/Γ2 ≃1 >∼(ρout/ρ)5∆Rwall−1∆t/4 for Γ0 ≫1. This could be anextremely important result, and would imply that the initial peculiar wall velocity couldnever be large enough to cause a void to expand for an arbitrarily long time.

However,because the second and third factors in Eqn (2.32) will change with time and depend onΓ, this is only a crude guess.As a concluding remark, we note that the wall of a void formed during first-orderinflation has an enormous outward peculiar velocity. This enormous velocity has beenthought to cause a void to expand “indefinitely”.

However, with such a large decelerationof the void wall, it will slow down in a finite (and possibly small) amount of time. Afuture paper will explore this numerically.34III.

Initial Conditions, Boundary Conditions, and Numerical TechniquesA. Initial ConditionsThe grid used in this code is initially equally spaced: ∆R(ti) ≡R(ti)j+1 −R(ti)j =constant, where the subscript j denotes the spatial grid point number and ranges fromj ∈[0, jB], where jB is its value at the outer boundary.

(Note from Eqns (2.19) and(2.20) that in general, ∆r ̸= constant then). The 0th and 1st-grid points are located atR0 = −∆R(ti)/2 and R1 = ∆R(ti)/2 respectively.

Thus, Rj(ti) = R1(ti)+(j −1)∆R(ti).The initial conditions for the functions at ti are determined as follows. The viscosityis set to zero: Q(ti, R) = 0.

The energy density ρ(ti, R) (or the “mass-energy” M(ti, R)),the “velocity” U(ti, R) (or Γ(ti, R)) and the specific internal energy, ǫ(ti, R), are chosenas functions of the radius R(ti). (If M is specified initially instead of ρ, we determineρ via Eqn.

(2.6)). For the cases run in this paper however, we will initially specify thefluid to be a polytrope (constant entropy on the initial time slice (II-C)) and ǫ/c2 ≫1or ǫ/c2 ≪1.

Using the relations found at the end of Section II-C, ǫ(ti, R) is determinedonce the internal energy ǫ(ti, RB) is specified at the outer boundary (RB ≡R(ti)jB):ǫ(ti, R)=ǫ(ti, RB) [ ρ(ti, R)/ρ(ti, RB) ] (γ−1)/γfor ǫ/c2 ≫1ǫ(ti, R)=ǫ(ti, RB) [ ρ(ti, R)/ρ(ti, RB) ] (γ−1)for ǫ/c2 ≪1(3.1)We then determine n and p by n(ti, R) = ρ(ti, R)/(1 + ǫ(ti, R)) and p(ti, R) = (γ −1)n(ti, R)ǫ(ti, R) respectively. Next, we find r (and M if it is not initially specified) byintegrating outward from r = 0 using the 4th-order Runge-Kutta method:r= 3Z R(ti,r)04πnR2dR/Γ!1/3 andM(ti, R) =Z R0 n(1 + ǫ/c2)R2dR!,(3.2)where we have omitted the (ti, R)’s for clarity.Finally, Φ is found by integratingEqn (2.26) inwards from the outer boundary (again using the 4th-order Runge-Kuttamethod) once Φ(ti, RB) is specified.We are interested in evolving a non-expanding, empty void embedded in a FRWuniverse.We will therefore initially set the energy density outside the void (ρout ≡ρ(ti, RB)) to be constant and equal to the spatially flat FRW value.

Then, H2out(ti) ≡H2(ti, RB) = (ξ/ti)2 = 8πGNc−2ρout(ti)/3, where aout(t) = a(ti)(t/ti)ξ is the cosmic scalefactor and cH−1out is the Hubble radius outside the void. (Note that for the k = 0 FRWuniverse, r = (4πn(ti))1/3R(ti, r), or a(ti) = (4πn(ti))−1/3).

For p = ρ/3 and p = 0,8

ξ = 1/2 and ξ = 2/3 respectively. We can quantify the initial size of a void by measuringits radius relative to the horizon size outside the void initially: c−1Rwall(ti)/H−1out(ti) =c−1ξRwall(ti)/ti, where Rwall(ti) is the void “radius”.Following past convention2, weloosely equate the Hubble radius with the horizon in the phrase “superhorizon-sized”.

(Horizon in this context is not to be confused with the particle horizon.) A void is definedto be superhorizon-sized if c−1Rwall(ti)/H−1out(ti) > 1 and subhorizon-sized if c−1Rwall(ti)/H−1out(ti) < 1.

In addition, Γout(ti) = 1 (see II-B) (or U =q2GNM/R ) and Φout(ti) ≡1(see II-A).In this paper, we consider voids which are initially either compensated or uncompen-sated in energy density and which have the following distributions. The “mass-energy”function for compensated voids is defined to beM(ti, R) = .5ρout(ti) [(1 + tanh x) + α(1 −tanh x)] R3(ti)/3,(3.3)where x ≡(R(ti) −Rwall(ti))/∆Rwall(ti), ∆Rwall(ti) is the wall thickness and α is a spec-ified constant less than or equal to 1.

Because M reaches its spatially flat FRW valueoutside the void, the energy density missing from the void has been put in the wall. Foruncompensated voids, the energy density is instead initially specified.

It isρ(ti, R) = .5ρout(ti)[(1 + tanh x) + α(1 −tanh x)]. (3.4)Here, the energy density missing from the void has not been put into the wall.

Therefore,the region outside a compensated void will always be a spatially flat (k = 0) FRWuniverse, whereas the region outside an uncompensated void is a negative spatially curved(k < 0) FRW universe.The inside of the void is chosen to be homogeneous initially. Since we want it to benonexpanding in the limit that it is empty (ρ →0), we choose the inside of the void to be aspatially flat (k = 0) FRW “mini” universe.

We reason as follows. Friedmann’s equationinside the void is H2in ≃8πGNc−2ρin/3 −kc−2r2/R2 (see II-B), where the subscript ‘in’denotes quantities inside the void, R = rain and Hin ≡Uin/R = Φ−1in ˙ain/ain.

In the limitthat ρ →0, the inside of the void is non-expanding (˙ain = 0) only if k = 0. Since wecannot numerically choose ρin = 0 (because Φin = ∞from Eqn(2.14)), we would like theinside of the void to not expand on time scales that the outside region expands in.

Thisis satisfied if ρin/ρout ≪1 because the Hubble time inside the void is much larger thanthat outside the void— H−1in /H−1out =qρout/ρin. (As a check, Eqns.

(2.22) and (2.25)show that as U →0, ǫ →0 and p →0, then ˙U →0 and ˙ǫ →0). Because the void isapproximately homogeneous, the “mass-energy” inside the void is M(ti, R) ≃c−2ρinR3/3.Finally, we set Γin(ti) = 1 inside the void since Γ =q1 −kr2/c2.

The “velocity is thenU =q2GNM/R ≃c−1Rq8πGNρ/3.We will only consider two types of initial velocity profiles in this paper. The first isU = UGRAV =q2GNM/R (or Γ(ti, R) = 1).

This specifies that the outward “velocity”of each particle is just large enough to compensate for the inward gravitational attraction(i.e. the peculiar velocity, UPEC, is zero).

The second isU(ti, R) = c−1Rq8πGNρ/3. (3.5)For R ≪Rwall(ti) −∆Rwall(ti) and R ≫Rwall(ti) + ∆Rwall(ti), Γ ≃1.

In the wall region,however, Γ > 1. This corresponds to an initial net outward “peculiar momentum” of theparticles in the wall.9

B. Boundary ConditionsWe set the outer boundary conditions for all times to be those for a homogeneousfluid.This specification works well in practice as long as the action is taking placeaway from this boundary. Thus we set n′ = ǫ′ = p′ = 0 and Q = 0 at j = jB, wherejB is the grid point number for the outermost comoving coordinate.

In addition, weset ΦjB(t) ≡1, its FRW value (see Eqn (2.3). Note that specifying ΦjB(t) completelyeliminates the arbitrariness of the time coordinate, as discussed after Eqn (2.2).

Thepresent definition sets t to be the FRW cosmic time outside the void. Thus if an initiallyinhomogeneous fluid becomes homogeneous, then from Eqn(2.26), Φ(t, r) = 1 everywhere,and t = constant hypersurfaces correspond to t = constant FRW homogeneous andisotropic hypersurfaces.It is possible to determine the outer boundary conditions for all t by solving theequations with Q = p′ = 0.

One is then left with two 1st−order ordinary differentialequations to solve. The boundary conditions determined this way, however, give largererrors than the ones shown below, and therefore were not used.

We instead integrate˙R=U˙U=−GNM/R2 + 4πGNpR/c2˙M=−4π pR2U/c2. (3.6)using the MacCormack method.

We then determine n, ǫ and p from Eqns (2.28) andEqn (2.18) by setting njB = njB−1, ǫjB = ǫjB(ti) [ njB(t)/njB(ti) ]γ−1 and pjB = (γ −1)njBǫjB.Finally, reflecting boundary conditions are used at the inner boundary: Rn0 = −Rn1,Un0 = −Un1 , pn0 = pn1, nn0 = nn1, ǫn0 = ǫn1 and Qn0 = Qn1, where the superscript n refers tovalues on the nth time slice.C. Numerical IntegrationThe ˙U,˙ǫ,˙R, ˙n and˙M equations are integrated using the 2-step MacCormackpredictor-corrector method.23 In Appendix B, we give the exact form for the differenceequations which allows inbound shocks to rebound offthe origin.

To illustrate the Mac-Cormack method, we show the predictor and corrector steps for ˙n as an example. Wecontinue to use the convention that nij is the value of n on the jth spatial grid point andat the ith time step.

Suppose we know all quantities on the ith time slice. We would liketo determine them on the (i + 1)th time slice.

First we predict the new quantities (withforward differencing) using the functional values on the ith time slice:npi+1j= nij −∆t nijΦijRij+12 U ij+1 −Rij2 U ijRij2(R ij+1 −R ij ). (3.7)After using similar methods to obtain Upi+1j , Rpi+1j , Mpi+1jand ǫpi+1j(and setting ppi+1j=(γ −1)npi+1jǫpi+1j) for all j, we integrate the Φ′ equation inwards from j = jB using the4th-order Runge-Kutta method with linear interpolations to determine Φpi+1jfor all j. Weintegrate again, using the predicted values obtained above (with backward differencing),and then average these with the previously predicted values:ni+1j= .5npi+1j + nij −∆t npijΦpijRpij2 Upij −Rpij−12 Upij−1Rpij2(Rpij −Rpij−1).

(3.8)10

We obtain the other values in a similar way, and then integrate again to find Φi+1j .As is well known,16 it is important to choose small enough time steps ∆t to satisfythe Courant condition. This condition requires ∆t to be smaller than the time taken forsound to cross from any one grid point to the next.

The speed of sound for relativisticfluids is cS =q(∂p/∂ρ)S.13,14 Using the fact that TdS = d(ρ/n) + pd(1/n), we findcS =r γpnw. (3.9)As c2 →∞, w = 1 and we obtain the usual expression for the speed of sound innonrelativistic fluids.

Note that for ǫ/c2 ≫1, cs = c√γ −1 = .57c for perfect fluids withγ = 4/3. The Courant condition requires the proper time per proper distance to be equalto C/cS, where C is the Courant number and is approximately 1/3 to 1/2 for strongshocks.

Since the infinitesimal proper time and proper distance are Φ∆t and Λ∆r =∆R/Γ respectively, the Courant condition becomes ∆tni = C (Rni+1 −Rni )/(Γni Φni (cS)ni ).14 In addition, when GN ̸= 0, regions of the universe expand and contract. We requirethe time steps to be small enough for the functions n, ǫ, and M to change sufficientlyslowly.We therefore set ∆tni = f (lni /˙lni ) to be the maximum time step allowed forl = n, ǫ and M, where f is a constant less than one.After the nth corrector step, ∆tni is calculated for all i, and ∆tn+1 is set to the smallestvalue obtained:(∆t)n+1 = min∆tmax,C Rni+1 −RniΓni Φni (cS)ni,"f nni˙nni,f ǫni˙ǫni,f Mni˙Mni#if GN̸=0,(3.10)where ∆tmax is a specified upper bound, if desired.25We apply a convergence test to the code for test problems where analytic solutionsare available.

The relative error in q, where q denotes any quantity, is defined to beei = |qi −eq(ri)|/eq(ri),(3.11)where eq(ri) is the exact solution and qi is the numerical solution. We obtain a globalmeasure of the error by definingL1 = 1NNXi=1ei,(3.12)where N is the total number of grid points.

This error is proportional to the grid spacingto some power: L1 ∝∆Rs, where s is the convergence rate.If s ≃2, the code issecond-order, as desired. These tools have been used previously to test codes in otherapplications.26IV.

1st-Crossing Time for Relativistic VoidsIn this section, we first review the standard lore for the evolution of superhorizon-sized voids, and then calculate the 1st-crossing time (the cosmic time taken for a photoninitially at the inner edge of the void wall to reach the origin) in this picture.Wethen calculate the actual 1st-crossing time.If at the 1st-crossing time the fluid wereapproximately homogeneous and isotropic, distortions from the original void would benegligible and a (nearly) FRW universe would result everywhere.11

It has been suggested that superhorizon-sized voids formed from first-order infla-tion would conformally expand with spacetime during the radiation-dominated period.6,8Thus, the size of a void at time t would be R = r0a(t), where r0 is the comoving coor-dinate of the void and a(t) is the cosmic scale factor. There are several justifications tosupport this belief.

First, small density perturbations conformally expand with space-time. Second, vacuum bubbles conformally expand during inflation after an initial shortgrowing period.

Third, the time taken for the void wall to slow down is expected to beenormous because the initial outward momentum of the void wall is enormous. However,this reasoning is not enough to conclude that a superhorizon-sized void conformally ex-pands with spacetime.

First, a void is not a small perturbation in spacetime. Thus linearresults can not be applied to the description of a void.

Second, it is the negative pressurethat causes vacuum bubbles to expand; a similar configuration having positive pressurewould instead acclerate inward.d Third, although part of the wall may still move out,the deceleration of the inner void wall is extremely large (see II-D), so that the void canstill collapse (see VII).If a superhorizon-sized void were to conformally expand in spacetime, then the earliesttime at which thermalization and homogenization can occur is when the horizon is of orderthe size of the void, since this is the expected 1st-crossing time. If the void has comovingsize r0 and the outside Hubble radius (that outside the void) is cH−1out(ti), this occurs whenr0a(t) = H−1out(t).

For evolution during the radiation-dominated epoch, H−1out(t) = 2t andp = ρ/3 so that the time is of order ∆tc ≡t −ti ≃H−1out(ti)(c−1Rwall(ti)/H−1out(ti))2. If thevoid is much larger than the Hubble radius outside the void (c−1Rwall(ti)/ H−1out(ti) ≫1),then the 1st-crossing time is very large: ∆tc/H−1out(ti) ≫1, independent of the “emptiness”of the void.

Other authors suggest that the void would fill in with radiation.9 If spacetimeat the void wall continues to expand as a(t) when the fluid diffuses into the void, thecomoving radius of the wall is roughly r ≃r0 −cR tti dt/a(t) orra(ti) ≃Rwall(ti) −cZ ttidt′qti/t′(4.1)in a radiation-dominated universe.The 1st-crossing time (i.e.the earliest thermal-ization time) then, is when r ≃0, or when ∆tc ≃.5H−1out(ti)(c−1Rwall(ti)/H−1out(ti))2,which is roughly the same as the time for the horizon to “engulf” a conformally ex-panding void.eBecause R = ra(t), the radius of the void would be given by R ≃(Rwall(ti)−cH−1out(ti)qt/ti )qt/ti, which for nearly all of the time is R ≃Rwall(ti)qt/ti; thevoid conformally stretches with spacetime. At time t ≃H−1out(ti)(c−1Rwall(ti)/H−1out(ti))2,the void radius decreases quite rapidly to zero.

Thus, although the qualitative void evo-lution is quite different in these two pictures, the quantitative 1st-crossing times are notbecause in both the void comoves with spacetime. The important point is that the ear-liest possible thermalization time in both pictures (i.e.

the 1st-crossing time) is thoughtto be∆tc ≃H−1out(ti)(c−1Rwall(ti)/H−1out(ti))2. (4.2)dIn the bubble wall, GN = 0, |p|′ > 0 and p < 0.

1 Using Eqn(2.22), we see that the acceleration ispositive, so that the bubble wall moves outward during inflation. But for normal positive pressure withp′ > 0, the acceleration is negative, so that the wall must accelerate inward.eWe thank Michael Turner for this explanation.12

We will show in this section that the actual 1st-crossing time is remarkably shorter thanEqn(4.2). In doing so we will show that Eqn(4.1) is fundamentally flawed.Consider two radially propagating photons A and B. Photon A starts at the inneredge of the void wall and propagates inward, and photon B begins at the outer edge ofthe void wall and moves outward.

We would like to calculate the distance each travels intime ∆t = t−ti, where ti is the initial time. Using Eqns (2.2) and (2.5), the infinitesimalcoordinate distance traveled by a photon in time dt is dr = cdtΦΓ/R′.

Define ρin(t)and ρout(t) to be the energy densities inside and outside the evolving wall region of thevoid. (Thus the subscripts “in” and “out” refer only to the undisturbed fluid).

Initially,ρ′in(ti) = ρ′out(ti) = 0, and Γ(ti, R) = 1. We also consider non-viscous, relativistic fluids,so that p = ρ/3 and Q(t, r) = 0.We will first calculate the distances photons A and B travel in special relativisticvoids (GN = 0).

From Eqn. (2.22), the fluid acceleration is zero inside and outside thevoid:˙Uin = ˙Uout = 0.

Therefore, R(t, r) = R(ti, r) inside and outside the void, so thatthe infinitesimal distance traveled by photons A and B in time dt is dR = cΦΓdt. FromEqn.

(2.10), we see that ρout and ρin are constant in time. Since Φout = 1, we find thatΦin(t, r) = Φin(ti, r) = (ρout(ti)/ρin(ti))1/4 from Eqn.

(2.14). In addition, since ˙Γ ∝p′(Eqn.

(2.30)), Γin(t) = Γout(t) = 1. Therefore in time ∆t ≡t −ti, photon B travelsoutward the distance ∆RB(t) ≡RB(t) −Rwall(ti) given by∆RB = cΦoutΓout∆t = c∆t,(4.3)while photon A travels inward the distance ∆RA(t) ≡Rwall(ti) −RA(t) given by∆RA(t) = cΦinΓin∆t = c [ρout(ti)/ρin(ti)]1/4 ∆t = c Tout(ti)/Tin(ti) ∆t > c∆t.

(4.4)Thus the emptier the void, the farther photon A moves relative to photon B! (It isimportant to note that this is a strictly relativistic effect; for nonrelativistic voids (II-C),Φ(t, r) ≃1 inside and outside this void so that the distance traveled by photons A andB are approximately the same: ∆RA ≃c∆t = ∆RB).

Define ∆tc to be the 1st-crossingtime (the time taken for photon A to reach the origin).Then from Eqn. (4.4) with∆RA = Rwall(ti), the 1st-crossing time is∆tc = c−1Rwall(ti) [ ρin(ti)/ρout(ti) ]1/4 .

(4.5)Since ρout(ti) ≥ρin(ti), ∆tc ranges from Rwall(ti)/c to zero. Therefore in the limit thatρin(ti)/ρout(ti) →0, a photon will reach the origin in zero cosmic time.

However, calling aregion a “fluid” if it is empty (ρin = 0) is incorrect. Suffice it to say that the 1st-crossingtime can be arbitrarily small.We note that if Eulerian, synchronous coordinates were used with the metric ds2 =−c2dT 2 + dR2 + R2dΩ2, the distance traveled by photon A or B in time ∆T would bethe same: ∆R = c∆T.

Thus the time ∆Tc for photon A to reach the origin is ∆Tc =c−1Rwall(ti). (We can also see this using Lagrangian coordinates, since the infinitesimalproper time measured by a comoving observer inside the void is dτ = Φindt = cdT,and therefore the time as measured by this observer for photon A to reach the originis ∆τc = cΦin∆tc = Rwall(ti)).

Thus, the quick 1st-crossing time is due to the choiceof comoving, synchronous coordinates, a choice we do not have for superhorizon-sizedgeneral relativistic voids (GN ̸= 0) embedded in a FRW universe.ff If the spatially flat FRW metric is transformed to the Eulerian gauge with synchronous coordinates,a coordinate singularity at the Hubble radius (see Appendix A). Thus this gauge and coordinate choicecannot be used to describe the evolution of superhorizon-sized general relativistic voids.13

We now calculate the distance traveled by photons A and B in the same cosmic timefor a general relativistic void.Again, the infinitesimal coordinate distance a photontravels in time dt is dr = cΦΓdt/R′. Because the pressure outside the void redshifts dueto Hubble expansion, Φin(t) decreases in time from Eqn(2.14) (since Φout(t) = 1), and∆tc consequently increases.

We first calculate the distance photon B travels. Outsidethe void, Γout(t) = 1, and R = ra(t) so that R′ = a(ti)qt/ti.

24 The comoving distancephoton B has traveled at time t is ∆r = cti/a(ti)(qt/ti −1), so that the distance photonB travels is ∆RB ≡RB −Rwall(ti) = (qt/ti −1)[Rwall(ti) + 1/2 cH−1out(ti)qt/ti].We now calculate the location of photon A. Assume that the energy density insidethe void changes negligibly: ρin(t) = constant.

(We will address this approximation in amoment). Then Φin(t) = (ρout(t)/ρin(ti))1/4.

Because the energy density outside the voidis redshifted as ρ ∝1/t2(see VI),24 Φin(t) ≃Φin(ti)qti/t. In addition, since ˙Γ ∝p′ (fromEqn.

(2.30)), Γin(t) = Γout(t) = 1. Integrating, we find that the cosmic time taken forphoton A to travel the distance ∆RA ≡Rwall(ti) −RA is∆t ≡t −ti = c−1∆RA ρin(ti)ρout(ti)!1/4 1 + c−1∆RA2H−1out(ti) ρin(ti)ρout(ti)!1/4,(4.6)and the location of photon A as a function of time isRA=Rwall(ti) −cΦin(ti)Z ttidt′qti/t′(4.7)=Rwall(ti) −cΦin(ti)H−1out(ti)qt/ti −1.

(4.8)As in the special relativistic case, ∆RA(t) > ∆RB(t), so that photon A travels fartherthan photon B in the same amount of cosmic time. The 1st-crossing time relative to theinitial outside Hubble time is then (∆RA = Rwall(ti))∆tcH−1out(ti) = c−1Rwall(ti)H−1out(ti) ρin(ti)ρout(ti)!1/4 1 + c−1Rwall(ti)2H−1out(ti) ρin(ti)ρout(ti)!1/4.

(4.9)In addition, if the more stringent condition c−1Rwall(ti)/H−1out(ti) < Φin(ti) holds, then∆tcH−1out(ti)<∼1whenc−1Rwall(ti)H−1out(ti)< ρout(ti)ρin(ti)!1/4(4.10)—the minimum thermalization time (i.e. the 1st-crossing time) is less than the initialHubble time outside the void!

(Note that the general and special relativistic results(Eqns(4.9) and (4.5)) are equal in this case, since the energy density outside the void isconstant during this time).Before discussing further implications, we find the condition for which Eqn(4.6) is sat-isfied; it was derived under the assumption that ρin(t) ≃ρin(ti). Using Eqn(2.10), insidethe void | ˙ρ/ρ| = 4(R2 ˙R)′/(3R2R′) = 4 ˙R/R (since R = ra(t)) so that if ρin ≃constant,then R ≃constant (in time) inside the void.

Because U2 = ( ˙R/Φ)2 = 2GNM/R ≃14

c−2(8πGNρin/3)R2 inside the void, the fractional change in the radius R over time scale∆t is roughlyf ≡∆RR =˙R∆tR=vuut ρin(t)ρout(t)Hout(t)Φin(t)∆t ≃Φ−1in (ti)stit∆tH−1out(ti),(4.11)where we have used the fact that H−1out(ti) = 2ti,qρin(t)/ρout(t)Hout(t) ≃Φ−2in (ti)Hout(ti)and Φin(ti) = (ρout(ti)/ρin(ti))1/4. We require f < 1.

Writing t = ∆t + ti, Eqn(4.11)becomes ∆t = f 2Φ2in(ti)H−1out(ti)[1 +q1 + 1/(fΦin(ti))2].Since Φin(ti) > 1, we havefΦin(ti) >∼1. The condition for which the density inside of the void remains approxi-mately constant during time ∆t then, is ∆t/H−1out(ti) ≤2f 2Φ2in(ti) ≃Φ2in(ti).

We now com-bine this with Eqn(4.6) to find the maximum allowed void size ∆RA/H−1out(ti) given Φin(ti).We find c−1∆RA/H−1out(ti) ≤.5Φin(ti)[−1 +q1 + 4(fΦin(ti))2] ≃fΦ2in(ti). Eqn (4.6) isthen valid when c−1∆RA/H−1out(ti) ≤fΦ2in(ti), and Eqn.

(4.9) is valid whenc−1Rwall(ti) / H−1out(ti) <∼qρout(ti) / ρin(ti). (4.12)Note that Eqn(4.10) is automatically satisfied.The quick 1st-crossing time might seem completely counterintuitive.

How can a pho-ton travel a distance much larger than the Hubble radius in less than a Hubble time? Theanswer lies in describing how one measures the size of an object which is not a small per-turbation in spacetime.

If size is measured circumferentially, then the void is enormousbecause its circumferential size is 2πRwall(ti) ≫H−1out(ti). (If spacetime were static, thenit would take a photon time ∆t ≃2πRwall(ti)/c to encircle the void).

However, if sizeis measured radially (the time taken for a photon to cross the object if spacetime werestatic), then the void is measured to be very small. If fact, an interesting comparisoncan be made to measuring the size of a black hole, an overdense region.

If one measuresits circumferential size, it is small (or at least finite), but its radial size is infinite.The 1st-crossing time for a superhorizon-sized void was previous calculated incorrectlybecause t was assumed to be the proper time outside and inside the void; in our notation,it was implicitly assumed that Φ(t, r) = 1. Comparing Eqns(4.1) and (4.7), we indeed seethat the factor Φin(ti) > 1 is missing from Eqn(4.1).

Note in addition that the factor ofqti/t in Eqn(4.7) does not come from spacetime expanding at the wall, but rather fromthe density outside the void redshifting, causing Φin(t) to decrease. In fact, spacetimeis roughly non-expanding at the inner wall edge.

If c−1Rwall(ti)/H−1out(ti) > Φin(ti) ≫1,the actual position of photon A is approximately constant in time until the last moment:RA ≃Rwall(ti) until t ≃ti + ∆tc (see Eqn(4.8)). Thus, spacetime at the inner edge of thevoid wall is neither expanding nor contracting.We note the interesting fact that if a photon inside the void is in thermal equilibriumwith average frequency ν, its frequency is Tin.

If it moves outside the void, its frequency isblue-shifted to ν′ = νΦin = ν(Tout/Tin) = Tout, the average frequency of thermal photonsoutside the void. Thus, this photon is automatically in thermal equilibrium outside thevoid, so that an outside observer could not detect the void’s initial presence unless non-thermal photons came out.g This is essentially because the fluid is relativistic, or that ρ,gBecause entropy is created at the shock (see VII), the photon with νTin would actuall have a slightlyhigher frequency than thermal outside the void.15

p, n and ǫ depend only on the temperature (II-C). We can understand why the photonis blue-shifted by way of comparison to an (overdense) black hole.

Suppose an observerfalls into a black hole ticking offphotons at a fixed frequency. As this observer crossesthe event horizon, the frequency of the photon emitted last gets redshifted to infinity asobserved by a stationary observer outside the black hole.12 The opposite effect happensfor a photon emitted from an underdense region.

Because a photon leaving a void enters aregion with a much larger gravitational potential, the frequency instead gets blue-shifted.In conclusion, the 1st-crossing time for a superhorizon-sized relativistic void embed-ded in a FRW expanding universe is given by Eqn (4.9) (if Eqn(4.12) is satisfied), anddepends sensitively on the quantity Rwall(ti)/H−1out(ti) (ρin(ti)/ρout(ti))1/4. We emphasizethe important point that if c−1Rwall(ti)/H−1out(ti) < (ρout(ti)/ρin(ti))1/4 = Tout(ti)/Tin(ti),then the 1st-crossing time is less than the outside Hubble time: ∆tc/H−1out(ti) < 1.V.

Pressureless Non-Viscous VoidsFor the special case when the pressure and viscosity are zero, the equations of motioncan be solved analytically. These give the Tolman-Bondi dust solutions,10,12 which wewill review here briefly.

It is important to study the pressureless case not only as a testproblem, but also to see how removing fluid forces affects the evolution of a void. (Becausep = Q = 0, the particles move only under gravitational forces: Φ−1 ˙U = −GNM/R2 (seeEqn (2.22))).Since Φ′ = 0 (Eqn (2.26)), we set Φ(t, r) ≡1.

The mass contained within r re-mains constant: M(r) =R r0 4πnR2dR/Γ, since˙M = 0 (see Eqn (2.9)).And fromEqn(2.30), Γ(t, r) = Γ(ti, r). For R′ ̸= 0, ρ = n can be found from Eqn(2.19): ρ(t, r) =ρ(ti, r) R(ti, r)′R(ti, r)2/[R(t, r)′R(t, r)2].

The generalized pressureless Friedmann equa-tion, Eq. (2.12), now becomes˙R2 = c2 hΓ(r)2 −1i+ 2GNM(r)/R.

(5.1)The quantity U2/2 −GNM/R = c2(Γ2 −1)/2 is conserved during evolution and can beinterpreted as the generalized total energy.Eqn (5.1) can be integrated for a shell ofradius r:R=c−2GNMΓ2 −1(cosh η −1) ,t = τ0(r) + c−3GNM(Γ2 −1)3/2 (sinh η −η) ,for Γ(r)2 > 1R=c−2GNM1 −Γ2(1 −cos η) ,t = τ0(r) + c−3GNM(1 −Γ2)3/2 (η −sin η) ,for Γ(r)2 < 1R=(9GNM/2)1/3 (t −τ0(r))2/3 ,for Γ(r) = 1. (5.2)These are the Tolman-Bondi solutions.In these models, shell-crossing can occur.

This happens when two adjacent shells(labeled by r and r +dr) occupy the same position so that R′ = 0 and ρ →∞. This maylead to a non-unique continuation of the solution, 27 and thus computations have to bestopped.

This problem is believed to occur because the pressure has been artificially setto zero. It is generally thought that adding pressure would prevent this situation fromoccurring.h As will be seen, when Γ(r) > 1 in the void wall region, shell-crossing doesoccur.

The addition of enough artificial viscosity can prevent this from happening, how-ever. Even though the viscosity given by Eqn (2.27) was designed to stabilize numericalshocks, it has been found to prevent shell-crossing.hThe author thanks T. Piran for this information.16

Table 1: Convergence test for Tolman-Bondi model∆R(ti)L1sL1asL1bs81.2 × 10−3...8.33 × 10−4...8.55 × 10−4...44.13 × 10−41.482.26 × 10−41.892.16 × 10−42.021.55 × 10−41.326.09 × 10−51.825.35 × 10−51.9716.56 × 10−51.241.81 × 10−51.751.41 × 10−51.92For the numerical simulations in this section, we set GN = 1, c = 1, ti = 1, ρ(ti, R) =0, ǫ(ti, RB) = 0, C = .3, f = .005 and γ = 5/3. Thus, the initial energy density andHubble radius outside the void are 4πρout(ti) = 2/3 and H−1out(ti) = 3/2, respectively (seeIII-A).We first examine the situation in which each mass shell’s velocity initially compensatesfor the gravitational attraction inwards: Γ(ti, R) = 1.

Then U(ti, R) =q2GNM/R.Figure 1 shows the energy density versus R R(ti, RB)/R(t, RB) (≡R RjB(ti)/RjB(t)) fora superhorizon-sized compensated void with c−1Rwall(ti)/H−1out(ti) = 333, ∆Rwall(ti) = 15,α = .001, ∆R(ti) = 2.5 and k2 = 0. We show the analytic and numerical results attimes t = 1, 10, 100, and 300 where RjB(t) = 1002, 4651, 2.16 × 104, and 4.49 × 104respectively.

The triangles and squares are the numerical and Tolman-Bondi solutions,respectively, although they are difficult to distinguish because the numerical results agreeso well with the analytic results. By t ∼300, the density everywhere is approximatelyconstant; there is hardly a trace of the void’s initial presence.Identical results areobtained for subhorizon-sized voids.

In addition, an initially uncompensated void evolvessimilarly.In Table 1, we show the results of a convergence test for q = ρ applied to the sameinitial conditions as in Figure 1, but for variable ∆R(ti) (see III-C). (We only do thistest for ρ, because the accumulated error in M and R are much smaller).

We set upanalytic conditions initially and integrate until t = 1.5. Because the inner grid pointor two ends up being the numerical culprit for non-second order convergence, we alsocalculate L1a ≡1N−2PNi=3 ei and L1b ≡1N−10PNi=10 ei, where ei is the relative error.

(Thisis because the code consistently underestimated ρ at the innermost few grid points). Tothe right of each global error estimate, the convergence rate is shown (Li ∝∆R(ti)s).We see that the convergence rate for L1 is less than second-order, whereas that for L1 isnearly second order.We have just seen that if U = UGRAV (Γ(ti, R) = 1), the void disappears.

Whathappens when U > UGRAV (Γ(ti, R) > 1) in the wall region? Recall that this correspondsto a net outward peculiar velocity (II-D).

Figure 2 shows the result for a compensatedvoid with c−1Rwall(ti)/H−1out(ti) = 333 , ∆Rwall(ti) = 15, α = .001, ∆R(ti) = 2.5, andk2 = 0. We choose the initial velocity to be given by Eqn(3.5).

Therefore, Γ(ti, R) = 1everywhere except in the wall region, where Γ(ti, R) > 1.Again, the triangles andsquares represent the numerical and analytic solutions, respectively.The initial time(ti) and t = 1.02, 1.048 are shown. Again, the difference between the numerical andanalytical results are small except at t = 1.048, where shell crossing occurs in bothsolutions, a good check on the code.

The comoving radius of the shell with the highest17

density, rshell(t), remains approximately constant in time. Because each shell in the wallhas constant total energy c2(Γ2 −1)/2, the wall expands outward.

(Since U > 0, Rincreases. But since the total energy is constant and M(r) is constant, U must increase).Identical results are obtained for subhorizon voids.Figures 3a and 3b show the density as a function of R R(ti, RB)/R(t, RB) for asubhorizon-sized, shallow, uncompensated and compensated void, respectively, at t =1, 50, 100, 400 and 1000.

Note that shell-crossing would have occurred at t = 116 and 2.1for the uncompensated and compensated voids, respectively. Here, c−1Rwall(ti)/H−1out(ti) =.0067, ∆Rwall(ti) = .001, ∆R(ti) = .0001, α = .5, and the velocity U(ti, R) is given byEqn (3.5).

In addition, for Figure 3a, k2 = 4 and RjB(t) = .035, .49, .79, 2.13 and 4.26,while for Figure 3b, k2 = 8 and RjB(t) = .035, .48, .76, 1.9 and 3.5. As can be seen, theinitial perturbation grows with time and eventually forms a thin, dense shell.

Again, thecomoving coordinate for the shell with the highest density is approximately constant intime after the shell has formed. As the shell travels outward, it pushes mass in front ofit, producing a shock.

This situation is similar to that of a fast car colliding with slowerones; although the faster car is not allowed to move through the slower cars, momentumis still transferred to them. Note that the initially compensated perturbation forms athick shell more quickly than the uncompensated perturbation, although it then proceedsto grow more slowly.As long as the voids formed remain subhorizon-sized, they will eventually grow ac-cording to a known similarity solution.28 An initially compensated (uncompensated) per-turbation in an expanding FRW p = 0 universe will eventually form a dense, thin shellthat expands outward as Rshell ≡R(t, rshell) ∝t4/5 (t8/9).

In Figure 3c, we show theposition of the void wall versus time for the results of Figures 3a and 3b. The trian-gles and squares connected by lines are the numerical solutions for the compensated anduncompensated cases, respectively, and the dashed and dotted lines are the self-similarsolutions for the compensated and uncompensated cases, respectively.

As can be seen, theinitially compensated perturbation approaches the similarity solution more quickly thanthe initially uncompensated perturbation, but for t >∼800, both solutions are self-similar.VI. FRW Homogeneous CosmologiesWe now test our code against the exact FRW homogeneous and isotropic solutionfor relativistic fluids with p = ρ/3.

As discussed in II-A,B and III-A, the FRW solutionis R(t, r) = ra(t) = R(ti, r)(t/ti)ξ where ξ = 1/2. In addition, 4πGNρ(t) = 3c2/(8t2)and U = R/(2t).For the numerical results shown in this section, we set GN = 1,c = 1, γ = 4/3, C = .3, ti = 1, ǫ(ti, RB) = 106 and c−1RB(ti)/H−1(ti) = 250 so that4πρ(ti, RB) = 3/8.Figure 4 shows the relative error in ρ (III-C) for a simulation with ∆R(ti) = .5 andk2 = 0.

The analytical solution was set up initially, and the code was run until t = 1.1.The solid, dotted and dashed lines are for f = .01, f = .005 and f = .0025. The relativeerror at the outer boundary is seen to be very sensitive to f; it is .05%, .15% and 2% forf = .0025, .005 and .01 respectively.

Note also the underprediction of ρ at the innermostfew grid points.Table 2 shows the accumulated error at time t = 1.1 for simulations with k2 = 4.iAgain, L1b ≡1N−10PNi=10 ei, where ei is the relative error for q = ρ. In the first 7 columns,we show the results for N = 25, 50, 100, 200, 400, 800 and 1600 (i.e.

∆R(ti) = N/250)iWhen evolving voids, viscosity is absolutely necessary. It is therefore important to see how it affectsthe solution where the fluid is approximately homogeneous and isotropic.

It turns out that it is virtuallyunaffected by Q ̸= 0, as it should be.18

Table 2: Ultrarelativistic Homogeneous convergence testNumber of grid points (N)L125501002004008001600f∆R →0sL1.016.010.0073.0058.0050.0046.01.004...L1b.0071.0056.0049.0045.0044.0043L1.012.0067.0038.0023.0016.0012.005.000852.00L1b.0015.0012.0010.00091.00088.00087L1.011.0059.0031.0017.00096.00060.00042.00025.000251.77L1b.00057.00041.00032.00027.00026.00025.00025and for f = .01, .005, and .0025. For a given value of f, as N increases, L1 and L1bapproach the same constant value even though L1 starts out much larger.

After thisconstant value has been reached, L1 and L1b remain unchanged when ∆R(ti) is furtherdecreased, even though the relative error near the origin improves.This is becausebeyond a certain point, all the accumulated error comes from the outer boundary, whichis immune to changes in ∆Rwall(ti) (see III-B).Knowing that a given value of f limits convergence of the code, we can calculate theconvergence as a function of the asymptotic value of L1. We assume that L1(∆R(ti) →0) ∝fs.Column 9 gives the estimated value for L1(∆R(ti) →0), and column 10estimates the value of s. We see that s is nearly 2, which means that convergence in f isnearly second order given a small enough value for ∆R(ti).VII.

Numerical Evolution of VoidsA. Nonrelativistic FluidsIn this subsection we examine the evolution of voids composed of nonrelativistic par-ticles in zero gravity (GN = 0).

If T is the fluid temperature and µ is a fluid particle’smass, T/µ ≪1. Since H−1 ∝GN−1, these voids are subhorizon-sized.

For the simu-lations in this subsection, we set GN = 0, c = 1010, ti = 1, C = .3, γ = 5/3, C = .3,4πρ(ti, RB) = 2/3 and ǫ(ti, RB) = 1.We start with the shock tube problem, a standard test of 1-D slab codes.29 In addition,it provides insight into the dynamics of collapsing voids. In a shock tube, the fluid isinitially at rest and is separated into two regions with different pressures and densities.The pressure discontinuity produces a shock wave which propagates into the low pressureregion and a rarefaction wave which propagates into the high pressure region.

An analyticsimilarity solution exits for a perfect fluid with slab geometry. It does not exist for thespherically symmetric geometry however.30 Far from the origin however, the sphericallysymmetric solution approaches the slab solution for small times and distances.14 We willthus set up a spherically symmetric shock tube by evolving an uncompensated void farfrom the origin, and compare the results to the exact slab similarity solution (brieflyreviewed in Appendix D).19

In Figure 5 we show the results for the shock tube problem with U(ti, R) = 0 (orΓ(ti, R) = 1), Rwall(ti) = 20, ∆Rwall(ti) = .01, ∆R(ti) = .01, α = .01 and k2 = 5. (Wedo not take the pressure gradient to be discontinuous initially, because the solution inthe original wall area is not as accurate then.j) We plot the pressure, number density,velocity and specific energy as a function of the position R at ti and at t = 1.15.

Thetriangles connected by lines is the numerical solution, and the dashed lines is the slabsimilarity solution. At t = 1.15, a strong shock wave (located at R ≃19.68) movesinward and a rarefaction wave (between 19.85 <∼R <∼20.14) moves outward.

Note thatthe shock is spread out over k2 ≃4−5 grid points. The numerical and analytical solutionsare seen to agree well in this limit, because the spherical geometrical effects are small((Rshock(t)−Rwall(ti))/Rwall(ti) ∼.4/20 = .02).

The distortion of the velocity distributionis due to the small geometrical effect; the shock gets slightly stronger directly behind theshock due to the smaller effective volume 4πR2∆R those mass shells occupy relative toshells further back. An important point to emphasize is that although initially the fluidis everywhere stationary (U(ti, R) = 0), it acquires a net momentum to the left in thewall region.In Figure 6a we show the long-term results for an initial configuration with Rwall(ti) =1 but otherwise identical to Figure 5.

The pressure, number density, velocity and specificenergy are shown at the initial time ti = 1 with dashed lines, and long after the collisionat t = 2.5 with triangles and connecting lines.Although it is not shown for clarity,the fluid configuration before the shock rebounds at the origin is similar to Figure 5: ashock heads toward the origin and a rarefaction wave moves away from the origin. Whenthe (spherical) shock crashes into the origin, the volume effect proves harsh as the fluidcollides with itself in a vanishingly small volume.

This causes the pressure at the originto become very large in order to repel the fluid (not pictured here). Note that the fluidat the shock as well as far behind it must be repelled, since it is all moving toward theorigin.k When the dust settles, we find that a weak shock has rebounded back.

Thisoutward-moving shock can be seen at R(t=2.5) ≃1.3. Due to the volume effect however,this shock will become weaker as it moves outward further.

Also seen at t = 2.5 is theoriginal (outgoing) rarefaction wave located between 1.3 <∼R <∼2.4. Note that at t = 2.5,the fluid is (and will approximately remain) at rest near the origin because the velocityis zero and the pressure is constant.

However, a large distortion has been left behind inthe fluid in the form of low density and high kinetic energy. This consists entirely of thefluid originally in the void.

The low-n, high-ǫ values for the first 5-6 grid points, however,is artificial. This is a consequence of using VonNeumann-type artificial viscosity called“wall-heating”, and is caused by the collision between 2 shocks.31In Figure 6b we show the initial and long-term pressure, number density, velocity andspecific energy as a function of the radius for a compensated, nonrelativistic void.

ThejThis is primarily because we calculate n via ˙n ∝(R2dR)−1 rather than ˙n ∝(r2dr)−1kIt is this reversal that requires a very robust difference scheme.All of the “obvious” differenceschemes failed after the inbound shocks reached the origin. Setting f = r2 in Eqn (2.20) and using thosedifference schemes listed in Appendix B are the very necessary requirements.20

Table 3: Nonrelativistic collapse timesVoid type∆Rwall(ti) = .02∆Rwall(ti) = .04∆Rwall(ti) = .08∆Rwall(ti) = .12Collapse Time ∆tc||Percent Change: (∆tc −∆tc(.02)un)/∆tc(.02)ununcompensated.38 ||....41 ||8%.45 ||18%.50 ||32%compensated.15 || −61%.20 || −47%.27 || −29%.32 || −16%initial distribution is identical to Figure 6a except for the energy density, and ∆Rwall(ti) =.02. Again, the fluid configuration at ti = 1 is shown as dashed lines and that at t = 2.5is shown as triangles connected by lines.

An inbound shock is again formed from theinitial pressure gradient at the inner edge of the void wall. Unlike the uncompensatedcase however, a weak outgoing shock wave is formed instead of a rarefaction wave.

Att = 2.5, it is located at R = 3.2. Like the uncompensated void, the pressure at the originbecomes very large after the shock collides with itself there, and a weak shock rebounds.At t = 2.5, this rebounded shock is located at R ≃2.3, which is farther out than that forthe uncompensated void (R ≃1.3).

This is because the inbound shock produced fromthe compensated case is much stronger than for the uncompensated case. And againthere is a large distortion left near the origin containing all the fluid initially in the void,although it is somewhat different spatially.We now compare the collapse times for uncompensated and compensated voids ofvarying wall thicknesses.

(The collapse time is defined to be the time taken for the shockto reach the origin). We set Rwall(ti) = 1, α = .01, k2 = 3 and ∆R(ti) = .01.

Table 3shows the results for ∆Rwall(ti) = .02, .04, .08 and .12, where we calculate the percentchange by comparing the collapse time with the uncompensated ∆Rwall(ti) = .02 collapsetime of ∆t = .38. As ∆Rwall(ti) increases, the collapse time increases.

In addition, thecollapse time for compensated voids is substantially smaller than for uncompensatedvoids.In conclusion, nonrelativistic voids with zero gravity collapse in the form of a shock,the strength of which depends on the details of the void wall. Some time after collapsing,the fluid is virtually at rest everywhere, with n and ǫ inhomogeneous near the origin.B.

Special Relativistic FluidsIn this subsection we consider the evolution of special relativistic voids with T/µ ≫1,where T is the temperature and µ is the mass of a fluid particle. Thus, ǫ/c2 ≫1, and weset GN = 0, c = 1, ti = 1, C = .3, γ = 4/3, ǫ(ti, RB) = 103, U(ti, R) = 0 (or Γ(ti, R) = 1)and 4πρ(ti, RB) = 3/8 in this subsection.In Figure 7a, we show M, Γ, 4πρ and Φ as a function of R for a relativistic shock tubeproblem.

Here Rwall(ti) = 1, k2 = 3, ∆Rwall(ti) = .02, ∆R(ti) = .01 and α = 10−4. As inthe nonrelativistic case, an inbound shock is formed.

(This is observed most easily in theplot of 4πρ versus R). However, an outgoing “rarefaction wave” is not observed, as it is in21

the nonrelativistic case (see Figure 5). If a wave were to propagate outward, it could goat most the distance a photon would travel.

But from Eqn(4.3), we know that a photonstarting from the outer wall edge only moves the distance ∆RB = .04 (.065) in time∆t = .04 (.065). Since this is of order the grid point thickness, if an outbound wave werepresent, it would not be observed at this time anyway.

On the other hand, an inboundphoton starting from the inner wall edge would travel the distance ∆RA = .4 (.65) (fromEqn(4.4)) in time ∆t = .04 (.065), since ρout(ti)/ρin(ti) ≃104. Reexamining Figure 7a,we now notice an important result; the inbound shock’s position is approximately equalto photon A’s location—the shock moves inward at roughly the speed of light.

This isactually not so surprising, because the speed of sound for a perfect fluid with p = ρ/3 is.57c (III-C).Consider next voids compensated in energy density. Because of relativistic-particlediffusion, we expect the inbound shock to again travel at approximately the speed of light.Since the value of Φin(ti) does not depend on the functional form of ρ in the void wall(Eqn (2.14)), the shock should move the same distance per time as for the uncompensatedvoid.

We ran a compensated void simulation with the same initial conditions as fromFigure 7a (except for ρ in the void wall), and found that at t = 1.04 and 1.065, thecompensated shock is ahead by only ∆R = .05. However, this can be accounted for bythe slightly different values of ρ initially at the inner edge of the wall.

Thus, specialrelativistic compensated and uncompensated voids collapse at approximately the samespeed, unlike nonrelativistic voids (see Table 3).In Figure 7b, we show the numerical results for an uncompensated void with Rwall(ti) =1, k2 = 3, α = 10−4 and ∆R(ti) = .01, and for ∆Rwall(ti) = .02, .04, .06 and .1. It isclear from this figure that in all four cases the shock reaches the origin at ∆tc ∼.08;∆tc is approximately independent of ∆Rwall(ti).

(Compare this with Table 3). UsingEqn (4.4), we estimate the time for light to reach the origin at ∆tc ≃Rwall(ti)/Φ ≃.09,in agreement with the simulations.l (During this time, photon B would only move out-ward the distance ∆RB = ∆tc ≃.09).

The value of ∆Rwall(ti) however, does influenceshock formation-time and strength. As ∆Rwall(ti) increases, the shock formation timeincreases and the shock strength decreases.We now examine the what happens to the void after the shock collides at the origin.In Figure 7c, we show the numerical results for an uncompensated void with Rwall(ti) = 1,∆Rwall(ti) = .02, α = 10−4, k2 = 8 and ∆R(ti) = .01.m We show p, n, U and ǫ asa function of R at the initial time ti and at times 1.04 and 3.0, where the second andthird times are before and after collision at the origin, respectively.

At t = 3, p′ ≃0 andU′ ≃0 near the origin; the fluid there is roughly at rest. Two weak outgoing waves arelBecause Φin actually increases slightly due to the dissipation of energy at the shock, the actual traveltime for the shock to reach the origin, ∆tc, is decreased slightly.mAfter colliding, a very weak shock rebounds.The generalized functional form for the artificialviscosity, however, was derived in the strong shock limit.Thus, a large value of k2 was needed tomaintain numerical stability after the collision.

A new functional form for Q will have to be used infuture simulations to stabilize the weak outgoing shock.3422

observed: the shock (at R ∼1) and a “rarefaction wave” (between 1.1 <∼R <∼2.1). Inaddition, in contrast with nonrelativistic voids, n′ ∼0 and ǫ′ ∼0 near the origin.n Thisresult is expected; since n ∝p3/4 ∝T 3 and ǫ ∝p1/4 ∝T (from the discussion followingEqns (2.28)) for a non-viscous fluid, if p′ ≃0, then it follows that n′ ≃0 and ǫ′ ≃0.This is a consequence of the fact that ρ, p, n and ǫ depend only on the temperature forǫ/c2 ≫1 and Q = 0.

We therefore find that after the collapse, there is only a small traceof the void’s initial presence!We conclude that a special relativistic void collapses in the form of a shock whichtravels at approximately the speed of light into the void. Thus, the collapse time is oforder the 1st-crossing time.

Some time after the collapse, the fluid becomes approximatelyhomogeneous and isotropic everywhere.C. General Relativistic FluidsIn this subsection we study the evolution of general relativistic voids for T/µ ≫1,where T is the temperature and µ is the mass of a fluid particle.

We set GN = 1, c = 1,ti = 1, C = .3 and γ = 4/3. Therefore, the outside Hubble radius and density (thatoutside the void) are cH−1out(ti) = 2 and 4πGNρ(ti, RB) = 3/8.In Figure 8, we show the pressure, number density, velocity and specific energyfor a general relativistic void at the initial time ti and for t = 1.04 and t = 1.065.The initial conditions are identical to those of Figure 7a, except that k2 = 4 andU(ti, R) = UGRAV =q2GNM/R (or Γ(ti, R) = 1).Because the void is subhorizon-sized (c−1Rwall(ti)/H−1out(ti) = 1/2), its evolution looks virtually the same as that forthe special relativistic void shown in Figure 7a; at the void wall, the gravitationalforce is M/R2 + 4πpR/c2 ≃8πρR/3 ≃10−1, which is much less than the fluid forceΓ2p′/(4pR′) ≃(4∆R(ti))−1 ≃12 (see Eqn.

(2.8)).In Figure 9a, we show the pressure as a function of R R(ti, RB)/R(t, RB) (≡R(t, r)RjB(ti)/RjB(t)) for a general relativistic superhorizon-sized void at ti, t = 2.0 and t = 8.0with Γ(ti, R) = 1, Rwall(ti) = 50, ∆Rwall(ti) = 1., k2 = 4, α = 10−4, ǫ(ti, RB) = 103 and∆R(ti) = .5. (This void is 50/2 = 25 times the outside Hubble radius).

In addition,RjB(t) = 100.2, 139.1, and 271.8. Even though GN ̸= 0 here, a strong inward shock stillforms.

This is because particles are diffusing into the void, having been accelerated awayfrom the high-pressure wall. (Note that at the wall, the acceleration due to gravity isGNM/R2 + 4πpR/c2 ≃10, whereas that due to the fluid force is only Γ2p′/(4pR′) ≃.25.This small relative amount however, is enough to form the shock).

Note that becauseof expansion, the pressure outside (and to a lessor extent inside) the void redshifts. Byt = 8, the pressure outside the void has redshifted from p ≃10−1 to 2×10−3, the expectedamount since p ∝1/t2 so that pout(8) ≃10−1/64 ≃2 × 10−3.

At the same time, thepressure on the inside has only redshifted from p ≃10−5 to 5 × 10−6, because the insideHubble time is larger than that outside: H−1in (ti)/H−1out(ti) =qρout(ti)/ρin(ti) = 100.From Eqn(4.9), the 1st-crossing time is ∆tc ≃50/10(1 + 25/(2 × 10)) = 11.25, in roughnFor the innermost 7-8 grid points, ǫ and n are too large and too small, respectively. This is againdue to “wall-heating”.23

agreement with the numerical collapse time of ∆tc ≃8.Thus, the collapse time isfound to be approximately equal to the 1st-crossing time—the the shock moves inwardat roughly the speed of light. This is not surprising however, because the speed of soundis .57c (III-C).Figures 9b and 9c show the pressure as a function of R R(ti, RB)/R(t, RB) for ageneral relativistic superhorizon-sized void with c−1Rwall(ti)/H−1out(ti) = 25, Γ(ti, R) = 1,k2 = 4, ∆Rwall(ti) = 1, ǫ(ti, RB) = 106 and ∆R(ti) = .5.In addition, α = 10−6and RjB(t) = 100.2, 125.1, and 148.7 for Figure 9b, and α = 10−10 and RjB(t) =100.2, 102.6, and 104.4 for Figure 9c.

The numerical collapse times in Figures 9b and 9care ∆tc = 1.3 and .09, respectively, which are both smaller than the outside Hubble time. (Note that ∆tc = .09 is 1/20th the outside Hubble time).

Since Φ(ti) equals 31.6 and316, respectively, the 1st-crossing times from Eqn(4.9) are ∆tc ≃50/31.6(1 + 25/31.6) ≃2.2 and 50/316(1 + 25/316) = .17, respectively. The 1st-crossing times are larger thanthe numerical collapse times because Φin(t) increases during the collapse due to thedissipation of energy at the shock.

Thus, the gain in entropy at the shock only makesthe collapse time shorter. For example, in Figure 9c, Φin(t) increases from its originalvalue of Φin = 316 to Φin ∼450 −500 during the collapse.

Using Φin = 500, we wouldpredict ∆tc ∼.1, which is roughly correct.We note however, that in Figure 9c (and to a lesser extent in Figure 9b), only part ofthe void has been filled in at the collapse time tc ≡ti + ∆tc. Thus, thermalization andhomogenization has not been achieved by tc.

We note from Figures 9a,b and c that asthe initial relative energy density inside the void decreases, the fraction of the void filledin at the collapse time decreases. However, since the energy density inside the somewhatfilled void is still much less than ρout(tc), Φ(tc) inside the “void” is still much greaterthan one.

And because Γ(tc) is also larger than one inside the “void”, the distance lightcan travel is still greater inside than outside the void. Thus, although the void has notyet homogenized, this may only take an additional small amount of time.Figures 10a and 10b show the pressure versus R R(ti, RB)/R(t, RB) for a superhorizon-sized uncompensated and compensated void, respectively.

Here, Γ(ti, R) = 1, Rwall(ti) =500, ∆Rwall(ti) = 10., k2 = 4, ǫ(ti, RB) = 106, ∆R(ti) = 5., and α = 10−10. (These voidsare 250 times the outside Hubble radius).

The pressure is shown at ti = 1 and t = 1.7for each void, and at t = 2.1 and t = 2.0 for the uncompensated and compensated voids,respectively. In addition, in Figure 10a, RjB(t) = 1002, 1288, and 1424, while in Figure10b, RjB(t) = 1002, 1307, and 1418.

The important point to note is that the shockreaches the origin at t ≃2.1 for both voids. This is due to the diffusion of particlesinto the void, and does not depend on the compensatedness of the void because Φin(ti)depends only on ρin(ti)/ρout(ti).

(This is similar to that for special relativistic relativisticvoids (VII-B)). This roughly agrees with the predicted collapse time ∆tc = 500/316(1 +250/(2 × 316)) ≃2.2 from Eqn(4.9).

Note also that at t ≃2, the difference betweenFigures 10a and 10b is small. For the compensated void, there is no perceptible outboundshock, and the original density “bump” in the void wall is stretched and damped out.

Infact, the fluid initially in the compensated void’s wall is moving inward at t ≃2.1.24

Up to this point, we have shown the evolution of general relativistic uncompensatedand compensated voids for the initial velocity profile U = UGRAV =q2GNM/R (orΓ(ti, R) = 1), where the initial velocity per shell just balances gravity. If the wall ini-tially has an outward peculiar velocity (as expected from first-order inflation, for exam-ple), then U(ti, Rwall(ti)) > UGRAV(ti, Rwall(ti)) (or Γ(ti, Rwall(ti)) > 1), as discussed inII-D.

In Figure 11, we show 4πp, M, Γ and Φ versus R R(ti, RB)/R(t, RB) for a com-pensated superhorizon-sized void with U/c = Rq8πGNρ/3 (Eqn(3.5)), Rwall(ti) = 500,∆Rwall(ti) = 10., k2 = 4, ǫ(ti, RB) = 106, ∆R(ti) = 2.5, f = .005 and α = 10−10.The voids are superhorizon-sized: Rwall(ti)/H(ti, RB)−1 = 250. In addition, RjB(t) =601.2, 783.9, and 911.8, and we show the configurations at ti and at times t = 1.7 andt = 2.3.

The void wall is initially moving outward with a very large peculiar velocity,since Γ(ti, Rwall(ti)) ≃900.We find the very interesting result that even though the wall has a large outwardpeculiar velocity, the inner part of the void still collapses. This is because the fluid nearthe base of the void wall gets accelerated into the void right away, pulling adjacent fluidwith it.

It is true however, that at t ∼2 the density profile looks different than thatfrom Figure 10b. This is because the fluid in the void wall takes more time to lose itsoutward velocity.

Since the numerical collapse time is ∆tc = 1.3 from Figure 11, we findthat the extra time taken for the void to collapse is approximately 1.3 −1.1 ≃.2, whichis still less than the outside Hubble time. (This follows qualitatively from our discussionin Section II-D, where it was argued that Γ would decrease very quickly at the inneredge of the wall, since ˙Γ ∝−Γ2p′).

These new results show that the collapse time of asuperhorizon-sized void with a large outward peculiar wall velocity can be of order the1st-crossing time. Since the minimum thermalization and homogenization time is the1st-crossing time, the time for thermalization and homogenization of this void may beshort.In Figure 12, we show the void radius versus cosmic time for a void with initialsize Rwall(ti) = c1023H−1out(ti) .

If the temperature outside the void initially is Tout(ti) =1014GeV, then the initial time is ti = 10−33 seconds.2 And because recombination oc-curs at trec ≃1012 seconds, trec/ti = 1045, which is near where the dashed lines in-tersect at the top of Figure 12. Therefore, the plot consists the radiation dominatedepoch in the early universe, where the scale factor outside the void is a(t) = a(ti)qt/ti.Because R(ti, RB)/R(t, RB) =qti/t, any comoving point (defined as R ∝qt/ti, notr = const) is a vertical line in this plot.The dashed line labeled rCM shows thevoid size if it were comoving with spacetime, as previously suggested (see IV).

TheHubble radius outside the void is cH−1out ≡RHOR = 2ct, and is shown as the dashedline labeled by rHOR.As can be seen, rHOR and rCM intersect at t/ti ≃1046 (or10−11 seconds), shortly after recombination. (This also follows from Eqn(4.2), since∆tc ≃H−1out(ti)(c−1Rwall(ti)/H−1out(ti))2 ≃1046).

If a superhorizon-sized void were to con-formally expand with spacetime, then this void would just barely be around to distortthe microwave background at recombination.825

The dash-dot lines show the position of the 1st-crossing photons (i.e. the inner voidwall) from Eqn(4.8).

We show Φin(ti) = 5×1011, 1015 and 1020 (since we require Φin(ti) >∼√1023 (from Eqn(4.12))), with the predicted 1st-crossing times of tc/ti ≃∆tc/ti = 4×1022,1016 and 106 respectively. If ∆tc/H−1out(ti) > 1, the radius of the inner edge of the wallis constant (R ≃const) until time t ≃tc ≃ti(Φ−1in (ti)Rwall(ti)/H−1out(ti))2, at whichpoint R rapidly decreases to zero (see Eqn(4.8)).

If the reheat temperature is TRH =Tout(ti) = 1014GeV, then the void with Φin(ti) = Tout(ti)/Tin(ti) = 1020 may thermalizeby t ≃10−27 seconds (or at T ≃TRHqti/t ≃1011GeV), far before recombination ornucleosynthesis. (Note that this value of Φin(ti) corresponds to an initial void temperatureof Tin(ti) = TRHΦin(ti)−1 = 10−6GeV = 1keV).In conclusion, we have seen that a compensated or uncompensated superhorizon-sizedvoid collapses as fluid in the wall diffuses into the void at the speed of light.

Thus, thecollapse time is found to be approximately equal to the 1st-crossing time calculated inSection IV, which is much smaller than previously thought.VIII. DiscussionIn this paper, the evolution of a superhorizon-sized void embedded in a Friedmann-Robertson-Walker (FRW) universe was studied by numerically evolving the sphericallysymmetric general relativistic equations in the Lagrangian gauge and synchronous co-ordinates.

(A superhorizon-sized void is a void larger than the Hubble radius outsidethe void: c−1Rwall(ti)/ H−1out(ti) > 1, where ti is the initial cosmic time, cH−1out(ti) is theHubble radius outside the void and Rwall(ti) is the radius of the void). The particles areassumed to be in local thermal equilibrium so that they can be described as a fluid, andthe perfect fluid equation of state is chosen.

The inside of the void is initially chosento be homogeneous and non-expanding. We are particularly interested in the evolutionof voids with T/µ ≫1, where T is the local temperature and µ is the mass of a fluidparticle.

This corresponds to a radiation-dominated period in the early universe. Wefind that for p = ρ/3, general relativistic voids collapse via a shock propagating inwardfrom the void wall.

The shock is formed from the steep pressure gradient at the inneredge of the wall. It moves at approximately the speed of light; this is not surprising,since the speed of sound is .57c.

Thus the void collapse time (the time taken for theshock to reach the origin) roughly equals the photon 1st-crossing time (the time takenfor a photon initially at the inner wall edge to reach the origin). At the time this shockreaches the origin, much of the fluid in the original wall area is moving toward the originbehind the shock.

The energy density of this fluid (i.e. the fluid behind the shock) isless than that outside the void.

At the collapse time then, only part of the void has beenfilled in with fluid. In particular, as the initial energy density inside the void decreases,the fraction of the void filled in at the collapse time decreases.Because the shock moves inward at roughly the speed of light, we can calculate theapproximate collapse time.

For p = ρ/3, the 1st-crossing time is∆tc = c−1Rwall(ti) ρin(ti)ρout(ti)!1/4 1 + c−1Rwall(ti)2H−1out(ti) ρin(ti)ρout(ti)!1/4(8.1)26

for Rwall(ti)/H−1out(ti) <∼qρout(ti)/ρin(ti), where ρin(ti) and ρout(ti) are the initial fluid en-ergy densities inside and outside the void, respectively, and ρ = T 4. If c−1Rwall(ti)/H−1out(ti) <(ρout(ti)/ρin(ti))1/4, then ∆tc/H−1out(ti) <∼1; the 1st-crossing time is less than or compara-ble to the initial Hubble time outside the void.

In fact, as the density inside the voidapproaches zero, the collapse time goes to zero!It may seem contradictory that the 1st-crossing time for a superhorizon-sized relativis-tic void can be less than the outside Hubble time (i.e. or that light travels comparativelyfarther inside a void than outside).

There are several points to make about this. First, ifwe choose Eulerian instead of Lagrangian synchronous coordinates to evolve special rela-tivistic voids (GN = 0), then the 1st-crossing time is not fast.

This is because spacetimeis sliced differently in time inside the void. Since a similar 1st-crossing time is obtainedfor general relativistic voids, it might be argued that the quick collapse time is a fig-ment of the Lagrangian gauge used.

However because the FRW metric has a coordinatesingularity at the Hubble radius when transformed to Eulerian synchronous coordinates,superhorizon-sized voids cannot be evolved in these coordinates. Any prior intuition fromspecial relativistic voids in these coordinates therefore, must be carefully applied.Second, calling a void “superhorizon-sized” is deceptive.

Although the 1st-crossingtime for a superhorizon-sized perturbation in the FRW universe is greater than the outsideHubble time, this does not imply that the same is true for a superhorizon-sized void. Thisis because the void is not a small perturbation in general (i.e.

(ρout −ρin)/ρout ≃1). Ifthe void’s size is defined to be the spacelike circumferential radius relative to the outsideHubble radius, then the void is superhorizon-sized if Rwall > cH−1out.However, if thevoids’s size is defined to be its radius relative to the Hubble radius inside the void, thenits size is smaller since c−1Rwall(ti)/H−1in (ti) = c−1Rwall(ti)/H−1out(ti)qρin(ti)/ρout(ti); it issubhorizon-sized if c−1Rwall(ti)/H−1out(ti)

(We call this latter size the“radial size”, because it is the time taken for a photon to cross a static void multipliedby c). This is because the Hubble radius is much larger inside than outside the void.Since the void is an underdense region, we expect the opposite relative size problem tooccur for overdense regions where the gravitational potential is large, not small.

As anexample, a black hole’s circumferential size is small (or at least finite), while its “radialsize” is infinite.The original motivation for studying this problem was to determine the evolution ofvoids formed from first-order (e.g. extended) inflation.

These voids are compensated inenergy density, and the walls initially have large outward peculiar velocities.34 Previousauthors estimated the minimum homogenization time by determining the 1st-crossingtime. They found it to be6,8,9∆tc ≃H−1out(ti) c−1Rwall(ti)H−1out(ti)!2.

(8.2)This estimate was based on the assumption that a superhorizon-sized void would con-formally expand with spacetime. The present work throws considerable doubt on this27

assumption and its implications (i.e. a long thermalization time for superhorizon-sizedvoids) for several reasons.

First, the 1st-crossing time is much shorter than this estimate.For ∆tc/H−1out(ti) > 1, this estimate is offby the factorqρin(ti)/ρout(ti) (see Eqn(8.1)).Since the minimum thermalization time is the 1st-crossing time, this implies that thermal-ization and homogenization may happen much quicker than suggested. This would haveprofound effects upon our understanding of the evolution of superhorizon-sized voids inthe early universe.

Second, the inner edge of the void collapses (not expands) in roughlythe 1st-crossing time, even (apparently) for voids with large outward peculiar veloci-ties. This collapse occurs because the wall decelerates from the positive wall pressure( ˙U ∝−p′ < 0).

If the wall pressure were negative instead (e.g. the wall pressure of avacuum bubble), the wall would accelerate outward.

Thus, fluid from at least part of thewall rushes into the void, partially “filling it” at the 1st-crossing time. Because the lighttravel distance is still large there, the void may still fill up in a relatively short amountof time.

Third, the void wall can cease expanding by many mechanisms. Because fluidrushes into the void, large velocity gradients are formed in the original wall area which“pull” on the wall and can slow it down.

In addition, when the peculiar wall velocity ismuch larger than the gravitational velocity, the deceleration of the wall is proportional tothe velocity squared which can be enormous. Finally, an outgoing wave will be dampedand slowed down in any case, by virtue of the volume effect caused by the sphericalgeometry.

Thus even if the wall initially expands outward conformally (or faster), it maystop doing so rather quickly, as the simulations appear to suggest. In any case, even ifthe outer part of the wall can carry out some mass for a long period of time, the interiorregion will continue to fill in and homogenize.

Since overdensities will most likely formnear the origin after collapse, perturbations, small waves and/or overdensities will mostlikely be around at recombination. These perturbations could have interesting ampli-tudes, and therefore could alter the standard Harrison-Zeldovich density spectrum.

Inany case, because the qualitative void evolution picture previously suggested is at leastpartially incorrect, it is reasonable to expect the estimated thermalization time to beerroneous; the “big-bubble problem” might not be a problem after all.It is important to point out that in order not to distort the microwave background,the thermalization time only needs to decrease somewhat from the estimate given byEqn(8.2). The largest superhorizon-sized void from minimalist inflation is roughly c−1Rwall(ti)/H−1out(ti) ≃1027.8 If the initial time after reheat is ti = 10−33 seconds (TRH =1014GeV), then at recombination (trec = 1012 seconds), this void would not have ther-malized according to Eqn(8.2).

Using these initial conditions with Eqn(8.2), it is easyto see that voids for which c−1Rwall(ti)/H−1out(ti) >∼1023 were traditionally troublesome.However, since 1023 is such an enormous number, a very slight change in the functionalform for the thermalization time, ∆tH, can cause thermalization to occur before recom-bination. For instance, if we suppose that the thermalization time can be written as∆tH ≃c−1H−1out(ti)(Rwall(ti)/H−1out(ti))p, then an acceptable range for p is p <∼1.6, which isnot dramatically different from 2.

The point is that the “big-bubble problem” can onlybe resolved by accurate knowledge of the thermalization time, because estimates offby28

a somewhat small amount can lead to erroneous conclusions.Of the many questions which remain unanswered, the most important is to deter-mine what happens to a superhorizon-sized void after it collapses. If it thermalizes andhomogenizes, at what time does this happen?

As mentioned previously, at the time of1st-crossing, the energy density within the former void region is greater than ρin(ti), butis still less than ρout(t). In addition, the bulk motion of much of the original wall fluidis inward so that spacetime is collapsing and not expanding there; thermalization andhomogenization has not occurred.

However, since the “void” is still relatively underdenseat this time, the light travel distance will be comparatively large, thereby allowing forthe possibility of relatively fast thermalization and homogenization. We have found thatspecial relativistic voids eventually nearly homogenize after creating an overdense regionat the origin.

In this case, all but the fluid near the origin becomes roughly homogeneousfairly quickly. The fluid near the origin however, is overdense for a relatively long time.If general relativistic voids behave similarly, then after collapsing, an overdense regionwould form at the origin.

This would take a relatively long amount of time to diffuseaway, and would eventually result in a nearly homogeneous and isotropic universe. In anycase, since this overdense region would form very quickly, it is unlikely that empty voidswould be around any time near recombination, as previously suggested.32,8,9 Instead,overdensities, perturbations and small waves with strange fluid velocities may be.

Futurework will study the consequent evolution of a general relativistic void after collapse,33as well as the evolution of a superhorizon-sized void from first-order (e.g.extended)inflation.34AcknowledgementsThe author wishes to thank her thesis advisor, Rocky Kolb, for suggesting this prob-lem and for helpful discussions. She also wishes to thank Andrea Malagoli, Rick Watkins,Sidney Bludman, Michael Turner, Marc-Mordecai Maclow, Paul Steinhardt and Tsvi Pi-ran for interesting and useful discussions.

She would also like to thank Stephen Arendtfor useful discussions and help with the manuscript. This work was supported in part byNASA under Grant NAGW-2381 at Fermilab and by the DOE at Chicago and Fermilab.Appendix A: Transformation of the FRW metric to Eulerian CoordinatesAs discussed in Section IIA., the metric in Eq.

(2.1) can be subjected to transforma-tions of the type t = f1(t′, r′) and r = f2(t′, r′) without altering its spherical symmetry.If we set r2 = k(t, r), then we can rewrite Eq. (2.1) asds2=−(N2 −S2)dt2 + 2Sdtdr + h2dr2 + r2dΩ2,(A.1)where we have introduced the lapse and shift functions as N(t, r) and S(t, r), respectively.This is the Eulerian gauge with synchronous coordinates if we choose S = 0.

As we shallsee, when the FRW spatially flat metric is transformed in this manner, the resultingmetric has a coordinate singularity at the horizon. Thus, the study of superhorizon-sizedvoids is not possible in these coordinates.Consider the FRW k = 0 metric given by Eqn (2.3).

Letting er = ra(t), we can rewrite29

this asds2 = − 1 −er2 ˙aa2!dt2 −2er ˙aaderdt + der2 + er2dΩ2,(A.2)where we set c = 1 everywhere in this appendix. This is the Eulerian asynchronousform of the FRW metric.

We notice that this takes the form of Eq. (A.1) with N = 1,S = er˙a/a and h = 1.

Since a(t) ∝tξ, where ξ = 1/2 and 2/3 in radiation (p = ρ/3)and matter (p = 0) dominated universes respectively, we can eliminate the derdt term bydefiningt = ξer/g(et, er)(A.3)and ensuring that the condition∂g∂er = g ξ + (1 −ξ)g2ξer(1 −g2)(A.4)hold. Defining κ ≡(1 −ξ)/ξ and α ≡1/[2(1 −ξ)], this last expression can be integratedto obtainh(et)er = g/(1 + κg2)α,(A.5)where h(et) is an arbitrary (integration) function of et only.

The metric then becomesds2 = − herg!4 (1 + κg2)2(α+1)1 −g2det2 +11 −g2der2 + er2dΩ2. (A.6)This metric has a coordinate singularity at g = 1, which from Eqn(A.3) occurs atr = t/a(t) ξ−1.But this is the comoving radius of the Hubble radius (rHOR), sinceH−1 = t/ξ = rHORa(t).

Therefore, we have shown that for the FRW models expressed inEulerian, synchronous coordinates, a coordinate singularity occurs at the Hubble radius.As an interesting example we take ξ = 1/2, which corresponds to radiation-domination(p = ρ/3). In addition, we choose h ≡ξ/et.

Using Eqs. (A.5) and (A.3) together wither = r a(t) = r a(ti)(t/ti)1/2, we find that α = κ = 1, et = t + (a(ti)r)2/(4ti) andg = (1−q1 −4(er/et)2)/(er/et).

In addition, the metric simplifies to ds2 = −(1−g2)−1det2+(1 −g2)−1der2 + er2dΩ2. Note that even though the energy density is spatially constanton a hypersurface t = constant, on a hypersurface of constant et, ρ ∝t−2 = 4g2er−2 =4et2(1 −q1 −4(er/et)2 )2/er4 which is not spatially constant.Appendix B: Differencing of the EquationsIn this appendix, we list the expressions as they are differenced in the code.

We usethe shorthand notation ∆G ≡Gji+1 −Gji or ∆G ≡Gji −Gji−1 for forward or backwarddifferencing, respectively. The equations involving spatial derivatives are˙U=−Φ GNMR2+ 4πGN(p + Q)Rc2!−3 ΓΦR2w∆(r3)4π∆(p + Q)˙n=−nΦ∆(R2U)R2∆R˙ǫ=−4πΦ(p + Q)∆(R2U)Γr2∆r,(B.1)30

We note that there are many other schemes which could potentially be used. (For ex-ample, one could write ∆r3/3 instead of r2∆r in the expression for ˙ǫ ).

Although thesevariations can produce and propagate the shock and rarefaction waves, they can not han-dle the bounce of the shock at the origin. The only set of difference equations which wehave found to do this are given above.

For the Tolman-Bondi pressureless dust models(Section V), it is found that more accurate solutions are obtained near the origin whendifferencing the ˙n equation as follows: ˙n = −(n2Φ∆(R2U))/(Γr2∆r). The Tolman-Bondifigures shown in this paper, however, use the difference equations as written in Eqns(B.1).Appendix C: General Relativistic Jump ConditionsIn this appendix, we derive the jump conditions for general relativistic shocks.

Thisclosely follows the derivation by May and White (1965)14. These conditions will then beused to derive a more general artificial viscosity expression.

In addition, the shock jumpconditions for ultra special relativistic shocks are derived.Let the variables a and b represent labels for comoving observers ahead and behindthe shock, respectively. If each observer measures the invariant interval separating thesame two events on the world surface of a shock, using Eqn(2.2), we obtain[ds2] = [−c2Φ2dt2 + Λ2dr2 + R2dΩ2] = 0,(C.1)where [G] ≡Ga −Gb.

Because the shocks are radial, we can choose the two events tohave dr = dt = 0 but dΩ̸= 0. Therefore, R is continuous across the shock: [R] = 0.Now consider two events with dr ̸= 0 and dt ̸= 0.

Then from Eqn(C.1),[c2Φ2 −M2s /n2] = 0,(C.2)where rs is the position of the shock in comoving coordinates, drs/dt is the shock “speed”,and Ms ≡f (drs/dt)/(4πR2) from Eqns(2.19) and (2.5). This is the first of the jumpcondition equations.It can be shown that the jump conditions for the Schwarzschild metric (but not forthe comoving metric) are 35,14[Tνµ∂g/∂xµ] = 0,(C.3)where g is the equation for the world surface of the shock.

Therefore we must first findthe jump conditions in the Schwarzschild frame and then transform back to the comovingframe. The Schwarzschild metric is ds2 = −c2A2dT 2 + B2dR2 + R2dΩ2, where R is nowthe “Eulerian” coordinate radius.

Since R is a coordinate, [R] = 0, and using the sameargument as above for the two observers,[c2A2 −B2S2] = 0,(C.4)where S is the shock “speed” in this frame: cS ≡dRs/dT. For the Schwarzschild metric,the solution for the metrics functions are well known: B2 = 1 −2MG/(Rc2) = A−2where M(R, T) = 4πc−2 R R0 ρR2dR.

As long as ρ is not infinite in the shock, the mass iscontinuous across the shock, [M] = 0. Then [B] = [A] = 0.

And using Eqn (C.4), we31

find that [S] = 0. We conclude that in the Schwarzschild frame, the metric componentsand the shock “speed” are continuous across the shock: [A] = [B] = [R] = [S] = 0.We now derive the jump conditions in the Schwarzschild frame.

If a shock is located atposition Rs(T) at time T, the equation for the world surface of a shock is g = Rs(T)−R =0. In addition, the perfect fluid stress-energy tensor in the Schwarzschild frameo isT ′µν = c−2(ρ + p)g′µλu′µu′λ + pg′µλg′νλ = nwg′µλu′µu′λ + pg′µλg′νλ,(C.5)and the conservation of mass equation is [nu′ν∂g/∂xν] = 0.

35,14 Using Eqns (C.3), thejunction conditions becomehT TT S −T RTi= c2 h(−A2(u′T)2nw + p)S + nwA2u′Tu′Ri=0(C.6)hT TR S −T RRi=hSnwB2u′Tu′R −(nwB2(u′R)2 + p)i=0(C.7)hnSu′T −nu′Ri=0. (C.8)We would like to express these equations in terms of the comoving metric functions.First, we can rewrite the shock velocity in the comoving frame as cS = (Rt +Rr ˙rs)/(Tt +Tr ˙rs).

In addition, we can relate the metric functions from the comoving frame to those inthe Schwarzschild frame by g′µν = (∂x′µ/∂xσ) (∂x′ν/∂xλ) gσλ. We then obtain (cA)−2 =(cΦ)−2T 2t −Λ−2T 2r , B−2 = −(cΦ)−2R2t + Λ−2R2r and (cΦ)−2TtRt = Λ−2TrRr.

In addition,the “mass” M(r, t) contained within comoving radius r is continuous across the shock,since [M] = [nR2(1 + ǫ/c2)]dR⇒dR→0 0. Using Eqn (2.12), we find that [Γ2 −U2/c2] = 0.We can use these relations to calculate the 4-velocity as measured in the Schwarzschildframe.

Since the 4-velocity in the comoving frame is uµ = (−cΦ−1, 0, 0, 0), and the ve-locity transforms as u′λ = (∂x′λ/∂xσ) uσ, u′λ = −Φ−1(c ˙T, ˙R, 0, 0), where ˙T = ∂T/∂tand ˙R = ∂R/∂t. Since U ≡˙R/Φ, and using the fact that the 4-velocity is normal-ized to u′λu′λ = −c2 = c2A2(u′T)2 −B2(u′R)2, we can rewrite the 4-velocity as u′λ =−(A−1√B2U2 + c2, U, 0, 0).

In the special relativistic limit (GN = 0), A = 1 and B = 1and this becomesu′λ = −((U2 + c2)1/2, U, 0, 0) = −(cΓ, Γv, 0, 0),(C.9)where we have defined v ≡U/Γ so that Γ = (q1 −(v/c)2)−1. In the special relativisticlimit, then, v is the fluid particle’s radial velocity, Γ is usual gamma-factor (i.e.

energyper particle mass), and U is the particle momentum per particle mass.We can now rewrite Eqn (C.6)-(C.8) and [S] = 0 in terms of the comoving metricfunctions using the fact that Tr/Tt = c−2UΛ/(ΓΦ) and [ ˙rs] = 0. We obtain[UMs/(nc2) + ΦΓ]=0(C.10)[c2MsΓ(1 + ǫ/c2) −pUΦ]=0(C.11)[MsU(1 + ǫ/c2) −pΓΦ]=0(C.12)[MsΓ/n + ΦU]=0.

(C.13)oWe denote quantities in the Schwarzschild frame by primes.32

Eqn (C.2) and Eqns (C.10)-(C.13) make up the required shock conditions. One of themhowever, is redundant.It can be shown that in the nonrelativistic limit, the aboveconditions reduce to the Lagrangian shock jump conditions given in Ref 21.The case considered by May and White (1967) is the penetration of a shock into anon-relativistic medium.

We set GN = 0 and ǫa = pa = Ua = 0 and take ǫ = p/[(γ −1)n].Using Eqns (C.11) and (C.12) with the fact that U2b = Γ2b −1, the energy behind theshock is found to be ǫb = Γb −1. Using Eqns (C.12) and (C.13), we eliminate Φb andfind that η ≡nb/na = (1 + Γbγ)/(γ −1).

From Eqn (C.12), we find thatM = (γ −1)nbǫbΦb/Ub,(C.14)and plugging this into Eqn (C.2), we determine that Φb = 1/(1 + γ(Γb −1)). Finally, weplug this into Eqn (C.14) and find Ms = cnaq(Γb −1)/(Γb + 1) (1+Γbγ)/[1+γ(Γb−1)].

(It turns out that Eqn (C.10) gives no new information). Thus, all quantities behind theshock can be expressed in terms of Γb, the “strength” of the shock.How does this affect the artificial viscosity needed to smooth the shock?Von-Neumann’s viscosity22 is Q = k2n(U′)2dr2 for U′ < 0 and 0 otherwise, with the shockbeing spread over roughly k2 grid points.

Over the k2 points, Q ∼n(Ub −Ua)2 = nU2b =c2n(Γ2b −1) = c2n(Γb −1)(Γb + 1). Using the fact that ǫb/c2 = Γb −1, Q ∼nǫbΓb ∼pΓwhich is greater than the pressure p by the factor of Γ ≥1.

Since Q should be of order thepressure in the shock region, we divide by Γ: Q = k2n(U′)2dr2/Γ for U′ < 0 and Q = 0otherwise. In the nonrelativistic limit, Γ = 1 and Von-Neumann’s artificial viscosity isobtained.In this paper, we are more interested in the case where the fluid on both sides ofthe shock is ultrarelativistic.

We first set GN = 0 and Ua = 0. Then Eqn(C.2) andEqns(C.10)-(C.13) can be manipulated to yieldUb=qΓ2b −1(C.15)M=cnaηqΓ2b −1 Φa / (Γbη −1)(C.16)η=Γb(1 + γǫb/c2) −(1 + ǫa/c2) / [ǫb/c2(γ −1)](C.17)ǫb/c2=−1 + Γb(1 + ǫa/c2) + (γ −1)(Γb −η−1)ǫa/c2(C.18)Φb=|η −Γb| Φa / (Γbη −1).

(C.19)(It can be shown that these equations reduce correctly in the ǫa →0 limit). Now, if weset γ = 4/3 and take the ultrarelativistic limit, ǫb/c2 > ǫa/c2 >> 1, Eqn (C.17) becomesη = 4Γb{1−3ǫa/(4Γbǫb)} ≃4Γb.

Then the equations can be approximately solved to giveǫb≃4ǫaΓb(1 −(4Γb)−2) / 3,(C.20)which to lowest order is ǫb ≃4ǫaΓb/3. Therefore, ǫb depends not only on Γb but also onǫa.

We generalize the artificial viscosity toQ=k2n1 + ǫ/ (Γc2)(U′)2dr2/ Γfor U′ < 0Q=0otherwise. (C.21)33

For a strong shock then, Q ∼nǫU2/(Γ2c2) ∼nǫ ∼p, as desired. And when c2 →∞, thisbecomes VonNeumann’s original expression.

Note that because this was derived in thelimit of strong shocks, it does not work as well for weaker special or general relativisticshocks.33We now write down the exact solution of Eqns (C.15)-(C.19) for conditions behind theshock in terms of Γb only and quantities in front of the shock. Because ρ = nǫ in this limit,we can rewrite Eqns (C.17) and (C.18) as η = 4ρbΓb/(3ρa+ρb) and η = (3ρb+ρa)/(4ρaΓb)respectively.

Setting them equal, we can solve for ρb. Then, all others quantities aredetermined.

They areUb=qΓ2b −1(C.22)ρbρa=16Γ2b −1061 +q1 −36 (16Γ2 −10)−2(C.23)M=± naΦa(3ρb/ρa + 1)4Γbq3ρb/ρa(C.24)η≡nbna=4Γb3(ρb/ρa)−1 + 1(C.25)ǫbǫa=3 + ρb/ρa4Γb(C.26)Φb=ǫaǫbΦa = 3 + (ρb/ρa)−14ΓbΦa. (C.27)We note from Eqn (C.27) that [Φw] = 0.Appendix D: Nonrelativistic Shock Tube ProblemIn this section, we sketch the derivation of the slab shock tube solution for non-relativistic fluids.29 We start with a fluid in which the pressure p and specific volumeV ≡1/n (where n is the number density) to the left and right of x0 are p1, V1 and p2, V2resectively—thus, p and V are discontinuous across x = x0.

We assume here that ourcoordinates are oriented so that p1 < p2 and V1 > V2. In addition, the fluid is initially atrest: v1 = v2 = 0.

Because there is no scale to the problem, the solution is a function ofx/t only, where x is the Eulerian coordinate. At time t later, there are 4 regions.

Region1 and 2 are at rest with p = p1, V = V1 and v = v1, and p = p2, V = V2 and v = v2respectively.Separating region 1 and 3 is a shock wave, with position xd. Region 3contains the fluid behind the shock, and region 3′ is wedged between regions 3 and 2 andcontains the rarefaction wave.

The boundary between region 3 and 3′ is called a contactdiscontinuity and is at position xc. The velocity and pressure are constant across thisboundary: p3′ = p3 and v3′ = v3.

The number density, however, is not. Finally theboundary between region 3′ and region 2 is located at xa.

We will calculate p, V and vfor regions 3 and 3′ at time t, assuming that ǫ = pV/(γ −1). We will solve for thefollowing unknowns: p3, v3, n3, n3′ and p, n, and v in the rarefaction wave.Across a shock front, v1−v3 =q(p3 −p1)(V1 −V3) and ǫ1−ǫ3+ 12(V1−V3)(p1+p3) = 0.34

Using ǫ = pV/(γ −1), they becomeV3=V1[(γ + 1)p1 + (γ −1)p3]/[(γ −1)p1 + (γ + 1)p3](D.1)v3=−(p3 −p1)q2V1/[(γ −1)p1 + (γ + 1)p3](D.2)The similarity solution for a rarefaction wave is 20p=p2[1 −(γ −1)|v|/(2c2)]2/(γ−1)(D.3)V=V2(p2/p))1/γ(D.4)v=2(c −c2)/(γ −1)(D.5)|v|=2(c2 −x/t)/(γ + 1),(D.6)where the local speed of sound is given by ci = √γpiVi. We can equate the speed of thefluid v3 to the fluid velocity in the rarefaction wave between regions 3 and 3′.

We findv3′ = 2c2[(p2/p3)(1−γ)/(2γ) −1]/(γ −1)(D.7)Combining Eqns (D.2) and (D.7), we find thatp3 = p1 +1γ −1s2γp2n1n2q(γ −1)p1 + (γ + 1)p3h1 −(p3/p2)(γ−1)(2γ)i. (D.8)This equation can now be solved numerically for p3.

We can then determine v3, n3, and n3′using Eqns (D.2), (D.1) and (D.5). In addition, we find p, n and v from Eqns (D.3),(D.4) and (D.5).

in the rarefaction wave.The final step to the solution is determining the location of the boundaries separatingthese regions.Since the rarefaction wave is moving at the sound speed in region 2,xa = tc2. Since the velocity at b is vb = v3′ from Eqn (D.7), we then use Eqn (D.6)to find its position: xb = tc2[1 −(γ + 1)/(γ −1) [1 −(p2/p3)(1−γ)(2γ)]].

The contactdiscontinuity moves with the fluid and its position is thus xc = v3 t. Finally, the shockposition can be determined by transforming to a frame in which the shock is constantand using mass conservation. We then transform back to find the the velocity of theshock to be vs = (n3/(n3 −n1))v3 so that xd = vst.REFERENCES1.

A. Guth, Phys. Rev.

D 23, 347 (1981); A. Guth and E. Weinberg, Nuc, Phys., B212, 321, (1983).2. E. W. Kolb and M. S. Turner, The Early Universe, (Addison-Wesley, RedwoodCity, Ca., 1990).35

3. D. La and P. J. Steinhardt, Phys.

Rev. Lett.

62, 376 (1989); Phys. Lett.

220B, 375(1989).4. S. Hawking, J. Stewart and I. Moss, Phys.

Rev. D 26, 2681 (1982); R. Watkins andL.

Widrow, Nucl. Phys.

B374, 446 (1992).5. Attempts have been made in the infinitely thin-wall limit to obtain solutions forthe evolution of very large bubbles during extended inflation.

(D. Goldwirth andH. Zaglauer, Phys.

Rev. Lett., 67, 3639, (1991).) They found that under certainconditions, large bubbles could collapse.

However, the very assumptions used in thecalculation have been found to be in conflict with the value of Newton’s constanttoday (A. Liddle and D. Wands, Phys. Rev.

D 46, 3655 (1992)).6. E. J. Weinberg, Phys.

Rev. D 40, 3950 (1989).7.

D. La, P. J. Steinhardt, and E. Bertschinger, Phys. Lett.

231B, 231 (1989)8. M. Turner, E. Weinberg, L. Widrow, Phys.

Rev. D 46, 2384 (1992).9.

A. Liddle and D. Wands, Mon. Not.

Roy. Aston.

Soc., 253, 637, (1991).10. R. Tolman, Proc.

Natl. Acad.

Sci., 20, 169, (1934); H. Bondi, Mon. Not.

R. Astron.Soc., 107, 410, (1947).11.

K. Maeda, Gener. Relat and Gravit, 18, 931, (1986);12.

L. Landau and E. Lifshitz, The Classical Theory of Fields, (Permagon Press, 4thedition, 1975).13. C. Misner and D. Sharp, Phys.

Rev , 136, B571, (1964); W. Hernandez and C.Misner, Astro. J.

143, 452, (1966); M. May and R. White Phys. Rev, 141, 1232,(1966).14.

M. May and R. White, Meth. Computat.

Phys, 73, 219 (1967).15. P. Schinder, S. Bludman and T. Piran, Phys.

Rev. D, 37, 2722 (1988).16.

R. Richmyer and K. Morton, Difference Methods for Initial-Value Problems, Inter-science Publishers, New York(1967).17. The temperature can be defined locally since each spacetime point is surroundedby a small region which is in approximate thermal equilibrium and for which theBolzmann distribution is approximately independent of direction.1818.

S. Weinberg, Gravitation and Cosmology, John Wiley & Sons, New York, (1972).19. B. Schutz, A first course in General Relativity, (Cambridge University Press, Cam-bridge, London, 1985)36

20. L. Landau and E. Lifshitz, Fluid Mechanics, Pergamon Press, London, (1959).21.

R. Courant and K. Friedrichs, Supersonic Flow and Shock Waves, Springer-Verlag,New York, (1985).22. J. VonNeumann and R. Richmyer, Jour.

of Applied Physics, 21, 232, (1950).23. D. Bernstein, D. Hobill and L. Smarr, Dynamical Spacetimes and Numerical Rela-tivity ; proceedings, ed.

by J. Centrella, Cambridge, New York, (1986).24. This is strictly only true for compensated voids, but is also approximately true foruncompensated voids for some amount of time.25.

Note that in order for the fluid to begin and remain in local thermal equilibrium,we must require the interaction length of the mediating particle to be smallerthan the grid spacing at any time t. If λ is the proper interaction length, thenλ < Λji∆r = ∆Rji/Γji. If a gauge boson is the mediating particle (e.g.

a photon),then λ ≃(α2T)−1 ≃g1/4∗/(α2ρ1/4), where√4πα is the gauge coupling, ρ = π2g∗/30and g∗is the number of effectively massless degrees of freedom.226. J.

Stone et. al., Astrophys.

J., 388, 415 (1992).27. R P A C Newman, Class.

Quantum Grav., 3, 527, (1986)28. E. Bertschinger, Ap.J.

Suppl., 58, 1, (1985)29. G. Sod, J. Comput.

Physics, 27, 1 (1978).30. Sedov, Similarity and Dimensional Methods in Mechanics, Academic Press, NewYork, 1959.31.

K. Winkler and M. Norman, Astrophysical Radiation Hydrodynamics, ed. by K.Winkler and M. Norman, D. Reidel Publishing Co., Boston, 1986.

(pg. 101).32.

D. La, Physics Letters, B 265, 232, (1991).33. S. Vadas, Numerical Evolution of General Relativistic Voids II, work in preparation.34.

S. Vadas, Evolution of SuperHorizon-Sized Voids from Extended Inflation, work inpreparation.35. J. Synge, Relativity: The General Theory, North-Holland Publ., Amsterdam (1960).FIGURE CAPTIONSFig.

1: The evolution of a compensated, pressureless void with radius Rwall(ti) =333cH−1out(ti) and for velocity U(ti, R) =q2GNM/R. The energy density as a function of37

RR(ti, RB)/R(t, RB) ≃R(ti/t)2/3 is shown at the initial time ti = 1 and at times t = 10,100 and 300.Fig. 2: Same as for Figure 1 but with velocity U(ti, R) = c−1Rq8πGNρ/3.

A dense,thin, outward-moving shell is formed. Shell-crossing occurs at t = 1.048.Fig.

3a: The evolution of a pressureless, uncompensated void with radius Rwall(ti) =.0067cH−1out(ti) and velocity U(ti, R) = c−1Rq8πGNρ/3. Artificial viscosity is used to pre-vent shell-crossing.

The energy density as a function of RR(ti, RB)/R(t, RB) ≃R(ti/t)2/3is shown at the initial time ti = 1 and at times t = 50, 100, 400 and 1000.Fig. 3b: Same as Figure 3a but for a compensated void.Fig.

3c: The radius of the shell versus time for the results pictured in Figures 3a and3b. The predicted self-similar solutions are overlaid as a dotted and dashed line for anuncompensated (Rshell ∝t8/9) and compensated (Rshell ∝t4/5) void respectively.Fig.

4: The relative errror in the energy density versus RR(ti, RB)/R(t, RB) ≃Rqti/t for an initially homogeneous and isotropic FRW universe with p = ρ/3 at time∆t/ti = .1. The results are for f = .01, .005, and .0025 shown with solid, dotted anddashed lines, respectively.Fig.

5: The nonrelativistic shock tube for Rwall(ti) = 20, ∆Rwall(ti) = .01 andGN = 0. The dashed lines represent the exact slab similarity solution.Fig.

6a: The nonrelativistic evolution of an uncompensated void with GN = 0. Thedashed line is the initial time ti = 1, and the triangles with connecting lines is at timet = 2.5 after the inbound shock has bounced and the fluid has begun to settle down.Fig.

6b: Same as Figure 6a but for a compensated void.Fig. 7a: Collapse of a special relativistic, uncompensated void for Rwall(ti) = 1,∆Rwall(ti) = .02, T/µ ≫1, GN = 0, U(ti, R) = 0 and α = 10−4.

Shown is M, Γ, 4πρand Φ versus R for the initial time ti = 1 and for times t = 1.04 and 1.065 (before theshock reaches the origin).Fig. 7b: We plot 4πp versus R for special relativistic voids with initial varying wallthicknesses at the initial time (unlabeled) and at later times t = 1.04 and 1.07, 1.08 or1.065.

The parameters are the same as in Figure 7a but with ∆Rwall(ti) = .02, .04, .06and .1.Fig. 7c: Collapse, thermalization and homogenization of a special relativistic void.The parameters are the same as in Figure 7a.

Shown is the initial time ti and timest = 1.04 (void collapsing) and t = 3.0 (after collapse).38

Fig. 8: Collapse of a general relativistic, uncompensated void for Rwall(ti) = 1,∆Rwall(ti) = .02, T/µ ≫1, U(ti, R) =q2GNM/R and α = 10−4 at the initial timeti = 1 and at times t = 1.04 and 1.065.Figs.

9a,b,c: Collapse of general relativistic, uncompensated voids for Rwall(ti)/H−1out(ti) = 25, ∆Rwall(ti) = 1, T/µ ≫1, U(ti, R) =q2GNM/R. Figures 9a, 9b and 9cshow the pressure versus RR(ti, RB)/R(t, RB) ≃Rqti/t for α = 10−4, α = 10−6 andα = 10−10, respectively.Figs.

10a,b: Collapse of general relativistic, uncompensated and compensated voidin 10a and b, respectively. Rwall(ti)/H−1out(ti) = 250, ∆Rwall(ti) = 10, T/µ ≫1, U(ti, R) =q2GNM/R and α = 10−10 for both.

Shown is the pressure versus RR(ti, RB)/R(t, RB) ≃Rqti/t.Fig. 11: Collapse of a general relativistic, compensated void for Rwall(ti)/H−1out(ti) =250, ∆Rwall(ti) = 10., T/µ ≫1, U(ti, R) = c−1Rq8πGNρ/3 and α = 10−10.

Shown isthe pressure versus RR(ti, RB)/R(t, RB) ≃Rqti/t.Fig. 12: 1st-crossing times versus Rqti/t for general relativistic voids with initialradius Rwall(ti) = c1023H−1out(ti) and outside temperature Tout = TRH, where TRH is thereheat temperature.

The dash-dot lines are for Φin(ti) = Tout(ti)/Tin(ti) = 5 × 1011, 1015and 1020 and the temperature at which each void potentially thermalizes is T = 500,106 and 1011GeV respectively. We also plot (long dashes) the evolution of a conformallystretched perturbation with the same initial size, rCM, and the Hubble radius (shortdashes) outside the void, rHOR.

If TRH = 1014GeV, then recombination occurs at t/ti =1045 (or at temperature Trec ≃3 × 10−9GeV), where rCM and rHOR intersect.39

---

출처: arXiv:9303.295 • 원문 보기