Uniqueness Criteria for the Fock Quantization of Dirac Fields and Applications in Hybrid Loop Quantum Cosmology

In generic curved spacetimes, the unavailability of a natural choice of vacuum state introduces a serious ambiguity in the Fock quantization of fields. In this review, we study the case of fermions described by a Dirac field in several cosmological spacetimes, and present recent results about well-motivated criteria that ensure the uniqueness in the selection of a vacuum up to unitary transformations. These criteria are based on two requirements. First, the invariance of the vacuum under the symmetries of the Dirac equations in the considered spacetime. Second, the unitary implementability of the Heisenberg dynamics of the annihilation and creation operators when the curved spacetime is treated as a fixed background. This last requirement not only permits the uniqueness of the Fock quantization but, remarkably, it also determines an essentially unique splitting between the phase space variables assigned to the background and the fermionic annihilation and creation variables. We first consider Dirac fields in cosmological spacetimes of 2+1 dimensions and then discuss the more relevant case of 3+1 dimensions. We use this analysis to investigate the hybrid loop quantization of a flat Friedmann-Lema\^itre-Robertson-Walker background cosmology coupled to a perturbative Dirac field. Among the Fock quantizations for the fermionic perturbations allowed by our criteria, we further restrict the choice of vacuum by demanding a finite fermionic backreaction and, moreover, by diagonalizing the fermionic contribution to the total Hamiltonian in the asymptotic limit of large wave numbers of the Dirac modes. Finally, we argue in support of the uniquess of the vacuum state selected by the extension of this diagonalization condition beyond the ultraviolet regime, proving that it picks out the standard Poincar\'e and Bunch-Davies vacua for fixed flat and de Sitter background spacetimes, respectively.


Introduction
Quantum Field Theory (QFT), namely the description of fields according to quantum rules, is one of the pillars of Modern Physics (see e.g. Refs. [1,2]). In this description, it is common to use a Fock formalism in which the physical processes are formulated in terms of creation and annihilation of field excitations around a vacuum state. Typically, these excitations are interpreted as particles (or antiparticles) of the field. This kind of description has been adopted successfully for all fundamental physical interactions, except for gravity: a fully satisfactory quantum field formulation of General Relativity, or more generally of the gravitational interaction, still remains as a defiant challenge.
In standard QFT, the particle excitations develop in flat spacetime. The symmetries of this background spacetime (regarded as classical) are employed in the selection of a natural state for the Fock representation of the field: Poincaré invariance fixes a unique vacuum in Minkowski spacetime [3]. However, it is well known that the generalization of the Fock quantization techniques to field theories in curved spacetimes is by no means straightforward, and in fact involves ambiguities. The most important of these ambiguities concerns the choice of a (quantum) representation of the algebra, obtained under Poisson brackets, of (the complex exponentiation of) the basic variables that describe the field, which are usually chosen to be canonical pairs. In traditional Non-Relativistic Quantum Mechanics, this algebra is known by the name of Weyl algebra (see e.g. Ref. [4]). Fortunately, in such case where the system has a finite number of degrees of freedom, the representation of the algebra in a Hilbert space is unique up to unitary equivalence, provided that certain mild conditions are imposed (including continuity). This uniqueness result is known as the Stone-von Neumann theorem [5]. This means that two different representations of the same Weyl algebra in the same Hilbert space are necessarily related by a unitary operator. This property ensures the robustness of the theoretical predictions of Quantum Mechanics, and in particular of those derived from the evolution of quantum states, something essential for the viability of the probabilistic interpretation of the quantum theory.
In QFT this uniqueness result for the quantum representation of the variables that describe the fields is generally no longer valid. It does not even apply to Fock representations of free field theories, which are typically governed by linear dynamical equations. In fact, it is well known (see e.g. [6]) that, for each given vacuum, there is an infinite number of linear canonical transformations, each of which provides a Fock representation of the field, that have no unitary correspondence in the Fock space. In short, this means that what in principle ought to be the corresponding unitary transformations just map the vacuum of the given representation to some vector which does not belong to its Hilbert space. Therefore, inequivalent quantum representations do exist in QFT, even in the simplest cases. This issue has been widely studied to be exclusive of systems, such as fields, with an infinite number of degrees of freedom, since a known mathematical criterion for the unitary implementability of a given linear canonical transformation involves a summability condition [7,8], which is trivially satisfied if the number of degrees of freedom is finite. This fundamental obstacle translates into the existence of an infinity of quantum descriptions of the same physical system that, in general, are not equivalent. However, one may introduce physically motivated requirements in order to reduce this ambiguity of the quantum description, and even to eliminate it completely, in certain situations. A remarkable example is given by fields that propagate in stationary spacetimes. This contains the case of fields in flat, Minkowski spacetime that we mentioned above. In these types of scenarios, the possible representations of the analogue of the Weyl algebra [commonly known as the field canonical commutation relations (CCRs) for bosons, or canonical anticommutation relations (CARs) for fermionic fields] are restricted to only a single representation if one imposes that the quantum theory incorporate the symmetry displayed by the background under timelike translations and that the evolution be generated by a positive Hamiltonian, that plays the role of an energy [6,9]. At the quantum level, this implies that the vacuum state of the field is stationary. There is therefore a natural unitary implementation of the dynamics.
The situation is notably more complicated when one considers fields that propagate in non-stationary spacetimes. Such systems describe scenarios of great physical interest, such as processes of star collapse or the cosmological evolution of the Universe (essentially since its very beginning). In these cases there is no timelike symmetry of the field equations that one can try and impose in order to restrict the admissible quantizations. The situation gets even worse if one takes into account that the quantum representations of the CCRs or CARs at different times are not necessarily unitarily equivalent. As a consequence, predictions based on the quantum evolution of the states lose robustness. For this reason, it is especially relevant to determine some physical criterion that allows us to remove the ambiguity in the choice of a Fock quantization in non-stationary spacetimes (or at least in a convenient subset of them) and, at the same time, regain a notion of unitary quantum dynamics for the fields.
Actually, this question has been investigated in recent years and the results indicate that the resolutions of the two problems are closely related. Indeed, for scalar fields in a multitude of non-stationary spacetimes of cosmological nature, it has been shown that a requirement of unitarity on the dynamics of the basic field operators in the Heisenberg picture (henceforth referred to as Heisenberg dynamics or evolution) can be used to guarantee the uniqueness of the Fock representation of the CCRs (up to unitary equivalence) if, in addition, one imposes invariance under the symmetries of the field equations [10,[16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34]. These symmetries include the spatial isometries. In many cases, these isometries suffice to reach the desired uniqueness when combined with the demand of a unitary Heisenberg evolution. Recent discussions on the topic of unitary dynamics in QFT in curved spacetimes from the canonical perspective can be found in Ref. [10] (see also Ref. [11] for related investigations on this issue). On the other hand, the nature of the vacuum state for fields in curved spacetimes and the Fock quantization of such fields from a covariant perspective have been widely investigated over the last decades, see e.g. Refs. [12][13][14][15].
In this review, instead, we focus our attention on fermionic fields. Many of the most abundant elementary particles in standard matter are fermions, and in this respect one can say that fermionic fields describe more realistic matter contents than scalar fields. In addition, although there exist well-founded results about the selection of Fock representations of fermionic fields in cosmology, the literature on this topic is not as prolific as in the case of scalar fields. In this sense, a review of the recent results obtained by us and our collaborators about uniqueness criteria for the Fock representation of fermionic fields, based on a unitary Heisenberg dynamics or other related properties, appears especially useful. For the sake of concreteness, most of our discussion is devoted to the particular case of a Dirac fermion field, to which we will henceforth refer simply as Dirac field.
The issue of determining a Heisenberg dynamics that can be realized as a unitary quantum transformation in fact involves a freedom in the splitting of the time dependence of the (fermionic) field. This time dependence can be separated in two parts: one that can be assigned to the quantum evolution of the creation and annihilation operators of the Fock representation and another that is due to the evolution of the background in which the field propagates. This second part can be treated as an explicit time dependence, via the background, when this is considered as a classical entity. Strictly speaking, this part of the evolution is not contained in the Heisenberg dynamics of the fermionic degrees of freedom. A fundamental idea in the search for a criterion to select a Fock quantization by imposing a notion of unitary dynamics is that the freedom in the splitting of the time dependence of the field can be employed to restrict the quantization in such a way that one ends up with a single family of equivalent representations while keeping nontrivial information about the fermionic evolution.
The idea of using the aforementioned freedom to arrive at a preferred class of Fock representations has proven to be very fruitful in frameworks that surpass the scheme of QFT in a curved classical spacetime. This is the case of fields with a dynamics that can be viewed as a propagation in an auxiliary background, or even quantum geometries that present regimes in which they can be treated effectively. For instance, this idea has been applied in the framework of hybrid loop quantum cosmology (hLQC), in which the spacetime is no longer a classical entity, but a quantum object [28][29][30]33]. For cosmological systems of notable physical interest, hLQC combines a loop quantization of the zero modes that (classically) describe the Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime that would correspond to a cosmological universe, with a Fock quantization of the field degrees of freedom that propagate in such a cosmological spacetime, typically viewed as perturbations [35,36]. The action of the gravitational system and its matter content is truncated at second order in these perturbations (generally assuming compact spatial sections). The combination of the loop and Fock techniques has to give rise to a consistent quantization of the total system, composed by the background cosmology and the perturbations. This consistency involves the imposition à la Dirac of a global Hamiltonian constraint, that interrelates the two types of representations as well as the evolution of the homogeneous universe and of its perturbations. For most of the relevant dynamical aspects of these perturbations, the information about the quantum geometry can be encapsulated in an effective geometry, that in many respects can be considered as emerging from a mean field approximation of the geometric degrees of freedom contained in the cosmological background. In this approximation, the fields that correspond to the (gauge invariant part of the) perturbations admit a QFT description where the background is the aforementioned effective geometry [35,36]. In this context, it is even clearer that building a formalism capable to maintain the unitarity of the Heisenberg dynamics in non-stationary spacetimes transcends the need to guarantee the robustness of the physical predictions of the theory. In hLQC, in particular, the additional advantage to pick out a unique family of equivalent Fock representations of the gauge invariant perturbations is that one can construct a Heisenberg evolution (with respect to some parameter of the complete quantum system) that behaves as a unitary quantum transformation in regimes where an effective background emerges.
Despite the attention that scalar (and tensor) perturbations deserve in cosmology, it is clear that a realistic matter content must include other types of fields, such as those that describe fermions, as we have already pointed out. Actually, in hLQC, recent works have introduced Dirac fermions and treated them as part of the perturbations [37]. The interest of contemplating the presence of these fields in the very early Universe goes beyond a formal question about the completeness of the description, because it is necessary to confirm that these fermionic fields do not affect substantially the otherwise well-established evolution of the primordial scalar perturbations, nor of the tensor ones. The results obtained in Refs. [37,38] support the expectation that the possible effects are ignorable.
It is worth remarking that there exists an inherent freedom to choose the splitting between the (fermionic) field variables and the degrees of freedom that describe the background, allowing one to change between different families of annihilation and creation variables. This can be done by means of transformations that mix all these degrees of freedom while preserving the canonical symplectic structure of the combined system, including the field canonical (anti-)commutation relations. If the fields are treated as perturbations, it suffices that the canonical symplectic structure is preserved at the level of the perturbative truncation adopted in the system. Instead of considering this freedom a nuisance, one can try and exploit it in order to define variables for which the Hamiltonian of the selected fermionic degrees of freedom has certain nice quantum properties, desirable from the viewpoint of a good physical and mathematical behavior [38].
In fact, when fermions were studied for the first time within the hybrid approach to loop quantum cosmology in Ref. [37], considering them as perturbations around a homogeneous and isotropic cosmological spacetime, the selection of fermionic variables for the corresponding Dirac field was restricted only by the requirements of invariance of the resulting Fock vacuum under the spatial isometries, a unitarily implementable Heisenberg evolution in the regime of QFT in a curved background, and a standard convention for particles and antiparticles. Nonetheless, with a rather reasonable choice made among the family restricted by these conditions, it was realized that the resulting Schrödinger equation for the fermionic degrees of freedom (after a sort of mean field approximation) involved ultraviolet divergences. To solve these divergences, one either has to appeal to a regularization scheme with subtraction of infinities or, alternatively, employ the remaining freedom in the choice of fermionic variables and restrict it even further by introducing additional requirements. Specifically, in Ref. [38] it was required that the fermionic backreaction be finite. In practice, this new restriction lowers the asymptotic order of the interaction part of the fermionic Hamiltonian at large wave numbers (defined for the Dirac field as the eigenvalues of the Dirac operator on the spatial sections of the background, in absolute value). As a consequence, the production of pairs of particles and antiparticles decreases for large wave numbers, becoming negligible asymptotically.
Actually, it is possible to go one step beyond in the same direction and, by taking advantage again of the freedom to split the degrees of freedom, demonstrate that one can absorb all the interaction terms of the fermionic Hamiltonian so as to make them identically zero order by order in the asymptotic regime of large wave numbers [39]. In this way, one clearly improves the quantum behavior of the (fermionic) field contribution to the Hamiltonian of the gravitational system. Moreover, one also greatly reduces the surviving ambiguity in the choice of Fock representation and vacuum for the field. Besides, since the resulting (fermionic) field Hamiltonian contribution is diagonal by construction on states with a definite number of particle (and antiparticle) excitations, at least asymptotically, the dynamics ruled by it is very simple. The vacuum of the naturally associated Fock representation changes only by a rotating phase. In this sense, one can interpret that this vacuum and the corresponding splitting of degrees of freedom in the hybrid quantization approach are those that get best adapted to the cosmological evolution. Finally, it is remarkable that this criterion of asymptotic diagonalization reproduces the standard choices of vacuum state in well-understood situations, within the scheme of QFT in a curved classical background [39,40]. This happens e.g. in the case of a flat spacetime, as well as for a de Sitter cosmology, scenario where the Bunch-Davies state is a natural vacuum [41].
The rest of this review is organized as follows. In Section 2 we provide the basics for the construction of a Fock representation of the CARs for a Dirac field, within the framework of QFT in a curved spacetime. We also summarize the results about the use of symmetries and of the unitarity of the Heisenberg evolution as criteria to select a preferred family of equivalent Fock representations. We then explicitly apply these criteria in Section 3, that deals with the case of a Dirac field in a non-stationary spacetime in 2+1 dimensions. The more interesting case of Dirac fields in a cosmological spacetime in 3+1 dimensions is reviewed in Section 4. After reviewing these aspects and results of QFT in curved spacetimes, in Section 5 we consider the use of canonical transformations to introduce a suitable splitting between the degrees of freedom of the cosmological background and of the Dirac field, treated in principle as a perturbation. In that section, we show how to use this freedom in the splitting to improve the quantum properties of the fermionic contribution to the Hamiltonian of the system, and in particular to make finite the backreaction that appears in it. The possibility of further employing this freedom to diagonalize the fermionic contribution to the Hamiltonian in the asymptotic limit of large wave numbers is reviewed in Section 6. There, we also explain that this diagonalization requirement can pick out a unique vacuum state under reasonable conditions, and that this state coincides with the natural one in situations of interest in QFT, like for Minkowski and de Sitter backgrounds. In addition, we also comment on the relation of adiabatic states with the vacuum selected by our criterion. Finally, we present the conclusions and some additional remarks in Section 7. We set the speed of light in vacuo, the Newton gravitational constant, and the reduced Planck constant equal to the unit.

Fock quantization of the Dirac field
This section contains some background material about the Fock quantization of a Dirac field in a curved spacetime. Special emphasis is put on the inherent ambiguity in the representation of the CARs associated with the infinitely many inequivalent complex structures available to construct the quantum theory, as well as on the combined criteria of symmetry invariance and of unitary implementability of the dynamics, that have been successfully employed to remove this ambiguity (and, even more, the ambiguity in the choice of basic field variables) in diverse, physically interesting fermionic (and bosonic) systems.
For the sake of clarity, let us begin our discussion by introducing the classical setting.

Background spacetime and Dirac equation
As backgrounds for the propagation of the field, we consider globally hyperbolic spacetimes, or just globally hyperbolic regions, (M ≈ I × Σ, g µν ), in either three or four dimensions. Here, I denotes an interval of the real line, and Σ is a Riemannian, Cauchy (hyper-)surface of dimension d, with d = 2, 3. For mathematical convenience, we restrict this surface to be topologically compact. In order to ensure that the spacetime (or region) admits a spin structure [42], we additionally require that M admit a global orthonormal frame [43]. So, the spacetime metric can be globally written as where η ab is the Minkowski metric in (d + 1)-dimensions, with signature {−, +, ..., +}, and e a µ is a(n orthonormal) coframe field, with dual frame e µ a . Throughout this work, Greek indices from the middle of the alphabet denote spacetime indices (µ, ν, ... = 0, ..., d), whereas Latin indices from the beginning of the alphabet account for the internal Lorentz gauge introduced by the frame (a, b, ... = 0, ..., d).
By employing an Arnowitt-Deser-Misner (ADM) decomposition of the considered spacetime (region) [44], we introduce a coordinate system in M, say {x µ } = {x 0 , x α }, with x 0 = t ∈ I being the time parameter and {x α } coordinatizing Σ. The line element in coordinates {x µ } reads where N and N α are, respectively, the lapse function and the shift vector, and h αβ is the induced metric on the Cauchy surface Σ. Let then Ψ be a free, complex, and anticommuting Dirac spinor with mass m, propagating in (M, g µν ). The dynamics of Ψ is governed by the first-order linear equation Here, the operator ∇ S µ stands for the spin lifting of the Levi-Cività covariant derivative [42], and γ a are the constant Dirac matrices that generate the Clifford algebra of a flat spacetime in (d + 1) dimensions: where I is the identity matrix and η ab is the (inverse of the) Minkowski metric. On account of the global hyperbolicity of (M, g µν ), the Dirac equation (3) has a well-posed Cauchy formulation [45]. So, given any smooth initial value Ψ( x) of the spinor field on a certain (compact) Cauchy surface Σ 0 , say at t = t 0 (where t 0 ∈ I is a fixed, but arbitrary, reference time), there exists a unique smooth solution Ψ(t, x) to Equation (3) which is defined on all of M and such that Ψ(t, x)| Σ 0 = Ψ( x). The solution Ψ(t, x), restricted to the domain of dependence of an arbitrary closed subset S of Σ 0 , depends only upon Ψ( x)| S . Henceforth, we fix Σ 0 (i.e., the section at t = t 0 ) as the Cauchy reference surface. Let S be the complex linear space of (smooth) solutions to the Dirac equation (3) which arises from the complex vector space P = {Ψ( x)} of (smooth) initial conditions at time t 0 . Note that, by construction, the map S Ψ → I t 0 (Ψ) = Ψ| Σ 0 is an isomorphism between the linear spaces S and P.

of 48
The space of Cauchy data P is naturally equipped with the product [45] (Ψ 1 , where γ a = η ab γ b , h is the determinant of h αβ , the dagger denotes the Hermitian adjoint, and n µ are the spacetime components of the future-directed unit normal to the Cauchy surface Σ 0 . The space of solutions S is endowed with an inner product of the form (5), though now at an arbitrary Cauchy surface Σ t ; this is so because of the independence of the mapping ( · , · ) D : S × S → C upon the spatial section on which it is evaluated. The exact meaning of what we understand by classical dynamical evolution in the space P is as follows. Given any solution Ψ(t, x) in S, the evaluation Ψ(t, x)| Σ t at a fixed time t defines the spinor Ψ ( x) = Ψ(t , x) on Σ 0 . Namely, we can take the induced spinor field on Σ t as initial condition on Σ 0 . Clearly, the mapping I t : S → P defined by is an isomorphism. Then, by considering the entire interval I, we get a one-parameter family of isomorphisms I t : S → P. The action of this family of mappings on a solution Ψ ∈ S gives the dynamical orbit of the Cauchy initial datum Ψ(t, x)| Σ 0 = Ψ( x) in P; that is, time evolution in the space of Cauchy data P is given by the one-parameter family of linear transformations T (t,t 0 ) = I t • I −1 t 0 .

Fock quantization, unitarity, and uniqueness
Let us next discuss the Fock quantization of Dirac fields using an approach based on the space of Cauchy data. It is worth remarking that an analogous approach constructed from the space of solutions is readily available, given the isomorphism I t 0 between both spaces.
We start by equipping the space of Cauchy data P with a complex structure, namely a real linear automorphism J : P → P with the property J 2 = −I, and such that it leaves the inner product (5) invariant. The complex structure allows for a natural splitting of P into two mutually complementary orthogonal subspaces (with respect to the considered inner product) P ± J = (P ∓ iJP )/2, that are eigenspaces of J with eigenvalue ±i. LetP be the complex conjugate of P, endowed with the complex where H J = H p J ⊕ H ap J is the one-particle Hilbert space associated with the complex structure J, and ⊗ a n H J denotes the n-fold antisymmetric tensor product of H J , with ⊗ a 0 H J = C [3,46]. Let {ψ p n ( x)} and {ψ ap n ( x)} be complete orthonormal bases for, respectively, H p J and H ap J . Then, the quantum field is (formally) represented in F J bŷ whereâ n is the annihilation operator associated with the spinorψ p n , whereasb † n is the creation operator associated with the spinor ψ ap n . The adjoint operatorsâ † n andb n correspond to the creation and annihilation operators of, respectively, particles and antiparticles. For a detailed discussion about the definition of fermionic annihilation and creation operators, see for instance Ref. [46]. For now, let us stress that the annihilation and creation operators here displayed correspond to the mode expansion projections of the (smeared) annihilation and creation operators specified in Ref. [46]. The basic operators {â n ,â † n ,b n ,b † n } satisfy the anticommutation relations, with the remaining anticommutators being null. The vacuum state of the theory corresponds to the unique (up to a phase) normalized cyclic vector in F J which vanishes under the action of all the annihilation operators,â n andb n . Let us emphasize that the choice of a complex structure for the Fock quantization determines the annihilation and creation operators of the theory. Thus, in general, different complex structures define different representations of the CARs. Furthermore, there exist infinitely many of these representations that fail to be unitarily equivalent [3]. This is where the ambiguity in the Fock representation of the Dirac field resides. Let us be more specific. Let F J and F J be two distinct Fock spaces, constructed from the different complex structures J and J . It can then be shown that, on the Fock space F J , the annihilation and creation operators defined by J (and which are naturally associated with the basis of F J ) are given by expressions of the form Here That is, the annihilation and creation operators defined by the two distinct complex structures J and J are related by a Bogoliubov transformation. By definition, unitary equivalence between the representations defined by J and J means that there exists a unitary operatorÛ : F J → F J intertwining the two representations, i.e. such that a n =Û −1â nÛ andb n =Û −1b nÛ . The transformation defined in Equations (9) and (10) is unitarily implementable then, in the sense that the transformation defined by the coefficients α f nm , β f nm , α g nm , and β g nm is a bona fide canonical transformation between the classical annihilation and creation variables corresponding to the considered operators.
A well-known result [8] states that unitary equivalence is achieved if and only if In general, given any two arbitrary complex structures, this condition is not satisfied. In fact, infinitely many inequivalent Fock representations of the CARs are possible, just as it happens with bosonic fields and their corresponding CCRs. The usual strategy to remove these type of ambiguities and to arrive at a (hopefully) unique Fock representation is to exploit the symmetries of the system. One typically requires that the complex structure (or the vacuum, in more physical terms) be invariant under some natural existing symmetries. As already mentioned in the Introduction, a crucial role is played here by time-translation invariance, and therefore that strict strategy fails to produce a unique representation in non-stationary scenarios, including very familiar and cosmologically relevant spacetimes.
Notice that a complex structure that remains invariant under time evolution immediately gives rise to a unitary implementation of (the canonical transformations generated by) the dynamics, therefore allowing the standard probabilistic interpretation of the quantum theory. It is therefore natural that, in non-stationary settings, one should try to preserve the unitary implementation of dynamical transformations, though giving up on (non-available) fully time-translation invariant complex structures, taking into account that this invariance is a sufficient, but by no means necessary, condition for such a unitary implementation. Let us make this more explicit. Suppose we are given a complex structure J on P, and construct the corresponding Fock representation, with associated operatorsâ n andb n as above. Since the field equations are linear, the transformations that correspond to time evolution from initial time t 0 to arbitrary time t are linear, and we obtain, for each t, new operators of the general form withâ † n (t) andb n (t) being supplied by the adjoint expressions of, respectively, Equations (13) and (14). It should be clear that, since the field evolution is a canonical transformation, the new operators satisfy the CARs, and we therefore have a family of new Fock representations. In fact, the new operators are simply those associated with the transformed complex structures T (t,t 0 ) JT −1 (t,t 0 ) [10], where T (t,t 0 ) is the evolution map introduced at the end of Section 2.1. Then, a unitary implementation of the dynamical transformations implies the unitary equivalence between all the new representations, for all t, and the original one defined by J. Of course, for a complex structure J that remains invariant under time evolution, the fulfillment of the unitary equivalence condition (12) is trivial, since all the nondiagonal beta coefficients of the associated Bogoliubov transformations are null.
On the other hand, there seems to be no compelling reason to relax the requirement of invariance under other natural remaining symmetries, such as isometries of the spatial manifold Σ, since invariant complex structures under these type of symmetries typically exist. So, the strategy that we adopt to deal with the ambiguity of the Fock quantization in non-stationary settings is the following. We require that the complex structure be invariant under the spatial isometries (and possibly other remaining symmetries of the system) and that it allows a unitary implementability of the dynamics. These combined criteria have been shown to be viable and effective in addressing the issue of the uniqueness of the quantization for a large class of field systems. The criteria were introduced for the first time in the context of midisuperspace models 1 , concretely to specify a unique preferred quantization of the inhomogeneous fields in Gowdy cosmological models [17][18][19][20][21]49], and since then they have been profusely and successfully employed to address the uniqueness of the quantization of (test) scalar fields in various, physically relevant cosmological backgrounds [22][23][24][25][26][27]29,30,32,34,50] (for a review, see Ref. [51]). Concerning fermionic fields and CARs, the same criteria have been successfully applied to single out a unique preferred quantum description for (test) Dirac fields in 2 + 1 dimensions [52] and in FLRW spacetimes [53][54][55][56], as we discuss in the next two sections.
It is worth pointing out that, as discussed in the Introduction, to achieve unitarily implementable dynamics in the type of non-stationary scenarios here considered, it is inevitable to explore the freedom in the splitting of the time dependence of the field between a genuine quantum Heisenberg evolution and an explicit dependence on the spacetime background. Typically, superimposed on the intrinsic dynamical evolution of the field variables, there is an explicitly time-dependent part coming from the non-stationary background spacetime itself. This last contribution to the total time dependence may be viewed as classical in nature, and effectively obstructs the possibility of a unitary quantum evolution. The solution is to extract the latter part by means of a time-dependent canonical transformation (performed at the classical level). This type of modification of the quantum notion of the field evolution is unavoidable in all of the cosmological systems analyzed so far, in order to 1 The term midisuperspace, originally introduced by K. Kuchař [47,48], refers to models that, even though possessing some symmetry, retain an infinite number of degrees of freedom, e.g. certain inhomogeneous models. recover a unitary implementability of the dynamics. Crucial in this approach is to pinpoint exactly the correct splitting between the intrinsic time dependence of the field and the time dependence coming from external factors, such as a non-stationary background. It is of the utmost importance to stress that this splitting is far from being arbitrary. It is guided, and to a great extent determined, by the requirement of a unitarily implementable dynamics.

Dirac fields in 2+1 dimensions
This section is devoted to discuss the applicability of the criteria of symmetry invariance and of unitary implementability of the dynamics in the Fock quantization of a concrete class of field systems, namely the case of a free Dirac field in 2 + 1-dimensional spacetimes which are conformally ultrastatic, with a time-dependent conformal factor. We show that, under rather non-stringent conditions on the time dependence of the cosmological background, and once a convention on the notions of particle and antiparticle has been established, a unique family of equivalent Fock representations is singled out by imposing (i) invariance under the unitary transformations that implement the symmetries of the equations of motion, and (ii) a nontrivial and unitarily implementable dynamics [52].

Dirac spinor in conformally ultrastatic spacetimes
Let us consider a fermionic field coupled to a globally hyperbolic, smooth manifold (or region) M, with the topology of I × Σ, where (as before) I ⊆ R is an interval of the real line, and Σ is a connected, compact, and orientable two-dimensional Riemannian manifold. Since, in particular, M is an orientable three-dimensional manifold, it is stably parallelizable [57]. We consider here conformally ultrastatic background geometries, so that the metric can be written as Up to the scale factor a(η), which contains the non-stationary information of the metric, 0 h αβ is the metric induced on the spatial surfaces Σ η defined at each fixed value of the conformal time η. The Dirac field couples to the geometry by means of the global (co)frame (1) defined, up to SO(2, 1) (orthochronous) gauge transformations, by the metric (15). Since in three dimensions any of the two irreducible complex representations of the Clifford algebra (4) are generated by 2 × 2 Dirac matrices, complex fermionic fields are locally represented by two-component spinors Ψ. In turn, we describe the components of these spinors by Grassmann variables, in order to encode the anticommuting nature of the fermionic field. We represent the Dirac matrices by The action for a fermionic field of mass m is given by where g is the determinant of the spacetime metric g µν and h.c. stands for Hermitian conjugate. Besides, using the spin connection one-form ω ab µ = 1 2 e νa ∂ µ e b ν + e νa e λb ∂ λ g µν − e νb ∂ µ e a ν − e νb e λa ∂ λ g µν , the spin lifting of the Levi-Cività covariant derivative is locally defined on the spinors as [42] It is convenient to partially fix the internal Lorentz gauge by choosing n µ e a µ = δ a 0 , where n µ is the future-directed unit vector field normal to the spatial surfaces Σ η . This leads to a reduction of the structure group of the bundle of oriented frames from SO(2, 1) (orthochronous) to SO(2) [58], restricting the spin structure to the double cover of the reduced frame bundle. For each of the considered spacetime manifolds and choice of spin structure on them, this restriction is well defined and provides an identical spin structure on each of all the two-dimensional leaves that foliate the spacetime manifold [42]. Thus, for each value of the time parameter, the field behaves as a spinor geometrically defined on each of the two-dimensional spatial manifolds Σ η that foliate the background. Equivalently, within the Cauchy data approach, the field can be described by one-parameter families of spinors, parametrized by the conformal time η, defined on an initial Cauchy reference surface Σ 0 , specified by η = η 0 . Let us denote by P the space of Cauchy data, namely the space of spinors at Σ 0 . The aforementioned one-parameter families on P are nothing but the result of evolving the Cauchy data in time (see Section 2). Thanks to the adopted gauge fixing, we can make a direct use of the spectral analysis of the Dirac operator defined on the two-dimensional Cauchy surface Σ 0 , instead of the Dirac operator on the whole Lorentzian geometry. In fact, in this gauge, the covariant derivative on spinors becomes simply where (2) ∇ S α is the spin covariant derivative on the spatial leaf with metric 0 h αβ . It can then be checked that where the prime stands for the derivative with respect to the conformal time η, and / D denotes the Dirac operator on Σ 0 .
Once the partial gauge fixing is performed, the inner product (5) on the space of Cauchy data P simplifies to where a 0 = a(η 0 ) and 0 h is the determinant of 0 h αβ . The Dirac operator / D is essentially self-adjoint with respect to the inner product (5) and, since Σ 0 is compact, it necessarily has a discrete spectrum, with eigenvalues ±ω n , labeled by natural numbers n, with ω n (≥ 0) growing with n [42]. Then, the space of Cauchy data P can be endowed with a basis formed by a set of eigenspinors of the Dirac operator. Let a −1 0 ρ np ( x) be the eigenspinors with positive eigenvalue ω n and orthonormal with respect to the inner product (5), where the index p accounts for the degeneracy. Since / D anticommutes with γ 1 γ 2 = γ 0 , we can chose as eigenspinors with negative eigenvalue −ω n those defined asσ np ( x) = −γ 0 ρ np ( x). Like ρ np ( x), these eigenspinors form an orthonormal set when the product is rescaled by the factor a −2 0 . With this rescaling, the set {ρ np ( x),σ np ( x)} provides a complete, orthonormal basis for P.
Let g n be the degeneracy of the eigenspace labeled by n, so that p = 1, ..., g n . The explicit form of g n depends on the spectral details of the Dirac operator / D and, consequently, on the particular 2-manifold considered. Nevertheless, for our purposes, we do not need the actual value of g n , but only to know its behavior in the ultraviolet regime of large eigenvalues ω n . So, let us introduce the counting function χ / D (ω) of the Dirac operator on a d-dimensional compact Riemannian manifold; that is, χ / D (ω) is the function that counts the number of positive eigenvalues of / D that are not greater than ω (including degeneracy). From the Weyl asymptotic formula [59], it follows that χ / D (ω) grows at most as ω d when ω goes to infinity. Using this result with d = 2, we conclude that the degeneracy behaves in the large n limit as g n = o(ω 2 n ), where the symbol o(ω 2 n ) means negligible with respect to ω 2 n . In terms of the basis {ρ np ( x),σ np ( x)}, the dynamical families of spinors in P (each one of them parametrized by the conformal time η) are given by Here, we use overlined symbols to indicate complex conjugation. Apart from the global factor a −1 (η), the time dependence (equivalently, the η-parameterization) is captured by the time-dependent coefficients s np andr np of ψ, which take care of the Grassmannian nature of the fermionic field. The auxiliary field ψ shows symmetric canonical Dirac brackets with its corresponding adjoint field that do not depend on the background [60,61]. We represent the algebra generated by these brackets in a Fock space, with the brackets being replaced with anticommutators , thus obtaining a Fock representation of the CARs and therefore a quantization of both the auxiliary and the original field, ψ and Ψ. By writing the field anticommutation relations in terms of the modes s np ,r np , and their complex conjugates, one finds that the only nonvanishing Dirac brackets are {s np ,s np } = −i and {r np ,r np } = −i, which are symmetric due to the anticommutativity of our Grassmann variables. In the quantum theory, they become anticommutators of the corresponding operators [60].
From Equations (17) and (24), the equations of motion for the fermionic modes are given by [52] s np = i(ω n + ima)r np , and their complex conjugates. These equations only couple the modes s np andr np (respectivelys np and r np ) with the same labels n and p, and do not depend on the degeneracy label p. They can be combined into the second-order differential equation where z np denotes either s np or r np . The general solution to this equation does not depend on the label p, except through the initial conditions, and is a linear combination of two complex independent solutions that we write in the form exp [(−1) l+1 iΘ l n (η)] with l = 1, 2. Let Θ l n (η 0 ) = Θ l n,0 and (Θ l n ) (η 0 ) = Θ l n,1 be the initial conditions at the initial reference time η 0 , and let us call Ω l n,0 = exp [(−1) l+1 iΘ l n,0 ]. A simple inspection shows that the integration constants of the general solution relate the initial conditions on Θ l n and their derivatives to the initial conditions s 0 np and r 0 np for the modes (and their complex conjugates), via Equations (25). One can then deduce that time evolution in the complex linear space of spinors ψ is dictated by the linear transformation [52] s np r np where the subindex η in column-vectors denotes evaluation at the given value of the conformal time, and the constants ∆ l n and ζ l n are wherel is the complementary of l in {1, 2}, namely {l,l} = {1, 2}.
In order to analyze whether the quantum theory admits a unitarily implementable dynamics, we do not really need to obtain the solution for V n (η, η 0 ). It is sufficient to know its behavior in the ultraviolet regime of large ω n . Under the mild condition that the scale factor be twice differentiable and with a second derivative that is integrable over each compact subinterval of the time domain I, a careful asymptotic analysis of the dynamics of the modes {z np } = {s np , r np } shows that two independent solutions to Equation (26) can be specified as follows [52]: for l = 1, 2, where ∆η = η − η 0 , and Λ l n (η) is a function with Λ l n (η 0 ) = 0 that, in the ultraviolet regime, is at most of order ω −1 n . With this choice, the constants (29) turn out to be

Fock quantization and unitary evolution
Let us now discuss the quantization of our fermionic system. Concretely, in this section we present the unique, preferred Fock quantization singled out by the criteria of symmetry invariance and of unitary implementability of the dynamics introduced in Section 2. The construction is performed in three steps. (1) We first focus on determining the family of invariant complex structures; namely those that commute with the action of the group of symmetries of the equations of motion for the modes. By construction, these complex structures lead to Fock vacuum states that are invariant under the unitary transformations generated by the symmetry group. We then consider time-dependent families of annihilation and creation variables associated with the invariant complex structures. By interpreting these families as dynamical trajectories, a specific redistribution of the implicit and explicit time dependence of the field is made. The dynamics that we wish to implement quantum mechanically is that of the implicitly time-dependent part, corresponding to the evolution of the annihilation and creation variables. (2) Each of the families of annihilation and creation variables defines an invariant Fock representation and a specific quantum evolution in the corresponding Fock space. We impose the criterion of unitary implementability of the dynamics, together with the requirement that the evolution be not trivialized. This leads us to set (or better said, characterize) all invariant Fock quantizations with a nontrivial and unitarily implementable dynamics. (3) Finally, we show that all such Fock quantizations turn out to be, in fact, unitarily equivalent, up to conventions in the notions of particles and antiparticles.
On account of Equation (25), it is clear that the field equations are invariant under the set of transformations that interchange eigenmodes of the Dirac operator with the same value of ω n . Since the Dirac operator is built from the spatial metric 0 h αβ , these symmetries include the isometries (if any) of the Cauchy surface Σ 0 . It should be clear that linear transformations commuting with the action of the symmetry group of the Dirac equation are composed of 2 × 2 blocks which, at most, can mix the modes s np andr np with the same value of p. In addition, using all the available symmetries one can reason that the blocks are necessarily the same for all the modes corresponding to the same eigenvalue of the Dirac operator (in norm) [52]. So, in particular, invariant complex structures are completely determined by a series of 2 × 2 matrices labeled by n ∈ N.
Let us remark that, given a complex structure, their associated annihilation and creation variables, namely the classical counterparts of the corresponding operators in an expansion of the type (7), diagonalize its action. Since invariant complex structures can mix only modes s np and r np with the same labels, the corresponding annihilation and creation variables must be linear combinations of these modes. The annihilation variables of particles and antiparticles are denoted by a np and b np , respectively, while the creation variables are their complex conjugatesā np andb np . These variables must satisfy the usual Dirac brackets for annihilation and creation sets, {a np ,ā np } = {b np ,b np } = −i and {a np , b np } = 0, which lead to the CARs (8) in the quantum theory. Now, let us consider time-dependent families of annihilation and creation variables associated with invariant complex structures. Specifically, In order that we get the required Dirac brackets for a np , b np , and their complex conjugates, the time-dependent functions f n l and g n l (l = 1, 2) must satisfy Combining these conditions, we can write where G n is a certain phase. Thus, it suffices just one complex function and two (real) phases for each n to characterize one family of annihilation and creation variables of the form (32). If we interpret the variables in each of these families as conforming to dynamical trajectories, their evolution differs (from each other, in general, and) from that of the modes s np andr np that dictate the evolution of the auxiliary spinor field ψ, because of the explicit dependence on η of the matrices F n . By substituting the inverse of Equation (32) in Equation (24), we obtain the expression of Ψ in terms of the introduced variables a np (η) andb np (η). The result makes it clear that the dynamics of Ψ is determined by the evolution of the annihilation and creation variables, and by the explicitly time-dependent contributions coming from F n and the scale factor. It is only the implicitly time-dependent part, that is, the part corresponding to the evolution of the annihilation and creation variables, the one that we want to implement as a quantum Heisenberg evolution. In the following, we restrict our attention to quantizations that are not only invariant, but also admit a unitary implementability of this dynamics.
For any of the allowed families of variables given in Equation (32), the dynamical evolution is a Bogoliubov transformation relating those variables at two different times, say the arbitrary time η and the initial time η 0 . By using Equations (27) and (32), we get where B n (η, η 0 ) = F n (η)V n (η, η 0 )F −1 n (η 0 ). This Bogoliubov transformation introduces the family of evolved complex structures J η = B n (η, η 0 )J η 0 B −1 n (η, η 0 ) on the space of initial Cauchy data, where J η is the invariant complex structure associated with (i.e., has a diagonal action on) the annihilation and creation variables at time η. As we have seen in Section 2, the transformations (35) are implementable as unitary operators on the Fock space defined by J η 0 if and only if Let us recall that the number of eigenstates of the Dirac operator with positive eigenvalue not greater than ω grows as ω 2 in the ultraviolet regime. Using this asymptotic behavior, it is possible to show [52] that, for large N 1 and arbitrary N 2 > N 1 , and with where K is some positive constant. This identity is important for the subsequent analysis because it implies that the sequence { √ g n ω −2 n } n∈N is square summable. Employing Equations (30) and (31), as well as relations (33) and (34), one can check that, in the asymptotic regime of large ω n , where the integrals are over conformal time from η 0 to η, the symbol O stands for asymptotic order, and h can be set equal to either f or g. In order to lighten the notation, we have omitted the dependence of these functions on η and denoted the evaluation at η 0 with the superscript 0, preceded by a comma. It is worth noticing that relations (33) and (34) guarantee that |β f n | = |β g n |. So, for the purpose of unitarity, it suffices to analyze just one of these types of coefficients.
It is convenient to employ again the notation {l,l} = {1, 2}. Then, given that |h ñ l | 2 + |h n l | 2 = 1, we get h ñ l = e iH ñ l 1 − |h n l | 2 , where H ñ l denotes some phase that may depend on time. It is worth remarking that, in addition to conditions (36) which severely restrict the behavior of the coefficients h n l , both in their mode and time dependence, one should naturally demand that h n l be such that the dynamics of the annihilation and creation variables is not trivialized when compared with the original Dirac evolution. Otherwise, the criterion of a unitary implementation of the dynamics would be useless, since it would pose no restriction on the Fock representation. Indeed, one may always extract all of the asymptotically dominant time dependence of the fermionic field by means of explicitly time-dependent canonical transformations, and trivialize in this way the requirement of a unitarily implementable evolution. More specifically, by examining Equation (38), it can be seen that the dominant contribution of the Dirac dynamics dictated by V n (η, η 0 ) is given by imaginary exponentials of the phases ±ω n ∆η. Thus, in order to avoid a trivial evolution, we rule out the possibility that these dynamical contributions are counterbalanced with a specific choice of (time and mode-dependent) phases in the linear combinations that determine the annihilation and creation variables.
All together, taking h as either f or g, the requirements of a nontrivial and unitarily implementable dynamics impose, as a necessary condition, that asymptotically [52] h n l = for a subset of the natural numbers, n ∈ N ± l , and with ϑ n h,l being some mode-dependent and time-dependent complex function that goes to zero in the limit of large ω n . Here, the union of the four subsets N ± l gives (up to a finite number of elements) the natural numbers, allowing for the possibility that up to three of these subsets be empty, and having identified h n l with h n 1 for n ∈ N ± 1 and with h n 2 for n ∈ N ± 2 , with the ± superscripts indicating the relative sign for h ñ l in the second and third identitites of Equation (39). By substituting this equation into relation (38), as well as using that the integral of (Σ 1 n − Σ 2 n ) behaves as m (a − a 0 ) /ω n + O(ω −2 n ) in the ultraviolet regime [52], one gets that the asymptotic behavior of the complex norm of the beta coefficients, for all n ∈ N ± l , is up to certain terms that are negligible compared with the largest order between ω −1 n and ϑ n h,l . So far, no condition has been set on ϑ n h,l other than it must go to zero in the ultraviolet limit. However, by adding the requirement that the sequence {g n |β h n | 2 } n∈N ± l be summable, one gets a restriction on how fast ϑ n h,l must tend to zero. The line of reasoning is the following. One assumes that the beta coefficients are square summable in the subsets N ± l , including degeneracy, and then one looks for the functions ϑ n h,l (if any) that solve the resulting conditions. From Equation (40), one finds that there are two different situations that lead to distinct conditions on ϑ n h,l , namely, either the sequence { √ g n ω −1 n } n∈N ± l is square summable, or not. In the first situation, it is not difficult to check that is not a square summable sequence), by using the implications of Equation (37) and recalling that any trivialization of the fermionic dynamics is excluded, one can see that the functions [52] must form a sequence that is square summable, including degeneracy, in the subsets N ± l where { √ g n ω −1 n } n∈N ± l fails to satisfy such square summability. In total, given an invariant family of complex structures J η characterized by the annihilation and creation variables (32), the necessary and sufficient conditions for the corresponding Fock representations of the CARs to be unitarily equivalent, and therefore to support a unitarily implementable nontrivial dynamics, are the following. (1) The functions h n l and h ñ l must be asymptotically of the form (39) for n ∈ N ± l . (2) The terms ϑ n h,l must be such that either is square summable, they form a sequence that is square summable (including over the degeneracy), or otherwise, (2b) the sequence { √ g nθ n h,l } n∈N ± l is square summable, withθ n h,l given in Equation (41).
Up to now, we have characterized all families of annihilation and creation variables that (i) share the symmetries of the equations of motion, and (ii) evolve according to nontrivial dynamics that are unitarily implementable at the quantum level. Each of these families determines an invariant Fock representation (e.g. that associated with the choice of an invariant complex structure at the initial time η 0 ) with a nontrivial, unitary quantum evolution in the corresponding Hilbert space. By referring to the combination of a Fock representation and a specific quantum dynamics as a Fock quantization of the system, the question now is whether the invariant Fock quantizations with unitarily implementable nontrivial dynamics are equivalent quantum theories or not. Before we address this issue of uniqueness, let us make some remarks about the preceding results.
Note that the criterion of unitary implementability, together with the requirement of a nontrivial dynamics, fix (up to phases) the leading order behavior of the coefficients h n l in the asymptotic regime of large ω n , as it is shown in Equation (39). In addition, for all those cases where { √ g n ω −1 n } n∈N ± l fails to be a square summable sequence, the imaginary part of e −iH n l ϑ n h,l must have its dominant asymptotic contribution of order ω −1 n , and equal to the function ma/( √ 2ω n ) in absolute value. This follows simply by realizing that the square summability of the sequence (41). Hence, for a nonzero fermionic mass m, the coefficients h n l asymptotically depend in a very specific way on the eigenvalue of the Dirac operator and on the mass of the field. Most importantly, there is also a specific dependence on time, by means of a dependence on the background where the field propagates. Finally, let us emphasize that the analysis about unitarity holds not only for a positive value of the mass m, but also when this mass is zero. That is, the families of annihilation and creation variables for the massless Dirac field, selected by the criteria of symmetry invariance and of unitarity, are fully characterized by coefficients (h n l , h ñ l ) with an asymptotic form given by Equation (39), and such that the sequences { √ g n ϑ n h,l } n∈N ± l are square summable. It is worth noticing that, within this family, the choice f n 1 = f n 2 = g n 1 = −g n 2 = 1/ √ 2 provides a representation which corresponds to the Fock quantization constructed from the (celebrated) conformal vacuum, i.e., the natural vacuum specified by imposing the conformal symmetry of the massless Dirac equation in the quantum theory.

Uniqueness of the quantization
Let us now address the issue of uniqueness. With this aim, we proceed as follows. First, as a reference, we adopt a certain Fock quantization that is invariant and possesses a nontrivial and unitarily implementable dynamics. Next, we consider any other invariant Fock quantization that admits a nontrivial, unitarily implementable dynamics, and we examine whether it is unitarily related with the reference Fock quantization or not. If the answer is in the affirmative, then the uniqueness is proven.
A simple choice of reference quantization is the Fock quantization characterized by annihilation and creation variables with Note that, in the case of the massless field, this choice defines the natural quantization with conformal vacuum. Let {ā np ,b np ,ã np ,b np } be any other choice of annihilation and creation variables that defines a Fock quantization with a nontrivial, unitarily implementable dynamics. The coefficientsf n l andg n l that characterize these variables then satisfy all of the conditions stipulated in the previous subsection. Given our reference Fock quantization and this other arbitrary one allowed by our criteria, it follows from Equation (32) that the annihilation and creation variables associated with them are related via the time-dependent Bogoliubov transformation K n (η) =F n (η)F −1 n (η), so that It is not difficult to check that the norm of the off-diagonal coefficients is |λ h n | = |h n 1 h n 2 −h n 2 h n 1 |. Besides, using Equation (34), one gets that |λ f n | = |λ g n |. Thus, the square summability conditions on λ f n and λ g n , that must be satisfied in order for the Bogoliuvov transformation to be unitarily implementable, turn out to be just one (and the same) condition. Then, the transformation determined by the sequence of matrices K n is implementable on the reference Fock space as a unitary operator if and only if where |λ In case this condition is satisfied, we can consider the two analyzed quantizations as physically equivalent. Note that the above condition ensures that the Fock representations defined for every value of the conformal time η are unitarily equivalent.
Taking e.g. h = f in the formulas of the preceding subsection, we have that the coefficientsf n l andf ñ l are, respectively, of the asymptotic form (39). Then, one gets that, for n ∈ N + l [52], up to terms O(|ϑ n h,l | 2 ) in the asymptotic limit of large ω n if the field is massless, or up to terms O(ω −1 n ) if the mass does not vanish and { √ g n ω −1 n } n∈N + l happens to be square summable. Alternatively, if this sequence is not square summable and the field is massive, one must have for n ∈ N + l . Since the sequences formed by ϑ ñ f ,1 and ϑ ñ f ,2 in the two first cases above, and byθ ñ f ,1 andθ ñ f ,2 in the last case, are square summable by hypothesis in their respective subsets, including degeneracy, it follows that the unitary equivalence condition (43) is immediately satisfied for all n ∈ N + l .
On the other hand, it can be shown that Equation (43) fails to be satisfied if any of the considered subsets of integers N − l has infinite cardinality. The reason for this failure is rooted at the difference in the relative sign of the coefficients in the pair (f n 1 ,f n 2 ) with respect to that in ( f n 1 , f n 2 ) [52]. Therefore, if one then insists and interchanges the roles off n l andg n l in the definition of the annihilation and creation variables for the subsets N − l , something that in practice amounts to an interchange between the relative signs of the pair (f n 1 ,f n 2 ) and the signs for the pair (g n 1 ,g n 2 ) [see Equation (34)], one gets that both pairs would display the same relative signs as the coefficients of the reference quantization, and then condition (43) would be satisfied for n ∈ N − l . According to Equations (32), the exchange off n l andg n l can be interpreted as a change in the convention of what are particles and what are antiparticles. In this sense, the inequivalence between quantizations before one performs the explained interchange can be understood as a spurious result coming from the the fact that we are just considering two Fock quantizations with the opposite convention for the concept of particle and antiparticle in an infinite number of modes.
In summary, the criterion of invariance under the symmetries of the equations of motion and the unitary implementability of a nontrivial quantum dynamics removes the ambiguities in the representation of the CARs, both for the massive and for the massless Dirac fields in 2 + 1 dimensions, selecting a unique family of unitarily equivalent Fock representations, together with a notion of quantum evolution, up to conventions about the concept of particles and antiparticles.

Fock quantization of Dirac fields in FLRW cosmologies
We now discuss the Fock quantization of Dirac fields in cosmological spacetimes of the FLRW type. In particular, in this section we show that one can achieve results about the uniqueness of the quantization very much like in the previous section. More precisely, we consider minimally coupled massive Dirac fields, propagating in homogeneous and isotropic FLRW spacetimes, with 3-dimensional spatial hypersurfaces that can be either spherical or toroidal. The spherical case was notably analyzed in great depth by D'Eath and Halliwell [62], within the context of the Wheeler-DeWitt approach to quantum cosmology. In that seminal treatment, a special time-dependent family of Fock representations was chosen for the Dirac field, by means of an instantaneous diagonalization of the Dirac Hamiltonian. In particular, such family is associated with vacua which are invariant under the symmetries of the Dirac equation. Moreover, it was shown in Ref. [62] that particle production over time remains finite for those vacua, a fact that is an exclusive characteristic of quantizations that admit an unitary implementation of the dynamics. We now show that this family of Fock representations is in fact uniquely selected, up to unitary equivalence, by the criteria of invariance under spatial symmetries and unitary implementability of the dynamics. Also, we present a similar result concerning the Dirac field in flat (compact) FLRW spacetime. In this case, besides spatial translations, the symmetry group includes as well helicity generated spin rotations.

Dirac spinors in FLRW cosmologies
As before, we consider spacetime manifolds with topology of the type I × Σ, where I is a connected interval of the real line and Σ is certain spatial Cauchy surface. In the case of a spherical universe, Σ is isomorphic to the 3-sphere, S 3 , whereas for universes with (compact) flat spatial sections Σ is isomorphic to the 3-torus, T 3 . In the following, we often refer to the above two situations as the spherical case and the flat case, for S 3 and T 3 respectively. The metrics associated with such universes can be written as where 0 h αβ is either the 3-dimensional spherical metric or the flat metric, and a(η) is the scale factor. It follows from previous comments in Section 2 that spin structures can always be defined in the above types of cosmological spacetimes [42,43]. Dirac fields Ψ, of mass m, correspond therefore to sections of the associated spinor bundle. Explicitly, we adopt the Weyl representation of the Dirac matrices where σ 0 = σ 0 is the identity 2 × 2 matrix and σ i = − σ i (with i = 1, 2, 3) are the Pauli matrices. Such representation of the generators of the Clifford algebra allows us to describe Dirac fields by means of a pair of two-component spinors φ A andχ A possessing well-defined and opposite chirality. We take φ A to be the left-handed projection of Ψ, whileχ A is the right-handed one. Moreover, we adopt the same conventions as in Ref. [62] concerning the treatment of spinor indices. In these cosmological spacetimes, the action for the Dirac field of mass m is where the spin covariant derivative ∇ S µ is given by relations (18) and (19), adapted to the models in question.
Let us perform again a partial fixing of the internal Lorentz gauge, with the purpose of providing a rigorous treatment of the spatial dependence of spinors, as well as to analyze their properties under the symmetry groups associated with the considered FLRW cosmologies. To be specific, the gauge group of the orthonormal and oriented frame bundle can be reduced from SO(3, 1) to SO(3) [58]. In the considered spacetime models, this procedure gives rise to a well-defined restriction of the spin structures to the double cover of the reduced bundles, having SU(2) as gauge group. Finally, this restriction provides one (and the same) spin structure on each of the spatial hypersurfaces that constitute the foliation of the considered cosmologies. In the case of S 3 , the spin structure turns out to be unique [63,64]. In the flat case, on the other hand, there are eight possible spin structures associated with distinct periodic conditions for the spinors in T 3 [65,66]. In practice, this partial gauge fixing is obtained again by imposing the conditions n µ e a µ = δ a 0 which, moreover, greatly simplify the Hamiltonian analysis of the system [61]. In particular, one can check that the Dirac operator (defined over left and right-handed spinors) on the reference spatial hypersurface Σ 0 takes the form where e αAA is the spinor version of the triad and (3) D α is the spin lifting of the Levi-Cività covariant derivative with respect to the metric 0 h αβ [62].
Once the gauge fixing has been performed, we introduce the following inner product on the space of left-handed spinors defined on the spatial hypersurface Σ 0 (as well as the corresponding definition for right-handed spinors): where summation over repeated indices is assumed, I AA denotes the components of the identity matrix, and φ A and χ A are arbitrary spinors. Since both the spherical and the toroidal hypersurfaces are geodesically complete [67], it follows that the Dirac operator (49) is essentially self-adjoint with respect to the inner product (50), with discrete spectrum in both cases. The spectrum of the Dirac operator on S 3 consists of the following sequence of eigenvalues [62,64]: each of which has a corresponding degeneracy g n = (n + 1)(n + 2). Using a notation similar to that employed in the previous section, the left-handed eigenspinors with eigenvalues ω n and −ω n are denoted by ρ np A andσ np A , respectively, where the label p = 1, ..., g n accounts for the degeneracy. Once normalized, the set of all such elements forms an orthonormal basis -with respect to the inner product (50) -for the space of left-handed spinors in S 3 . An analogous basis for right-handed spinors is readily obtained by complex conjugation.
In the spherical case, the chiral projections of the Dirac field can then be written as Concerning now the flat FLRW model, the spectrum of the corresponding Dirac operator -with compactification period l 0 -consists of the following sequence of eigenvalues [65,66]: where the three v j 's form the standard orthonormal basis for the lattice Z 3 , and the three numbers j ∈ {0, 1} characterize each of the possible choices of spin structure 2 . Given any such spin structure, 2 These spin structures determine the periodicity or antiperiodicity of the Dirac field in each of the orthogonal directions that define T 3 . If, in harmony with spatial isotropy, one imposes the same global behavior for the field in all these directions, the choice of spin structure is restricted to either j = 0 or j = 1 for all j.
we identify the label k in ω k (or equivalently in −ω k ) with the norm of any of the wave vectors k ∈ Z 3 corresponding to ω k . The degeneracy g k associated with each eigenvalue ω k (or −ω k ) does not possess in this case a closed expression. However, well-known results in Riemannian geometry allow us to conclude that g k grows asymptotically as O(ω 2 k ), for unboundedly large ω k [29,68]. Let us then fix a spin structure on T 3 and choose a set of triads associated with the flat metric. The eigenspinors of the Dirac operator then form a basis of the space of spinors in T 3 , basis that can be made orthonormal with respect to the inner product (50). In particular, if one chooses a diagonal triad such that the spin connection 1-form becomes null, the Dirac operator turns out to be the standard one for flat Euclidean space. Its (left-handed) eigenspinors corresponding to eigenvalues ±ω k are where u are some x-independent two-component spinors, subject to the condition that the eigenvalue equation must hold. They can be chosen so that the spinors w are normalized in the inner product (50) and such that Σ 0 for all k, k = − j v j /2. Summation over repeated indices is assumed 3 . Finally, the constants C  (52) and (53), with corresponding Grassmann variables m k , r k , t k , and s k .
Returning to the spherical FLRW case, notice that upon introduction of the mode expansions (52) and (53) in the Dirac action, and once the second-class constraints of the fermionic system are solved [69], one ends up with the following symmetric Dirac brackets for the mode variables [60,61]: where the ordered pair (x np , y np ) stands for (m np , s np ) or (t np , r np ). Using Grassmann variational derivatives and requiring stationarity of the action, one obtains Dirac equations for the modes: as well as the complex conjugate versions. One can combine these dynamical equations to obtain decoupled second-order equations, that are actually the same for all modes {x np , y np } with the same label n. Denoting the mode variables generically by {z np }, the resulting equation for given ω n and (nonzero 4 ) mass m is the following: The general solution to this equation is a linear combination of two independent complex solutions, which we write again in the form exp[iΘ 1 n (η)] and exp[−iΘ 2 n (η)]. The general expression of the fermionic modes at arbitrary time η can then be written as a linear transformation of the 3 Except for the index k on the right-hand side of the second relation in Equation (56). 4 The case m = 0 is slightly different, although straightforward to handle.
corresponding initial values, that assigns a different weight to the two independent solutions of Equation (59): where ∆ l n and ζ l n , l = 1, 2, are constants that depend on the initial conditions of the independent solutions exp[iΘ 1 n (η)] and exp[−iΘ 2 n (η)], and of their derivatives, at the reference time η 0 (see Ref. [53] for explicit expressions and details).
Just as in the previous section, one can obtain the asymptotic behavior of the linear transformations (60) in the ultraviolet regime of large ω n , and subsequently study the unitary implementability of such transformations at the quantum level. For that matter, let us impose the initial conditions Θ l n (η 0 ) = 0. One obtains the following expressions for the (exponents of the) solutions to Equation (59) [53]: where a 0 = a(η 0 ), ∆η = η − η 0 , and the functions Λ l n are solutions of an equation of Riccati type which become negligible in the limit of unboundedly large ω n . In particular, the functions Λ l n have a behavior of the type O(ω −1 n ). The asymptotic values obtained for the constants ∆ l n and ζ l n are the following [53]: The analysis concerning the fermionic dynamics in the flat case is quite similar, leading to mode solutions that are the analogues of Equations (60): The corresponding asymptotic expressions, in the limit of large ω k , are and where the functions Λ l k are solutions of a Riccati equation with an asymptotic behavior of the type O(ω −1 k ).

Fock quantization and unitary evolution
Considering both the spherical and the flat FLRW cosmologies, we proceed now to characterize all Fock representations for the Dirac field which satisfy the following requirements. First, the associated vacua must be invariant under the action of the natural symmetries of the system, among them the spatial isometries (and helicity generated spin rotations, in the flat case). Secondly, the Fock quantizations are required to admit a (nontrivial) unitary implementation of the dynamics at the quantum level.
Let us start by analyzing the behavior of Dirac spinors in the spherical case, when the isometry transformations of S 3 are applied. The transformation group is then SO(4), or equivalently the double cover Spin(4) = SU(2) × SU(2), with action in S 3 defined by means of a Clifford multiplication. As a consequence, Spin(4) acts on the cross-sections of the spinor bundle on S 3 [64]. Notice that this action, when viewed on the four-component Dirac spinor, is reducible to two blocks: the action of Spin(4) over spinors φ A , and the complex conjugate action over spinorsχ A . Both such representations of Spin(4) are unitary with respect to the inner product (50), and it follows that each block is further decomposable in a direct sum of irreducible representations.
Invariant vacua are associated with invariant complex structures, and we therefore seek complex structures which commute with the action of Spin(4) over spinors. To begin with, given the decompositions (52) and (53), any complex structure can be seen as an infinite dimensional matrix in the basis formed by the modes {m np ,r np , t np ,s np }. From this point on, we follow the analysis carried out in Ref. [64], concerning the eigenspaces of the Dirac operator on S 3 . Consider the action of Spin(4) on spinors of the type φ A . Using Frobenius theorem [70], it is shown in Ref. [64] that each eigenspace corresponds exactly to a representation space of one of the irreducible representations of Spin(4) on spinors. Moreover, each such irreducible component shows up in the direct sum with multiplicity equal to one. This means that the representation spaces generated for each n by the sets simply reproduce those generated by the sets (67), and therefore each irreducible representation associated, for each n, with one of the sets (68) is unitarily equivalent to a corresponding one coming from the sets (67). All this information (combined with an inspection of the dynamical mode equations [54]) allows us to apply Schur's lemmas [71] to conclude that a complex structure that commutes with the action of the group of isometries of S 3 over the space of Dirac spinors cannot mix modes m np ,r np , t np , ands np with different values of n. Moreover, within the subspace associated with a fixed value of n, such complex structures cannot mix the modes {m np ,s np } with {t np ,r np }, since the two sets provide inequivalent irreducible representations of Spin(4). Let us consider first, for each given n, the subspace generated by the modes {m np ,s np }. The restriction of an invariant complex structure to any such subspace can then be characterized by means of four linear maps, relating in all possible pairings the two subspaces generated by {m np } and by {s np }. Each of these four maps has to be proportional to the identity, as ensured by Schur's lemma. Finally, similar considerations can be applied to the subspace generated by the modes {t np ,r np } [54]. We now turn to the flat FLRW model, and consider the action of isometries of T 3 on the Dirac field. Restricting our attention to continuous transformations, the isometry group is generated by constant translations along each of the orthogonal directions of the 3-torus. A general translation on the torus is thus x → x + θ, where for each component we have 2πθ α /l 0 ∈ S 1 . For any given choice of spin structure on T 3 , one can easily check that a general translation simply results into the following transformation: in each of the elements (55) of the basis of (left-handed) eigenspinors. For each k ∈ Z 3 there are therefore two copies of the same 1-dimensional complex irreducible representation, with different k's giving rise to inequivalent representations [55]. Again, one can perform the same analysis for the spinors of opposite chirality. Then, taking into account the mode decomposition of the Dirac field Ψ, and using again Schur's lemma, one concludes the following. A complex structure that commutes with the action of translations on T 3 can at most mix the modes (m k ,s − k−2 τ , t − k−2 τ ,r k ) among themselves, for each fixed k ∈ Z 3 , and is trivial otherwise. In the flat case there is an additional symmetry of the Dirac system, following from the conservation of helicity in the evolution of the Dirac field in conformal time η. In fact, one can consider the projection of the spin angular momentum in the direction of the linear momentum of the particle, a projection that (except for the subspace generated by the modes with ω k = 0) defines the helicity operator h [72]: Here, ∇ is the standard 3-dimensional Euclidean gradient and σ denotes a vector with components given by the Pauli matrices. Eigenspinors of h with eigenvalues +1 or −1 are said to have positive or negative helicity, respectively. With the choice of gauge such that the spin connection vanishes, it turns out that the matrix blocks of h (apart from the factor [− ∇ 2 ] −1/2 ) correspond exactly to the Dirac operator on T 3 . One can then check that the positive helicity part of the Dirac field Ψ is generated by the coefficients m k ands k , for all k ∈ Z 3 different from τ. On the other hand, the negative helicity contribution is generated by the modes t k andr k , for all k ∈ Z 3 different from τ. A simple inspection of the equations of motion (58) shows that helicity is indeed a conserved quantity. Therefore, one can include, as an additional symmetry of the fermionic system, the 1-parameter group of spin rotations generated by helicity, by means of the complex exponentiation of h/2 multiplied by the angle of rotation. Such group is immediately unitary with respect to the inner product (5), since the operator h is essentially self-adjoint. It follows that the unitary implementation of this symmetry at the quantum level is ensured whenever the complex structure that defines the quantization does not mix positive helicity modes with negative helicity ones. The complex structures characterized above define the sets of creation and annihilation operators that provide invariant Fock representations of the CARs for the Dirac field in the considered homogeneous and isotropic scenarios. In the spherical case, let us denote the classical counterparts of the annihilation operators for particles and antiparticles by a , respectively for particles and antiparticles. We recall that the pairs (x, y) (with the appropriate labels) denote any of the ordered pairs of mode coefficients (m, s) or (t, r). In the following, we consider all the possible (time-dependent) families of fermionic creation and annihilation variables selected by invariant complex structures.
In the case of S 3 , the creation and annihilation variables in question can then be written as Once more, the label η denotes dependence on conformal time. Notice that the time-dependent functions f n l and g n l (with l = 1, 2) may differ for the pairs of modes (m np ,s np ) and (t np ,r np ), although this is not explicit in the notation. The following relations must again be satisfied: such that anticommutators of the type (8) are obtained at the quantum level.
Turning to the flat case, the relations (71) are replaced with where conditions analogous to Equation (72) again apply. When evaluated at different times, the sets of variables (71) are related to each other by means of dynamical Bogoliubov transformations. The general form of such linear transformations can be obtained taking into account Equations (60) for the evolution of the fermionic modes in the spherical case. In general terms, the Bogoliubov transformation relating the annihilation and creation variables at the initial time η 0 with those at any time η is given by a sequence of blocks B n such that The absolute values of the coefficients β where the integrals are over conformal time from η 0 to η, h denotes either f or g, and the superscript 0 stands for evaluation at the initial time.
Analogous considerations, of course, apply to the flat case, with dynamical Bogoliubov transformations between variables of the type (73) that are characterized by matrices B k , with corresponding coefficients β f k and β g k , for which the explicit expressions can be found in Ref. [55]. Considering for instance the Fock representation defined by the annihilation and creation variables at the initial time η 0 (or equivalently by the associated complex structure), the dynamical Bogoliubov transformations (74) can be implemented on the corresponding Fock space by means of unitary quantum operators if and only if [7,8] Actually, it is sufficient to ensure the above condition for either h = f or h = g, since it follows again from relations (72) that |β g n (η, η 0 )| = |β f n (η, η 0 )| [53]. Let us then fix h to be either f or g. The unitary implementability of the dynamics therefore depends on the asymptotic behavior of the beta coefficients in the limit of large ω n , that in turn depends on the behavior of the sequences h n l . A detailed analysis, carried out in Ref. [54], shows that, apart from an uninteresting alternative which would effectively trivialize the quantum dynamics in a similar way as it was discussed in the previous section, the fulfillment of the unitarity condition (76) requires that the functions h n l behave asymptotically as where {l,l} is the set {1, 2}. Moreover, the sequences ϑ n h,l must be square summable, including degeneracy.
Before continuing with our discussion, a comment is in order. Since we have not put any restriction on the global asymptotic behavior of the sequences h n l for fixed l, it is possible that neither h n 1 nor h n 2 actually converges over N. In fact, the sum (76) can be made finite with h n 1 taking the role of h n l in Equation (77) for n in a subset of N and h n 2 taking that role over a complementary subset (modulo finite subsets of N). Hence, the above behavior of h n l is required only for n ∈ N l ⊂ N, with N 1 ∪ N 2 = N modulo finite subsets (including the possibility of one of the subsets N l being empty).
The analysis concerning the flat case follows similar lines, apart from a careful handling of the already mentioned issue of the accidental degeneracy of the Dirac eigenspaces (for full details, see Ref. [55]). The conclusion is that the requirement of (nontrivial) unitary implementation of the dynamics completely fixes the explicit dependence on time in the dominant part of the Dirac field, in a very similar way as in the spherical case [see Equation (77)], for an infinite number of modes in the asymptotic sector of large ω k . Hence, the time-dependent scaling that must be introduced in the mode dynamics of the Dirac field, such that the remaining part of the evolution can be unitarily implemented, is completely fixed (in absolute value) at dominant order. In particular, the scaling factor for each mode is similarly determined in terms of ma/(2ω k ), and also includes the global term a −3/2 , introduced in the mode expansions (52) and (53) adapted to the case of T 3 .

Uniqueness of the quantization
In the previous subsection we have characterized the Fock quantizations of the Dirac field which allow a unitary implementation of the dynamics and which possess vacua that are invariant under the natural symmetries of the considered cosmological models. We now show that, for each of these models, all such quantizations are unitarily equivalent. Direct consequences are the removal of quantization ambiguities typically present in QFT and the selection of a very specific, well-defined notion of quantum evolution.
Let us start with the case of spherical sections. One of the simplest choices of functions that satisfy the conditions derived above for a unitarily implementable quantum dynamics is for all n ∈ N and for both pairs (m np ,s np ) and (t np ,r np ). We take this choice as defining our reference family of complex structures. Let us then consider any other family of invariant complex structures allowing a unitary implementation of the dynamics. Such family is defined by certain annihilation and creation variables,ã (x,y) np andb (x,y) np , with coefficientsh n l (for h identified with f or g and for l = 1, 2) that have the asymptotic behavior described at the end of the preceding subsection. In particular, the subdominant sequences ϑ ñ h,l appearing in Equation (77) are square summable (degeneracy included) over subsequences N l . The relation between this family and the reference one is given by a Bogoliubov transformation determined by a sequence of matrices K n such that One can check that the absolute values of the nondiagonal elements of these matrices have the following expression [53,54]: where h stands again for both types of functions f and g. The two Fock quantizations, namely the reference one associated with the family (78) and the one defined by the coefficientsh n l , are unitarily equivalent if and only if the Bogoliubov transformation (79) is itself unitarily implementable, a demand that is equivalent to the fulfillment of the conditions ∑ n g n |λ f n (η)| 2 < ∞ and ∑ n g n |λ g n (η)| 2 < ∞, (81) for all η of interest. Once more, one of the above conditions is redundant, since the equality |λ g n | = |λ f n | is again ensured by relations (72).
Let us then focus on the first condition in Equation (81) and consider the situation such that the asymptotic behavior (77) applies to the functionsf n l . The alternative, with Equation (77) applying to the functionsg n l , can be treated in a completely analogous way. One can show from Equation (80) that, in the limit of unboundedly large ω n , the coefficients λ f n have the following behavior on the subsequence N 1 : Since by hypothesis ϑ ñ f ,1 is square summable (including the degeneracy g n ) over N 1 , it follows immediately that the condition for unitary equivalence is satisfied if the set N 2 is finite (or empty). Suppose now that N 2 is an infinite set. It is clear that λ f n behaves asymptotically like O(1) for n ∈ N 2 , and therefore the summability required in Equation (81) cannot be attained, leading to apparently inequivalent quantizations. However, this inequivalence stems, as it happened in the previous section, from the fact that the quantization associated with the coefficientsh n l defines a convention for the concept of particles and antiparticles which is completely the opposite, for an infinite number of modes, of the convention corresponding to the reference quantization. Once both conventions are reconciled, the quantizations are seen to be physically equivalent. In fact, suppose that, in our reference quantization (78), we switch the convention concerning particles and antiparticles for all modes corresponding to N 2 . This redefinition is effectively attained with the interchange f n l ↔ g n l (n ∈ N 2 ), as follows from the definition (71). Then, the behavior of the new coefficients λ f n would no longer be O(1) on N 2 , but they would behave (in norm) as |ϑ ñ f ,2 | + O(ω −2 n ) instead. Since, also by the hypothesis of a unitarily implementable evolution, ϑ ñ f ,2 is square summable (degeneracy included) over N 2 , one concludes that conditions (81) are now satisfied, confirming therefore the unitary equivalence between the two Fock quantizations, after the two conventions concerning particles and antiparticles have been harmonized.
The analysis of the uniqueness of the Fock quantization of the Dirac field in the flat FLRW case proceeds in a similar fashion. One can again choose a reference complex structure allowing a unitary implementation of the dynamics, namely the one characterized by the following matrix elements in Equation (73), for all ω k = 0: Applying the same type of arguments as above, one can prove [55] that, as it happened in the case of S 3 , once a convention concerning particles and antiparticles is fixed, the condition that there exist a nontrivial unitary implementation of the dynamics is sufficient to ensure the unitary equivalence of all Fock quantizations associated with invariant vacua. Let us conclude with a brief comment on a key difference between the current study and previous analogous works on quantum scalar fields in FLRW cosmologies [22,24,25,29,30,50]. Like in the scalar field case, our criteria uniquely fix not only the Fock representation once a suitable set of variables for the field has been chosen, but actually they greatly reduce as well the ambiguity in this choice of variables that arises from time-dependent linear redefinitions. Considering, for instance, a spherical spatial topology, this affects the global scaling introduced in the decompositions (52) and (53), as it does for the scalar field, but now it affects also the scalings in the particle and antiparticle contributions that are induced by the time dependence of the functions f n l and g n l subject to the unitarity conditions (77). As such, the Dirac field presents specific and different time-dependent scalings in its particle and antiparticle parts, that are also different for each of the two chiralities, introducing an aspect which is absent in the scalar field analysis. For the scalar field, the requirement of unitary dynamics imposes a global scaling of the original field variable, such that the scaled field in practice behaves like a conformally coupled field, and one might wrongly believe that the possibility of attaining a unitarily implementable dynamics is somehow constrained by the availability of a conformal symmetry in the scaled theory (at least in the ultraviolet regime). The work on Dirac fields definitely puts aside such type of misconceptions.

Hamiltonian backreaction of Dirac perturbations in hLQC
The previous uniqueness results about the Fock quantization of Dirac fields in FLRW cosmologies are of special importance beyond the context of QFT. Actually, they can be used (and further developed) to confer great robustness to the full quantization of a homogeneous and isotropic cosmological spacetime perturbed with small matter inhomogeneities described by a Dirac field, in the context of hybrid quantum cosmology. As mentioned in the Introduction, this framework for the quantum description of cosmological systems employs techniques from a theory of quantum cosmology (here we focus our attention on the case of loop quantum cosmology) for quantizing the spatially homogeneous zero modes of the geometry, while the inhomogeneous fields are quantized using standard Fock representations. In this setting, the features of the resulting quantum theory and its predictions are strongly affected by the precise knowledge obtained not only about a unique preferred Fock space for the fields, but also about which part of their degrees of freedom displays a genuine unitary quantum evolution when one reaches regimes where the cosmological background behaves classically. This knowledge serves in hLQC to separate in a specific way this homogeneous background from the variables that describe the fermionic perturbations, and such splitting can be refined by further imposing some physically sound properties on the full quantum system, such as a proper definition of the fermionic part of the Hamiltonian operator. Along these lines, in this section we are going to revisit the main results in the hybrid loop quantization of a flat homogeneous and isotropic cosmology with fermionic perturbations, in particular in what concerns the study of the Hamiltonian operator and the consequences on the quantum backreaction of the fermionic matter on the cosmological background.

Fermionic perturbations in flat FLRW: Splitting of the phase space
Let us start by considering the Einstein-Hilbert action restricted to symmetry-reduced universes with a metric given in Equation (46), taking 0 h αβ as the Euclidean metric in coordinates adapted to the spatial homogeneity, and letting the lapse function, that we call N 0 , be arbitrary. We particularize again the discussion to a topology of the flat spatial sections equal to the compact T 3 -topology. In order to include standard inflationary scenarios in our system, we minimally couple a homogeneous scalar field φ with potential V(φ), that plays the role of an inflaton. In the canonical ADM framework, the degrees of freedom of this cosmological model can be described with the scale factor a, the inflaton φ, and their canonical momenta, respectively denoted by π a and π φ . On classical FLRW solutions, these variables are subject only to one constraint, arising from the zero mode of the Hamiltonian constraint, that generates global time reparameterizations: where we recall that l 0 is the compactification length of T 3 .
To include inhomogeneous fermionic content in this cosmological model, we minimally couple a Dirac field Ψ of mass m and treat it entirely as a perturbation (including its purely homogeneous part, if it had any). For physical completeness, one may also introduce purely inhomogeneous and anisotropic perturbations of the metric and the inflaton field 5 . Within this perturbative hierarchy, we conduct the analysis at the lowest nontrivial order and thus we truncate the whole action of the system (and its associated symplectic structure) at quadratic order in all of the perturbations. Since the Dirac action is precisely quadratic in the fermionic field, at this order it only couples with the homogeneous sector of the cosmology 6 . Therefore, in practice, we can treat its associated spinor bundle as if it were defined on a pure (spatially flat) FLRW universe. Furthermore, this coupling implies that the Dirac field is a perturbative gauge invariant at our order of truncation, namely it is independent of perturbatively linear coordinate redefinitions that respect the manifest homogeneity of the background.
In order to take advantage of the previous results about the uniqueness of the Fock quantization of the Dirac field, and to allow for a convenient spatial mode expansion of the fermionic fields, we partially fix again the local Lorentz gauge imposing the condition n µ e a µ = δ a 0 on the tetrad of our background. This allows us to generically perform a decomposition of the two chiral components of the fermionic perturbations as that given in Equations (52) and (53) (after replacing the eigenspinors of the Dirac operator on S 3 by their analogues for T 3 ). Recall that these behave under local gauge transformations as SU(2) spinors defined in T 3 , with a spin structure that is given by the choice of the vector τ [c.f. Equation (54)]. Of equal importance is the fact that this partial gauge fixing eliminates all the nontrivial canonical brackets between the homogeneous FLRW geometry and the (rescaled) Dirac field [61]. Then, the only nonvanishing brackets between the fermionic mode coefficients m k ,r k , t k , ands k (and their complex conjugates) are where (x, y) again denotes any of the two possible ordered pairs of mode coefficients (m, s) or (t, r) (omitting their associated wave vector labels). The introduction of this fermionic mode expansion in the Dirac Hamiltonian coupled to our considered background gives rise to the following contribution to the total Hamiltonian of the system: As we have explained above, N 0 is the homogeneous lapse function of the FLRW background. It follows that the Dirac perturbations, at our considered order of quadratic truncation in the action, contribute only to the global zero mode of the Hamiltonian constraint of the entire system. Explicitly, if one ignores the rest of perturbative fields in our model (which do not couple to the fermionic field at quadratic order), the zero mode of the Hamiltonian constraint is given by the sum of H |0 and H D . At this point in the discussion, it is worth remarking that there exists an inherent freedom in the description of the cosmological model at hand. Treating the system formed by the FLRW universe 5 The zero modes of the metric and scalar field can be conveniently isolated and accounted for in the scale factor and homogeneous inflaton, owing to the compactness of the spatial sections. and its perturbations as a whole entity, one can always mix the different sectors of the phase space by means of canonical transformations. Even if this mixing does not affect the physical behavior of the system at the classical level at the end of the day, the freedom in identifying the sets of basic variables that describe each sector can strongly affect the properties of the hybrid quantization, given that the perturbative fields are quantized with a different type of representation (à la Fock) than the homogeneous background (with loop techniques). Focusing the attention on the choice of splitting between the FLRW sector and the fermionic one, one can understand this as the specific assignment of how each of these two types of degrees of freedom contributes to the dynamics of the entire system. To take into account this panorama, and with the aim put on characterizing the Fock representation for the fermionic perturbations, we introduce general families of annihilation and creation variables of the form (73) where the coefficients f k 1 , f k 2 , g k 1 , and g k 2 are now taken as functions of the canonical variables that describe the geometric sector of the cosmological background, understood hereafter as the sector that corresponds to the scale factor and its momentum 7 . Furthermore, these functions are subject to the conditions where G k and F k 2 are phases given by generic real functions of the FLRW geometry. These conditions are nothing but the equivalent of Equation (72) (adapted to the case of T 3 ). Therefore, they guarantee that the definition of the fermionic annihilation and creation variables is a canonical transformation in the fermionic sector of the phase space (namely, if one freezes the FLRW background variables and ignores their Poisson brackets).
Each of these new families of fermionic variables can codify in different ways the possible splittings in the dynamical behavior, between the homogeneous cosmological geometry and the fermionic perturbations, that preserve the linearity in the perturbative sector. However, when the cosmological system is viewed as a whole dynamical entity, and unless the coefficients f k 1 , f k 2 , g k 1 , and g k 2 are constant, the new fermionic variables do not form a canonical set with the scale factor and its canonical momentum. These variables must be modified if one wishes to restore the canonical algebra fulfilled by the original basic set for the description of the system. In other words, the canonical transformation that started with the above definition of families of fermionic annihilation and creation variables must be completed. This can be readily done, at our order of quadratic perturbative truncation, by demanding that the symplectic potential of the FLRW and fermionic sectors remain unchanged after the transformation, up to contributions that are of higher perturbative order than quadratic. This procedure leads to the new corrected variables for the scale factor and its momentum: Here, the partial derivatives affect only the explicit dependence of the fermionic modes on the cosmological background geometry, via the functions f k 1 , f k 2 , g k 1 , and g k 2 . Thus, the new FLRW variablesã,π a , φ, and π φ form a canonical set with the fermionic annihilation and creation variables defined by means of Equations (73) and (89). 7 One can generalize the analysis to coefficients that are functions also of the inflaton and its momentum, thus allowing for a dependence on all the degrees of freedom of the FLRW background. However, this generalization is not necessary for our discussion, except at some point in Section 6, where we comment on it explicitly.
In the following, for convenience, we restrict our attention to families of annihilation and creation quantization of the entire system, it is therefore of the utmost importance to adhere to physical criteria and restrict the allowed families of fermionic annihilation and creation variables. In view of the results presented in the previous sections, a first reasonable condition to impose is that, when a QFT regime in a classical background spacetime is recovered, the quantum Heisenberg evolution of the fermionic annihilation and creation operators can be implemented unitarily. As we have seen, once we set a convention for particles and antiparticles, this criterion on the considered families of variables completely fixes the (symmetry invariant) quantum representation of the field, up to unitary equivalence. In other words, it allows us to determine the Fock space for the representation of the fermionic field. For concreteness, let us set the particle/antiparticle convention to be such that it corresponds to the standard one in Minkowski spacetime as the mass of the fermions goes to zero (situation in which the rescaled Dirac field propagates as if it were in flat spacetime, in conformal time) [55]. Then, the families of annihilation and creation variables restricted by the criterion of a unitarily implementable dynamics are those such that, in the asymptotic limit of large ω k : The above condition greatly restricts the asymptotic behavior of the allowed families of annihilation and creation variables. Furthermore, all such choices lead to unitarily equivalent representations. However, there is still much freedom left, even in the asymptotic regime of large ω k , as the function ϑ k is only constrained to be square summable (including degeneracy). Actually, this freedom can be used to refine the phase space splitting in such a way that the part of the hybrid quantization that concerns the fermionic degrees of freedom displays several desirable properties. In particular, it appears physically sound to demand that the Fock quantization of the contributionH D to the Hamiltonian constraint result into a well-defined operator on the fermionic vacuum state (and thus on the dense Fock subspace of its associated n-particle states). Such an operator can be simply obtained by promoting the variables a , on the other hand, respectively to annihilation and creation operators in the fermionic Hamiltonian given in Equation (93). As these variables commute with the (perturbatively corrected) ones that describe the homogeneous FLRW background, the well-definiteness ofH D on the fermionic vacuum is insensitive to whether these FLRW variables are fixed as classical or promoted to quantum operators as well (within the hybrid loop scheme). In turn, if one imposes normal ordering on the products of annihilation and creation operators, it is fairly easy to see that this property on the Fock quantization ofH D depends exclusively on the asymptotic dependence on ω k of the terms h k I defined in Equation (96). These provide the interactive part of the fermionic Hamiltonian, that is responsible for the annihilation and creation of pairs of particles and antiparticles, and the resulting operator is well defined (with normal ordering) on the vacuum state if and only if they form a square summable sequence, including the degeneracy of the Dirac eigenvalues. Specifically, if one introduces condition (97), together with relation (89), in h k I , the dominant contribution to this function in the asymptotic regime of large ω k is cancelled out. From a Hamiltonian perspective, this cancelation is the ultimate responsible for the unitarity of the Heisenberg evolution (in the context of QFT in curved spacetimes). However, it does not guarantee that h k I forms a square summable sequence over all k. The necessary and sufficient condition for this to happen, and thus for the Fock quantization ofH D to be well defined on the vacuum, turns out to be that asymptotically [38] It is worth noticing that this last condition on the allowed families of fermionic annihilation and creation variables restricts even further the admissible phase space splittings between the fermionic sector and the background. Namely, it specifies even further how the assignment of the dynamical content of each sector should be made.
Finally, hereafter we restrict the discussion exclusively to phases G k such that {G k , H |0 } = 0, in order not to introduce any artificial asymmetry between the dynamical behavior of the fermionic variables that describe particles and antiparticles [see Equation (93)].

Hybrid quantization: Fermionic backreaction
Let us next summarize the main steps that must be followed for the hybrid quantization of our cosmological system, formed by an inflationary FLRW background coupled to a perturbative Dirac field. As mentioned above, one can freely include metric and inflaton perturbations as well, and discuss a similar phase space splitting and choice of Fock representation for them. We nonetheless ignore them in this review, as they do not necessarily affect the fermionic sector at the considered order of truncation in the action. For details on the treatment of metric and inflaton perturbations in the context of hLQC, we refer the reader to Refs. [35,36,73].
The first step towards the hybrid quantization of the full system is to specify a concrete splitting of phase space with physically good properties, along the aforementioned lines. In other words, one starts by identifying a specific set of fermionic annihilation and creation variables [by means of Equations (73) and (89)] within the family restricted by the asymptotic conditions (97) and (98). In addition to the phase space splitting, such choice fixes also the Fock representation for the fermionic degrees of freedom. Indeed, for their quantization, one just needs to promote these fermionic variables to the corresponding annihilation and creation operators, that in turn completely specify the Fock space F D from their associated cyclic vacuum state. On the other hand, once a preferred phase space splitting has been identified, the FLRW background spacetime is described by the set of variablesã, φ, and their canonical momenta [see Equations (90) and (91)], that, as mentioned above, Poisson commute with the chosen set of annihilation and creation variables, at the classical level. We then adopt a (discrete) loop quantum cosmology representation for the canonical variables that describe the homogeneous background geometry, with operators defined on a Hilbert space that we call H grav kin . For specific details on this representation, see Refs. [37,73]. In this review, we only recall its main features when they are relevant for the quantum dynamics of the fermionic sector. In addition, a standard Schrödinger representation is adopted for the homogeneous inflaton and its momentum, with Hilbert space given by L 2 (R, dφ), such that the inflaton acts by multiplication and the momentum is represented as −i∂ φ . The total representation space for the hybrid quantization of the full canonical set of basic variables of our cosmological system is then the tensor product of all the introduced individual spaces, namely H grav kin ⊗ L 2 (R, dφ) ⊗ F D . This tensor product space is often called the kinematical space, and it is not the fully physical one. Indeed, the whole system is subject to the zero mode of the Hamiltonian constraint, that can be found in Equation (92) and classically generates global time reparameterizations. We implement this symmetry at the quantum level by demanding that physical states be annihilated by the representation of the constraint as an operator on the kinematical Hilbert space 8 . Actually, for mathematical convenience, one often rather imposes a rescaled version of this constraint, obtained through multiplication by the volumeṼ =ã 3 l 3 0 of the homogeneous FLRW sector. In order to find physical states in the system, it is therefore necessary to specify the corresponding constraint operator and deal with the ambiguities that its representation involves, as it is not a linear function of the canonical variables. In this construction, in particular, we impose normal ordering in the fermionic contributionH D to the constraint. For details about the remaining ambiguities and how to reasonably fix them in the context of loop quantum cosmology, we refer the reader to Ref. [37].
Let a specific quantum representation of the zero mode of the Hamiltonian constraint (with normal ordering for the fermionic operators) be provided in this way. In order to search for states of physical interest, such that the influence of the perturbations on the FLRW background can be made controllably small, we look for solutions to the quantum constraint starting from the following ansatz for the allowed wavefunctions Ξ: Here,Ṽ denotes dependence of the wavefunction on the geometric sector of the homogeneous FLRW background, and N D is a generic label indicating the occupancy numbers in the fermionic Fock space. Hence, in our states we can regard Γ as the partial wavefunction that describes the behavior of the FLRW cosmology, whereas ψ D is the partial state for the fermionic perturbations. We notice that both of these contributions are allowed to depend on the inflaton. It is also worth noting that, in case inhomogeneous perturbations of the metric and inflaton were included, this ansatz can be correspondingly generalized, separating the dependence of the total wavefunction in each different type of perturbation. In what concerns the partial FLRW state Γ, we further impose as part of our ansatz that it evolves unitarily with respect to its variation in φ (so that we can normalize it in the Hilbert space H grav kin ), with a generator that we callĤ 0 such that We take this generator to be a positive self-adjoint operator and furthermore impose that it gives rise to partial states that are close to exact quantum solutions (i.e., annihilated by the action) of the (rescaled) constraint operator of the homogeneous FLRW background. Explicitly, let us write the constraint operator that we would have for our FLRW model in the absence of perturbations in the form where −Ĥ (2) 0 is the operator that represents the (rescaled) contribution of the inflaton potential and of the homogeneous FLRW geometry [see Equation (84)]. Our restriction onĤ 0 translates into demanding that the action of (Ĥ 0 ) 2 −Ĥ (2) 0 − i[∂ φ ,Ĥ 0 ] on Γ be small, let us say at most comparable to terms of quadratic order in the perturbative parameter of our system.
For the wave profiles selected by the above ansatz to potentially become physical states, we must impose that they be annihilated by the (rescaled) constraint operator. Namely, the action of the quantum representation of the function (92), rescaled with the homogeneous volumeṼ, on states of the form (99) (and satisfying the aforementioned conditions) must be zero. The resulting constraint equation can be greatly simplified, and viewed as an evolution equation on the partial fermionic wavefunction, if the following approximations are applied [37,73]: i) The partial state Γ is such that one can ignore transitions in the FLRW geometry mediated by the zero mode of the Hamiltonian constraint. Then, one can apply a kind of mean-field approximation and take the inner product of the constraint equation with Γ, with respect to the Hilbert space H grav kin . ii) The contribution ∂ 2 φ ψ D can be neglected when compared with Ĥ 0 Γ ∂ φ ψ D . In other words, the contribution to the inflaton momentum of the fermionic partial state is negligible compared with the contribution of the homogeneous FLRW state. The self-consistency of this approximation can be explicitly checked using the constraint equation [73].
If these approximations hold within our perturbative treatment, then the constraint equation can be recast as the following Schrödinger-like one: where the hat denotes the hybrid loop representation of the underlying function, and [37,38] C (Γ) This function C (Γ) D encodes, in mean value, how much the partial FLRW state Γ differs from being an exact solution of the unperturbed system. In this sense, it can be understood as a quantum backreaction term between the fermionic sector and the homogeneous background. Besides, it is worth mentioning that the generator of the fermionic evolution dictated by this Schrödinger equation is a φ-dependent operator that only acts on F D , and depends on the homogeneous background geometry by means of expectation values of FLRW operators on the partial state Γ. Furthermore, this operator is automatically well defined on the vacuum, given the preferred Fock quantization adopted for the fermionic degrees of freedom. Therefore, the introduced approximations in the hybrid quantum constraint equation allow one to recover a QFT regime in a quantum spacetime, with good physical properties.
Solutions to the Schrödinger equation (102) can be obtained by considering its associated Heisenberg dynamics. For that purpose, it is convenient to introduce the following change of evolution parameter: of our hybrid quantization. In addition, d η Γ denotes the derivative with respect to the parameter η Γ . Besides, where the diagonal mode coefficients h k D and the interaction mode coefficients h k I of the fermionic Hamiltonian are respectively given in Equations (94) and (96) (up to the elimination of any dependence on the inflaton momentum by identifying H |0 with zero, as we have explained). It is worth noting that, strictly speaking, these Heisenberg equations for the fermionic modes can be derived without the second approximation (ii) introduced above. In fact, given only the first approximation (i), it suffices that there exists a regime for all the modes in which the annihilation and creation operators find a straightforward counterpart in the Grassmann variables that they represent [36,73].
The resulting evolution from time η 0 to time η of the annihilation and creation operators is a Bogoliubov transformation that can be easily seen to take the form [37] a (x,y) k where α k (η 0 , η 0 ) = 1, β k (η 0 , η 0 ) = 0, and, for all η, In particular, recall that we are focusing our attention on phase space splittings and Fock representations of the fermionic degrees of freedom characterized by Equations (97) and (98) in the asymptotic regime of large ω k . These lead to interacting contributions to the fermionic Hamiltonian such that the beta coefficients of this transformation have the asymptotic behavior , where the function Max[., .] picks out the argument of dominant asymptotic order [38]. Taking into account the second condition in Equation (98), it is clear that these beta coefficients form a square summable sequence, over all k. It then follows that our quantum Heisenberg dynamics is unitarily implementable on the fermionic Fock space.
The unitary operator that implements the considered Heisenberg evolution can be explicitly constructed [37,38]. Furthermore, it can be checked that its action on the fermionic vacuum state picked out by the hybrid quantization provides a solution to the Schrödinger equation (102) if the backreaction function is of the form Here, c (x,y) k are arbitrary real phases of the evolution operator and is certain function with an asymptotic behavior which depends on that of the interaction contribution G (Γ) k (or, equivalently, on that of β k ). For our allowed families of Fock representations of the fermionic degrees of freedom, one can see that this function is of dominant asymptotic order O[Max(ω k |θ k | 2 , ω −5 k )]. Therefore, our backreaction term relating the FLRW background and the fermionic perturbations in the regime of QFT in a quantum spacetime, attained by the hybrid quantization of the system, turns out to be an absolutely convergent quantity. Thus, one needs not perform any regularization or substraction of infinities to render this backreaction term finite, something that could be done by using the arbitrary phases c (x,y) k that one can freely add to the evolution operator, but that would imply an unjustified adjustment of an infinite number of quantities.
We end this section with the following remark. The first work [37] that was carried out about the introduction of fermionic perturbations in hLQC employed a choice of annihilation and creation variables analogous to that introduced by D'Eath and Halliwell in Ref. [62]. As mentioned in the previous section, this choice is within the family of unitarily equivalent quantizations that allow for a unitarily implementable dynamics, in the context of QFT in curved spacetimes. The analysis in Ref. [37] shows that, with this choice of fermionic variables, the backreaction term C (Γ) D in the hybrid quantum theory fails to be an absolutely convergent quantity and needs to be regularized. This is an example of how the requirement of a unitarily implementable evolution alone is not enough to guarantee nice properties of the hybrid quantization of the full system. In fact, the conditions that are needed for the absolute convergence of C (Γ) D are very similar, but slightly weaker, than those for a proper definition of the fermionic Hamiltonian on the vacuum. Indeed, it is enough that Equation (98) is fulfilled and ω k |θ k | 2 is a summable sequence over all k [38]. In this respect, any choice of fermionic variables that leads to a good definition of the Hamiltonian operator on the vacuum automatically guarantees that the quantum fermionic backreaction is finite, without the need of regularization.

Fermionic Hamiltonian diagonalization: Choice of vacuum state
The restrictions on the definition (73) of fermionic annihilation and creation variables in hybrid quantum cosmology to guarantee: (i) unitarity of the QFT evolution [Equation (97)], and (ii) a well-defined Hamiltonian on the vacuum [Equation (98)], allow us to considerably reduce the possible choices of such variables, at least in the asymptotic regime of large ω k . However, and even though all of the resulting Fock representations are unitarily equivalent, there is still much freedom in identifying a particular set of annihilation and creation variables (even in the asymptotic regime). In other words, there remains ambiguity in the complete characterization of the phase space splitting between the FLRW background and the fermionic sector, as well as of the fermionic vacuum state of the theory. In this section we motivate and adhere to the physical criterion of Hamiltonian diagonalization to try and fix this remaining choice, following a procedure that is specifically adapted to the spatially local structure of the fermionic dynamics.

Hamiltonian diagonalization in hLQC
The previously imposed criteria on the allowed choices of fermionic annihilation and creation variables [Equations (97) and (98)] are tailored to have the net effect of diminishing the dominant asymptotic order, in the regime of large ω k , of the interaction parts h k I of the fermionic Hamiltoniañ H D [given in Equation (93)]. Ultimately, this fact is responsible for the unitarity of the fermionic evolution, as well as for a proper definition of the fermionic Hamiltonian on the vacuum. As we have already commented, these interaction parts annihilate and create pairs of particles and antiparticles in the quantum fermionic dynamics. From a physical perspective, one would think that a choice of phase space splitting where the assignment of the dynamical contribution of the FLRW background to the system is such that the fermionic states undergo no annihilation and creation of pairs, would be one that is naturally adapted to the dynamics of the entire cosmological system. This choice should then be such that h k I = 0, so that the Fock quantization of the resulting fermionic HamiltonianH D would have a diagonal action on the n-particle states associated with the selected set of annihilation and creation operators.
According to this line of reasoning, we refer to any choice of fermionic annihilation and creation variables that lead to vanishing interaction terms h k I , for all k, as variables for the Hamiltonian diagonalization. In general, restricting to cases with f k 2 = 0 (which include all those choices that respect the standard convention for particles and antiparticles in the massless limit [55]), one can readily check that the diagonalization condition h k I = 0 is fulfilled for all k if and only if where f k 1 = f k 2 ϕ k . This is a semilinear partial differential equation which has locally unique solutions, as long as the section of initial conditions is transversal to the flow of the Hamiltonian vector field {., H |0 } [74]. Naturally, there are several possible families of such solutions for all k, and each of them completely characterizes (up to the two phases G k and F k 2 ) a different set of fermionic variables for the Hamiltonian diagonalization, in virtue of Equation (89). Explicitly, it holds that Let us point out that, in fact, the structure of the differential Equation (110) allows for solutions that, besides on the FLRW geometry, can depend also on the inflaton and its canonical momentum. Remarkably, nonetheless, the definition of fermionic annihilation and creation variables resulting from any such choice of coefficients for the Hamiltonian diagonalization can still be completed to become canonical in the entire cosmological system, following an analogous procedure to that discussed in the previous section. We adopt this extended framework for the definition of the fermionic variables in this and the next subsection.
Finally, it is worth noticing that relation (111)

Asymptotic diagonalization
In the following, we try and fix a preferred solution to Equation (110), using our previous knowledge on the restrictions that unitary evolution and a proper definition of the fermionic Hamiltonian impose, and by looking into the details of this Hamiltonian in the asymptotic regime of large ω k . For fermionic variables that admit a unitarily implementable evolution, namely when Equation (97) holds, the interaction coefficients h k I of the fermionic HamiltonianH D behave asymptotically as where the second summand results from the second Poisson bracket in Equation (96). As we pointed out above, we explicitly see here that demanding condition (98) for a well-defined action of the Fock quantization of the Hamiltonian on the fermionic vacuum is equivalent to diminishing the dominant asymptotic order, in inverse powers of ω k , of the interaction coefficients. Once this condition is imposed, one arrives at an asymptotic behavior for h k I with an analogous structure as in Equation (113). More concretely, its dominant contribution is given by 2ω k θ k /a plus certain specific terms that are proportional to ω −2 k . One can cancel this contribution again by conveniently fixing the dominant asymptotic behavior of θ k . The resulting interaction coefficient h k I displays, once more, a similar asymptotic structure, but with the role of θ k played by its subdominant contributions and the remaining summands having a lower asymptotic order in powers of ω k . Owing to the asymptotic structure of the Hamiltonian, this pattern repeats itself at each asymptotic order in inverse powers of ω k , if one admits an asymptotic expansion of this form and imposes the cancellation of the previous dominant terms in the interaction coefficient h k I . Motivated by these properties of the fermionic Hamiltonian, and taking into account that asymptotically if the unitary evolution condition (97) is fulfilled, we propose the following asymptotic series as an ansatz for a Hamiltonian diagonalization in the regime of large ω k : Here, the symbol ∼ indicates the equality of the asymptotic expansions, and γ n are functions of the homogeneous FLRW background canonical variables. These are completely fixed in an iterative way if one introduces our ansatz in the interaction coefficients h k I of the fermionic Hamiltonian, and imposes that each contribution in inverse powers of ω k be equal to zero. Specifically, one then obtains [39,75] γ n+1 = a H |0 , γ n + ma where γ −n ≡ 0 for all n > 0. This is a deterministic recurrence relation that can be used to fix all of the functions γ n , starting from the initial datum γ 0 . It is worth mentioning that the proposed function ϕ k for the asymptotic diagonalization of the fermionic Hamiltonian, given by Equations (115) and (116), provides a very specific solution to the general equation (110), in the sector of unboundedly large ω k . Such an asymptotic solution can be thought of as a physically preferred solution, inasmuch as it has been obtained by exclusively adhering to local features of the fermionic Hamiltonian. However, despite the strong asymptotic restriction that our requirement sets on the admissible solutions to Equation (110), there may exist many such solutions for all k that, viewed as functions of ω k , display the same, preferred asymptotic behavior. It seems therefore most convenient to investigate whether the imposition of certain smoothness conditions on the dependence of any such ϕ k on ω k (e.g., continuity or analiticity) can allow us to fix this solution completely. Indeed, if we were able to ensure this uniqueness, by taking into account relations (89), we would solve the last remaining ambiguity in the choice of fermionic variables for the hybrid description of the system, up to the phases G k and F k 2 . In other words, once these phases were chosen, we would succeed in specifying a Fock representation (with its associated vacuum state) of the fermionic degrees of freedom, together with a particular splitting of the fermionic and FLRW sectors of the phase space and of their contribution to the dynamics of the entire system. Actually, we have already restricted the phase G k in a substantial way by demanding it to be a dynamically irrelevant constant in the classical linearized system, after imposing {G k , H |0 } = 0 in order to eliminate asymmetries in the evolution of particles and antiparticles. As for the other phase, F k 2 , it can be naturally selected by demanding that the original dynamics that is extracted from the Dirac field by means of the background-dependent transformation (73) be minimal, along the lines detailed in Ref. [40].
We end this subsection by noting that, if one specifies a preferred choice of fermionic annihilation and creation variables for the Hamiltonian diagonalization, then one is not only fixing the relevant structures for the hybrid quantization of the cosmological system, but also a concrete Fock representation of the Dirac field, within the framework of QFT in curved spacetimes. This regime is attained when the homogeneous background obeys the Friedmann equations, whereas the fermionic annihilation and creation variables evolve with the dynamics dictated by the fermionic contribution to the Hamiltonian (92). If the interaction terms in this Hamiltonian are zero, then this fermionic dynamics can be straightforwardly solved, namely, the annihilation and creation variables just evolve via multiplication of their initial data (at an arbitrary initial time) by a complex phase. With these solutions at hand, if one takes the inverse of the defining transformation of these variables given by Equation (73), and introduces it in the mode expansion of the Dirac field, one immediately obtains a complete basis of solutions for the Dirac equation, in the sense of Equation (7). The constant coefficients of the elements of this basis in the expansion of the field are the annihilation and creation initial data, that select a unique Fock representation (and its associated vacuum state) once they are promoted to operators.

Uniqueness of the vacuum: Minkowski and de Sitter spacetimes
In this subsection we focus our discussion on the aforementioned regime of QFT in curved spacetime, and explain how our ansatz for asymptotic diagonalization succeeds in the selection of natural vacuum states in Minkowski and de Sitter spacetimes. In fact, the asymptotic expansion given in Equations (115) and (116) allows us to determine a complete basis of solutions for the Dirac equation [as in relation (7)] that turns out to correspond to the choice of the Poincaré or the Bunch-Davies vacuum, respectively, when the background cosmology is fixed as the Minkowski or the de Sitter spacetime. Taking into account that, when the homogeneous background obeys the Friedmann equations, we have that a{., H |0 } is simply the derivative with respect to the conformal time, in the considered situations in QFT the recurrence relation (116) becomes Let us start by considering the case of a background given by the classical Minkowski spacetime. This particularization is easily implemented by setting the scale factor as the unit constant, the inflaton as an arbitrary constant, and its potential equal to zero. Then, we immediately have that γ 0 = m, while any other γ n , determined by the recursion relation (117), has a vanishing time derivative. That iterative equation can be solved by introducing the generating function G(x) = ∑ ∞ n=0 γ n x n , that leads to a quadratic equation with only one solution consistent with the initial datum γ 0 = m: Around x = 0, this is an analytic function with power series in x characterized by the coefficients γ n , by construction. Then, comparing this series to our ansatz (115) and employing the uniqueness of the asymptotic expansion, we can directly identify the function ϕ k that leads to an asymptotic diagonalization with the following analytic function: Using the relations (89) and (111), that arise from the requirement that the fermionic annihilation and creation variables be defined by means of a canonical transformation, one eventually obtains where ξ k = ω 2 k + m 2 . Since the selected function ϕ k that permits the Hamiltonian diagonalization is completely independent of time in this case, according to our comments above, it is natural to demand that the phase F k 2 be simply an arbitrary constant, as well as G k . From Equation (112), one can straightforwardly check that the diagonal coefficients of the resulting fermionic Hamiltonian are then h k D = ξ k . A simple inspection of this result, together with Equation (120), immediately reveals that our criterion of asymptotic diagonalization selects indeed the basis of solutions to the Dirac equation in Minkowski spacetime that corresponds to the Poincaré Fock representation of the field. Namely, it is the quantization that, with a standard convention for particles and antiparticles, separates between positive and negative mass-shell frequencies ξ k .
Let us show, in addition, how our criterion of asymptotic diagonalization recovers the common notion of Bunch-Davies vacuum for the Dirac field in de Sitter spacetime. In a flat slicing, this spacetime can be understood as a cosmological solution of Friedmann equations obtained by setting the inflaton potential equal to the constant 3H 2 Λ /(8π) and the inflaton momentum equal to zero. Here, H Λ is the constant Hubble parameter. In conformal time, the expanding scale factor then behaves as The analysis of the restrictions that the asymptotic diagonalization imposes on the choice of a fermionic vacuum is easier if one first considers the general differential equation for ϕ k in our de Sitter background. It reads The general solution of this equation can be found by introducing a mode-dependent complex time T k = −2iω k η and the following change of variables [75]: The function v k turns out to satisfy a confluent hypergeometric equation in the complex variable T k , that has the general solution where A and B are integration constants, and 1 F 1 (.; .; z) is the hypergeometric function of type (1, 1), that is absolutely convergent for all values of its complex argument z [76]. The asymptotic expansion determined by our criterion of Hamiltonian diagonalization, given in Equations (115) and (117), actually picks out a specific solution of the mentioned hypergeometric equation, up to an irrelevant multiplicative constant, namely a particular ratio between the integration constants A and B. Indeed, the iterative relation (117), particularized to our de Sitter background, can be seen to lead to coefficients γ n such that where C n are constants that are completely specified by a complicated nonlinear recurrence relation, with the initial value C 0 = mH −1 Λ [75]. We need not solve this relation, since the T k -dependence of the above asymptotic series, together with the known value of C 0 , is enough to restrict the associated expansion of v k in Equation (123) By introducing this ansatz for the asymptotic behavior of v k in the confluent hypergeometric equation that it must satisfy, one can determine the coefficients v n of the asymptotic expansion exclusively in terms of v 0 , yielding where 2 F 0 (., .; −; z) is the hypergeometric function of type (2, 0), that has a zero radius of convergence. Even though it formally diverges, its series is known to provide the asymptotic expansion of a very particular type of solution to the confluent hypergeometric equation, namely the Tricomi solution [76]. In fact, using the asymptotic properties of the hypergeometric functions in Equation (124), it is possible to prove that the Tricomi solution is the unique one that has an asymptotic expansion of the form (127). Therefore, after substituting this solution in relation (123), we immediately see that our procedure of asymptotic Hamiltonian diagonalization allows us to obtain again a unique solution to the general equation (110), for all wave vectors k.
In more detail, the Tricomi solution for v k selected by our criterion of asymptotic diagonalization leads, after several manipulations, to the following function ϕ k in Equation (123) [75]: where H (1) ν denotes the Hankel function of the first kind [77]. In turn, using relations (89) and (111), this result determines the coefficients f k 1 , f k 2 , g k 1 , and g k 2 that define the fermionic annihilation and creation variables, up to the phases F k 2 and G k . We recall that the latter is fixed as an arbitrary constant. As for the former, namely F k 2 , one can select it following the ideas put forward in the previous subsection, so that it minimizes the amount of dynamics extracted by the time-dependent canonical transformation (73). In any case, the details about the time dependence of this phase are irrelevant for the basis of solutions to the Dirac equation that ϕ k selects (in the context of QFT in a fixed curved spacetime), as one can straightforwardly check using Equations (73) and (112). Then, after taking into account the diagonal evolution of the annihilation and creation variables dictated by h k D in Equation (112), the resulting complete set of solutions for the Dirac equation in which the field decomposes is given by very specific linear combinations of Hankel functions of the first (in the case of antiparticles) and second (in the case of particles) kinds, different for each chirality and helicity [75]. We recall that any such basis decomposition fixes a particular Fock representation of the field. In this case, the basis of solutions turns out to be precisely the one that has been naturally associated in the literature with the fermionic analog of the Bunch-Davies vacuum in de Sitter spacetime [75,78,79]. Therefore, our criterion of asymptotic Hamiltonian diagonalization provides again the vacuum state that is physically accepted as preferred, in this case in a de Sitter background.
We end this section with a final remark. Beyond the well-known background spacetimes analyzed here, namely Minkowski and de Sitter, quantum fields in FLRW cosmologies suffer from the lack of a natural choice of vacuum state, at least if one appeals only to the symmetries of the system to select it. In the case of scalar fields, a common approach to mitigate this issue is the introduction of adiabatic states (see e.g. Refs. [12,[80][81][82]). Their construction is based on an iterative procedure to solve the field equations such that the resulting quantization displays certain local Poincaré-like features. In the case of Dirac fields, there have been at least two notable attempts to generalize the notion of adiabatic states [83,84]. The proposal in Ref. [84] follows closely the construction procedures previously established for scalar fields. In particular, this work introduces an algorithm to iteratively solve the Dirac equation that, at each consecutive step, is able to approach the actual mode solutions up to contributions that are more and more subdominant in the asymptotic regime of large ω k . An adiabatic state of nth-order is defined by truncating this procedure at the nth-step of the iteration and setting the value of the resulting approximate mode solutions, at an arbitrary time, as the initial data that specify the basis of solutions for the Fock representation of the field. Ref. [85] analyzed the relation between these adiabatic states and our family of unitarily equivalent quantizations of the Dirac field, including those that satisfy the condition of asymptotic diagonalization. It was shown that the representations associated with adiabatic states of all orders belong to the same equivalence class. In particular, the zeroth-order state already corresponds to a representation that admits unitarily implentable evolution, once the time dependence attributed to the FLRW background has been conveniently extracted. In fact, this unitarity guarantees that any two states defined with adiabatic initial data at different times are unitarily equivalent, so that the choice of initial time for the definition of the adiabatic states is not a relevant ambiguity. Furthermore, the first-order adiabatic state directly leads to a Fock quantization of the Dirac field in the family selected by imposing that the fermionic Hamiltonian for the annihilation and creation variables be well defined on the vacuum. The question of whether higher-order adiabatic states give rise to representations that behave, in the asymptotic regime of large ω k , increasingly closer to the one(s) selected by our criterion of asymptotic Hamiltonian diagonalization is yet an open issue.

Conclusions
In this work, we have reviewed some recent investigations, carried out by us and our collaborators, about the physical motivation and use of certain criteria capable to ensure the uniqueness of the Fock quantization of fields in cosmological systems, specialized to the case of fermions described by Dirac fields. The presented results have been applied to the study of the hybrid quantization of the primordial Universe with perturbations, that contain all the fermionic degrees of freedom described by a Dirac field (and may also include other matter field perturbations and metric perturbations).
We have first considered the Fock quantization of the CARs for Dirac fields in conformally ultrastatic three-dimensional spacetimes, as well as in cosmological FLRW spacetimes in four dimensions, with spherical or toroidal spatial hypersurfaces. We have characterized the set of vacua that are invariant under the physical symmetries of the Dirac equation in these spacetimes. These symmetries include the continuous isometries of the spatial hypersurfaces, enlarged with the spin rotations generated by the helicity in the case of FLRW cosmologies with sections of toroidal topology.
For all the Fock representations associated with the above set of invariant vacua, we have proven that there exists a subset that admits unitary implementability of the dynamics on the Fock space. This evolution comes from the Dirac equation, after extracting from the fermionic field some of its time variation that can be attributed to the dependence on the variables that describe the spacetime background. In the Heisenberg picture, the extracted part is regarded as explicitly time dependent, and therefore is not included in the proper quantum dynamics of the annihilation and creation operators. In fact, this extraction is necessary in order to achive the unitary implementability of the quantum evolution. It must be restricted, nonetheless, by the condition that the evolution remaining from the original Dirac equation be not trivialized.
After determining all the Fock representations that are allowed by the criteria of invariance under the symmetries of the equations of motion and of a nontrivial unitary implementability of the dynamics, we have shown that all these representations are unitarily equivalent for each of the spacetime scenarios that we have considered, provided that one fixes a convention to distinguish between particles and antiparticles of the Dirac field. In other words, our well-motivated conditions of unitarity and invariance guarantee the uniqueness of the Fock representation, up to unitary equivalence.
This uniqueness result has the immediate consequence of specifying also a unique concept of quantum dynamics for the fermionic annihilation and creation operators, modulo unitary redefinitions. Indeed, our analysis allows us to fully characterize the functions of the spacetime background that need to be removed from the time dependence of the fields, at least in the ultraviolet limit of large eigenvalues of the Dirac operator on the spatial hypersurfaces. This characterization can alternatively be understood as the determination of which field excitations are the particles and antiparticles that preserve their coherence over time.
We have provided an optimal description, with an eye to its quantization, of the phase space of a homogeneous and isotropic cosmology coupled to a homogeneous scalar field (that acts as an inflaton in General Relativity) and with fermionic perturbations, when the Einstein-Dirac action is truncated at quadratic perturbative order. For the fermionic sector of the phase space, we have used our previous result about the Fock representation of a Dirac field in order to select a quantization of the fermionic degrees of freedom, up to unitary modifications. Hence, for the fermionic field, we have chosen certain annihilation and creation variables that are related with the Dirac modes through a canonical transformation that depends on the homogeneous and isotropic background, and that supports a unitarily implementable Heisenberg dynamics when the background is viewed as a fixed entity. This leads to a specific splitting of the phase space between the background degrees of freedom and the fermionic content. A particular consequence is the modification of the contribution to the global Hamiltonian constraint associated with the fermionic perturbations. We have taken advantage of this modification and, going beyond the criterion of unitary dynamics for the selection of the fermionic variables, we have employed the remaining freedom in the background dependence of this choice to obtain other desirable properties in the quantization of our system. One property that we have investigated is a proper definition of the fermionic Hamiltonian operator on the set of finite particle/antiparticle states constructed from the vacuum. On the other hand, the discussed splitting of the phase space implies also a change in the canonical variables that describe the homogeneous background, so as to preserve the symplectic canonical structure of the system at the perturbative order of our truncation. The corresponding change in the background variables amounts to the correction of the original ones with terms that are quadratic in the perturbations.
Using the above description of the phase space, the only nontrivial constraint that needs to be imposed quantum mechanically is the zero mode of the Hamiltonian constraint. This global constraint interrelates the different physically relevant sectors of the phase space, namely, the geometric FLRW sector, for which a loop representation is adopted, the inflaton, with a Schrödinger-like representation, and the fermionic perturbations, for which one takes a Fock representation in the selected family (in addition, it is possible to include scalar and tensor perturbations, described by perturbative gauge invariants, with Fock representations that can be picked out as well with our proposed criteria). In order to single out this preferred family of Fock representations, in addition to the invariance under the symmetries of the field equations and the unitary implementability of the Heisenberg dynamics, we have chosen a convention for the distinction between particles and antiparticles which smoothly connects with the standard convention of QFT when the mass of the field vanishes.
We have shown how to impose the operator that represents the zero mode of the Hamiltonian constraint of our perturbed cosmological system on states for which the dependence on the different sectors of the phase space, except the inflaton, becomes separable. In this way, we have been able to find mild conditions under which the imposition of the constraint turns out to be essentially equivalent to a certain master constraint equation on the fermionic perturbations. Given our perturbative hierarchy, these conditions amount to have negligible FLRW geometry transitions mediated by the zero mode of the Hamiltonian for the partial wave function that describes such an homogeneous geometry in our state. The resulting equation is special inasmuch as the dependence on the FLRW geometry only persists by the inclusion of expectation values over that geometry. With an additional approximation on the variation of the partial wave function of the fermionic perturbations with respect to the inflaton, that can be checked at least for self-consistency, one can deduce from this master constraint on the perturbations a Schrödinger equation for the partial state that describes the fermionic content. This Schrödinger equation involves the quantum backreaction that the fermionic perturbations produce on the FLRW geometry.
Besides, within our approximations, it is possible to solve the quantum dynamics dictated by the commented master contraint equation on the perturbative fermionic modes. We have reviewed how such dynamics can be implemented on our Fock space. The resulting evolution depends, in particular, on the FLRW geometry, but this is so exclusively through expectation values, that turn out to be different for each fermionic mode and that are well defined thanks to the loop representation adopted in the scheme of hLQC. These facts, together with ultraviolet properties, guarantee the unitary implementability of the fermionic Heisenberg dynamics. One can construct the associated unitary evolution operator, that is generated by the fermionic Hamiltonian that appears in the Schrödinger equation that has been derived. Actually, we have seen that the requirement that this Hamiltonian be well defined on the vacuum is enough to guarantee a finite fermionic backreaction on the FLRW background, without the need of any regularization. On the other hand, the unitarity of the fermionic dynamics translates into a finite production of pairs of particles and antiparticles in the evolved vacuum.
We have gone one step beyond and employed the still remaining freedom in the determination of the Fock representation of the fermionic degrees of freedom, and their splitting from the FLRW sector of the phase space, to demand an additional feature in the fermionic Hamiltonian, namely, that it become diagonal in terms of the fermionic annihilation and creation variables in the asymptotic region of large wave numbers of the modes, in the sense that it do not contain interactions in that region that produce pairs of particles and antiparticles. We have seen that this condition indeed fixes asymptotically the choice of vacuum state. Furthermore, we have argued in favor of the uniqueness of the vacuum selected by means of this asymptotic Hamiltonian diagonalization when extended to all wave numbers by suitable smoothness conditions. In this respect, we have demonstrated the uniqueness in the case of standard QFT in Minkowski and de Sitter spacetimes, treated as fixed backgrounds, showing in addition that the vacua that are picked out by the diagonalization procedure are the Poincaré and the Bunch-Davies vacua, respectively. For more general backgrounds, either of classical or quantum nature, our proposal can potentially serve to attain a well-defined and unique choice of vacuum state with especially good physical and mathematical properties.
Finally, we have commented on the relation between adiabatic states and the vacua selected by our criteria. For iterative constructions of fermionic adiabatic states, all of them turn out to belong to the unitary equivalence class of Fock states that incorporate symmetry invariance and allow for a unitarily implementable dynamics. Besides, states of first or higher adiabatic order belong to the family of Fock states picked out by the additional requirement of a fermionic Hamiltonian with a well-defined action on the dense set of finite particle/antiparticle states, and therefore it is ensured that they lead also to a finite fermionic backreation. As for the issue of Hamiltonian diagonalization, it is an open question whether higher-order adiabatic states give rise to representations in which the vacuum state increasingly approaches our choice in the asymptotic regime of large wave numbers.