Fermion production at the boundary of an expanding universe: a cold-atom gravitational analogue

We study the phenomenon of cosmological particle production of Dirac fermions in a Friedman-Robertson-Walker spacetime, focusing on a (1+1)-dimensional case in which the evolution of the scale factor is set by the equations of Jackiw-Teitelboim gravity. As a first step towards a quantum simulation of this phenomenon, we consider two possible lattice regularizations, which allow us to explore the interplay of particle production and topological phenomena in spacetimes with a boundary. In particular, for a Wilson-type discretization of the Dirac field, the asymptotic Minkowski vacua connected by the intermediate expansion corresponds to symmetry-protected topological groundstates, and have a boundary manifestation in the form of zero-modes exponentially localized to the spatial boundaries. We show that particle production can also populate these zero modes, which contrasts with the situation with a na\"ive-fermion discretization, in which conformal zero-mass fields exhibit no particle production. We present a scheme for the quantum simulation of this gravitational analogue by means of ultra-cold atoms in Raman optical lattices, which requires real-time control of the Raman-beam detuning according to the scale factor of the simulated spacetime, as well as band-mapping measurements.

1 Introduction Quantum field theory (QFT) provides a unifying language to describe quantum many-body systems at widely different scales. For instance, observed phenomena in high-energy physics can be accounted for by the standard model of particle physics [1], a QFT of fermions coupled to scalar and vector bosons. Here, Poincaré invariance determines the arena for such fields: the flat Minkowski spacetime of special relativity. At much smaller energy scales, within the realm of condensed-matter systems, non-relativistic QFTs are routinely used to explain various collective phenomena [2]. Interestingly, relativistic QFTs analogous to those of particle physics also arise in coarse-grained descriptions of certain phase transitions [3], or in materials such as graphene [4], Weyl semimetals [5], and topological insulators and superconductors [6]. Here, Poincaré invariance and an effective speed of light emerge at long wavelengths [7], such that the lowenergy excitations can be described by quantum fields in an effective Minkowski spacetime. As first realized in the context of the propagation of quantized sound waves in fluids [8], there are also situations in which the emergent invariance is related to general coordinate transformations, i.e. diffeomorphisms [9]. One then obtains emergent QFTs in a curved spacetime, leading to condensed-matter analogues of phenomena studied within the realm of general relativity [10].
Note that, in this so-called analogue gravity [11], the emerging spacetime metric g µν (x) is typically a classical field corresponding to a particular solution of Einstein's field equations [9]. Accordingly, these condensed-matter analogue systems do not aim at mimicking a full quantum theory of gravity in the laboratory, but rather at reproducing characteristic phenomena of QFTs under classical background gravitational fields [12], exploring the interplay of gravitation and quantum physics well below the Planck scale. The study of QFTs in curved spacetimes has lead to important predictions of this interplay, with paradigmatic examples being (i) the evaporation of black holes due to quantum effects [13], which elucidates on the thermodynamic nature of such objects; and (ii) particle production during inflation, which is crucial to understand the large-scale behaviour of the universe [14]. One of the attractive features of analogue gravity is that one can mimic these phenomena, which are notoriously difficult to observe in a real gravitational context, in a controlled tabletop experiment. In fact, the range of most applications of QFTs in curved spacetimes is believed to lie far away from any experimental probe [9]. One such example is the aforementioned evaporation of black holes [13], where the emission of thermal Hawking radiation from stellar-size black holes leads to vanishingly-small temperatures in comparison to the cosmic microwave background (i.e. 10nK versus T CMB ∼ 2.7K, and its observed inhomogeneities δT CMB ∼ 10µK). In the context of analogue systems, on the contrary, Hawking radiation of bosonic fields has already been observed in the propagation of either light in non-linear media [15][16][17][18], or sound in Bose-Einstein condensates [19][20][21][22]. The related phenomenon of the Unruh effect [23] has also been observed with Bose-Einstein condensates [24]. These experiments are leading to a paradigmatic shift: while, for many decades, gedanken experiments have been crucial to understand the interplay of general relativity and quantum mechanics, it is nowadays pos-sible to turn them into real experiments in analoguegravity labs.
In this manuscript, we will focus on cosmological particle production in an expanding universe [25,26], which can find analogues in the quantized sound waves of Bose-Einstein condensates [27,28] and in trapped-ion crystals [29]. In fact, the essence of quantum fields in an expanding spacetime has been recently observed in these two experimental platforms [30,31]. Further progress along these lines has allowed to implement various specific metrics of expanding curved spacetimes in the lab [32][33][34], opening the route to very promising future advances. Although most of the recent progress has focused on bosonic fields, these atomic experiments can also be performed with Fermi gases [35], and it would be very interesting to observe gravity analogues of the more elusive fermion production in expanding spacetimes [36]. In fact, one could go a step further and realise emerging spacetimes with exotic geometries and topologies that, despite being allowed by the theory, do not have any clear observational pathway in a real gravitational context. In fact, some of the above analogues have been realised in ring-shaped condensates [30], such that the bosonic fields are defined on spacetimes with a non-trivial spatial topology R × R → R × S 1 , where S 1 is a circle. In spacetimes with D = (d + 1) dimensions, one may consider R d → R d−n × S n , where there are n compactified spatial dimensions. Such toroidal topologies, originally addressed in the context of Kaluza-Klein compactification of extra dimensions [37,38], can lead to interesting consequences for the quantum fields, such as topological mass generation and topological symmetry restoration [39]. Other boundary conditions can also play a role in general relativity, such as in the context of black hole thermodynamics [40]. In QFTs, one may consider base spaces R d → R d−1 ×I, where I is a finite interval of the real line in which one imposes Dirichlet boundary conditions on the fields, leading to analogues of the Casimir effect in both its static [41] and dynamical [42] incarnations. We note that gauge field theories in manifolds with a boundary have also been explored [43][44][45], and can lead to an interesting bulk-boundary correspondence.
Coming back to the notion of extra dimensions, rather than compactifying them, one may instead interpret them as the bulk of certain lattice models displaying non-trivial topology in reciprocal space [46]. This bulk topology has a boundary manifestation in the form of field solutions that are exponentially localized within the boundaries. Remarkably, the lattice field theories describing these boundary degrees of freedom can display properties that cannot exist in the absence of a bulk, such as certain quantum anomalies [47][48][49]. This connects directly to the aforementioned condensed-matter experiments with topological insulators and superconductors [6]. We believe it would be interesting to explore the interplay of these effects with gravitation in analogue systems, where the corresponding QFTs with non-trivial topologies evolve in real time under a background curved spacetime.
Let us finally note that the essence of most phenomena of QFTs in curved spacetimes is exemplified by studying free quantum fields evolving in a background curved metric, although interesting effects can also arise when exploring such real-time dynamics for interacting quantum fields. In the spirit of quantum simulations [50], one could exploit some of these analogue gravity systems, those in which one can prepare various initial states and measure relevant observables after a controllable real-time evolution, as quantum simulators that address dynamical and nonperturbative problems of QFTs in curved spacetimes, going beyond the capabilities of current numerical simulations based on classical computers.
The goal of our work is to explore a specific system where all of the above points can be addressed. First of all, we will focus on analogue gravity for fermions, exploring the cosmological production of Dirac fermions in an expanding Friedmann-Robertson-Walker spacetime in D = (1 + 1) dimensions. We will show that specific lattice regularizations of these QFTs lead to the aforementioned nontrivial topologies in reciprocal space, which are manifested by the existence of zero modes localized at the spatial boundaries of the expanding spacetime. We show how bulk fermions reproduce exactly the continuum QFT prediction for the production of particleantiparticle pairs as a consequence of the accelerated cosmological expansion. In particular, in the zeromass limit, one reaches a conformally-invariant situation where the fermion production in the bulk is zero. This contrasts with the situation within the spacetime boundaries, where fermions bound to the spatial edges can be created at finite rates despite having strictly zero energy due to a protecting symmetry. We discuss how, by working in conformal time, this phenomenon could be observed in experiments of ultra-cold fermions in Raman optical lattices. This experimental realization brings a very interesting perspective for future work, as one could explore how this real-time dynamics is affected by non-perturbative phenomena such as dynamical mass generation in a Gross-Neveu QFT [51].
This article is organised as follows. In Sec. 2, we revise the theoretical background of particle creation in a (1+1)-dimensional Friedmann-Robertson-Walker universe during a de Sitter expansion phase. This expansion is characterized by an exponentially-growing scale factor. We obtain both analytical and numerical results for the spectrum of created particles, and its number density. In Sec. 3, we start by reviewing some topics related to symmetry-protected topological phases, and their relevance in lattice models Pictorial representation of the phenomenon of particle production in an expanding spacetime. The dynamics of the background metric produces excitations of the field in pairs, which are interpreted as particles (represented in red) and antiparticles (represented in orange). For certain lattice discretizations of the field, the spatial boundaries of the spacetime can actually host zero-energy modes that are exponentially localized to the boundaries, and thus propagate in a lower effective dimension. We explore the phenomenon of particle production for such zero modes, which is a landmark of the interplay of topological vacua and the curved spacetime.
within condensed matter and high-energy physics. We then introduce two different discretization schemes for the theory of Dirac fields in a FRW background, discussing how they deal with the phenomenon of fermion doubling, and how this can affect the description of particle production for periodic boundary conditions. We then move to study the effect of imposing open boundary conditions in Wilson's scheme of discretization, and observe the appearance of topological phases connected to the asymptotic vacua of the lattice field theory. These vacua are characterized by a non-zero topological invariant in reciprocal space, which has a boundary correspondence in the form of zero-modes exponentially localized to the spatial edges of the system. We conclude that these topological modes are also produced as a consequence of the expansion of the universe. In Sec. 4, we propose a detailed experimental scheme for the quantum simulation of this phenomenon in systems of ultra-cold fermionic atoms in Raman optical lattices. Finally, in Sec. 5, we present the conclusions and outlook of the presented results.
Robertson-Walker (FRW) universe. The details on how to describe a QFT of Dirac fermions in curved spacetimes can be found in Appendix A. The FRW metric forms the basis of the standard model of cosmology, capturing the large-scale structure of the universe. The line element is where a(t) is a dimensionless scale factor, and dΣ 2 depends on the metric of the spatial slices, which corresponds to that of a maximally-symmetric manifold [9]. Accordingly, the spatial metric has spherical symmetry, and a uniform curvature that can be either positive, negative, or zero. In FRW spacetimes, the scale factor a(t) determines how big the spatial slices are and, thus, how the universe expands, being its specific time dependence determined by the nature of the stress-energy tensor that sources Einstein's field equations. Assuming that this source corresponds to a perfect fluid, which is homogeneous and isotropic, one can derive the so-called Friedmann equations that determine how a(t) evolves in time. For instance, when the Einstein field equations are sourced by the vacuum energy (i.e. a positive cosmological constant), one obtains an exponentially-growing scale factor, leading to a de Sitter expansion, which will be discussed in more detail below. In D = 1 + 1, the symmetric manifold would correspond to the line dΣ = dx, such that the spatial curvature vanishes and the discussion is, in principle, greatly simplified. However, regardless of the nature of the stress-energy tensor, the Einstein tensor vanishes identically for such reduced dimensionality [52] and, thus, matter cannot act as a source of curvature following Einstein's theory. Moreover, including a positive cosmological constant in Einstein's equations, which leads to the aforementioned de Sitter expansion in D = 3 + 1 dimensions, now implies that the volume element of the metric vanishes, which is incompatible with the desired FRW spacetime (1). The details on how to overcome these problems in D = 1 + 1 can be found in Appendix B.
For the moment, we consider an arbitrary time dependence for the scale factor, and develop the formalism by introducing the conformal time .
Let us also note that, in D = 1 + 1 dimensions, the flat gamma matrices are proportional to the Pauli matrices, such that one can work with two-components spinors as mentioned previously. In this manuscript, we make the following choice In this reduced dimensionality, there is a single Lorentz transformation Λ with a boost of speed v = atanhξ, which acts on the spinors via Eq. (75) with the generators S 01 = −S 10 = σ y /2, and depends on the boost rapidity via Ω 01 = −Ω 10 = ξ. The covariant derivative (76) also depends on these generators via the connection field (77) which, in this case, only depends on conformal time We then substitute all these expressions in the generic action (83), and find that the resulting dynamics amounts to that of Dirac fermions with a timedependent mass in a static Minkowski spacetime with coordinates x = (η, x), and is dictated by the following action Here, we have rescaled the spinor field as ψ(x) → χ(x) = a(η)ψ(x), such that the principle of stationary actions yields the following Dirac equation Therefore, all the effects of the expanding spacetime are encoded in a dynamical multiplicative renormalization of the bare mass, which will depend on the specific time dependence of the scale factor that is itself determined by the underlying classical field equations of the gravity model. Let us note that, regardless of the particular scale factor, one can see that, for conformal QFTs in which the bare fermion mass vanishes, m = 0, there is no effect of the expanding spacetime apart from a trivial rescaling, and thus no particle production.
Before moving on, let us emphasize again that the fermionic fields are not the source of the specific expansion of the scale factor, which should have some other origin as will be discussed in Sec. 2.3. Before moving there, we discuss the physics of particle production for a generic scale factor in the following subsection, recalling that we do not consider any backaction from the field onto the metric, and thus we work in the QFT in curved spacetime formalism. In any case, this effect should be negligible in light of the Fermi-Dirac statistics of the fields.

Bogoliubov transformations and particle creation in asymptotic Minkowski vacua
Let us momentarily go back to QFT in flat spacetimes. Upon quantization, the vacuum of a Dirac field in a static Minkowski spacetime can be uniquely defined. In fact, the vacuum remains the same at any instant of time, such that the notion of particle and antiparticle excitations is unambiguous [1]. Introducing additional interactions in the Dirac action (74) brings in interesting dynamical effects, since these excitations can scatter, and particle-antiparticle pairs can be created from the vacuum. This situation changes for curved spacetimes since, as advanced in the introduction, there are specific situations where the curvature/dynamics of the universe can lead to particle creation [13,14] even in the absence of interactions with other quantum fields. For the problem at hand, where the fermions are coupled to an expanding FRW universe, this occurs when the field absorbs the required energy to create excitations from the gravitational background [26,36,53]. This is particularly transparent in the conformal-time description of the action, Eq. (7), where the effect of the expanding universe is encoded in a time-dependent mass, suggesting that the energy of the Dirac fields shall not be conserved and may result in particle creation.
However, a precise interpretation of particle production in curved spacetimes is complicated by the fact that, in contrast to the Minkowski spacetime, the vacuum is coordinate-dependent and, thus, not uniquely defined [53,54]. Therefore, the typical notion of particles and antiparticles that arises in the canonical quantization of field theories in flat spacetimes [1] must be re-addressed with some care. Let us discuss the canonical quantization of the Dirac theory (7) at a fixed instant of time x = (η ⋆ , x), assuming thus a specific inertial frame [1]. One imposes canonical anti-commutation relations for χ(x) and its conjugate momentum Π χ (x) =χ(x)γ 0 = iχ † (x), which are then upgraded to field operators and denoted using a hat χ(x), Π χ (x) →χ(x),Π χ (x) fulfilling equal-time canonical anti-commutation relations. Since the Dirac equation (8) is linear, we can expand the field operators in the complete basis of its positive-and negativefrequency solutions, so-called mode functions. Accounting for the previous rescaling with the scale factor, we find that the Dirac field at a fixed instant of time can be expanded aŝ Here, the mode functions are obtained by the product of spinor solutions with Fourier components e ±ik x In the expansion (9), we have chosen a normalization that depends on ω k = (k 2 + m 2 a 2 (η ⋆ )) 1/2 , and is consistent with the anti-commutation algebra of the creation-annihilation operators while the remaining anti-commutators vanish. The spinor solutions are normalized according to a standard choice [1], where the dot here represents the matrix-vector multiplication. In addition, they ful- In the case of a flat and static spacetime, the scale factor is trivial a(η) = 1, andũ k (η ⋆ ) = u k e −iω k η⋆ ,ṽ k (η ⋆ ) =ṽ k e +iω k η⋆ , such that the mode functions (10) can be expressed in terms of plane waves e ±ikx , where kx = η µν k µ x ν and the 2momentum is defined on mass shell k = (ω k , k). One can then find the single-particle states |k⟩ as unitary irreducible representations of the Poincaré group x → x ′ = Λx + d. Accordingly, these states transform as |k⟩ → U (Λ, d) |k⟩ = e −ikd |Λk⟩, where the operators yield a representation of the group fulfilling U † (Λ, d)U (Λ, d) = I and U (Λ 1 , d 1 )U (Λ 2 , d 2 ) = U (Λ 1 Λ 2 , Λ 1 d 2 + d 1 ), and are generated by acting with rescaled creation operators on the QFT vacuum We recall that the vacuum fulfillsâ k |0⟩ =b k |0⟩ = 0, and is the only state left invariant under transformations within the Poincaré group [55]. Using this invariance, one can change the inertial frame used to define the canonical momenta at any other instant of time η = η ⋆ , and extend the notion of the vacuum to any other instant of time, arriving in this way to an unambiguous notion of particles and antiparticles. In particular, the modes do not mix under Poincaré transformations and, moreover, the individual number operators can be shown to be Poincaré invariant such that any inertial observer would agree on the specific particle/antiantiparticleparticle content of the state. In fact, Eq. (9) is valid for any instant of time letting x = (η ⋆ , x) → (η, x), which shows that the evolution of the operators is trivialâ k (η) = a k e −iω k η ,b k (η) =b k e −iω k η , and there is no particle/antiparticle creation from the initial vacuum, nor scattering between different particles, unless additional interactions are incorporated in the QFT.
A similar philosophy can be followed for curved spacetimes, although, as advanced previously, more profound conceptual challenges arise. In the case of an expanding spacetime a(η) ̸ = 1, Poincaré invariance is superseded by diffeomorphism invariance, the vacuum becomes coordinate dependent, and particles can no longer be associated with unitary irreducible representations of the Poincaré group [55]. Since there is, in principle, no preferable coordinate system, the notion of vacuum and particle becomes ambiguous for curved spacetimes. As discussed in [54,56], a reasonable approximation to discuss particle production is that of adiabatic vacua, which connect to the so-called Bunch-Davies vacuum in de Sitter spacetimes [57] after an extended period of inflation [53]. In the case of Dirac fields, recursive methods to construct such adiabatic vacua have been recently discussed in [58]. Ultimately, one may adopt an operational philosophy, where the notion of particles/antiparticles is related to specific local detectors that click by absorbing energy from the field, the so-called Unruh-DeWitt detectors [56]. In this article, however, we will be interested in situations amenable to analogue-gravity experiments, in which the simple and unambiguous notions of the vacuum and particles in flat spacetimes are still useful. Paralleling the in-out formalism of scattering in interacting QFTs mentioned above, we consider scale factors that tend adiabatically to constant values in the remote past and distant future. Therefore, the metric (3) tends to asymptotically-flat Minkowski spacetimes, where the vacuum and particle/antiparticle states have a well-defined meaning. Accordingly, we can use the field decomposition in Eq. (10) for those distant times η ≈ η 0 and η ≈ η f , using their corresponding mode functions and creationannihilation operators. However, although customary in flat spacetime, there are some problems with the mode expansion as defined in Eq. (9) when a(η) ̸ = 1. In particular, the above normalization of the modes u k (η) andṽ k (η) would imply that their norm is timedependent, and thus their evolution cannot be unitary. We thus define new modes normalized to 1 as = 0 remain the same. Using these normalized modes, the expansion of the Dirac field (9) becomeŝ where the mode functions are defined in analogy to Eq. (10), and evolve in conformal time following the Dirac equation (8) in momentum space This convention is discussed in more detail in Appendix E. Let us then define the asymptotic vacua |0 0 ⟩ and |0 f ⟩, which are annihilated by the corresponding operatorsâ k (η 0 ),â k (η f ), andb k (η 0 ),b k (η f ). Note that, due to the intermediate expansion, the initial vacuum |0 0 ⟩ may not evolve into the instantaneous groundstate of the QFT at any later time, nor to the asymptotic distant-future vacuum |0 f ⟩ after the complete expansion. In contrast to the flat spacetime, in which the time-evolution of the creation-annihilation operators was a trivial complex phase, the corresponding evolution in the expanding spacetime leads to the following canonical, so-called Bogoliubov, transformationâ where α k (η f ) and β k (η f ) are known as the Bogoliubov coefficients [59,60]. To maintain the anticommutation algebra (11), these dimensionless coefficients fulfill Under such a Bogoliubov transformation, the annihilation operators in the distant future become a linear superposition of both the creation and annihilation operators of the remote past, which is the key to account for particle production. One can readily see that the number of particles in the far future, , upon which all inertial observers agree, is given by N a = n a δ(0). Here, we have divided by the scale factor to take into account the rescaling of the field (9). The divergent factor δ(0) in the particle number corresponds to the infinite spatial volume of the FRW universe, such that the mean density of produced particles, in a fiducial volume cell, is Given that this particle density is obtained by integrating over the spatial momentum, which has inverse length units, whereas the scale factor is dimensionless, one finds that the Bogoliubov coefficient |β k (η f )| 2 is proportional to the number of produced particles for a specific mode. Additionally, according to Eq. (17), the mean density of antiparticles is equal n b = n a , which is a consequence of the fact that the expansion conserves the total charge.

Fermion production for a de Sitter expansion
Once the formalism to understand the phenomenon of particle creation has been discussed, we present the details on how to calculate these quantities for a specific expansion of a FRW spacetime (3). In particular, we consider an exponentially-growing scale factor where we note that the conformal and cosmological times (2) are related via η = −exp{−Ht}/H, such that η ∈ (−∞, 0) for t ∈ (−∞, +∞). This evolution displays a constant rate of expansion a −1 da/dt = H that is commonly known as the Hubble parameter.
Let us now discuss how this specific expansion arises in the two approaches to define a Dirac QFT in a FRW spacetime of D = (1 + 1) dimensions introduced in Sec. 2.1.
First of all, we consider the situation in which this arises as an effective QFT when the fermions are forced to move along a single spatial section of the D = (3 + 1)-dimensional FRW spacetime. The scale factor (20) can be obtained from the standard Einstein field equations (21) which are expressed in terms of the Riemann curvature tensor a positive cosmological constant Λ, and the scalar curvature Using the specific expressions of these quantities for the (3 + 1)-dimensional FRW spacetime [9], one can then derive the so-called Friedmann equations for the scale factor from Eq. (21) which, in this case, have a simple exponentially-growing solution (20) with H = Λ/3. This de Sitter expansion, which is an exact solution of Einstein's equations, is actually used to model the inflationary epoch of the early universe in the standard model of cosmology, focusing on a slowroll regime in which the Hubble parameter is approximately constant [9]. By including a term proportional to the stress-energy tensor in the right-hand side of Eq. (21), one can also obtain other evolutions of the scale factor associated to matter-and radiationdominated universes. Indeed, from this perspective, the positive cosmological constant can be interpreted as the result of a vacuum energy acting as a source of Einstein's equations.
Let us now move to the second alternative, where the Dirac fields move in a (1 + 1)dimensional FRW spacetime that evolves according to JT gravity. As noted above, the problem with Einstein's equations (21) in this reduced dimensionality is that G µν (x) = 0. JT gravity constructs an alternative field equation using directly the curvature scalar. In the presence of a positive cosmological constant, the constant-curvature Jackiw-Teitelboim equation simply reads For the FRW spacetime, the scalar curvature is R = 2 a d 2 a dt 2 , and one can easily obtain the simple exponentially-growing solution (20) with H = Λ/2, which differs from the dependence in higher dimensions. In the context of JT gravity, one can also derive the analogue of the Friedmann equation when a term proportional to the stress-energy scalar is added in the right-hand side of Eq. (24). This leads to other time evolutions of the scale factor a(t) for a matterand radiation-dominated universes [61], which differ from those of Einstein's gravity. In any case, since we are interested in a de Sitter expansion, and both approaches lead to the same exponential growth with the same Hubble constant, any of the interpretations of the origin of the (1 + 1)-dimensional Dirac QFT in the expanding FRW spacetime will be valid. We can thus carry on with the phenomenon of particle production.
To connect with our previous discussion of particle production, this scale factor must be connected to the asymptotic flat spacetimes limits, which requires adiabatically ramping up/down the scale factor. One way of doing this [62] is by means of the following factor of expansion which is used to interpolate smoothly between each of the three following regimes This scale factor is regulated by a parameter λ. The greater λ is, the flatter are the side regions and the better approximated is the de Sitter phase of expansion. We expect that, for λ sufficiently large, the vacuum at η 0 ≪ η in is adiabatically connected to the instantaneous groundstate at η in , such that no excitations are produced during the ramping-up phase. Then, the de Sitter expansion between (η in , η out ) will cause non-adiabatic effects, such that the timeevolved state shall no longer coincide with the instantaneous groundstate at η out , leading to particle production. Conversely, after the de Sitter phase of expansion, the instantaneous groundstate will be adiabatically connected to the vacuum state of the asymptotic future η f ≫ η out , such that no extra particles are produced during the ramping-down phase.
In this way, the particle creation has a well-defined interpretation considering the flat-spacetime asymptotic limits, and is essentially caused by the period of de Sitter expansion which, as will be shown below, admits a closed analytical expression. We note that other hyperbolic-tangent scale factors also allow for closed analytical expressions for the production of Dirac fermion in FRW spacetimes [63][64][65], although the specific expansions cannot be connected to the Einstein or JT gravity field equations sourced by a simple cosmological constant. We now discuss how to calculate the Bogoliubov coefficients for this expansion. For reasons that will become clear when discussing the cold-atom analoguegravity implementation, we follow the diagonalization method [66,67], which uses the instantaneous eigenstates of the single-particle Hamiltonian as the mode functions, and consequently introducing at each time annihilation operatorsâ k (η),b k (η) which are used to define the vacuum state. The details of this method can be found in Appendix E. This approach starts by noticing that the dynamics of the rescaled spinor field in the FRW spacetime can be described by a singleparticle Hamiltonian H k (η) with instantaneous eigenvalues ±ω k (η) = ± k 2 + m 2 a 2 (η), and normalized eigenvectors v ± k (η). These eigenvectors correspond, up to a normalization factor, to the spinor solutions introduced above, Eq. (9), for the specific time in- , which will now be set to the asymptotic remote past η ⋆ = η 0 . At later times, the evolution of the Dirac field is described by Eq. (13), but the normalised mode functions u k (η), v k (η) generally depart from these instantaneous spinor solutions v ± k (η). This departure arises as one leaves the asymptotic remote past η ≈ η in ≫ η 0 and enters in a region with non-adiabatic changes in the scale factor a(η). In this period of expansion, the mode functions must be found by solving the set of coupled ordinary differential equations (ODEs) in Eq. (14), which for each component and in our particular representation reads where u k,1 (u k,2 ) is the upper (lower) spinor component. These ODEs can be expressed as i∂ η u k = H k (η)u k , where the aforementioned single-particle Dirac Hamiltonian reads One finds that the mode solution u k (η) at any instant of time is related to the instantaneous eigenstates of the single-particle Hamiltonian via which is the manifestation of the Bogoliubov transformation in Eqs. (16)- (17) at the level of single-particle solutions. Neglecting the adiabatic changes in the asymptotic regions (η 0 , η in ) and (η out , η 0 ), and assuming a purely de Sitter phase (20), we solve Eq. (27) for u k (η), choosing as the initial condition the instantaneous eigenstate so that β k (η in ) = 0. Then, after the expansion, the Bogoliubov coefficient can be obtained from the overlap of the evolved mode function with the negativeenergy instantaneous eigenstate From this expression, it becomes clear that the production of particles will be negligible if the adiabatic theorem [68] holds, since the mode solutions will remain instantaneous eigenstates of the single particle Hamiltonian at any latter times η > η in , and so the overlap in (32) will be zero. For this particular Hamiltonian and for a purely de Sitter expansion phase, the adiabatic theorem holds as long as the parameters satisfy which is consistent with previously proposed adiabatic parameters [69]. We are thus interested in situations where the expansion does not satisfy condition (33), because it is the non-adiabaticity in the expansion of the universe which induces particle production. As shown in Appendix C, the analytical solution for the mode functions is found by decoupling the ODEs (27) into a pair of Bessel differential equations, whose solutions can be expressed in terms of Hankel functions [70] Here, the four integration constants C i,j are not independent, as we started from two first-order ODEs in Eq. (27). We thus need two initial conditions, which are given by Eq. (31). We note that similar expressions in terms of Hankel functions can be found in the literature for the (3 + 1)dimensional case [71], where differences arise due to the helicity of the spinor solutions, and also for scalar fields [72], where the order of the Hankel functions is real ν ∈ R. As discussed in the Appendix C, our solution in terms of Hankel functions can also be related to previously-found solutions that make use of the lessfamiliar cylinder functions [73]. Before moving on, let us note that this analytical solutions rests on the assumption that no particles will be produced on the adiabatic switching regions, the validity of which will be explored below numerically for specific switchings. Let us now comment on a simple analytical expression for the Bogoliubov coefficient (32), and thus the density of produced particles (19) after an infinitelylong phase of expansion, namely where the only restriction is that of a non-vanishing bare mass m ̸ = 0. The eigenvalues in those limits adopt (34) are determined by the initial condition (31), such that where ) . As discussed in Appendix D, using the asymptotic behaviour of the Hankel functions for z = kη out → 0 − , which assumes the existence of a cutoff for the spatial momenta k ≤ Λ c , we arrive at the result |β k (0 − )| 2 = 1/(e 2πm H + 1), which is reminiscent of a Dirac-Fermi distribution at an effective temperature T = H/2π. Let us note, however, that the instantaneous energy dispersion ω k (η) = (k 2 + m 2 a 2 (η)) 1/2 does not appear in the expression, which would thus yield an infinite density of produced particles when integrated. Altogether, our result is where θ(m) is Heaviside's step function, θ(m) = 1 if m > 0 and zero elsewhere. Note that, in the limit of large masses m ≫ H, particle production is exponentially suppressed.
Let us now assess the validity of these results by taking into account specific parameters for the adiabatic switching regions. We will no longer use the approximations η in → −∞ and η out → 0 − , and so we will not use the limiting forms that we stated before, but calculate numerically the instantaneous eigenstates in the corresponding adiabatic regions using the specific scale factor (25). A numerical benchmark of our result can be found in Fig. 3, where we present the results of a numerical calculation of |β k (η f )| 2 , which involves solving the system of ODEs (27) with the specific scale factor (25), as a function of the bare mass. We consider various finite values of the final expansion time, and fix the parameter of the adiabatic switching to λ = 30, which ensures a smooth and slow connection to the asymptotic Minkowski vacua. We observe that, as η out approaches zero, |β k (η f )| 2 indeed tends to the previous Fermi-Dirac-like function (36), confirming the validity of our analytical treatment and, in particular, the adiabatic switching that leads to a faithful interpretation of the particle production. Let us also note that we have also implemented other switching profiles, which lead to a similar agreement with the analytical prediction. We now explore how particle production depends on the spatial momentum of the fermions. We see in Fig. 2 that, for masses within the sub-Hubble regime (m < H), there is a peak in the spectrum for k ≪ H, and also a non-trivial shape for k ≈ 0. Again, the production is lower for heavier particles, although the distribution is broader. For super-Hubble masses (m > H), we see that the peak occurs at higher spatial momentum k. Analogously, its height lowers and its width broadens as the mass increases, although not as much as in the sub-Hubble case.
Finally, if we integrate this spectrum in momentum space, we can obtain the density of created particles via Eq. (19). If we calculate this integral for different masses, we obtain the plot in Fig. 4. We can see that, for lighter fermions, the total number of produced particles grows as the mass increases (i.e. although the height of the peak of the spectrum lowers as mass increases, it also broadens in such a way that the total area increases, yielding a higher density of produced particles). This occurs until a certain value of the mass is reached, where the broadening of the spectrum is not sufficient to compensate for the lowering of the peak, and the density of produced fermions starts to decrease with the bare mass. The physical reason for this decrease is that the creation of heavy energetic fermions is suppressed as the gravitational background does not have enough energy to produce them. On the other hand, if the mass is very small, we are close to a conformal-invariant expansion where particle production is also suppressed, as displayed in the figure. Altogether, there is a maximum for the density of the produced fermions with an intermediate mass that balances between these two effects. This occurs for fermions with a mass m ≈ 0.274H, which differs with respect to the scalar field case with m = H [72].

Lattice regularization and spacetime boundaries
In this section, we consider two different lattice discretizations of the Dirac QFT in a curved spacetime. In high-energy physics, the lattice is an artificial scaffolding for the fields that serves to regularize the QFT, allowing to treat interacting problems with ultraviolet divergences beyond the perturbative renormalization group [74]. In the context of flat Minkowski spacetimes, lattice field theories (LFTs) are routinely used for this purpose in quantum chromodynamics [75][76][77]. As advanced in the introduction, there are certain lattice discretizations [48,49,78] that display a nontrivial topology in reciprocal space [46,79], and con- nect these QFTs [80][81][82][83][84][85][86][87] to the physics of topological insulators and superconductors in condensed matter [6]. In the condensed-matter context, the lattice is actually physical, and one is not only interested in recovering a continuum limit devoid of lattice artifacts that can appear around certain phase transitions, but also in charting the full phase diagram in which the specific lattice discretization can play a key role.
A celebrated example is that of the quantum anomalous Hall (QAH) effect [88] in both the honeycomb [89,90], and square [91] lattices. These models can be connected to a Hamiltonian formulation of Dirac fermions in (2+1) dimensions with a Wilson-type discretization [78], and host groundstates that cannot be understood from the paradigm of Landau's theory of spontaneous symmetry breaking [92,93]. In fact, characterizing these states requires introducing an invariant known as the Chern number, which characterizes the topology in reciprocal space, and is responsible for the robustness of the quantized transverse conductance [94,95]. The QAH phases are a specific type of the so-called symmetryprotected topological (SPT) phases [96]. In general, different SPT phases can be found within the same symmetry class, and cannot thus be described by a local order parameter and connected via a symmetrybreaking mechanism. Characterizing these phases and the phase transitions requires instead the use of other topological invariants that can only change via a gap-closing phase transition that is not associated to symmetry breaking. For fermionic models, there are various SPT phases for different symmetry classes and spatial dimensionalities [97,97,98], which can also be connected to discretizations of Dirac QFTs [99], and turn out to be robust to perturbations that do not explicitly break the specific symmetry.
In (1+1) dimensions, an archetype of SPT physics is the so-called Su-Schrieffer-Heeger model in the limit of a classical lattice dimerization [100][101][102]. Here, one can distinguish a topological phase in the symmetry class BDI from a trivial band insulator using a topological invariant known as the Zak's phase φ Zak [103], which is defined as the integral of the Berry connection A n (k) [104,105] in reciprocal space. The model allows for a topological phase transition as one changes the microscopic couplings, through which the energy gap to excitations vanishes, and the Zak's phase changes abruptly φ Zak : 0 → π. This marks the onset of topological effects, which have a boundary manifestation: the existence of charge fractionalization in zero-energy excitations that are exponentially localized to the edges of the lattice, i.e. boundary zero modes. If the lattice dimerization is dynamical, such boundary modes can be localized to solitons interpolating between two different groundstates, paralleling the Jackiw-Rebbi mechanism of charge fractionalization in QFTs [106].
Although the Su-Schrieffer-Heeger model can be rigorously connected to a specific regularization of a continuum Dirac QFT with boundary zero modes [107,108], it is not a standard discretization in the Hamiltonian formulation of LFTs [109,110]. As discussed in [83], it can be related to the aforementioned Wilson discretization [78] for a specific choice of gamma matrices and microscopic parameters. This Wilson discretization can also be depicted as a lattice model of fermions hopping on a crosslink two-leg ladder and subjected to an external πflux, the so-called Creutz ladder [111,112]. This model can also host boundary zero modes exponentially localized to the left-and right-most boundaries of the ladder, signaling the occurrence of an SPT phase in the symmetry class AIII [113][114][115]. The role of finite temperatures [116,117], dynamical quenches [118,119], and charge pumping [120] has also been discussed. The interplay of interactions and topology has been explored recently, including repulsive and attractive Hubbard-type interactions [82,83,85,[121][122][123], many-body [124] and disorder [125,126] induced localization, and interactions mediated by a discrete gauge field [127]. Additionally, by exploring regimes away from the external π-flux limit, the continuum limit connects to Lorentz-violating QFTs, allowing for the characterization of topological phenomena via persistent groundstate currents [128].
Let us note that all of the above phenomena correspond to discretizations in a flat Minkowski spacetime. To the best of our knowledge, the introduction of a curved metric remains largely unexplored in the context of SPT phases. Consequently, it is very natural to ask ourselves if the topological properties can have an interplay with the typical phenomena studied within the realm of QFTs in curved spacetimes. This is not only interesting from the theoretical perspective, but might also be fruitful for a dynamical manifestation of topological effects in possible quantum simulation experiments, as discussed below.
The most intuitive way to discretize this theory is by defining a chain of N evenly-spaced sites, separated by a length a, and defining fermionic creation and annihilation operators on those sites. The spatial derivatives appearing in the Hamiltonian must be replaced by a finite difference, such that The discretized Hamiltonian is then obtained by direct substitution on Eq. (38), and by discretizing the spatial integral where we have introduced the chiral gamma matrix which is anti-Hermitian for the mostly-plus metric. In LFTs, one usually confines the fields in a box, and then allows its size to diverge. Since one expects the fields' amplitudes to decay sufficiently fast, boundary effects are usually neglected, and so periodic boundary conditions (PBC) and a basis of plane waves are generally used. We follow this approach now, and relegate the study of other boundary conditions to the following subsection, which will allow us to study boundary effects related to the aforementioned SPT phases of matter. For now, we assume PBC, and useχ n = 1 √ N kχ k e ikan . Since this is a discrete theory with spatial periodicity a, only those values of the crystal momentum within the Brillouin zone k ∈ BZ = {− π a + 2πn N a , n ∈ Z N } will be allowed, leading to the Hamiltonian in momentum spacê This way of discretizing the theory is usually known as the naïve discretization, which is afflicted by the socalled fermion doubling [76,77,110]. In d = 1 spatial dimension, this implies that, when taking the continuum limit a → 0, one recovers twice as many fermions as there were in the original continuous theory. In general, for Hamiltonian field theories in D = d + 1 dimensional spacetimes, one would encounter 2 d doublers. This presents a problem when one is interested in particle production, since it means that the discrete theory has additional low frequency excitations. For our particular model (45), this will result in an overproduction of particles. In order to prove that, we need to calculate the Bogoliubov coefficients β k (η), noting that the previous system of differential equations (27) gets modified due to the discretization The instantaneous eigenvectors v ± k (η), which we used to impose the initial condition (31) and calculate the final density of produced particles (32), are also modified, since the single-particle Dirac Hamiltonian (28) changes due to the discretization. This can be directly read from Eq. (45), yielding for our representation With these new equations of motion and instantaneous eigenstates, we calculate numerically the fermion production. One would expect to recover a good approximation to the continuum results for k ≪ 1/a, where the dispersion relation becomes similar to that of Dirac fermions. As the momentum increases, however, the differences between the continuum equations and the naïvely-discretized ones becomes more important, until reaching the edge of the Brillouin zone k ≈ π/a, where the fermion doubler lies. Since the dispersion relation there is again similar to that of a Dirac fermion, a high contribution to the production of particles can again take place. As can be seen in Fig. 5, the spectrum of particle production at small momenta (solid line) reproduces accurately the continuum result (dots), but we observe spurious creation of particles caused by the fermion doubling when approaching the boundary of the Brillouin zone. One can clearly see that the distortion of the particle-production spectrum is fully symmetric, as both low-energy fermions are equally affected by the time-dependent mass.
Let us now discuss how the total number of produced fermions in this discretized theory is related to that in the continuum theory. Let us recall that, in the original continuum QFT, the spectra in Figs.  2 (a) and (b) peak at a certain momentum k that decreases when the bare mass m is lowered with respect to the Hubble parameter H. If such peaks correspond to the region in the Brillouin zone in which the dispersion relation of the discretized model faithfully approximates the continuum one, we expect to get twice the density of produced particles. We can check this numerically by calculating the continuum density of produced particles, and then comparing that value with the one obtained after the discretization is done. We show the results of doing this in Fig.  6, where the ratio between the continuum density of created particles and the one obtained after the dis- cretization n Naive n Cont is presented. As shown in this figure, as aH = H/Λ c ≪ 1, which corresponds to a regime where the peak of particle production lies well-below the lattice cutoff Λ c = 1/a, the total density of produced fermions in the discretized model is twice the one estimated from the continuum QFT.
There are several proposals to get rid of these doublers. For example, in the staggered fermion approach, due to Kogut and Susskind [109], it is proposed to reduce the number of degrees of freedom by using a single component field in each site of the lattice, which halves the doublers. An alternative that makes connection to topological phases of matter is the so-called Wilson's fermion approach [129], which we develop below in the context of gravitational particle creation.

Wilson fermions and reciprocal-space topology
The idea proposed by Wilson [129] is to include an additional term in the Hamiltonian that acts as a momentum-dependent mass, known as the Wilson mass. This mass sends all of the spurious doublers to the cutoff of the theory, apparently removing all of their effects in the long wavelength properties of the continuum limit. The objective is to leave the mass equal to the bare mass around k = 0, but making it very heavy around the edges of the Brillouin zone, where the doublers lie. The Wilson term that must be added to the naïve Hamiltonian (43) iŝ where r is the so-called Wilson parameter, a dimensionless parameter that is typically set to r = 1, although it can take other real values for more generality. This term is the discretized version of a second derivative, which is a priori irrelevant in a renormalization-group (RG) sense for the continuum QFT [129]. However, allowing for negative bare masses can change the physics considerably, as it can lead to topological phases with boundary zero modes that would not be affected by the coarse-graining and rescaling of the RG, and also a non-zero topological invariant that is preserved under the RG flow [85]. Let us discuss these properties by considering the total HamiltonianĤ =Ĥ N +Ĥ W .
Again, working with PBC and a basis of plane waves, the new single-particle Hamiltonian is which corresponds to a substitution of the bare mass term for a momentum-dependent one, the so-called Wilson mass This momentum-dependent mass behaves as it should: one recovers the original mass around k ≈ 0, whereas the mass of the doubler becomes very heavy in the continuum limit a → 0 When calculating the particle-production spectrum of the discrete theory in Wilson's approach, we use this new single-particle Hamiltonian to obtain the instantaneous eigenstates, and modify the equations of motion (46) by substituting the bare mass with the Wilson one ma(η) → m W (k, η). The numerical solution leads to the spectral distribution of Fig. 7, where the solid line corresponds to Wilson's discretization and the dots stand for the continuum predictions. Comparing this figure with the naïve-fermion case in Fig. 5, we see that the peak of the contribution of the doubler on the edge of the Brillouin zone has become negligible. For the parameters chosen, this is consistent with our previous results, as we concluded in Sec. 2.3 that production of very heavy particles is highly suppressed. Accordingly, in order to recover the correct continuum results, we need to impose aH ≪ 1 such that the dynamical change of the bare mass does not interfere with Wilson's prescription: the doublers must remain at the cutoff of the QFT. We must consider the combined effects of the modification of the dispersion relation away from the doublers, which can be neglected by choosing values of m and a such that the tail of the spectrum of created particles is negligible around k ≈ ± π 2a , and the effect of the new Wilson mass term around the center of the BZ, whose correction to the original mass is of order O(k 2 a 2 ). Once more, we should choose small values of the mass m, so that its spectrum presents negligible amplitudes for large values of the momentum, and values of a which enlarge the Brillouin zone, making the correction of the mass term smaller. In Fig. 8 we represent the ratio between the continuum density of created particles and the one obtained after the Wilson's discretization. We can see that the regime aH = H/Λ C ≪ 1 corresponds to a situation where the peak of production lies well-below the lattice cutoff Λ c = 1/a, and one recovers the continuum results without the effects of the fermion doubler.

Topological fermion production at the boundary
Now that we have our discrete Hamiltonian field theory in Eqs. (43) and (48), we discuss how its vacuum is related to zero-temperature SPT phases, and how this can modify the spectrum of particle production depending on the boundary conditions. We stated above that the characterization of SPT phases is based on topological invariants which, in the (1+1)-dimensional case correspond to the Zak's phase φ Zak [103]. This topological invariant is defined as the integral of the Berry connection [104,105] over the Brillouin zone, φ Zak (η) = BZ dkA k (η). The Berry connection is defined as where we recall that v − k (η) is the single-particle negative-energy instantaneous eigenstate. Since we have a dynamical mass in our problem, the dependence with η is here treated as parametric, and we consider the topological invariant associated to each one of these instantaneous groundstates. The Zak's phase will differentiate SPT phases ( φ Zak π ∈ Z) from topologically-trivial ones ( φ Zak 2π ∈ Z), which can be rephrased in terms of a Z 2 -valued gauge-invariant Wilson loop W Zak = e iφ Zak ∈ {−1, +1}. In Wilson's approach to the discretized Dirac QFT [83], the Zak's phase is given by which amounts to a sign difference of the mass of the Dirac fermions, those at the center and edges of the BZ. Accordingly, the instantaneous groundstate displays a non-trivial Zak's phase φ Zak (η) = ±π when the Dirac fermion and its spatial doubler have opposite masses. Since the expansion of the spacetime is embodied in a time-dependent bare mass, there may be certain parameters of the theory for which the evolution of the spacetime itself will induce a topological phase transition of the instantaneous groundstates.
As an example, let us suppose a bare mass m = −2/a, and a Wilson parameter that will henceforth be set to r = 1. We also consider a Hubble constant H = 1/a, and assume that the spacetime expansion begins at η in = −2a and finishes at η out = −0.5a. We consider a purely de Sitter expansion a(η) = −1/Hη, neglecting the asymptotic flat regions for this example. Then, at η = η in , we have m W (0, η in ) = −1/a and m W (π/a, η in ) = +1/a, i.e. the sign of the masses is opposite and, thus, we are initially in a SPT phase with φ Zak = −π. At η = η out , after the spacetime has undergone an expansion, the instantaneous Wilson masses evolve into m W (0, η out ) = − 4 a and m W (π/a, η out ) = − 2 a , acquiring the same sign, such that the Zak's phase vanishes in this case φ Zak (η out ) = 0. Thus, the system would start in a topological phase and end in a trivial one, which can only occur through an intermediate gap-vanishing phase transition. This situation reminds of the adiabatic dynamical quenches mentioned in the introduction [118,119], where a quantum system crosses a critical point by the external modification of a microscopic parameter. Interestingly, gravity is responsible for such external modification of the parameters, leading to a topological phase transition via the de Sitter expansion. Note that, due to the Kibble-Zurek mechanism [130][131][132], the crossing of the phase transition requires a breakdown of the adiabatic approximation, and can lead to excitations that are not connected to the particle production of the continuous QFT described in the previous section. We will thus avoid this situation and explore those parameters for which the instantaneous groundstates remain in a SPT phase during the whole de Sitter expansion. The whole particleantiparticle production will thus be a consequence of the breakdown of Poincaré invariance in the expanding spacetime, and the change of the notion of vacua in the asymptotically-flat spacetimes.
So far, we have only looked at bulk properties by imposing PBC. However, we know that SPT phases have a bulk-boundary correspondence [6] manifested in the appearance of zero-energy states, which are exponentially localized to the boundaries of the chain, i.e. to the spatial boundaries of spacetime. The first step to do this is to change our boundary conditions from periodic (PBC) to open (OBC). As a consequence, we will no longer be able to use momentum as a good quantum number (49). In turn, our strategy will be to directly diagonalize the total Hamiltonian in position space. It is very intuitive to understand this Hamiltonian as a tight-binding-like model [133], with tunnelings and self-energies shown in Figs. 9 (a)  and (b).
Thus, when looking for a single-particle Hamiltonian with OBC, we can assume a chain of N sites and spacing a, and express this Hamiltonian as a 2N × 2N matrix H n,n ′ (η) considering also the two internal degrees of freedom of the Dirac spinor. Let us note that this matrix is block-tridiagonal, having 2 × 2 matrices in the diagonal that depend on the bare and Wilson mass discretization. The upper and lower blocks correspond to 2 × 2 matrices that include the naïve and Wilson tunnelings, bringing the singleparticle excitation to the nearest-neighboring sites. Then, the procedure is to diagonalize this matrix, obtaining N pairs of positive-and negative-semidefinite eigenenergies, ω j,± (η) with their associated eigenvectors, v ± j (η), where j is an index that labels the eigenvalues, and plays the role of the spatial momentum in the translationally-invariant situation for PBC. After this, the procedure to obtain the Bogoliubov coefficient β j is analogous to that followed in previously. The system of coupled ODEs (46) now reads which must be solved numerically after imposing an initial condition analogous to Eq. (31). In our case, we consider ∀j ∈ {1, · · · , N } : ω j+ ≥ 0. This amounts to an initial situation in which the asymptotic Dirac sea is obtained by filling all the negative-energy solutions and only one of the two topological states of zero en- ergy, which is localized to one edge of the chain. After the de Sitter expansion, we calculate the Bogoliubov coefficient by the analogue of Eq. (32), but the matrix vector products are now ∀j ∈ {1, · · · , N } : ω j− ≤ 0. This is to be interpreted as the overlap between the evolved states and the eigenstates of the final Hamiltonian. For the propagating modes, this can be understood as an excitation to positive-energy eigenstates, and for the topological modes, as a directed flow from the edge that was initially populated towards the other edge of the chain or bulk states. We thus recover the particle-production spectrum by associating each |β j (η out )| 2 parameter with its corresponding energy ω j+ (η out ). Accordingly, the particle-production spectrum will no loner depend on momentum, but rather on the energy of the particles and antiparticles in the asymptotically-flat spacetimes. Although conceptually equivalent, the process with OBC is more involved, as it requires solving N coupled ODEs for each different value of j, so a full spectrum calculation requires solving N coupled ODEs N times. This difficulty is due to the lost of periodicity, but it is worthy since we can now look for dynamical manifestations of the SPT phases due to the gravitational production of particles.
Since there are zero-energy states exponentially localized to the spatial edges in the SPT phase, we are interested in distinguishing bulk and boundary contributions to the production caused by the de Sitter expansion. This is automatically allowed by the previous spectrum, as the boundary contribution stemming from the edge states must appear at ω j+ = 0 energy. If we solve numerically the system of ODEs, and perform the corresponding numerical diagonalizations, we obtain the result shown in Fig. 10. In this figure, we see that there is in fact particle creation for this topological zero-modes when we impose OBC, which becomes manifest via a non-zero fermion number inside the energy gap of the spectrum. This is an interesting result that could be observed in a quantum simulation experiment, as discussed in the following section. Since the production of zero-modes is intensive (as can be seen in Fig. 10, it is inside the energy gap due to the mass), its contribution to the total density of particles (19) will be negligible in comparison to the contribution of the bulk, which is extensive. On the other hand, if one has frequency resolution, or spatial resolution to localize a probe to the boundaries, the effect of these topological edge states should be distinguishable from the bulk, as we have just shown. It is interesting to highlight that, although there is no particle production for zero-energy bulk propagating modes, as this would require a massless conformallyinvariant limit, the edge states change this paradigm allowing for zero-energy particle production.

Analogue gravity in ultracold Fermi gases
In the previous sections, we have discussed the phenomenon of fermion production during a de Sitter phase of expansion by solving the real-time dynamics of Dirac fermions in a D = (1 + 1)-dimensional Friedmann-Robertson-Walker spacetime. After that, we have regularized the problem on a lattice, which has allowed us to change the topology of the basis domain of the fields from R × R to R × I, with I a finite subinterval of the real line associated to open boundary conditions. In this case, topological effects in reciprocal space have a boundary manifestation in the form of zero-energy modes exponentially localized to the boundaries of I. We have demonstrated that this topological modes, although having zero energy, can accommodate for particles-antiparticles being created by the de Sitter expansion of the universe, which contrasts with the case of zero-energy massless bulk modes. In this section, we discuss a path for the quantum simulation of the phenomenon of fermion production in expanding universes using table-top experiments, proposing an experimental scheme that employs ultra-cold atomic gases in optical lattices [35,[134][135][136]. This would allow us to test the results presented in the previous sections and, more interestingly, would open new perspectives for the study of non-perturbative effects in real-time dynamics, for example by adding a four-Fermi interaction in the form of a Gross-Neveu-type model which, among other interesting phenomena, can lead to chiral symmetry breaking in this real-time scenario Before presenting the details of the scheme, let us note that previous proposals for the quantum simulations of Dirac QFTs in a curved spacetime can be found in the literature, e.g. [137][138][139][140][141]. In this case, there has been a certain focus in spacetimes of reduced dimensionality with spatial inhomogeneities, such as Rindler ones, where a generic formulation based on the previous formalism of vielbeins and the spin connections generally require to implement an inhomogeneous and non-unitary tunneling of the fermionic atoms in the optical lattice. For the specific case of the FRW spacetimes considered in our work, working with conformal time simplifies things considerably, as the inhomogeneities only occur in the temporal direction, and shall amount to a specific real-time modulation of the experimental parameters. In this sense, our scheme can directly exploit the progress in the quantum simulation of lattice field theories in flat spacetimes reviewed in [142][143][144][145][146][147], alleviating some of the difficulties associated to the quantum simulation of more generic curved spacetimes. In particular, we will show that the schemes presented in [128,148,149], which are based on the idea of Raman optical lattices for implementing an effective spin-orbit coupling [150][151][152][153], can be modified minimally, including a real-time modulation of a single experimental parameter, to allow for a quantum simulation of Dirac fields under the de Sitter expansion. For the sake of completeness, we present a selfcontained discussion of the various ingredients of this proposal, emphasizing where the differences with respect to [128,148,149] would arise.

Raman optical lattices and expanding spacetimes
Our objective is to simulate the dynamics of the discretized system shown in Figs. 9 (a) and (b). To do so, we consider a gas of fermionic atoms, such as the alkaline-earth 87 Sr atoms [154]. In this case, the total electronic angular and spin momentum vanishes in the ground-state manifold, which is composed of the Zeeman sub-levels associated to the nuclear spin F = I = 9/2, so there are 10 Zeeman sub-levels M F ∈ {−9/2, −7/2, · · · , 9/2}, which can be split by applying a weak magnetic field B ex = Be z (we choose e z as the quantization direction). Since we are interested in simulating the two spinor states of the Dirac field, we focus only on two of those Zeeman sub-levels, which we shall denote |↑⟩ = | 1 S 0 , F, M ↑ ⟩ and |↓⟩ = | 1 S 0 , F, M ↓ ⟩. Note that one must choose M ↑ and M ↓ such that the selection rules allow for two-photon Raman transitions between them. The interest of working with these atomic species is that, in addition to the internal SU (N ) symmetry of their scattering [154], they have ultra-narrow optical transitions that allow to minimize the residual photon scattering associated to these Raman transitions.
This atomic gas is under the influence of three counter-propagating laser beams, as depicted in Fig.  11. This set-up generates a blue-detuned threedimensional optical lattice [35], with an optical potential of the form Here, j ∈ {1, 2, 3} ≡ {x, y, z} denote each one of the spatial directions, and k j = k j e j is the wave-vector of the laser beams, with mutually-orthogonal polarizations ε j , and V 0,j is the amplitude of the ac-Stark shift experienced by the states in the groundstate manifold | 1 S 0 , F, M F ⟩. We choose V 0,1 ≪ V 0,2 , V 0,3 , so that the dynamics are effectively frozen in the directions y and z, and the system can simulate our (1 + 1)dimensional original problem. We also consider the use of an additional laser beam in a traveling-wave configuration to drive two-photon Raman transitions between |↑⟩ and |↓⟩ by means of off-resonant couplings to states within an excitedstate manifold |e⟩ ∈ {| 3 P 1 , F ′ , M ′ F ⟩}, and will thus be referred to as the Raman beam. The wavevector of this extra laser beam will be denoted byk 4 , and its polarization byε 4 . We choose them to satisfỹ k 4 · k 1 = 0 andε 4 · ε 1 = 0 (e.g. k 1 = k 1 e 1 , k 4 = k 4 e 2 ), adjusting the polarizations in a way that the respective selection rules allow for two-photon transitions between the two Zeeman sub-levels. By virtually populating the excited states {|e⟩}, transitions between |↑⟩ and |↓⟩ can occur, which will be used to simulate the spin-dependent tunnelings of Figs. 9 (a) and  (b). We consider large detunings from these excited states which, in addition to the narrow linewidth of these transitions, allow to minimize the heating mechanism due to the residual spontaneous photon emission. Moreover, the additional ac-Stark effect induces a non-linear shift with respect to the Zeeman quantum number, which can be exploited to set the laser-beam frequencies such that the two-photon processes only involve the states |↑⟩ and |↓⟩ [155], provided that the Fermi gas has been polarized to one of those states using a preliminary stage of optical pumping. We note that, due to the traveling-wave configuration of the Raman beam, the two-photon processes that absorb a photon from the weaker standing wave and emit it in the Raman beam (and viceversa), lead to an optical Raman potential [150] with a period that is twice the one of the optical-lattice potential cos(k 1 r 1 )e i(k4r2−∆ωt) σ + + h.c., (58) whereṼ 0 is the amplitude of the two-photon Raman transition, and we have introduced σ + = |↑⟩ ⟨↓|, and Figure 11: Scheme for the quantum simulation with ultra-cold atoms in an optical Raman potential: (a) A cloud of atoms is subjected to three pairs of counter-propagating laser beams in a standing-wave configuration with mutually orthogonal linear polarizations (blue arrows), leading to a cubic optical lattice. The optical-potential depths in the y and z directions are much larger than in the x direction, leading to effective one-dimensional dynamics. A fourth laser beam in a traveling-wave configuration (orange arrow) is used to induce two-photon Raman transitions in combination with the standing wave along the x axis, which require their beatnote to be tuned close to the resonance δω ≈ E ↑ − E ↓ , as depicted in (b). In this way, the Raman terms oscillate with twice the period in comparison to the ac-Stark shifts that lead to the optical lattice. The expanding background is encoded in the time-dependence of the Raman detuning ∆ω = ω 4 − ω 1 . When tuned close to the resonance, i.e. ∆ω ≈ ω 0 = E ↑ −E ↓ , this term can drive the aforementioned spin-flip transitions by virtually populating the excited state. This process involves absorbing a photon from the traveling-wave Raman beam and re-emitting it to the standing wave, while simultaneously exciting the atom. However, since the period of the Raman potential (58) is exactly twice the one of the optical-lattice potential (57), atoms standing in the minima of the latter see a vanishing Raman-beam intensity and, therefore, no local spin-flips are driven. Instead, only spin-flipping tunnelings are induced by the Raman potential, which can be used to generate spin-orbit coupling [150][151][152][153], as demonstrated in landmark experiments with bosonic and fermionic gases [155][156][157][158]. In [128,148,149], slight modifications of this Raman lattice scheme were considered for the quantum simulation of relativistic Dirac QFTs in flat spacetimes with four-Fermi interactions. Let us now discuss how, using conformal time (3), the quantum simulation of Dirac fields in an FRW spacetime describing an expanding universe can also be realised with specific modifications of the scheme. The Hamiltonian field theory of the 87 Sr Fermi gas reads, in second quantization [136], as followŝ whereΦ † σ (r) andΦ σ (r) are fermionic field operators that create/annihilate an atom at position r in one of the two selected Zeeman sub-levels σ ∈ {↑, ↓}, V (r) = V latt (r)+V Ram (r), m is the mass of the atoms, and a s is the s-wave scattering length for collisions, which can be controlled via optical Feshbach resonances [154]. Although these can interfere with the SU(N ) symmetry of interactions, we note that only a couple of levels is required in this work, so there is no problem in controlling the scattering by admixing with excited states that have a non-vanishing electronic angular momentum.
We want to map this Hamiltonian to the lattice model in Eqs. (43) and (48). First, in order to obtain a lattice description of the ultra-cold atomic model, we consider the regime of deep optical lattices by imposing |V 0,j | ≫ E R,j , where E R,j = k 2 j /2m is the recoil energy. In this limit, atoms are tightly confined to the minima of the cubic optical lattice, which forms a periodic crystal at r 0 n = j λj 2 (n j + 1/2)e j , with n j ∈ Z Nj , and λ j = 2π/k j being the laser wavelengths. A better description of this lattice model is obtained by working in the so-called Wannier basis [159,160], which leads to the transformation Here, w(r − r 0 n ) represents the Wannier function localized around r 0 n , andf n,σ ,f † n,σ are dimensionless creation-annihilation operators on the corresponding lattice site, which satisfy a fermionic algebra {f n,σ ,f † n ′ ,σ ′ } = δ n,n ′ δ σ,σ ′ . Using this basis, the Hamiltonian (59) can be split into a spin-conserving term and a spin-flipping one. Since the resulting microscopic parameters will be expressed as integrals of the Wannier functions, which are tightly confined to the lattice sites, we can treat the problem as a lattice model with only nearest-neighbor couplings. The spin-conserving terms of the Hamiltonian read wheren n,σ =f † n,σfn,σ is the number operator, and we have introduced t j as the tunneling strength along the e j direction, and U as the Hubbard interaction. The explicit expression of these parameters is obtained by performing the corresponding integrals of Wannier functions [35] and, assuming λ j = λ = 2π k ∀j, they read As aforementioned, by setting the optical-lattice depths V 0, 1 ≫ V 0, 2 , V 0, 3 , the tunnelings along the y and z directions become negligible compared to that along the x direction, t 1 ≫ t 2 , t 3 . Moreover, we will also consider that t 1 ≫ U , and neglect the effect of the Hubbard interactions. In the following, we will show how the effect of the additional Raman potential can modify the (1+1)-dimensional tight-binding model in a way that connects to the two different discretizations of Dirac fields in FRW spacetimes discussed in the previous section. We now discuss how the spin-flipping tunnelings of Fig. 9 (a) can be mediated by the Raman potential. One proceeds with the second-quantized Raman term (59) in a similar way, expanding the fields in the Wannier basis (9), and performing the corresponding overlap integrals. As advanced previously, since the Raman potential vanishes at each minima of the lattice, the contributions to local spin-flipping terms vanish. To be more precise, the doubled period of the Raman potential makes the integrand of the corresponding Wannier integral to be an odd function over a symmetric interval of integration, which must thus vanish. The leading-order contributions are then nearest-neighbor laser-assisted tunnelings of strength t along the direction of the standing wave which, simultaneously, change the Zeeman sub-level. Let us now discuss the connection to the two possible discretizations: (i) Wilson-fermion scheme.-In this case, we also need the spin-dependent tunnelings of Fig. 9 (b). The idea is to allow for a detuning δ in the Raman beam, such that it drives slightly off-resonant two-photon transitions between |↑⟩ and |↓⟩, such that the beatnote frequency in Eq. (58) is set to with δ ≪ ω 0 . Moving into a rotating frame, the detuned Raman potential contributes to the lattice Hamiltonian witĥ In this expression, the approximate form oft after applying a Gaussian approximation around the minima of each optical-lattice minima reads The missing step is that, after a U (2) gauge transformation and a rescaling to obtain the correct units for the lattice field operators (41), one can map the cold-atom creation-annihilation operators to those of the Dirac spinor field where we have used the notationχ n1 = (χ n1,u ,χ n1,d ) t for the lattice spinor field. After this transformation, one can see that the site-dependent phase of the Raman tunneling (65) disappears, and we get exactly the nearest-neighbor tunneling used for the discretization of the Dirac kinetic term (43). In addition, this transformation also affects the spin-conserving tunneling in Eq. (61), turning it into a spin-dependent tunneling that can be mapped exactly onto the Wilson mass term (48) when considering also the Raman detuning. In summary, we recover the lattice field theory in Eqs. (43) and (48) with the correspondence between parameters In comparison to previous schemes for the quantum simulation of Dirac fields in flat spacetimes [128,148,149], we see that the required ingredients when using conformal time for the Dirac fields under a FRW spacetime are exactly the same. One of the differences of the mapping is that the transformation in Eq. (67) has been modified with respect to those of previous works [128,148,149], which is a consequence of the different convention of the metric signature. Additionally, this choice also changes the sign of the detuning in Eq.(68), which will require using Raman beams that are blue-detuned with respect to the transition between the Zeeman sub-levels. Finally, the most important difference with respect to the quantum simulation of Dirac fields in flat spacetimes [128,148,149] is that simulating the expansion of the FRW universe requires using a time-dependent detuning of the Raman beam δ → δ(η).
-Contrary to what would be expected, finding a quantum simulation scheme for the naïve-fermion discretization (43) requires additional experimental complexity in comparison to the Wilson-fermion one. First of all, the spindependent tunneling that lead to the Wilson mass is no longer required. This tunneling can be inhibited [161] by exploiting a linear gradient of the onsite energies which can be achieved by lattice acceleration, requiring a linear drift of the optical-lattice beams detuning with time [162], or the application of a magnetic-field gradient. In both cases, the effective Hubbard model in Eq. (61) receives a correction where δω n1,σ = n 1 ∆ σ , and ∆ σ is the aforementioned gradient that can depend on the internal state if it arises from a magnetic field. Provided that |t 1 | ≪ ∆ σ , the tunneling mediated by the standing wave becomes energetically penalized, and can be neglected up to leading order. In general, this gradient can also inhibit the Raman-mediated tunneling. However, one can modify the Raman-beam frequency in Eq. (64), such that the Raman beams provide the required energy to overcome the gradient penalty during the tunneling. In particular, if one considers a spin-independent gradient ∆ ↑ = ∆ ↓ , and imposes one obtains a direct mapping to the naïve-fermion discretization (43) with parameters

Conformal time in the laboratory
In this subsection, we start with a small digression to emphasize the simplifications that arise from using conformal time to describe the properties of the Dirac field, specially in light of the requirements for its quantum simulation. Since quantum simulators require a Hamiltonian formulation, one must be aware that the canonical quantization of Dirac fields in nonstatic curved spacetimes can present important subtleties [163]. In particular, if we follow the standard quantization route discussed for conformal time below Eq. (37), but considering now the cosmological time t (2), we would get a canonical momentum Π ψ (x) = − √ −gψ(x)γ 0 (x) = ia(t)ψ † (x) and, upon quantization of the field and its canonical momentum, arrive directly to the Hamiltonian field theorŷ In this expression, one clearly notices that the last term in brackets leads to a non-Hermitian operator. As discussed in [163], this non-Hermitian contribution is generated precisely by the spin connection (82), and is a generic consequence of the covariant derivative of fermionic problems in non-static metrics such as the FRW spacetime. Accordingly, one should be more careful in defining a correct process of canonical quantization, which entails a rescaling of the above field operatorsψ(x),ψ † (x) [138]. This leads to more complicated Hamiltonians that, upon a naive-fermion discretization similar to Eq. (43), are described by tight-binding models with both spin-conserving and spin-flipping tunnelings. The strengths of these tunnelings will generally depend on the spacetime coordinates [138]. In our specific situation, the spinconserving and spin-flipping tunnelings would depend on the scale factor a(t), and thus become timedependent.
On a technical level, one would have to implement the specific time dependence on both type of tunnelings by modulating the intensity of the opticaland Raman-lattice potentials, which contrast with the simpler modulation of the detuning that is required when working with conformal time (70). On a more fundamental level, we see that the naïve discretization using cosmological time already requires a combination of both spin-conserving and spin-flipping tunnelings, whereas only spin-flipping tunnelings where required in the conformal-time case. (43). As a consequence, it is not clear how one would proceed to get a Wilson-fermion discretization which, in the conformal-time case, exploited a momentumdependent mass term (50) that comes from spinconserving tunneling processes as well. Altogether, performing a cosmological-time quantum simulation of fermion production in the boundary of the FRW spacetime would considerably increase the complexity of the scheme and, quite likely, forbid a direct application of the Raman-lattice toolbox that is currently being used in several experiments [155][156][157][158].
According to our proposal, when performing the quantum simulation, we should interpret the real time of the experiment as representing the conformal time. We recall that the conformal time is negative for the de Sitter expansion (−∞ < η < 0), but the lapse between an initial instant η 0 and a final one η f is actually positive. We should then identify the real time of the experiment, which starts as t = 0, with the conformaltime lapse t = η − η 0 during which we want to simulate a period of expansion of the universe, leading to a bare mass that increases from its initial value ma(η 0 ) to ma(η f ). This mass is controlled by the detuning of the Raman beam, and its explicit expression is that of (68) for Wilson fermions and (71) for naïve fermions. Accordingly, we must tune the Raman detuning as a function of the experimental time by simply shifting the desired profile that is set by the scale factor In both cases, the time-dependence of the Ramanbeam detuning changes in time according to ∆δ(t) = 2ma(t − η 0 ), which depends on the bare mass m and is proportional to the shifted scale factor a(t − η 0 ). We recall that, in the numerical simulations, we have considered a particularly-smooth adiabatic switching (25), but other simpler profiles can also be explored in the laboratory. Finally, in the case of Wilson fermions, the detuning will also incorporate a static part δ 0 = 4t 1 , which depends on the recoil energy and the optical lattice depth (62).
In summary, the dynamics of a Dirac field in a de Sitter phase of expansion can be simulated through a specific time-dependent control of the Raman-beam detuning. In recent experiments with Fermi gases in two-dimensional Raman lattices [155], the consequences of changing the value of this detuning in the quench dynamics of the fermions has been explored, which can actually be used to infer the value of a topological invariant that captures the essence of an SPT groundstate, and even pinpoint the appearance of topological phase transitions. For the situation studied in our work, it is not sufficient to change the value of the detuning prior to the time evolution, but one rather needs to change it as time evolves. This allows to connect to the physics of quantum fields in expanding spacetimes, and could allow for a direct observation of the gravitational analogue of particle production, including the gravitational creation of topological modes in the spacetime boundary.
It should be noted that we have considered that the Hubbard interactions U (63) vanish via optical Feshbach resonances, since we have focused on a free LFT in a cosmological background. However, the same experimental scheme opens the path for the study of different interacting QFTs in a variety of background metrics. For instance, one of the most direct extensions of our model would be the Gross-Neveu model [51], which only differs from our original Hamiltonian by a four-Fermi interaction term that connects to the Hubbard interactions in the case of two-component Dirac spinors. This interaction term can be switched on by increasing the value of the scattering length, and would allow us to study non-perturbative phenomena such as chiral symmetry breaking or dynamical mass generation, and their interplay with particle production.

Measurement and fermion production
The remaining ingredient for the quantum simulation of fermion production in an expanding spacetime is to discuss how to measure the key observables (i.e. the spectrum of produced particles) in real time. To do so, one can take advantage of the various detection methods in ultra-cold atoms [164], such as the so-called time-of-flight (TOF) measurements, and the bandmapping technique. In TOF measurements [165], one abruptly removes all the applied fields that trap the atoms, letting the gas to expand freely. After this sudden turn off, if the atoms expand ballistically, there is a relation between their initial momentum and their final position ℏk = M x/t. Accordingly, after a certain time, absorptive imaging is used to obtain the spatial distribution n(x) of the atoms, which gives information about the momentum distribution prior to the release. In this absorptive imaging, photons from an incoming resonant laser are absorbed by the atoms, which consequently cast a shadow that can be recorded by a CCD-camera, giving one access to the so-called columnar -integrated-density. This technique can also be done in a spin-resolved manner by using laser beams with different frequencies, addressing each of the internal states. On the other hand, the band-mapping technique [166][167][168] turns off the external fields adiabatically, such that the band structure of the many-body system is slowly transformed into a free-particle dispersion relation. During this ramping-down, the quasi-momentum is approximately conserved, and Bloch states on the nth band are mapped onto free states with linear momentum on the nth Brillouin zone, giving thus direct access to the population of the different bands prior to the ramping down of the external fields and to its quasimomentum distribution.
Following the discussion in the previous sections, the quantum simulator can provide a gravitationalanalogue for particle creation by considering an initial half-filling condition that first populates the lower band of the lattice models. Then, a period of de Sitter expansion between two asymptotic flat vacua is simulated by the real-time evolution of the Fermi gas loaded in the Raman lattice with a time-dependent Raman detuning (73). During the expansion, not only the state of the Fermi gas will change, but also the Hamiltonian itself, and consequently the band structure of the system. Hence, if a fermionic atom is initially in a specific Bloch state, after the expansion, it will generally be left in a superposition of Bloch states corresponding to the two different bands, which is the analogue of Eq. (17) provided that a particle-hole transformation is applied. The number of produced particles for a certain quasi-momentum is then given by the probability for an atom with a certain quasimomentum to get excited to the higher band. The |β k (t + η 0 )| 2 coefficient can be thus obtained from this probability. The details of this result can be found in Appendix F. Since we can measure the population of each band for each quasi-momentum with the bandmapping technique [169], this Bogoliubov coefficient can thus be accessed experimentally. Going further, combining the band-mapping technique with spinresolved measurements [170], one could even characterize the band topology of the expanding lattice field theory [171].

Conclusions and outlook
In this work, we have developed the theory of particle production for a fermionic Dirac field in a (1+1)dimensional FRW spacetime, both from the perspective of the usual continuous theory and for two types of discretizations on the lattice. We have shown that the phenomenon of fermion production for a de Sitter expansion admits an exact solution in terms of a pair of decoupled Bessel differential equations. To avoid problems with the interpretation of the instantaneous vacua, we added an adiabatic switching that connects the de Sitter expansion to a pair of asymptotic Minkowski spacetimes, both of which have a well-defined notion of vacuum. We have numerically shown that the extra switching periods do not modify the particle production, which can still be described in terms of the analytically-solvable Bessel equation.
To pave the way for a quantum simulation of this phenomenon, we have considered two possible lattice discretizations of the Hamiltonian field theory associated to Dirac fermions in a curved spacetime: a naïve-and a Wilson-fermion discretization, which allow us to discuss universes with spatial boundaries. Focusing on the bulk of the lattice, we have shown that the naïve-fermion discretization reproduces the spectrum of particle production predicted by the continuum QFT for small momenta. However, as one approaches the edge of the Brillouin zone, the phenomenon of fermion doubling leads to a mirror-image of the spectrum, resulting in a doubled number of the produced fermions with respect to the continuum prediction. Turning into the Wilson-fermion discretization, which sends the fermion doubler to the ultraviolet cutoff of the QFT, we have shown that the spectrum and total number of particles produced at the bulk matches nicely the continuum expressions.
On the other hand, the asymptotic vacua of the Wilson-fermion lattice field theories can actually correspond to a couple of SPT groundstates characterized by a non-zero topological invariant, which has a boundary correspondence in terms of the groundstate degeneracy by the appearance of topological zeroenergy modes localized to the spatial boundaries of the FRW universe. This has allowed us to study the role of topology in the phenomenon of particle production and, in particular, to show that there can be gravitational particle production of zero-energy fermions localized to the edges of the spacetime. A similar phenomenon would also occur in situations where the mass of the Dirac fermion has some solitonic profile that changes its sign, leading to particle production in the form of domain-wall fermions. By numerically solving a set of coupled differential equations for the case of open boundary conditions, we have calculated the spectrum of fermion production, which is now a function of the energy of the produced particles as momentum is no longer a good quantum number in the absence of translational symmetry. We have shown that this spectrum allows one to identify clearly the production of fermions at the boundary of the expanding FRW universe, as one finds production of particles for energies below the mass gap.
To conclude, we have shown that current experiments of ultra-cold alkaline-earth Fermi gases in Raman optical lattices would be an ideal platform for the quantum simulation of this curved quantum field theories, provided that one exploits the simplifications that arise when working with conformal time. In connection to those experiments, our model requires to simulate the expanding spacetime by encoding its effect in a time-dependent mass that depends on the scale factor of the expansion, which corresponds in the experiment to changing the Raman-beam detuning as a function of time. This would allow for an experimental simulation of fermion production in a FRW spacetime, including its interplay with topology and the boundaries of such an effective universe.

Quantum field theories of Dirac fermions in a (d + 1)dimensional curved spacetime
In this Appendix, we review the formulation of Dirac QFTs in a curved spacetime of D = d + 1 dimensions. One typically starts from the flat-spacetime limit, in which the events are described by D-vectors x = (t, x) in a Minkowski spacetime of metric g µν (x) → η µν = diag(−1, +1, · · · , +1), where µ, ν ∈ {0, 1, · · · , d} label the spacetime coordinates. Here, we have chosen the mostly-plus metric [172], which is customarily used in treatments of general relativity [9]. The underlying Poincaré symmetry, which consists of spacetime translations and Lorentz transformations, has a specific representation for Dirac fermions that requires introducing spinor fields ψ(x) and the so-called gamma matrices γ µ , which obey a Clifford algebra {γ µ , γ ν } = 2η µν . Note that the choice of the mostlyplus metric exchanges the Hermitian (anti-Hermitian) nature of the temporal (spatial) gamma matrices with respect to the mostly-minus metricγ µ , which is the typical choice in particle physics [1]. Both matrices can be related by a simple global phase γ µ = iγ µ . In this manuscript, we stick to the mostly-plus metric in which the action of massive Dirac fields, which corresponds to the simplest field bi-linear that is a scalar under the Poincaré group [173], reads as follows Here, ψ(x) = iψ † (x)γ 0 is the adjoint field, which again differs with the standard choice for the mostly-minus metric [1], ∂ µ = ∂/∂x µ , and m is the bare mass. In this section, we use natural units ℏ = c = 1 and the repeated-index summation. When considering curved spacetimes, the metric is not necessarily flat, and one needs to exchange η µν → g µν (x) in the action (74). In addition, the integral measure must now include the volume form of the Lorentzian manifold associated to the curved spacetime. Therefore, one must substitute d D x → d D x √ −g, where g = det(g µν (x)), such that different regions of spacetime are weighted in a way that is invariant under diffeomorphisms. Finally, the partial derivatives in Eq. (74), which connect fields defined on nearby spacetime points of the flat Minkowski spacetime, must also be generalized in the presence of curvature, which leads to the covariant derivative Here, Ω ab parametrizes the specific Lorentz transformation, where we note that the changes with respect to Ref. [1] are due to the choice of the mostly-plus metric, and the consequent change in the metric and gamma matrices. To connect nearby spinor fields, the covariant derivative also contains a correction due to the spin connection ω ab µ (x), which leads to where we have introduced a connection field The specific form of the spin connection is found by introducing D-beins [174,175], which form a basis of vector fields with components e µ a (x) that allows us to express the curved metric in terms of the flat Minkowski one [12,54,56], namely g µν (x) = η ab e a µ (x)e b ν (x).
On the one hand, D-beins allow one to generalise the Clifford algebra to the curved spacetime γ µ →γ µ (x) {γ µ (x),γ ν (x)} = 2g µν (x), where the new set of curved gamma matrices is On the other hand, D-beins also play a key role in the spin connection. Making use of the Christoffel symbols Γ ν τ µ (x) = 1 2 g νσ (x)(∂ µ g στ (x) + ∂ τ g σµ (x) − ∂ σ g τ µ (x)) (81) for the curved spacetime, the spin connection becomes Here, the overall minus sign with respect to Ref. [175] comes from the different signature of the metric in the mostly-plus and -minus conventions. Note that, as customary in the context of curved Dirac fields, latin indexes a, b are raised (lowered) via the flat metric η ab (η ab ), whereas greek ones µ, ν require g µν (x)(g µν (x)).
Equipped with all these tools, Dirac fermions in curved spacetimes can be finally described by the action B Dynamical gravity in (1 + 1) dimensions In this Appendix, we review two formalisms that allow for a dynamical description of gravity in (1 + 1) dimensions, where Einstein's equations are not dynamical. There exist different alternative formulations for a consistent theory of gravity in D = 1+1 dimensions, such as the so-called Jackiw-Teitelboim model [176,177], which reproduces several phenomena characteristic of the higher-dimensional Einstein gravity. This model can also incorporate source terms that induce spacetime curvature, provided that one adds a term proportional to the trace of the stress-energy tensor in the classical field equations [61,178]. The solution of these classical field equations determines the corresponding metric g µν (x), and can lead to certain analogues of Einstein's gravity, including gravitational collapse, black-hole physics and, more relevant for the topic of this work, analogues of matter-and radiation-dominated FRW spacetimes [61]. We recall that, in this manuscript, we are concerned with a semi-classical approach that neglects back-action and treats the metric and the spin connection as classical background fields. At this level, we can use the solutions of JT gravity, and calculate their effect on the properties of the low-dimensional Dirac field. We thus avoid the need of treating the Dirac fields in conjunction with the auxiliary dilaton fields that allow for a covariant action of JT gravity [52,177].
Another possibility is to consider that the lowdimensional QFT arises as an effective field theory in situations in which the fermions are forced to propagate only along a two-dimensional section of an underlying curved spacetime in D = 3 + 1 dimensions [179,180]. This can be achieved, for instance, by considering an anisotropic mass of the fermion that is very large along the two remaining spatial directions. Alternatively, one may consider solitonic profiles for the mass, which connect to the aforementioned domain-wall constructions [46,48,49,79], and can also be used to impose constraints on the propagation of the fermions. From this perspective, one can use the specific time-dependence of the scale factor in a vacuum-dominated FRW spacetime in four dimensions, which are determined by the Friedmann equations [9], and incorporate it in the (1 + 1)-dimensional QFT of Dirac fermions in a curved metric. Note, however, that the number of spinor components is four in D = 3 + 1, whereas one works with two-component spinors in D = 1+1. Luckily, when considering a twodimensional section of the FRW spacetime (1), the required gamma matrices and the resulting connection field respect a block structure, such that the 4 spinor components decouple in two disconnected pairs that evolve in time independently. Accordingly, one can use a QFT with two-component Dirac spinors in the FRW spacetime with a time-dependence of the scale factor determined by the higher-dimensional Einstein equations. Both approaches lead to the same de Sitter expansion. We remark that this does not occur for generic metrics and stress-energy tensors.
The explicit expression of |ψ S ⟩ in terms of the Bogoliubov parameters can be obtained by noting that the fermionic operators in (117) realize an su(2) algebra, and so we can use the disentangled form of the two-mode squeezing operator [185] U k (ζ k ) = exp e iθ k tan(r k )â † k (η 0 )b † −k (η 0 ) · exp − ln(cos(r k )) n a k (η 0 ) +n b −k (η 0 ) − 1 · exp e −iθ k tan(r k )â k (η 0 )b −k (η 0 ) , wheren a k (η 0 ) =â † k (η 0 )â k (η 0 ) andn b −k (η 0 ) = b † −k (η 0 )b −k (η 0 ). Inserting this expresion in (124) and recalling that the Schrödinger and Heisenberg creation/annihilation operators, as well as the vacuum states, coincide at η = η 0 , prior to any dynamical evolution, one arrives at where the interpretation of particle creation is manifest: particles and antiparticles are created in pairs with oposite momentum as a consequence of the dynamics induced by the expansion of the universe. Note that if the adiabatic theorem holds, β k (η f ) would be negligible and thus the evolution would keep the state of the system in the instantaneous groundstate, |ψ S ⟩ ∼ |0⟩. This picture is particularly useful for the analysis of the analogue experiment with ultra-cold atoms. The experimental set-up described in Section 4 comprises a positive-and a negative-energy band which are symmetric around zero energy. The groundstate of this system is obtained by filling the lower band with atoms, and leaving the upper band empty. This situation, upon a particle-hole transformation, represents the vacuum state |0⟩ of the Dirac QFT. On the other hand, a situation where all the atoms are excited within the upper band represents a state with maximal content of pairs of particle/antiparticle, Thus, to simulate the phenomenon of particle production, the system must be initially prepared at halffilling in the groundstate, with all atoms within the lower band, so the initial state is |0⟩. Then, the dynamics induced by the expansion of the universe are simulated by means of the optical potential, as described in Section 4 and, consequently, at the end of the expansion atoms will no longer be in a well-defined Bloch state, but rather delocalized in both bands according to (126). Conversely, quasi-momentum is conserved. Since we can access to the energy band of each atom by means of the band-mapping technique, the coefficients |α k (η f )| 2 and |β k (η f )| 2 will be given by the probability of finding an atom with quasi-momentum k in a certain Bloch state. In particular, for a certain quasi-momentum, the probability of finding it in the lower (↓) or in the upper (↑) band will be given by Thus, the relation between the Bogoliubov parameters and the probability of excitation of the atoms is made explicit within the Schrödinger picture.