Vorticity generation in cosmology and the role of shell crossing

There is no source for cosmic vorticity within the cold dark matter cosmology. However, vorticity has been observed in the universe, especially on the scales of clusters, filaments, galaxies, etc. Recent results from high-resolution general relativistic N-body simulation show that the vorticity power spectrum dominates over the power spectrum of the divergence of the peculiar velocity field on scales where the effective field theory of large-scale structure breaks down. Incidentally, this scale also corresponds to the scale where shell-crossing occurs. Several studies have suggested a link between shell crossing in the dark matter fluid and the vorticity generation in the universe, however, no clear proof of how it works within general relativity exists yet. We describe for the first time how vorticity is generated in a universe such as ours with expanding and collapsing regions. We show how vorticity is generated at the boundary of the expanding and collapsing regions. Our result indicates that the amplitude of the generated vorticity is determined by the jump in gradients of the gravitational potential, pressure and the expansion rate at the boundary. In addition, we argue that the presence of vorticity in the matter fields implies a non-vanishing magnetic part of the Weyl tensor. This has implications for the generation of Maxwell's magnetic field and the dynamics of clusters. The impact of accelerated expansion of the universe and the existence of causal limit for massive particles are discussed

The standard cosmological model also known as the Lambda-Cold Dark Matter (ΛCDM) model has been successful in explaining some of the observed large-scale features of the Universe, for example, the observed anisotropies in the cosmic microwave background radiation [1,2].The model assumes that despite the inhomogeneous distribution of structures visible to an observer, the universe is well-approximated by the Friedmann-Lemaıtre-Robertson-Walker (FLRW) spacetime on all length scales.It asserts that large-scale structures can be described as small perturbations on top of the background homogeneous, isotropic FLRW spacetime.These perturbations can be decomposed into three modes: scalars, vectors and tensors.At the linear level, these modes propagate independently [3,4].The scalar perturbations cannot induce a rotational part of the peculiar velocity.The vector perturbations can if there is a source initially but even if there is a source it decays very fast as the universe expands [5].The tensor perturbations can also induce the cosmic vorticity field but the amplitude is very small [6].At non-linear order, the evolution of the scalar perturbations can generate vector and tensor perturbations [5,7], however, their amplitude cannot explain the observed vorticity in large-scale structures [8,9].Vorticity can also be sourced by the entropy perturbations but the adiabatic perturbation appears to be preferred by the current observation [10].
The cosmic vorticity field has been observed in galaxy clusters, filaments, galaxies, etc [11].It is well known that most galaxies rotate and that the angular velocities of neighbouring galaxies are correlated [12].They play an important role in determining the observed galaxy spin and alignment [13].On the solar system's scales, it dominates the dynamics of weather patterns [14], yet the evolution of the cosmic vorticity has no known source within the ΛCDM model of the universe.Within the general relativistic N-body simulation, the vorticity power spectrum was recently estimated and it was found that it dominates over the power spectrum of the divergence of the peculiar velocity field of the matter field on scales where the effective field theory of large-scale structure breaks down [15,16].In this N-body simulation, the vector and tensor perturbations were turned off.Also, it was shown conclusively in the paper that the measured vorticity field is in the peculiar velocity of the matter as measured by the Eulerian observer [17].These details are important for the discussion that will follow.
The work of [18] was the first to provide insights on how the generation of vorticity may be related to shell crossing of the matter fields.Using the N-body simulation, they showed that the vorticity field tends to peak in the outskirts of virialised structures.This result was confirmed by [19] using a different suit of N-body simulation.Hints of this were earlier discussed in [20], where the amount of vorticity generated after the first shell crossing in large-scale structure caustics was done.There is no theoretical understanding of the connection between the shell-crossing of the mater fluid element and the generation of vorticity in the universe.This lack of understanding motivated [21] to consider whether the pair-wise velocity of galaxies can explain the observed vorticity.Similarly, the contribution of higher cumulant of the phase-space distribution function was studied [22,23].Their conclusion of [21] appears to show that the contribution of the pair-wise velocity is insufficient to explain the measured vorticity field.
It is this gap that this paper plans to fill.We describe for the first time how cosmic vorticity field may be generated in the neighbourhood of "shell crossing singularity".We clarify the role of caustics or shell-crossing singularity in the generation of vorticity in cosmology.To achieve this, we develop a model of the universe that describes more consistently the expanding and collapsing regions of universe.The standard cosmology model, describes only the expanding regions neglecting the gravitationally bound regions that have decoupled from the Hubble expansion.We show that a consistent treatment of both regions is crucial to the generation of vorticity in the universe.The amplitude of the generated cosmic vorticity field depends on the difference between gradients of the gravitational potential, pressure and convergence of flow lines between expanding and collapsing shells of matter.
The consequences of nonzero vorticity are enormous and we explored a couple of them.Firstly, the existence of vorticity in the matter fields implies a non-vanishing magnetic part of the Weyl tensor associated with the matter flow velocity.The magnetic part of the Weyl tensor vanishes in Newtonian gravity, hence its measurement will constitute another test of general relativity.The nonzero magnetic part of the Weyl tensor could have implications for the existence of dark matter.This was studied in [24,25].In addition, a non-zero cosmic vorticity field could also have implications for the generation and propagation of Maxwell's magnetic field in clusters.Finally, we show that differences in the expansion and contraction rates of the expanding and collapsing regions respectively could help explain the late-time accelerated expansion of the universe.
The rest of this paper is organised as follows: we describe particle trajectory, fluid flow and shell crossing singularity in Newtonian gravity in section II.We extended the same treatment to general relativity where we describe the geodesic of a massive particle on curved spacetime.We also identify the point where a geodesic ceases to be geodesic using the focusing theorem in this section.We show how the existence of an apparent horizon or the causal limit allows changing to a more appropriate coordinate since the Jacobian determinant at the causal limit is non-zero.The change of coordinate is possible because of the inverse function theorem.We introduce the model of the universe that describes the expanding and collapsing regions with appropriate boundary conditions consistently in section III.We describe vorticity generation in section IV and discuss other obvious implications of the model in section V and conclude in section VI.

II. DYNAMICS OF MASSIVE PARTICLES ON CURVED SPACETIMES
In Newtonian gravity, the concept of massive particle motion is formulated in terms of the force the particles feel.In general relativity, massive particles travel along time-like geodesics.Geodesics are not globally valid, especially in curved spacetime.A curve could start out as geodesic but could cease to be geodesic in a finite time.Thus, the concept of geodesics on curved spacetime is only locally defined.In this section, we shall describe in detail how shell crossing singularity or caustics in cosmology is related to the breakdown of the geodesic.We also discuss how the formation of the causal horizon before caustic formation allows it to change to a more appropriate coordinate thereby avoiding the shell-crossing singularity.

II.1. Fluid flow in Newtonian gravity
The modelling of the large-scale structures of the universe depends crucially on the solution of the geodesic equation.The initial conditions for the Newtonian N-body simulations are set using the solution of the geodesic equation [26].The force among particles distributed all over the universe is calculated using geodesic equation [27].One key point that is usually not mentioned when this approach is adopted is that a curve is geodesic only locally.A geodesic can cease to be a geodesic within a finite time.This is usually not studied as a breakdown of geodesics in cosmology, however, it is a big area of research in Differential geometry [28].Particles moving with non-relativistic velocities in the weak gravitational field regime in Newtonian gravity, its trajectory is given by where ∇ i is a spatial derivative on Euclidean space, Φ is the gravitational potential, τ is proper time, x i is the position of the particle or fluid element at τ .x i is related to the initial position, q i , of the particle according to where Ψ i is known as the displacement field.It is customary to describe x i and q i as the Eulerian and Lagrangian coordinates respectively.Initially, that is at τ = 0, Ψ i (q, τ ) = 0.In the expanding universe, it is beneficial to work in comoving coordinates r i = x i /a(η), where a is the scale factor of the universe, η is the conformal time, it is related to the proper time according to dη = dτ /a(η).The gravitational potential is related to the mass density, ρ through the Poisson equation ∇ 2 Φ = 4πa 2 δρ m , where δρ is a perturbation in the mass-density with respect to a background FLRW value ρ.Imposing conservation of mass-density, that is ρ(τ )d 3 q = ρ(τ, x i )d 3 x, leads to where δ ≡ δρ/ρ is the density contrast or fluctuation in the mass-density and J = det [δ ij + Ψ i,j ] is the Jacobian of the transformation.According to Zel'dovich [29], the leading order approximation to J is given by where α, β and γ are eigenvalues and D 1 is the matter density growth function.The caustics occurs, i.e J → 0 whenever 1 − αD 1 (τ ) → 0. According to equation (3), the density contrast diverges at a caustics 1 + δ(x, η) → ∞.This is sometimes interpreted as an indication of a breakdown of Zel'dovich approximation.However, going beyond the Zeldovich approximation by adding higher-order perturbation theory terms does not resolve or remove the caustics [30].
In fact, it is more of an indication of a breakdown of the 'one-parameter' cosmological perturbation theory.It is important to mention that it is more of an indication of a breakdown of one-parameter cosmological perturbation theory because cosmological perturbation theory can be formulated in many-parameters by dynamically switching to a most suitable background spacetime [31][32][33].In the language of fluid dynamics, it is an indication of a breakdown of the single-stream approximation [34].It is also clear in this language that when the single-steam approximation breaks down, the natural progression towards progress is to introduce the two-stream approximation [35,36].
There have been several attempts towards developing a consistent two-stream approximation or two-parameters perturbation theory in cosmology [31-33, 37, 38].These attempts focus on the impact of the coupling between small and large-scale dynamics on the large-scale features of the universe.The results so far appear to show that the coupling could be important for the gravitational waves [39].Consistent treatment of the boundary condition for the scalar perturbations which dominate dynamics on small scales is missing.One surprising attribute of these approaches is that they assume a priori that Newtonian approximation is valid on small scales [38,40].This is one of the crucial points that we highlight here.Newtonian gravity does not appear to describe some of the critical events preceding the formation of caustics.For example, it has been pointed out earlier that a massive particle comoving with the Hubble expansion cannot influence dynamics within a virialised local environment if there exists a causal horizon [41].Newtonian gravity lacks the tools to describe the formation of causal horizons [42].There have been works attempting to justify the use of Newtonian gravity in cosmology on all scales [43].These claims are yet to account for some degrees of freedom in general relativity such as the magnetic part of the Weyl tensor that is fundamentally not contained in Newtonian gravity [44,45].The magnetic part of the Weyl tensor is non-zero if the vorticity is non-zero [46].
As a result of these, the rest of our discussion will be based on general relativity.It is possible to realise some of these features that precede the formation of caustics in Newtonian gravity after the fact.Furthermore, general relativity is consistent with the principle of least action.The principle of least action is central to our model building because it is straightforward to derive consistent boundary conditions for geodesics using the principle of least action.Some of the ideas we discuss here are similar to those employed in the study of gravitational memory effect but our approach is more fundamental [47].

II.2. Geodesic of relativistic massive particle and its evolution equation
The action of a relativistic massive particle is minus the rest energy times the change in time S = −E∆τ = −m∆τ , where ∆τ = τ f − τ i , where τ i is the initial time when the seed was created and τ f is the future time.The particle follows a maximal geodesic as it evolves from τ ini to τ f .The most essential point here is that the concept of a maximal geodesic is only locally defined in curved spacetime.Therefore, it is essential to ascertain the range of validity of geodesic.
Let γ denote a smooth time-like curve defined within an interval [τ ini , τ f ] on Pseudo-Riemann manifold M 4 in 4 dimensions.The massive particle action is given by where g ab is the metric of the spacetime on M 4 , L( γ(τ ), γ(τ )) is the Lagrangian and γ = dγ/dτ .In the second equality, we introduced x a , i.e the coordinates of the points on the manifold.We will drop the −m in equation ( 5) for the rest of the presentation to reduce clutter.Equation ( 5) is invariant under reparametrisation.For γ to be a geodesic within [τ ini , τ f ], it has to be the critical point of the infinitesimal variation of γ(τ ).Let x a (τ, s) = x(τ ) = sδx(τ ) be a variation, where s parameterises nearby curves s ∈ (−ϵ, ϵ).The central geodesic is given by xa (τ ) = x a (τ, 0).A variation is proper when all the nearby curves converge at the end-points On curved spacetime, a family of nearby curves could converge before the endpoint.When this happens a geodesic will no longer maximax the action if it is extended beyond this point.The point where a family of geodesics converge before the end-point is called a conjugate point.A conjugate point is a caustic since the Jacobian vanishes there.Let's introduce a time-like 4-velocity, u a and a deviation vector, ξ a that tracks the propagation of the nearby family of curves: It is well known that the first variation of an action is equivalent to taking the functional derivative of the action.In our case, the action in (5) has some amazing symmetry by the Noether theorem, i.e the Lagrangian does not depend explicitly on x a .This implies that ξ a can be Lie dragged along u a : where we have imposed the orthogonality condition ξ a u a = 0 in the last line.And in the third line, we imposed the conditions for proper variation of the action.One can check that the result in equation ( 10) is correct by simply taking the functional derivative of an action with an arbitrary Lagrangian which gives Putting the Lagrangian in equation ( 5) into equation (11) gives the same result.The critical point of these equations (i.e dS/ds| s=0 = 0 gives the geodesic equation u c ∇ c u a = 0.

II.3. Validity of geodesics, focusing theorem and horizon
To understand when a geodesic ceases to be a geodesic, we need a second derivative test.This is also called the second variation of the action.The second variation measures how fast nearby geodesics are expanding or contracting towards the central geodesics γ 0 (τ ).The second variation of equation ( 5) may be obtained by simply taking the derivative of equation ( 8) where we have introduce a shorthand notation for directional derivatives ξ c ∇ c = ∇ ξ and u c ∇ c = ∇ u .The index position on the first term can be switched with the help of the Ricci identity After some simplification, it gives Performing integration parts leads to where R a def is the Riemann tensor and Dξ a /Dτ = u a ∇ a ξ b .The equation of motion resulting from the first variation has been imposed.Further imposing the condition for proper variation ξ a (p) = ξ a (q) = 0, equation ( 16) reduces to The critical point (i.e d 2 S/ds 2 = 0) gives the geodesic deviation equation.
We can also obtain the same geodesic deviation equation in equation ( 18) using the following Lagrangian in equation (11) and setting ξ b ∇ b ξ c = 0: This connection will be important later.It is easier to extract information from equation ( 18) by decomposing the spacetime into temporal and the spatial part, the most consistent way to do this is to consider foliations where ξ a Lie dragged along the integral curves of u a [42].
where B ab = ∇ b u a .B ab measures the deformation of the curved spacetime in comparison to flat space.It can be decomposed into irreducible coordinate independent components [46,48] where A a is the acceleration, Θ describes the expansion of the one-parameter family of geodesics if Θ > 0 and contraction or collapsing of one-parameter family of geodesics if Θ < 0. σ ab called the shear tensor, describes the rate of change of the deformation of a one-parameter family of geodesics when compared to flat spacetime.ω ab is the vorticity tensor.it is an anti-symmetric tensor, It describes the twisting of a one-parameter family of nearby geodesics.We define also the scalar invariant of these tensors as follows, for the shear tensor σ 2 = σ ab σ ab /2 and the vorticity tensor ω 2 = ω c ω c /2 = ω ab ω ab /2), where ω a = 1 2 ε abc ω bc is a vorticity vector.h ab is the metric on the hypersurface.It is defined in terms of the metric g ab and u a where is the alternating tensor for the full spacetime [46].Note that B ab is related to the extrinsic curvature tensor according to Putting equation (20) into equation (18), we find that these geometric quantities satisfy the following equations [46,[49][50][51]: where C acbd is the Wely tensor and R ab is the Ricci tensor.Equations ( 23), ( 24) and ( 25) can also be derived using the Ricci identity.General relativity is needed to relate the Ricci tensor in equation ( 23) to the energy-momentum tensor.We make minimal assumptions about the form of the energy-momentum tensor.The Weyl curvature tensor may be decomposed further into the electric E ab and the magnetic part H ab with respect to u a where E ab and H ab are defined as E ab = C acbd u c u d and H ab = ε a cd C cdbe u e /2.E ab and H ab live on the hypersurface E ab u b = 0 = H ab u b .Note that E ab describes the impact of the tidal forces due to local massive distribution, while the magnetic part describes the tidal forces due to the twisting or stretching of spacetime along different directions.In the Newtonian limit, E ij ≈ ∂ i ∂ j Φ, where Φ is the gravitational potential and and H ab vanishes in the Newtonian limit.The most well-known consequence of the above equation is the prediction that the vorticity vanishes exactly.This can be seen by expressing the derivative in equation (25) in terms of the Lie derivative In the gravitational rest frame A a is a gradient of a scalar A a = ∇ a Φ, hence D [a A b] = 0. Therefore, irrespective of the coordinate system, ω ab vanishes if the initial vorticity is zero.The vanishing of the vorticity or the existence of the vorticity-free congruence implies u a is hypersurface orthogonal.This also means that u a maybe derived from the covariant of a scalar field S: u a = −∇ a S/||∇S||, where ||∇S|| is a normalization factor and u a is pointing in the future direction.The shear propagation equation may be written in a coordinate-independent form as The shear is sourced by the electric part of the Weyl tensor and the trace-free part of the Ricci tensor.
Finally, from equation (23) we can obtain the time-like geodesic version of the focussing theorem [42,52]From equation (29), it is clear that σ ab σ ab ≥ 0, for zero vorticity and assuming that the weak energy condition holds that is Integrating with respect to τ gives 1/Θ ≥ 1/Θ 0 + τ /3, where Θ 0 is the initial value of the expansion.Equation (30) describes the features of geodesics that must collapse to caustics.Such geodesics must be collapsing initially.We go into greater detail in sub-section II.4 to describe how this happens in cosmology.

II.4. Inevitability of more than one-parameter family of curves in a universe like ours
Within the standard model of cosmology, the cosmological inflation models are usually built on an FLRW background spacetime.The model predicts the initial conditions for the large-scale structures of the universe.The current observation suggests that the seeds of large-scale structure formation follow a Gaussian distribution [54].General relativity is a deterministic theory, hence, it is possible to probe models of cosmological inflation using the observations in the late universe [55,56].One could extend this argument to the focusing theorem(equation ( 30)), that the family of geodesics that collapses to form clusters, galaxies we see today are those that found themselves within over-dense regions at the initial time, while those that found themselves in under-dense regions evolve to form voids. Evolution histories of these regions are different as we will show in section III.The concept of tracer bias in modelling the clustering of large-scale structures is based on this idea [57,58].The discussion below goes into greater detail to describe the distinction between the one-parameter family of geodesics that collapses to clusters and the one-parameter family of geodesics that evolves to form voids.
Equation ( 20) is a linear first-order differential equation, hence, the solution at present time is related to its values at some point q in the past according to ξ i (τ, x) = J ij (τ, x)ξ j (τ ini , q) , where J ij is a Jacobi matrix and i, j indicates components.Putting ξ i (x) = J ij (x)ξ j (q) , in equation (20) gives where, K ab = h b c ∇ c u a is the extrinsic curvature of the hypersurface orthogonal to u a , Although the relationships between various components of J i j and K ik are important(see [59] in the case null geodesics), our interest at the moment is on the determinant which is given by where det[J] is the determinant of the Jacobian.In order to obtain this equation, we made use of the Jacobi formula [60], which expresses the derivative of the determinant of any matrix A whose inverse exists in terms of the adjugate of A and the derivative of A. In relation to equation (31), we are assuming that there are no caustics det[J](τ, x) ̸ = 0 .Integrating equation (32) gives At Θ = 0, the Jacobian becomes det[J](τ, x dp ) = det[J](τ ini , x dp ), where x dp defines a spatial location where Θ = 0.At this location, the property of the fluid element changes; it is incompressible fluid at x dp [14,61].Furthermore, equation (32) has the form of an autonomous differential equation (dy/dτ = F (τ, y)), hence one could argue that Θ = 0 is a critical point.In the observed universe, there are locations where Θ(τ, x dp ) = 0.This is usually measurable in peculiar velocity surveys [62][63][64].For the local group observer, this location is known as the zero-velocity surface and its radius has been determined precisely [62,65,66].The consequences of this for the supernova absolute magnitude tension were explored in [67,68].Our discussion here is more general, we treat these as critical location or causal horizon that indicates an end to a collapsing region of the universe and the beginning of an expanding region for an observer in a virialised region such as ours.The causal horizon can easily be determined by splitting Θ into two parts: where Θ H denotes the expanding part, i.e Θ H = 3H(this is determined by the background FLRW spacetime) and Θ L describes the local component.The expanding component is always positive, Θ H > 0, while the local component could be positive, negative and zero.Θ L = 0 implies a universe without large-scale structures or the large-scale structures have zero peculiar velocity relative to Hubble expansion, Θ L > 0, implies an equally expanding local regions.The relative dominance of Θ H and Θ L divides the universe into expanding and collapsing regions.The locations in the universe where the gravitational field of gravitationally bound structures dominate Θ L > Θ H (e.g.haloes, etc) define the collapsing regions.The star and galaxy formation happen within this region [41].Within this region, the oneparameter family of geodesics are converging to a singularity/caustic according to the focusing theorem, however, nonlocal effects may intervene to prevent a singularity formation [42].The locations where Θ L < Θ H (e.g.void, vacuum) are expanding and will continue to expand.These are the expanding regions.This dichotomy may be better understood by calculating Θ in a perturbed FLRW spacetime in comoving synchronous gauge where ψ and E are perturbed metric variables respectively.Note that we neglect the vector and tensor perturbations since their propagation is null-like.To the leading order, Θ is given by where "˙" is the derivative with respect to proper time, Θ H = 3H(τ ) and [69].Again, the expansion or contraction of nearby geodesics is determined by the relative dominance of Θ H and Θ L .Within the Halo model [70], it is possible to estimate x dp assuming spherical symmetry, i.e x dp = r dp r (r is a unit vector).So we can express the time derivative in terms of the radial derivative On the FLRW background spacetime, ∂r/∂τ can easily be evaluated for null geodesics.In the spherically symmetric LTB model, there exists a clear relationship between the time and the areal radius [71].This is a general property of inhomogeneous spacetime.However, we consider a much-simplified approach since the vanishing of Θ(τ, x)) is time independent, we can estimate ∂r/∂τ using the null geodesic relation.In this limit r dp is interpreted as the comoving distance to the zero-velocity surface for an observer centred at r = 0. Implementing this in the equation (37) gives where r is the comoving radial distance, c is the speed of light, δ m ≡ δρ/ρ = (ρ − ρ)/ρ and ρ is the matter density.Θ(τ, r dp ) = 0, when r dp = −cd ln ρ/(3Hd ln r).Given any halo density profiles, it is straightforward to calculate d ln ρ/d ln r [72].The plot of Θ as function of r is shown in figure 1, We made use of the NFW (Navarro-Frenk-White) dark matter density profile with the FLRW exterior [73].M⊗is the mass of the sun.Note that the causal horizon is much greater than the splashback radius.Both of these radii can in principle be measured.
There are three essential parts of Θ according to figure 1: • Expanding region(+): There are regions with a comoving distance greater than the causal limit r > r dp for a given gravitationally bound cluster.Within this region, the global Hubble expansion of spacetime dominates.
A typical example is void.
• Collapsing region (-): These are regions ( r < r dp ) that have gravitationally decoupled from the Hubble expansion because they are moving with a slower velocity to catch up with the Hubble expansion.Within this region, a one-parameter family of nearby geodesics are collapsing with respect to an observer in the expanding region.Also with respect to the observer in the expanding region, the geodesics within this region would appear to be converging to caustics.
• Boundary: This is a thin shell located at r = r dp .Within the spherical collapse model, It is related to the turn-around radius.It is the critical point of equation (32).We refer to it as the causal horizon for massive particles with velocities less than three times the global Hubble rate.Analysis of several observations and the N-body simulation of the local universe indicates that this scale exists and it is fundamental [74][75][76].
Finally, one key point to note here is that at Θ(τ, r dp ) = 0, the determinant of the Jacobian is a non-zero constant.In the next sub-section, we study the dynamics of a one-parameter family of time-like geodesics in the neighbourhood of the causal horizon.

II.5. Caustics and inverse function theorem
The existence of the causal horizon, i.e Θ(τ, r dp ) = 0 divides a family of time-like geodesics that starts at the same time in the past into two regions.At τ = τ ini hypersurface, one-parameter family of geodesics within r < r dp ) are converging while one-parameter family of geodesics in r > r dp are expanding.In this sub-section, we study the dynamics of a one-parameter family of time-like geodesics in the neighbourhood of the causal horizon by perturbing the geodesics around Θ(τ, r dp ) = 0 surface: where ∆τ is an infinitesimally small difference between τ and τ ′ .Under the infinitesimal perturbation, the Jacobian and the expansion scalar change according to Substituting in equation (32) and keeping only terms that are of linear order in ∆τ gives The time evolution of Θ is given by the Raychaudhuri equation ( 23) and imposing Θ(τ ′ , x dp ) = 0, leads to Assuming that the weak energy condition holds R αβ u α u β ≥ 0 and we know that σ ab σ ab is positive definite, then the second derivative must be negative indicating that x dp is a local maximum.This is a typical example of a ball rolling down a hill.Any slight perturbation in the particle position causes it to roll down hill.The ball accelerates as it rolls down the hill because the gravitational force is pulling it downwards.Most crucial lesson here is that the global Hubble expansion breaks down at x dp , therefore, an expanding coordinate system cannot be extended beyond the zero-velocity surface.Extending the geodesics that started out in an expanding spacetime beyond x dp will end up in a singularity or a caustic immediately after x dp , hence x dp is a boundary.The fact that the determinant of the Jacobian (det[J](τ, x)) is a non-zero constant at x dp provides hints on how to proceed.For example, by the inverse function theorem, it indicates that we can find another more suitable set of coordinates in the immediate neighbourhood of x dp on the collapsing region and join it seamlessly to the FLRW spacetime at x dp .We show in section IV that conditions for joining two families of geodesics across the zero-velocity surface.

III. MODEL OF THE UNIVERSE WITH COLLAPSING AND EXPANDING REGIONS
The analysis in sub-section II.5 shows that the one-parameter family of geodesics which describes the dynamics of massive particles in the expanding regions of the universe cannot be extended beyond the causal horizon.At the causal horizon, the Jacobian determinant is constant and by inverse function theorem, we can define another oneparameter family of geodesics to describe the dynamics in the collapsing region since the extension of the geodesics in the expanding universe into the collapsing region leads to a caustics.We describe in detail how to do this in the remainder of this section.

III.1. Dynamics of geodesics in both regions and their junction conditions
The diffeomorphism symmetry of general relativity in 4 dimensions allows freedom to choose coordinates.We want to find a smooth coordinate transformation that will smoothly join the expanding and collapsing coordinates.Therefore, we require that the spatial coordinates satisfy the following condition at the boundary where x a + is the coordinates in the expanding region, x a − is the coordinates in the collapsing region and Σ is the spatial hyoersurface.We parameterise the geodesics by a translated time parameter t = τ − τ ini (r dp ), where τ ini (r dp ), perhaps is related to the bang time when collapsing and expanding regions of the universe were delineated.This parameterisation allows us to place the boundary hypersurface Σ at t = 0. Thus, geodesics with positive t > 0 (τ > τ ini (r dp )) are in the expanding spacetime while the geodesics with t < 0 (τ < τ ini (r dp )) are in the collapsing spacetime.Therefore, we define the 4-velocities in the collapsing and expanding spacetimes as respectively.Note that these vectors are time-like u a ± u ± a = −1 in there respective regions.The two regions are modelled as oriented Lorentzian manifolds denoted as M ± = M ± ∪ ∂M ± .The boundary lies on the hypersurface Σ of both spacetimes Σ ∈ ∂M ± .We consider a situation where both spacetimes can be combined into an ambient spacetime (M, g), whose manifold is the union of the manifolds of the individual parts such that This setup is better understood using the language of the distribution function.The Heaside function H(t) is used to constrain evolutionary history in both manifolds.The Heaside function is normalised such that it is equal to +1 if t > 0, 0 if t < 0 and indeterminate if t = 0, with the following properties where δ(t) is the Dirac distribution.The metric of the ambient spacetime is related to metrics in M + and M − as where last term δg ab denotes the metric at the boundary.For the smooth joining of the metrics at the boundary, we require that δg ab = 0 vanishes and the metrics join smoothly at the boundary [77] In order to reduce clutter, we drop H in the subsequent discussion.We define the action of the curve between point p to point q in the ambient spaceetime as a sum of the actions of the smooth curves in the two manifolds where L − is the Lagrangian for the smooth curves in the collapsing region and L + is the Lagrangian for smooth curves in the expanding region.Prime indicates a derivative with respect to the argument.For a smooth curve γ − is geodesic within the range [p, t dp ] and γ + is geodesic within the range [t dp , q].The critical point of the total action with respect to the infinitesimal variation as described in sub-section II.2 corresponds to Then varying both actions following the steps described in section II.2 without imposing the proper variation condition to zero gives Now imposing proper variation at the endpoints of the geodesic ξ i − (p) = ξ i + (q) and not at the boundary gives Given the smoothness condition for the coordinates given in equation ( 44), we impose that the curves are piecewise smooth, that is d dt Inserting the Lagrangian introduced in equation ( 5) into equations ( 55) and ( 55) gives the corresponding geodesic equations in both spacetime regions.And at the boundary, we have the following condition that must hold It will become clearer shortly that equation ( 57) is the generalised Israel junction conditions [77].The boundary conditions for the 4-vector for the geodesic equation are obtained by plugging in tthe geodesic equation given in equation ( 5) The geodesic deviation equation associated with the piecewise geodesics can be obtained by performing a second variation of equation (49) as discussed in sub-section (II.3).However, since equations ( 55) and ( 56) are simply Euler-Lagrange equations of motion we can obtain the respective geodesic deviation equation using the Lagrangian for the geodesic deviation equation given in equation (19).Furthermore, putting the Lagrangian for the geodesic deviation equation given in equation (19) in the generalised Junction condition gives the junction condition for the second fundamental form where we made use of equation (20) in the third equality.Equations ( 58) and ( 59) define the junction conditions that allow glueing spacetimes together across a boundary hypersurface Σ = ∂M + ∩ ∂M − via a thin shell.Note that equation ( 59) is the second Israel Junction condition [77].

III.2. Gravity and the validity of fluid approximation
We showed in the previous sub-section that the general conditions for joining two families of geodesics at a given boundary.We derived the junction conditions that the 4-velocity vector and the second fundamental form that the two families of geodesics must satisfy.We have not explicitly made use of any specific theory of gravity.What we have derived so far will apply to any theory that respects the principle of least action.Now, we need to derive corresponding conditions for a given theory of gravity.Here, we consider the Einstein general theory of relativity, it relates the geometry of spacetime to the matter content of the universe: where T ab is the energy-momentum tensor and R is the Ricci scalar.We have so far derived the generalised junction conditions and equation of motion for the trajectory of massive particles on curved spacetime.Equation ( 60) may be obtained from the Einstein-Hilbert action.For oriented manifolds such as those described in sub-section III.1, the Einstein-Hilbert action has non-vanishing boundary terms [78,79].One could vary the Einstein-Hilbert action similarly as we did in equation ( 50) and avoid setting the tangential derivatives to zero at the boundary [? ] to obtain the corresponding equation of motion.A slightly different way to obtain the same result is to recall that Einstein's theory of gravity is a second-order partial differential equation.In this approach, one finds that the momentum constraint component of the Einstein field equations contains the tangential derivative of the extrinsic curvature tensor which does not need to vanish at the boundary.This constraint can easily be derived using the Gauss-Codazzi identity and equation (58).With respect to the Einstein field equation, this is interpreted as stress-energy tensor at the boundary [80] Note that the violation of the [[K ab ]] = 0 implies that the spacetime is not smooth at Σ.This has sound physical interpretation because it indicates that the surface layer has a non-vanishing stress-energy tensor.In cosmology for Gaussian initial conditions, one expects fluctuations of order 10 −5 in the matter density field on large scales.On small scales, the fluctuations are much larger.
The observed universe with a characteristic size L1 is modelled as a fluid that contains N1 fluid elements, i.e the green balls.The green balls are gravitationally bound clusters, they are virialised.The fluid description description of the evolution of the universe breaks down when the interaction length between the fluid elements is of the order of L2.That is the condition for validity of fluid approximation is given by N1 ≫ N2 ≫ 1 and L1 ≫ L2.We can extend this analogy because in cosmology the red balls are not fundamental particles.We consider them as stars, hence we can also describe the dynamics within the green ball using fluid approximation where the red balls are fluid elements.
In physics, the elementary particles of nature are leptons, quarks, and gauge bosons.These particles are quantum mechanical in nature.It is not clear yet how to fit quantum mechanism and general relativity together.So we cannot associate elementary particles to the trajectory of massive particles we have derived so far.In cosmology, however, fundamental interactions between these elementary particles are not important, rather the dynamics of planets, stars, galaxies, clusters, etc are important depending on the length scale of interest.On the Giga-parces scales, for example, one could consider, galaxies as fluid parcels or fluid elements and assign each fluid parcel a geodesic.This is known as a fluid approximation.It is assumed that on Giga-parces scales, the internal dynamics of a galaxy are not important, hence, interactions within it are averaged over.Even in the N-body simulation in cosmology, a similar approximation is made but it is interpreted as mass resolution.In this case, a fluid element which in principle is made up of many dark matter particles is assigned a mass and a geodesic [81].The N-body simulation evolves the fluid element and not individual dark matter particles.
The fluid approximation breaks down when the internal dynamics within the fluid element become important.In our case with clusters as fluid elements, the fluid approximation will break down at the causal limit or the zero-velocity surface.This indicates that internal dynamics of clusters cannot be ignored.Therefore, in order to describe what happens within clusters, we could consider stars as fluid elements and assign a different one-parameter family of geodesics to each state.The fluid approximation will apply but in a different one-parameter family of geodesics.The conditions for the validity of fluid approximation are described in detail in figure III.2.The essential point is that fluid element before shell crossing is different from the fluid element after shell crossing.

III.3. Fluid rest frame and the far-away observer 4-velocity
In the end, the total energy-momentum tensor in the ambient spacetime would include the stress-energy tensor at the boundary of the different fluid approximations where T + ab and T − ab are the energy-momentum tensor in the expanding and collapsing regions.S ab is the stress-energy momentum tensor due to a jump discontinuity in the Riemann tensor.The physical interpretation of S ab is given in terms of the energy-momentum tensor.Within the standard cosmological model, for example, the late universe is dominated by the cosmological constant and dust.In this limit T + ab may be decomposed as where ρm± is the matter density field ρm± = T ±ab ũa ± ũb ± .Note that one can work with perfect fluid or fluid with non-vanishing anisotropic stress.For the stress-energy tensor due to the jump in the Riemann tensor, it can also be decomposed in a similar way [? ] Note that the extrinsic curvature is related to the covariant derivative of where ÃSa = ũb S ∇ b ũSa is the acceleration in the rest frame of u Sa , ΘS is expansion of the geodesics ũSa , σSab and ωSab are the corresponding shear and vorticity respectively.Furthermore, it is instructive to interpret contributions to S ab (i.e equation ( 64)) in a similar fashion as T ab ± .This can be done by decomposing S ab into irreducible units with respect to ũa S and h±ab where ρS , PS , qSa and π⟨ab⟩ S are the corresponding boundary layer energy density, pressure, energy flux vector and stress-energy tensor respectively: Using equation (65) in equation ( 64), we can calculate these observables ρS = 0 , PS = 1 2π As we shall see later, these fluid variables (i.e.PS , qSa and π⟨ab⟩ S ) are generated by the relative motion between adjacent fluid elements in the neighbourhood of the boundary.The relative motion induces internal friction (viscosity) at the boundary even if it were a perfect fluid in the bulk.To capture this effect, we parametrise π ⟨ab⟩ + in terms of the bulk and shear viscosity components π⟨ab⟩ where ξ ± > 0 and η ± > 0 are bulk and shear viscosity respectively.In the second equality, we have approximated π⟨ab⟩ ± with the shear viscosity only for simplicity.It is straightforward, to include the additional contribution due to the bulk viscosity.
The energy-momentum tensor for each of the fluid specie is measured in its rest frame is given by T ab tot = ρm+ ũa where we have replaced S ab with The energy-momentum tensor that goes into the Einstein field equation in the ambient spacetime is the sum of the energy-momentum tensor of the fluid element decomposed in terms of the threading 4-vector.The matter 4-velocity is related to the fundamental(threading time-like) 4-velocity (we refer to this as the 4-velocity of the observer) according to where v a + is the relative velocity between the matter and observer frames and γ + is the Lorentz boost factor.The projection tensor to matter hypersurface is given by g The full covariant decomposition of spacetime covariant derivative of ũa + is given in equation (65).The decomposition of the full spacetime covariant derivative of v a + with respect to u a + is given by Up to the leading order in v a + , the inverse of equation ( 73) is given by u a + = γ ũa + + va + ≈ ũa + + va + , va + = −v a + , where va + ũ+a = 0, and va + v +a = v a + v +a .At the leading order in v a + these observable quantities in both frames are related according to [82] The transformation between the components of the energy-momentum tensor between these frames is given by [49,82] qa In General relativity, the fundamental 4-velocity is curl free, i.e D [b A +c] = 0, therefore from equation ( 25), we have that ω +ab = 0, hence, ω+ab = D [a v +b] .Also, on the homogenous background σ +ab = 0, therefore σ+ab = D ⟨a v +b⟩ .Note that we can obtain a similar set of relations and equations for the collapsing region.

IV. VORTICITY GENERATION AT THE BOUNDARY
In this section, we describe how vorticity is generated at the boundary of the two fluid elements due to the viscosity or relative friction at the boundary.The gradient of the pressure, the gravitational potential and the rate of expansion in the immediate neighbourhood of the boundary play a very important role.Note that the presence of vorticity invalidates the focusing theorem argument.

IV.1. Continuity and Euler equations in a dust dominated oriented universe
Using the Bianchi identity ∇ [a R bc]d c = 0 and contracting it twice leads to ∇ a G ab = 0. Taking the divergence of equation ( 60) leads to the conservation equation for the total energy-momentum tensor, i.e. ∇ b T ab tot = 0.The time and spatial components of ∇ b T ab tot = 0 at leading order in relative velocity is given by ρm+ + (ρ m+ + P + ) [Θ + + Divv + ] = 0 , (79) where Divσ +a = D b σ +ab is the divergence of the rate of shear deformation tensor and D b is a spatial derivative on the corresponding hypersurface.Again, a similar equation exists for the collapsing region.We split the boundary stressenergy tensor into collapsing and expanding parts (see equation ( 64) for details).The motivation for this approach is that there exists physical processes such as diffusion that transfer of energy/information from the boundary to the bulk.This is slightly different from cases where the conservation of the boundary stress-energy tensor is treatment as a separate unit [83].The approach we discuss here is common in the field of hydrodynamics [84] Equation ( 79) is the continuity equation or mater conservation equation and equation ( 80) is the Euler equation.The Ṗ+ comes from the time derivative of equation ( 78) since qa vanishes in the rest-frame of the fluid qa The Euler equation agrees with the results obtained in [82] in the limit of vanishing viscosity.In addition, we have also set terms such as A b + σ +b a , σ +bc σ bc + v a + , to zero since they will only contribute at high order.Going beyond the linear approximation in v a will be straight-forward.
Acting on equation ( 80) with a spatial derivative operator and taking the antisymmetric part gives the vorticity propagation equation The boundary pressure, P + acts like pressure associated with a Barotropic fluid (P + ∝ ρ +m ).This can be seen from the definition of boundary pressure given in equation ( 68): ∝ ρ m+ (Θ + is related to the matter density field through the continuity equation).Thus, it is consistent to neglect its contribution to vorticity from terms such as D [a P + D b] ρ +m /(ρ m+ + P + ) 2 → 0 .Similarly, we can use equation (79) to relate Ṗ+ to the square of sound speed, ∂P + /∂ρ +m = c 2 s+ : + leads to vorticity diffusion term.This can be seen by acting on the Ricci identity with ε abc [85] Taking the spatial derivative of equation ( 84 where ω+a = ε +abc ω[bc] + /2.The second term on the RHS is a vorticity diffusion term.The two terms on the LHS are derivative terms, this becomes clear when the Lie derivative operator is used to bring them together L ũ ω+a = ω+a + Θ + ω+a /3.

IV.2. Vorticity generation and the line of sight
Vorticity is one of the most important physical quantities in fluid dynamics.It is the most crucial observable for weather forecasting in the local environment [86].Vorticity gives a microscopic measure of the rotation at every point in the fluid flow.As mentioned earlier, there is no true source term for vorticity in the vorticity propagation equation for the fluid (i.e equation ( 85)).The vorticity generation mechanism we describe here is built on the boundary layer theory developed by Morten (1984) [84] for incompressible fluids.Morton showed that the source for all vorticity in the fluid flow emanates from the boundary layers and the rate at which vorticity enters the fluid from the boundary is determined by the viscous diffusion term.The rate of viscous diffusion is determined by the conditions imposed on circulation at the boundary layer.
We discussed in sub-section III.2 how vorticity is generated at the boundary between two different fluid elements in Cosmological context.The Morten boundary layer theory fits naturally into our problem because Θ vanishes at the boundary, therefore, the fluids act much like an incompressible fluid in the neighbourhood of the boundary.Our target here is to use this idea to show how vorticity is generated at the boundary between the expanding and collapsing regions of the universe.
Firstly, we would like to show how the radial component of the vorticity is generated.This component is very crucial for observational purposes because, for an observer, the radial component of vorticity could be compared to the Kaiser redshift distortion term [87].To accomplish this, we decompose every observable into 1 + 1 + 2 irreducible units.The first '1' denotes the time direction., the second '1' denotes the radial direction while the '2' denotes a closed 2-surface or the screen space.For the radial direction, we use radial unit vector n a with the following normalisation n a + n +a = 1.n a + is orthogonal to u a + : u a + n +a = 0.The metric on the 2-surface is given by N +a b orthogonal to n a + and u a + : n a + N +ab = 0 = u a + N +ab .It is used to project tensors onto 2-surfaces (N +a a = 2).The spatial derivative of n a + can be decomposed into irreducible form [88]: where ϕ + denotes the expansion/contraction of the 2-surface, ζ +ab is the shear distortion of the 2-surface, ξ + denotes the twisting of the 2-surface and na + = β b + is the acceleration of n a + .The propagation equations for these terms are given in [88].We do not need their explicit form for the discussion that will follow.Furthermore, any 3-vector can be decomposed into radial components and tangential components.For example, we decompose the vorticity further as ωa + = ω+∥ n a + + ωa +⊥ , where ω+∥ ≡ ω+a n a + , and ωa The acceleration vector A a is decomposed accordingly.The shear tensor on the spatial hypersurface is decomposed as where we have introduced tensors that transforms as scalar, vector and tensor on the 2-surface

Rate of change of scalar circulation
The vorticity is related to the circulation through the Stokes theorem where dS a is an oriented sheet orthogonal to the spacelike vector n a : dS a = n a dS.The surface S is bounded by a closed contour dl a and t a is a 2-vector on the sheet.The circulation defined in equation ( 92) is more appropriately known as scalar circulation since it probes only the component of the vorticity parallel to n a .The total circulation within the universe with two different regions as shown in figure 3 is given by 3. Consider a projected sheet of two fluids in the region S = −, + within the sheet AS .The sheet is bounded by two curves parallel to the interface between the collapsing and expanding regions.The circulation is calculated in the limit where the two curves approach the interface between the two regions/sheets.
where Γ − and Γ + is the circulation in A − and A + respectively.
B A γ ds is the circulation contained in the interface region and ta is the density of circulation contained within the boundary region due to a junction in the velocity vector.A and B are the limits of the integrals.The full spacetime decomposition of the covariant derivative of u a in 1 + 3 is given in equation ( 21) and in 1 + 1 + 2 is given by The time derivative of circulation in each time slice can be performed with the help of the fundamental theorem of calculus: where A = A a dx a is a 1-form and its Lie derivative is given by L Note that u a is the threading 4-vector.We compute the time rate of change of circulation along the threading 4-vector where the Lie derivative of vorticity is given in equation ( 85) and dS a = n a dS.To perform the integral over a closed 2-sphere, we need to further decompose ωa into components living on a closed 2-sphere with one component parallel to n a .Further decomposition of the spatial Laplacian of ωa is given by Note that ϕ is just an inverse of the fixed radial distance.We require that both the first and second radial derivatives of ωa on closed 2-sphere vanish, thereby leading to a much-simplified expression Now we are in a position to evaluate the integrals in equation (97).Note that the total circulation around the whole loop is the sum of the circulations around the two loops: where the circulation around Γ − is the sum of an integral along Γ A and along Γ BA .Similarly, the circulation around Γ + is the sum of two parts, one along Γ B and the other along Γ AB .The integral along Γ AB differs from Γ BA by an opposite sign.This is because the direction of travel is opposite.The radial component of the vorticity propagation equation becomes Note that ϕ vanishes at the boundary since Θ vanishes.Putting equation (101) in equation ( 97) while remembering to use equation ( 100), we can perform the integration using equation ( 92) and Divergence theorem gives [89]: where Σ − = η − t a ∇ ⊥a ω−∥ and Σ + = −η + t a ∇ ⊥a ω+∥ is the diffusive vorticity flux of vorticity from the fluid in − and + respectively.To evaluate the rate of change of the γ, we have to, first of all, separate the terms We can make use of equation ( 95) to evaluate the time rate of change of the integral.Note that to do this consistently we will need a 1 + 1 + 1 + 1 decomposition equivalent of equation ( 95).However, since we are considering only the linear order approximation, we can proceed as in equation ( 97) The circulation within the interface is generated by the relative acceleration of fluid elements on both sides of the boundary.The Lie derivative of closed 2-sphere projected relative velocity (i.e L u v + a ), is given by the momentum conservation equation where L u v a ⊥ = va ⊥ + Θv a ⊥ /3 is a coordinate independent acceleration term.Note that at leading order and assuming the FLRW background ṅa ≈ 0. Contracting equation ( 105) with a 2-vector t a gives where ε ab ∇ b ⊥ ω∥ ≡ ∇ ⊥a ω∥ and at leading order A ⊥a = ∇ ⊥a Φ.We introduce another directional directive t a ∇ ⊥a = ∂/∂s to improve clarity.Identifying the last term in equation ( 106) with the diffusive flux introduced in equation ( 102) and making both terms the subject of the expression leads to Combing equations ( 107) and (108) give where the boundary vorticity flux is related to the boundary acceleration through the momentum equation.This provides a source for the vorticity propagation equation (101).We can write the jump in the viscous acceleration of fluid elements on each side of the closed 2-sphere as The jump in viscous acceleration is sourced by the angular gradients in the gravitational potential, pressure and expansion.Plugging equation (110) into equation ( 104) and performing the arc-length integration gives The circulation within the boundary is generated by the relative acceleration of fluid elements on both sides of the boundary which is sourced by the differences in gravitational potential, the sum of pressure and matter density and relative expansion.The last two terms on the second line vanishes at the boundary since Θ = 0. Neglecting these terms, we recover the standard result in boundary layer vorticity generation [89,90] In this limit, circulation within the interface is generated by the differences in gravitational potential and the ratio of pressure to the matter density.Although there is a flux of vorticity term in equation ( 112), it does not contribute to the net rate of circulation The flux diffusion term transport vorticity from the boundary to the interior of the fluid, they do not generate vorticity.The vorticity is generated by the difference in the gravitational and pressure across the boundary.Moreover, if one imposes the smoothness of the gravitational potential at the boundary [[Φ B ]] = [[Φ A ]] = 0, then the vorticity can only be generated by the differences in pressure across the interface, otherwise the circulation is globally conserved.
A mutual generation and annihilation of vorticity could still happen in the neighbourhood of the boundary.Finally, if we further impose no-slip condition at the boundary, we find a simple expression for the vorticity flux across the interface Here the vorticity flux is sourced exclusively by the pressure, gravitational potential and expansion gradients, this is in agreement with [90,91].

IV.3. Vorticity generation at the boundary and the rate of change of vector circulation
In this sub-section, we extend the treatment to the vorticity vector on the hypersurface.The basic building block of this is the curl theorem.It relates vorticity at a given instant in time to the relative velocity [89] where Γ a is called a vector circulation, it gives the macroscopic picture of rotation within a given local region and V is the volume of a given hypersurface, ∂V denotes the boundary of the hypersurface, The total vector circulation of the universe at a given time slice is given by the sum of vorticity in the fluids in the two regions separated by a shell and the vorticity within the shell [89] Γ Note that n a is pointing in the direction of increasing radial distance.The time derivative of the vector circulation is given by where − is the density of circulation contained within the boundary region due to a possible jump in the velocity vector.Using equation (95), we evaluate the time derivative over the integral and substituting for the Lie derivative of ωa using equation ( 82) leads to Performing the integration in equation ( 117) while taking into account the possible discontinuity across the shell gives where we have introduced the viscous vorticity flux and the acoustic flux in the neighbourhood of the shell For the integrals containing the viscous terms in equation ( 119), we made use of the Divergence theorem to perform the integration.To evaluate the rate of change of the vector circulationl, we have to first of all separate out the terms for clarity The directional derivative of ε ab along u a is given by where A c ⊥ is a closed 2-sphere projected acceleration vector.The directional derivative of n a along u a is given by [88] ṅa = A ∥ ua + α a where α a ≡ ṅā and A = n a ua . (124) These terms vanish on an FLRW background, similarly, εab and σ c [a ε b]c vanish, hence the Lie derivative of γ a at linear order in relative velocity is given by Using the momentum constraint equation, we find that the Lie derivative of v a is given by Note that we can now define the projected angular derivatives and projected tensors The last two terms in equation ( 126) can be identified with the terms defined in equations ( 120) and (120).Making these terms the subject of the formulae leads to where we take the following approximation ε ab L u v b = L u v ⊥a .Combing these two equations leads to the expression of the diffusive flux vector in terms of the acceleration term, angular pressure gradient, gravitational potential gradients and the gradient in the expansion rates Equation ( 130) indicates that the effective vorticity flux out of a given shell is equal to the difference in acceleration of fluid elements on both sides of the shell.Imposing the no-slip Junction conditions for the geodesic (i.e equation ( 58) The gradients in pressure, gravitational field and expansion across the shell act as sources for vorticity in the universe.The momentum constraint equation at shell crossing is given by Putting equation (132) into equation ( 125) and performing the angular integration gives where we set P = wρ m with w = constant for simplicity.Substituting equation (133) into equation ( 119) and performing further algebraic simplification leads to The circulation is generated by the relative acceleration between fluid elements on each side of the boundary.For [[Θ]] = 0, equation (134) reduces to Again the difference in the integral of the gravitational potential and pressure at the boundary of the thin shell are responsible for the generation of circulation in an initially irrotational fluid.This is caused by the relative acceleration due to the viscosity between fluid elements on each side of the boundary.The viscous forces on their own do not generate circulation on the boundary, rather their role is to transfer the circulation between the boundary and the fluid interior.Finally, the integral of the gravitational potential and pressure at the boundary generates circulation, vorticity is obtained from circulation.

V. IMPLICATIONS FOR THE STANDARD MODEL OF COSMOLOGY
Although the stated aim of this project was to describe a possible mechanism for the generation of vorticity in the universe.Our approach to the problem builds on the knowledge from the N-body simulation where vorticity generation is linked to shell crossing singularity [15,92].Vorticity generation and shell-crossing singularity in cosmology are complex topics on their own, putting them together to come up with a unified picture of vorticity generation is a much more complicated task.We have been able to put this together into a consistent picture.The key result we present here is that vorticity generation at the boundary of the causal horizon is essential for avoiding the shellcrossing singularity that appears when a given coordinate is extended beyond the causal horizon.In addition to this key result, we would like to discuss some obvious features and predictions that emerged from this model.

V.1. Existence of a causal limit for massive particles
The key feature of the universe model described here is the existence of a "causal horizon" (boundary) for time-like geodesics.The causal horizon divides the universe into two regions: expanding and collapsing regions.Using the principle of least action we derive the equations of motion and boundary conditions applicable in each case.With these equations and boundary conditions, we showed how the vorticity is sourced by the gradients in gravitational potential and pressure across the boundary.The vorticity flux and acoustic flux generated at the boundary then diffuse away from the boundary.This is a unique model of the universe that does not introduce any new free parameter, rather it makes a prediction about the existence of a causal horizon for massive particles in an expanding universe.
In general, causal horizons are determined by the dynamics of time-like and light-like geodesics [42].Within the standard model of cosmology, the well-known causal horizon is the particle horizon.It is determined by the dynamics of the light-like geodesic, it is important to note that the particle horizon indicates the maximum distance light from particles could have travelled to the observer in the age of the universe.Ellis and Stoeger had earlier argued that there must exist a causal horizon in our universe that is determined by the dynamics of time-like geodesics [41].This causal horizon is a unique feature of general relativity.We showed how the causal limit could be determined given a halo model for a gravitationally bound cluster.Regions within the causal horizon are collapsing while regions outside are expanding.Ellis and Stoeger argued that only the collapsing region contributes most significantly to the dynamics of our local environment and that the dominant interaction within this region is not mediated by massless particles, instead, they are mediated by massive particles that travel at very low speeds relative to the cosmic rest frame.The vector and tensor perturbations on an FLRW background spacetime have negligible impact on the dynamics collapsed region.It is the differences in speed that cause the gravitationally bound cluster to decouple from the Hubble expansion since it cannot keep up with the speed of the expanding cosmic rest frame [41,93].The dynamics of null geodesics in the presence of the time-like causal limit will be discussed elsewhere.

V.2. Generation of of the magnetic part of the Weyl tensor
One other unique feature of the model is the emphasise on the distinction between the frames of reference of the observer and the fluids or matter fields.The observer frame is tilted with respect to the fluid frame according to ũa ≈ u a + v a , where u a is the observer four velocity, ũa is the fluid four-velocity and v a is the relative velocity between them.The dynamics of the fluid element in the expanding region are different from the dynamics of the fluid element in the collapsing region.The friction between the fluid elements at the boundary creates the necessary condition for the generation of vorticity at the boundary.Time and ruler are determined by the observer with u a .The vorticity we describe is associated with ũa (see equation ( 77)).This is crucial because we can make use of the Ricci identity 2∇ [c ∇ d] ũa = R e dca ũe to relate the fluid vorticity to the existence of the magnetic part of the Weyl tensor.The Weyl tensor is the trace-free part of the Riemann tensor, it describes the curvature of spacetime.The fluctuations in the magnetic part of the Weyl tensor have been associated with the generation of gravitational waves [94].To see how the existence of vorticity is connected to existence of the magnetic part of the Weyl tensor, we project the two free indices of the Ricci identity with ε cdb .The further simplification of the results gives For only scalar perturbations ε acd Dc σd b = 0, then equation (136) reduces to Hab + 2 D⟨a ωb⟩ = 0.This implies that for non-zero vorticity in the fluid, the magnetic part of the Weyl tensor is non-vanishing.This has some consequences: • Gravitational wave can be generated in the limit of the scalar perturbations if the vorticity generated at the boundary is non-zero.
• The non-zero magnetic part of the Weyl tensor could have a consequence for the theory of dark matter.This was discussed in detail in [25].
• The nonzero magnetic part of the Weyl tensor could also have implications for the generation and propagation of Maxwell's magnetic field in clusters [85].
Finally, non-zero vorticity implies the existence of coherent helicity in the fluid interior.Helicity is usually defined as the scalar product of the fluid velocity and the vorticity vector at each point in the fluid.The helicity measures the extent to which a flow field carries vorticity in a specific direction.

V.3. Non-local variance in Hubble rate contributing to acceleration in the universe
In section III, we introduced a model of the universe that treats the expanding and collapsing regions of the universe consistently with suitable junction conditions.One unique property of this model is that both regions do not exchange mass rather the expanding regions with modes greater than the size of the collapsing region simply rescales the size of the collapsing regions [58].In terms of the energy-momentum tensor, we showed in sub-section III.2 that the jump in the Riemann tensor manifests as viscosity in the fluid element at the boundary of the collapsing and expanding regions.In this sub-section, we argue that an observer associated will infer an accelerated expanding universe even if the rate of expansion of the expanding and collapsing regions of the universe are decelerating.This can be seen by calculating the volume average of expansion, i.e Θ.To see how this works, we start with the definition of an integral of a spacetime scalar S in a manifold with an ambient metric g ab [95,96] where W M4 is window function that selects slicing and foliation hypersurface and x b is the adapted coordrinate system.
With respect to the model described in section III, (M 4 , g) corresponds to the ambient spacetime, which is a union of the manifolds M = M − ∪ M + describing the expanding and collapsing regions respectively.A full treatment of this system will include a consistent implementation of the boundary conditions at the level of the volume integration, however, in sub-section III.2 we consider a limit where the effect of the boundary conditions is treated as effective fluid.This is equivalent to decomposing equation (138) into disjoint set: The average of a scalar on an arbitrary manifold is defined as where . With respect to equation (139), the average of a scalar decomposes where and ⟨S⟩ W M − are average of a scalar defined on M + and M − respectively.
The window function selects a slice and foliation W A0,A,B0,B = (V a ∇ a A) δ (D) (A 0 − A) (B 0 − B) , where A defines the foliation and B defines the radial extent [97].We focus on the foliation defined by the 4-velocity of the fluid (ũ a ) introduced in equation ( 73) as seen by a tilted observer with 4-velocity u a .The effective scale factor in the expanding and collapsing regions is defined as where h is the metric on the hypersurface of the fluid velocity(i.e equation ( 74)).The Hubble rate in both regions is given by The acceleration of the ambient spacetime is then given by A full treatment of this involves the use of the Einstein field equations and spatial averaging techniques [98][99][100], however, the point here is that the model we describe does not only describe a mechanism for the generation of vorticity in the universe, it also gives a natural explanation for accelerated expansion in the universe without the need for any assumption about the energy content of the universe.Here, the acceleration of the ambient spacetime volume, d 2 a W Σ 3 /dt 2 > 0, is easily realised for known standard matter source, the accelerated expansion is simply due to the fact that Hubble rate in the collapsing region has a sign opposite to the Hubble rate in the expanding region H W Σ − ∝ −H W Σ + .

VI. CONCLUSION
Our understanding of the evolution of large-scale structures in the universe is built on an expanding FLRW background spacetime.The general relativistic 1 + 3 covariant decomposition provides tools for studying evolution of spacetimes by simply looking at how one-parameter family of geodesics propagates.One key property of geodesics in general relativity is that they can cease to be geodesics in a finite time or affine parameter [101].The geodesic, the path a particle takes in a gravitational field is determined by the curvature of spacetime, which in turn is determined by the distribution of matter and energy.The changes in the nature of matter and energy distribution are reflected on the propagation of geodesics.The study of validity range of geodesics have received very little attention in cosmology even when geodesics constitute the bedrock for studying the growth of structures using the N-body simulations [26,27,102].Geodesics are studied more diligently in the field of differential geometry, especially in the blackholes physics where it is well understood that presence of horizons is inevitable in general relativity [42,103].
In this paper, we have studied the validity range of time-like geodesics in a universe like ours which has overdense(gravitationally bound region, not undergoing Hubble flow) and under-dense regions (undergoing Hubble flow).We show that geodesics on expanding region of the universe cannot be be extended to the over-dense region if singularities are to be avoided.We show that a causal horizon forms at the zero-velocity surface which then serve as a boundary separating families of one-parameter family of geodesics that are causally disconnected.This distinction was pointed earlier by Ellis and Stoeger [41].We argue that the fact that the determinant is finite at the causal horizon allows to define consistent junction conditions for the coordinates, metrics and extrinsic curvature tensors at the boundary.We derived a generalised boundary condition starting from the principle of least action.The generalised boundary condition we found reduces naturally to the Israel junctions conditions [77].
The physical picture of what we have described mathematically can be easily be visualised as follows.The existence of a causal horizon for time-like geodesics divides the universe into expanding and collapsing regions.The over-dense region such as the gravitationally bound systems like clusters are decoupled from the Hubble expansion.Dynamics within the gravitational bound regions are more crucial for the formation of sub-structures such as galaxies and they are causally disconnected from the dynamics in the expanding region.The causal horizon exists because massive particles within the over-dense region are moving with a slower velocity such that they cannot catch up with the Hubble flow [41].
One of the key feature of the scenario we described is that it provides a mechanism for the generation of vorticity in the universe.The jump in the Riemann tensor at the boundary of both regions could be interpreted in terms of the stress-energy tensor.We show that the stress-energy tensor at the boundary has a natural physical interpretation in terms of the effective energy-momentum tensor and therefore, can be decomposed with respect to the fundamental four vector which then allows to apply Morton's boundary layer theory for the generation of vorticity [84].The components of the boundary effective energy-momentum tensor include a non-vanishing boundary anisotropic stress tensor, pressure and the momentum flux vector.These observables exist even if the matter in the universe is purely dust.We then showed that the non-vanishing contribution of these observables leads to a non-vanishing scalar and vector circulation sourced by the jump in gradients of gravitational potential, pressure and expansion.
The results of this study is very crucial towards understanding the local universe because the vorticity has been observed around filaments, clusters, etc [11].Also, some high-resolution N-body simulations have been able to detect vorticity in the outskirts of the large-scale structures but the mechanism for its formation has eluded a clear analytic understanding [15,92].Although there are earlier works based on the N-body simulations that established a link between vorticity generation and caustics formation or shell-crossing singularity in cosmology [20,104], This is the first paper to describe in detail how vorticity could be generated at the boundary layer between the over-dense and under-dense regions in the universe in line with the Morton boundary layer theory [84].Furthermore, we showed that the vorticity flux is sourced by the jump in gradients of the gravitational potential, pressure and rate of expansion.The vorticity flux is generated at the boundary, then gradually diffuses towards the fluid interior.
Finally, we discussed other obvious predictions/extensions of our work in section V. Some of these predictions/extensions include the existence of a causal horizon for massive particles in the universe.The implications for the existence of causal horizon in relation to the possible role the magnetic part of the Weyl tensor could play in the local universe.There is also a possibility that non-zero vorticity could act as a source for the magnetic fields in clusters.We also discussed how the non-local variance in the Hubble rate between the collapsing and expanding regions could manifest as an apparent accelerated expansion of the ambient spacetime.A comprehensive discussion of these and many more would be provided elsewhere.
The perturbation theory computations in this paper were done with the help of tensor algebra software xLightcone and xPand [105], These apps are based on xPert and xTensor [106].I made use of COLOSSUS (A python toolkit for cosmology, large-scale structure, and dark matter halos) developed by Benedikt Diemer for computations involving the dark matter halo profiles [72].
of massive particles on curved spacetimes II.1.Fluid flow in Newtonian gravity II.2.Geodesic of relativistic massive particle and its evolution equation II.3.Validity of geodesics, focusing theorem and horizon II.4.Inevitability of more than one-parameter family of curves in a universe like ours II.5.Caustics and inverse function theorem III.Model of the universe with collapsing and expanding regions III.1.Dynamics of geodesics in both regions and their junction conditions III.2.Gravity and the validity of fluid approximation III.3.Fluid rest frame and the far-away observer 4-velocity IV.Vorticity generation at the boundary IV.1.Continuity and Euler equations in a dust dominated oriented universe IV.2.Vorticity generation and the line of sight IV.2.1.Rate of change of scalar circulation IV.3.Vorticity generation at the boundary and the rate of change of vector circulation I. INTRODUCTION

FIG. 1 .
FIG. 1.In the left panel, we show the plot of d ln ρ/d ln r vs the comoving radius.The position of the sharpest drop in density indicates the location of the halo boundary otherwise known as the splashback radius.On the right panel, we show the expansion as a function of the comoving radius.The thick horizontal line corresponds to Θ = 0. We considered the following halo masses {M1, M2, M3, M4} = 1 × 10 11 , 1 × 10 12 , 1 × 10 13 , 1 × 10 14 M⊗ and fixed the halo concentration to cvir = 7.M⊗is the mass of the sun.Note that the causal horizon is much greater than the splashback radius.Both of these radii can in principle be measured.

a ũ+ b and at the leading order in v a + is given by
) However, [[Θ]] vanishes at the boundary, therefore, we set the contribution of Ṗ+ to zero.Using the Ricci identity for the matter 4-vector: 2∇ [c ∇ b] ũ+a = R d abc ũ+d , we show that η + D [a Divσ b] ) and extracting the anti-symmetric part gives ε abc D b Divσ c + = 2D b D b ω+a .Putting all these back into equation (82) gives the propagation equation for the vorticity ω+a + 1 3 Θ+ ω+a = η + D b D b ω+a .( h ab + Hab + 4 Ã(a ωb) + 2 D⟨a ωb⟩ = 0 (136) Similarly, contracting the Ricci identity with ε cda gives another constraint equation for the divergence of the vorticity vector Dc ωc = Ãb ωb .In the rest frame of matter fields, Ãb = 0, therefore D c ωc vanishes, hence 2ε acd Dc σd b + Hab + 2 D⟨a ωb⟩ = 0 .(137) at the boundary, therefore, the Euler Lagrange equation holds in the separate spacetimes d dt