On the propagation across the big bounce in an open quantum FLRW cosmology

The propagation of a scalar field in an open FLRW bounce-type quantum spacetime is examined, which arises within the framework of the IKKT matrix theory. In the first part of the paper, we employ general-relativity tools to study null and timelike geodesics at the classical level. This analysis reveals that massless and massive non-interacting particles can travel across the big bounce. We then exploit quantum-field-theory techniques to evaluate the scalar field propagator. In the late-time regime, we find that it resembles the standard Feynman propagator of flat Minkowski space, whereas for early times it governs the propagation across the big bounce and gives rise to a well-defined correlation between two points on opposite sheets of the spacetime.

rics. The three-dimensional submanifold of the spacetime where the metric is degenerate represents a spacetime defect and it has been recently suggested that its origin can be explained within the IIB matrix model [5][6][7][8]. Universes having a bouncing-like behavior have been addressed in the context of f (R) theories in Refs. [9][10][11] and in Gauss-Bonnet modified gravity in Ref. [12]. Furthermore, the authors of Refs. [13][14][15] have constructed classical nonsingular bouncing models by means of generalized cubic Galileon theories, which permit to realize a pattern where the null energy condition is violated without introducing ghost and gradient instabilities. Another viable theory is represented by the matter bounce scenario [16]. Last, the bounce mechanism has been investigated in the framework of loop quantum cosmology [17][18][19][20] and string cosmology [21,22]. For a review on bouncing cosmologies we refer the reader to Refs. [23][24][25][26].
Recently, solutions of matrix models have been found [5][6][7], which can be interpreted as 3 + 1-dimensional quantum geometries 1 describing an effective FLRW cosmology with a big bounce (BB). The underlying model is known as IKKT model [33], which has been proposed as a constructive definition of (some corner of) string theory. These solutions or backgrounds have an intrinsic quantum structure, with spacetime uncertainty or "fuzzyness" akin to quantum mechanical phase space. In this framework, a classical spacetime geometry is recovered in the semi-classical or IR regime, while the quantum structure of geometry becomes important only in the UV regime, i.e. at very short distances. In particular, the singularity of classical geometry at the BB is completely under control. Fluctuations of such a background lead to fields, including scalar fields, gauge fields, and gravitons; in fact for the geometry given in [7], a whole tower of higherspin gauge fields arises, which is described by a ghost-free higher-spin gauge theory [34].
In particular, the framework of matrix models allows to study the physics on such backgrounds with a BB. This study was initiated in [35] on a 1 + 1-dimensional toy model, which allowed to compute the propagator for the global geometry including the BB. Furthermore, the Bogoljubov coefficients which govern the asymptotic properties of fields propagating in and out of the BB were obtained.
In the present paper, we extend the results of [35] to the case of 3+1 dimensions. We obtain explicitly the fluctuation modes of a scalar field on the 3 + 1-dimensional FLRW background in the semi-classical regime, paying special attention to the asymptotic regimes of late times and close to the BB. This then allows to compute the propagator explicitly, using a path integral approach provided by the underlying matrix model framework. Remarkably, the propagator is found to be regular at the BB, and has the local structure of a Feynman propagator at late times. This means that physical modes can 1 See also e.g. [27][28][29][30][31][32] for related work.
propagate across the BB singularity in a well-defined way, and some of their structure will survive.
The main message of this paper is that matrix models provide a suitable framework for quantum geometry, which allows to address and resolve the singularities which arise in the framework of general relativity. Moreover, the present example demonstrates how a time evolution and a 3 + 1dimensional causal structure can emerge from the underlying matrix model, which has no a priori notions of space and time. It should be emphasized that the dynamics is governed here by the matrix model, which is different from general relativity, at least at the classical level. As demonstrated in [36], the Einstein-Hilbert action does arise upon including 1-loop effects in the IKKT model, under suitable assumptions. This will of course affect some of the results of the present paper, however we expect that the qualitative features of the present classical analysis will also apply after including quantum effects.
The paper is morally divided into two parts, a classical and a quantum one. The two parts are consistent with each other. In the classical part, we study some aspects of the present FLRW geometry, with special emphasis on the near-BB regime. In particular, we elaborate the geodesics, and show that they extend smoothly across the BB. The BB singularity is found to be rather "mild" in a sense that will be explained below. The propagator is obtained in the quantum setting, by computing explicitly the (free) path integral of modes as defined by the matrix model. At late times, we recover again the standard Feynman propagator with the appropriate iε structure. At early times near the BB, the propagator also turns out to be well-defined, and allows to study the propagation of scalar particles across the BB. This result agrees with the classical analysis regarding null and timelike geodesic.

The background geometry
We recall [7] that the background M 3,1 under consideration can be described semi-classically as a projection of fuzzy H 4 n , which is obtained from five matrices X a ∼ x a interpreted as quantizations of five embedding functions where a = 0, . . . , 4. A convenient parametrization of this four-dimensional hyperboloid is as follows for η ∈ R. Note that χ can be restricted to be positive. Projecting this along the x 4 axis leads to a 2-sheeted cover of the following region where the upper sheet ("post-BB", corresponding to x 4 > 0) is covered by η > 0, while the lower sheet ("pre-BB", corresponding to x 4 < 0) is covered by η < 0. The BB separates these sheets, and corresponds to x μ x μ = −R 2 . This leads to the following parametrization of M 3,1 Note that the flow of time will be along increasing η on both sheets; this arises from the iε regularization discussed in Sect.

5.
In principle, we can restrict χ to be either positive or negative on either sheet. However, we will see in Sects. 4 and 5 that it is convenient to choose χ > 0.

Effective metric
To understand the effective metric on M 3,1 , we recall [7] that the background solution T μ of the matrix model leads to the following kinetic term which governs all fluctuations in the matrix model. Using the semi-classical relation [T μ , .] ∼ i{t μ , .} in terms of Poisson brackets (here t μ represents the semiclassical limit of T μ ) and recalling {t μ , x ν } = sinh(η)δ μ ν , this can be rewritten uniquely in the standard form [7,37] where [7] G μν = | sinh(η)| −3 γ μν , γ αβ = sinh 2 (η)η αβ , dropping some irrelevant constant (here γ μν is an auxiliary metric which is relevant for the torsion). This metric can be recognized as SO(3, 1)-invariant FLRW metric, where is the invariant length element on the space-like hyperboloids . Equivalently, we can write with r = sinh χ.
From Eq. (8), we can read off the cosmic scale parameter a(η) and the relation linking the differentials dt and dη, i.e., For late times, we have This is a reasonable FLRW cosmology given the simplicity of the model, which is asymptotically coasting at late time with a(t) ∼ 3 2 t, cf. Refs. [38,39]. Note that it arises directly from the matrix model, without using or assuming general relativity. Last, it is worth mentioning that the scale factor (12) can be either odd or even (the odd solution, in particular, may be of interest as discussed in Ref. [40]).

Classical analysis of the FLRW spacetime
In this section, we perform a classical investigation of the FLRW geometry (8). Curvature invariants are considered in Sect. 3.1, whereas null and timelike geodesics are studied in Sect. 3.2. We conclude the section with the analysis of some cosmological observables (see Sect. 3.3).

Curvature invariants
In order to describe the behaviour of the spacetime near the BB, we begin our analysis with the investigation of some curvature invariants. Starting from the metric (8), we find that the Kretschmann scalar is the squared Ricci tensor reads as whereas the (topological) Euler invariant is [41][42][43] the star indicating the duality operation. It is clear that the scalars (15)- (17) blow up at the BB, i.e., at η = 0.

Null and timelike geodesics
In this section, we investigate null and timelike geodesics of the FLRW geometry having θ and ϕ constant.
Null geodesics can be described in terms of the function χ(η). We seek a solution which reaches the BB at η = 0 after having travelled toward it for η < 0 and away from it when η > 0. Therefore, it follows from Eq. (8) that the motion in the outward χ -direction is parametrized by the differential equation which, with the boundary condition χ(η = 0) = 0, leads to The behaviour of the solution (19) is shown in Fig. 1, whereas Fig. 2 represents the plot obtained by means of the embedding functions (2). Having obtained a continuous geodesic solution which can be extended uniquely at η = 0, we can conclude that light (and hence the physical information) can travel across the BB, despite the singularity occurring in the invariants (15)- (17).
It is interesting to work out the solution χ = χ(t) of null geodesics (with θ and ϕ constant) for early times, i.e., when t → t 0 (with t 0 ∈ R) or, equivalently η → 0. First of all,  The inversion of the above function yields where we note that η(t = t 0 ) = 0.
We are now ready to obtain the expression of the cosmic scale factor valid near the BB. Indeed, by means of Eqs. (12) and (21), we have which vanishes at t = t 0 . We note that a(t) is positive (resp. negative) for t > t 0 (resp. t < t 0 ). The sign of a(t) drops out in the metric (8), but the above choice has no cusp. The equation of null geodesics (with θ and ϕ constant) can be written in the equivalent form which in view of (22) gives the desired early-time solution see Fig. 3. The above solution can also be obtained if we first expand Eq. (19) about η = 0, and then exploit Eq. (21).
At this stage, we analyze the timelike geodesics of a noncomoving observer. Hence, let denote the unit timelike four-velocity vector of such observer, whose proper time is indicated with τ . If we suppose, like before, that the observer moves along the direction of constant θ and ϕ, then from the normalization condition we obtain Furthermore, if we consider the χ -translational Killing vector ξ α , the conserved momentum of the non-comoving observer is given by The above equation gives a relation to express dχ(τ )/dτ in terms of and a(η), which can be exploited in Eq. (27). In this way, after some manipulations, we end up with the geodesic equation in terms of the function χ(η), i.e., We have solved numerically Eq. (29) with the boundary condition χ(η = 0) = 0. As it is clear from Fig. 4, the non-comoving observer can travel across the BB, likewise a massless particle.
We can provide the analytic solution for the timelike geodesic of the non-comoving observer moving in the χdirection in the regime of small times. First of all, it is easy to see that the equation of the timelike geodesic of the noncomoving observer in terms of the cosmic time variable is Bearing in mind Eq. (22) and considering, for simplicity, t 0 = 0, R = 1, and = 1, the solution of Eq. (30) with the The plot of the function (31) is given in Fig. 5.
Although we have shown that both null and timelike geodesics have no pathological behaviour, it should be noted that the corresponding velocities blow up at the BB, as it is clear from Eqs. (23) and (30).

Singularity or not? "Mild singularity" or a new kind of singularity
In the general relativity framework, singularity theorems are based on the criterion that timelike and null geodesic completeness are minimum conditions for a spacetime to be considered singularity-free [44][45][46]. However, these theorems do not prove that singularities of spacetime are necessarily related to unboundedly large curvature. Indeed, the characterization of singularities via the divergent behaviour of the curvature can be inadequate in some situations (e.g., the case of conical singularity). Our model is not framed in general relativity, but emerges in the semiclassical limit of the matrix theory. In other words, the cosmic scale factor occurring in the metric (8) is not a solution of Einstein field equations. In our analysis, a(η) vanishes at the BB, where, in addition, the invariants (15)-(17) blow up. Furthermore, the effective metric (8) is zero (and hence is degenerate) at η = 0. Despite that, no pathological behaviour occur in the analysis of null and timelike geodesics, as we have shown before. Therefore, if we are to make a comparison with the recipes of general relativity, we could say that our FLRW solution features a "mild singularity" 2 which does not prevent both massless particles and non-comoving observers from crossing the BB, but introduces some diverging curvature scalars. Moreover, owing to the degenerate character of the metric (8), another option is possible, where we could conclude that we are dealing with a new kind of singularity, which, in principle, might be present also in general relativity, as discussed by Hawking in Ref. [49].

Cosmological observables
In this section, we evaluate some cosmological observables of our FLRW geometry. In Sect. 3.3.1 we consider the particle horizon, while in Sect. 3.3.2 we deal with the luminosity distance.

The particle horizon
In terms of the time coordinate, the particle horizon D hor (t) at t > 0 reads as [50] where we note that, for the background under consideration, we are free to choose −t arbitrarily negative. In general relativity settings, this operation would be correct only if the spacetime were singularity-free. However, following our discussion of Sect. 3.2.1, it makes sense to consider such a limit in our model, since we are dealing with a new kind of singularity which cannot be totally explained with the standard tools of Einstein theory.
Bearing in mind Eqs. (12) and (13), the particle horizon where η ≡ η(t ). Since D hor (η) diverges, it is possible to receive signals from any coming particle of the spacetime. This can be interpreted as a hint that our model has no need for inflation.

Luminosity distance
Let us consider a light signal emitted by a comoving source at time t e > 0 which is then detected by a comoving observer at time t 0 > t e . If we further suppose that the signal moves keeping θ and ϕ constant, the luminosity distance d L (t e ; t 0 ) is given by [50] where the coordinate distance r e travelled by the signal is expressed by the relation In our model, it follows from Eqs. (11)-(13) that the luminosity distance can be expressed as where η e ≡ η(t e ) > 0 and η 0 ≡ η(t 0 ) > η e . It is easy to see that which represents an expected result since the cosmic scale factor a(η) vanishes at the BB.

The scalar modes on M 3,1
In this second part of the paper, we would like to compute the propagator for a scalar field on the above background, and see if there are any interesting effects due to the BB. A priori it is not evident how to define the propagator, due to the boundary provided by the BB. We take as starting point the definition as 2-point function defined by a Gaussian integral in the matrix model (or rather its semi-classical limit). Schematically, where Z being the generating functional. The necessary details for this computation will be provided below.

Relevant Operators
In this section, we introduce some relevant wave operators which will play a crucial role in our forthcoming analysis.

The effective d'Alembertian G
The metric (8) is encoded in the "matrix" d'Alembertian which governs the propagation of a scalar fields φ, and is related to the metric d'Alembertian through This can be seen by rewriting the action (6) as follows is the SO(4, 1)-invariant volume form on H 4 in Cartesian and hyperbolic coordinates (4), respectively (hereafter, we suppose that χ > 0). Explicitly, one finds (cf. (2.32) in [51]) where (3) η is the Laplacian on the spacelike three-dimensional hyperboloid H 3 where we have defined We note that to derive Eq. (45) we have exploited that the metric on the space-like H 3 is

The Laplacian operator G on H 4
The induced metric on the four-dimensional hyperboloid H 4 is (cf. Eq. (2)) where the length element d 2 on a spatial standard threedimensional hyperboloid H 3 has been given in Eq. (9). By means of the metric (48), the Laplacian operator G on a generic function φ = φ(η, χ, θ, ϕ) reads as where

Relations between G and
η and G and The Laplacian G on the four-dimensional hyperboloid H 4 and the Laplacian η on the spacelike three-dimensional hyperboloid H 3 are related by (cf. Eqs. (45) and (49)) On the other hand, the effective d'Alembertian G is related to (3) η as (cf. Eq. (44)) Last, we observe that the effective d'Alembertian G is related to the Laplacian G on H 4 as follows (cf. Eqs. (44) and (49)) An inspection of Eqs. (50)-(52) reveals that the eigenfunctions of the operators G , (3) , and share the same building blocks. The eigenfunctions of d'Alembertian operator are derived in the next section, while those of G and (3) can be found in Appendix A.

Eigenfunctions of the d'Alembertian operator
It is convenient to work out the eigenfunctions of the operator instead of G for at least two reasons. First of all, the d'Alembertian is self-adjoint with respect to the symplectic volume form (43); secondly, the operator admits handier eigenfunctions than G .
Bearing in mind Eqs. (41a), (44), and (45), the d'Alembertian operator reads as The eigenfunctions of the d'Alembertian operator are defined by the equation whose resolution can be tackled via the separation ansatz where the spherical harmonic functions Y m l (θ, ϕ) of degree l and order m (with l |m|) satisfy the well-known property [52] Bearing in mind Eqs. (55) and (56), the eigenvalue problem (54) gives where we have divided both sides byφ(η, χ)Y m l (θ, ϕ). After having performed some manipulations, the above equation can be solved through the method of separation of variables, yielding the following two ordinary differential equations: β being a real-valued constant. The solution of Eq. (58) and the eigenfunctions of the d'Alembertian operator (53) will be provided in the next sections.

The time-like equation
It is instructive to work out the details of the solution of Eq. (58a). If we introduce the variable then the derivative operators read as and hence Eq. (58a) becomes If we write then we end up with the general Legendre equation [53,54] ( which can be solved via the associated Legendre functions of the first and second kind P μ ν (w) and Q μ ν (w), respectively, having degree ν and order μ given by Therefore, bearing in mind Eqs. (59) and (62), the solution of Eq. (58a) in terms of the variable η is c 1 and c 2 being integration constants. Here Q μ ν can equivalently be replaced by P −μ ν , which will be done in the following.

The radial equation
The solution of (58b) can be written in terms of the Gaussian or ordinary hypergeometric function [53] as where c 3 and c 4 are integration constants and It is possible to write Eqs. (58b) and (66) in a more convenient form. Indeed, if we introduce the variable (recall that we are considering χ > 0) then we have and hence Eq. (58b) becomes The solution of the above equation reads as where we have exploited Eq. (68) and The physically meaningful solutions can be identified by considering the boundary conditions. Recall that χ ∈ [0, ∞) plays the role of a radial variable. In the limit χ → 0, we have (B.8) since coth χ ∼ χ −1 , while Pμ l (coth χ) is divergent for l 0, and finite for l = 0. Therefore to have non-singular solutions, we should choose the Qμ l solutions and reject the Pμ l .

The eigenmodes on the spacetime
Thanks to the solutions (65) and (71), the eigenfunctions (55) of the d'Alembertian operator (53) are known. Moreover, Eq. (64b) permits obtaining the explicit expression of the eigenvalue λ, which reads as In order to have oscillatory (square-integrable) solutions, we should assume that the order (64b) of the solution (65) is purely imaginary. Therefore, we set where Similarly, the order (72) of the solution (71) should be purely imaginary, i.e.
for q ∈ R. This leads to which implies that From Eqs. (76) and (78), the eigenvalue (74) can be written equivalently as As a consequence of Eq. (79), we see that the degree (64a) of the solution (65) is complex and reads as In order not to overcount equivalent solutions, we observe that the P μ ν (tanh η) coincide for the two choices of ν = − 1 2 ± iq. We can therefore make the convention that recalling that χ > 0. The above-defined eigenmodes are regular functions in the limit χ → 0. Indeed, due to the relation (B.8), we have Now consider the limit χ → +∞, which means that coth χ → 1 + . The functions Q iq l (x) are oscillatory but bounded when x → 1 + as it can be inferred from Eq. (B.9). In this way, we obtain , as χ → +∞.

Flat regime
Now consider the following "flat" regime 3 where "FR" stands for "flat regime", and q will be a typical momentum. Then the Q iq l reduces to the spherical Bessel functions, This should be expected, since the eigenfunctions (83) should reduce to the standard ones on R 3,1 . For χ → 0, the relation (87) is guaranteed by the standard expansion formulas for Q iq l (coth χ) and j l (qχ); 4 here, we have verified that it works very well in the range χ ∈ (0, 1) for q l. Moreover, it also holds that and hence we can conclude that Due to the relation [55] we easily obtain Thus, bearing in mind Eqs. (89) and (91), the eigenmodes (83) for large times (i.e., η → +∞) and in the flat regime (86) become 3 We will consider this regime also in the computation of the propagator for η → 0, where the geometry is strictly speaking no longer flat. 4 The expansion of Q iq l (coth χ) when χ → 0 can be read off from Eq. (B.8); for j l (qχ) we refer the reader to Ref. [52].
the symplectic volume form being given in Eq. (43). Explicitly, we obtain where we have exploited the orthogonality of the spherical harmonics [52] and the ensuing Kronecker delta factor δ ll . We first consider the integral involving the χ variable. If we define ξ = tanh χ , then we obtain which, via the substitution x = 1/ξ , yields The above integral can be worked out with the help of the results of Appendix B, where we have shown that The above relation leads to which in turn means that, for any real q, q , where we have exploited (see Eqs. (B.24)-(B.26)) By virtue of the identity (99), Eq. (94) gives We recall that the expression of the associated Legendre function of the first kind with a generic degree ν and argument x lying in the interval (−1, 1) is [53,55] P is 2 F 1 (a, b; c; x) being, like before, the Gaussian (or ordinary) hypergeometric function and (x) the gamma function. In our case, the degree ν of P is ν (tanh η) assumes the form given in Eq. (82). For this choice of ν, we find where we have employed the substitution y = tanh η. This last integral can be calculated with the help of following result: which has been proved in Ref. [56] (see Appendix B for further details). Therefore, bearing in mind the above formula along with Eq. (82), from Eq. (104) we finally obtain the sought-after orthogonality relations where a(q, s) b(q, s) = 2 sinh(π s) s 1 + cosh 2 (πq) sinh 2 (π s) = b(q, −s). (107b) To prove this, it suffices to multiply this equation with Q iq l (coth χ) and integrate over +∞ 0 dχ using (99). This completeness relation allows to represent any (squareintegrable) function φ(χ) defined in the interval (0, +∞) as superposition of Q iq l modes, whereφ(q) is given bŷ The factor q sinh (πq) e −2πq (π/2) 2 can of course be absorbed via a suitable redefinition ofφ(q).

Quantization and matrix version
The above discussion regarding the eigenmodes of the d'Alembertian (53) is completely classical, even though we are claiming to work in the framework of matrix models. This may seem suspect, but it can be fully justified as follows. Since the underlying fuzzy hyperboloid H 4 n is a quantized coadjoint orbit, there is an SO(4, 1)-equivariant quantization map 5 which establishes an isometric equivalence between commutative and the fuzzy "functions", and maps the matrix d'Alembertian to the above semi-classical version. Therefore all the eigenmodes and eigenvalues computed in the classical case carry over without corrections to the fuzzy case, in the free theory. This justifies the classical computations in this work. 5 Strictly speaking, the underlying space is twistor space CP 1,2 , and we are restricting ourselves to the lowest sector of scalar modes here [7].

Path integral quantization
In this section we compute the propagator of a scalar field φ(x) in the FLRW geometry (8). After having investigated the action in Sect. 5.1, the propagator of φ(x) is explicitly worked out in Sect. 5.2.

The action
In the semi-classical limit, the action for a scalar field φ(x) having mass m can be written as where the expression of can read from Eq. (43) and, as usual, ε is a small positive number which should be let tend to zero after integration. The exact matrix version of the action has the same form with the trace Tr replacing the integral d , and the iε term ensures that the matrix path integral Dφ e i S is well-defined [57]. We can evaluate S ε [φ] by employing the following decomposition of φ(x) in the basis of eigenmodes (83): φ + s,q,l,m , φ − s,q,l,m being the coefficients of such a decomposition. In order to ease the notation, in our forthcoming calculations we will set Bearing in mind Eqs. (80) and (113), the action (112) becomes Thanks to the orthogonality relations (106), the above equation gives where in the last passage we have exploited the fact that all terms proportional to δ(s + s ) give a vanishing contribution, since both s and s are positive (cf. Eq. (76)). If we introduce the square matrix then we can write and hence, in conclusion, the action assumes the form

The propagator in momentum space and position space
In the following, we evaluate the propagator of the massive scalar field φ(x). We first provide the general expressions for the propagator both in momentum space and position space. The late-time propagator is then elaborated in Sect. 5.2.1, recalling also a related topic -i.e. the expansion of plane waves in terms of spherical harmonics -in Sect. 5.2.2. Finally, in Sect. 5.2.3 we compute the propagator across the BB.
Starting from the action (119) and adopting the compact notation (cf. Eq. (114)) the propagator in momentum space reads as and (see Eq. (107)) Therefore, bearing in mind Eq. (121), the local propagator in position space reads as

The propagator in the flat regime and with η → +∞
We are mainly interested in the local propagator for distances far below the curvature scale, but keeping the oscillating nature of the modes in the late-time regime η → +∞. This is the flat regime FR defined in (86), where the eigenmodes (83) reduce to (92). It then follows from Eq. (124) that the late-time local propagator can be written as the sum of a leading piece and a subleading part, i.e., where "L" and "SL" stand for "leading" and "subleading", respectively. The leading term reads as 6 The above equation can be simplified by means of the identity which can be proved via the recurrence relation (1 + z) = z (z) [53] along with Eq. (B.2a). Moreover, the formula (127) leads to (cf. (86)) q sinh(πq) | (iq + l + 1)| 2 q l ∼ π q 2l+2 .
The restriction to q l means that we ignore the extreme IR regime of the propagator, which is justified for the typical applications of (quantum) field theory.
Therefore, by exploiting (128) and after some calculation, Eq. (126) becomes 6 Notice that if we would keep here only the asymptotic behavior of the spherical Bessel functions j l rather than their full form, the propagator would be ill-defined, as the sum over l, m would lead to a singular dependence on θ, ϕ.
where we have taken into account that the integrand is an even function of q due to the property j l (−x) = (−1) l j l (x).
Up to the η-dependent normalization factor which reflects the cosmic expansion, we recover precisely the Feynman propagator on a flat 3 + 1-dimensional Minkowski space, including the appropriate iε prescription which ensures local causality and determines the arrow of time along growing η.
The subleading contribution occurring in Eq. (125) is By exploiting Eq. (128) and after some arrangement, we find The last expression depends on the sum η+η and hence it is rapidly oscillating at late times (i.e., for η, η → +∞). This justifies the fact that Eq. (131) gives a subleading contribution with respect to Eq. (129), which, on the contrary, depends on the difference η − η . However, a crucial consistency check consists in proving that the leading term (129) has, apart from the normalization factor, the form of the usual local Feynman propagator on a flat four-dimensional spacetime. This task is fulfilled in the Sect. 5.2.2.

Plane-wave expansion in terms of the spherical harmonics
The scalar plane waves can be expanded in R 3 in terms of spherical harmonics via the Rayleigh equation [58] e i q· x = 4π where q = | q|, r = | x|,x = x r ,q = q q , and j l (qr) are the spherical Bessel functions [53]. In the above equation, the complex conjugation can be interchanged between the two spherical harmonics due to the symmetry of the scalar product q · x. Therefore, we have recalling that Y m l (−x ) = (−1) l Y m l (x ). Now we integrate this relation over all q with fixed radius q: using the orthogonality of the spherical harmonics. Therefore, we can write e.g. the propagator in the form The analogy with Eq. (129) should be noted.

The propagator in the flat regime and with η → 0
The behaviour of the scalar field (cf. Eq. (112)) near the BB can be inferred from the form assumed by its propagator (124) for small times, i.e., when η → 0. Here, we will also employ the flat regime (86). The definition of the associated Legendre function of the first kind P jointly with Eq. (82), leads to the expression P ±is ν (tanh η) where we have exploited the relation 2 The expansion about η = 0 of the hypergeometric function occurring above can be written as Thus, from the last relation jointly with Eq. (89), the earlytime expression of the eigenfunctions (83) evaluated in the flat regime (86) become As a consequence of Eqs. (124) and (141), we find that the early-time propagator can be written as where "ET" means "early-time" propagator. The first term reads as After some calculations, the leading contribution becomes where we have exploited Eq. (128) and the fact that the integrand is an even function of q. If we exploit the following relations [53] (2z) = 2 2z−1 √ π (z) (z + 1/2), then we can simplify the last term occurring in Eq. (144) as and hence, after some algebra, we get We can obtain an estimate for the last term appearing in the above equation. Indeed, taking into account the relation [59] | (x + iy)| 2 2 (x) we obtain the inequality and, upon taking into account Eq. (128) jointly with the fact that the integrand is an even function of q, we obtain, after some algebra, Therefore, if we sum the contributions (147) and (151), the early-time propagator (142) becomes The last term occurring in the above equation can be arranged as follows and hence, for η, η → 0, we get the final formula The early-time propagator (154) turns out to be bounded and well-defined owing to the upper bound (149). It gives rise to a more-or-less standard correlation between points x, x lying on the same sheet of the projected spacetime M 3,1 . More remarkably, it also gives rise to a well-defined correlation and in fact a propagation between two points located on opposite sheets of M 3,1 near the BB; for larger |η| and |η | this is suppressed compared to the case of two points on the same sheet, due to the oscillating term e is(η−η ) .
The observation that a scalar particle can cross the BB agrees at least qualitatively with the classical analysis regard-ing null and timelike geodesics of Sect. 3.2. A more detailed analysis of this fascinating result is left for future work.

Conclusions
In this paper, we elaborated in detail the scalar modes and the propagator on a quantum version of a 3 + 1-dimensional FLRW spacetime with BB, in the semi-classical limit. The underlying framework of matrix models provides a clean setup to work with a spacetime which is singular as a classical manifold, but well-defined as a quantum geometry.
The most interesting conclusion is that the physics of scalar fluctuations is perfectly well-defined even at or near the classical singularity, and it is possible to relate the preand post-BB eras in a meaningful way. This implies some intriguing correlation between the two sides of the BB, which remain to be worked out in detail. The framework of matrix models provides a clear prescription how to define and compute the propagator, and the causal structure of a Feynman propagator is obtained automatically, at least at late times. This is remarkable, since time emerges on the same footing as space from the underlying matrix model, and there is no a priori notion of classical time.
This work should be seen in the context of emergent spacetime and gravity within the IKKT matrix model, which is closely related to string theory [33]. The present background under considerations then arises as a classical solution [7], and a gravitational action arises from one-loop effects [36], cf. Ref. [37]. Fluctuations are governed by a noncommutative field theory, which for scalar fields reduces to the one elaborated here. This suggests that analogous properties can be expected for all fluctuations in the model. In particular, the resulting theory is expected to satisfy the standard causality and unitarity requirements in quantum field theory, given the emergence of a Feynman propagator and the maximal supersymmetry of the underlying IKKT model. Even though many aspects of such a theory remain to be understood, the results of the present paper should provide sufficient motivation for further studies.
The present analysis is restricted to non-interacting test particles on the background geometry. This is of course not entirely satisfactory, since the density of matter is expected to be singular at the BB, which would lead to modifications of the background and hence of the metric. The inclusion of such effects as well as the induced Einstein-Hilbert action would be very desirable, but is beyond the scope of this paper.
Finally, it is interesting to point out that recent numerical simulations of the bosonic IKKT model with mass term produced evidence for an emergent spacetime, with features reminiscent of a bouncing 1 + 1-dimensional cosmology [60]. One may therefore hope to eventually relate the matrix background under consideration to numerical simulations of the model.