Spherical and nonspherical models of primordial black hole formation: exact solutions

We construct spacetimes which provide spherical and nonspherical models of black hole formation in the flat Friedmann--Lemaitre--Robertson--Walker (FLRW) universe with the Lemaitre--Tolman--Bondi solution and the Szekeres quasispherical solution, respectively. These dust solutions may contain both shell-crossing and shell-focusing naked singularities. These singularities can be physically regarded as the breakdown of dust description, where strong pressure gradient force plays a role. We adopt the simultaneous big bang condition to extract a growing mode of adiabatic perturbation in the flat FLRW universe. If the density perturbation has a sufficiently homogeneous central region and a sufficiently sharp transition to the background FLRW universe, its central shell-focusing singularity is globally covered. If the density concentration is sufficiently large, no shell-crossing singularity appears and a black hole is formed. If the density concentration is not sufficiently large, a shell-crossing singularity appears. In this case, a large dipole moment significantly advances shell-crossing singularities and they tend to appear before the black hole formation. In contrast, a shell-crossing singularity unavoidably appears in the spherical and nonspherical evolution of cosmological voids. The present analysis is general and applicable to cosmological nonlinear structure formation described by these dust solutions.

We construct spacetimes which provide spherical and nonspherical models of black hole formation in the flat Friedmann-Lemaitre-Robertson-Walker (FLRW) universe with the Lemaitre-Tolman-Bondi solution and the Szekeres quasispherical solution, respectively. These dust solutions may contain both shell-crossing and shell-focusing naked singularities. These singularities can be physically regarded as the breakdown of dust description, where strong pressure gradient force plays a role. We adopt the simultaneous big bang condition to extract a growing mode of adiabatic perturbation in the flat FLRW universe. If the density perturbation has a sufficiently homogeneous central region and a sufficiently sharp transition to the background FLRW universe, its central shell-focusing singularity is globally covered. If the density concentration is sufficiently large, no shell-crossing singularity appears and a black hole is formed. If the density concentration is not sufficiently large, a shell-crossing singularity appears. In this case, a large dipole moment significantly advances shell-crossing singularities and they tend to appear before the black hole formation. In contrast, a shell-crossing singularity unavoidably appears in the spherical and nonspherical evolution of cosmological voids. The present analysis is general and applicable to cosmological nonlinear structure formation described by these dust solutions.

Introduction
Black holes may have formed in the early Universe. These black holes are called primordial black holes (PBHs) [1][2][3]. Since PBHs can act as sources of Hawking radiation, strong gravitational fields and gravitational radiation, current observations can place a stringent upper limit on the abundance of PBHs [4]. Generally speaking, PBHs of mass M were formed at an epoch when the total mass enclosed within the Hubble radius was of order M , if we neglect critical behavior in gravitational collapse, mass loss due to Hawking evaporation, and mass gain due to accretion. The abundance of PBHs of mass M depends primarily on the primordial fluctuations of mass M and the state of the Universe at the formation epoch.
Since the Hawking radiation and its final outcome depend on quantum gravity, PBHs can also provide observable phenomena of quantum gravity, even if no positive signal has been detected. Thus, PBHs can be seen as a probe into the early Universe, high energy physics, and quantum gravity. In particular, since PBHs can convey information about primordial A singularity which is not behind a black hole horizon, called a naked singularity, may appear in the course of gravitational collapse with regular initial data in the LTB solution [29][30][31][32][33][34]. This is also the case in the Szekeres solution [35][36][37][38]. Clearly this departure from spherical symmetry cannot avoid the formation of naked singularities. Goncalves [39] studied the cosmic censorship and the curvature strength for the shell-focusing singularities in the quasispherical Szekeres solutions with dust and a cosmological constant based on the radial null geodesics, which do not exist in general. Hellaby and Krasinski [40] studied the condition for the avoidance of shell-crossing singularities in the Szekeres solutions for the quasispherical, quasi-pseudospherical and quasiplanar cases in the very general formulations. The same authors [41] presented the geometrical interpretation of the Szekeres solutions for the quasi-pseudospherical and quasiplanar cases. Krasinski and Bolejko [42] defined an absolute apparent horizon in the quasispherical Szekeres solutions, discussed its difference from an apparent horizon, which is commonly used, and concluded that the apparent horizon can be regarded as the true horizon. Vrba and Svitek [43] rewrote the condition for the occurrence of shell-crossing singularities in terms of the maximum, minimum, and average density of the shell at the moment of occurrence. They also discussed the time evolution of the solutions, restricting themselves to the marginally bound case which is not relevant to the cosmological growing perturbation.
It is important to see the condition for the occurrence and non-occurrence of naked singularities and its physical interpretation in the context of the nonlinear evolution of cosmological primordial fluctuations. We adopt the simultaneous big bang condition to extract a growing mode of adiabatic perturbations. The formation and evolution of cosmological nonlinear perturbations have been analyzed on the same grounds [44,45]. For these nonlinear cosmological perturbations, we investigate the formation of black holes and cosmological voids against the formation of shell-focusing and shell-crossing singularities and see how the deviation from spherical symmetry affects these physical situations for the first time. We also discuss the link between the current result and the physical interpretation by Khlopov and Polnarev [12] that shell-focusing naked singularities may be physically regarded as the onset of strong pressure gradient force. This paper is organized as follows. In Sect. 2, we present and interpret the Szekeres solution. In Sect. 3, we describe the dynamics of the Szekeres solution. In Sect. 4, we analyze the shell-crossing and shell-focusing singularities under the simultaneous big bang condition. There, we find a necessary and sufficient condition for the Szekeres solution to describe the spherical and nonspherical formation of PBHs without suffering naked singularity formation. In Sect. 5, we propose several concrete models that describe the spherical and nonspherical formation of PBHs without naked singularities. Section 6 is devoted to our conclusions. In Appendix A, we briefly review the derivation of the Szekeres solution in order for the paper to be self-contained. We use geometrized units, where c = G = 1.

Presentation of the metric
We present here a brief overview of Szekeres's quasispherical dust solution. For completeness, the derivation of this solution is given in Appendix A. The line element is given by r, x, y), and φ = φ(t, r). The prime denotes the ordinary derivative with respect to the argument and hence is only used for functions of one variable such as A(r), B 1 (r), B 2 (r), C(r), W (r), and S(r). For functions of more than one variable, the partial derivatives are denoted by a comma followed by the index of the differentiating variable. Equation (6) can be integrated to give where H(r) is an arbitrary function of r. Thus, the solution contains seven arbitrary functions of r. With one scaling freedom and one algebraic constraint taken into account, the solution contains five arbitrary functions. Equation (5) implies that neither A(r) nor C(r) can vanish, while B 1 (r) or B 2 (r) can. We assume A(r) > 0 and C(r) > 0 without loss of generality. Equation (4) can be transformed into the following form: which implies P ≥ 1/[4A(r)] and hence P too is positive definite. It is clear that the solution reduces to the LTB solution if A(r) = C(r) = 1/2 and B 1 (r) = B 2 (r) = 0. See Refs. [46,47] for the LTB solution, where W 2 (r) − 1, S(r), and H(r) correspond to the energy, mass, and time functions, respectively, and φ(t, r) gives the areal radius of the two-sphere of constant t and r. The coordinates (x, y) on the two-sphere can be understood in terms of the stereographic projection. One of the three arbitrary functions corresponds to the gauge freedom. For example, H(r) can be fixed if we fix the radial coordinate at t = t 0 so that φ(t 0 , r) = r. The simplest choice to recover the Minkowski spacetime is S(r) = 0, W (r) = 1, and H(r) = 0.

4/28
Equation (8) describing the evolution of shells can be integrated explicitly. For W (r) = 1, we have For 0 < W (r) < 1, we have For W (r) > 1, we have The above solutions can be summarized into the following compact form: where and where Note that for 0 < W (r) < 1 the expanding and collapsing branches are combined into a single complete solution with both big bang and big crunch. Since we are interested in the cosmological solutions, we focus on the branches which possess a big bang, so that G = G + for W (r) ≥ 1, while G = G ± for 0 < W (r) < 1. Figure 1 shows the relevant branches of G for −∞ < Y ≤ 1. We note that G − (Y ) for 0 < W (r) < 1 admits the following expansion: q are defined by

Comoving coordinates and the dipole moment
and then from (p, q) to (θ, ϕ) based on the stereographic projection, where p = cot θ 2 cos ϕ and q = cot The function P can now be expressed as and the metric induced on the two-surface S t,r is given by the standard form of the two-sphere whereas the line element in the spacetime contains off-diagonal terms drdθ and drdϕ in general. Thus, we establish the following interesting relation: where M is the proper mass andρ is given byρ := ρe λ φ 2 /r 2 and is calculated to givẽ Therefore, the mass contained within the volume element dV = r 2 sin θdrdθdϕ is constant in time. This means the coordinates (r, θ, φ) play the role of the comoving spherical coordinates andρ can be interpreted as the conserved mass density. We have seen thatρ gives the conserved mass density in terms of the comoving coordinates (r, θ, ϕ) or (x,ỹ,z). This suggests that withρ we can define a conserved dipole moment. In the comoving spherical coordinates, P can be transformed to the standard form given in Eq. (21). 6/28 Thus, defining vectors β and n as and n := (sin θ cos ϕ, sin θ sin ϕ, cos θ) = 1 r (x,ỹ,z), (26) respectively, we find P ,1 P = − n · β r (27) in the coordinate basis of the comoving Cartesian coordinates (x,ỹ,z). Since the nonspherical dependence of the density distributionρ, as well as ρ, appears only through this combination, the matter distribution has monopole and dipole moments only and the vector β is proportional to the dipole moment. The absolute value β(r) of the vector β(r) is given by It can be easily shown that β 2 (r) is positive definite. The definition of β used here is the same as in Szekeres [35]. We can rewrite the expression for the energy density, Eq. (24), in the following form: whereρ Thus, we can naturally define the nondimensional mass dipole moment localized in the shell labeled r with d.
On the other hand, the physical density ρ can also be written in the following form: where ρ 0 := S 8πφ 2 φ ,1 and χ := We can interpret ρ 0 as the spherical part of the density field, which is identical to that of the reference LTB solution, and χ as the deviation from it. The nonsphericity χ does depend on time through φ ,1 /φ. Although χ cannot be simply interpreted as the deviation due to the dipole moment, it is closely related to the vector β. If β = 0, then χ vanishes identically. Conversely, if χ vanishes identically, then β = 0 or φ ,1 /φ = S /(3S).

Spherically symmetric and axially symmetric subclasses
As is seen in Eq. (7), the energy density depends only on t and r, if is satisfied [48]. From Eqs. (5), Eq. (33) implies that all of A(r), B 1 (r), B 2 (r), and C(r) are constant. Then, the spacetime becomes spherically symmetric and the solution reduces to 7/28 the LTB solution. In this case, we can see β = 0. This supports our interpretation of β as a dipole moment. Equation (34) holds if and only if φ is separable as we can show from Eq. (6). In this case, φ = a(t)S 2/3 (r) and therefore ρ = 3/[8πa 3 (t)], i.e., the density is homogeneous. In this case, since the spatial component of the metric is written as g ij = a(t)γ ij (x i ), γ ij can be shown to be that of the constant curvature space in three dimensions and therefore this is the FLRW solution [22]. The above discussion also implies that the Szekeres solution can be matched to the Schwarzschild solution at any radius r = r m . In fact, we can always choose W (r), S(r), H(r), A(r), B 1 (r), B 2 (r), C(r) to be constant for r > r m . The mass parameter of the Schwarzschild black hole is given by M = S(r m )/2. Then, the region for r > r m is spherically symmetric vacuum and hence the Schwarzschild solution by Birkhoff's theorem.
Next, let us consider a subclass where B 1 (r) = B 2 (r) = 0 identically but A(r) is still allowed to be a function of r. Then, Eqs. (4) and (5) reduce to In the comoving spherical coordinates (r, θ, ϕ), we can easily find Since the metric components in the comoving spherical coordinates do not depend on ϕ, the spacetime is axially symmetric. In the comoving Cartesian coordinates (x,ỹ,z), we find Thus, β is directed along thez axis. This also supports our interpretation of β as a dipole moment.

Shell-crossing singularities
Following Szekeres [35], we quote two lemmas for quadratic forms.
where c is a constant and θ = θ(r).
If we denote the time of the shell-crossing singularity by t = t SC (r, x, y), where t SC depends not only on r but also on x and y in general. Thus, shell-crossing singularities are affected by nonsphericity. If we fix t, a shell-crossing singularity occurs not on the twosurface S t,r but on a different two-surface in general. From Eq. (32), we can find that if P ,1 = 0 and φ ,1 /φ = S /(3S), then χ diverges to infinity at shell-crossing singularities. In 8/28 other words, if the shell is nonspherical, it is the nonspherical rather than the spherical part of the density field that diverges at the shell-crossing singularities. Therefore, generic shell-crossing singularities in the Szekeres spacetime are essentially nonspherical. Applying the lemmas presented above to the function φ ,1 P − φP ,1 , we can show that the two-surface S t,r possesses a shell-crossing singularity if and only if which can be rewritten as directly linking the shell-crossing condition with nonsphericity parameter β. We define t SC (r) as the time of the earliest occurrence of shell-crossing singularity on S t,r so that the quadratic form P φ ,1 − φP ,1 with fixed r begins to have a zero at t = t SC (r). This implies

Regularity condition
We assume the existence of a regular initial data surface (t = t 0 ). Since the areal radius of the two-surface S t,r is φ(t, r), we also assume that φ 0 (r) := φ(t 0 , r) is an increasing function of r and scale r so that φ 0 (r) = O(r). For the center to be locally Minkowski, Eq. (2) implies W (0) = 1. For φ ,0 to be bounded, Eq. (6) implies S(r) = O(r). Then, for ρ 0 (r) := ρ(t 0 , r) to be bounded as r → 0, Eq. (7) implies S (r) = O(r 2 ). This implies We further assume that the metric is C 2− in the comoving Cartesian coordinates. For rP ,1 /P to be C 2− in terms of (x,ỹ,z), β(r) = O(r), and hence the nonspherical functions A(r), B 1 (r), B 2 (r), and C(r) are continuous and differentiable at r = 0.
At t = t 0 , the initial data is shell-crossing free for all r(> 0), x, and y if, and only if In turn, applying Lemma 3.2 to the function S P − 3SP ,1 , we can see that ρ is positive definite at t = t 0 if, and only if Thus, regularity imposes the following condition on β: The above condition is automatically satisfied near r = 0, since β(r) = O(r).

9/28
3.3. Trapped surfaces, and apparent horizons Following Szekeres [35], we consider a trapping condition for a spacelike two-surface S t,r . Since the tangent space of S t,r is spanned by ∂/∂x, and ∂/∂y, any normal vector n µ to S t,r should satisfy n 2 = n 3 = 0. Thus, if we consider the congruence of null geodesics with tangent vector k µ normal to S t,r , we find with k 2 = k 3 = 0 on S t,r . Assuming k 0 > 0, we can identify the null geodesics of k 1 > 0 (< 0) with outgoing (ingoing) ones. We can choose k 1 = 1 (−1) on S t,r by choosing the scaling of the affine parameter. The sign of the expansion coincides with the sign of k µ ;µ , which is calculated to give where we put k 1 = ,, and = ±1. Differentiating Eq. (47) with respect to t, we find In Eq. (48) with µ = 1, we find We eliminate k 1 ,0 from Eq. (50) by Eq. (51),, and eliminate k 0 ,0 + k 1 ,1 from Eq. (49) by the resultant equation. Then, we find Using Eqs. (2), (3), and (6), we can transform Eq. (52) to where ι = sign(φ ,0 ). Equation (44) implies that the first factor is positive because it cannot change the sign without encountering a shell-crossing singularity. Thus, we establish that if the dust is collapsing (ι = −1), the outgoing ( = 1) null normal can have vanishing expansion if S = φ, while if the dust is expanding (ι = 1), it is the ingoing ( = −1) null normal that can have vanishing expansion if S = φ.
Here we identify a marginally trapped two-surface S t,r with an apparent horizon according to Krasinski and Bolejko [42]. Then we call an apparent horizon with vanishing outgoing (ingoing) null expansion a future (past) apparent horizon. Note that this is somewhat different from the notion of an apparent horizon defined by Hawking and Ellis [49]. Although a black hole horizon is usually identified with an event horizon in the asymptotically flat spacetime, such an identification is not so strongly motivated in the cosmological spacetime because of the teleological nature of an event horizon and because the asymptotic condition at null infinity is less physically meaningful in cosmology with a finite particle horizon. Although the local and dynamical definition of black hole horizon is ambiguous, we can identify a future apparent horizon with a local black hole horizon in this paper for not only its physical reasonableness but also its simplicity and usefulness for the analysis.
In the dust models, a noncentral shell-focusing singularity is always future trapped because φ(t, r) = 0 and S(r) > 0 there. This means that noncentral singularity is causally disconnected from a distant observer. 10/28

Six critical events
As for the dynamics of the Szekeres solution, there are six important events: big bang, past apparent horizon, maximum expansion, future apparent horizon, big crunch and shellcrossing singularity. We denote the times of the occurrence of these six events at each shell with t BB (r), t PH (r), t ME (r), t FH (r), t BC (r), and t SC (r), respectively.
The order of these events is trivial except for shell-crossing singularity. This is because, except for shell-crossing singularity, the events are determined solely by each shell and its dynamics is identical to the FLRW spacetime.
For W (r) ≥ 1, we can find t BB (r) < t PH (r) and there is no maximum expansion, no future apparent horizon and no big crunch. For the big bang time, Eq. (15) implies Since φ = S at apparent horizons, we find For 0 < W (r) < 1, we can find t BB (r) < t PH (r) < t ME (r) < t FH (r) < t BC (r). The big bang time and the time of past apparent horizon are given by Eqs. (54) and (55). The time of future apparent horizon is given by For the big crunch time, we find and the result is The time of maximum expansion is given by φ = S(r)/[1 − W 2 (r)], as seen from Eq. (6). Therefore, from Eqs. (15) and (16), we find The time of shell-crossing singularity is highly nontrivial because t SC (r) is determined not solely by the dynamics of the shell labeled r, but that of the neighboring shells. We will study it in detail in the next section.

Szekeres solution as nonlinear cosmological perturbations 4.1. Simultaneous big bang condition
The simultaneous big bang condition is often used in the cosmological context and suitable for a nonlinear growing mode of adiabatic perturbation. This condition implies H(r) = −t BB (r) = const. in Eq. (54). If we assume this condition, we find from Eqs. (15) and (16) in the limit of t → t BB , and hence the solution approaches the flat FLRW solution. This implies that the simultaneous big bang condition corresponds to extracting a purely growing 11/28 mode of adiabatic perturbation from the flat FLRW universe. This condition in the LTB solution has been adopted for the construction of the primordial black hole formation model in the flat FLRW universe in Harada et al. [24]. The LTB solution as an exact model of black holes in the evolving Universe is also discussed from a very broad scope in Sect. 18.9 of Ref. [46].

Shell-crossing singularities
We should note that very general treatments of the occurrence of shell-crossing singularities for the LTB solution and the Szekeres solution are given in Sects. 18.10 and 19.7.4 of Ref. [46], respectively. Our aim in the next few subsections is to analyze the occurrence of shell-crossing singularities and get a physical insight into the nature of shell-crossing occurrence in the present cosmological setting.
4.2.1. Expression for ψ = rφ ,1 /φ. In order to see the occurrence of shell-crossing singularities, it is important to have the explicit form of rφ ,1 /φ. In the general case, this can be found by differentiating Eq. (15) with respect to r as follows: where the derivative of G is given by the following form: With the simultaneous big bang condition H(r) = −t BB =const., Eq. (61) is reduced to Putting we can rewrite ψ in the following form: The function X(Y ) is plotted in Fig. 2. For W (r) = 1, we can see that the expression (65) is very singular. For this case, regularity is obvious in the following expression for ψ: In the limit W → 1 where 0 < S < ∞ and 0 < φ < ∞ are fixed, we find Y → 0 and d ln G/dY → 3/10. Therefore, for W (r) = 1, we find ψ(r,X) =ã 1 (r)X + a 0 (r), whereX and a 0 (r) is the same as in Eq. (64). From Eq. (41), the two-surface S t,r possesses a shell-crossing if and only if ψ ≤ β(r). According to Eq. (65) for W = 1 and Eq. (67) for W = 1, we can now analyze rigorously the occurrence of shell-crossing singularity for given t and r under the simultaneous big bang condition.
We should note that Eq. (43) implies a 1 = O(r 2 ) for W = 1 andã 1 = O(r 2 ) for W = 1. As for a 0 (r), Eq. (43) implies a 0 (r) → 1 as r → 0. Since β(r) = O(r) implies a 0 (r) > β(r) for sufficiently small r, there is no shell-crossing singularity for sufficiently small r. From Eqs. (65) and (67), since X for W = 1 andX for W = 1 both begin with 0 at the big bang, no shell-crossing singularity appears at least for a sufficiently short time interval after the big bang.
We now discuss in general the cases in terms of the shell labeled r with 0 < W (r) < 1, W (r) > 1, and W (r) = 1, separately. Hereafter, for simplicity, we choose the scaling of r so that S(r) = S 3 r 3 in accordance with Eq. (43), where S 3 is a positive constant. The result does not depend on the choice of the scaling. Recall that Eq. (46) holds for β(r).

4.2.2.
Bound shell: 0 < W (r) < 1. For 0 < W (r) < 1, as the shell labeled r begins with a big bang, reaches maximum expansion, and ends in a big crunch, Y monotonically increases 13/28 from 0, takes value 1 at maximum expansion, and then monotonically decreases to 0. Meanwhile, X = X + monotonically increases from 0 to 1, which is its value at maximum expansion, and then switches to the X = X − branch, which monotonically increases from 1 to ∞ as time proceeds. We should note that X − admits the following expansion: So, if a 1 (r) > 0, ψ monotonically increases from a 0 (r) to ∞ as time proceeds. Therefore, we can conclude that there appears no shell-crossing singularity. This is also the case if a 1 (r) = 0, where ψ is constant in time. If a 1 (r) < 0, ψ monotonically decreases from a 0 (r) to −∞.
Since β(r) ≥ 0, this means that a shell-crossing singularity necessarily appears before the big crunch. In this case, it is important whether or not the shell-crossing singularity appears before the maximum expansion and the future apparent horizon. A maximum expansion is characterized by Y = 1 or X = 1. Since ψ is a monotonically decreasing function of time, we can conclude that a shell-crossing singularity at r appears before or coinciding with a maximum expansion if and only if a 1 (r) + a 0 (r) ≤ β(r). A future apparent horizon is characterized by φ = S in the collapsing branch. The value of ψ on the future apparent horizon is given by ψ = a 1 (r)X − (1 − W 2 (r)) + a 0 (r). Since ψ is a monotonically decreasing function of time, we can conclude that a shell-crossing singularity at r appears before or coinciding with a future apparent horizon if and only if If |1 − W 2 |/S 2/3 is a monotonically decreasing function of r, a 1 (r) > 0 and hence there is no shell-crossing singularity. Otherwise, there exists r(> 0) for which a 1 (r) < 0 and a shellcrossing singularity appears. Then, Eq. (70) determines whether or not the shell-crossing singularity appears before or coinciding with a future apparent horizon. We can see that a large dipole moment can promote and advance the occurrence of shell-crossing singularity before the future apparent horizon.

4.2.3.
Unbound shell: W (r) > 1. If W (r) > 1, the shell labeled r begins with a big bang and expands forever. In this case, Y monotonically decreases from 0 to −∞ and X monotonically decreases from 0 to −1/2. So, if a 1 (r) < 0, ψ monotonically increases from a 0 (r) to −a 1 (r)/2 + a 0 (r) as time proceeds. Therefore, we can conclude that no shell-crossing singularity appears. This is also the case if a 1 (r) = 0, where ψ is constant in time. If a 1 (r) > 0, ψ monotonically decreases from a 0 (r) to −a 1 (r)/2 + a 0 (r). Since β(r) ≥ 0, this means that a shell-crossing singularity eventually appears in the course of expansion if and only if −a 1 (r)/2 + a 0 (r) ≤ β(r), i.e., This means that if W has an extremum which is greater than 1, there necessarily appears a shell-crossing singularity, irrespective of the value of β(r). A cosmological void in the asymptotically flat FLRW spacetime is characterized by W (r) > 1 near the center and W (r) → 1 in the asymptotic region r → ∞. Because of Eq. (43), W must have a maximum which is greater than 1 in this case. Therefore, Eq. (71) implies that a shell-crossing singularity necessarily appears. This also applies to PBH formation which has a region where W > 1. Therefore, for PBH formation without shell-crossing singularities before a future apparent horizon, W ≤ 1 is necessary for all r ≥ 0. 14/28 4.2.4. Marginally bound shell: W (r) = 1. For W (r) = 1, the shell labeled r begins with a big bang and expands forever. For this case, Eq. (67) is relevant.X monotonically increases from 0 to ∞ as time proceeds. Ifã 1 (r) > 0, ψ monotonically increases from a 0 (r) to ∞ and hence there is no shell-crossing if Eq. (46) is satisfied. This is also the case ifã 1 (r) = 0, where ψ is constant in time. Ifã 1 (r) < 0, ψ monotonically decreases from a 0 (r) to −∞. In this case, a shell-crossing singularity appears in the course of expansion.
Therefore, the condition for the occurrence of shell-crossing singularity is given by W (r) < 0. If the region inside the shell with W = 1 is unbound (W > 1), and the region outside is bound (0 < W < 1), a shell-crossing singularity appears on this shell.

Criterion in terms of the density distribution.
We can define the following quantities:ρ We can identifyρ with the density averaged over the two-surface S t,r , and ρ with the density averaged over the three-ball of which the surface is given by S t,r . The explicit expressions of the above quantities for the Szekeres solution are given bȳ Recall that the regularity at t = t 0 requires β(r) < 1. We can find that the ratio of the ball-averaged density to the shell-averaged density is calculated to give According to the above analysis, we conclude that if there exists t at which ( ρ /ρ)(t, r) ≥ 1, then no shell-crossing singularity appears at r. For 0 < W (r) < 1, even if ( ρ /ρ)(t, r) < 1, no shell-crossing singularity appears before a future apparent horizon if ( ρ /ρ)(t FH (r), r) > β(r). For W (r) > 1, ( ρ /ρ)(t, r) ≥ 1 for some t is a sufficient condition for the absence of shell-crossing singularity. Even if ( ρ /ρ)(t, r) < 1, no shell-crossing singularity appears if Eq. (71) is not satisfied. For W (r) = 1, ( ρ /ρ)(t, r) ≥ 1 is a necessary and sufficient condition 15/28 for the absence of shell-crossing singularity. Note that ifρ(t, r) is a monotonically decreasing function of r for fixed t, ( ρ /ρ)(t, r) ≥ 1 holds for all r > 0 but not vice versa. Interestingly, Szekeres [35] has proven that for zero-energy collapse and time-symmetric collapse, the condition for shell-crossing singularities to appear before the future apparent horizon is the same as we have obtained in the simultaneous big bang collapse, where the zero-energy collapse means W (r) = 1 for all r ≥ 0 and the time-symmetric collapse means that 0 < W (r) < 1 and t ME (r) =const. for all r ≥ 0. The simultaneous big bang condition is a physical condition, which is independent from the zero-energy condition and timesymmetric condition. If we assume the simultaneous big bang condition and zero-energy condition simultaneously, the solution reduces to the flat FLRW solution and hence there is no perturbation from it. Also, if we assume the simultaneous big bang condition and timesymmetric condition simultaneously, the solution reduces to the closed FLRW solution and again there is no perturbation from it.

Shell-focusing singularity
Let us move on to the possibility that a shell-focusing singularity is naked. As we have seen, noncentral shell-focusing singularities are covered by a future apparent horizon. Therefore, we focus on central shell-focusing singularity as a result of gravitational collapse and hence we assume 0 < W < 1 in the neighborhood of the center. In accordance with Eq. (43), we assume If W 2 = 0, the central shell-focusing singularity does not appear in finite proper time and hence we assume W 2 > 0. From Eqs. (18), (56), and (58), we find Then, if W 4 ≥ 0, the central shell becomes singular only after its neighborhood gets trapped and hence there is no null geodesic which emanates from the central singularity. So, we only have to study the case where W 4 < 0. Equations (1), (2), and (3) imply that, along the outgoing radial null tangent vector, the following equation must be satisfied: Using Eqs. (27) and (63), if there is an outgoing radial null tangent from the center, the following equation must be satisfied where we have used Eqs. (6) and (17). Let us put Y = Yr α and transform Eq. (85) to where α is constant. Since 1 − W 2 = O(r 2 ) for regularity, we can assume 0 ≤ α ≤ 2. If α > 2, the center is trapped and hence there is no outgoing null curve from the center. We are 16/28 interested in the limit to the center, while Y approaches a positive finite value. Let us write the right-hand side of Eq. (86) as F (r, Y). Because of Eq. (81), we can expand a 1 (r) in the form and ψ admits the expansion where we have used Eq. (69) and left the possible lowest-order terms only. Note that the smaller the value of α is, the earlier the outgoing radial null curve is. Since we are interested in the earliest one, we first assume 0 < α < 2. Then, the possible lowest-order terms of F are given by Let us put α = 4/3. In this case, we can write F (r, Y) as Noting that W 2 > 0 and W 4 < 0, F (r, Y) can be rewritten as where Thus, we can see that This means there exists a single outgoing radial null curve which emanates from the central shell-focusing singularity. For α = 4/3 but 0 < α < 2, we can find that F has no root for a positive and finite Y. The existence of the radial null tangent ensures the existence of an outgoing null geodesic which emanates from the singularity. In summary, there appears an outgoing null curve with radial null tangent at the center if and only if W 2 > 0 and W 4 < 0. The global visibility in the Szekeres spacetime in the sense of an event horizon is not analytically tractable in general because there is no geodesic which remains radial [38]. It should again be noted that the global visibility in the sense of an event horizon is not physically clearly motivated in the cosmological context as in the present case.

Condition for PBH formation models
In the context of cosmological perturbation in the flat FLRW universe, we cannot generally expect a monotonically decreasing density profile because the density perturbation should be given randomly according to some probability distribution function. For example, in order to have a model where the deviation from the flat FLRW solution falls off sufficiently fast in the asymptotic region, the mass excess in the overdense perturbation must be compensated by the surrounding underdense perturbation. This is called a compensated model [11]. Therefore, it is clear that the requirement of a monotonically decreasing density profile is too restricted to discuss the evolution of nonlinear cosmological perturbation in the Universe. Primordial black hole formation in the asymptotically flat FLRW spacetime is characterized by 0 < W < 1 near the center and W → 1 in the asymptotic region r → ∞. However, the shell labeled r with W (r) < 1 eventually collapses. Therefore, even if W → 1 as r → ∞, the mass of the black hole becomes infinite. However, in this situation, t FH (r) → ∞ as r → ∞ and hence the infinitely massive black hole is unphysical. To avoid this unphysical situation, we introduce a cut-off scale r s > 0 so that we assume W (r) = 1 for r > r s . Khlopov and Polnarev [12] argue that in the context of the GUT phase transition, the occurrence of caustics of massive particles will prevent collapse to a black hole because modeling by a dust fluid is no longer valid, where the matter's pressure gradient force cannot be neglected in such a high-density region. More precisely, they assume in the LTB solution that if a central shell-focusing singularity is not behind a future apparent horizon, then it can be regarded to prevent collapse to a black hole because the resultant high density implies the breakdown of the dust approximation of the real matter field and the collapse of the shell is impeded by a strong pressure gradient force. We extend their assumption to the Szekeres solution. Also, at shell-crossing singularities, the density grows infinitely large and hence the dust description should break down. On one hand, the pressure gradient cannot be neglected there and the free-fall collapse will be altered. On the other hand, it is not so clear whether or not it impedes black hole formation.
We can construct a model that is free from shell-crossing singularity if and only if S and W satisfy 0 < W (r) ≤ 1, lim r→∞ W (r) = 1, and −r ln We should note that Eq. (94) is equivalent to It is interesting that once we assume the first condition (93), no dependence on nonsphericity appears in the second condition even though the model can be highly nonspherical. However, shell-crossing singularities formed after the formation of a future apparent horizon do not affect the formation of primordial black holes. This means that in the present context we can relax the condition to that for the absence of shell-crossing singularities before the future apparent horizon formation, so that Eq. (94) is replaced by

Global visibility
In the previous section, we have studied whether or not shell-crossing and shell-focusing singularities at r are formed before the shell labeled r is trapped by a future apparent horizon. This corresponds to the question whether or not singularities are locally naked. We can also ask whether or not the singularity is globally naked. For a singularity, if a light ray which emanates from the singularity can get out of the black hole, we call this singularity globally naked. If there is no such light ray, then we call this singularity globally covered.
In asymptotically flat spacetimes, we can define globally naked singularities in terms of future null infinities. In the cosmological setting this is not so well motivated, as we have already seen. Here we deal with this problem by specifying the mass of the black hole.
If we fix the mass of the black hole to M , we can identify r s which gives the black hole mass, i.e., S(r s ) = 2M . To investigate the global visibility, we need to track null geodesics which emanate from the singularity. However, this problem is very difficult to tackle not only because null geodesic equations cannot be integrated analytically in general but also because null geodesics cannot be kept radial even though they started as radial null ones momentarily. For this reason, we study a sufficient condition for the singularity to be globally covered. This can be done by requiring that the singularity which occurs at t is surrounded by a trapped sphere at the same time t. For the central shell-focusing singularity, this condition yields For the shell-crossing singularity at r, this condition yields Equation (97) is equivalent to while Eq. (98) has no compact expression in terms of the arbitrary functions. The left-hand side of Eq. (98) decreases if β(r) increases, while the right-hand side does not depend on β(r).

Exact models of PBH formation
Here we construct exact models of PBH formation. We fix the scaling of r so that S(r) = r 3 in this subsection. Then, only the choice of W (r) determines the reference LTB solution.

Model A.
First, we consider the following natural choice of W : where r c is a positive constant and is constant. We call this choice of W (r) "model A". For consistency, we need ≤ 1. In this model, as we have seen, the central shell-focusing singularity is naked if > 0, while it is covered if ≤ 0. However, the global visibility is not so trivial. Let us choose r c = 1 and = 0, 0.1, and 0.2. Then, the time of future apparent horizon, t FH (r), is plotted in Fig. 3. We can see that the model with = 0.1 has a central 19/28 shell-focusing singularity which is locally naked but globally covered. For ≥ 0, since |1 − W 2 |/S 2/3 is a monotonically decreasing function of r, no shell-crossing singularity appears. It should be noted that the function W (r) is discontinuous at r = r c . With the above choice, the model is only perturbed within the comoving radius r s = r c and the region outside it is identical to the flat FLRW universe. This is also the case irrespective of the choice of the nonspherical functions because with the simultaneous big bang condition, the function φ becomes identical to that in the flat FLRW spacetime, and this reduces to the flat FLRW spacetime as we have already seen. So, in this model the compensation is perfectly realized at r s = r c .  Next, we present another example free from shell-crossing singularity and shell-focusing naked singularity as follows: where r c and r w are positive constants, and the following inequality must be satisfied for W 2 > 0: r c /r w > 2 6 /3 9/2 . We call this choice of function W (r) "model B". In this case, the metric function is C 2− and we have neither naked shell-focusing singularity nor shell-crossing singularity because W 4 = 0 and |1 − W 2 |/S 2/3 is a monotonically decreasing function of r. Also, this model is compensated at r s = r w , irrespective of the choice of the nonspherical functions.

Model C.
We can consider another example, for which where r c , r w (< r c ), and r f are positive constants, and the following inequality must be satisfied for W 2 > 0: We call this choice of W (r) "model C". This model is also free from shell-focusing naked singularity and shell-crossing singularity. For 0 < r < r w , from the simultaneous big bang condition, the function φ is identical to that in the closed FLRW spacetime. Therefore, the region inside r w is identical to the closed FLRW spacetime irrespective of the choice of the nonspherical functions. Since the fall-off of the function in the asymptotic region r → ∞ is sufficiently fast, this model is compensated at infinity [11], while the mass of the black hole becomes infinite because r s = ∞. We should note that the above choices of the energy function are similar to those presented in Ref. [24] for the PBH formation with the LTB solution but with a slightly different scaling of r.

5.3.4.
Model D. Now, we propose an example in which shell-crossing appears but only after the future apparent horizon forms. Here, we consider a polynomial function of W (r) parametrized as follows: where r c and r w are positive constants and n 1 and n 2 are even integers, and for consistency, the following inequality must be satisfied: r c /r w > √ f max , where f max is the maximum value of the function f (x) = x 2 (1 + x n1 − 2x n2 ) 4 for 0 < x < 1. If n 1 > 2 and n 2 > 2, shell-focusing singularity is censored. We call this choice of W (r) "model D". This model is compensated at r s = r w , irrespective of the choice of the nonspherical functions. We put n 1 = 8, n 2 = 10, r c = 0.8, and r w = 1 and plot ψ(r, 1) and ψ(r, 1 − W 2 ) as functions of r in In this case, since a 1 (r) is not positive definite, ( ρ /ρ)(t, r) ≥ 1 does not hold for some r. Noting that ψ ≤ β(r) implies shell-crossing singularities, we can understand the evolution of this model. Notice that ψ(r, 0) = 1, ψ(r, 1), and ψ(r, 1 − W 2 ) are the values of ψ at the big bang, maximum expansion and future apparent horizon, respectively. ψ is a monotonically increasing function of time if ψ(r, 1) > 1 and monotonically decreasing function of time if ψ(r, 1) < 1 because ψ(r, Y ) = a 1 (r)X(Y ) + 1. We should note that in our scaling, β(r) < 1 must be satisfied from Eq. (46). So, if β(r) = 0, which is the LTB model, no shell-crossing singularity appears before the future apparent horizon forms. This is also the case if β(r) is sufficiently small for 0 < r < r w . But if there is an r at which ψ(r, 1 − W 2 (r)) < β(r) < 1, shell-crossing singularity appears before the formation of a future apparent horizon. Moreover, if there is an r at which ψ(r, 1) < β(r) < 1, shell-crossing singularity appears at r even before the maximum expansion. This shows that a large dipole moment can promote shell-crossing singularity occurrence before the formation of an apparent horizon: if shellcrossing singularity is to occur, the larger the dipole moment is and the earlier the time of shell-crossing singularity becomes.

Conclusion
We have constructed exact models of spherical and nonspherical formation of primordial black holes with the LTB solution and Szekeres's quasispherical solution of the Einstein equation. The matter content is restricted to being a pressureless fluid or dust. The LTB solution contains three arbitrary functions of one variable with one scaling freedom remaining and hence describes a general spherically symmetric and inhomogeneous dust spacetime. The Szekeres solution additionally contains four arbitrary functions of one variable with one algebraic constraint equation. We interpret the Szekeres solution as the nonspherical deformation of the LTB solution by adding dipole moment distribution. These solutions may contain both shell-crossing and shell-focusing singularities. We use these solutions to model the evolution of cosmological nonlinear fluctuations. In this context, these singularities 22/28 are regarded as the breakdown of the model with a pressureless fluid, where strong pressure gradient force instead is at work.
Based on the Szekeres solution, we have analyzed how inhomogeneous effects and nonspherical effects affect the formation of naked singularities from cosmological fluctuations. If the perturbation has a sufficiently homogeneous central region and a sufficiently sharp transition in the matching with the background flat FLRW universe, then the central shellfocusing singularity is globally covered. Moreover, if the central density concentration is sufficiently large, no shell-crossing singularity appears from regular initial data, irrespective of the dipole moment distribution. However, if the central density concentration is not sufficiently large, shell-crossing singularity appears. In this case, large dipole moment distribution significantly advances shell-crossing singularities and they tend to appear before a future apparent horizon. This is the first exact approach to nonspherical effects on primordial black hole formation in full general relativity. The current analysis on the Szekeres solution is important in the cosmological epoch when the effective equation of state is very soft, which can be realized in the ending phase of inflation.
Hereafter we redefineφ as φ, and so on. Defining W (r) := 1/h(r), we can write k(r) in terms of W (r) as k(r) = K − W 2 (r). Thus, Eq. (A32) becomes The solutions are called quasispherical, quasiplanar, and quasi-pseudospherical, respectively, for K = 1, 0, and −1. Equation (A35) is integrated to give where H(r) is an arbitrary function of r. Thus, the solution contains seven arbitrary functions with an algebraic constraint. Substituting the explicit expressions for λ and ω into G00 = 8π , after a rather lengthy calculation we can obtain the following compact expression for : 8π = S P − 3SP ,1 φ 2 (φ ,1 P − φP ,1 ) .
Substituting the expressions for λ and ω, we can show that Eqs. (A13)-(A17) identically vanish and hence all the components of the Einstein equations are satisfied.