Evolution of confined quantum scalar fields in curved spacetime. Part I

We develop a method for computing the Bogoliubov transformation experienced by a confined quantum scalar field in a globally hyperbolic spacetime, due to the changes in the geometry and/or the confining boundaries. The method constructs a basis of modes of the field associated to each Cauchy hypersurface, by means of an eigenvalue problem posed in the hypersurface. The Bogoliubov transformation between bases associated to different times can be computed through a differential equation, which coefficients have simple expressions in terms of the solutions to the eigenvalue problem. This transformation can be interpreted physically when it connects two regions of the spacetime where the metric is static. Conceptually, the method is a generalisation of Parker's early work on cosmological particle creation. It proves especially useful in the regime of small perturbations, where it allows one to easily make quantitative predictions on the amplitude of the resonances of the field, providing an important tool in the growing research area of confined quantum fields in table-top experiments. We give examples within the perturbative regime (gravitational waves) and the non-perturbative regime (cosmological particle creation). This is the first of two articles introducing the method, dedicated to spacetimes without boundaries or which boundaries remain static in some synchronous gauge.


Introduction
Quantum field theory in curved spacetime is the theory that studies the evolution of quantum fields which propagate in a classical general relativistic background geometry. To this date, the theory has been able to yield many quantitative theoretical predictions for physically relevant problems within its scope, such as cosmological particle creation [1,2], the Unruh effect [3] or Hawking radiation [4]. The success in the study of those concrete well-known problems can be related to the use of different mathematical techniques and simplifications, adapted to the specific problem and which yield the computations tractable. For instance, in the study of phenomena such as Hawking radiation or the Unruh effect, the presence of horizons (or "would-be horizons" [5,6]) allows to simplify the computations leading to a thermal spectrum for the quantum field, at least within most of the usual approaches to these phenomena. Also, the presence of symmetries like homogeneity or isotropy is convenient to obtain quantitative results on particle creation in cosmological models.
In the recent years, the interest to verify the theory experimentally has increased significantly [7][8][9][10][11][12][13]. Thus far, an analogous of Hawking radiation in Bose-Einstein condensates has already been experimentally demonstrated [14] (although the result is not free of controversy, see for example [15]). In particular, and due to the great improvement of the precision of quantum measurements in table-top experiments, it is expected that the theoretical predictions on quantum fields confined in cavities, and under the effect of small changes in the background geometry or the non-inertial motion of the cavity, will be tested experimentally in the near future [16][17][18][19]. Consequently, new open problems within the framework of the theory are arising, related to the new concrete systems under study, and therefore new arXiv:1811.10507v6 [quant-ph] 25 Oct 2021 adapted mathematical techniques will be necessary to approach them.
In this work, we introduce a general method for computing the evolution of a confined quantum scalar field in a globally hyperbolic spacetime, by means of a time-dependent Bogoliubov transformation. As it will be clearly shown along the article, the method proves to be especially useful to address the new kind of problems that we just mentioned, related to quantum fields confined in cavities and undergoing small perturbations, since it provides very simple recipes for computing the resonance spectrum and sensibility of the field to a given perturbation of the background metric or the boundary conditions. However, we wish to emphasise that the method can be applied to any confined quantum scalar field, both real or complex, evolving as a result of any kind of regular changes in the background metric and/or the boundary conditions, just under some minor assumptions. This generality will also be illustrated through a concrete example of an application of the method outside the scope of small perturbations (in cosmological particle creation) that will be considered in the article.
The method is based on the foliation of the spacetime in spacelike Cauchy hypersurfaces using a time coordinate. Since the method is intended for confined fields, these hypersurfaces will be compact, although not necessarily with boundaries. The core idea is to construct a basis of modes naturally associated to each hypersurface, and then provide a differential equation in time for the Bogoliubov transformation between the modes associated to two hypersurfaces at different times. This way, the evolution of the field in time is not obtained by solving the Klein-Gordon equation, but rather by solving a differential equation for the time-dependent Bogoliubov transformation.
The core idea of transferring the dynamics from the mode functions to the Bogoliubov coefficients appears for the first time in the pioneer work by Parker [2]. In this sense, our work is a conceptual extension of Parker's work to the general case of any confined quantum scalar field (with minor assumptions). However, the key mathematical construction has been developed as a generalisation of the technique introduced in [20]. That previous technique, although with the strong limitation of being applicable only in Minkowski spacetime (the influence in the evolution of the field coming from the motion of the cavity), has already proven to be of special practical use in the study of quantum fields confined in cavities. Benefiting from conformal invariance, extensions to a few other metrics were also considered [17,[21][22][23][24]. Yet, the scope of that method was notoriously reduced, something which the present work overcomes.
We want to emphasise that our way of proceeding makes the method completely different to the Hamiltonian diagonalisation approach introduced for the problem of cosmological particle creation [25,26]. In such approach, the construction of bases of modes associated to each time is used to diagonalise the Hamiltonian of the quantum field by a time-dependent Bogoliubov transformation in the Fock representation, such that the Hamiltonian is diagonal in terms of particle operators at every moment; and then this diagonalised Hamiltonian is used to compute the evolution of a state in time directly in the Fock representation. Therefore, this approach introduces a different notion of particles at each instant of time. Problems with this approach appear because the Fock representations for different times can be unitarily inequivalent [27][28][29]; that is, because the particle interpretation in this framework is in general wrong. Conversely, the method introduced in this article does not intend to quantise the field with a Fock representation at all times during its evolution. The construction of the modes at each time will be just an intermediate mathematical tool used to compute a time-dependent Bogoliubov transformation. Such Bogoliubov transformation has the usual physical interpretation in terms of effects on quantum particles (mode-mixing and particle creation) only when it relates regions of spacetime where the Fock quantisation has a valid physical interpretation in terms of particles, in particular due to the presence of a timelike Killing field. The Bogoliubov transformation then unambiguously describes the evolution of the field between the two regions, since it allows to compute the scattering matrix between the bases of the two different Fock quantisations [29][30][31]. We provide strong physical arguments to explain why the Fock representations that we consider are unitary equivalent and physically meaningful.
Following the above discussion, we take the opportunity to stress that, for the quantisation of the field, we simply rely on canonical Fock quantisation in the regions where it is clearly physically meaningful. In other words, the method is of general practical applicability to compute the scattering matrix between two physically valid Fock representations of the field, and it is indeed an extremely powerful tool for such purpose; while, on the other hand, we have no intention of providing a fully consistent field quantisation in any region of a general spacetime (see the approaches to such problem in e.g. [32][33][34][35][36]).
This article constitutes the first of two articles introducing the method. In this first part we consider spacetimes without boundaries, or with timelike bound-aries that remain static in some synchronous gauge, as explained in the next section. Spacetimes with timelike boundaries which do not remain static in any synchronous gauge require an specific treatment, which due to its extension we leave for a second forthcoming article. The present article is organised as follows. In Sect. 2 we state the general physical problem for which we will construct the method, introducing the background metric, the field theory and the different assumptions that we consider; and also define two important mathematical objects that we will use. The nuclear part of the article is Sect. 3. In this section we construct the basis of modes associated to each hypersurface of the foliation of the spacetime, and formally compute the time-dependent Bogoliubov transformation between the modes of two different hypersurfaces, for which we give a differential equation and a formal solution. We also discuss the physical meaning of both the modes and the transformations. In Sect. 4 we consider the particularly important case of small perturbations and resonances, obtaining especially simple recipes for its solution. As an example of application, we consider the problem of a confined quantum field in the presence of a gravitational wave. In Sect. 5 we apply the method to an example of cosmological particle creation. Finally, in Sect. 6 we present the summary and conclusions. In addition, in Appendix A we discuss in detail the different conditions that we require for the method to work. In Appendix B we describe how the method handles the problems in which modes with exponential time dependence or zero-frequency modes are present. In Appendix C we prove that the bases of solutions to the Klein-Gordon equation on which we build our method are complete and orthonormal. In Appendix D we provide the detailed computation of the differential equation provided in Sect. 3, and prove that the Bogoliubov transformations obtained satisfy the Bogoliubov identities. In Appendix E we further support and extend the results obtained on small perturbations. In Appendix F we provide for convenience a summary of the useful formulae for the application of the method.

Statement of the problem
We consider a globally hyperbolic spacetime (M, g) of dimension N + 1, possibly with timelike boundary ∂M [37]. In this spacetime we define a Klein-Gordon scalar field Φ satisfying the equation where m ≥ 0 is the rest mass of the field, g µν is the spacetime metric, R its scalar curvature and ξ ∈ R is a coupling constant (we use natural units = c = 1). In the case of the existence of boundaries, we impose one of the following four possible boundary conditions to the field: a) Dirichlet vanishing boundary conditions b) Neumann vanishing boundary conditions where n µ (x) is the normal vector to ∂M . c) Robin vanishing boundary conditions where γ ∈ R \ {0} is the boundary parameter. Later on, we will introduce the possibility that this parameter is replaced by a function.

d) Mixed vanishing boundary conditions
In this case, Dirichlet, Neumann and Robin vanishing boundary conditions apply each to complementary regions of the boundary. Being just a combination of the other cases, we mention them here as a possibility, but for simplicity we will not consider them explicitly when constructing the method.
The types of boundary conditions considered are by far the most common ones in physical problems. We do not discard that the method could also apply to other less used conditions, but we simply do not approach such possibility within this work.
Thanks to the global hyperbolicity, it is always possible to construct a Cauchy temporal function t in the full spacetime [37], which provides a foliation in Cauchy hypersurfaces Σ t of constant time. We introduce the Klein-Gordon inner product between two solutions of (1), given by which, for convenience, we already evaluated at a given Cauchy hypersurface Σt, with dVt being its volume element. Both in the absence of boundaries, or under the boundary conditions (2)(3)(4), this inner product is independent of Σt.
Finally, we introduce the three conditions on the Cauchy hypersurfaces and the temporal function that we need to ensure the applicability of the method. These conditions are: A. The Cauchy hypersurfaces Σ t must be compact. B. For any Cauchy hypersurface Σ t , the Cauchy problem for the Klein-Gordon equation (1) must be wellposed; that is, given as initial conditions the value of the field and of its first derivative with respect to t at Σ t , there exists an unique solution to the Klein-Gordon equation in the whole spacetime satisfying these conditions. 1 C. Using the temporal function as a coordinate, the metric should be written as where h ij (t) is a regular Riemannian metric. 2 This is called a synchronous gauge.
The necessity of each condition will become clear when constructing the method. A detailed discussion on their physical meaning and the limitations they introduce can be found in Appendix A. In this first article introducing the method, we are going to consider the problems for which M either has no boundary, or its boundary remains parallel to the gradient of the temporal function t, that is, it remains static in some synchronous gauge. Problems for which there exists no adequate temporal function (satisfying the required conditions) such that the boundary remains parallel to its gradient will be considered in the second article.
2.2 Space of functions over Σ t and self-adjoint operator Let us introduce two mathematical objects that are pivotal for the method. First, we define Γ t as a subspace of the space of square integrable smooth functions over a Cauchy hypersurface, If M (and therefore Σ t ) has no boundary, then Γ t = C ∞ (Σ t )∩L 2 (Σ t ). If there are boundaries, then Γ t is the restriction of C ∞ (Σ t ) ∩ L 2 (Σ t ) to functions satisfying the boundary conditions in ∂Σ t compatible with the 1 Under the presence of boundaries, these initial conditions must be compatible with the boundary conditions at the intersection ∂Σ t = Σ t ∩ ∂M . 2 For simplicity in the notation, from here on we will omit the dependence of h ij (t) on the spatial coordinates, but we are considering the most general case in which this metric can have any (regular) dependence on the coordinates. We will not write the explicit dependence on the spatial coordinates for other quantities either, except when necessary. Also, in (6) for simplicity we are assuming that for each hypersurface Σ t there is one coordinate chart (x 1 , . . . , x N ) that completely covers it. This might not be the case, but considering several coordinate charts would be straightforward and would not affect the construction of the method, as far as for all of them the metric can be written as in (6). boundary conditions in ∂M ; which, since the boundaries are parallel to ∂ t , read for x ∈ ∂Σ t , where n( x) is the normal vector to the boundary ∂Σ t and ∇ h(t) is the covariant derivative corresponding to the spatial metric h ij (t). At this point, we mention that the boundary parameter γ could be also considered a function γ( x), but it must be independent of t. 3 Second, we define the operatorÔ(t) on Γ t aŝ where ∇ 2 h(t) is the Laplace-Beltrami differential operator, R h (t) the scalar curvature corresponding to the spatial metric h ij (t) and F (t) ≥ 0 a time-dependent non-negative quantity. The quantity F (t) is given by the condition that the operatorÔ(t) at each time must be positive definite in the corresponding space Γ t (any values of F (t) that meet this requirement would be valid). Notice that, thanks to condition A and eventually to the boundary conditions, this operator is selfadjoint [with respect to the usual scalar product in L 2 (Σ t )] and has a discrete spectrum with no accumulation points. 4 Moreover, if we fixed F (t) = 0 the operator would always have a minimum eigenvalue. Consequently, a finite value for F (t) such that the operator is positive definite can always be found. We refer to Appendix B for a discussion on the role of the quantity F (t), related to the presence of exponential time dependence of the solutions and of zero-frequency modes. 5 In many problems (such as the examples provided in this article) it is possible to fix F (t) = 0 while keepingÔ(t) positive definite, in which case the quantity can simply be ignored. 3 The justification for this apparently arbitrary limitation can be found in Appendix C. 4 Assertions about the self-adjoint nature of the operator and about the properties of its spectrum should strictly be done over the extension of the operator defined on Γ t to the full Hilbert space L 2 (Σ t ) (of which Γ t is a dense subspace). This is a well-known mathematical procedure in the analysis of partial differential equations with elliptic operators and boundary value problems, which is the context in which we will make use of the operator. Therefore, we shall not get into the details of it. Another possibility is to directly consider the operator as acting on a Sobolev space over the Riemannian manifold. See for example [38]. 5 The discussion in Appendix B cannot however be fully understood until reading Sect. 3.
Once we have introduced the operatorÔ(t), and making use of condition C, we shall simplify the Klein-Gordon equation (1) taking into account the form of the metric in (6), obtaining where with h(t) being the determinant of the spatial metric h ij (t) (it is a factor which depends on the change of the metric of the spacelike hypersurfaces with time), andR(t) := R(t) − R h (t) is the part of the full scalar curvature of g µν which depends on time derivatives, given bȳ The key role of condition C has been to yield equation (11), in which all the spatial derivatives present are those in the Laplace-Beltrami operator contained inÔ(t).
3 Construction of the method

Construction of the bases of modes
For each Cauchy hypersurface Σt we will construct a set of modes {Φ n (t)} fulfilling the following two properties: I. They will form, together with their complex conjugates {Φ n (t) * }, a complete orthonormal basis of the space of global solutions to the Klein-Gordon equation (1) with respect to the inner product (5). We stress that each mode of the basis is defined in the whole spacetime, the label [t] meaning only that we associate it to the corresponding hypersurface. Specifically, it is in this hypersurface that we will set its initial conditions. II. In the regions where ∂ t behaves like a Killing field around the hypersurface Σt [that is, h ij (t) remains constant for a long enough period around t =t], and where the function F (t) introduced in (10) can be made to vanish, the modes will be those with positive frequency with respect to t in that region (the corresponding complex conjugates will be those with negative frequency). 6 When this is not the case, the way we choose the modes can be thought just as an unambiguous recipe to associate an orthonormal basis of solutions to each hypersurface, even if no notion of positive frequency modes exists.  n (t)} to Cauchy hypersurfaces Σt.
In Fig. 1 we provide a graphical depiction of this association of bases of modes to Cauchy hypersurfaces, which shall be helpful when following the construction of the method. Since we impose that the modes are solutions to the Klein-Gordon equation, then because of condition B the only quantities left to fully determine them are the initial conditions; which, as mentioned before, we are going to fix at Σt. That is, we need to fix the quantities Φ  n (t)| t=t for each mode. The evaluations of the modes at t =t, which will be functions in Γt, are going to be given by the eigenvectors of the operatorÔ(t) in (10): Notice thatt in (14) is just a parameter: This is a partial differential equation in the spatial coordinates of the given Σt. 7 SinceÔ(t) is positive definite, we can define the real and positive quantities ω n ) 2 out of the eigenvalues of equation (14). We also impose the following orthogonality and normalisation condition to the solutions: Finally, the first time derivative of the modes evaluated at t =t will be given by Equations (14)(15)(16) fully determine the initial conditions, and therefore their time evolution determines the modes {Φ n ) * }. In Appendix C we prove that this 6 Regions in which F (t) cannot be made to vanish, but can be made arbitrarily small, may also be considered. We refer to Appendix B for further details on this condition. 7 As we commented previously for the operatorÔ(t), this equation should also be posed in the full Hilbert space L 2 (Σ t ). It is known, however, that the functions representing the eigenvectors can be taken to belong to the dense subspace Γt [38].

set of modes constitutes a complete orthonormal basis [in the Klein-Gordon inner product (5)] of the space of solutions to the Klein-Gordon equation (satisfying the possible boundary conditions).
We have then proven that the modes fulfil property I. Let us now consider property II. One can easily realise that, when ∂ t behaves like a Killing field in some spacetime region S around Σt, then q(t) =R(t) = 0 in S. Also, because of the time symmetry and the requisite that F (t) = 0 the eigenvalue problem (14) will have the same solutions for all the hypersurfaces Σ t within S. Taking into account these facts, it is easy to check that the modes of the form satisfy both the initial conditions (14)(15)(16) in Σt and the Klein-Gordon equation (11) in the whole region S. Therefore, these modes correspond to the modes that we are actually assigning to the hypersurface Σt. Now, if the interval of time that S embraces is large enough so as to explore the minimum frequency in the spectrum (now we can really call ω n a frequency), then it is clear that we can talk about the modes (17) as forming an orthonormal basis of modes with well-defined positive frequency with respect to t. Property II is therefore also fulfilled.

Time-dependent Bogoliubov transformation
Let us write down the Bogoliubov coefficients between any two orthonormal bases of modes, associated to the hypersurfaces Σt 0 and Σt: Since, according to condition B, all the modes are well defined globally, these coefficients exist and are unique as functions of time for all times. The objective of the method is to provide a differential equation in time that allows for their computation, without having to explicitly compute the time dependence of the modes. For convenience, we write the Bogoliubov transformation in the matrix notation as In Appendix D.1 we prove that this time-dependent Bogoliubov transformation satisfies the differential equation where we define the quantities With the initial condition U (t 0 ,t 0 ) = I, equation (21) has the formal solution where T denotes time ordering.
Equations (21) and (25) are the main result of this work: They switch from the time evolution of the modes to the time evolution of the transformation between bases. We also stress that one of the strengths of the method is that all the quantities appearing in (22)(23)(24), which are the coefficients of our differential equation, are known just by constructing the orthonormal basis of modes solutions to the eigenvalue equation (14), for which the timet is just a parameter.
We can get rid of the trivial phase in (21) and write the evolution in a more compact form. Apart from the phase introduced by the term iΩ(t), one can check that the diagonal elementsα nn (t) are purely imaginary [see relation (D.23) in Appendix D.2]. Thus, we can get rid of the full trivial phase by defining the diagonal matrix and writing the evolution in terms of a new operator Q(t,t 0 ) defined by Replacing (27) in (21), we get the differential equation which, again with the initial condition Q(t 0 ,t 0 ) = I, has the formal solution In Appendix D.2 we check that the Bogoliubov transformations (25) and (29) obtained with this method satisfy the Bogoliubov identities.

Physical interpretation
The time-dependent Bogoliubov transformation that the method provides contains all the information necessary to compute the evolution of the field in time. However, as we advanced in the Introduction, we do not pretend to give a quantisation for each and every basis of modes that we have constructed at each time. It is only in those regions S with time symmetry, that we described at the end of the Subsect. 3.1, where we can proceed to the usual Fock quantisation of the field; 8 that is, defining the corresponding Fock space with its vacuum state and creation and annihilation operators (ones the adjoints of the others in the case of a real field) associated to the mode decomposition given by the method, which in such region would be a decomposition in modes with well-defined frequency (for brevity, we do not expose this whole construction explicitly, see e.g. [29]). We rely on the fact that, under the existence of a timelike Killing vector field, Fock representation gives the correct physical description of a field in terms of particles perceived by the family of observers following the integral trajectories of the timelike Killing vector field. The "time-dependent Bogoliubov coefficients" that we constructed in the previous Subsection can be interpreted as true Bogoliubov coefficients relating the annihilation and creation operators of different mode decompositions only when they connect two regions in time where this valid quantisation procedure can be done. When this is not the case, within the scope of this article they should be considered just as an intermediate computational tool. Whether there is still some physical interpretation of the intermediate times between two static regions will be studied in future works. The Bogoliubov coefficients at those intermediate times may well be related to field observables (such as the stress-energy tensor) instead of particle observables, since those remain physically meaningful in the absence of timelike Killing fields.
A concern that can still be raised is that, even if two Fock quantisations done in two different regions S 1 and S 2 , in which ∂ t is a Killing field, are both valid and physically correct (which we know to be the case), they could eventually still be unitarily inequivalent. Unfortunately, we do not have a rigorous mathematical demonstration against such possibility occurring, but we can provide a solid physical argument for the cases considered with our method, taking into account that we only consider regular metrics. The condition for unitary inequivalence in the case of two quantisations related by a Bogoliubov transformation is that n,m |β nm | 2 → ∞. This implies a divergent number of particles in one of the quantisations for any state of the field with a finite number of particles in the other quantisation. If the energies of the particles in each quantisation were not limited from below, a divergent number of particles could happen without a divergent amount of energy (the socalled infrared catastrophe [39]). But thanks to condition A and the positive definiteness ofÔ(t) a minimum non-zero energy always exists. Therefore, for the problems considered with this method unitarily inequivalent quantisations would imply that, for states with finite energy in one of the quantisations, the energy diverges according to the other quantisation. However, for physically valid quantisations this possibility is also discarded by the regularity of the metric and the compactness of the spatial hypersurfaces, since these two facts imply that only finite amounts of energy can be interchanged with the field in finite times. Therefore, for the problems considered with this method, physically valid quantisations must be also unitarily equivalent.

Small perturbations and resonances
Let us consider now the case in which the spatial metric h ij (t) only changes in time by a small perturbation around some constant metric h 0 ij ; that is, where ε 1. Since, if boundaries exist, they remain parallel to ∂ t , all the spacelike hypersurfaces will be identical, and for convenience we will call them Σ 0 := Σt. We additionally require that F (t) remains O(ε), so that the solutions to the problem for ε = 0 are modes with well-defined frequency. Later on, we will check that such value of F (t) actually does not contribute at all to the physically relevant result of resonances, and thus the function F (t) can be fully ignored in the small perturbations regime. Staying to first order in ε, we can write the relevant quantities for computing the Bogoliubov coefficients as , where Φ 0 n and ω 0 n are the modes and frequencies solution to (14) with h 0 ij as the metric. It is immediate thatα nm (t) andβ nm (t) in (23) and (24) are O(ε), since the quantities dΦ where dV 0 is the volume element of the metric h 0 ij . Later on, we will further simplify these expressions for practical purposes.
Let us obtain the Bogoliubov coefficients to order ε in terms of these expressions. In the case of small perturbations, it is more convenient to consider the evolution as expressed by the Bogoliubov transformation Q(t f ,t 0 ) in (29). Noticing thatK(t) is O(ε), we only need to consider the zeroth order in the matrix Θ(t).
To first order in ε, the transformation reads where We can show the resonance behaviour of the field in a clear way if we write explicitly the expressions for the Bogoliubov coefficients: One can see that, in general, the Bogoliubov transformation will differ from the identity just by terms of first order in ε, except for the cases where there are resonances. That is, if the perturbation of the metric contains some characteristic frequency ω p , then the same frequency will usually be also present in the quantities ∆α nm (t) and ∆β nm (t); and if this frequency coincides with some difference between the frequencies of two modes, ω p = ω 0 n − ω 0 m (it is in resonance), then the corresponding coefficient α nm (t f ,t 0 ) will grow linearly with the time differencet f −t 0 , and after enough time it will overcome the O(ε). Respectively, if the characteristic frequency coincides with some sum between the frequencies of two modes, ω p = ω 0 n + ω 0 m , then the corresponding coefficient β nm (t f ,t 0 ) will grow linearly in time and eventually overcome the O(ε). For completeness, in Appendix E.1 we show that the resonances remain stable under small deviations of the frequency of the perturbation from the exact resonant frequency.
If the Fourier transform F of ∆α nm (t) [respectively ∆β nm (t)] exists as a well-defined function, which necessarily implies that the perturbation vanishes fast enough in the asymptotic past and future, another way to consider the resonances is by taking the limitst 0 → −∞ andt f → ∞ in (36) and (37) and writing That is, the Bogoliubov coefficients between the asymptotic past and future are proportional to the Fourier transforms evaluated at the corresponding subtraction (respectively addition) of frequencies. Evidently, resonances will occur if the frequency spectrum is peaked around one or more of these values.
It is important to remark that, since there is no strict time symmetry [neither F (t) vanishes in general], strictly speaking there would be no notion of particles associated to modes with well-defined frequency, except if the perturbation really vanishes in certain regions. However, even while the perturbation is ongoing, the deviation from the time symmetry is of order ε. Thus, in the resonance regime, when one or more coefficients grow beyond this order, this lack of symmetry can be neglected with a solid physical argument: The deviation from time symmetry is small, and in particular much smaller than the effects measured. Therefore, in the case of resonance one can safely say that the physical effects are taking place for particles associated to modes with well-defined frequency.
In Appendix E.2, we further support the interpretation in the previous paragraph by proving that the results on resonances do not depend on the particular choice of the basis of modes to order ε that we have used to compute the evolution. In order to give such proof, in Appendix E.2 we also prove a result that is useful at this point to further simplify the expressions in (32) and (33): Any contribution of the form respectively where X(t) can be any function of time, appearing in the expression of ∆α nm (t), respectively ∆β nm (t), will not affect the resonance behaviour of the corresponding modes. In simple words, this happens because, if X(t) contains the correct resonant frequency in its Fourier expansion, the contribution of the corresponding term of the expansion to (40) [resp. (41)] vanishes. We refer to Appendix E.2 for the details. Taking into account this fact, the expressions in (32) and (33) can be simplified to where the symbol '≡' denotes the equivalence relation "gives the same resonances as". Since resonances are the only physically meaningful result to be obtained from this computation, given a concrete problem one can always use these simplified expressions to compute the sensibility of the field to each resonance, by identifying the corresponding term of the Fourier expansion or value of the Fourier transform.
It is useful to notice that, although in principle we have assumed the orthonormality of the modes Φ n (t), for using the expressions in (42) and (43) only the nonperturbed modes Φ 0 n need to be strictly orthogonal [with respect to the usual inner product in L 2 (Σ 0 )] and with the correct normalisation given by (15) to zeroth order. The modes Φ n (t) need only to satisfy the unavoidable orthogonality already guaranteed by being correct solutions (to order ε) of the eigenvalue problem (14); that is, they will necessarily be orthogonal (to order ε) unless they have the same eigenvalue (ω It is easy to check that any change of order ε in the normalisation of the modes (or in the inner product between modes with the same eigenvalue) would not contribute to (42) or (43), because the corresponding factor (ω 0 n ) 2 − (ω 0 m ) 2 would make it to vanish. Therefore, when solving a specific problem, one only needs to explicitly ensure orthogonality and properly normalise the modes to zeroth order.
Finally, we give other alternative simplified formulas for the computation of ∆α nm (t) and ∆β nm (t), which justification is given in Appendix E.3. These formulas have the great advantage that they only imply the modes Φ 0 n and eigenvalues ω 0 n of the static problem: where∆(t) is a linear operator defined by its action on the basis {Φ 0 n } aŝ with ε∆Ô(t) being the first order term in ε of the oper-atorÔ(t). For problems in which it is easy to compute the eigenvalues and eigenfunctions to first order in ε, it might be easier to apply (42) and (43); while for other problems the expressions in (44)(45)(46) are definitely more convenient. Finally, as we already advanced, we can check in these last expressions that the quantity ∆F (t) never contributes to the resonances, since in (46) its direct contribution to∆(t) cancels with its contribution through ∆Ô(t) [recall the definition in (10)]. The resonance can be consistently described in the regime of duration of the perturbation ∆t such that 1 ω p ∆t 1/ε; since one needs the period of time to be reasonably larger than the inverse of the frequency being described, but on the other hand, one should keep the second order term in ε that we dropped in (34) significantly smaller than the first order term that we kept, so that the first order expansion of the exponential in (29) remains valid. Notice that, because of these limitations, the net effect after integrating in time under resonance might be of order greater than ε, but still should be small as compared to unity. However, even the effect of such small Bogoliubov transformation could in principle be measured using the adequate probe states as the initial states of the field and quantum metrology techniques (see [17,21,[40][41][42]). Thus, even in these situations, in which the total energy involved was too small to produce a pair of particles or to promote an existing particle to a higher energy state, this would not mean that there are no measurable effects at all. Moreover, in this work we restrict ourselves to an idealised system in which the evolution of the field is perfectly unitary. More realistic systems may also present decoherence any time quanta is created or promoted, something that can then lead to cumulative effects.

A concrete example
In order to show the method in action, we will consider an example of application within the perturbative regime. In particular, we consider a field trapped in a cavity which is perturbed by a gravitational wave. This is a very interesting problem in the growing research field of gravitational wave detection with Bose-Einstein condensates. In [21] the authors solve a similar problem by using the restricted version of the method in [20]. They introduce a simplified toy model in one spatial dimension, taking advantage then of the conformal invariance in order to "translate" the problem into an equivalent one in Minkowski spacetime. With the general method presented here, it is possible to go beyond that scheme and solve the problem in three dimensions. The three-dimensional problem was also considered in [42] with a very different approach. We will reproduce and extend the result in that work with very simple and straightforward calculations. As in [21] and [42], we consider our field to be the phonon field of a trapped Bose-Einstein condensate, therefore a real scalar massless quantum field. Following the work in [43], one can see that, in the case that the condensate remains stationary, the phonon field experiences an effective metric (with minimal coupling) which is the gravitational wave metric with the speed of light appearing in the component g 00 replaced by the speed of sound in the condensate c s . This speed of sound is the one that we normalise for this problem, c s = 1. We work in the TT-gauge and consider a wave with amplitude ε and frequency Ω propagating in the z-direction and with polarisation in the xy-directions. This yields the metric where we have simplified sin[Ω(t − z/c)] → sin(Ωt), as we have that c ΩL z (being L z the size of the condensate in the z-direction), because of the orders of magnitude between the speed of light and the speed of sound. For simplicity, we consider that the condensate is trapped in a rectangular prism of lengths L x , L y and L z aligned with the directions of propagation and polarisation of the wave. The boundaries are free-falling, and we impose Dirichlet vanishing boundary conditions.
For the metric in (47) and for m = ξ = F (t) = 0 the eigenvalue equation (14) where it is clear that we will need three quantum numbers. Since the boundaries are free-falling, the boundary conditions read Φ(x = 0) = Φ(x = L x ) = 0, and equivalently for the other dimensions. The solutions to order ε, with the normalisation (15) to zeroth order, are where k x n := πn/L x , and equivalently for the other dimensions; with the eigenvalues to order ε given by Notice that, for this particular problem, the free-falling boundaries make the eigenvalues to depend on time while the modes are time-independent to order ε. We 9 When constructing the integration method, we have been using the different notations t andt for time, in order to distinguish between the evolution of modes in time and the choice of a concrete basis of modes associated to a Cauchy hypersurface Σt. This was indeed necessary for constructing the method. But for its application to a concrete problem we will not need to explicitly consider the evolution in time t of a basis of modes anymore. Thus, we can relax the notation and replacet → t.
plug the solutions into (42) and (43), together with ∆q(t) = 0 computed from (12), obtaining (n , m , ) = (n, m, ); By looking at these expressions, and considering the resonance behaviour found in (37), it is straightforward to obtain the conclusions: When Ω = 2ω 0 nm , the corresponding diagonal β-coefficient will grow linearly in time due to resonance as while the rest of the coefficients will remain trivial in any case.
In [42], the authors consider the more realistic case in which the gravitational wave lasts some approximate time τ and vanishes asymptotically, by replacing sin(Ωt) → e −t 2 /τ 2 sin(Ωt) in the metric (47). We can easily handle this perturbation with our method by using the expression with the Fourier transform (39), obtaining (for the only non trivial coefficients) Equation (54) exactly reproduces, through a very simple calculation (just a trivial eigenvalue problem and a trivial Fourier transform), the relation found in [42] (equation 2.21). When Ω 2ω 0 nm , it reproduces the resonance behaviour for the range of durations 1/Ω τ 1/|Ω − 2ω 0 nm |, as discussed in Appendix E.1. We can obtain some further conclusions on this simple model. For example, it is straightforward to check that imposing Neumann vanishing boundary conditions leads to the same results. One would need to temporarily introduce a small rest mass term ∆m in the field equation, in order to dodge the presence of a zerofrequency mode, as described in Appendix B; but it is easy to check that taking the limit ∆m → 0 after solving the problem is in this case trivially well defined. Also, we can discuss what happens if we include the cross-polarisation in the metric. This polarisation would add a term proportional to ∂ x ∂ y in the operator in (48). By using the expressions (44)(45)(46), it is easy to check that this does not contribute to any Bogoliubov coefficient. Therefore, for waves coming perpendicular to any of the faces, the field is sensitive only to the polarisation which is aligned with the other faces. Again by using those expressions, it would not be difficult to even consider a wave propagating in an arbitrary direction, by applying the corresponding rotation matrices with the Euler angles to the metric, but for brevity we will not approach this problem here. In future publications [44], the authors will apply the method introduced in this article to study in detail the resonance properties of this and other models of trapping cavities, in what will be a fruitful application of the method to the study of gravitational wave detection with Bose-Einstein condensates.

Application to cosmological particle creation
In this section, we give an example of the application of the method beyond the perturbative regime. In particular, we use it to analyse a massive, minimally coupled Klein-Gordon field in a one-dimensional FLRW flat spacetime. The line element of this spacetime is given by For the metric in (55), where h xx (t) = a(t) 2 , and for ξ = F (t) = 0, the eigenvalue equation (14) takes the form 10 As it is customary [2,29], we consider the spatial dimension to be a 1-torus, for which purpose we introduce periodic boundary conditions in space, so that Φ(x = 0) = Φ(x = L), with L being the length of the torus. The solutions to (56), with the normalisation in (15), are given by n |a(t)|L e iknx , n ∈ Z; where ω [t] n = k 2 n a(t) 2 + m 2 and k n := 2πn L . From (12) we get q(t) =ȧ(t)/a(t). Computing (23) and (24) with the field modes (57) we obtain α nm (t) = 0;β nm (t) = 0, m = −n; Observe that if the scale factor is constant, then we obtainβ n(−n) (t) = 0 and the transformation is trivial. The transformation is also trivial when the scalar field is massless, which agrees with the fact that particles are created in (1 + 1)-dimensions only if a non-zero mass breaks the conformal invariance [29].
With the result obtained in (58), the differential equation (28) written directly for the Bogoliubov coefficients readṡ By replacing n → −n in the second equation, we realise that the system of equations can be decoupled in independent systems of two differential equations for each pair {α nm , β (−n)m }. Moreover, since the equations are homogeneous and the initial conditions are α nm (t 0 , t 0 ) = δ nm and β nm (t 0 , t 0 ) = 0, only the systems for which m = n will have a non-trivial evolution, while the rest of the coefficients will all stay equal to zero. For the non-trivial systems, we can write down the two first order differential equations more conveniently as the following second order differential equation for α nn (t, t 0 ): (60) with the initial conditions given by α nn (t 0 , t 0 ) = 1 andα nn (t 0 , t 0 ) ∝ β (−n)n (t 0 , t 0 ) * = 0. Once computed the α nn (t, t 0 ) coefficients, computing another integral would give the β (−n)n (t, t 0 ) coefficients. However, it is much easier to compute them (up to an unimportant phase) out of the Bogoliubov identity m (|α nm | 2 − |β nm | 2 ) = 1, which in this case simplifies to The mathematical results up to here are completely general, there is no assumption made on the form of the scale factor. As it was already discussed, the Bogoliubov coefficients computed will only have a clear physical interpretation when they connect regions in which ∂ t is a Killing field. We already reproduce the known result [2,29] that in the problem considered there is never mode-mixing (the only non-zero α's are the diagonal ones), but only particle creation in pairs of fully entangled particles with the same frequency and equal but opposite momentum. Finally, we point out that we could easily take the limit of an unconfined quantum field (L → ∞) by simply replacing the role of the quantum number n by the wave number now raised to a continuum quantity, k n → k ∈ R; and the discrete frequencies by the corresponding dispersion relation, ω [t] n → ω [t] (k).

A concrete example
To illustrate the previous general result, let us use it to compute numerically the Bogoliubov coefficients for a well-known toy model of cosmological expansion; in particular, the toy model appearing in [29], pp. 59-62, and which was first proposed in [45]. In this model, the expansion is explicitly written in terms of the conformal time η, defined by dη/dt = 1/a(t). In this time, the scale factor is given by where A, B and ρ are constants. One can see that, in the asymptotic past and future (η, t → ±∞), the spacetime approaches the Minkowskian flat spacetime. In [29] the Bogoliubov coefficients between the asymptotic regions α nn (∞, −∞) and β (−n)n (∞, −∞) are obtained analytically by first solving the Klein-Gordon equation in conformal time through separation of variables.
Unfortunately, for the given expression of a(η) we cannot obtain an explicit analytic expression for a(t) (one could expect this to happen, since the toy model is prepared ad hoc to be solved in conformal time). Nonetheless, we can solve this problem numerically using the differential equation (60) that we have obtained. We first compute the functions appearing in the coefficients of the equation, and finally solve the differential equation itself. The asymptotic initial conditions are α nn (−∞, −∞) = 1 andα nn (−∞, −∞) = 0. Being trivially related by (61), it is equivalent to follow the evolution of the α's or the β's. For convenience, we consider the second ones, since they directly give the interesting result of particle creation for the quantum field. In Fig. 2, we plot a graphic with an example of numerical results for |β (−n)n (t, −∞)| 2 for a concrete choice of the parameters in the problem, comparing them to the asymptotic values obtained analytically in [29].
One can clearly see that the asymptotic values in [29] are correctly reproduced. We also claim the attention about the fact that, even if the intermediate values of the coefficients cannot be assigned an immediate interpretation in terms of particle creation, the evolution of the coefficients in time is perfectly smooth and monotonically increasing, with a time dependence qualitatively similar to that of the scale factor itself. This "well-behaved" result is remarkable, and makes us to suspect that those intermediate values may have some physical meaning in terms of some sort of approximation, although the discussion of this idea lies beyond the scope of this article.

Summary and conclusions
We have developed a method for computing the evolution of a confined quantum scalar field under general changes of the spacetime metric, under the assumptions given in Subsect. 2.1. In order to keep track of the evolution of the field, instead of solving the Klein-Gordon differential equation, the method computes a time-dependent Bogoliubov transformation between the bases of modes associated to different spacelike hypersurfaces. It provides the way to obtain the initial conditions of the basis of modes associated to each spacelike hypersurface as an eigenvalue problem in that hypersurface [equation (14)]. Once obtained these initial conditions, the method provides a system of homogeneous first order differential equations in time [equation (21) or (28)] for the Bogoliubov transformation between the bases of modes associated to different times. The coefficients of such system of equations are computed out of the initial conditions of the modes themselves [equations (22)(23)(24)]. Finally, we can also arrive to the formal generic solution of the system of equations [equation (25) or (29)]. In this article we have considered the geometries without boundaries and those for which the boundaries remain static in some synchronous gauge, leaving the case of "moving boundaries" for a second article.
The time-dependent Bogoliubov transformation obtained with this method cannot, in general, be given a direct physical interpretation in terms of mode-mixing and particle creation phenomena, simply because the modes constructed for each spacelike hypersurface do not correspond (in general) to modes with well-defined frequency with respect to some time, and therefore they do not yield a physically unambiguous quantisation in terms of particles. However, if there is a region of spacetime (long enough in time) in which the time translation along the chosen time coordinate is a symmetry, 11 then automatically the modes as constructed with the method correspond to modes with well-defined frequency in that region. In such cases, the quantisation in these modes has the adequate interpretation in terms of particles, and the Bogoliubov transformation between two such regions has also its usual physical meaning in terms of mode-mixing and particle creation. In the regions where this is not the case, the modes and Bogoliubov transformations constructed can be considered just as an "integration tool" for computing the physically meaningful cases, although their potential meaning as some sort of approximation will be analysed in future works.
We stress that the method is of general applicability, just under the considered assumptions. However, as much as we stress this generality, we should also emphasise that the method proves to be of special practical use in the case of small perturbations of the metric (or the boundary conditions, within the next forthcoming article on the method) around a static situation. For clarity, we number here the reasons of this utility: -The perturbations of the modes can be obtained out of the eigenvalue equation (14) just with perturbation theory to first order. In the same way, for computing the quantities ∆α(t) and ∆β(t) one should only keep the first order in the perturbation, taking advantage of the equivalence relation with respect to resonances described in Appendix E.2 to further simplify the expressions [equations (42) and (43)]. Moreover, one can alternatively use the expressions (44)(45)(46) to obtain the results directly from the solutions to the corresponding static problem. -There is no need to explicitly solve the differential equation (28) or to compute the exact transformation (29). In (36) and (37) we can already see that, to first order, there is just a resonance behaviour. Thus, writing down the Fourier expansion of the perturbations ∆α(t) and ∆β(t) automatically gives which Bogoliubov coefficients will stay to order ε (the order of the perturbation) and which (if any) will grow linearly in time due to resonance, and at which rate. Alternatively, when it is possible one can also take the Fourier transform to find the Bogoliubov coefficients between the asymptotic regions using (38) and (39).
-Since the problem consist of a small perturbation of order ε around a static situation, in which there would be an exact notion of particles with welldefined frequency, the indefiniteness in such a notion will stay to order ε. Thus, for the modes which are in resonance (the interesting ones), for which the Bogoliubov coefficients will grow beyond the order ε, the physical interpretation in terms of mode-mixing and particle creation is perfectly valid. -In Appendix E we prove that the result on resonances is stable under small deviations from the resonant frequency. We also verify that the resonances obtained are the same regardless of which bases of modes to order ε one uses to integrate the evolution, which further supports that they are the physically meaningful result.
As an example of the application of the method, we have considered the perturbation of a confined field by gravitational waves in a realistic 3 + 1 dimensional background. Through very simple calculations, we obtain the Bogoliubov coefficients, which agree with those given in [42] for the same problem, and further discuss the role of the boundary conditions and the polarisation. The perturbative method could also prove its utility in other problems which are now under study, in which quantum systems are perturbed by small gravitational effects [16,46,47].
Finally, we have included an example of an application of the method to cosmological particle creation. We have obtained a differential equation that allows to compute the Bogoliubov coefficients for any onedimensional FLRW flat spacetime, and applied this differential equation to an already known toy model of cosmological expansion, reproducing the correct results for the Bogoliubov coefficients between the asymptotic regions. In this way, we checked that the method can also be useful in problems beyond the perturbative regime. Out of the results obtained for the toy model, we highlight the smooth and monotonic behaviour in time of the Bogoliubov coefficients computed, even when they do not necessarily have an immediate physical interpretation. As we already mentioned, we leave for future work the analysis of the possible meaning of those values as some sort of approximation. We also leave for future work potential new applications of the integration method introduced in this article to other cosmological scenarios, as well as further applications in other problems in quantum field theory in curved spacetime. In this Appendix we provide some further discussion on the conditions A, B and C introduced in Subsect. 2.1.
Condition A reflects the fact that the method applies only for confined fields. It plays its role in ensuring that the operatorÔ(t) in (10) has a discrete spectrum, something crucial for the matrix formalism that we employ. In some situations (such as the example we give in Sect. 5), it might be possible to obtain results for unconfined fields by taking the limit of some typical length of the confinement region going to infinity after solving the problem [2,29]. However, in most of the situations this limit may not be well defined, or may not commute with a limit previously taken when solving the problem, such as a late time limit.
Condition B is introduced to ensure that the Cauchy problem for the scalar field is well posed, so that we can relate one-to-one initial conditions and solutions. Notice that this condition amounts to the predictability of the behaviour of the field in the given geometry. The aim of the method presented in this article is not to prove such predictability, but rather to compute the evolution of the field, for which it is reasonable to assume the predictability as given. If M has no boundaries, then condition B is automatically ensured by the well-known theorems on existence and uniqueness of global solutions to partial differential equations. Unfortunately, there is to our knowledge no extension of these theorems in the presence of timelike boundaries. For Dirichlet boundary conditions, a proof of local well-posedness of the Cauchy problem for the Klein-Gordon equation is given in [48] (where the authors consider the uniform Lopatinski boundary condition, which includes Dirichlet as a specific case). Global wellposedness has been proven for the Dirac equation under MIT-boundary conditions, which are the analogous to Neumann boundary conditions for a scalar field [49]. In the case of Robin boundary conditions, global wellposedness is proven in [50] for the Klein-Gordon equation in static spacetimes of bounded geometry (for a wider class of boundary conditions including the Robin case). It is probable that global well-posedness of the Cauchy problem for the Klein-Gordon equation can be proven for the boundary conditions considered in this article as an extension of the results cited, but approaching such task is not in the aim of the present work, which goals are mainly computational. Therefore, lacking (so far) definitive theorems, we prefer to explicitly impose global well-posedness as a requisite. Nonetheless, we strongly believe that it will be satisfied, in the worst case, except for too stilted geometries, most probably without any physical interest.
Condition C, which may seem unmotivated, is crucial for being able to write the Klein-Gordon equation as in (11), with all the spatial derivatives collected in the self-adjoint operatorÔ(t). Without this, the construction of the bases of modes and the time-dependent Bogoliubov transformation would not be possible, at least in the form we provide it. The existence of a gauge in which the crossed time-space terms in the metric cancel is guaranteed by the global hyperbolicity, with or without boundaries [37]. In most of the cases of physical interest, it should be also possible to make the lapse function −g tt equal to one globally. In any case, this is always possible locally, so even if no global synchronous gauge exists, one may still apply the method "by segments". However, we cannot provide a general recipe on how to approach the joint of the results of the different segments.
Appendix B: Function F (t), exponential time dependence and zero-frequency modes The quantity F (t) has been introduced in the opera-torÔ(t) in (10) for a purely technical reason; namely, in order to have a positive-definite operator, so that all the eigenvalues obtained from the equation (14) are positive, something necessary for the construction of bases of modes and the computation of Bogoliubov coefficients to hold as it has been exposed. Since, while there is no timelike Killing field, the construction of bases and Bogoliubov transformations is purely instrumental, there is no need to further discuss the role of F (t) in such case. However, when stating the properties of the bases constructed in Subsect. 3.1, in property II, together with ∂ t being a timelike Killing field in a region of spacetime, we require also that F (t) = 0 in such region so that the modes constructed are modes with well-defined frequency (and we prove that this is the case). We can better understand this requisite, and in general the role of F (t), if we analyse what happens when ∂ t is a timelike Killing field but F (t) > 0. There are two possibilities: -F (t) must take some finite positive value, so that the operatorÔ(t) is positive definite as required.
If we insisted in setting F (t) = 0, we would obtain an operator with some negative eigenvalues. This would imply the existence of imaginary values for the quantities ω n obtained from the eigenvalue equation (14). It is not difficult to check that in such case there exist fundamental solutions of the Klein-Gordon equation that have an exponential dependence on t, by considering equation (17) with imaginary ω [t] n . 12 Under this circumstance, it is not possible to proceed to the Fock quantisation of the field, even under the existence of a timelike Killing field. However, the method can still keep track of the evolution of the field in that region, if not for quantising it. By introducing the positive function F (t), we build bases of solutions associated to each hypersurface, although none of these solutions correspond to the fundamental solutions with exponential behaviour. They are just bases of solutions as instrumental as in the cases in which there is no timelike Killing field. But, precisely because they are not fundamental solutions, the Bogoliubov transformations between these different basis, obtained from the dynamics of the field, are not trivial, since F (t) is present in (23) and (24). In other words, the exponential behaviour that is "artificially dropped" from the modes by the introduction of F (t) is picked back and reproduced by the Bogoliubov transformation between them.
-F (t) could be made arbitrarily small, but not zero, in order for the operatorÔ(t) to be positive definite. In this case, setting F (t) = 0 would yield a positive semi-definite operator, which implies the existence of a fundamental zero-frequency mode. If one is not interested in the quantisation of the field in the region, one can proceed by giving F (t) some fixed value and compute the corresponding Bogoliubov coefficients. However, if one wishes to attempt quantisation in such region, one could instead add a small increment to the rest mass ∆m, so that the operatorÔ(t) is positive definite even with F (t) = 0. Yet, after solving the problem one should take the limit ∆m → 0 in order to recover the original field. This limit might not be trivial, and the method cannot give an a priori prescription on how to take it. Even so, in Subsect. 4.1 we give an example where this procedure proves to be successful in a trivial way. Notice however, that this difficulty with a zerofrequency mode is intrinsic to the field properties and not peculiar to the method.
We can identify some circumstances in which one can always set F (t) = 0 and forget about the quantity. In particular, let us consider the important case of a minimally coupled field (ξ = 0). In such case, nonzero rest mass and/or Dirichlet boundary conditions or Robin boundary conditions with γ > 0 [52] automatically make the operatorÔ(t) to be positive definite; while for massless fields without boundaries or with Neumann boundary conditions the operator is positive semidefinite and could be handled with the procedure described in the second point above.

Appendix C: Completeness and orthonormality of the bases of modes
Let us prove that the set {Φ n ) * } constructed in Subsect. 3.1 is a complete basis of the space of solutions to the Klein-Gordon equation which satisfy the possible boundary conditions. Because of condition B and the linearity of the Klein-Gordon equation, it is equivalent to prove that the set formed by the initial conditions for each mode in that set at Σt is itself a basis of the space of initial conditions at Σt. If we recall (16), this set of initial conditions of the modes is given by the following set of pairs of values of the field and of its first time derivative at Σt, (Φ(t), ∂ t Φ(t)| t=t ): It is clear that the space of initial conditions at Σt [pairs of functions (Φ(t), ∂ t Φ(t)| t=t )] is L 2 (Σt)⊕L 2 (Σt). We will prove in short that the set in (C.1) is indeed a basis of such space. However, we also need to ensure that the elements in the basis satisfy the possible boundary conditions at ∂Σt compatible with the boundary conditions in ∂M , when they exist. For the first component of the elements in (C.1) (the evaluation of the mode at Σt), this compatibility means being in Γt, since this is the definition of such subspace [together with the smoothness condition that simply picks up a representative function of the vector in the Hilbert space L 2 (Σt)]. Such condition is clearly satisfied. Let us consider the second component of the elements in (C.1) (the time derivative of the mode at Σt). If we take the total time derivative of any of the boundary conditions in (2)(3)(4) along the boundary, since the boundary itself remains static, we will easily find that the same boundary condition that applies to Φ(t) must apply also to ∂ t Φ(t). Therefore, we have that the second component of the elements in (C.1) must also belong to Γt, something which is also clearly satisfied.
Observe that the fact that both the evaluation of the mode and of its first time derivative at some Cauchy hypersurface Σt satisfy the same boundary conditions is what allows for the construction of the modes as it has been done in this article [that is, imposing equation (16)]. This construction was therefore feasible thanks to the condition that the boundaries shall remain static in the synchronous gauge considered. If this is not the case and the boundaries are not parallel to the time flow ∂ t , the boundary conditions at ∂M as given in (2-4) (or their derivatives along the boundary itself) mix spatial and time derivatives of the mode in the same equation, and the construction of initial conditions as presented in this article cannot be done any more. This is the situation that we will consider in the second forthcoming article on the method. 13 Finally, let us prove that the set (C.1) is a basis of L 2 (Σt) ⊕ L 2 (Σt). For convenience, we temporarily simplify the notation for the functions and eigenvalues by replacing Φ n → ω n . We notice first that clearly both sets of eigenvectors {Φ n } and {Φ * n } are complete bases of L 2 (Σt). 14 This means that the modes of the second basis can be expanded in the first basis with unique coefficients: since the subspaces generated by eigenfunctions of the same operator with different eigenvalues are orthogonal, and p nm is clearly invertible. Also, any two arbitrary functions Φ A and Φ B in L 2 (Σt) can be expanded in the basis {Φ n } with unique coefficients: Using the expansions in (C.2) and (C.4) we will prove that the generic pair (Φ A , Φ B ) ∈ L 2 (Σt)⊕L 2 (Σt), that is, a generic initial condition, can be expanded in the set of pairs in (C.1) with unique coefficients, which proves that such set is a basis of the space of initial conditions. If we write Multiplying (C.6) by ω m and subtracting both equations, we arrive at where for obtaining the last equation we have used (C.3). Since p nm is invertible, the coefficients d n exist and are unique. Using (C.6) we also obtain unique coefficients c n . This completes the proof. With respect to the orthonormal character of the bases, using the initial conditions and the normalisation given in (15) one can easily compute the Klein-Gordon inner product (5) between any two modes of the basis, obtaining: where in (C.10) we have used that the integral is nonzero only when the eigenvalues coincide. Therefore, we have proven the completeness and orthonormality of the bases of modes that we are constructing for each hypersurface.
Because of relations (16) and (D.16) (replacingt → t + ∆t), we have that This is the local evolution of the needed quantities to first order in ∆t, and therefore we are ready to compute the derivatives in (D.14) and (D.15). If we plug (D.17) and (D.18) into (D.14) and (D.15), we get, with some manipulation, where the quantities are those defined in (23) and (24). This completes the proof of equation (21).

Appendix D.2: Bogoliubov identity checking
Any Bogoliubov transformation B [written in matrix notation as in (20)] must satisfy the following identity: The Bogoliubov transformations that we have found in (25) and (29) both have the form B = e D (the temporal ordering only affects the limits of the integrals and clearly does not play any role). In such a case, noticing that J is invertible, the condition (D.21) implies the following condition for D: Thus, noticing that Ω(t), Θ(t) and A(t) are diagonal, it is clear that both U (t,t 0 ) and Q(t,t 0 ) will satisfy the Bogoliubov identity (D. 21) if and only if identity (D. 22) is satisfied for D = K(t). This condition is equivalent to the relationŝ α nm (t) = −α mn (t) * ,β nm (t) =β mn (t). (D.23) Let us prove these relations by taking the derivative with respect tot of the inner products given in (C.9) and (C.10). Since these inner products are constant, with some manipulation we have that also need thatt f −t 0 1/|ω 0 n − ω 0 m | for the resonance to take place, this means that, as far as the frequency is sufficiently close to the resonant one, there is a wide window of time between 1/|ω 0 n − ω 0 m | and 1/|∆ω| in which resonance occurs. Analogous arguments apply to the resonance for the β-coefficients.
Appendix E.2: Invariance under small changes in the basis of modes In the perturbative regime, the modes with physical interpretation, those that can be used for quantisation, are the resonance modes of the static configuration; that is, the static solutions Φ 0 n . The only physically meaningful effect that can be described are the resonances between them, since anything else remains as tiny as the perturbation, and therefore as tiny as the degree of arbitrariness in the choice of the solutions to order ε that we use to keep track of the evolution. Here, we prove that the results on resonances do not depend on the concrete perturbed bases used to integrate the evolution of the field (as far as they are bases of solutions correctly computed to order ε).
Consider that, instead of the bases {Φ    n (t)} which differ just by order ε of the original ones (notice that we do not even require that "half" of the modes are the complex conjugate of the others). Since both are complete bases of solutions, they have to be related by      1 χ  Consider now the differential equation (21). Taking into account that K(t) has no zeroth order terms in ε, the corresponding differential equation for V (t,t 0 ) correct to order ε will be d dt V (t,t 0 ) = iΩ(t) + K(t) with Ω 0 := Ω(t)| ε=0 = diag(ω 0 1 , ω 0 2 , . . . , −ω 0 1 , −ω 0 2 , . . .). One can see that, in principle, there is a non-zero contribution to the differential equation due to the change of bases. Let us now write down the matrix T (t) in the block form This is just for notational purposes, there is no need that the coefficients satisfy any relation. We consider for example the coefficient X nm (t), which would contribute to the quantity ∆α nm (t) in the new transformation with the amount which remains O(ε) at any time, and therefore does not contribute to the resonances. If we had considered the Bogoliubov coefficient as given by the Fourier transform in (38), the conclusion is even more obvious, since the Fourier transform of (E.31) identically vanishes when evaluated at ω 0 n − ω 0 m . Analogous arguments apply to the other entries of the transformation.
Notice that, if the transformation I + εT (t) is not a Bogoliubov transformation to order ε, then the Bogoliubov coefficients obtained will not fulfil the Bogoliubov identities to order ε. This happens because the bases of perturbed solutions used in this case are not orthonormal to order ε. However, this is completely irrelevant, since these bases are just an intermediate mathematical tool to correctly compute the evolution. The modes in the bases need to be correct solutions to order ε and coincide with the static solutions to zeroth order, nothing else. The final Bogoliubov transformation should not be interpreted physically to order ε, but only beyond it; that is, only for the resonances, which as we have proven remain the same, no matter which transformation (to order ε) is done to the bases.
It is straightforward to connect the computation done in this Appendix to the claim done in Sect. 4, which states that the quantities of the form in (40) appearing in ∆α nm (t) [respectively of the form in (41) in ∆β nm (t)] will not contribute to the resonances: As one can see from (E.31), adding or subtracting these quantities is equivalent to an irrelevant change in the choice of bases to order ε.

Appendix E.3: Further simplification of the expressions
Let us obtain expressions (44)(45)(46) from (38) and (39). First, we write the operatorÔ(t) to first order in ε aŝ O(t) ≈Ô 0 + ε∆Ô(t). (E.34) Zeroth and first orders in ε of equation (14) then read, respectively,Ô Using these two expressions and the fact thatÔ 0 is selfadjoint, we can write Finally, we notice that, because of the equivalence relation with respect to resonances, any diagonal term in ∆α nm (t) is always meaningless. We can thus ignore the second term in the l.h.s., obtaining (44) and (46) for ∆α nm (t). An identical derivation with the replacement (Φ 0 m ) * → Φ 0 m yields (45) for ∆β nm (t), in which case the diagonal terms are relevant.
Appendix F: Summary of formulae for the application of the method In this Appendix, we provide Tables 1 and 2 with all the formulae necessary in order to apply the method to a concrete problem. As we indicated in the examples in Subsect. 4.1 and Sect. 5, once the expressions for the method have been found, there is no need to consider the explicit evolution in time of the modes constructed any more. Therefore, the notation for the "time label" of the different modes can be simplified replacingt → t. In the tables we use this simplification.