# Four-mode squeezed states: two-field quantum systems and the symplectic group $Sp(4, R)$

Thomas Colas, Julien Grain, Vincent Vennin

► To cite this version:

Thomas Colas, Julien Grain, Vincent Vennin. Four-mode squeezed states: two-field quantum systems and the symplectic group $Sp(4, R)$. Eur.Phys.J.C, 2022, 82 (1), pp.6. 10.1140/epjc/s10052-021-09922-y. hal-03229409

HAL Id: hal-03229409

https://hal.science/hal-03229409

Submitted on 20 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Four-mode squeezed states: two-field quantum systems and the symplectic group Sp(4, $\mathbb{R}$)

Thomas Colas<fnref n="1" />,<fnref n="2" />,ᵃ ID, Julien Grain<fnref n="1" />,ᵇ, Vincent Vennin<fnref n="2" />,ᶜ

Received: 29 October 2021 / Accepted: 8 December 2021 / Published online: 5 January 2022

© The Author(s) 2022

**Abstract** We construct the four-mode squeezed states and study their physical properties. These states describe two linearly-coupled quantum scalar fields, which makes them physically relevant in various contexts such as cosmology. They are shown to generalise the usual two-mode squeezed states of single-field systems, with additional transfers of quanta between the fields. To build them in the Fock space, we use the symplectic structure of the phase space. For this reason, we first present a pedagogical analysis of the symplectic group Sp(4, $\mathbb{R}$) and its Lie algebra, from which we construct the four-mode squeezed states and discuss their structure. We also study the reduced single-field system obtained by tracing out one of the two fields. This procedure being easier in the phase space, it motivates the use of the Wigner function which we introduce as an alternative description of the state. It allows us to discuss environmental effects in the case of linear interactions. In particular, we find that there is always a range of interaction coupling for which decoherence occurs without substantially affecting the power spectra (hence the observables) of the system.

## Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Quantum phase space of two free fields</td><td>2</td></tr><tr><td>3</td><td>Sp(4, $\mathbb{R}$) toolkit</td><td>5</td></tr><tr><td>4</td><td>Quantum dynamics</td><td>10</td></tr><tr><td>5</td><td>Quantum state</td><td>13</td></tr><tr><td>6</td><td>Decoherence</td><td>19</td></tr><tr><td>7</td><td>Conclusion</td><td>23</td></tr></table>

<table><tr><td>Appendix A: Commutation relations in the Lie algebra sp(4, $\mathbb{R}$)</td><td>25</td></tr><tr><td>Appendix B: Evolved vacuum state</td><td>26</td></tr><tr><td>Appendix C: Reduced density matrix</td><td>28</td></tr><tr><td>Appendix D: Equivalence between tracing out in Hilbert space and marginalisation in the phase space</td><td>29</td></tr><tr><td>References</td><td>31</td></tr></table>

## 1 Introduction

Two-mode squeezed states [1–3] have been widely studied in the past for the important role they play in quantum optics (see e.g. Refs. [4,5] for reviews), but also in the cosmological context where they describe primordial density perturbations, amplified by gravitational instability from the vacuum quantum fluctuations [6–12]. In general, they characterise the quantum state of linear single-field systems, where each pair of Fourier modes is placed in a two-mode squeezed state [13].

When more degrees of freedom are present however, two-mode squeezed state are insufficient and the squeezing formalism needs to be generalised to higher numbers of modes. This is why, in this work, we construct the four-mode squeezed states, which describe two linearly-coupled quantum scalar fields. The motivation behind this analysis is twofold. First, there are a number of situations where two-field systems are directly relevant, for instance during the inflationary phase our primordial universe underwent. Even if current cosmological data is consistent with single-field setups [14], from a theoretical point of view, inflation takes place in a regime that is far beyond the reach of accelerators, and most physical setups that have been proposed to embed inflation contain extra scalar fields, for instance in the string-theoretic context [15–19]. Those additional degrees of freedom are usually associated to entropic perturbations. Four-mode squeezed states would then naturally appear in

<footnote n="1">¹ Institut d'Astrophysique Spatiale, Université Paris-Saclay, CNRS, 91405 Orsay, France</footnote>

<footnote n="2">² Laboratoire Astroparticule et Cosmologie, Université Denis Diderot Paris 7, 10 rue Alice Domon et Léonie Duquet, 75013 Paris, France</footnote>

<footnote>ᵃ e-mail: thomas.colas@universite-paris-saclay.fr (corresponding author)</footnote>

<footnote>ᵇ e-mail: julien.grain@universite-paris-saclay.fr</footnote>

<footnote>ᶜ e-mail: vincent.vennin@apc.in2p3.fr</footnote>
two-field inflation models, and provide insight about multi-field cosmology in general. Second, this setup provides a way to investigate environmental effects in the case of linear interactions, by tracing over one of the two fields. More precisely, when the system of observational relevance couples to unobserved degrees of freedom (referred to as the “environment”), quantum entanglement builds up between the system and the environment. This affects observational predictions and also leads to the quantum decoherence of the observable sector [20–22]. This phenomenon is usually investigated by the means of effective methods that only provide results that are perturbative in the interaction strength and that rely on additional assumptions, see e.g. Ref. [23]. By considering that one of the two fields represents the observed system and the other field stands for the environment, the formalism we develop will allow us to go beyond those methods and present exact results.

Let us stress that the explicit construction of the squeezed quantum states, especially in the Fock’s space, is not only of formal interest. As we will explain, it provides important insight into the physical mechanisms at play in the dynamics of those states and in the emergence of peculiar properties such as quantum entanglement. Furthermore, it is required in a number of concrete computations (see for instance Refs. [24–28]).

Although we are inspired by problems formulated in the context of cosmology, it is worth mentioning that the formalism we develop here is generic and broad in applicability. It does not require prior knowledge of the concepts and tools relevant in cosmology (which will only be mentioned in our concluding remarks for illustrative purpose), to which we plan to apply our results in separate publications.

Let us now describe how this article is organised, and highlight its main results. In Sect. 2, we introduce the physical setup describing two free scalar fields, both at the classical and quantum levels, and we highlight the symplectic structure that underlies its phase space. This leads us to introducing the symplectic group in four dimensions, Sp(4, $\mathbb{R}$), which we formally describe in Sect. 3. This section reviews material that may also be found in other references on the same topic (apart from the fully factorised form of group elements, Eq. (3.22), which, up to our knowledge, is a new result), see e.g. Refs. [29–33]. As a consequence, it may be skipped by those readers already familiar with the use of symplectic groups in quantum mechanics. Otherwise, it provides a self-contained presentation of the techniques employed in the rest of the paper. The Hamiltonians leading to the four-mode squeezed states are then built in Sect. 4, where we also construct the evolution operator using these results. In particular, we comment on the physical interpretation of the generators of the Lie algebra, when acting in the Hamiltonian. As an example, we also briefly apply our formalism to describe two massless fields in a cosmological background. Finally, we use the

previous results to write a tractable expression for the evolution operator from which four-mode squeezed states can be obtained. Section 5 is devoted to the explicit construction of the four-mode squeezed states in the Fock basis, see Eqs. (5.6) and (5.7), which constitute one of the main results of this paper. An expansion around the limit where the two fields are uncoupled further allows us to discuss the physical interpretation of these formulas, to interpret the four-mode squeezed states in terms of particle transfer and to link their structure to the relevant microphysical parameters. We also derive the Wigner function of the system, which provides an alternative description of the state in the phase space. Although equivalent to the Fock-space description, its simple Gaussian form, built out of the power spectra of the configuration fields, makes some calculations simpler. Finally, in Sect. 6, we investigate environmental effects by tracing out one of the two fields and studying how the reduced quantum state of the first field is affected. This is done both at the level of the density matrix and of the Wigner function. Two independent calculations of the purity are thus performed, and then expanded in the small-coupling limit where they are shown to lead to the same result. Section 7 presents our conclusions, and the paper ends by four appendices to which various technical aspects of the calculations presented in the main text are deferred.

## 2 Quantum phase space of two free fields

### 2.1 Two scalar fields in a homogeneous and isotropic background

Let us consider two real-valued scalar fields $\phi_1(t, \vec{x})$ and $\phi_2(t, \vec{x})$, with conjugate momenta $\pi_1(t, \vec{x})$ and $\pi_2(t, \vec{x})$. These phase-space coordinates can be arranged into the four-dimensional vector $z(\vec{x}) = (\phi_1(\vec{x}), \phi_2(\vec{x}), \pi_1(\vec{x}), \pi_2(\vec{x}))^\mathrm{T}$ where “T” stands for the transpose and explicit time dependence is dropped for notational convenience. In this work, for simplicity, we focus on the case of free fields, for which the Hamiltonian $H$ is a local quadratic form,

$$H = \frac{1}{2} \int d^3\vec{x} z^\mathrm{T}(\vec{x}) H(t) z(\vec{x}). \qquad (2.1)$$

In the cosmological context, one may view $\phi_1$ and $\phi_2$ as two test fields (i.e. they do not backreact on the background geometry), or as the perturbations of some cosmological fields where cosmological perturbation theory is carried out at leading order. In Eq. (2.1), $H(t)$ is a four-by-four real symmetric matrix, which we assume does not depend on the spatial coordinate $\vec{x}$, so the background on which the fields evolve is homogeneous. Note that $H(t)$ may however involve the
gradient operator $\partial/\partial\vec{x}$ (though to positive powers only, to be compatible with the locality assumption).

Phase space is equipped with the Poisson bracket

$$ \{F, G\} = \int d^3\vec{x} \left[ \frac{\delta F}{\delta\phi_1(\vec{x})} \frac{\delta G}{\delta\pi_1(\vec{x})} - \frac{\delta F}{\delta\pi_1(\vec{x})} \frac{\delta G}{\delta\phi_1(\vec{x})} + \frac{\delta F}{\delta\phi_2(\vec{x})} \frac{\delta G}{\delta\pi_2(\vec{x})} - \frac{\delta F}{\delta\pi_2(\vec{x})} \frac{\delta G}{\delta\phi_2(\vec{x})} \right], \quad (2.2) $$

which can be written in matricial form for the phase-space coordinates as

$$ \{z(\vec{x}), z^T(\vec{y})\} = \Omega \delta^3(\vec{x} - \vec{y}). \quad (2.3) $$

In this expression, $\delta$ is the Dirac distribution, the matrix $\Omega$ is given by

$$ \Omega = \begin{pmatrix} 0 & I_2 \\ -I_2 & 0 \end{pmatrix}, \quad (2.4) $$

where $I_n$ is the $n \times n$ identity matrix, and the notation in Eq. (2.3) has to be understood as $\{z_\mu(\vec{x}), z_\nu(\vec{y})\} = \Omega_{\mu\nu}\delta^3(\vec{x} - \vec{y})$ with $\mu, \nu = 1 \dots 4$. The time evolution of any function $F$ of the phase-space field variables is given by

$$ \dot{F}(z) = \{F(z), H\}, \quad (2.5) $$

where a dot means differentiating with respect to the time variable.

## Fourier space

Since the background is homogeneous, it is useful to work in Fourier space and to introduce

$$ z_{\vec{k}} = \begin{pmatrix} \phi_{1,\vec{k}} \\ \phi_{2,\vec{k}} \\ \pi_{1,\vec{k}} \\ \pi_{2,\vec{k}} \end{pmatrix} = \int \frac{d^3\vec{x}}{(2\pi)^{3/2}} z(\vec{x}) e^{-i\vec{k}\cdot\vec{x}}. \quad (2.6) $$

The condition for the fields to be real-valued, $z(\vec{x}) = z^*(\vec{x})$, translates into $z_{\vec{k}}^* = z_{-\vec{k}}$. This means that half of the Fourier modes are enough to parametrise the entire (now complex) phase space. In practice, any integral over $k \in \mathbb{R}^3$ can be split into an integral over $\mathbb{R}^{3+} \equiv \mathbb{R}^2 \times \mathbb{R}^+$ and $\mathbb{R}^{3-} \equiv \mathbb{R}^2 \times \mathbb{R}^-$, where the latter can be related to the former by a simple change of integration variable $\vec{k} \to -\vec{k}$ and using the relation $z_{\vec{k}}^* = z_{-\vec{k}}$. For the Hamiltonian (2.1), this leads to

$$ H = \int_{\mathbb{R}^{3+}} d^3\vec{k} z_{\vec{k}}^{\dagger} \mathbf{H}_k(t) z_{\vec{k}}, \quad (2.7) $$

which avoids double counting the degrees of freedom of the theory. In Eq. (2.7), $\mathbf{H}_k(t)$ corresponds to $\mathbf{H}(t)$ where the spatial gradient $\partial/\partial\vec{x}$ is replaced with $i\vec{k}$. If we further assume the background to be isotropic, $\mathbf{H}_k(t)$ only depends on the norm $k = |\vec{k}|$ of $\vec{k}$, hence the notation.

Plugging Eq. (2.6) into Eq. (2.2), one can compute

$$ \{z_{\vec{k}}, z_{\vec{q}}^{\dagger}\} = \Omega \delta^3(\vec{k} - \vec{q}), \quad (2.8) $$

which is of the same form as Eq. (2.3). This shows that the Poisson brackets are preserved when going to Fourier space (in the language that will be introduced in Sect. 2.2, the Fourier transform is a “symplectic transformation”), so phase space can be equivalently parametrised with the Fourier coordinates, and the Poisson bracket (2.2) can also be written as

$$ \{F, G\} = \int_{\mathbb{R}^3} d^3\vec{k} \left( \frac{\delta F}{\delta\phi_{1,\vec{k}}} \frac{\delta G}{\delta\pi_{1,\vec{k}}} - \frac{\delta F}{\delta\pi_{1,\vec{k}}} \frac{\delta G}{\delta\phi_{1,\vec{k}}} + \frac{\delta F}{\delta\phi_{2,\vec{k}}} \frac{\delta G}{\delta\pi_{2,\vec{k}}} - \frac{\delta F}{\delta\pi_{2,\vec{k}}} \frac{\delta G}{\delta\phi_{2,\vec{k}}} \right). \quad (2.9) $$

Finally, applying the equation of motion (2.5) to the Fourier phase-space coordinates yields

$$ \dot{z}_{\vec{k}} = (\Omega \mathbf{H}_k) z_{\vec{k}}, \quad (2.10) $$

where the explicit time dependence of $\mathbf{H}_k$ has been dropped for notational convenience. This equation shows that, for free fields in a homogeneous and isotropic background, Fourier modes decouple and evolve independently. This implies that each Fourier sector can be studied separately, which greatly simplifies the analysis. From now on, we therefore focus on a single Fourier sector, i.e. $\vec{k}$ in $\mathbb{R}^{3+}$ is fixed hereafter.

## Quantisation

So far, linear Hamiltonian systems have been described at the classical level. We now follow the canonical quantisation prescriptions and promote the phase-space coordinates to quantum operators acting on a Hilbert space, $\hat{z}_{\vec{k}} := (\hat{\phi}_{1,\vec{k}}, \hat{\phi}_{2,\vec{k}}, \hat{\pi}_{1,\vec{k}}, \hat{\pi}_{2,\vec{k}})^T$ (where, from now on, hats denote quantum operators). These operators satisfy the canonical commutation relations

$$ [\hat{z}_{\vec{k}}, \hat{z}_{\vec{q}}^{\dagger}] = i \Omega \delta^3(\vec{k} - \vec{q}), \quad (2.11) $$

which is the quantum analogue of Eq. (2.8) and where, hereafter, we work with $\hbar = 1$. The Hamiltonian operator reads $\hat{H} = \int_{\mathbb{R}^{3+}} d^3\vec{k} \hat{z}_{\vec{k}}^{\dagger} \mathbf{H}_k \hat{z}_{\vec{k}}$, and the dynamics of any function of the quantum phase-space variables is given by the Heisenberg equation $\dot{\hat{F}}(\hat{z}_{\vec{k}}) = -i[F(\hat{z}_{\vec{k}}), \hat{H}]$. In particular, for the field variables themselves, this gives rise to $\hat{z}_{\vec{k}} = (\Omega \mathbf{H}_k) \hat{z}_{\vec{k}}$, which directly transposes Eq. (2.10).

Creation and annihilation operators are defined in the usual way, i.e.

$$ \hat{a}_{j,\vec{k}} = \frac{1}{\sqrt{2}} \left( \sqrt{k} \hat{\phi}_{j,\vec{k}} + \frac{i}{\sqrt{k}} \hat{\pi}_{j,\vec{k}} \right) \quad \text{for } j = 1, 2, \quad (2.12) $$
where the prefactors $\sqrt{k}$ and $1/\sqrt{k}$ are introduced for dimensional reasons. This can be written in matricial form as

$$ \hat{a}_{\vec{k}} = \begin{pmatrix} \hat{a}_{1,\vec{k}} \\ \hat{a}_{2,\vec{k}} \\ \hat{a}_{1,-\vec{k}}^{\dagger} \\ \hat{a}_{2,-\vec{k}}^{\dagger} \end{pmatrix} = U D_k \hat{z}_{\vec{k}}, \quad (2.13) $$

where “†” stands for the conjugate transpose, and the matrices $U$ and $D_k$ are defined as

$$ \begin{gathered} U = \frac{1}{\sqrt{2}} \begin{pmatrix} I_2 & i I_2 \\ I_2 & -i I_2 \end{pmatrix} \quad \text{and} \\ D_k = \begin{pmatrix} \sqrt{k} I_2 & 0 \\ 0 & I_2/\sqrt{k} \end{pmatrix}. \end{gathered} \quad (2.14) $$

One can check that $U$ is a unitary matrix, i.e. $UU^\dagger = U^\dagger U = I_4$. In principle, $D_k$ may be replaced with $M_k D_k$, where $M_k$ is any (dimensionless) symplectic matrix (formally defined below in Sect. 2.2). It only leads to a different definition of the vacuum state (i.e. the state that is annihilated by the annihilation operators) [13]. Let us note that the ordering in $\hat{a}_{\vec{k}}$ is different than in $\hat{z}_{\vec{k}}$, since in $\hat{a}_{\vec{k}}$ the first two entries concern the $\vec{k}$ sector and the last two entries the $-\vec{k}$ sector.<fnref n="1" />

The dynamics of the creation and annihilation operators is generated by the quadratic Hamiltonian $\hat{\mathcal{H}} = \int_{\mathbb{R}^{3+}} d^3\vec{k} \hat{a}_{\vec{k}}^\dagger \mathcal{H}_k \hat{a}_{\vec{k}}$, where $\mathcal{H}_k$ reads [13]

$$ \mathcal{H}_k = U \left[ (D_k^{-1})^\mathrm{T} \mathcal{H}_k D_k^{-1} + (D_k^{-1})^\mathrm{T} \mathbf{\Omega} \dot{\mathbf{D}}_k^{-1} \right] U^\dagger. \quad (2.15) $$

The canonical commutation relations are given by

$$ \begin{gathered} [\hat{a}_{\vec{k}}, \hat{a}_{\vec{q}}^{\dagger}] = i \mathcal{J} \delta^3(\vec{k} - \vec{q}) \\ \text{where } \mathcal{J} = U \mathbf{\Omega} U^{\dagger} = -i \begin{pmatrix} I_2 & 0 \\ 0 & -I_2 \end{pmatrix}, \end{gathered} \quad (2.16) $$

and any function of the creation and annihilation operators evolves according to $\dot{F}(\hat{a}_{\vec{k}}) = [F(\hat{a}_{\vec{k}}), \hat{\mathcal{H}}]$. In particular, for the creation and annihilation operators themselves, one obtains

$$ \dot{\hat{a}}_{\vec{k}} = (\mathcal{J}\mathcal{H}_k)\hat{a}_{\vec{k}}. \quad (2.17) $$

## 2.2 Symplectic structure of the phase space

In Sect. 2.1, we saw that the Fourier transform preserves the Poisson brackets of the phase-space variables. Another example of a transformation that preserves the Poisson brackets

is provided by the Hamiltonian evolution itself. Indeed, the equation of motion (2.10) can be solved as

$$ z_{\vec{k}}(t) = \mathcal{G}_k(t, t_{\mathrm{in}}) z_{\vec{k}}(t_{\mathrm{in}}), \quad (2.18) $$

where $\mathcal{G}_k$ is a $(4 \times 4)$-real matrix called the Green's matrix and that satisfies $\dot{\mathcal{G}}_k = \mathbf{\Omega} \mathcal{H}_k \mathcal{G}_k + I_4 \delta(t - t_{\mathrm{in}})$, with initial condition $\mathcal{G}_k(t_{\mathrm{in}}, t_{\mathrm{in}}) = I_4$. Note that, as $\mathcal{H}_k$, $\mathcal{G}_k$ only depends on the wavenumber $k$. One can also check that $\mathcal{G}_k$ satisfies

$$ \mathcal{G}_k^\mathrm{T} \mathbf{\Omega} \mathcal{G}_k = \mathbf{\Omega}. \quad (2.19) $$

This is indeed obviously the case at initial time, and by plugging the equation of motion for $\mathcal{G}_k$ in the time derivative of the left-hand-side of Eq. (2.19), one obtains a vanishing result after using that $\mathcal{H}_k$ is symmetric and that $\mathbf{\Omega}^\mathrm{T} \mathbf{\Omega} = -\mathbf{\Omega}^2 = I_4$.

In general, real matrices satisfying Eq. (2.19) are called symplectic, and they form the symplectic group $\mathrm{Sp}(4, \mathbb{R})$. They describe all possible reparametrisations of phase space through linear canonical transformations, i.e. transformations that preserve the Poisson brackets. Indeed, consider two phase-space coordinates $z_{\vec{k}}$ and $\tilde{z}_{\vec{k}}$, related through a linear transformation

$$ \tilde{z}_{\vec{k}} = \mathcal{M}_k z_{\vec{k}}. \quad (2.20) $$

One can check that the Poisson brackets are preserved, i.e. $\{\tilde{z}_{\vec{k}}, \tilde{z}_{\vec{k}}^{\dagger}\} = \{z_{\vec{k}}, z_{\vec{k}}^{\dagger}\} = \mathbf{\Omega}$, if and only if $\mathcal{M}_k \in \mathrm{Sp}(4, \mathbb{R})$ [i.e. $\mathcal{M}_k$ satisfies Eq. (2.19)]. This ensures that the Poisson bracket between two arbitrary phase-space functions is the same when calculated with the $z_{\vec{k}}$-variables or with the $\tilde{z}_{\vec{k}}$-variables. One can check that the equations of motion for the new set of canonical variables $\tilde{z}_{\vec{k}}$ are then given by Hamilton equations with the new Hamiltonian kernel

$$ \tilde{\mathcal{H}}_k = (\mathcal{M}_k^{-1})^\mathrm{T} \mathcal{H}_k \mathcal{M}_k^{-1} + (\mathcal{M}_k^{-1})^\mathrm{T} \mathbf{\Omega} \dot{\mathbf{M}}_k^{-1}, \quad (2.21) $$

and the dynamics is solved by $\tilde{z}_{\vec{k}}(t) = \tilde{\mathcal{G}}_k(t, t_{\mathrm{in}}) \tilde{z}_{\vec{k}}(t_{\mathrm{in}})$ where $\tilde{\mathcal{G}}_k(t, t_{\mathrm{in}}) = \mathcal{M}_k(t) \mathcal{G}_k(t, t_{\mathrm{in}}) \mathcal{M}_k^{-1}(t_{\mathrm{in}})$. Symplectic transformations thus constitute a fundamental symmetry of the Hamiltonian phase space, and this is why the symplectic group for linear scalar-field systems is the main topic of the present work.

Let us note that in this framework, the dynamical evolution is nothing but a particular symplectic transformation since, as stressed above, the Green's matrix is symplectic. One can also check that for the transformation that goes from $z_{\vec{k}}(t)$ to $z_{\vec{k}}(t_{\mathrm{in}})$, generated by the matrix $\mathcal{G}_k^{-1}(t, t_{\mathrm{in}})$, the new Hamiltonian vanishes, as can be shown by plugging the equation of motion for $\mathcal{G}_k$ into Eq. (2.21). This is consistent with the fact that the $z_{\vec{k}}(t_{\mathrm{in}})$ variables are indeed time independent, and in that case the dynamics is entirely contained in the canonical transformation that relates them with the primary variables $z_{\vec{k}}(t)$.

<footnote n="1">¹ This is because the classical version of Eq. (2.13) gives $z_{-\vec{k}} = D_k^{-1} U^\dagger a_{-\vec{k}}$ and $z_{\vec{k}}^* = D_k^{-1} U^\mathrm{T} a_{\vec{k}}^*$. The reality condition $z_{\vec{k}}^* = z_{-\vec{k}}$ thus entails $a_{-\vec{k}} = U U^\mathrm{T} a_{\vec{k}}^*$, where $UU^\mathrm{T} = \begin{pmatrix} 0 & I_2 \\ I_2 & 0 \end{pmatrix}$, explaining the structure of the $\hat{a}_{\vec{k}}$ vector.</footnote>
Finally, when one works with the creation and annihilation operators introduced in Sect. 2.1, a similar description applies. The equation of motion given in Eq. (2.17) can be solved in terms of the Green’s matrix,

$$ \hat{a}_{\vec{k}}(t) = \mathcal{G}_k(t, t_{\text{in}}) \hat{a}_{\vec{k}}(t_{\text{in}}), \qquad (2.22) $$

where $\mathcal{G}_k(t, t_{\text{in}}) = U D_k \mathcal{G}_k(t, t_{\text{in}}) D_k^{-1} U^\dagger$ such that $\mathcal{G}_k^\dagger \mathcal{J} \mathcal{G}_k = \mathcal{J}$ and $\det(\mathcal{G}_k) = 1$. Similarly, one can check that any generic canonical transformation $\hat{z}_{\vec{k}} \rightarrow M_k \hat{z}_{\vec{k}}$ gives rise to $\hat{a}_{\vec{k}} \rightarrow M_k \hat{a}_{\vec{k}}$, where $M_k = U D_k M_k D_k^{-1} U^\dagger$ satisfies $M_k^\dagger \mathcal{J} M_k = \mathcal{J}$ and $\det(M_k) = 1$.

### 3 Sp(4, $\mathbb{R}$) toolkit

In the previous section, we have seen how symplectic transformations naturally arise in the phase-space description of linear Hamiltonian systems, both as a fundamental reparametrisation symmetry and as a way to generate the dynamics. This is why in this section, we further study the mathematical structure of the symplectic group in four dimensions, which is relevant to discuss the physics of two scalar fields. Readers already familiar with the use of symplectic groups in quantum mechanics can easily skip this section, which mostly consists in a review of the mathematical tools employed in the rest of the paper. It may otherwise serve as a pedagogical introduction to the techniques employed in the subsequent calculations, and set out our main notations.

#### 3.1 Generators and Lie algebra

As explained around Eq. (2.19), the group of symplectic $(4 \times 4)$-matrices, denoted $\text{Sp}(4, \mathbb{R})$, is defined as

$$ \text{Sp}(4, \mathbb{R}) = \{ \mathcal{M} \in \mathcal{M}_4(\mathbb{R}) : \mathcal{M}^\mathrm{T} \Omega \mathcal{M} = \Omega \}, \qquad (3.1) $$

where $\mathcal{M}_n(\mathbb{R})$ is the set of $(n \times n)$-real matrices. Since $\Omega^\mathrm{T} \Omega = -\Omega^2 = I_4$, one can show that, if $\mathcal{M} \in \text{Sp}(4, \mathbb{R})$, then $\mathcal{M}^\mathrm{T} \in \text{Sp}(4, \mathbb{R})$ (this is because Eq. (3.1) leads to $\mathcal{M}^\mathrm{T} = -\Omega \mathcal{M}^{-1} \Omega$, which implies that $\mathcal{M} \Omega \mathcal{M}^\mathrm{T} = \Omega$). One can also readily check that $\text{Sp}(4, \mathbb{R})$ is indeed a group, and it follows from Eq. (3.1) that all symplectic matrices have unit determinant [34].

$\text{Sp}(4, \mathbb{R})$ is a Lie group, that is a continuous group whose multiplication and inversion operations are differentiable, so one can investigate global properties of the group by looking at its local or linearised version, given in terms of its so-called Lie algebra. The Lie algebra is a vector space under the bracket operation $[X, Y] = XY - YX$ that completely captures the local structure of the group. When analysing Lie groups, finding a basis of this vector space, i.e. a set of so-called “generators”, and deriving their commutators, is of paramount importance. The number of generators spec-

fies both the dimension of the Lie group and its associated Lie algebra. In what follows, the Lie algebra of $\text{Sp}(4, \mathbb{R})$ is denoted $\text{sp}(4, \mathbb{R})$.

The exponential map allows us to connect the Lie algebra to its corresponding Lie group. For any $X \in \text{sp}(4, \mathbb{R})$, $\mathcal{M} = \exp(X)$ is an element of $\text{Sp}(4, \mathbb{R})$ (see Ref. [33] for a recent complete characterisation of the exponential map of $\text{Sp}(4, \mathbb{R})$). Note that the exponential map is not surjective, meaning that there exist elements of $\text{Sp}(4, \mathbb{R})$ that cannot be written as $\exp(X)$. By plugging $\mathcal{M} = \exp(X)$ into Eq. (3.1), and upon writing the obtained formula as an expansion in $X$, one finds that the Lie algebra is defined according to

$$ \text{sp}(4, \mathbb{R}) = \{ X \in \mathcal{M}_4(\mathbb{R}) : \Omega X + X^\mathrm{T} \Omega = 0 \}. \qquad (3.2) $$

This allows one to write a generic element of the Lie algebra as

$$ X = \begin{pmatrix} A & B \\ C & -A^\mathrm{T} \end{pmatrix} \qquad (3.3) $$

where $A, B$ and $C$ are three $(2 \times 2)$ real matrices, $B$ and $C$ being symmetric. This implies that $\text{sp}(4, \mathbb{R})$ is a ten-dimensional vector space, since $A$ contains four degrees of freedom while $B$ and $C$ being symmetric, they contain three degrees of freedom each. Denoting by $K_i, i = 1 \cdots 10$, a basis of 10 generators, a generic element of the Lie algebra can be decomposed as

$$ X = \sum_{i=1}^{10} \alpha_i K_i, \quad \alpha_i \in \mathbb{R}. \qquad (3.4) $$

Our next step is to exhibit such a basis.

To that end, we first introduce the Kronecker product. For two $(2 \times 2)$-matrices

$$ A = \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix} \quad \text{and} \quad B = \begin{pmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{pmatrix}, \qquad (3.5) $$

the Kronecker product $A \otimes B$ is a $(4 \times 4)$-matrix defined as

$$ \begin{aligned} A \otimes B &= \begin{pmatrix} a_{11}B & a_{12}B \\ a_{21}B & a_{22}B \end{pmatrix} \\ &= \begin{pmatrix} a_{11}b_{11} & a_{11}b_{12} & a_{12}b_{11} & a_{12}b_{12} \\ a_{11}b_{21} & a_{11}b_{22} & a_{12}b_{21} & a_{12}b_{22} \\ a_{21}b_{11} & a_{21}b_{12} & a_{22}b_{11} & a_{22}b_{12} \\ a_{21}b_{21} & a_{21}b_{22} & a_{22}b_{21} & a_{22}b_{22} \end{pmatrix}. \end{aligned} \qquad (3.6) $$

The reason why this construction is useful is because, within the two-field system at hand, each field is individually described by an $\text{sp}(2, \mathbb{R})$ algebra. It is therefore natural for expect that $\text{sp}(4, \mathbb{R})$ contains products of elements of $\text{sp}(2, \mathbb{R})$ with themselves. The generators of $\text{sp}(2, \mathbb{R})$ can be simply written in terms of the Pauli matrices [13],

$$ \text{sp}(2, \mathbb{R}) = \text{Span}\{\sigma_x, i\sigma_y, \sigma_z\}, \qquad (3.7) $$
<table><caption>**Table 1** Generators of sp(4, R) in the fundamental representation</caption><thead><tr><th>Squeezing</th><th>Rotation</th><th>Boost</th></tr></thead><tbody><tr><td><b>K</b><sub>1</sub> = σ<sub>z</sub> ⊗ σ<sub>z</sub> = <span style="font-size: 2.5em; vertical-align: middle;">(</span><table style="display: inline-table; vertical-align: middle; text-align: center;"><tbody><tr><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0</td><td>-1</td><td>0</td><td>0</td></tr><tr><td>0</td><td>0</td><td>-1</td><td>0</td></tr><tr><td>0</td><td>0</td><td>0</td><td>1</td></tr></tbody></table><span style="font-size: 2.5em; vertical-align: middle;">)</span></td><td><b>K</b><sub>3</sub> = iσ<sub>y</sub> ⊗ <b>I</b><sub>2</sub> = <span style="font-size: 2.5em; vertical-align: middle;">(</span><table style="display: inline-table; vertical-align: middle; text-align: center;"><tbody><tr><td>0</td><td>0</td><td>1</td><td>0</td></tr><tr><td>0</td><td>0</td><td>0</td><td>1</td></tr><tr><td>-1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0</td><td>-1</td><td>0</td><td>0</td></tr></tbody></table><span style="font-size: 2.5em; vertical-align: middle;">)</span></td><td><b>K</b><sub>7</sub> = σ<sub>x</sub> ⊗ <b>I</b><sub>2</sub> = <span style="font-size: 2.5em; vertical-align: middle;">(</span><table style="display: inline-table; vertical-align: middle; text-align: center;"><tbody><tr><td>0</td><td>0</td><td>1</td><td>0</td></tr><tr><td>0</td><td>0</td><td>0</td><td>1</td></tr><tr><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0</td><td>1</td><td>0</td><td>0</td></tr></tbody></table><span style="font-size: 2.5em; vertical-align: middle;">)</span></td></tr><tr><td><b>K</b><sub>2</sub> = σ<sub>z</sub> ⊗ <b>I</b><sub>2</sub> = <span style="font-size: 2.5em; vertical-align: middle;">(</span><table style="display: inline-table; vertical-align: middle; text-align: center;"><tbody><tr><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td>0</td><td>0</td><td>-1</td><td>0</td></tr><tr><td>0</td><td>0</td><td>0</td><td>-1</td></tr></tbody></table><span style="font-size: 2.5em; vertical-align: middle;">)</span></td><td><b>K</b><sub>4</sub> = iσ<sub>y</sub> ⊗ σ<sub>z</sub> = <span style="font-size: 2.5em; vertical-align: middle;">(</span><table style="display: inline-table; vertical-align: middle; text-align: center;"><tbody><tr><td>0</td><td>0</td><td>1</td><td>0</td></tr><tr><td>0</td><td>0</td><td>0</td><td>-1</td></tr><tr><td>-1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0</td><td>1</td><td>0</td><td>0</td></tr></tbody></table><span style="font-size: 2.5em; vertical-align: middle;">)</span></td><td><b>K</b><sub>8</sub> = σ<sub>x</sub> ⊗ σ<sub>z</sub> = <span style="font-size: 2.5em; vertical-align: middle;">(</span><table style="display: inline-table; vertical-align: middle; text-align: center;"><tbody><tr><td>0</td><td>0</td><td>1</td><td>0</td></tr><tr><td>0</td><td>0</td><td>0</td><td>-1</td></tr><tr><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0</td><td>-1</td><td>0</td><td>0</td></tr></tbody></table><span style="font-size: 2.5em; vertical-align: middle;">)</span></td></tr><tr><td><b>K</b><sub>5</sub> = <b>I</b><sub>2</sub> ⊗ iσ<sub>y</sub> = <span style="font-size: 2.5em; vertical-align: middle;">(</span><table style="display: inline-table; vertical-align: middle; text-align: center;"><tbody><tr><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td>-1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0</td><td>0</td><td>1</td><td>0</td></tr><tr><td>0</td><td>0</td><td>-1</td><td>0</td></tr></tbody></table><span style="font-size: 2.5em; vertical-align: middle;">)</span></td><td><b>K</b><sub>9</sub> = σ<sub>z</sub> ⊗ σ<sub>x</sub> = <span style="font-size: 2.5em; vertical-align: middle;">(</span><table style="display: inline-table; vertical-align: middle; text-align: center;"><tbody><tr><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0</td><td>0</td><td>0</td><td>-1</td></tr><tr><td>0</td><td>0</td><td>-1</td><td>0</td></tr></tbody></table><span style="font-size: 2.5em; vertical-align: middle;">)</span></td></tr><tr><td><b>K</b><sub>6</sub> = iσ<sub>y</sub> ⊗ σ<sub>x</sub> = <span style="font-size: 2.5em; vertical-align: middle;">(</span><table style="display: inline-table; vertical-align: middle; text-align: center;"><tbody><tr><td>0</td><td>0</td><td>0</td><td>1</td></tr><tr><td>0</td><td>0</td><td>1</td><td>0</td></tr><tr><td>0</td><td>-1</td><td>0</td><td>0</td></tr><tr><td>-1</td><td>0</td><td>0</td><td>0</td></tr></tbody></table><span style="font-size: 2.5em; vertical-align: middle;">)</span></td><td><b>K</b><sub>10</sub> = σ<sub>x</sub> ⊗ σ<sub>x</sub> = <span style="font-size: 2.5em; vertical-align: middle;">(</span><table style="display: inline-table; vertical-align: middle; text-align: center;"><tbody><tr><td>0</td><td>0</td><td>0</td><td>1</td></tr><tr><td>0</td><td>0</td><td>1</td><td>0</td></tr><tr><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td>1</td><td>0</td><td>0</td><td>0</td></tr></tbody></table><span style="font-size: 2.5em; vertical-align: middle;">)</span></td></tr></tbody></table>

with

$$ \sigma_x = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad \sigma_y = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}, \quad \sigma_z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}, \qquad (3.8) $$

so sp(2, R) is of dimension 3. This implies that over the 10 generators of sp(4, R), 6 provide two copies of sp(2, R) and describe the two sectors separately, and 4 are related to the coupling between the two sectors. For $\sigma_a, \sigma_b \in \{I_2, \sigma_x, i\sigma_y, \sigma_z\}$, there are 16 combinations $\sigma_a \otimes \sigma_b$ for $a, b \in \{0, \dots, 3\}$, namely $I_4$ and the 15 Dirac matrices [29]. Among them, only 10 are of the form (3.3),<fnref n="2" /> which are listed in Table 1, where they are organised in three subsets. The so-called squeezing generators are diagonal. Through the exponential map, they give rise to group elements of the form $\exp(d_1 K_1) = \text{diag}(e^{d_1}, e^{-d_1}, e^{-d_1}, e^{d_1})$ and $\exp(d_2 K_2) = \text{diag}(e^{d_2}, e^{d_2}, e^{-d_2}, e^{-d_2})$ with $d_1$ and $d_2$ two real parameters. Therefore, they elongate one phase-space direction while contracting the other, hence their name. The two other kinds of generators are called rotations and boosts. When squared, rotations give $-I_4$ and boosts give $I_4$. Therefore, once exponentiated, rotations generate group elements of the form $\exp(\theta_i K_i) = \cos \theta_i I_4 + \sin \theta_i K_i$ for $i \in \{3, 4, 5, 6\}$ with $\theta_i \in \mathbb{R}$, while boosts generate $\exp(\alpha_i K_i) = \cosh \alpha_i I_4 + \sinh \alpha_i K_i$ for $i \in \{7, 8, 9, 10\}$ with $\alpha_i \in \mathbb{R}$, hence their name.

Now that we have explicitly obtained the generators of the Lie algebra, let us identify their respective role. The two sp(2, R) algebras are given by $\{(K_1 + K_2)/2, (K_3 + K_4)/2, (K_7+K_8)/2\}$ and $\{(K_1-K_2)/2, (K_3-K_4)/2, (K_7-K_8)/2\}$, each of them being composed of a squeezing, a rotation and a boost. They act on each sector separately and would be enough to describe two non-interacting degrees of free-

dom. Formally, they generate the Sp(2, R) × Sp(2, R) subgroup of Sp(4, R). The 4 remaining generators are associated with the coupling between the two sectors. They correspond to the rotations $K_5$ and $K_6$ and the boosts $K_9$ and $K_{10}$, which entangle the two sectors, as will be made clear in Sect. 3.4. One can indeed check that these 4 generators (and only them) have non-vanishing off-diagonal elements within the $(2 \times 2)$ blocks, which, from the ordering of the phase-space variables in Eq. (2.6), implies that they mix the two fields. This is why, hereafter, they will be referred to as the coupling generators.

To complete our description of the Lie algebra, let us finally provide the commutators between its generators. They can be obtained from the multiplication rule of the Kronecker product for square matrices,<fnref n="3" />

$$ (A \otimes B) (C \otimes D) = (AC) \otimes (BD), \qquad (3.9) $$

together with the formula

$$ \sigma_i \sigma_j = \delta_{ij} I_2 + i\epsilon_{ijk} \sigma_k \quad \text{for } \{i, j, k\} \in \{x, y, z\}, \qquad (3.10) $$

where $\delta_{ab}$ is the Kronecker symbol and $\epsilon_{abc}$ is the Levi-Civita symbol. For instance, one has $[K_6, K_4] = [i\sigma_y \otimes \sigma_x, i\sigma_y \otimes \sigma_z] = (i\sigma_y)^2 \otimes \sigma_x \sigma_z - (i\sigma_y)^2 \otimes \sigma_z \sigma_x = -I_2 \otimes (-i\sigma_y) + I_2 \otimes (i\sigma_y) = 2K_5$. All other commutators are presented in Appendix A, where we also derive the various subalgebras. We do not reproduce these formulas here for display convenience, and will simply refer to Appendix A when needed.

<footnote n="2">² Formally, the 15 Dirac matrices form the $\mathfrak{so}(3, 3)$ algebra. sp(4, R) is isomorphic to $\mathfrak{so}(3, 2)$, which is a subalgebra of $\mathfrak{so}(3, 3)$ [29].</footnote>

<footnote n="3">³ Other properties of the Kronecker product that will be used in the following are that it is bilinear and associative, and that for square matrices, $(A \otimes B)^{-1} = A^{-1} \otimes B^{-1}$, $(A \otimes B)^* = A^* \otimes B^*$ and $(A \otimes B)^T = A^T \otimes B^T$.</footnote>
## 3.2 Bloch–Messiah decomposition

The explicit derivation of the 10 generators of the Lie algebra allows one to decompose any element of the group (within the exponential map) onto these generators, upon exponentiating Eq. (3.4). This parametrisation of the exponential map is performed in Ref. [33] (see also Refs. [31,32]). When all 10 parameters are non vanishing, it however leads to expressions that may be cumbersome to manipulate, and which, as mentioned above, do not reach all the elements of the group. This is why, in this section, we turn our attention to an alternative decomposition, the so-called Bloch–Messiah decomposition (also sometimes called Euler decomposition) [35], which allows one to write any symplectic matrix as

$$ M(\theta, d, \varphi) = R(\theta)Z(d)R(\varphi). \qquad (3.11) $$

Here, $R(\theta)$, $R(\varphi) \in \text{Sp}(4, \mathbb{R}) \cap \text{SO}(4)$ are constructed from the four rotation generators and $Z(d)$ from the two squeezing generators, i.e.

$$ R(\theta) = \exp(\theta_3 \mathbf{K}_3 + \theta_4 \mathbf{K}_4 + \theta_5 \mathbf{K}_5 + \theta_6 \mathbf{K}_6) \qquad (3.12) $$

and a similar expression for $R(\varphi)$ with $\varphi = (\varphi_1, \varphi_2, \varphi_3, \varphi_4)$, and

$$ Z(d) = \exp(d_1 \mathbf{K}_1 + d_2 \mathbf{K}_2). \qquad (3.13) $$

The 8 parameters contained in $\theta$ and $\varphi$ are called the rotation parameters, while $d$ contains the so-called squeezing parameters. Note that the parameters associated with the coupling generators are $\theta_5, \varphi_5, \theta_6$ and $\varphi_6$, which thus control the mixing between the two sectors.

One may note that only 6 out of the 10 generators of $\text{Sp}(4, \mathbb{R})$ are involved in the Bloch–Messiah decomposition. It therefore provides a factorised expression of the group elements that is slightly more convenient to manipulate. Moreover, the three blocks of the decomposition can be further factorised down. Indeed, in Appendix A, it is shown that $\mathbf{K}_1$ and $\mathbf{K}_2$ commute, so

$$ Z(d) = \exp(d_1 \mathbf{K}_1) \cdot \exp(d_2 \mathbf{K}_2). \qquad (3.14) $$

Regarding the rotation operators, still in Appendix A, it is shown that $\mathbf{K}_3$ commutes with the other three rotation generators, namely $\mathbf{K}_4, \mathbf{K}_5$ and $\mathbf{K}_6$ [see Eq. (A.2)], so $\mathbf{K}_3$ generates a separate $\text{U}(1)$ Lie group; while $\mathbf{K}_4, \mathbf{K}_5$ and $\mathbf{K}_6$ generate a $\text{SU}(2)$ Lie group [see Eq. (A.20)]. As a consequence,

the four rotation generators form$^4$

$$ \text{Sp}(4, \mathbb{R}) \cap \text{SO}(4) \cong \text{U}(2) \cong \text{SU}(2) \times \text{U}(1), \qquad (3.15) $$

where “$\cong$” indicates group isomorphisms. This leads to factorising

$$ R(\theta) = \exp(\theta_3 \mathbf{K}_3) \cdot \exp(\theta_4 \mathbf{K}_4 + \theta_5 \mathbf{K}_5 + \theta_6 \mathbf{K}_6), \qquad (3.16) $$

and a similar expression for $R(\varphi)$. Finally, one can use the Baker–Campbell–Haussdorf formula [36–39] to further factorise the remaining $\text{SU}(2)$ part. This can be done by first introducing the complexified Lie algebra

$$ \begin{aligned} S_z &= \frac{1}{2i} \mathbf{K}_4, & S_+ &= \frac{1}{2i} (\mathbf{K}_5 - i \mathbf{K}_6), \\ S_- &= \frac{1}{2i} (\mathbf{K}_5 + i \mathbf{K}_6), & & \end{aligned} \qquad (3.17) $$

where the notation is purposely reminiscent of spin physics. One can check that $S_z^\dagger = S_z$ and $S_+^\dagger = S_-$, and from the $\text{SU}(2)$ commutation relations derived in Eq. (A.20), one has

$$ \begin{aligned} [S_z, S_+] &= S_+, & [S_z, S_-] &= -S_-, \\ [S_+, S_-] &= 2S_z. & & \end{aligned} \qquad (3.18) $$

Expanding $\theta_4 \mathbf{K}_4 + \theta_5 \mathbf{K}_5 + \theta_6 \mathbf{K}_6$ onto $S_z, S_+$ and $S_-$, the Baker–Campbell–Haussdorf formula thus leads to [39]

$$ \begin{aligned} & \exp(\theta_4 \mathbf{K}_4 + \theta_5 \mathbf{K}_5 + \theta_6 \mathbf{K}_6) \\ &= \exp(p_+ \mathbf{S}_+) \exp(p_z \mathbf{S}_z) \exp(p_- \mathbf{S}_-), & \end{aligned} \qquad (3.19) $$

where

$$ \begin{aligned} p_z &= -2 \ln \left( \cos \theta - i \frac{\theta_4}{\theta} \sin \theta \right), \\ p_- &= \frac{-\tau^* \sin \theta}{\theta \cos \theta - i \theta_4 \sin \theta}, \\ p_+ &= \frac{\tau \sin \theta}{\theta \cos \theta - i \theta_4 \sin \theta}, & \end{aligned} \qquad (3.20) $$

and where we have introduced

$$ \tau = -\theta_6 + i\theta_5 \quad \text{and} \quad \theta = \sqrt{\theta_4^2 + \theta_5^2 + \theta_6^2}. \qquad (3.21) $$

A similar expression can be found for $R(\varphi)$, where $q_z, q_-$, $p_+$ denote the parameters analogue to $p_z, p_-$, $p_+$. Combining

<footnote>$^4$ The role played by the 4 rotation generators can be further understood as follows. As explained in Sect. 3.1, $(\mathbf{K}_3 + \mathbf{K}_4)/2$ and $(\mathbf{K}_3 - \mathbf{K}_4)/2$ generate separate rotations in the first and second sectors respectively. Thus $\mathbf{K}_3$ induces the same rotation in both sectors and $\theta_3$ can be thought of as a “coherent phase” (which is why $\mathbf{K}_3$ decouples from the other generators), while $\mathbf{K}_4$ generates opposite rotations in the two sectors and $\theta_4$ can be thought of as a “phase shift”. For the two coupling generators, $\mathbf{K}_5$ operates the same rotation in the position plane and in the momentum plane, so it can be understood as a field redefinition; while $\mathbf{K}_6$ operates a rotation that mixes positions and momenta of the two sectors.</footnote>
the above results, any element of Sp(4, $\mathbb{R}$) can be decomposed according to

$$ M = \underbrace{[\exp(p_+S_+) \cdot \exp(p_zS_z) \cdot \exp(p_-S_-) \cdot \exp(\theta_3 K_3)]}_{R(\theta)} \cdot \underbrace{[\exp(d_1K_1) \cdot \exp(d_2K_2)]}_{Z(d)} \cdot \underbrace{[\exp(q_+S_+) \cdot \exp(q_zS_z) \cdot \exp(q_-S_-) \cdot \exp(\varphi_3 K_3)]}_{R(\varphi)} \quad (3.22) $$

This fully factorised form will be of particular convenience when it comes to characterising the quantum states of the system in Sect. 5.

### 3.3 Helicity basis

In Sect. 2.1, the creation and annihilation operators have been introduced as an equivalent parametrisation of phase space, called the helicity basis. In the quantum mechanical context, it leads to the convenient occupation-number representation, which is why we now translate the above considerations into that basis.

We recall that when applying a canonical transformation, the helicity variables transform via matrices of the form $M_k = UM_kU^\dagger$ for $M_k \in \text{Sp}(4, \mathbb{R})$, which satisfy $M_k^\dagger J M_k = J$ and $\det(M_k) = 1$, see the discussion below Eq. (2.22). In particular, this is the case for the Green matrix $G_k(t, t_{\text{in}})$. Those two conditions define the SU(2, 2) group, but given that SU(2, 2) is a fifteen dimensional Lie group, Sp(4, $\mathbb{R}$) cannot be isomorphic to the whole group and instead constitutes a ten dimensional subgroup, which we denote $\mathcal{S}p(4, \mathbb{R})$. More precisely, decomposing a generic matrix $M_k \in \mathcal{M}_4(\mathbb{R})$ into blocks according to

$$ M_k = \begin{pmatrix} A & B \\ C & D \end{pmatrix} \quad (3.23) $$

with $A, B, C, D \in \mathcal{M}_2(\mathbb{R})$, the condition $M_k = UM_k(t)U^\dagger$ can be written as

$$ M_k = \begin{pmatrix} A & B \\ B^* & A^* \end{pmatrix}, \quad (3.24) $$

where $A = (1/2)[(A+D)+i(C-B)]$ and $B = (1/2)[(A-D)+i(C+B)]$. As explained below Eq. (3.1), if $M_k \in \text{Sp}(4, \mathbb{R})$ then $M_k^T \in \text{Sp}(4, \mathbb{R})$, which implies that $M_k^\dagger$ is also symplectic, hence $M_k J M_k^\dagger = J$. This leads to the two conditions

$$ AA^\dagger - BB^\dagger = I_2 \quad \text{and} \quad AB^T - BA^T = 0. \quad (3.25) $$

Upon expanding the matrices $A$ and $B$ in terms of the so-called Bogolyubov coefficients,

$$ A = \begin{pmatrix} \alpha_{11} & \alpha_{12} \\ \alpha_{21} & \alpha_{22} \end{pmatrix} \quad \text{and} \quad B = \begin{pmatrix} \beta_{11} & \beta_{12} \\ \beta_{21} & \beta_{22} \end{pmatrix}, \quad (3.26) $$

Eq. (3.25) leads to the four conditions

$$ |\alpha_{11}|^2 + |\alpha_{12}|^2 - |\beta_{11}|^2 - |\beta_{12}|^2 = 1, \quad (3.27) $$

$$ |\alpha_{21}|^2 + |\alpha_{22}|^2 - |\beta_{21}|^2 - |\beta_{22}|^2 = 1, \quad (3.28) $$

$$ \alpha_{11}\alpha_{21}^* + \alpha_{12}\alpha_{22}^* - \beta_{11}\beta_{21}^* - \beta_{12}\beta_{22}^* = 0, \quad (3.29) $$

$$ \alpha_{11}\beta_{21} + \alpha_{12}\beta_{22} - \alpha_{21}\beta_{11} - \alpha_{22}\beta_{12} = 0. \quad (3.30) $$

This fixes 2 real and 2 complex combinations out of the 8 complex Bogolyubov coefficients, that is 6 out of the 16 real parameters, and one recovers the 10 degrees of freedom of Sp(4, $\mathbb{R}$).

One may note that the Green matrix $G_k(t, t_{\text{in}})$ being an element of $\mathcal{S}p(4, \mathbb{R})$, it can also be written in terms of Bogolyubov coefficients according to Eqs. (3.24) and (3.26). In that case, the relations (3.27)–(3.30) translate the fact that the Poisson brackets (or their quantum analogue, the commutation relations) are preserved by the dynamical evolution. Moreover, when applying $\hat{a}_k(t) = G_k(t, t_{\text{in}})\hat{a}_k(t_{\text{in}})$, one notices that $A$ and $A^*$ transform annihilation into annihilation operators, and creation into creation operators, respectively. Therefore, they maintain the overall excitation number while reshuffling the excitations between the different sectors. On the contrary, $B$ and $B^*$ convert annihilation into creation operators and conversely, so they create or annihilate new excitations. We finally note that $\alpha_{11}, \alpha_{22}, \beta_{11}$ and $\beta_{22}$ act on each sector separately while $\alpha_{12}, \alpha_{21}, \beta_{12}$ and $\beta_{21}$ mix the two sectors. As a consistency check, one can verify that when these mixing Bogolyubov coefficients vanish, Eqs. (3.27)–(3.30) reduce to the Sp(2, $\mathbb{R}$)-constraint on the Bogolyubov coefficients [13], namely $|\alpha_{11}|^2 - |\beta_{11}|^2 = |\alpha_{22}|^2 - |\beta_{22}|^2 = 1$.

Let us now study the infinitesimal properties of $\mathcal{S}p(4, \mathbb{R})$, as we did for Sp(4, $\mathbb{R}$) in Sect. 3.1. We observe that the ten generators of $\mathcal{S}p(4, \mathbb{R})$ can be found by simple correspondence with the ten generators of Sp(4, $\mathbb{R}$), upon introducing $\exp(\alpha_i L_i) \equiv U \exp(\alpha_i K_i) U^\dagger = \exp(\alpha_i U K_i U^\dagger)$, where $i = 1 \cdots 10$ and where we have used the fact that $U$ is unitary. Writing $K_i = \sigma_{a_i} \otimes \sigma_{b_i}$ and $U = u \otimes I_2$ with

$$ u = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & i \\ 1 & -i \end{pmatrix}, \quad (3.31) $$

the generators of $\mathcal{S}p(4, \mathbb{R})$ can be calculated by means of Eq. (3.9) and one has $L_i = U K_i U^\dagger = (u \otimes I_2)(\sigma_{a_i} \otimes \sigma_{b_i}) (u^\dagger \otimes I_2) = u \sigma_{a_i} u^\dagger \otimes \sigma_{b_i}$. Given that $u \sigma_x u^\dagger = -\sigma_y$, $u \sigma_y u^\dagger = -\sigma_z$ and $u \sigma_z u^\dagger = \sigma_x$, one concludes that the generators of $\mathcal{S}p(4, \mathbb{R})$ are merely a reshuffling of those of Sp(4, $\mathbb{R}$), where the detailed correspondence is given in Table 2. In particular, one observes that rotations in the helic-
<table><caption>**Table 2** Generators of the helicity basis. *Sp*(4, **R**) in the fundamental representation, $L_i = UK_iU^\dagger$, where the $K_i$ generators are given in Table 1</caption><thead><tr><th>Squeezing</th><th>Rotation</th><th>Boost</th></tr></thead><tbody><tr><td>$L_1 = K_8$</td><td>$L_3 = -i K_2$</td><td>$L_7 = i K_3$</td></tr><tr><td>$L_2 = K_7$</td><td>$L_4 = -i K_1$<br>$L_5 = K_5$<br>$L_6 = -i K_9$</td><td>$L_8 = i K_4$<br>$L_9 = K_{10}$<br>$L_{10} = i K_6$</td></tr></tbody></table>

ity basis are block diagonal (which generalises the Sp(2, **R**) result where rotations in the helicity basis are (2×2) diagonal matrices [13]). Note that the commutation relations between the $L_i$ operators directly follow from those between the $K_i$ operators given in Appendix A.

The Bloch–Messiah decomposition (3.11) can also be performed in the helicity basis,

$$ \mathcal{M}(\theta, d, \varphi) = \mathcal{R}(\theta)\mathcal{Z}(d)\mathcal{R}(\varphi), \quad (3.32) $$

where $\mathcal{Z} = UZU^\dagger$ and $\mathcal{R} = URU^\dagger$ are the squeezing matrix and the rotation matrix in the helicity basis. This expression allows us to connect the Bogolyubov coefficients with the squeezing and rotation parameters. Using the above results, an explicit calculation yields

$$ \alpha_{11} = e^{-i(\theta_3+\varphi_3)} \left( \cos\theta - i \frac{\theta_4}{\theta} \sin\theta \right) \left( \cos\varphi - i \frac{\varphi_4}{\varphi} \sin\varphi \right) \cosh r_1 \\ - e^{-i(\theta_3+\varphi_3)} \left[ (\theta_5-i\theta_6) \frac{\sin\theta}{\theta} \right] \left[ (\varphi_5+i\varphi_6) \frac{\sin\varphi}{\varphi} \right] \cosh r_2, \quad (3.33) $$

$$ \alpha_{12} = e^{-i(\theta_3+\varphi_3)} \left( \cos\theta - i \frac{\theta_4}{\theta} \sin\theta \right) \left[ (\varphi_5 - i\varphi_6) \frac{\sin\varphi}{\varphi} \right] \cosh r_1 \\ + e^{-i(\theta_3+\varphi_3)} \left[ (\theta_5 - i\theta_6) \frac{\sin\theta}{\theta} \right] \left( \cos\varphi + i \frac{\varphi_4}{\varphi} \sin\varphi \right) \cosh r_2, \quad (3.34) $$

$$ \alpha_{21} = -e^{-i(\theta_3+\varphi_3)} \left[ (\theta_5+i\theta_6) \frac{\sin\theta}{\theta} \right] \left( \cos\varphi - i \frac{\varphi_4}{\varphi} \sin\varphi \right) \cosh r_1 \\ - e^{-i(\theta_3+\varphi_3)} \left( \cos\theta + i \frac{\theta_4}{\theta} \sin\theta \right) \left[ (\varphi_5+i\varphi_6) \frac{\sin\varphi}{\varphi} \right] \cosh r_2, \quad (3.35) $$

$$ \alpha_{22} = -e^{-i(\theta_3+\varphi_3)} \left[ (\theta_5 + i\theta_6) \frac{\sin\theta}{\theta} \right] \left[ (\varphi_5 - i\varphi_6) \frac{\sin\varphi}{\varphi} \right] \cosh r_1 \\ + e^{-i(\theta_3+\varphi_3)} \left( \cos\theta + i \frac{\theta_4}{\theta} \sin\theta \right) \left( \cos\varphi + i \frac{\varphi_4}{\varphi} \sin\varphi \right) \cosh r_2, \quad (3.36) $$

and

$$ \beta_{11} = e^{-i(\theta_3-\varphi_3)} \left( \cos\theta - i \frac{\theta_4}{\theta} \sin\theta \right) \left( \cos\varphi + i \frac{\varphi_4}{\varphi} \sin\varphi \right) \sinh r_1 \\ - e^{-i(\theta_3-\varphi_3)} \left[ (\theta_5 - i\theta_6) \frac{\sin\theta}{\theta} \right] \left[ (\varphi_5 - i\varphi_6) \frac{\sin\varphi}{\varphi} \right] \sinh r_2, \quad (3.37) $$

$$ \beta_{12} = e^{-i(\theta_3-\varphi_3)} \left( \cos\theta - i \frac{\theta_4}{\theta} \sin\theta \right) \left[ (\varphi_5 + i\varphi_6) \frac{\sin\varphi}{\varphi} \right] \sinh r_1 \\ + e^{-i(\theta_3-\varphi_3)} \left[ (\theta_5 - i\theta_6) \frac{\sin\theta}{\theta} \right] \left( \cos\varphi - i \frac{\varphi_4}{\varphi} \sin\varphi \right) \sinh r_2, \quad (3.38) $$

$$ \beta_{21} = -e^{-i(\theta_3-\varphi_3)} \left[ (\theta_5 + i\theta_6) \frac{\sin\theta}{\theta} \right] \left( \cos\varphi + i \frac{\varphi_4}{\varphi} \sin\varphi \right) \sinh r_1 \\ - e^{-i(\theta_3-\varphi_3)} \left( \cos\theta + i \frac{\theta_4}{\theta} \sin\theta \right) \left[ (\varphi_5 - i\varphi_6) \frac{\sin\varphi}{\varphi} \right] \sinh r_2, \quad (3.39) $$

$$ \beta_{22} = -e^{-i(\theta_3-\varphi_3)} \left[ (\theta_5 + i\theta_6) \frac{\sin\theta}{\theta} \right] \left[ (\varphi_5 + i\varphi_6) \frac{\sin\varphi}{\varphi} \right] \sinh r_1 \\ + e^{-i(\theta_3-\varphi_3)} \left( \cos\theta + i \frac{\theta_4}{\theta} \sin\theta \right) \left( \cos\varphi - i \frac{\varphi_4}{\varphi} \sin\varphi \right) \sinh r_2, \quad (3.40) $$

where $\theta = \sqrt{\theta_4^2 + \theta_5^2 + \theta_6^2}$ and $\varphi = \sqrt{\varphi_4^2 + \varphi_5^2 + \varphi_6^2}$ have already been defined in Eq. (3.21), and where we have introduced

$$ r_1 = d_1 + d_2 \quad \text{and} \quad r_2 = d_2 - d_1. \quad (3.41) $$

One can check that if the rotation parameters associated with the coupling generators vanish, i.e. if $\theta_5 = \theta_6 = \varphi_5 = \varphi_6 = 0$, then the mixing Bogolyubov coefficients vanish too, i.e. $\alpha_{12} = \alpha_{21} = \beta_{12} = \beta_{21} = 0$, and the link between the Bogolyubov coefficients and the squeezing and rotation parameters reduces to the one obtained in Sp(2, **R**) [13].

Finally, a fully factorised form for the elements of *Sp*(4, **R**) can be obtained from transposing Eq. (3.22), which gives rise to

$$ \mathcal{M} = \underbrace{\overbrace{\left[ \exp(p_+ L_+) \cdot \exp(p_z L_z) \cdot \exp(p_- L_-) \cdot \exp(\theta_3 L_3) \right]}^{\mathcal{R}(\theta)}}_{\mathcal{Z}(d)} \\ \cdot \underbrace{\overbrace{\left[ \exp(d_1 L_1) \cdot \exp(d_2 L_2) \right]}^{\mathcal{R}(\varphi)}}_{\mathcal{R}(\varphi)} \\ \cdot \underbrace{\overbrace{\left[ \exp(q_+ L_+) \cdot \exp(q_z L_z) \cdot \exp(q_- L_-) \cdot \exp(\varphi_3 L_3) \right]}^{\mathcal{R}(\varphi)}}_{\mathcal{R}(\varphi)}, \quad (3.42) $$

where $L_z = US_zU^\dagger = L_4/(2i)$, $L_+ = US_+U^\dagger = (L_5 - iL_6)/(2i)$ and $L_- = US_-U^\dagger = (L_5 + iL_6)/(2i)$.

### 3.4 Quantum representation

Since the Green's matrix is an element of the symplectic group, in order to describe the dynamics in the occupation-number representation, one first needs to derive the quantum representation of the elements of *Sp*(4, **R**). This can be done by following the procedure outlined in Refs. [40–42] and presented in details in Appendix B of Ref. [13]. It consists
in first linearising the elements of the *Sp*(4, **R**) Lie group,

$$ \mathcal{M}_k \approx \mathcal{I}_4 + \sum_{a=1}^{10} \epsilon_k^a \mathcal{L}_a, \qquad (3.43) $$

where $\epsilon_k^a$ are (small) real parameters. This generates a canonical transformation on the helicity variables, given by

$$ \hat{\mathbf{a}}'_{\vec{k}} \approx \hat{\mathbf{a}}_{\vec{k}} + \sum_{a=1}^{10} \epsilon_k^a \mathbf{L}_a \hat{\mathbf{a}}_{\vec{k}}. \qquad (3.44) $$

Our goal is to find a set of operators $\hat{\mathcal{L}}_a^{\vec{k}}$ such that this can be written as a unitary transformation $\hat{\mathbf{a}}'_{\vec{k}} = \hat{\mathcal{M}}^{\dagger} \hat{\mathbf{a}}_{\vec{k}} \hat{\mathcal{M}}$, with

$$ \hat{\mathcal{M}} \approx \hat{\mathbb{I}} + \int_{\mathbb{R}^3} d^3 \vec{q} \sum_{a=1}^{10} \epsilon_{\vec{q}}^a \hat{\mathcal{L}}_a^{\vec{q}}. \qquad (3.45) $$

Note that $\hat{\mathcal{L}}_a^{\vec{q}}$ has to be anti-hermitian for $\hat{\mathcal{M}}$ to be unitary. By expanding this transformation in $\epsilon_{\vec{q}}^a$, one obtains $\hat{\mathbf{a}}'_{\vec{k}} \approx \hat{\mathbf{a}}_{\vec{k}} + [\hat{\mathbf{a}}_{\vec{k}}, \int_{\mathbb{R}^3} d^3 q \sum_{a=1}^{10} \epsilon_{\vec{q}}^a \hat{\mathcal{L}}_a^{\vec{q}}] (\text{where the commutator is performed on each entry of } \hat{\mathbf{a}}_{\vec{k}})$. Since this formula should match Eq. (3.44), and given the commutation relations (2.16), this suggests to look for $\hat{\mathcal{L}}_a^{\vec{q}}$ operators that are quadratic in the creation and annihilation operators,

$$ \hat{\mathcal{L}}_a^{\vec{q}} = \hat{\mathbf{a}}_{\vec{q}}^{\dagger} \mathcal{Q}^a \hat{\mathbf{a}}_{\vec{q}}, \qquad (3.46) $$

where $\mathcal{Q}^a$ is anti-hermitian since $\hat{\mathcal{L}}_a^{\vec{q}}$ is. Making use of Eq. (2.16), this leads to $[\hat{\mathbf{a}}_{\vec{k}}, \int_{\mathbb{R}^3} d^3 q \sum_{a=1}^{10} \epsilon_{\vec{q}}^a \hat{\mathcal{L}}_a^{\vec{q}}] = i \mathcal{J} \sum_{a=1}^{10} \epsilon_k^a \mathcal{Q}^a \hat{\mathbf{a}}_{\vec{k}}$. By identification with Eq. (3.44), this gives $\mathcal{L}_a = i \mathcal{J} \mathcal{Q}^a$, which can be inverted as

$$ \mathcal{Q}^a = i \mathcal{J} \mathcal{L}_a. \qquad (3.47) $$

This yields the generators $\hat{\mathcal{L}}_a^{\vec{q}}$ of the quantum representation listed in Table 3.

Another set of generators that is obtained from straightforward linear combinations of the $\hat{\mathcal{L}}_a^{\vec{k}}$ operators is given by

$$ \hat{a}_{i, \vec{k}}^{\dagger} \hat{a}_{i, -\vec{k}}^{\dagger}, \quad \hat{a}_{i, \vec{k}} \hat{a}_{i, -\vec{k}}, \quad \hat{a}_{i, \vec{k}}^{\dagger} \hat{a}_{i, \vec{k}} + \hat{a}_{i, -\vec{k}}^{\dagger} \hat{a}_{i, -\vec{k}} + 1 \qquad (3.48) $$

with $i = 1, 2$,

$$ \hat{a}_{i, \vec{k}}^{\dagger} \hat{a}_{j, -\vec{k}}^{\dagger} + \hat{a}_{j, \vec{k}}^{\dagger} \hat{a}_{i, -\vec{k}}^{\dagger}, \quad \hat{a}_{i, \vec{k}} \hat{a}_{j, -\vec{k}} + \hat{a}_{j, \vec{k}} \hat{a}_{i, -\vec{k}}, \\ \hat{a}_{i, \vec{k}}^{\dagger} \hat{a}_{j, \vec{k}} + \hat{a}_{i, -\vec{k}}^{\dagger} \hat{a}_{j, -\vec{k}} \quad \text{with } i, j = 1, 2, \qquad (3.49) $$

which correspond to the coupling generators. The operators in Eqs. (3.48) generate two-mode creation, two-mode annihilation and a number counting operation respectively, acting on each sector separately. The first and second mixing generators in Eq. (3.49) originate from the two mixing boosts ($\hat{\mathcal{L}}_9^{\vec{k}}$

<table><caption>**Table 3** Quantum representation of the generators of *Sp*(4, **R**)</caption><tbody><tr><td>Squeezing</td><td>$$ \begin{aligned} \hat{\mathcal{L}}_1^{\vec{k}} &= (\hat{a}_{1, \vec{k}}^{\dagger} \hat{a}_{1, -\vec{k}}^{\dagger} - \hat{a}_{1, \vec{k}} \hat{a}_{1, -\vec{k}}) \\ &\quad - (\hat{a}_{2, \vec{k}}^{\dagger} \hat{a}_{2, -\vec{k}}^{\dagger} - \hat{a}_{2, \vec{k}} \hat{a}_{2, -\vec{k}}) \\ \hat{\mathcal{L}}_2^{\vec{k}} &= (\hat{a}_{1, \vec{k}}^{\dagger} \hat{a}_{1, -\vec{k}}^{\dagger} - \hat{a}_{1, \vec{k}} \hat{a}_{1, -\vec{k}}) \\ &\quad + (\hat{a}_{2, \vec{k}}^{\dagger} \hat{a}_{2, -\vec{k}}^{\dagger} - \hat{a}_{2, \vec{k}} \hat{a}_{2, -\vec{k}}) \\ \hat{\mathcal{L}}_3^{\vec{k}} &= -i (\hat{a}_{1, \vec{k}}^{\dagger} \hat{a}_{1, \vec{k}} + \hat{a}_{1, -\vec{k}}^{\dagger} \hat{a}_{1, -\vec{k}} + 1) \\ &\quad -i (\hat{a}_{2, \vec{k}}^{\dagger} \hat{a}_{2, \vec{k}} + \hat{a}_{2, -\vec{k}}^{\dagger} \hat{a}_{2, -\vec{k}} + 1) \\ \hat{\mathcal{L}}_4^{\vec{k}} &= -i (\hat{a}_{1, \vec{k}}^{\dagger} \hat{a}_{1, \vec{k}} + \hat{a}_{1, -\vec{k}}^{\dagger} \hat{a}_{1, -\vec{k}} + 1) \\ &\quad +i (\hat{a}_{2, \vec{k}}^{\dagger} \hat{a}_{2, \vec{k}} + \hat{a}_{2, -\vec{k}}^{\dagger} \hat{a}_{2, -\vec{k}} + 1) \\ \hat{\mathcal{L}}_5^{\vec{k}} &= (\hat{a}_{1, \vec{k}}^{\dagger} \hat{a}_{2, \vec{k}} + \hat{a}_{1, -\vec{k}}^{\dagger} \hat{a}_{2, -\vec{k}}) \\ &\quad - (\hat{a}_{2, \vec{k}}^{\dagger} \hat{a}_{1, \vec{k}} + \hat{a}_{2, -\vec{k}}^{\dagger} \hat{a}_{1, -\vec{k}}) \\ \hat{\mathcal{L}}_6^{\vec{k}} &= -i (\hat{a}_{1, \vec{k}}^{\dagger} \hat{a}_{2, \vec{k}} + \hat{a}_{1, -\vec{k}}^{\dagger} \hat{a}_{2, -\vec{k}}) \\ &\quad -i (\hat{a}_{2, \vec{k}}^{\dagger} \hat{a}_{1, \vec{k}} + \hat{a}_{2, -\vec{k}}^{\dagger} \hat{a}_{1, -\vec{k}}) \\ \hat{\mathcal{L}}_7^{\vec{k}} &= i (\hat{a}_{1, \vec{k}}^{\dagger} \hat{a}_{1, -\vec{k}} + \hat{a}_{1, \vec{k}} \hat{a}_{1, -\vec{k}}) \\ &\quad +i (\hat{a}_{2, \vec{k}}^{\dagger} \hat{a}_{2, -\vec{k}} + \hat{a}_{2, \vec{k}} \hat{a}_{2, -\vec{k}}) \\ \hat{\mathcal{L}}_8^{\vec{k}} &= i (\hat{a}_{1, \vec{k}}^{\dagger} \hat{a}_{1, -\vec{k}}^{\dagger} + \hat{a}_{1, \vec{k}} \hat{a}_{1, -\vec{k}}) \\ &\quad -i (\hat{a}_{2, \vec{k}}^{\dagger} \hat{a}_{2, -\vec{k}}^{\dagger} + \hat{a}_{2, \vec{k}} \hat{a}_{2, -\vec{k}}) \\ \hat{\mathcal{L}}_9^{\vec{k}} &= (\hat{a}_{1, \vec{k}}^{\dagger} \hat{a}_{2, -\vec{k}}^{\dagger} + \hat{a}_{2, \vec{k}}^{\dagger} \hat{a}_{1, -\vec{k}}^{\dagger}) \\ &\quad - (\hat{a}_{1, \vec{k}} \hat{a}_{2, -\vec{k}} + \hat{a}_{2, \vec{k}} \hat{a}_{1, -\vec{k}}) \\ \hat{\mathcal{L}}_{10}^{\vec{k}} &= i (\hat{a}_{1, \vec{k}}^{\dagger} \hat{a}_{2, -\vec{k}}^{\dagger} + \hat{a}_{2, \vec{k}}^{\dagger} \hat{a}_{1, -\vec{k}}^{\dagger}) \\ &\quad +i (\hat{a}_{1, \vec{k}} \hat{a}_{2, -\vec{k}} + \hat{a}_{2, \vec{k}} \hat{a}_{1, -\vec{k}}) \end{aligned} $$</td></tr><tr><td>Rotation</td><td></td></tr><tr><td>Boost</td><td></td></tr></tbody></table>

and $\hat{\mathcal{L}}_{10}^{\vec{k}}$). They correspond to the creation and annihilation of entangled pairs of particles with opposite momenta in the two sectors. The third mixing generator in Eq. (3.49) originates from the two mixing rotations ($\hat{\mathcal{L}}_5^{\vec{k}}$ and $\hat{\mathcal{L}}_6^{\vec{k}}$). It does not lead to net particle creation, but rather transfer excitations from one sector to the other. One can check that all these operations preserve momentum conservation.

## 4 Quantum dynamics

In Sect. 2, we explained how the symplectic group Sp(4, **R**) naturally appears in the description of physical systems made of two free fields. Having then studied its mathematical structure in Sect. 3, we now want to apply these tools to describe the dynamics of two-field systems. To this end, we explore the structure of the system's Hamiltonian in Sect. 4.1, before applying our findings to an explicit example. Then, we provide in Sect. 4.2 a tractable expression of the evolution operator.
## 4.1 Hamiltonian structure

In this section, we express the Hamiltonian in the occupation-number representation, to see which interactions are allowed in the theory on generic grounds. The only constraint we have on the Hamiltonian is that it stems from a real symmetric kernel, i.e. $H_k$ is a real symmetric matrix. From Eq. (2.15), $\mathcal{H}_k$ can thus be written as $\mathcal{H}_k = UN_kU^\dagger$, where $N_k = (D_k^{-1})^\mathrm{T}H_kD_k^{-1} + (D_k^{-1})^\mathrm{T}\Omega\dot{D}_k^{-1}$ is a real symmetric matrix (the first term is obviously symmetric, while for the second term, this can be shown by differentiating with respect to time the symplectic relation (2.19) satisfied by $D_k^{-1}$, and by using that $\Omega^\mathrm{T} = -\Omega$), which we parametrise as

$$N_k = \begin{pmatrix} A & C \\ C^\mathrm{T} & B \end{pmatrix}, \qquad (4.1)$$

with $A, B, C \in \mathcal{M}_2(\mathbb{R})$ and $A^\mathrm{T} = A, B^\mathrm{T} = B$. This leads to

$$\mathcal{H}_k = \begin{pmatrix} a & b \\ b^* & a^* \end{pmatrix} \qquad (4.2)$$

with $a = (A+B)+i(C^\mathrm{T}-C)$ and $b = (A-B)+i(C^\mathrm{T}+C)$. Inverting those formulas, one obtains that $a+a^* = 2(A+B)$ and $a-a^* = 2i(C^\mathrm{T}-C)$, so the real part of $a$ should be symmetric while its imaginary part should be antisymmetric; and $b+b^* = 2(A-B)$ and $b-b^* = 2i(C^\mathrm{T}+C)$, so both the real and imaginary parts of $b$ should be symmetric. This imposes that $a$ and $b$ are of the form

$$\begin{aligned} a &= \begin{pmatrix} F_{1,k} & F_{1 \leftrightarrow 2,k} e^{-i\varphi_k} \\ F_{1 \leftrightarrow 2,k} e^{-i\varphi_k} & F_{2,k} \end{pmatrix}, \\ b &= \begin{pmatrix} R_{1,k} e^{i\Theta_{1,k}} & R_{1 \leftrightarrow 2,k} e^{i\xi_k} \\ R_{1 \leftrightarrow 2,k} e^{i\xi_k} & R_{2,k} e^{i\Theta_{2,k}} \end{pmatrix}, \end{aligned} \qquad (4.3)$$

where $F_{1,k}, F_{1 \leftrightarrow 2,k}, F_{2,k}, \varphi_k, R_{1,k}, R_{1 \leftrightarrow 2,k}, R_{2,k}, \Theta_{1,k}, \Theta_{2,k}$ and $\xi_k$ are ten real parameters, that in general depend on time, and which indeed saturate the ten degrees of freedom contained in $N_k$. In fact, since the Green matrix is generated from the (integrated) Hamiltonian, those ten parameters fully determine the squeezing and rotation parameters, or equivalently the Bogolyubov coefficients, that characterise the dynamics.

The Hamiltonian $\hat{\mathcal{H}} = \int_{\mathbb{R}^{3+}} d^3\vec{k} \hat{a}_{\vec{k}}^\dagger \mathcal{H}_k \hat{a}_{\vec{k}}$ thus contain three terms,

$$\hat{\mathcal{H}} = \int_{\mathbb{R}^{3+}} d^3\vec{k} (\hat{\mathcal{H}}_{1,\vec{k}} + \hat{\mathcal{H}}_{2,\vec{k}} + \hat{\mathcal{H}}_{1 \leftrightarrow 2,\vec{k}}), \qquad (4.4)$$

with

$$\begin{aligned} \hat{\mathcal{H}}_{i,\vec{k}} = F_{i,k} & \begin{pmatrix} \hat{a}_{i,\vec{k}}^{\dagger} \hat{a}_{i,\vec{k}} + \hat{a}_{i,-\vec{k}}^{\dagger} \hat{a}_{i,-\vec{k}} + 1 \\ + R_{i,k} \begin{pmatrix} e^{i\Theta_{i,k}} \hat{a}_{i,\vec{k}}^{\dagger} \hat{a}_{i,-\vec{k}}^{\dagger} + \text{h.c.} \end{pmatrix} \end{aligned} \quad \text{for } i=1,2, \qquad (4.5)$$

$$\begin{aligned} \hat{\mathcal{H}}_{1 \leftrightarrow 2, \vec{k}} = F_{1 \leftrightarrow 2, k} e^{i \varphi_k} & \begin{pmatrix} \hat{a}_{1, \vec{k}}^{\dagger} \hat{a}_{2, \vec{k}} + \hat{a}_{1, -\vec{k}}^{\dagger} \hat{a}_{2, -\vec{k}} \\ + R_{1 \leftrightarrow 2, k} e^{i \xi_k} & \begin{pmatrix} \hat{a}_{1, \vec{k}}^{\dagger} \hat{a}_{2, -\vec{k}}^{\dagger} + \hat{a}_{2, \vec{k}}^{\dagger} \hat{a}_{1, -\vec{k}}^{\dagger} \end{pmatrix} \end{aligned} + \text{h.c.} \qquad (4.6)$$

The Hamiltonian can thus be decomposed onto the quantum generators listed in Table 3 (this is because $\mathcal{J}\mathcal{H}_k$ is an element of the Lie algebra), hence we can benefit from their physical interpretation discussed in Sect. 3.4. The components $\hat{\mathcal{H}}_1$ and $\hat{\mathcal{H}}_2$ drive each sector separately and are made of two terms: a harmonic part controlled by $F_{i,k}$ that does not induce particle creation, and a parametric part controlled by $R_{i,k}$ that changes the particle content. Physically, this parametric amplification is related to the presence of an external field (the electric field for the Schwinger effect [43], the gravitational field for the Hawking effect, space-time curvature for the physics of cosmological perturbations, etc.). The component $\hat{\mathcal{H}}_{1 \leftrightarrow 2}$ is an interaction term between the two sectors, and also contains two contributions: a “transferring” part, controlled by $F_{1 \leftrightarrow 2, k}(t)$ and built from the two mixing rotation operators, that transfers particles from one sector to the other; and an “entangling” part, controlled by $R_{1 \leftrightarrow 2, k}$ and built from the two mixing boost operators, that creates or annihilates joint pairs of particles in the two sectors.

Let us note that even when the entangling part is absent, i.e. when $R_{1 \leftrightarrow 2, k} = 0$, entanglement between the two sectors can still indirectly arise from first creating particles in the two sectors separately with the parametric terms, and then transferring those particles between sectors by means of the transferring term.

### Example of two massless fields in a cosmological background

Let us illustrate the formalism introduced above with a simple example. We consider two massless scalar fields $\phi$ and $\chi$ on a Friedmann–Lemaître–Robertson–Walker spatially flat geometry,

$$ds^2 = a^2(\eta) (-d\eta^2 + \delta^{ij} dx_i dx_j), \qquad (4.7)$$

where $a$ is the scale factor and $\eta$ is the conformal time. The Ricci scalar $R$ of this metric is given by $R = 6a''/a^3$, where a prime denotes derivation with respect to conformal time. The action is given by

$$S = - \int d^4 x \sqrt{-\det g} \left( \frac{1}{2} g^{\mu\nu} \partial_{\mu} \phi \partial_{\nu} \phi + \frac{R}{2} \zeta \phi^2 + \frac{1}{2} g^{\mu\nu} \partial_{\mu} \chi \partial_{\nu} \chi + \frac{R}{2} \zeta \chi^2 + \lambda^2 \phi \chi \right), \qquad (4.8)$$

where $\zeta$ is the conformal coupling constant and $\lambda$ is a coupling parameter that has dimension of a mass. Expanding the scalar fields into Fourier modes as in Eq. (2.6) and making use of the reality prescription $\phi_{\vec{k}}^* = \phi_{-\vec{k}}$ and $\chi_{\vec{k}}^* = \chi_{-\vec{k}}$, with
the metric (4.7) the action reads

$$ S = \int d\eta L = - \int d\eta \int_{\mathbb{R}^{3+}} d^3\vec{k} \left[ a^2 \phi'_{-\vec{k}} \phi'_{\vec{k}} + (k^2 a^2 + R \zeta a^4) \phi_{-\vec{k}} \phi_{\vec{k}} + a^2 \chi'_{-\vec{k}} \chi'_{\vec{k}} + (k^2 a^2 + R \zeta a^4) \chi_{-\vec{k}} \chi_{\vec{k}} + \lambda^2 a^4 (\phi_{-\vec{k}} \chi_{\vec{k}} + \chi_{-\vec{k}} \phi_{\vec{k}}) \right], \quad (4.9) $$

which also defines the Lagrangian density $L$. The conjugate momenta can be identified as $p_{\vec{k}}^{\phi} = a^2 \phi'_{-\vec{k}}$ and $p_{\vec{k}}^{\chi} = a^2 \chi'_{-\vec{k}}$, and upon performing a Legendre transform, one obtains a Hamiltonian of the form (2.7), with

$$ \mathbf{H}_k = \begin{pmatrix} k^2 a^2 + R \zeta a^4 & \lambda^2 a^4 & 0 & 0 \\ \lambda^2 a^4 & k^2 a^2 + R \zeta a^4 & 0 & 0 \\ 0 & 0 & 1/a^2 & 0 \\ 0 & 0 & 0 & 1/a^2 \end{pmatrix} \quad \text{and} \\ z_{\vec{k}} = \begin{pmatrix} \phi_{\vec{k}} \\ \chi_{\vec{k}}^{\phi} \\ p_{\vec{k}}^{\chi} \\ p_{\vec{k}}^{\phi} \end{pmatrix}. \quad (4.10) $$

The momentum sector of the Hamiltonian (i.e. the bottom-right block) can be simplified by performing a canonical transformation of the form (2.20), $\tilde{z}_{\vec{k}} = M_k z_{\vec{k}}$, with

$$ M_k = \begin{pmatrix} a & 0 & 0 & 0 \\ 0 & a & 0 & 0 \\ a' & 0 & 1/a & 0 \\ 0 & a' & 0 & 1/a \end{pmatrix}. \quad (4.11) $$

One can check that $M_k$ satisfies the symplectic relation (2.19), and that the new Hamiltonian kernel, given by Eq. (2.21), reads

$$ \tilde{\mathbf{H}}_k = \begin{pmatrix} k^2 + R \zeta a^2 - a''/a & \lambda^2 a^2 & 0 & 0 \\ \lambda^2 a^2 & k^2 + R \zeta a^2 - a''/a & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}. \quad (4.12) $$

The $\tilde{z}_{\vec{k}}$ variables are usually referred to as the Mukhanov-Sasaki variables, and in the helicity basis, Eq. (2.15) gives rise to a Hamiltonian kernel of the form (4.2), with

$$ k R_{i,k} = \frac{a^2}{2} R \zeta - \frac{a''}{2a}, \quad k F_{i,k} = k^2 + k R_{i,k}, \quad (4.13) \\ F_{1 \leftrightarrow 2,k} = R_{1 \leftrightarrow 2,k} = \frac{\lambda^2 a^2}{2k}, $$

and where all the phases vanish, i.e. $\varphi_k = \Theta_{1,k} = \Theta_{2,k} = \xi_k = 0$. This allows one to relate the four Hamiltonian contributions identified in Sect. 4.1, namely the harmonic, parametric, transferring and entangling terms, to the microphysical

parameters of the problem. In particular, one can see that the parametric term is generated by the external “field” $a(t)$ (since when the scale factor is a constant, $R_{i,k} = 0$), and that the transferring and entangling terms are controlled by the coupling parameter $\lambda$, in agreement with the discussion in Sect. 4.1.

Since $R = 6a''/a^3$, the parametric term is given by $kR_{i,k} = (3\zeta - 1/2)a''/a$. Therefore, if one sets $\zeta = 1/6$, the parametric term vanishes, and the harmonic term simply becomes $F_{i,k} = k$, as in flat space time. This is because, in that case, the fields are conformally coupled to the metric, which removes the effect of space-time expansion that otherwise generates parametric amplification. This implies that, in the absence of direct coupling (i.e. if $\lambda = 0$), there is no particle creation, and the above system is made of uncoupled harmonic oscillators. However, in the presence of direct coupling, both the transferring and the entangling terms become non vanishing, which guarantees that entangled pairs of particles are created and exchanged between the two fields. This is because the interaction term breaks conformal invariance, as noticed in Ref. [44].

## 4.2 Evolution operator

We now investigate the integrated dynamics. As explained in Sect. 3.3, the Green’s matrix in the helicity basis, $G_k(t, t_{in})$, belongs to $Sp(4, \mathbb{R})$. Its quantum analogue, the so-called evolution operator $\hat{U}_{\vec{k}}(t, t_{in})$, thus lies in the quantum representation of $Sp(4, \mathbb{R})$, which was studied in Sect. 3.4. In this section, we make use of the formal results derived above to derive a tractable expression for the evolution operator.

Since the Green matrix is separable in Fourier space for free fields, the evolution operator is also factorisable as

$$ \hat{U}(t, t_{in}) = \prod_{\vec{k} \in \mathbb{R}^{3+}} \hat{U}_{\vec{k}}(t, t_{in}), \quad (4.14) $$

where $\hat{U}_{\vec{k}}(t, t_{in}) \in Sp(4, \mathbb{R})$. Making use of the Bloch-Messiah decomposition presented in Sects. 3.2 and 3.3, one can write

$$ \hat{U}_{\vec{k}}(t, t_{in}) = \hat{R}_{\vec{k}}(\theta_k) \cdot \hat{Z}_{\vec{k}}(d_k) \cdot \hat{R}_{\vec{k}}(\varphi_k), \quad (4.15) $$

see Eq. (3.32), where the squeezing and rotation parameters depend only on the norm of the wavevector, since this is also the case for the Green’s matrix. Further factorisation can be obtained by the procedure outlined in Sect. 3.3 and leading to Eq. (3.42), and replacing the generators in Eq. (3.42) by their expression in the quantum representation given in Table 3, one finds that the evolution operator involves three types of operation only, namely

$$ \hat{R}_{i,\vec{k}}(\theta_k) = \exp \left[ -i\theta_k \left( \hat{a}_{i,\vec{k}}^{\dagger} \hat{a}_{i,\vec{k}} + \hat{a}_{i,-\vec{k}}^{\dagger} \hat{a}_{i,-\vec{k}} + 1 \right) \right], \quad (4.16) $$
$$ \hat{Z}_{i,\vec{k}}(r_k) = \exp \left[ r_k \left( \hat{a}_{i,k}^{\dagger} \hat{a}_{i,-\vec{k}}^{\dagger} - \hat{a}_{i,\vec{k}} \hat{a}_{i,-\vec{k}} \right) \right], \quad (4.17) $$

$$ \hat{R}_{i \to j, \vec{k}}(p_k) = \exp \left[ i p_k \left( \hat{a}_{j,\vec{k}}^{\dagger} \hat{a}_{i,\vec{k}} + \hat{a}_{j,-\vec{k}}^{\dagger} \hat{a}_{i,-\vec{k}} \right) \right]. \quad (4.18) $$

The operators $\hat{R}_{i,\vec{k}}(\theta_k)$ are constructed from the rotation generators and induce global phase shifts in sector $i$, without changing the particle content. The squeezing operators $\hat{Z}_{i,\vec{k}}(r_k)$ create pairs of entangled particles in each sector separately, and $\hat{R}_{i \to j, \vec{k}}(p_k)$ transfers particles from one sector to the other without changing their overall number. In order to express the operators $\hat{R}_{\vec{k}}(\theta_k)$, $\hat{R}_{\vec{k}}(\varphi_k)$ and $\hat{Z}_{\vec{k}}(d_k)$ that appear in Eq. (4.15) in terms of $\hat{R}_{i,\vec{k}}$, $\hat{R}_{i \to j, \vec{k}}$ and $\hat{Z}_{i,\vec{k}}$, we adopt a diagrammatic representation analogous to quantum circuits:

## 5 Quantum state

We are now in a position to write down the quantum state of the system in the occupation-number basis. To that hand, we first need to equip our Hilbert space with a Fock-space structure. In practice, we impose that the creation and annihilation operators are ladder operators for the Hamiltonian at initial time $t_{in}$, i.e.

$$ [\hat{H}, \hat{a}_{i, \pm \vec{k}}^{\dagger}] = c_{i,k} \hat{a}_{i, \pm \vec{k}}^{\dagger} \quad \text{and} \quad [\hat{H}, \hat{a}_{i, \pm \vec{k}}] = -c_{i,k} \hat{a}_{i, \pm \vec{k}} \quad (5.1) $$

at $t_{in}$, for $i = 1, 2$ and where $c_{i,k}$ are real parameters. Among the various terms (4.5) and (4.6) that the Hamiltonian may contain, only number counting operators, $\hat{N}_{i, \pm \vec{k}} \equiv \hat{a}_{i, \pm \vec{k}}^{\dagger} \hat{a}_{i, \pm \vec{k}}$, give commutators of the form (5.1). This imposes

<figure><img src="image_3.png" /></figure>

In those graphical representations, the top line stands for operations performed on the first sector, the bottom line on the second sector, and entangling operations are displayed with joint boxes. To run a “circuit”, one successively applies the operations from the left to right (along the direction shown with the arrows). The parameters entering the circuits were introduced in Sect. 3.2. The squeezing parameters $r_1^k = d_1^k + d_2^k$ and $r_2^k = d_2^k - d_1^k$ control the two-mode creation in each sector, while the mixing parameters $p_z^k, p_\perp^k, p_+^k$ and $q_z^k, q_\perp^k, q_+^k$ control the entanglement between the two sectors. This shows that dynamical evolution can be seen as successive applications of phase rotations within each sector, creations of particles with opposite momenta in each subspace, and particle transfers between the two sectors.

that only the harmonic terms, i.e. those controlled by $F_{i,k}$, can be present at initial time, leading to $c_{i,k} = F_{i,k}(t_{in})$; while one must have $R_{i,k}(t_{in}) = R_{1 \leftrightarrow 2,k}(t_{in}) = F_{1 \leftrightarrow 2,k}(t_{in}) = 0$. This is for instance the case in the example discussed in Sect. 4.1 if the expansion is initially accelerating, since one can check that, then, in the asymptotic past (i.e. in the limit $a \to 0$), only $F_{i,k} \simeq k$ survives in Eq. (4.13). More generally, this is true in inflating backgrounds where Fourier modes get blue-shifted below the Hubble scale at early time, and hereafter we will assume that this condition is indeed satisfied.

We then build the Fock basis of the initial Hilbert space,

$$ \mathcal{E}(t_{in}) = \prod_{\vec{k} \in \mathbb{R}^{3+}} \mathcal{E}_{\vec{k}}^{(1)}(t_{in}) \otimes \mathcal{E}_{-\vec{k}}^{(1)}(t_{in}) \otimes \mathcal{E}_{\vec{k}}^{(2)}(t_{in}) $$
$$ \otimes \mathcal{E}_{-\vec{k}}^{(2)}(t_{\text{in}}), \qquad (5.2) $$

which is a quadripartite system. The four-mode vacuum state is the one annihilated by all four annihilation operators, denoted by

$$ |\Phi(t_{\text{in}})\rangle = \prod_{\vec{k} \in \mathbb{R}^{3+}} |0_{\vec{k}}^{(1)}, 0_{-\vec{k}}^{(1)}, 0_{\vec{k}}^{(2)}, 0_{-\vec{k}}^{(2)}\rangle(t_{\text{in}}), \qquad (5.3) $$

and the rest of the Fock space can be built by successive applications of creation operators, which leads to the Fock states

$$ \begin{aligned} |m_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}, s_{\vec{k}}^{(2)}, t_{-\vec{k}}^{(2)}\rangle(t_{\text{in}}) \\ = \frac{(\hat{a}_{1,\vec{k}}^{\dagger})^{m}}{m!} \frac{(\hat{a}_{1,-\vec{k}}^{\dagger})^{n}}{n!} \frac{(\hat{a}_{2,\vec{k}}^{\dagger})^{s}}{s!} \frac{(\hat{a}_{2,-\vec{k}}^{\dagger})^{t}}{t!} \\ |0_{\vec{k}}^{(1)}, 0_{-\vec{k}}^{(1)}, 0_{\vec{k}}^{(2)}, 0_{-\vec{k}}^{(2)}\rangle(t_{\text{in}}). \end{aligned} \qquad (5.4) $$

From now on, we drop the argument $t_{\text{in}}$ to make the notation lighter.

## 5.1 Four-mode squeezed state

The evolved vacuum is obtained by application of the evolution operator on the initial vacuum<fnref n="5" />

$$ |\Phi(t)\rangle = \hat{U}(t, t_{\text{in}})|\Phi(t_{\text{in}})\rangle. \qquad (5.5) $$

In this section, we present its explicit expression in the Fock space, where we derive the most generic form of a four-mode squeezed state.

In Appendix B, we apply the operators appearing in the circuits sketched in Sect. 4.2 onto the vacuum state (5.3) one after the other, and derive the following expression for the evolved vacuum state $|\Phi(t)\rangle = \prod_{\vec{k} \in \mathbb{R}^{3+}} |\Phi_{\vec{k}}(t)\rangle$,

$$ |\Phi_{\vec{k}}(t)\rangle = \sum_{n,m=0}^{\infty} \sum_{s,t=-n}^{m} c_k(n, m, s, t) \\ |(n+s)_{\vec{k}}^{(1)}, (n+t)_{-\vec{k}}^{(1)}, (m-s)_{\vec{k}}^{(2)}, (m-t)_{-\vec{k}}^{(2)}\rangle, \qquad (5.6) $$

with

$$ \begin{aligned} c_k(n, m, s, t) = & \frac{e^{-2i[\varphi_3^k(n+m+1)+\varphi_3^k]}}{\cosh r_1^k \cosh r_2^k} (-1)^{n+m} e^{p_z^k(m-n)} \\ & \times \tanh^n(r_1^k) \tanh^m(r_2^k) (ip_+^k)^{-s-t} \frac{m!}{n!} \\ & \times \sqrt{\frac{(m-s)!(m-t)!}{(n+s)!(n+t)!}} \\ & \times \sum_{i=\max(0,s)}^{m} \frac{\left(p_-^k p_+^k e^{-p_z^k}\right)^i (n+i)!}{i!(i-s)!(m-i)!} \\ & \times \sum_{j=\max(0,t)}^{m} \frac{\left(p_-^k p_+^k e^{-p_z^k}\right)^j (n+j)!}{j!(j-t)!(m-j)!} \end{aligned} \qquad (5.7) $$

(and where the summation index $t$ should not be confused with the time parameter). This generic expression of four-mode squeezed states expanded in the Fock basis is one of the main results of this paper. It features an infinite tower of entangled states, characterised by non-trivial numbers of excitations. More precisely, there are 6 indices being summed over, which can be interpreted as follows. The indices $i$ and $j$ are internal and can be resumed within each Fock state separately, their physical interpretation will be made clearer below.<fnref n="6" /> The indices $n$ and $m$ are related to the creation of entangled pairs of excitations with opposite wavevector inside each sector separately, and arise from the squeezing operation (4.17). The indices $s$ and $t$ are related to the transfer of excitations between the two sectors, with wavevector $\vec{k}$ and $-\vec{k}$ respectively. They arise from the transferring operation (4.18) and entangle the two sectors. They are negative when excitations go from the first sector to the second sector, positive otherwise, and their bounds guarantee that the num-

$$ \begin{aligned} c_k(n, m, s, t) = & \frac{e^{-2i[\varphi_3^k + \varphi_3^k(n+m+1)]}}{\cosh r_1^k \cosh r_2^k} (-1)^{n+m} e^{p_z^k(m-n)} \\ & \times \tanh^n(r_1^k) \tanh^m(r_2^k) (p_-^k p_+^k e^{-p_z^k})^{\max(0,s)+\max(0,t)} \\ & \times (ip_+^k)^{-s-t} \frac{m!}{n!} \sqrt{\frac{(m-s)!(m-t)!}{(n+s)!(n+t)!}} \frac{1}{|s|!|t|!} \\ & \times \frac{[n+\max(0,s)]![n+\max(0,t)]!}{[m-\max(0,s)]![m-\max(0,t)]!} \\ & {}_2F_1[-m+\max(0,s), 1+n+\max(0,s), 1-s+2\max(0,s), -p_-^k p_+^k e^{-p_z^k}] \\ & {}_2F_1[-m+\max(0,t), 1+n+\max(0,t), 1-t+2\max(0,t), -p_-^k p_+^k e^{-p_z^k}]. \end{aligned} \qquad (5.8) $$

<footnote n="5">⁵ Note that, in place of the initial vacuum, one may consider an initial state of the form $|\Psi(t_{\text{in}})\rangle = \prod_{\vec{k}} \hat{A}_{\vec{k}}(t_{\text{in}})|\Phi(t_{\text{in}})\rangle$, where $\hat{A}_{\vec{k}}(t_{\text{in}}) \in \mathcal{S}p(4, \mathbb{R})$. Such states are often referred to as “alpha vacua” [45] (where further restrictions on $\hat{A}_{\vec{k}}(t_{\text{in}})$ sometimes apply). Since the product of two elements of $\mathcal{S}p(4, \mathbb{R})$ still lies in $\mathcal{S}p(4, \mathbb{R})$, the evolved alpha vacua are still of the form (5.5), so our discussion encompasses this possibility.</footnote>

<footnote n="6">⁶ The internal indices $i$ and $j$ can be resumed in terms of the hypergeometric function $_2F_1$, leading to</footnote>
ber of excitations remains non negative. One should also note that $c_k(m, n, s, t)$ is invariant when swapping $n, m$ and $s, t$, i.e. $\vec{k}$ and $-\vec{k}$, which is a consequence of statistical isotropy. The form of the expansion (5.6) is therefore a direct consequence of the symmetries of the problem, although the precise expression of $c_k(m, n, s, t)$ given in Eq. (5.7) required a non-trivial calculation.

Let us stress that only seven out of the ten Sp(4, $\mathbb{R}$) squeezing and rotation parameters enter Eq. (5.7). This is because, as explained in Appendix B, the rotation $\hat{\mathcal{R}}_{\vec{k}}(\varphi_k)$ only adds a global phase to the initial vacuum state, since it is invariant under rotations (in practice, this global phase is irrelevant, which reduces the number of effective parameters down to six). As a consequence, $q_{-}^k, q_{+}^k$ and $q_z^k$ are ineffective (had we started from a different initial state, those parameters would have entered the final result, see footnote 5). In practice, once the Hamiltonian of the system is specified, the dynamics can be integrated, which yields the seven relevant squeezing and rotation parameters, $r_1^k, r_2^k, \varphi_3^k, \theta_3^k, p_-^k, p_+^k$ and $p_z^k$, and thus fully determines the quantum state of the system at any time.

As a consistency check, one may verify that in the decoupled limit, the product of two two-mode squeezed states is recovered. Setting the two mixing parameters $\theta_5^k$ and $\theta_6^k$ to zero, Eq. (3.21) leads to $\theta^k = |\theta_4^k|$ and $\tau_k = 0$, so Eq. (3.20) gives $p_-^k = p_+^k = 0$ and $p_z^k = 2i\theta_4^k$. In Eq. (5.7), the overall factor $(p_-^k)^{i+j}$ is therefore non vanishing only when $i + j = 0$, i.e. $i = j = 0$ since $i$ and $j$ are non negative. Then, the overall factor $(p_+^k)^{i+j-s-t}$ is non vanishing only when $s + t = 0$, but since $s \le i = 0$ and $t \le j = 0$, this implies that $s = t = 0$, so only one term remains. This gives rise to

$$|\varphi_{\vec{k}}(t)\rangle = e^{-2i(\theta_3^k + \varphi_3^k)} \sum_{n=0}^{\infty} \frac{(-1)^n e^{-2in(\theta_3^k + \theta_4^k)}}{\cosh r_1^k} \tanh^n(r_1^k) |n_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}\rangle \\ \times \sum_{m=0}^{\infty} \frac{(-1)^m e^{-2in(\theta_5^k - \theta_4^k)}}{\cosh r_2^k} \tanh^m(r_2^k) |m_{\vec{k}}^{(2)}, m_{-\vec{k}}^{(2)}\rangle, \quad (5.9)$$

which defines the two-mode squeezed states coefficients $c_{i,k}(n)$. Up to a global irrelevant phase, this is indeed the product of two, uncoupled and disentangled, two-mode squeezed states [13].

### Expansion around the uncoupled limit

In order to gain some physical insight into the structure of the four-mode squeezed states, and having in mind possible comparisons with perturbative techniques that expand in the amplitude of the interaction Hamiltonian (i.e. in the terms in the Hamiltonian that mix the two sectors), let us now expand the evolved vacuum state (5.7) around the uncoupled limit (5.9). From Eq. (3.21), one can see that the cou-

pling parameters, i.e. $\theta_5^k, \varphi_5^k, \theta_6^k$ and $\varphi_6^k$, only appear through the combination $\tau_k = -(\theta_6^k - i\theta_5^k) = |\tau_k|e^{i\arg(\tau_k)}$, since $\theta^k = \sqrt{(\theta_4^k)^2 + |\tau_k|^2}$, and given that $\varphi_5^k$ and $\varphi_6^k$ are irrelevant, as explained above. As a consequence, an expansion around the uncoupled limit is an expansion in $|\tau_k|$. Upon expanding Eq. (5.7) up to quadratic order in $|\tau_k|$, one obtains

$$|\varphi_{\vec{k}}(t)\rangle = \sum_{n,m=0}^{\infty} c_{1,k}(n)c_{2,k}(m) \left\{ \begin{aligned} & |n_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}, m_{\vec{k}}^{(2)}, m_{-\vec{k}}^{(2)}\rangle \\ & + |\tau_k| \left[ \mathcal{F}_k(n, m+1) |1_{1\to 2}\rangle - \mathcal{F}_k^*(n+1, m) |1_{2\to 1}\rangle \right] \\ & + \frac{|\tau_k|^2}{2} \left[ \mathcal{F}_k(n, m+2)\mathcal{F}_k(n-1, m+1) |2_{1\to 2}\rangle \right. \\ & \qquad \left. + \mathcal{F}_k^*(n+2, m)\mathcal{F}_k^*(n+1, m-1) |2_{2\to 1}\rangle \right. \\ & \qquad \left. + 2\mathcal{F}_k^2(n, m+1) |1-1_{1\to 2}\rangle + 2(\mathcal{F}_k^*)^2(n+1, m) |1-1_{2\to 1}\rangle \right. \\ & \qquad \left. - 2\mathcal{F}_k(n, m)\mathcal{F}_k^*(n+1, m+1) |2_{1\leftrightarrow 2}\rangle \right. \\ & \qquad \left. + 2\mathcal{G}_k(n, m) |n_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}, m_{\vec{k}}^{(2)}, m_{-\vec{k}}^{(2)}\rangle \right] \end{aligned} \right\}, \quad (5.10)$$

where the Fock states that appear in this expansion are labeled by the number of particles being transferred from one sector to the other (while respecting statistical isotropy), and are given by

$$\begin{aligned} |1_{1\to 2}\rangle &\approx |(n-1)_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}, (m+1)_{\vec{k}}^{(2)}, m_{-\vec{k}}^{(2)}\rangle \\ &\quad + |n_{\vec{k}}^{(1)}, (n-1)_{-\vec{k}}^{(1)}, m_{\vec{k}}^{(2)}, (m+1)_{-\vec{k}}^{(2)}\rangle \\ |1_{2\to 1}\rangle &\approx |(n+1)_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}, (m-1)_{\vec{k}}^{(2)}, m_{-\vec{k}}^{(2)}\rangle \\ &\quad + |n_{\vec{k}}^{(1)}, (n+1)_{-\vec{k}}^{(1)}, m_{\vec{k}}^{(2)}, (m-1)_{-\vec{k}}^{(2)}\rangle \\ |2_{1\to 2}\rangle &\approx |(n-2)_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}, (m+2)_{\vec{k}}^{(2)}, m_{-\vec{k}}^{(2)}\rangle \\ &\quad + |n_{\vec{k}}^{(1)}, (n-2)_{-\vec{k}}^{(1)}, m_{\vec{k}}^{(2)}, (m+2)_{-\vec{k}}^{(2)}\rangle \\ |2_{2\to 1}\rangle &\approx |(n+2)_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}, (m-2)_{\vec{k}}^{(2)}, m_{-\vec{k}}^{(2)}\rangle \\ &\quad + |n_{\vec{k}}^{(1)}, (n+2)_{-\vec{k}}^{(1)}, m_{\vec{k}}^{(2)}, (m-2)_{-\vec{k}}^{(2)}\rangle \\ |2_{1\leftrightarrow 2}\rangle &\approx |(n+1)_{\vec{k}}^{(1)}, (n-1)_{-\vec{k}}^{(1)}, (m-1)_{\vec{k}}^{(2)}, (m+1)_{-\vec{k}}^{(2)}\rangle \\ &\quad + |(n-1)_{\vec{k}}^{(1)}, (n+1)_{-\vec{k}}^{(1)}, (m+1)_{\vec{k}}^{(2)}, (m-1)_{-\vec{k}}^{(2)}\rangle \\ |1-1_{1\to 2}\rangle &\approx |(n-1)_{\vec{k}}^{(1)}, (n-1)_{-\vec{k}}^{(1)}, (m+1)_{\vec{k}}^{(2)}, (m+1)_{-\vec{k}}^{(2)}\rangle \\ |1-1_{2\to 1}\rangle &\approx |(n+1)_{\vec{k}}^{(1)}, (n+1)_{-\vec{k}}^{(1)}, (m-1)_{\vec{k}}^{(2)}, (m-1)_{-\vec{k}}^{(2)}\rangle \end{aligned} \quad (5.11)$$

with the weighting functions

$$\mathcal{F}_k(n, m) = i \frac{\sin \theta_4^k}{\theta_4^k} e^{i[\theta_4^k + \arg(\tau_k)]} \sqrt{nm}, \quad (5.12)$$
$$G_k(n, m) = -2\mathcal{F}_k(n+1, m)\mathcal{F}_k^*(n+1, m) - i(n-m)\frac{\theta_4^k - e^{i\theta_4^k} \sin \theta_4^k}{(\theta_4^k)^2}. \quad (5.13)$$

These expressions can be interpreted as follows. At linear order in $\tau$, the only effect of the interaction is to add contributions from states where one particle has been exchanged between the two sectors [those are displayed in the second line of Eq. (5.10)]. The amplitude of these additional states is controlled by $\mathcal{F}_k$, which thus measures the rate at which particles transfer. It involves the phase $\theta_4^k$ through a function of order one, and the product of the numbers of particles in the two sectors, which mostly determines its amplitude. At quadratic order in $\tau$, the new states that appear in the expansion are obtained by exchanging two particles: either two particles with the same wavevector transfer from one sector to the other [third line of Eq. (5.10)], or two particles with opposite wavevector transfer from one sector to the other [fourth line of Eq. (5.10)], or two particles with opposite wavevector and from opposite sectors change sector [first term in the fifth line of Eq. (5.10)], or, finally, one particle changes sector and then moves back to its original sector [last term in the fifth line of Eq. (5.10)]. The amplitude of those states are now controlled by squared powers of $\mathcal{F}_k$, the only exception being the last state, controlled by $G_k$, which also involves the second term in Eq. (5.13). This can be interpreted as follows. When a particle transiently visits the opposite sector and is then reinstated, its journey to the other side modifies its state if the two sectors evolve differently. This is the case if $n \neq m$, which is why the second term in Eq. (5.13) is controlled by $n-m$.

One concludes that an expansion in the amplitude of the interaction, around the uncoupled limit, is essentially an expansion in the number of particles being exchanged: at order $|\tau_k|^p$, $p$ particles are transferred, which allows one to predict the form of the new states that appear in the expansion (e.g. $|p_1 \to 2\rangle = |(n-p)_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}, (m+p)_{\vec{k}}^{(2)}, m_{-\vec{k}}^{(2)}\rangle + |n_{\vec{k}}^{(1)}, (n-p)_{-\vec{k}}^{(1)}, m_{-\vec{k}}^{(2)}, (m+p)_{-\vec{k}}^{(2)}\rangle$, etc.). In particular, one can see that the diagonal elements, i.e. the amplitudes associated to the states of the form $|n_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}, m_{\vec{k}}^{(2)}, m_{-\vec{k}}^{(2)}\rangle$, where, within each sector, there are the same number of excitations in each wavevector, receive contributions from even powers of $|\tau_k|$ only. More generally, one can show that $c_k(n, m, s, t) = O(|\tau|^{s+|t|})$, and that $c_k(n, m, s, t)$ only receives contributions of order $|\tau|^{s+|t|+2q}$, where $q$ is a positive integer number counting the number of particles being sent away and then sent back to their original sector. This is described by the sums over $i$ and $j$ in Eq. (5.7), where $i$ counts the number of quanta travelling from sector 1 to sector 2 and travelling back, and $j$ counts the number of quanta travelling

along the opposite journey. This completes the physical interpretation of all six summation indices appearing in Eq. (5.7).

## 5.2 Phase-space representation

Although the quantum state derived in the Fock basis in Sect. 5.1 is enough to fully characterise the system, and in spite of the clear physical interpretation it yields, the seemingly complicated structure of Eq. (5.7) may call for alternative, simpler representations of the state. In this section, we provide such an alternative description by means of the Wigner function [40,41,46].

### Gaussian state

As we will see, the reason why the Wigner function is a convenient tool is that, the dynamics being linear, the evolved vacuum state is a Gaussian state, which is fully characterised by its two-point function

$$\Sigma_{\vec{k}, \vec{q}}(t) = \langle \emptyset | \hat{z}_{\vec{k}}(t) \hat{z}_{\vec{q}}^{\dagger}(t) | \emptyset \rangle, \quad (5.14)$$

which here is expressed in the Heisenberg picture. Using the statistical isotropy and homogeneity we have already invoked, different Fourier modes are uncoupled and one has

$$\Sigma_{\vec{k}, \vec{q}}(t) = \Sigma_k(t) \delta^3(\vec{k} - \vec{q}). \quad (5.15)$$

Our first goal is to relate $\Sigma_k$ to the Bogolyubov coefficients introduced in Sect. 3.3, since those describe the dynamical evolution of the system. In the helicity basis (2.13), Eq. (5.14) gives rise to

$$\Sigma_k(t) = D_k^{-1} U^\dagger \langle \emptyset | \hat{a}_{\vec{k}}(t) \hat{a}_{\vec{k}}^{\dagger}(t) | \emptyset \rangle U D_k^{-1}, \quad (5.16)$$

where the ladder operators evolve according to Eq. (2.22),<fnref n="7" /> i.e. by means of the Green matrix $G_k(t, t_{in})$. Since the Green matrix is an element of $\mathcal{S}p(4, \mathbb{R})$, it can be written down in terms of Bogolyubov coefficients as in Eqs. (3.24) and (3.26). Moreover, the ordering of the ladder operators in Eq. (2.13) is such that

$$\langle \emptyset | \hat{a}_{\vec{k}}(t_{in}) \hat{a}_{\vec{k}}^{\dagger}(t_{in}) | \emptyset \rangle = \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} \otimes I_2, \quad (5.17)$$

<footnote n="7">⁷ In the quantum representation constructed in Sect. 3.4, the evolution of the ladder operators is rather given by $\hat{a}_{\vec{k}}(t) = \hat{U}_{\vec{k}}^{\dagger}(t, t_{in})\hat{a}_{\vec{k}}(t_{in})\hat{U}_{\vec{k}}(t, t_{in})$, but one can show that this is strictly equivalent to Eq. (2.22). Indeed, upon writing the evolution operator $\hat{U}_{\vec{k}}$ in the fully factorised form derived in Sect. 4.2, computing the commutators between the ten elementary operations and the $\hat{a}_{\vec{k}}(t_{in})$ operators, and relating the Bogolyubov coefficients to the squeezing and rotation parameters using Eqs. (3.33)–(3.37), one can check that the same result (5.19)–(5.21) is obtained.</footnote>
so one obtains

$$ \Sigma_k(t) = \underbrace{\left( \begin{array}{cc} \mathrm{Cov}_k^{(\phi\phi)} & \mathrm{Cov}_k^{(\phi\pi)} \\ \mathrm{Cov}_k^{(\phi\pi)\mathrm{T}} & \mathrm{Cov}_k^{(\pi\pi)} \end{array} \right)}_{\mathrm{Cov}_k(t)} + \frac{i}{2} \Omega, \qquad (5.18) $$

where the second term, $i\Omega/2$, stems from the quantum commutator (2.11), and which defines the covariance matrix $\mathrm{Cov}_k(t)$. It is a real, symmetric, positive and semi-definite 4 × 4 matrix, where $\mathrm{Cov}_k^{(\phi\phi)}$, $\mathrm{Cov}_k^{(\pi\pi)}$ and $\mathrm{Cov}_k^{(\phi\pi)}$ are 2 × 2 matrices containing the position-position spectra, momentum-momentum spectra and position-momentum spectra, reading

$$ \mathrm{Cov}_{12,k}^{(\phi\phi)} = \frac{1}{2k} \mathrm{sinc}(\theta^k) |\tilde{\tau}_k| |\tau_k| \left\{ \begin{aligned} & \sin(\mathrm{arg}\tau_k + \mathrm{arg}\tilde{\tau}_k) \\ & \times [\cosh(2r_2^k) - \cosh(2r_1^k)] \\ & + \sin(2\theta_3^k - \mathrm{arg}\tau_k + \mathrm{arg}\tilde{\tau}_k) \sinh(2r_1^k) \\ & + \sin(2\theta_3^k + \mathrm{arg}\tau_k - \mathrm{arg}\tilde{\tau}_k) \sinh(2r_2^k) \end{aligned} \right\}. \qquad (5.25) $$

The momentum-momentum power spectra are given by

$$ \mathrm{Cov}_k^{(\phi\phi)} = \frac{1}{k} \left( \mathrm{Re} \left[ \beta_{11} (\beta_{21}^* + \alpha_{21}) + \beta_{12} (\beta_{22}^* + \alpha_{22}) \right] \mathrm{Re} \left[ \beta_{11} (\beta_{21}^* + \alpha_{21}) + \beta_{12} (\beta_{22}^* + \alpha_{22}) \right] \right) \quad (5.19) $$

$$ \mathrm{Cov}_k^{(\pi\pi)} = k \left( \mathrm{Re} \left[ \beta_{11} (\beta_{21}^* - \alpha_{21}) + \beta_{12} (\beta_{22}^* - \alpha_{22}) \right] \mathrm{Re} \left[ \beta_{11} (\beta_{21}^* - \alpha_{21}) + \beta_{12} (\beta_{22}^* - \alpha_{22}) \right] \right) \quad (5.20) $$

$$ \mathrm{Cov}_k^{(\phi\pi)} = \left( \mathrm{Im} \left[ \beta_{11} (\beta_{21}^* + \alpha_{21}) + \beta_{12} (\beta_{22}^* + \alpha_{22}) \right] - \mathrm{Im} \left[ \beta_{11} (\beta_{21}^* - \alpha_{21}) + \beta_{12} (\beta_{22}^* - \alpha_{22}) \right] \right) \cdot \mathrm{Im} \left( \alpha_{21} \beta_{21} + \alpha_{22} \beta_{22} \right) \quad (5.21) $$

An expression of $\mathrm{Cov}_k(t)$ in terms of the squeezing and rotation parameters can also be obtained from the Bloch–Messiah decomposition of $\mathcal{G}_k(t, t_{in})$ given in Eq. (3.32), where one notes that the squeezing and rotation parameters entering that decomposition are the same as those appearing in the evolution operator $\hat{U}_k^*(t, t_{in})$ given in Eq. (4.15), see footnote 7. In order to write the result in a compact form, we introduce the complex parameter

$$ \tilde{\tau}_k \equiv e^{-\frac{\mu z}{2}} = \cos(\theta^k) + i\theta_4^k \mathrm{sinc}(\theta^k), \qquad (5.22) $$

in terms of which the position-position power spectra read

$$ \mathrm{Cov}_{11,k}^{(\phi\phi)} = \frac{1}{2k} \times \left\{ |\tilde{\tau}_k|^2 \left[ \cosh(2r_1^k) + \cos(2\theta_3^k + 2\arg\tilde{\tau}_k) \sinh(2r_1^k) \right] + \sin^2(\theta^k) |\tau_k|^2 \left[ \cosh(2r_2^k) - \cos(2\theta_3^k + 2\arg\tau_k) \sinh(2r_2^k) \right] \right\}, \quad (5.23) $$

$$ \mathrm{Cov}_{22,k}^{(\phi\phi)} = \frac{1}{2k} \left\{ |\tilde{\tau}_k|^2 \left[ \cosh(2r_2^k) + \cos(2\theta_3^k - 2\arg\tilde{\tau}_k) \sinh(2r_2^k) \right] + \sin^2(\theta^k) |\tau_k|^2 \left[ \cosh(2r_1^k) - \cos(2\theta_3^k - 2\arg\tau_k) \sinh(2r_1^k) \right] \right\}, \quad (5.24) $$

$$ \mathrm{Cov}_{11,k}^{(\pi\pi)} = \frac{k}{2} \times \left\{ |\tilde{\tau}_k|^2 \left[ \cosh(2r_1^k) - \cos(2\theta_3^k + 2\arg\tilde{\tau}_k) \sinh(2r_1^k) \right] + \sin^2(\theta^k) |\tau_k|^2 \left[ \cosh(2r_2^k) + \cos(2\theta_3^k + 2\arg\tau_k) \sinh(2r_2^k) \right] \right\}, \quad (5.26) $$

$$ \mathrm{Cov}_{22,k}^{(\pi\pi)} = \frac{k}{2} \times \left\{ |\tilde{\tau}_k|^2 \left[ \cosh(2r_2^k) - \cos(2\theta_3^k - 2\arg\tilde{\tau}_k) \sinh(2r_2^k) \right] + \sin^2(\theta^k) |\tau_k|^2 \left[ \cosh(2r_1^k) + \cos(2\theta_3^k - 2\arg\tau_k) \sinh(2r_1^k) \right] \right\}, \quad (5.27) $$

$$ \mathrm{Cov}_{12,k}^{(\pi\pi)} = \frac{k}{2} \mathrm{sinc}(\theta^k) |\tilde{\tau}_k| |\tau_k| \times \left\{ \begin{aligned} & \sin(\mathrm{arg}\tau_k + \mathrm{arg}\tilde{\tau}_k) \\ & \times [\cosh(2r_2^k) - \cosh(2r_1^k)] \\ & - \sin(2\theta_3^k - \mathrm{arg}\tau_k + \mathrm{arg}\tilde{\tau}_k) \sinh(2r_1^k) \\ & - \sin(2\theta_3^k + \mathrm{arg}\tau_k - \mathrm{arg}\tilde{\tau}_k) \sinh(2r_2^k) \end{aligned} \right\}, \quad (5.28) $$
and the position-momentum power spectra by

$$ \text{Cov}_{11,k}^{(\phi\pi)} = \frac{1}{2} \left[ -|\tilde{\tau}_k|^2 \sin(2\theta_3^k + 2\arg\tilde{\tau}_k) \sinh(2r_1^k) + \text{sinc}^2(\theta^k) |\tau_k|^2 \sin(2\theta_3^k + 2\arg\tau_k) \sinh(2r_2^k) \right], \quad (5.29) $$

$$ \text{Cov}_{22,k}^{(\phi\pi)} = \frac{1}{2} \left[ -|\tilde{\tau}_k|^2 \sin(2\theta_3^k - 2\arg\tilde{\tau}_k) \sinh(2r_2^k) + \text{sinc}^2(\theta^k) |\tau_k|^2 \sin(2\theta_3^k - 2\arg\tau_k) \sinh(2r_1^k) \right], \quad (5.30) $$

$$ \text{Cov}_{12,k}^{(\phi\pi)} = \frac{1}{2} \text{sinc}(\theta^k) |\tilde{\tau}_k| |\tau_k| \left\{ \begin{aligned} & \cos(\arg\tau_k + \arg\tilde{\tau}_k) \\ & \times [\cosh(2r_1^k) - \cosh(2r_2^k)] \\ & + \cos(2\theta_3^k - \arg\tau_k + \arg\tilde{\tau}_k) \sinh(2r_1^k) \\ & + \cos(2\theta_3^k + \arg\tau_k - \arg\tilde{\tau}_k) \sinh(2r_2^k) \end{aligned} \right\}, \quad (5.31) $$

$$ \text{Cov}_{21,k}^{(\phi\pi)} = \frac{1}{2} \text{sinc}(\theta^k) |\tilde{\tau}_k| |\tau_k| \left\{ \begin{aligned} & \cos(\arg\tau_k + \arg\tilde{\tau}_k) \\ & \times [\cosh(2r_2^k) - \cosh(2r_1^k)] \\ & + \cos(2\theta_3^k - \arg\tau_k + \arg\tilde{\tau}_k) \sinh(2r_1^k) \\ & + \cos(2\theta_3^k + \arg\tau_k - \arg\tilde{\tau}_k) \sinh(2r_2^k) \end{aligned} \right\}, \quad (5.32) $$

where the remaining spectra are obtained from the symmetry of $\text{Cov}_k(t)$. We stress that, contrary to what was done around Eq. (5.10), there is no expansion in $\tau_k$ in these expressions. In practice, solving the dynamics of the system yields the Bogolyubov coefficients, or equivalently the squeezing and rotation parameters, from which the covariance matrix can be obtained from the above expressions. Finally, let us note that in the absence of interaction, $\tilde{\tau}_k = e^{i\theta_4^k}$ and $\tau_k = 0$, and Eqs. (5.23), (5.26) and (5.29) reduce to the power spectra describing a two-mode squeezed state (see e.g. Eqs. (6.18), (6.19) and (6.20) of Ref. [13]). This allows us to highlight the great similarity in the structure of the power spectra obtained from the two-mode and the four-mode squeezed states. For instance, in the expression of the diagonal elements of $\text{Cov}_k^{(\phi\phi)}$, $\text{Cov}_k^{(\pi\pi)}$ and $\text{Cov}_k^{(\phi\pi)}$, the first line is given by $|\tilde{\tau}_k|^2$ multiplied the power spectrum of a two-mode squeezed state in the sector of interest (position-position, momentum-momentum and position-momentum respectively). Since $|\tilde{\tau}_k|^2 \le 1$, this parameter can be interpreted as describing a loss of power in a given sector induced by its couplings to the other sector. On the contrary, the second line in these expressions is given by $|\tau_k|^2$ multiplied by the power spectrum of a two-mode squeezed state

in the other sector, and thus corresponds to an increase of power in a given sector coming from the opposite sector.

### Wigner function

Let us now introduce the Wigner function. By inverting Eq. (2.12), one can see that the field position and momentum operators, $\hat{\phi}_{i,\vec{k}}$ and $\hat{\pi}_{i,\vec{k}}$, involve creation and annihilation operators with opposite wavevectors. For convenience, it is useful to treat the two sectors $\vec{k}$ and $-\vec{k}$ separately, and to introduce the new position and momentum variables [12]

$$ \hat{\mathfrak{q}}_{i,\vec{k}} = \frac{1}{\sqrt{2k}} (\hat{\mathfrak{a}}_{i,\vec{k}} + \hat{\mathfrak{a}}_{i,\vec{k}}^{\dagger}), \quad (5.33) $$

$$ \hat{\mathfrak{p}}_{i,\vec{k}} = -i \sqrt{\frac{k}{2}} (\hat{\mathfrak{a}}_{i,\vec{k}} - \hat{\mathfrak{a}}_{i,\vec{k}}^{\dagger}), \quad (5.34) $$

with $i = 1, 2$, which do not mix opposite Fourier modes. When going from $\hat{\phi}_{i,\vec{k}}$ and $\hat{\pi}_{i,\vec{k}}$ to $\hat{\mathfrak{a}}_{i,\vec{k}}$ and $\hat{\mathfrak{p}}_{i,\vec{k}}$, one simply decomposes the four complex field variables, $z_{\vec{k}} = (\phi_{1,\vec{k}}, \phi_{2,\vec{k}}, \pi_{1,\vec{k}}, \pi_{2,\vec{k}})^{\mathrm{T}}$ into eight real field variables, $\mathfrak{q}_{\vec{k}} = (\mathfrak{q}_{1,\vec{k}}, \mathfrak{q}_{1,-\vec{k}}, \mathfrak{q}_{2,\vec{k}}, \mathfrak{q}_{2,-\vec{k}}, \mathfrak{p}_{1,\vec{k}}, \mathfrak{p}_{1,-\vec{k}}, \mathfrak{p}_{2,\vec{k}}, \mathfrak{p}_{2,-\vec{k}})^{\mathrm{T}}$, according to

$$ \phi_{i,\vec{k}} = \frac{1}{2} \left[ (\mathfrak{q}_{i,\vec{k}} + \mathfrak{q}_{i,-\vec{k}}) + \frac{i}{k} (\mathfrak{p}_{i,\vec{k}} - \mathfrak{p}_{i,-\vec{k}}) \right], \quad (5.35) $$

$$ \pi_{i,\vec{k}} = \frac{1}{2i} \left[ k (\mathfrak{q}_{i,\vec{k}} - \mathfrak{q}_{i,-\vec{k}}) + i (\mathfrak{p}_{i,\vec{k}} + \mathfrak{p}_{i,-\vec{k}}) \right]. \quad (5.36) $$

This transformation can be readily inverted, and is canonical.<fnref n="8" /> The helicity basis can be described by the eight-dimensional version of Eq. (2.13), i.e. $\hat{\mathfrak{a}}_{\vec{k}} \equiv (\hat{\mathfrak{a}}_{1,\vec{k}}, \hat{\mathfrak{a}}_{1,-\vec{k}}, \hat{\mathfrak{a}}_{2,\vec{k}}, \hat{\mathfrak{a}}_{2,-\vec{k}}, \hat{\mathfrak{a}}_{1,\vec{k}}^{\dagger}, \hat{\mathfrak{a}}_{1,-\vec{k}}^{\dagger}, \hat{\mathfrak{a}}_{2,\vec{k}}^{\dagger}, \hat{\mathfrak{a}}_{2,-\vec{k}}^{\dagger})^{\mathrm{T}}$, which is related to the vector $\hat{\mathfrak{q}}_{\vec{k}}$ via

$$ \hat{\mathfrak{a}}_{\vec{k}} = [(\mathfrak{U} \mathfrak{D}_k) \otimes \mathfrak{I}_2] \hat{\mathfrak{q}}_{\vec{k}}, \quad (5.37) $$

and whose commutators are given by

$$ [\hat{\mathfrak{a}}_{\vec{k}}, \hat{\mathfrak{a}}_{\vec{k}}^{\dagger}] = i (\mathcal{J} \otimes \mathfrak{I}_2) \delta^3(\vec{k} - \vec{q}). \quad (5.38) $$

The real variables $\mathfrak{q}_{\vec{k}}$ can be used to introduce the Wigner-Weyl transform [47–49], which allows one to connect quantum operators $\hat{\mathcal{O}}_{\vec{k}}$ to classical functions in the phase space $\hat{\mathcal{O}}_{\vec{k}}(\mathfrak{q}_{\vec{k}})$, according to

$$ \hat{\mathcal{O}}_{\vec{k}}(\mathfrak{q}_{\vec{k}}) = \frac{1}{(2\pi)^4} \int_{\mathbb{R}^4} dx_1 dx_2 dy_1 dy_2 e^{-i\mathfrak{p}_{1,\vec{k}}x_1 - i\mathfrak{p}_{2,\vec{k}}x_2 - i\mathfrak{p}_{1,-\vec{k}}y_1 - i\mathfrak{p}_{2,-\vec{k}}y_2} $$

<footnote n="8">⁸ This can be seen by introducing $\hat{\phi}_{i,\vec{k}}^{\mathrm{R}} = (\hat{\phi}_{i,\vec{k}} + \hat{\phi}_{i,\vec{k}}^{\dagger})/\sqrt{2}$, $\hat{\phi}_{i,\vec{k}}^{\mathrm{I}} = (\hat{\phi}_{i,\vec{k}} - \hat{\phi}_{i,\vec{k}}^{\dagger})/(\sqrt{2}i)$ and similarly for $\hat{\pi}_{i,\vec{k}}^{\mathrm{R}}$ and $\hat{\pi}_{i,\vec{k}}^{\mathrm{I}}$, and by showing that the $8 \times 8$ matrix that relates those variables to $\mathfrak{q}_{\vec{k}}$ satisfies the symplectic relation (2.19) in 8 dimensions.</footnote>
$$ \begin{aligned} & \langle \mathbf{q}_{1,\bar{k}} + \frac{x_1}{2}, \mathbf{q}_{2,\bar{k}} + \frac{x_2}{2}, \mathbf{q}_{1,-\bar{k}} + \frac{y_1}{2}, \mathbf{q}_{2,-\bar{k}} + \frac{y_2}{2} | \hat{O}_{\bar{k}} \\ & \quad | \mathbf{q}_{1,\bar{k}} - \frac{x_1}{2}, \mathbf{q}_{2,\bar{k}} - \frac{x_2}{2}, \mathbf{q}_{1,-\bar{k}} - \frac{y_1}{2}, \mathbf{q}_{2,-\bar{k}} - \frac{y_2}{2} \rangle. \end{aligned} \quad (5.39) $$

In general, the quantum state of a system can be equivalently described in terms of its density matrix $\hat{\rho}(t) = |\phi(t)\rangle\langle\phi(t)| = \prod_{\bar{k}\in\mathbb{R}^{3+}} \hat{\rho}_{\bar{k}}(t)$, or in terms of its Wigner function $W(t) = \prod_{\bar{k}\in\mathbb{R}^{3+}} W_{\bar{k}}(\mathbf{q}_{\bar{k}}, t)$, which is the Wigner-Weyl transform of the density matrix [50,51].

The inverse Wigner-Weyl transform, which allows one to go from the Wigner function to the density matrix, is given in Appendix D, see Eq. (D.1). The Wigner functions $W_{\bar{k}}(\mathbf{q}_{\bar{k}}, t)$ can be interpreted as quasi distribution functions in the sense that expectation values of quantum operators can be expressed as

$$ \begin{aligned} \langle\phi_{\bar{k}}(t)|\hat{O}_{\bar{k}}|\phi_{\bar{k}}(t)\rangle &= \operatorname{Tr}[\hat{\rho}_{\bar{k}}(t)\hat{O}_{\bar{k}}] \\ &= (2\pi)^4 \int d\mathbf{q}_{\bar{k}} \tilde{O}_{\bar{k}}(\mathbf{q}_{\bar{k}}) W_{\bar{k}}(\mathbf{q}_{\bar{k}}, t), \end{aligned} \quad (5.40) $$

where the $(2\pi)^4$ prefactor can be absorbed in the normalisation of the Wigner function if needed. For Gaussian states, as shown e.g. in Appendix G of Ref. [12], the Wigner functions are Gaussian functions (hence the statement that the state is “Gaussian”) and read

$$ W_{\bar{k}}(\mathbf{q}_{\bar{k}}, t) = \frac{1}{(2\pi)^4 \sqrt{\det \mathrm{Cov}_{8\times 8, k}}} e^{-\frac{1}{2} \mathbf{q}_{\bar{k}}^{\mathrm{T}} \mathrm{Cov}_{8\times 8, k}^{-1} \mathbf{q}_{\bar{k}}}, \quad (5.41) $$

where $\mathrm{Cov}_{8\times 8, k}$ is the $8 \times 8$ covariance matrix in the field basis $\mathbf{q}_{\bar{k}}$, defined through

$$ \langle\phi|\hat{\mathbf{q}}_{\bar{k}}(t)\hat{\mathbf{q}}_{\bar{k}}^{\dagger}(t)|\phi\rangle = \mathrm{Cov}_{8\times 8, k}(t) + \frac{i}{2}\mathbf{\Omega} \otimes \mathbf{I}_2. \quad (5.42) $$

It can be computed in a similar way as what was done around Eq. (5.16), and Eq. (5.37) gives rise to

$$ \mathrm{Cov}_{8\times 8, k}(t) = \mathrm{Cov}_k(t) \otimes \mathbf{I}_2, \quad (5.43) $$

where the inverse is simply given by $\mathrm{Cov}_{8\times 8, k}^{-1} = \mathrm{Cov}_k^{-1} \otimes \mathbf{I}_2$, see footnote 3.

The advantage of working with the Wigner functions $W_{\bar{k}}(\mathbf{q}_{\bar{k}}, t)$ rather than the quantum state or the density matrix in the Fock basis, lies in the simplicity of the Gaussian function (5.41). It describes entirely the quantum state of the system from the knowledge of its two-point correlation functions, and can therefore greatly simplify some calculations, as will be made explicit in the next section. The drawback of this approach is that the entangled structure, otherwise easily interpretable from Eq. (5.7), is now hidden in the details of the power spectra of the system.

# 6 Decoherence

Having determined in Sect. 5 the quantum state of a linear two-field system, both in the Fock basis where one obtains direct products of four-mode squeezed states and at the level of the Wigner function, we now consider the case where one of the two fields is unobservable and can be traced over. This corresponds to situations where measurements can only be performed on the first field, dubbed the “system”, while the second field, dubbed the “environment”, cannot be directly accessed. When the two fields become entangled, this leads to the concept of quantum decoherence [22], and we choose to illustrate the usefulness of the various tools introduced above by studying this notion for the systems at hand.

## 6.1 Reduced density matrix

Let $\hat{O} = \hat{O}_1 \otimes \hat{\mathbb{I}}_2$ be a quantum operator describing an observable of the first field only, where $\hat{\mathbb{I}}_2$ denotes the identity acting in the second-field sector (i.e. the environment). Its expectation value is given by $\operatorname{Tr}(\hat{\rho}\hat{O}) = \operatorname{Tr}_1\operatorname{Tr}_2[\hat{\rho}(\hat{O}_1 \otimes \hat{\mathbb{I}}_2)] = \operatorname{Tr}_1[\operatorname{Tr}_2(\hat{\rho})\hat{O}_1]$, see Eq. (5.40), where $\operatorname{Tr}_i$ denotes the trace over the degrees of freedom contained in the field $i$. The expectation value of $\hat{O}$ can therefore be obtained from the reduced density matrices

$$ \begin{aligned} \hat{\rho}_{\bar{k}, \mathrm{red}} &= \operatorname{Tr}_2(\hat{\rho}_{\bar{k}}) \\ &= \sum_{n,m=0}^{\infty} [\mathbf{I}_2^{(1)} \otimes \langle n_{\bar{k}}^{(2)}, m_{-\bar{k}}^{(2)} |] \hat{\rho}_{\bar{k}} \\ & \quad \times [\mathbf{I}_2^{(1)} \otimes |n_{\bar{k}}^{(2)}, m_{-\bar{k}}^{(2)}\rangle] \end{aligned} \quad (6.1) $$

by tracing over in the first-field sector (i.e. the system). In Eq. (6.1), for explicitness, the trace over the environmental degrees of freedom has been expanded in the Fock basis. The reduced density matrix thus contains all accessible information about the system.

If the two fields are entangled, the reduced density matrix follows non-unitary evolution, and describes a mixed (as opposed to pure) state. This can be described by the so-called purity

$$ \gamma_{\bar{k}} = \operatorname{Tr}(\hat{\rho}_{\bar{k}, \mathrm{red}}^2). \quad (6.2) $$

For a pure state, $\rho^2 = \rho$ and $\gamma = 1$, while mixed states have $1/d \le \gamma \le 1$ in general, where $d$ is the dimension of the Fock space of the system (which is infinite in the present case) and the limit $\gamma \to 1/d$ corresponds to a maximally decohered state. Therefore, purity provides a measure of the information loss into the environment, hence of decoherence. Note that it is simply related to the linear entropy $S_{\bar{k}, \mathrm{lin}} = 1 - \gamma_{\bar{k}}$, which itself provides a lower bound to the entanglement entropy $S_{\bar{k}, \mathrm{ent}} = -\operatorname{Tr}(\rho \ln \rho)$, and which, as the entangle-
ment entropy, characterises one's ignorance about the state of a system.

For the system at hand, one can check that if the two fields are uncoupled, then the reduced density matrix is pure. Indeed, if the quantum state is given by Eq. (5.9), the reduced density matrix reads

$$ \begin{aligned} \hat{\rho}_{\vec{k},\text{red}} &= \text{Tr}_2 \left( \sum_{n,m,\vec{n},\vec{m}=0}^{\infty} c_{1,k}(n)c_{1,k}^*(\vec{n})c_{2,k}(m)c_{2,k}^*(\vec{m}) \right. \\ & \qquad \left. \left| n_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}m_{\vec{k}}^{(2)}, m_{-\vec{k}}^{(2)} \right\rangle \left\langle \vec{n}_{\vec{k}}^{(1)}, \vec{n}_{-\vec{k}}^{(1)}\vec{m}_{\vec{k}}^{(2)}, \vec{m}_{-\vec{k}}^{(2)} \right| \right) \\ &= \sum_{m=0}^{\infty} |c_{2,k}(m)|^2 \sum_{n,\vec{n}=0}^{\infty} c_{1,k}(n)c_{1,k}^*(\vec{n}) \left| n_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)} \right\rangle \\ & \qquad \left. \times \left| \vec{n}_{\vec{k}}^{(1)}, \vec{n}_{-\vec{k}}^{(1)} \right| \right\rangle, \end{aligned} \quad (6.3) $$

where we have performed the trace in the Fock basis as in Eq. (6.1). From the expression of the $c_{i,k}(n)$ coefficients given in Eq. (5.9), one can easily check that $\sum_n |c_{i,k}(n)|^2 = 1$, hence the reduced density matrix is the one of a (pure) two-mode squeezed state, and using again the identity $\sum_n |c_{i,k}(n)|^2 = 1$, it has purity $\gamma_{\vec{k}} = 1$.

In general, for the evolved vacuum state given in Eq. (5.6), the density matrix $\hat{\rho}_{\vec{k}}(t) = |\phi_{\vec{k}}(t)\rangle\langle\phi_{\vec{k}}(t)|$ is written down explicitly in Appendix 1, where it is shown that the tracing-out procedure of Eq. (6.1) gives rise to

$$ \hat{\rho}_{\vec{k},\text{red}} = \sum_{n=0}^{\infty} \sum_{n'==-\infty}^{\infty} \sum_{s,t=-\min(n,n')}^{\infty} \Xi_k(n, n', s, t) \\ \qquad |(n+s)_{\vec{k}}^{(1)}, (n+t)_{-\vec{k}}^{(1)}\rangle\langle(n'+s)_{\vec{k}}^{(1)}, (n'+t)_{-\vec{k}}^{(1)}| \quad (6.4) $$

with

$$ \begin{aligned} \Xi_k(n, n', s, t) &= \sum_{m=\max(0, s, t)}^{\infty} \sum_{m'=-m}^{n'} c_k(n, m, s, t) c_k^* \\ &= (n' - m', m + m', s + m', t + m'). \end{aligned} \quad (6.5) $$

This expression differs from the uncoupled case (6.3) due to the presence of the s and t indices, related to the transfer of excitations from one sector to the other. One can check that the invariance of the $c_k$ coefficients under exchanging s and t, which was noted below Eq. (5.7) as a manifestation of the statistical isotropy of the state, guarantees that the $\Xi_k$ coefficients are also invariant under swapping s and t, hence the reduced state is also statistically isotropic.

Let us also recall that if one performs a perturbative expansion in the interaction Hamiltonian, that is, as argued

in Sect. 5.1, an expansion in $\tau_k$, then the corrections to $c_k(n, m, s, t)$ are of order $O(|\tau_k|^{s+|t|+2q})$, where $q$ is a non-negative integer number that stands for the number of particles that change sectors and then change back. As a consequence, the corrections to $\Xi_k(n, n', s, t)$ are of order $O(|\tau_k|^{2(|s|+|t|+Q)})$, where $Q = m' + q + q'$, where $q$ and $q'$ are associated with the two $c_k$ coefficients appearing in Eq. (6.5). This implies that, while we saw that the diagonal elements of the full quantum state received even corrections in $|\tau_k|$ only, see the discussion at the very end of Sect. 5.1, we have now showed that all the entries of the reduced density matrix receive even corrections in $|\tau_k|$ only. As a consequence, the leading correction to observables performed on the system is always quadratic in the coupling constants, in agreement with the results of Ref. [23]. From Eq. (6.4), the purity (6.2) can also be obtained, and in Appendix C it is shown that

$$ \begin{aligned} \gamma_{\vec{k}}(t) = & \sum_{n=0}^{\infty} \sum_{n', u=-\infty}^{\infty} \sum_{s,t=-\min(n,n')}^{\infty} \Xi_k(n, n', s, t) \\ & \Xi_k(n' - u, n - u, s + u, t + u). \end{aligned} \quad (6.6) $$

Let us stress that the above formulas are not perturbative and allow one to compute the reduced density matrix and its purity up to arbitrary order in the interaction terms.

### Small-coupling limit

In order to gain more insight, one can however derive the leading-order result in the interaction parameter $|\tau_k|$. After a lengthy though straightforward calculation, the reduced density matrix derived from the perturbed evolved vacuum state (5.10) is given by

$$ \begin{aligned} \hat{\rho}_{\vec{k},\text{red}}(t) = & \sum_{n,n'=0}^{\infty} c_{1,k}(n)c_{1,k}^*(n') |n_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}\rangle\langle n'_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}| \\ & + |\tau_k|^2 \sum_{n,n',m=0}^{\infty} c_{1,k}(n)c_{1,k}^*(n') \\ & \times \left\{ |c_{2,k}(m)|^2 \mathcal{F}_k(n, m+1) \mathcal{F}_k^*(n', m+1) \right. \\ & \qquad \left. \times \left[ |(n-1)_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}\rangle\langle(n'-1)_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}| \right. \right. \\ & \qquad \left. \left. + |n_{\vec{k}}^{(1)}, (n-1)_{-\vec{k}}^{(1)}\rangle\langle n'_{\vec{k}}^{(1)}, (n'-1)_{-\vec{k}}^{(1)}| \right] \right. \\ & \qquad \left. + |c_{2,k}(m)|^2 \mathcal{F}_k^*(n+1, m) \mathcal{F}_k(n'+1, m) \right. \\ & \qquad \left. \times \left[ |(n+1)_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}\rangle\langle(n'+1)_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}| \right. \right. \\ & \qquad \left. \left. + |n_{\vec{k}}^{(1)}, (n+1)_{-\vec{k}}^{(1)}\rangle\langle n'_{\vec{k}}^{(1)}, (n'+1)_{-\vec{k}}^{(1)}| \right] \right. \\ & \qquad \left. - c_{2,k}(m)c_{2,k}^*(m+1) \right. \\ & \qquad \left. \times \mathcal{F}_k(n, m+1) \mathcal{F}_k(n'+1, m+1) \right\} \end{aligned} $$
$$
\begin{align}
& \times \left[ \begin{aligned}[t]
& |(n-1)_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}\rangle \langle n'_{\vec{k}}^{(1)}, (n'+1)_{-\vec{k}}^{(1)}| \\
& + |n'_{\vec{k}}^{(1)}, (n-1)_{-\vec{k}}^{(1)}\rangle \langle n'_{-\vec{k}}^{(1)}, (n'+1)_{-\vec{k}}^{(1)}|
\end{aligned} \right] \nonumber \\
& - c_{2,k}(m) c_{2,k}^*(n-1, m) \mathcal{F}_k^*(n', m) \nonumber \\
& \times \left[ \begin{aligned}[t]
& |(n+1)_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)}\rangle \langle n'_{\vec{k}}^{(1)}, (n'-1)_{-\vec{k}}^{(1)}| \\
& + |n'_{\vec{k}}^{(1)}, (n+1)_{-\vec{k}}^{(1)}\rangle \langle n'_{-\vec{k}}^{(1)}, (n'-1)_{-\vec{k}}^{(1)}|
\end{aligned} \right] \nonumber \\
& + c_{2,k}(m) c_{2,k}^*(m-1) \left[ \mathcal{F}_k^*(n', m) \right]^2 \nonumber \\
& \times \left| n'_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)} \right\rangle \langle n'_{\vec{k}}^{(1)}, (n'-1)_{-\vec{k}}^{(1)} \right| \nonumber \\
& + c_{2,k}(m) c_{2,k}^*(m+1) \left[ \mathcal{F}_k(n, m+1) \right]^2 \nonumber \\
& \times \left| (n-1)_{\vec{k}}^{(1)}, (n-1)_{-\vec{k}}^{(1)} \right\rangle \langle n'_{\vec{k}}^{(1)}, n'_{-\vec{k}}^{(1)} \right| \nonumber \\
& + c_{2,k}(m) c_{2,k}^*(m+1) \left[ \mathcal{F}_k(n'+1, m+1) \right]^2 \nonumber \\
& \times \left| n'_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)} \right\rangle \langle n'_{\vec{k}}^{(1)}, (n'+1)_{-\vec{k}}^{(1)} \right| \nonumber \\
& + c_{2,k}(m) c_{2,k}^*(m-1) \left[ \mathcal{F}_k^*(n+1, m) \right]^2 \nonumber \\
& \times \left| (n+1)_{\vec{k}}^{(1)}, (n+1)_{-\vec{k}}^{(1)} \right\rangle \langle n'_{\vec{k}}^{(1)}, n'_{-\vec{k}}^{(1)} \right| \nonumber \\
& + |c_{2,k}(m)|^2 \left[ \mathcal{G}_k^*(n', m) + \mathcal{G}_k(n, m) \right] \nonumber \\
& \times \left| n'_{\vec{k}}^{(1)}, n_{-\vec{k}}^{(1)} \right\rangle \langle n'_{\vec{k}}^{(1)}, n'_{-\vec{k}}^{(1)} \right| . \tag{6.7}
\end{align}
$$

One can see that the leading correction in $|\tau_k|$ is indeed of quadratic order [despite the full state (5.10) having linear contributions], in agreement with the above discussion. This expression allows one to evaluate the purity (6.2) and one obtains

$$
\begin{equation}
\begin{split}
\gamma_{\vec{k}}(t) = 1 + 4 |\tau_k|^2 \sum_{n,m=0}^{\infty} \Re \left\{ |c_{1,k}(n)|^2 |c_{2,k}(m)|^2 \mathcal{G}_k(n,m) \right. \\
\left. \qquad + 2c_{1,k}(n+1)c_{1,k}^*(n)c_{2,k}(m)c_{2,k}^*(m+1) \right. \\
\left. \qquad \left[ \mathcal{F}_k(n+1,m+1) \right]^2 \right\},
\end{split}
\tag{6.8}
\end{equation}
$$

where the relation $\sum_n |c_{i,k}(n)|^2 = 1$ has been used. Making use of Eqs. (5.12) and (5.13), the sums appearing in the expression can be carried out explicitly, and in terms of the squeezing and rotation parameters, one obtains

$$
\gamma_{\vec{k}} = 1 - 4 |\tau_k|^2 \sin^2(\theta_4^k) \left[ \sinh^2(r_1^k - r_2^k) + \cos^2(\theta_4^k - \arg \tau_k) \sinh(2r_1^k) \sinh(2r_2^k) \right]. \quad (6.9)
$$

One can check that, as $|\tau_k|$ increases away from 0, $\gamma_{\vec{k}}$ decreases away from 1, since the term inside the squared braces in Eq. (6.9) is always positive. One can also see that larger squeezing amplitudes $r_1^k$ and $r_2^k$ lead to more efficient decoherence, as usually encountered [11]. It is finally worth pointing out that in the specific configuration where $r_1^k = r_2^k$

and $\theta_4^k = \arg \tau_k \pm \pi$, the leading-order correction to the purity vanishes.

## 6.2 Reduced Wigner function

As argued in Sect. 5.2, a complementary (and sometimes simpler from a computational standpoint) tool to analyse the four-mode squeezed state is provided by the Wigner function. Let us first establish how the Wigner function of the reduced system can be obtained from the one of the full two-field setup.

We consider again an operator of the form $\hat{O}_{\vec{k}} = \hat{O}_{1,\vec{k}} \otimes \hat{1}_2$. From Eq. (5.39), its Wigner-Weyl transform is simply given by $\tilde{O}_{\vec{k}} = \tilde{O}_{1,\vec{k}}/(2\pi)^2$, where $\tilde{O}_{1,\vec{k}}$ is the Wigner-Weyl transform of $\hat{O}_{1,\vec{k}}$ within the first-field sector,

$$
\begin{equation}
\begin{split}
\tilde{O}_{1,\vec{k}}(\mathbf{q}_{1,\vec{k}}) = \int_{\mathbb{R}^2} \frac{dx_1 dy_1}{(2\pi)^2} e^{-i\mathbf{p}_{1,\vec{k}}x_1 - i\mathbf{p}_{1,-\vec{k}}y_1} \\
& \qquad \left\langle \mathbf{q}_{1,\vec{k}} + \frac{x_1}{2}, \mathbf{q}_{1,-\vec{k}} + \frac{y_1}{2} \right| \hat{O}_{1,\vec{k}} \left| \mathbf{q}_{1,\vec{k}} - \frac{x_1}{2}, \mathbf{q}_{1,-\vec{k}} - \frac{y_1}{2} \right\rangle,
\end{split}
\tag{6.10}
\end{equation}
$$

with $\mathbf{q}_{1,\vec{k}} \equiv (\mathbf{q}_{1,\vec{k}}, \mathbf{q}_{1,-\vec{k}}, \mathbf{p}_{1,\vec{k}}, \mathbf{p}_{1,-\vec{k}})^{\mathrm{T}}$. Plugging this result into Eq. (5.40), the expectation value of $\hat{O}_{\vec{k}}$ is given by

$$
\begin{equation}
\begin{split}
\langle \emptyset_{\vec{k}}(t) | \hat{O}_{\vec{k}} | \emptyset_{\vec{k}}(t) \rangle ={}& (2\pi)^2 \int d\mathbf{q}_{1,\vec{k}} \, \tilde{O}_{1,\vec{k}}(\mathbf{q}_{1,\vec{k}}) \\
& \qquad W_{\vec{k},\mathrm{red}}(\mathbf{q}_{1,\vec{k}}, t),
\end{split}
\tag{6.11}
\end{equation}
$$

where we have defined the reduced Wigner function

$$
W_{\vec{k},\mathrm{red}}(\mathbf{q}_{1,\vec{k}}, t) = \int_{\mathbb{R}^4} d^4 \mathbf{q}_{2,\vec{k}} \, W_{\vec{k}}(\mathbf{q}_{1,\vec{k}}, \mathbf{q}_{2,\vec{k}}, t). \quad (6.12)
$$

Comparing Eq. (6.11) with Eq. (5.40), one can see that $W_{\vec{k},\mathrm{red}}$ can be used to compute quantum expectation values of observables in the first-field space. Therefore, it corresponds to the Wigner function in the reduced phase space (where the different powers of $2\pi$ come from the different dimensions of the phase spaces). In other words, $W_{\vec{k},\mathrm{red}}$ is the Wigner-Weyl transform of $\hat{\rho}_{\vec{k},\mathrm{red}}(t)$, which is shown explicitly in Appendix D. This is why partial trace in the Hilbert space is equivalent to partial integration in the phase space.

If the Wigner function is Gaussian, this partial integration can be easily done, and marginalisation over a phase-space variable simply corresponds to removing the associated lines and columns in the covariance matrix. From Eq. (5.41), this implies that

$$
W_{\vec{k},\mathrm{red}}(\mathbf{q}_{1,\vec{k}}, t) = \frac{1}{(2\pi)^2} \sqrt{\det \mathrm{Cov}_{\vec{k},\mathrm{red}}} e^{-\frac{1}{2} \mathbf{q}_{1,\vec{k}}^{\mathrm{T}} \mathrm{Cov}_{\vec{k},\mathrm{red}}^{-1} \mathbf{q}_{1,\vec{k}}}, \quad (6.13)
$$
where $\text{Cov}_{\tilde{k},\text{red}}$ is obtained from Eq. (5.43) by removing the lines and columns related to the second sector (i.e. the third, the fourth, the seventh and the eighth lines and columns).

The purity of the reduced system can then be computed from the reduced Wigner function, by using the property of the Wigner-Weyl transform [46]

$$ \mathrm{Tr} (\hat{A}\hat{B}) = (2\pi)^2 \int_{\mathbb{R}^4} d^4 \mathbf{q}_{1,\tilde{k}} \tilde{A} \tilde{B}, \qquad (6.14) $$

where $\hat{A}$ and $\hat{B}$ are two quantum operators acting on the first-field system and $\tilde{A}$ and $\tilde{B}$ are their Wigner-Weyl transforms. With $\hat{A} = \hat{B} = \hat{\rho}_{\tilde{k},\text{red}}$, Eq. (6.14) gives rise to the following expression for the purity (6.2),

$$ \gamma_{\tilde{k}}(t) = (2\pi)^2 \int_{\mathbb{R}^4} d^4 \mathbf{q}_{1,\tilde{k}} W_{\tilde{k},\text{red}}^2(\mathbf{q}_{1,\tilde{k}}, t). \qquad (6.15) $$

Plugging Eq. (6.13) into that formula, one obtains

$$ \begin{aligned} \gamma_{\tilde{k}}(t) &= (16 \det \mathrm{Cov}_{\tilde{k},\text{red}})^{-1/2} \\ &= \frac{1}{4} \left[ \mathrm{Cov}_{11,k}^{(\phi\phi)} \mathrm{Cov}_{11,k}^{(\pi\pi)} - \left( \mathrm{Cov}_{11,k}^{(\phi\pi)} \right)^2 \right]^{-1}. \end{aligned} \qquad (6.16) $$

As a consequence, for a Gaussian state, the purity of the system can be directly evaluated from the knowledge of the power spectra in the observable sector [52]. More precisely, the power spectra appear through a specific combination, i.e. the determinant of the (reduced) covariance matrix, which makes the purity invariant under canonical transformations [13]. Let us indeed recall that in a two-dimensional (i.e. single-field) system, there is a single symplectic invariant [52,53], the so-called symplectic eigenvalue $\sigma_{\tilde{k}}(t)$, which is such that the eigenvalues of $\text{Cov}_{\tilde{k},\text{red}}\Omega$ are given by $\pm i\sigma_{\tilde{k}}(t)$. This leads to $\sigma_{\tilde{k}} = 1/(2\sqrt{\gamma_{\tilde{k}}})$, and the condition $\gamma_{\tilde{k}} < 1$ is equivalent to $\sigma_{\tilde{k}} > 1/2$. In passing, let us note that, for Gaussian states, the entanglement entropy introduced below Eq. (6.2) is related to the symplectic eigenvalue by [54] $S_{\tilde{k},\text{ent}} = (\sigma_{\tilde{k}} + 1/2) \log_2(\sigma_{\tilde{k}} + 1/2) - (\sigma_{\tilde{k}} - 1/2) \log_2(\sigma_{\tilde{k}} - 1/2)$. This allows one to relate the linear and entanglement entropies, and check the above statement that the former provides a lower bound to the latter.

If the two fields are uncoupled, the reduced system is in a pure state and one has [13] $\text{Cov}_{11,k}^{(\phi\phi)} \text{Cov}_{11,k}^{(\pi\pi)} - (\text{Cov}_{11,k}^{(\phi\pi)})^2 = 1/4$. This implies that $\gamma_{\tilde{k}}(t) = 1$, $\sigma_{\tilde{k}}(t) = 1/2$ and $S_{\tilde{k},\text{lin}}(t) = S_{\tilde{k},\text{ent}}(t) = 0$. Otherwise, $\gamma_{\tilde{k}} < 1$ signals the presence of decoherence. An obvious, yet crucial consequence of Eq. (6.16) is that decoherence, i.e. the reduction of $\gamma_{\tilde{k}}$ away from unity, cannot be achieved without modifying the power spectra. In other words, for any system that undergoes decoherence, the observational predictions are necessarily altered, and an important question that will be addressed below is whether or not decoherence can proceed while keeping this alteration negligible [23,55].

Let us finally stress that unlike the approach presented in Sect. 6.1, which relies on a detailed analysis of the mathematical structure of Sp(4, $\mathbb{R}$) and leads to a formula, Eq. (6.6), that involves nine infinite sums; the Wigner function formalism only makes use of Gaussian integrals. It can therefore be straightforwardly generalised to higher-dimension systems (i.e. containing more fields), while the approach of Sect. 6.1 would require further analyses of the groups Sp(2n, $\mathbb{R}$). By plugging Eqs. (5.23), (5.26) and (5.29) into Eq. (6.16), one can finally obtain a fully non-perturbative expression for the purity in terms of the squeezing parameters

$$ \gamma_{\tilde{k}}(t) = \left\{ |\tilde{\tau}_k|^4 + \mathrm{sinc}^4(\theta^k) |\tau_k|^4 + 2\mathrm{sinc}^2(\theta^k) |\tau_k|^2 |\tilde{\tau}_k|^2 \right. \\ \left. \times \left[ \cosh(2r_1^k) \cosh(2r_2^k) + \cos(2\arg\tilde{\tau}_k - 2\arg\tau_k) \right] \times \sinh(2r_1^k) \sinh(2r_2^k) \right\}^{-1}, \qquad (6.17) $$

where we recall that $\tilde{\tau}_k$ is defined in Eq. (5.22), and $\tau_k$ and $\theta^k$ in Eq. (3.21).

### Small-coupling limit

Similarly to what was done at the end of Sect. 6.1, let us now further expand the purity in the limit where the system and the environment fields are weakly coupled. The power spectra in the observable sector are given by Eqs. (5.23), (5.26) and (5.29). Expanding these expressions up to second order in the interaction parameter $|\tau_k|$ [recalling that $(\theta_6^k - i\theta_5^k) = |\tau_k|e^{i\arg\tau_k}$ and $\theta^k = \sqrt{(\theta_4^k)^2 + |\tau_k|^2}$], one obtains

$$ \begin{aligned} \mathrm{Cov}_{11,k}^{(\phi\phi)} = \frac{1}{2k} \left( \left[ \cosh(2r_1^k) + \cos(2\theta_3^k + 2\theta_4^k) \sinh(2r_1^k) \right] \right. \\ \left. + |\tau_k|^2 \left\{ - \frac{\sin^2 \theta_4^k}{\theta_4^k} \cosh(2r_1^k) \right. \right. \\ \left. \left. + \frac{1}{2(\theta_4^k)^2} \left[ \cos(2\theta_3^k) - \cos(2\theta_3^k + 2\theta_4^k) \right. \right. \\ \left. \left. - 2\theta_4^k \sin(2\theta_3^k + 2\theta_4^k) \right] \sinh(2r_1^k) \right. \right. \\ \left. \left. + \frac{\sin^2 \theta_4^k}{(\theta_4^k)^2} \left[ \cosh(2r_2^k) - \cos(2\theta_3^k + 2\arg\tau_k) \right] \right. \right. \\ \left. \left. \times \sinh(2r_2^k) \right] \right\} \end{aligned} \qquad (6.18) $$

$$ \begin{aligned} \mathrm{Cov}_{11,k}^{(\pi\pi)} = \frac{k}{2} \left( \left[ \cosh(2r_1^k) - \cos(2\theta_3^k + 2\theta_4^k) \sinh(2r_1^k) \right] \right. \\ \left. + |\tau_k|^2 \left\{ - \frac{\sin^2 \theta_4^k}{\theta_4^k} \cosh(2r_1^k) \right. \right. \\ \left. \left. \times \sinh(2r_2^k) \right] \right\} \end{aligned} $$
$$ -\frac{1}{2(\theta_4^k)^2} \left[ \begin{aligned} & \cos(2\theta_3^k) - \cos(2\theta_3^k + 2\theta_4^k) \\ & -2\theta_4^k \sin(2\theta_3^k + 2\theta_4^k) \end{aligned} \right] \sinh(2r_1^k) \\ & + \frac{\sin^2 \theta_4^k}{(\theta_4^k)^2} \left[ \begin{aligned} & \cosh(2r_2^k) + \cos(2\theta_3^k + 2\arg \tau_k) \\ & \times \sinh(2r_2^k) \end{aligned} \right] \end{aligned} \right) \quad (6.19) $$

$$ \mathrm{Cov}_{11,k}^{(\phi\pi)} = \frac{1}{2} \left( \begin{aligned} & -\sin(2\theta_3^k + 2\theta_4^k) \sinh(2r_1^k) \\ & + |\tau_k|^2 \left\{ \begin{aligned} & \frac{\sin^2 \theta_4^k}{(\theta_4^k)^2} \sin(2\theta_3^k + 2\arg \tau_k) \sinh(2r_2^k) \\ & + \frac{1}{2(\theta_4^k)^2} \left[ \begin{aligned} & -\sin(2\theta_3^k) + \sin(2\theta_3^k + 2\theta_4^k) \\ & -2\theta_4^k \cos(2\theta_3^k + 2\theta_4^k) \end{aligned} \right] \sinh(2r_1^k) \end{aligned} \right\} \end{aligned} \right). \quad (6.20) $$

In the limit where the two fields are uncoupled, $\tau_k = 0$, one recovers the result obtained for two-mode squeezed states in Ref. [13]. One can also check that in agreement with the discussion of Sect. 6.1, the leading correction to the power spectra is of quadratic order in $|\tau_k|$. Let us stress that those corrections involve parameters that describe the environment sector, such as $r_2^k$, hence observations carried out on the system alone can a priori lead to indirect information about the microphysical evolution of the traced over field(s). One should also note that, through the dynamical evolution, the presence of the interaction modifies all Bogolyubov coefficients, hence all squeezing and rotation parameters. This is why formally, in the above expressions, the leading terms (i.e. the ones before $|\tau_k|^2$) need also be expanded.

By plugging Eqs. (6.18), (6.19) and (6.20) into Eq. (6.16), one can finally derive an expression for the purity expanded at quadratic order in $|\tau_k|^2$, and by doing so one exactly recovers Eq. (6.9). This is an important consistency check as the two methods employed to derive this result are completely independent (and, as already argued, the approach based on the Wigner function is computationally less heavy).

This calculation also makes explicit that as one increases the interaction strength, one decreases the purity and hence makes decoherence more efficient, but one also induces larger corrections to the observable power spectra. This allows one to answer the question asked above, namely whether or not decoherence can proceed without affecting too much the power spectra. As noticed below Eq. (6.9), decoherence becomes more efficient as quantum squeezing increases, and in the large-squeezing limit, the correction to $\gamma_k = 1$ is controlled by $|\tau_k|^2 e^{2r}$. However, from Eqs. (6.18), (6.19) and (6.20), one can see that the relative correction to the power spectra is rather controlled by $|\tau_k|^2$ in that limit.

As a consequence, if the interaction strength is such that

$$ e^{-r} \ll |\tau_k| \ll 1, \quad (6.21) $$

decoherence takes place while keeping corrections to observable predictions tiny.

## 7 Conclusion

In this work, we have performed a detailed study of the quantum dynamics of two scalar fields, quadratically coupled, and embedded in a homogeneous and isotropic background. Their dynamics is generated by a quadratic Hamiltonian with time-dependent coefficients. Evolution of such systems (either classical or quantum) is obtained by applying elements of the symplectic group $Sp(4, \mathbb{R})$ to the initial configuration, so we have first investigated the mathematical structure of this group by presenting various descriptions of it. In particular, using the Bloch-Messiah decomposition that we further developed using the commutation relations of the Lie algebra, we have derived fully factorised expressions for the group elements of $Sp(4, \mathbb{R})$. The ten parameters entering these expressions are dubbed the squeezing and rotation parameters, as per the three parameters entering the decomposition of $Sp(2, \mathbb{R})$. Alternatively, this group is described by Bogolyubov coefficients, which we explicitly related to the squeezing and rotation parameters.

We then provided the quantum representation of the Lie algebra of $Sp(4, \mathbb{R})$, from which the quantum Hamiltonian can be easily expressed and interpreted. Couplings between the two fields manifest themselves through exchanges of quanta from one sector to the other (hence preserving the total number of quanta), and through direct productions of pairs in which each quantum belongs to a different sector. The latter provides a direct way to entangle the two sectors and those particles add up to the direct pair production occurring within each sector separately. The former also leads to entanglement but in an indirect way by transferring quanta which have been previously created in a given sector. Using this group-theoretic approach, we then showed that the evolution operator can be interpreted as the successive application of three blocks of quantum operations on the initial state. This sequence of operations schematically consists in first exchanging quanta between the two sectors, then creating pairs within each sector separately, and finally mixing again these newly-created quanta between the two sectors.

Applying the evolution operator to the vacuum state allowed us to derive the most general expression for the four-mode squeezed states, see Eqs. (5.6) and (5.7), which, to our knowledge, has not been presented in the literature so far. It can be viewed as the copy of two two-mode squeezed states, one for each sector, which then exchange quanta according to
their couplings. Its mathematical structure exhibits a power expansion in the coupling between the two fields, in which the power at each order gives the number of transfers between the two sectors. As an example, we provided explicit formulas for an expansion truncated at second order, i.e. including up to two exchanges of particles between the two fields. We finally described the four-mode squeezed states in terms of their Wigner functions, which were shown to be Gaussian. Their covariance is built from all the cross-spectra and we expressed them using either the Bogolyubov coefficients or the ten squeezing and rotation parameters.

In cases where one of the two sectors is unobserved (say the second sector, which we referred to as the “environment”), entanglement between the two fields leads to quantum decoherence in the first sector (dubbed the “system”), as well as modifications of its observable predictions. We studied this mechanism by first computing the reduced density matrix starting from the four-mode squeezed state. We showed that the environment induces corrections to observables that are necessarily of even power in the coupling strength. This is because, in order to preserve statistical isotropy, any particle transfer between the two sectors must be compensated by the inverse transfer of a particle with the same wavenumber, or by the transfer of a particle with opposite wavenumber, so the number of transfers is even. The purity of the system, $\gamma_k$, was also calculated from the reduced density matrix. We then investigated decoherence using the Wigner function and found that it substantially simplifies calculations. Indeed, we showed that tracing out the environment is readily obtained in the phase-space representation by marginalising the Wigner function over the phase-space of the environment. Since the Wigner function of a four-mode squeezed state is Gaussian, this operation is trivial. This allowed us to obtain a non-perturbative expression of the purity in terms of the power spectra of the system. We have finally expanded the result at second-order in the coupling parameters in both approaches (where we have checked that the same result is obtained).

The fact that the purity can be expressed in terms of the power spectra of the system entails that decoherence, i.e. the decrease of the purity, cannot proceed without affecting the observable power spectra. However, we have shown that in the large squeezing limit, there exists a regime, given by Eq. (6.21), where the interaction strength is large enough to make the system decohere but small enough to keep the observables mostly unchanged, shedding some light on the results of Refs. [23,55]. This also confirmed that squeezed states are more easily subject to decoherence [11].

Though limited to quadratic coupling, let us stress that our approach to decoherence does not rely on any approximation scheme, since the full quantum state of the joint system-plus-environment setup is first derived exactly, before tracing out the environment. It thus provides an ideal case study to

frame the range of applicability of approximate approaches to decoherence, such as the Lindblad formalism (at least in the simple situation of quadratic couplings). One may indeed compare the exact results obtained in this work with the ones derived from those effective methods, and this is the topic of a future work.

As mentioned in Sect. 1, our results are directly relevant for cosmology, in order to describe scalar fields in flat and non-flat Friedmann-Lemaître-Robertson-Walker geometries, quantum fields in curved spaces possibly with derivative couplings, and in the context of primordial cosmology, adiabatic and isocurvature perturbations in multiple-field inflation scenarios. But more generally, our results are relevant for any time-dependent, quadratic Hamiltonian that couples two degrees of freedom, regardless of the origin of these degrees of freedom. They thus offer a wide range of applications. In particular, the phase-space approach can be readily extended to cases where the environment is made of more than one scalar field. Suppose indeed that $N$ scalar fields, initially in their vacuum state, compose the environment. The Wigner function of the system-plus-environment setup is described by an $(N+1)$-dimensional Gaussian. Tracing out (i.e., in phase space, marginalising) over $N$ scalar fields is thus as trivial as tracing out over one single field. This may be used to understand how isocurvature modes can lead to the decoherence of the adiabatic sector, and should be compared with effective-field theory approaches for such systems [56].

Finally, since squeezed states feature quantum entanglement, they are an interesting playground to discuss possible setups for Bell and Leggett-Garg inequality violations in continuous systems, see Refs. [25, 27, 57, 58]. These analyses have been carried out for two-mode squeezed states, where only one type of entanglement can be harvested, namely the one between modes $\bar{k}$ and $-\bar{k}$ of the same field. Since the two subsystems $\bar{k}$ and $-\bar{k}$ are not locally distinct in real space, this necessarily restricts the analysis to those Bell inequalities where locality is not part of the assumptions being tested (such as, for instance, the Leggett-Garg or the temporal-Bell inequalities). The situation is however different for multiple-field systems, where one may chose to measure a field $\phi_1$ at position $\bar{x}_1$ and another field $\phi_2$ at position $\bar{x}_2$. This is because, on top of the entanglement between modes $\bar{k}$ and $-\bar{k}$ of the same field, one now has entanglement between quanta in different fields. This thus opens up the possibility to test for a wider class of Bell inequalities. Let us also mention that the present work would additionally allow one to assess how decoherence affects the ability to test Bell inequalities with squeezed states.

**Data Availability Statement** This manuscript has no associated data or the data will not be deposited. [Authors' comment: Data sharing not applicable to this article as no datasets were generated or analysed during the current study.]
**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Funded by SCOAP³.

## Appendix A: Commutation relations in the Lie algebra sp(4, R)

Following the procedure exemplified at the end of Sect. 3.1, which relies on combining the multiplication rule of the Kronecker product, Eq. (3.9), with the commutation relations between the Pauli matrices, Eq. (3.10), the commutators between the generators of sp(4, R), as listed in Table 1, are given by

$$ \text{Sq./Sq. : } [\mathbf{K}_1, \mathbf{K}_2] = 0; \qquad (\text{A.1}) $$

$$ \text{Rot./Rot. : } \begin{aligned} [\mathbf{K}_3, \mathbf{K}_5] &= 0, & [\mathbf{K}_3, \mathbf{K}_6] &= 0, \\ [\mathbf{K}_3, \mathbf{K}_4] &= 0, & [\mathbf{K}_5, \mathbf{K}_6] &= 2\mathbf{K}_4, & [\mathbf{K}_6, \mathbf{K}_4] &= 2\mathbf{K}_5, \\ [\mathbf{K}_4, \mathbf{K}_5] &= 2\mathbf{K}_6; & & & (\text{A.2}) \\ & & & & (\text{A.3}) \end{aligned} $$

$$ \text{Boost/Boost : } \begin{aligned} [\mathbf{K}_7, \mathbf{K}_8] &= 0, & [\mathbf{K}_7, \mathbf{K}_{10}] &= 0, \\ [\mathbf{K}_9, \mathbf{K}_8] &= 0, & [\mathbf{K}_9, \mathbf{K}_7] &= 2\mathbf{K}_6, & [\mathbf{K}_8, \mathbf{K}_{10}] &= 2\mathbf{K}_5, \\ [\mathbf{K}_9, \mathbf{K}_{10}] &= 2\mathbf{K}_3; & & & (\text{A.4}) \\ & & & & (\text{A.5}) \end{aligned} $$

$$ \text{Sq./Rot. : } \begin{aligned} [\mathbf{K}_1, \mathbf{K}_6] &= 0, & [\mathbf{K}_2, \mathbf{K}_5] &= 0, & (\text{A.6}) \\ [\mathbf{K}_1, \mathbf{K}_3] &= 2\mathbf{K}_8, & [\mathbf{K}_1, \mathbf{K}_4] &= 2\mathbf{K}_7, \\ [\mathbf{K}_1, \mathbf{K}_5] &= 2\mathbf{K}_9, & & & (\text{A.7}) \\ [\mathbf{K}_2, \mathbf{K}_3] &= 2\mathbf{K}_7, & [\mathbf{K}_2, \mathbf{K}_4] &= 2\mathbf{K}_8, \\ [\mathbf{K}_2, \mathbf{K}_6] &= 2\mathbf{K}_{10}; & & & (\text{A.8}) \end{aligned} $$

$$ \text{Sq./Boost : } \begin{aligned} [\mathbf{K}_1, \mathbf{K}_{10}] &= 0, & [\mathbf{K}_2, \mathbf{K}_9] &= 0, & (\text{A.9}) \\ [\mathbf{K}_1, \mathbf{K}_7] &= 2\mathbf{K}_4, & [\mathbf{K}_1, \mathbf{K}_8] &= 2\mathbf{K}_3, \\ [\mathbf{K}_1, \mathbf{K}_9] &= 2\mathbf{K}_5, & & & (\text{A.10}) \\ [\mathbf{K}_2, \mathbf{K}_7] &= 2\mathbf{K}_3, & [\mathbf{K}_2, \mathbf{K}_8] &= 2\mathbf{K}_4, \\ [\mathbf{K}_2, \mathbf{K}_{10}] &= 2\mathbf{K}_6; & & & (\text{A.11}) \end{aligned} $$

$$ \text{Rot./Boost : } \begin{aligned} [\mathbf{K}_4, \mathbf{K}_9] &= 0, & [\mathbf{K}_4, \mathbf{K}_{10}] &= 0, & (\text{A.12}) \\ [\mathbf{K}_5, \mathbf{K}_7] &= 0, & [\mathbf{K}_6, \mathbf{K}_8] &= 0, & (\text{A.13}) \\ [\mathbf{K}_3, \mathbf{K}_7] &= 2\mathbf{K}_2, & & & (\text{A.14}) \\ [\mathbf{K}_3, \mathbf{K}_8] &= 2\mathbf{K}_1, & & & \\ [\mathbf{K}_4, \mathbf{K}_7] &= 2\mathbf{K}_1, & & & \end{aligned} $$

$$ [\mathbf{K}_4, \mathbf{K}_8] = 2\mathbf{K}_2, \qquad (\text{A.15}) $$

$$ \begin{aligned} [\mathbf{K}_5, \mathbf{K}_9] &= 2\mathbf{K}_1, \\ [\mathbf{K}_6, \mathbf{K}_{10}] &= 2\mathbf{K}_2, & (\text{A.16}) \end{aligned} $$

$$ \begin{aligned} [\mathbf{K}_3, \mathbf{K}_{10}] &= 2\mathbf{K}_9, \\ [\mathbf{K}_9, \mathbf{K}_3] &= 2\mathbf{K}_{10}, & (\text{A.17}) \end{aligned} $$

$$ \begin{aligned} [\mathbf{K}_5, \mathbf{K}_{10}] &= 2\mathbf{K}_8, \\ [\mathbf{K}_8, \mathbf{K}_5] &= 2\mathbf{K}_{10}, & (\text{A.18}) \end{aligned} $$

$$ \begin{aligned} [\mathbf{K}_6, \mathbf{K}_7] &= 2\mathbf{K}_9, \\ [\mathbf{K}_9, \mathbf{K}_6] &= 2\mathbf{K}_7. & (\text{A.19}) \end{aligned} $$

One can identify various subalgebras, which is particularly useful when it comes to factorising down elements of the group. Since all subalgebras are three dimensional, we can use Bianchi classification to sort them. In the following equation, each line corresponds to a subalgebra, i.e. a set of closed generators by the adjoint operation:

$$ \text{Type IX : } \begin{aligned} [\mathbf{K}_5, \mathbf{K}_6] &= 2\mathbf{K}_4, & [\mathbf{K}_6, \mathbf{K}_4] &= 2\mathbf{K}_5, \\ [\mathbf{K}_4, \mathbf{K}_5] &= 2\mathbf{K}_6; & & & (\text{A.20}) \end{aligned} $$

$$ \text{Type VIII : } \begin{aligned} [\mathbf{K}_1, \mathbf{K}_3] &= 2\mathbf{K}_8, & [\mathbf{K}_3, \mathbf{K}_8] &= 2\mathbf{K}_1, \\ [\mathbf{K}_8, \mathbf{K}_1] &= -2\mathbf{K}_3; & & & (\text{A.21}) \\ [\mathbf{K}_1, \mathbf{K}_4] &= 2\mathbf{K}_7, & [\mathbf{K}_4, \mathbf{K}_7] &= 2\mathbf{K}_1, \\ [\mathbf{K}_7, \mathbf{K}_1] &= -2\mathbf{K}_4; & & & (\text{A.22}) \\ [\mathbf{K}_1, \mathbf{K}_5] &= 2\mathbf{K}_9, & [\mathbf{K}_5, \mathbf{K}_9] &= 2\mathbf{K}_1, \\ [\mathbf{K}_9, \mathbf{K}_1] &= -2\mathbf{K}_5; & & & (\text{A.23}) \\ [\mathbf{K}_2, \mathbf{K}_3] &= 2\mathbf{K}_7, & [\mathbf{K}_3, \mathbf{K}_7] &= 2\mathbf{K}_2, \\ [\mathbf{K}_7, \mathbf{K}_2] &= -2\mathbf{K}_3; & & & (\text{A.24}) \\ [\mathbf{K}_2, \mathbf{K}_4] &= 2\mathbf{K}_8, & [\mathbf{K}_4, \mathbf{K}_8] &= 2\mathbf{K}_2, \\ [\mathbf{K}_8, \mathbf{K}_2] &= -2\mathbf{K}_4; & & & (\text{A.25}) \\ [\mathbf{K}_2, \mathbf{K}_6] &= 2\mathbf{K}_{10}, & [\mathbf{K}_6, \mathbf{K}_{10}] &= 2\mathbf{K}_2, \\ [\mathbf{K}_{10}, \mathbf{K}_2] &= -2\mathbf{K}_6; & & & (\text{A.26}) \\ [\mathbf{K}_9, \mathbf{K}_6] &= 2\mathbf{K}_7, & [\mathbf{K}_6, \mathbf{K}_7] &= 2\mathbf{K}_9, \\ [\mathbf{K}_7, \mathbf{K}_9] &= -2\mathbf{K}_6; & & & (\text{A.27}) \\ [\mathbf{K}_8, \mathbf{K}_5] &= 2\mathbf{K}_{10}, & [\mathbf{K}_5, \mathbf{K}_{10}] &= 2\mathbf{K}_8, \\ [\mathbf{K}_{10}, \mathbf{K}_8] &= -2\mathbf{K}_5; & & & (\text{A.28}) \\ [\mathbf{K}_9, \mathbf{K}_3] &= 2\mathbf{K}_{10}, & [\mathbf{K}_3, \mathbf{K}_{10}] &= 2\mathbf{K}_9, \\ [\mathbf{K}_{10}, \mathbf{K}_9] &= -2\mathbf{K}_3. & & & (\text{A.29}) \end{aligned} $$

The type IX Bianchi algebra is related to a su(2) subalgebra that contains three rotation generators while the type VIII Bianchi algebra is isomorphic to sl(2, R) subalgebras. The first six type-VIII subalgebras contain a rotation, a squeezing and a boost generator, so they can be related to sp(2, R) ≅ sl(2, R), the last three type-VIII subalgebras are made of a rotation and two boost generators, making their interpretation less obvious.
## Appendix B: Evolved vacuum state

This appendix presents the computation that leads to the expression of the evolved vacuum state in the occupation-number representation given in Eqs. (5.6) and (5.7). Starting from the initial vacuum state (5.3), we follow the different operations displayed in the quantum circuits of Sect. 4.2. Let us first consider $\hat{\mathcal{R}}_{\bar{k}}(\varphi_k)$, which for convenience we reproduce here:

<figure><img src="image_4.png" /></figure>

Since the operators contained in $\hat{\mathcal{R}}_{\bar{k}}(\varphi_k)$ do not create particles, the components $\hat{\mathcal{R}}_{2 \to 1}$ and $\hat{\mathcal{R}}_{1 \to 2}$ leave the initial vacuum state unchanged. Moreover, the operator $\hat{\mathcal{R}}_1(-iq_z^k/2)$ generates an overall factor $e^{-iq_z^k/2}$, which is exactly compensated by the operator $\hat{\mathcal{R}}_2(iq_z^k/2)$ that generates an overall factor $e^{iq_z^k/2}$. As a consequence, only the overall phase shift controlled by $\varphi_3^k$ remains, and the evolved vacuum state can be obtained as

$$ \hat{\mathcal{Z}}_2(r_2^k) [\chi_{\bar{k}}^{(1)} \otimes |0_{\bar{k}}^{(2)}, 0_{-\bar{k}}^{(2)}\rangle] \quad (B.2) $$

$$ = \chi_{\bar{k}}^{(1)} \otimes \frac{1}{\cosh r_2^k} \sum_{m=0}^{\infty} (-1)^m \tanh^m r_2^k |m_{\bar{k}}^{(2)}, m_{-\bar{k}}^{(2)}\rangle, \quad (B.3) $$

where $\chi_{\bar{k}}^{(i)}$ is any vector belonging to $\mathcal{E}_{\bar{k}}^{(i)} \otimes \mathcal{E}_{-\bar{k}}^{(i)}$ for $i = 1, 2$ and the squeezing parameters $r_1^k$ and $r_2^k$ control the two-mode creation in each sector. At this stage the state is thus given by

<figure><img src="image_5.png" /></figure>

where we recall that $\hat{\mathcal{R}}_i$, $\hat{\mathcal{Z}}_i$ and $\hat{\mathcal{R}}_{i \to j}$ with $i, j = 1, 2$ are defined in Eqs. (4.16), (4.17) and (4.18). Let us see how these operators act one after the other.

The operators $\hat{\mathcal{R}}_1(\varphi_3^k)$ and $\hat{\mathcal{R}}_2(\varphi_3^k)$ simply add a global phase factor $e^{-2i\varphi_3^k}$, so the initial vacuum state is first transformed according to

$$ e^{-2i\varphi_3^k} |0_{\bar{k}}^{(1)}, 0_{-\bar{k}}^{(1)}, 0_{\bar{k}}^{(2)}, 0_{-\bar{k}}^{(2)}\rangle. \quad (B.1) $$

To derive the action of the squeezing operators $\hat{\mathcal{Z}}_i$, we recall that the two squeezing generators commute, and that they act on each sector separately. As a consequence, we can use the result for two-mode squeezed states in Sp(2, $\mathbb{R}$), as derived e.g. in Ref. [13], and write

$$ \hat{\mathcal{Z}}_1(r_1^k) [0_{\bar{k}}^{(1)}, 0_{-\bar{k}}^{(1)} \otimes \chi_{\bar{k}}^{(2)}] = \frac{1}{\cosh r_1^k} \sum_{n=0}^{\infty} (-1)^n \tanh^n r_1^k |n_{\bar{k}}^{(1)}, n_{-\bar{k}}^{(1)}\rangle \otimes \chi_{\bar{k}}^{(2)}, $$

$$ \frac{e^{-2i\varphi_3^k}}{\cosh r_1^k \cosh r_2^k} \sum_{n,m=0}^{\infty} (-1)^{n+m} \tanh^n r_1^k \tanh^m r_2^k |n_{\bar{k}}^{(1)}, n_{-\bar{k}}^{(1)}, m_{\bar{k}}^{(2)}, m_{-\bar{k}}^{(2)}\rangle. \quad (B.4) $$

Then comes the contribution from $\hat{\mathcal{R}}_i(\theta_3^k)$, which simply involves number counting operators, see Eq. (4.16), so the state becomes

$$ \frac{e^{-2i(\varphi_3^k + \theta_3^k)}}{\cosh r_1^k \cosh r_2^k} \sum_{n,m=0}^{\infty} (-1)^{n+m} e^{-2i\theta_3^k(n+m)} \tanh^n r_1^k \tanh^m r_2^k \\ \times |n_{\bar{k}}^{(1)}, n_{-\bar{k}}^{(1)}, m_{\bar{k}}^{(2)}, m_{-\bar{k}}^{(2)}\rangle. \quad (B.5) $$

The action of $\hat{\mathcal{R}}_{2 \to 1}$ is more involved since it transfers excitations from one sector to the other. We first note that, in $\hat{\mathcal{R}}_{2 \to 1}$ and $\hat{\mathcal{R}}_{1 \to 2}$, the domains $\bar{k}$ and $-\bar{k}$ can be factorised out since $[\hat{a}_{1,\bar{k}^*}^\dagger \hat{a}_{2,\bar{k}}^\dagger, \hat{a}_{1,-\bar{k}}^\dagger \hat{a}_{2,-\bar{k}}^\dagger] = [\hat{a}_{2,\bar{k}}^\dagger \hat{a}_{1,\bar{k}}^\dagger, \hat{a}_{2,-\bar{k}}^\dagger \hat{a}_{1,-\bar{k}}^\dagger] = 0$.
This implies that

$$
\begin{align}
\hat{\mathcal{R}}_{2 \to 1}(-p_-^k) &= \exp \left[ -ip_-^k (\hat{a}_{1,\vec{k}}^\dagger \hat{a}_{2,\vec{k}}) \right] \nonumber \\
&\quad \times \exp \left[ -ip_-^k (\hat{a}_{1,-\vec{k}}^\dagger \hat{a}_{2,-\vec{k}}) \right] \tag{B.6} \\
&= \left[ \sum_{i=0}^{\infty} \frac{(-ip_-^k)^i}{i!} (\hat{a}_{1,\vec{k}}^\dagger \hat{a}_{2,\vec{k}})^i \right] \nonumber \\
&\quad \times \left[ \sum_{j=0}^{\infty} \frac{(-ip_-^k)^j}{j!} (\hat{a}_{1,-\vec{k}}^\dagger \hat{a}_{2,-\vec{k}})^j \right], \tag{B.7}
\end{align}
$$

where the exponentials have been Taylor expanded, and a
similar expression for $\hat{\mathcal{R}}_{1 \to 2}(p_+^k)$ can be written down for
future use, namely

$$
\begin{equation}
\begin{split}
\hat{\mathcal{R}}_{1 \to 2}(p_+^k) = & \left[ \sum_{k=0}^{\infty} \frac{(ip_+^k)^k}{k!} (\hat{a}_{2,\vec{k}}^\dagger \hat{a}_{1,\vec{k}})^k \right] \\
& \times \left[ \sum_{\ell=0}^{\infty} \frac{(ip_+^k)^{\ell}}{\ell!} (\hat{a}_{2,-\vec{k}}^\dagger \hat{a}_{1,-\vec{k}})^{\ell} \right].
\end{split}
\tag{B.8}
\end{equation}
$$

The action of $\hat{a}_{i,\pm\vec{k}}^\dagger \hat{a}_{j\mp\vec{k}}$ is to transfer an excitation from sector $j$ to sector $i$, so upon applying Eq. (B.7) onto Eq. (B.5), with the schematic normalisation $\hat{a}|n\rangle = \sqrt{n}|n-1\rangle$ and $\hat{a}^\dagger|n\rangle = \sqrt{n+1}|n+1\rangle$, the state becomes

$$
\begin{equation}
\begin{split}
& \frac{e^{-2i(\phi_3^k + \theta_3^k)}}{\cosh r_1^k \cosh r_2^k} \sum_{n,m=0}^{\infty} (-1)^{n+m} e^{-2i\theta_3^k(n+m)} \\
& \quad \times \tanh^n r_1^k \tanh^m r_2^k \sum_{i=0}^m \frac{(-ip_-^k)^i}{i!} \sum_{j=0}^m \frac{(-ip_-^k)^j}{j!} \\
& \quad \times \sqrt{\frac{m!}{(m-i)!}} \sqrt{\frac{m!}{(m-j)!}} \sqrt{\frac{(n+i)!}{n!}} \sqrt{\frac{(n+j)!}{n!}} \\
& \quad \times |(n+i)_{\vec{k}}^{(1)}, (n+j)_{-\vec{k}}^{(1)}, (m-i)_{\vec{k}}^{(2)}, (m-j)_{-\vec{k}}^{(2)}\rangle.
\end{split}
\tag{B.9}
\end{equation}
$$

The next step is to apply $\hat{\mathcal{R}}_1(-\frac{i}{2}p_z^k)$ and $\hat{\mathcal{R}}_2(\frac{i}{2}p_z^k)$, which add imaginary phases (i.e. exponential modulation) to each term that depend on their number of particles, and the state (B.9) becomes

$$
\begin{equation}
\begin{split}
& \frac{e^{-2i(\phi_3^k + \theta_3^k)}}{\cosh r_1^k \cosh r_2^k} \sum_{n,m=0}^{\infty} (-1)^{n+m} e^{-2i\theta_3^k(n+m)} \tanh^n r_1^k \tanh^m r_2^k \\
& \quad \times \sum_{i=0}^{m} \frac{(-ip_-^k)^i}{i!} \sum_{j=0}^{m} \frac{(-ip_-^k)^j}{j!} \\
& \quad \times \sqrt{\frac{m!}{(m-i)!}} \sqrt{\frac{m!}{(m-j)!}} \sqrt{\frac{(n+i)!}{n!}} \sqrt{\frac{(n+j)!}{n!}} \\
& \quad \times e^{-\frac{p_z^k}{2}(2n+i+j)} e^{\frac{p_z^k}{2}(2m-i-j)} \\
& \quad \times |(n+i)_{\vec{k}}^{(1)}, (n+j)_{-\vec{k}}^{(1)}, (m-i)_{\vec{k}}^{(2)}, (m-j)_{-\vec{k}}^{(2)}\rangle.
\end{split}
\tag{B.10}
\end{equation}
$$

Finally, the application of $\hat{\mathcal{R}}_{1 \to 2}(p_+^k)$ can be done using Eq. (B.8), and one obtains

$$
\begin{equation}
\begin{split}
& \frac{e^{-2i(\phi_3^k + \theta_3^k)}}{\cosh r_1^k \cosh r_2^k} \sum_{n,m=0}^{\infty} (-1)^{n+m} e^{-2i\theta_3^k(n+m)} \\
& \quad \times \tanh^n r_1^k \tanh^m r_2^k \sum_{i=0}^{m} \frac{(-ip_-^k)^i}{i!} \sum_{j=0}^{m} \frac{(-ip_-^k)^j}{j!} \\
& \quad \times \sqrt{\frac{m!}{(m-i)!}} \sqrt{\frac{m!}{(m-j)!}} \sqrt{\frac{(n+i)!}{n!}} \sqrt{\frac{(n+j)!}{n!}} \\
& \quad \times e^{-\frac{p_z^k}{2}(2n+i+j)} e^{\frac{p_z^k}{2}(2m-i-j)} \\
& \quad \times \sum_{k=0}^{n+i} \frac{(ip_+^k)^k}{k!} \sum_{\ell=0}^{n+j} \frac{(ip_+^k)^{\ell}}{\ell!} \sqrt{\frac{(n+i)!}{(n+i-k)!}} \sqrt{\frac{(n+j)!}{(n+j-\ell)!}} \\
& \quad \times \sqrt{\frac{(m-i+k)!}{(m-i)!}} \sqrt{\frac{(m-j+k)!}{(m-j)!}} \\
& \quad \times |(n+i-k)_{\vec{k}}^{(1)}, (n+j-\ell)_{-\vec{k}}^{(1)}, (m-i+k)_{\vec{k}}^{(2)}, \\
& \quad \times (m-j+\ell)_{-\vec{k}}^{(2)}\rangle.
\end{split}
\tag{B.11}
\end{equation}
$$

This expression can be slightly simplified by replacing the
sum over $k$ and $\ell$ with a sum over $s \equiv i - k$ and $t = j - \ell$,
and one obtains

$$
\begin{equation}
\begin{split}
|\phi_{\vec{k}}(t)\rangle ={}& \frac{e^{-2i(\phi_3^k + \theta_3^k)}}{\cosh r_1^k \cosh r_2^k} \\
& \times \sum_{n,m=0}^{\infty} (-1)^{n+m} e^{-2i\theta_3^k(n+m)} e^{p_z^k(m-n)} \\
& \times \tanh^n r_1^k \tanh^m r_2^k \\
& \times \sum_{i,j=0}^{m} \frac{(-ip_-^k e^{-p_z^k})^{i+j}}{i!j!} \frac{(n+i)!(n+j)!m!}{(m-i)!(m-j)!n!} \\
& \times \sum_{s=-n}^{i} \sum_{t=-n}^{j} \frac{(ip_+^k)^{i+j-s-t}}{(i-s)!(j-t)!} \sqrt{\frac{(m-s)!(m-t)!}{(n+s)!(n+t)!}} \\
& \times |(n+s)_{\vec{k}}^{(1)}, (n+t)_{-\vec{k}}^{(1)}, (m-s)_{\vec{k}}^{(2)}, (m-t)_{-\vec{k}}^{(2)}\rangle.
\end{split}
\tag{B.12}
\end{equation}
$$

In order to first sum over the indices appearing in the num-
ber of particle eigenstates, one can flip the ordering of
the sums over $i, j$ and $s, t$, using that $\sum_{i=0}^{m} \sum_{s=-n}^{i} =$
$\sum_{s=-n}^{m} \sum_{i=\max(s,0)}^{m}$ and $\sum_{j=0}^{m} \sum_{t=-n}^{j} =$
$\sum_{t=-n}^{m} \sum_{j=\max(t,0)}^{m}$, leading to

$$
|\phi_{\vec{k}}(t)\rangle = \frac{e^{-2i(\phi_3^k + \theta_3^k)}}{\cosh r_1^k \cosh r_2^k} \times \sum_{n,m=0}^{\infty} (-1)^{n+m} e^{-2i\theta_3^k(n+m)} e^{p_z^k(m-n)}
$$
$$ 
\begin{aligned}
& \times \tanh^n r_1^k \tanh^m r_2^k \frac{m!}{n!} \\
& \times \sum_{s,t=-n}^{m} \sqrt{\frac{(m-s)!(m-t)!}{(n+s)!(n+t)!}} (ip_+^k)^{-s-t} \\
& \times \sum_{i=\max(s,0)}^{m} \sum_{j=\max(t,0)}^{m} \frac{(p_-^k p_+^k e^{-p_z^k})^{i+j}}{i!j!(i-s)!(j-t)!} \\
& \times \frac{(n+i)!(n+j)!}{(m-i)!(m-j)!} \\
& \times |(n+s)_{\bar{k}}^{(1)}, (n+t)_{-\bar{k}}^{(1)}, (m-s)_{\bar{k}}^{(2)}, (m-t)_{-\bar{k}}^{(2)}\rangle.
\end{aligned}
\quad (B.13) $$

This is the result presented in Eqs. (5.6) and (5.7) in the main text.

## Appendix C: Reduced density matrix

In this appendix, we explicitly trace out the environmental degrees of freedom in Fock space as discussed in Sect. 6.1, leading to the expression for the reduced density matrix and the purity given in Eqs. (6.4) and (6.6). For the evolved vacuum state given in Eq. (5.6), the density matrix $\hat{\rho}_{\bar{k}}(t) = |\Phi_{\bar{k}}(t)\rangle\langle\Phi_{\bar{k}}(t)|$ reads

$$ 
\begin{aligned}
\hat{\rho}_{\bar{k}} = & \sum_{n,m=0}^{\infty} \sum_{s,t=-n}^{m} \sum_{n',m'=0}^{\infty} \sum_{s',t'=-n'}^{m'} \\
& \times c_k(n, m, s, t) c_k^*(n', m', s', t') \\
& \times |(n+s)_{\bar{k}}^{(1)}, (n+t)_{-\bar{k}}^{(1)}, (m-s)_{\bar{k}}^{(2)}, (m-t)_{-\bar{k}}^{(2)}\rangle \\
& \times \langle(n'+s')_{\bar{k}}^{(1)}, (n'+t')_{-\bar{k}}^{(1)}, (m'-s')_{\bar{k}}^{(2)}, (m'-t')_{-\bar{k}}^{(2)}|.
\end{aligned}
\quad (C.1) $$

Following Eq. (6.1), one can trace out the environmental degrees of freedom in the Fock basis. When doing so, the only non-vanishing terms are such that $m - s = m' - s'$ and $m - t = m' - t'$. This allows one to fix $s'$ and $t'$, and since the conditions $-n' \le s'$, $t' \le m'$ imposed by the sum boundaries in Eq. (C.1) imply that $m - m' - n' \le s, t \le m$, the reduced density matrix reads

$$ 
\begin{aligned}
\hat{\rho}_{\bar{k}, \text{red}}(t) = & \sum_{n,m,n',m'=0}^{\infty} \sum_{s,t=\max(-n,m-m'-n')}^{m} \\
& \times c_k(n, m, s, t) c_k^* \\
& \times (n', m', s+m'-m, t+m'-m) \\
& \times |(n+s)_{\bar{k}}^{(1)}, (n+t)_{-\bar{k}}^{(1)}\rangle \\
& \times \langle(n'+s+m'-m)_{\bar{k}}^{(1)}, (n'+t+m'-m)_{-\bar{k}}^{(1)}|.
\end{aligned}
\quad (C.2) $$

Let us now replace $n'$ and $m'$ by the new indices $N = n' + m' - m$ and $M = m' - m$. This gives rise to

$$ 
\begin{aligned}
\hat{\rho}_{\bar{k}, \text{red}}(t) = & \sum_{n,m=0}^{\infty} \sum_{M=-m}^{\infty} \sum_{N=M}^{\infty} \sum_{s,t=\max(-n,-N)}^{m} \\
& \times c_k(n, m, s, t) \\
& \times c_k^*(N-M, m+M, s+M, t+M) \\
& \times |(n+s)_{\bar{k}}^{(1)}, \\
& (n+t)_{-\bar{k}}^{(1)}\rangle\langle(N+s)_{\bar{k}}^{(1)}, (N+t)_{-\bar{k}}^{(1)}|.
\end{aligned}
\quad (C.3) $$

One can note that the indices $m$ and $M$ do not appear explicitly in the elements of the Fock basis, which is the reason why we have performed the above change of indices. In order to rewrite the sums over $m$ and $M$ as internal sums, similarly to what was done between Eqs. (B.12) and (B.13), one can re-order the various indices and write

$$ 
\begin{aligned}
\hat{\rho}_{\bar{k}, \text{red}}(t) = & \sum_{n=0}^{\infty} \sum_{N=-\infty}^{\infty} \sum_{s,t=-\min(n,N)}^{\infty} \\
& \times |(n+s)_{\bar{k}}^{(1)}, (n+t)_{-\bar{k}}^{(1)}\rangle\langle(N+s)_{\bar{k}}^{(1)}, (N+t)_{-\bar{k}}^{(1)}| \\
& \times \underbrace{\sum_{m=\max(0,s,t)}^{\infty} \sum_{M=-m}^{N} c_k(n, m, s, t) c_k^*(N-M, m+M, s+M, t+M)}_{\Xi_k(n,N,s,t)}
\end{aligned}
\quad (C.4) $$

which defines the coefficients $\Xi_k$, and which matches Eqs. (6.4) and (6.5) given in the main text (where the indices have been renamed for notational convenience).

From this expression, the purity of the system can also be computed. Squaring Eq. (C.4), one obtains

$$ 
\begin{aligned}
\hat{\rho}_{\bar{k}, \text{red}}^2 = & \sum_{n,n'=0}^{\infty} \sum_{N,N'=-\infty}^{\infty} \sum_{s,t=-\min(n,N)}^{\infty} \\
& \times \Xi_k(n, N, s, t) \Xi_k \\
& \times (n', N', s+N-n', t+N-n') \\
& \times |(n+s)_{\bar{k}}^{(1)}, (n+t)_{-\bar{k}}^{(1)}\rangle \\
& \times \langle(N'+N-n'+s)_{\bar{k}}^{(1)}, (N'+N-n'+t)_{-\bar{k}}^{(1)}|,
\end{aligned}
\quad (C.5) $$

from which the purity (6.2) can be expressed as

$$ 
\begin{aligned}
\gamma_{\bar{k}}(t) = & \sum_{n=0}^{\infty} \sum_{N,N'=-\infty}^{\infty} \sum_{s,t=-\min(n,N)}^{\infty} \Xi_k(n, N, s, t) \\
& \times \Xi_k(N+N'-n, N', s+n-N', t+n-N').
\end{aligned}
\quad (C.6) $$

This expression can be slightly simplified by replacing the sum over $N'$ by a sum over $u \equiv n - N'$, and one obtains

$$ 
\gamma_{\bar{k}}(t) = \sum_{n=0}^{\infty} \sum_{N,u=-\infty}^{\infty} \sum_{s,t=-\min(n,N)}^{\infty} \Xi_k(n, N, s, t)
$$
$$ \times \Xi_k(N-u, n-u, s+u, t+u). \qquad (C.7) $$

This matches Eq. (6.6) given in the main text, where the indices have been slightly renamed for notational convenience.

## Appendix D: Equivalence between tracing out in Hilbert space and marginalisation in the phase space

In Sect. 6.2, we have shown that a partial trace of the density matrix in Hilbert space is associated to a partial integration over the corresponding degrees of freedom of the Wigner function in phase space. This proof was however indirect, and relied on the equivalence between different ways of computing the expectation value of quantum operators. In this appendix, for completeness, we provide a direct proof that the reduced Wigner function is nothing but the Wigner–Weyl transform of the reduced density matrix.

For Hermitian operators, the Wigner–Weyl transform introduced in Eq. (5.39) can be inverted, and the classical phase-space function $\tilde{O}_{\vec{k}}(\mathbf{q}_{\vec{k}})$ gives rise to the quantum operator $\hat{O}_{\vec{k}}$ through

$$ \hat{O}_{\vec{k}} = \frac{1}{(2\pi)^4} \int_{\mathbb{R}^8} d^8 \zeta_{\vec{k}} \hat{D}(\zeta_{\vec{k}}) \times \int_{\mathbb{R}^8} d^8 \mathbf{q}_{\vec{k}} \exp (i \zeta_{\vec{k}} \cdot \mathbf{q}_{\vec{k}}) \tilde{O}_{\vec{k}}(\mathbf{q}_{\vec{k}}). \qquad (D.1) $$

It is obtained by integrating the four-mode displacement operator $\hat{D}(\zeta_{\vec{k}})$ against the Fourier transform of $\tilde{O}_{\vec{k}}(\mathbf{q}_{\vec{k}})$, where $\zeta_{\vec{k}}$ is the conjugate vector to $\mathbf{q}_{\vec{k}}$, expanded as

$$ \zeta_{\vec{k}} = (\zeta_{1,\vec{k}}, \zeta_{1,-\vec{k}}, \zeta_{2,\vec{k}}, \zeta_{2,-\vec{k}}, \kappa_{1,\vec{k}}, \kappa_{1,-\vec{k}}, \kappa_{2,\vec{k}}, \kappa_{2,-\vec{k}})^{\mathrm{T}}, \qquad (D.2) $$

and the displacement operator $\hat{D}(\zeta_{\vec{k}})$ is defined as

$$ \hat{D}(\zeta_{\vec{k}}) = \hat{D}_{1,\vec{k}}(\zeta_{\vec{k}}) \hat{D}_{1,-\vec{k}}(\zeta_{\vec{k}}) \hat{D}_{2,\vec{k}}(\zeta_{\vec{k}}) \hat{D}_{2,-\vec{k}}(\zeta_{\vec{k}}), \qquad (D.3) $$

where $\hat{D}_{i,\vec{k}}(\zeta_{\vec{k}})$ are the one-mode displacement operators

$$ \hat{D}_{i,\pm\vec{k}}(\zeta_{\vec{k}}) = \exp \left( \gamma_{i,\pm\vec{k}} \hat{a}_{i,\pm\vec{k}}^{\dagger} - \gamma_{i,\pm\vec{k}}^{*} \hat{a}_{i,\pm\vec{k}} \right) \qquad (D.4) $$

for $i = 1, 2$ and

$$ \gamma_{i,\pm\vec{k}} = -\frac{1}{\sqrt{2}} (\kappa_{i,\pm\vec{k}} + i\zeta_{i,\pm\vec{k}}). \qquad (D.5) $$

The reason why $\hat{D}(\zeta_{\vec{k}})$ can be factorised in the form (D.3) is because the creation and annihilation operators commute across the four modes.

Since the Wigner function $W_{\vec{k}}(\mathbf{q}_{\vec{k}}, t)$ is the Wigner–Weyl transform of the density matrix $\hat{\rho}_{\vec{k}}(t)$, one can use Eq. (D.1)

to extract the density matrix from the Wigner function,

$$ \hat{\rho}_{\vec{k}}(t) = \frac{1}{(2\pi)^4} \int_{\mathbb{R}^8} d^8 \zeta_{\vec{k}} \hat{D}(\zeta_{\vec{k}}) \int_{\mathbb{R}^8} d^8 \mathbf{q}_{\vec{k}} \times \exp (i \zeta_{\vec{k}} \cdot \mathbf{q}_{\vec{k}}) W_{\vec{k}}(\mathbf{q}_{\vec{k}}, t). \qquad (D.6) $$

Similarly, the reduced density matrix must be connected to the reduced Wigner function (that we here aim at determining) through a relation of the form

$$ \hat{\rho}_{\vec{k},\mathrm{red}}(t) = \frac{1}{(2\pi)^2} \int_{\mathbb{R}^4} d^4 \zeta_{1,\vec{k}} \hat{D}_1(\zeta_{1,\vec{k}}) \int_{\mathbb{R}^4} d^4 \mathbf{q}_{1,\vec{k}} \times \exp (i \mathbf{q}_{1,\vec{k}} \cdot \zeta_{1,\vec{k}}) W_{\vec{k},\mathrm{red}}(\mathbf{q}_{1,\vec{k}}, t), \qquad (D.7) $$

where $\hat{D}_1(\zeta_{1,\vec{k}}) = \hat{D}_{1,\vec{k}}(\zeta_{1,\vec{k}}) \hat{D}_{1,-\vec{k}}(\zeta_{1,\vec{k}})$ is the first-sector displacement operator. Let us recall that $\hat{\rho}_{\vec{k},\mathrm{red}}$ is defined in Eq. (6.1) as the partial trace of the density matrix $\hat{\rho}_{\vec{k}}$ over the environmental degrees of freedom. By plugging Eq. (D.6) into Eq. (6.1), and expanding $\hat{D}(\zeta_{\vec{k}})$ as in Eq. (D.3), one finds

$$ \hat{\rho}_{\vec{k},\mathrm{red}}(t) = \frac{1}{(2\pi)^2} \int_{\mathbb{R}^4} d^4 \zeta_{1,\vec{k}} \hat{D}_1(\zeta_{1,\vec{k}}) \times \int_{\mathbb{R}^4} d^4 \mathbf{q}_{1,\vec{k}} \exp (i \mathbf{q}_{1,\vec{k}} \cdot \zeta_{1,\vec{k}}) \times \left[ \frac{1}{(2\pi)^2} \int_{\mathbb{R}^4} d^4 \zeta_{2,\vec{k}} \sum_{u,v=0}^{\infty} \right. \\ \times \langle u_{\vec{k}}^{(2)}, v_{-\vec{k}}^{(2)}(t) | \hat{D}_2(\zeta_{2,\vec{k}}) | u_{\vec{k}}^{(2)}, v_{-\vec{k}}^{(2)}(t) \rangle \\ \times \left. \int_{\mathbb{R}^4} d^4 \mathbf{q}_{2,\vec{k}} \exp (i \mathbf{q}_{2,\vec{k}} \cdot \zeta_{2,\vec{k}}) W_{\vec{k}}(\mathbf{q}_{1,\vec{k}}, \mathbf{q}_{2,\vec{k}}, t) \right]. \qquad (D.8) $$

By comparing this expression with Eq. (D.7), one can read off

$$ W_{\vec{k},\mathrm{red}}(\mathbf{q}_{1,\vec{k}}, t) = \frac{1}{(2\pi)^2} \int_{\mathbb{R}^4} d^4 \mathbf{q}_{2,\vec{k}} W_{\vec{k}}(\mathbf{q}_{1,\vec{k}}, \mathbf{q}_{2,\vec{k}}, t) \mathcal{N}(\mathbf{q}_{2,\vec{k}}), \qquad (D.9) $$

where

$$ \mathcal{N}(\mathbf{q}_{2,\vec{k}}) = \int_{\mathbb{R}^4} d^4 \zeta_{2,\vec{k}} \exp (i \mathbf{q}_{2,\vec{k}} \cdot \zeta_{2,\vec{k}}) \times \sum_{u,v=0}^{\infty} \langle u_{\vec{k}}^{(2)}, v_{-\vec{k}}^{(2)}(t) | \hat{D}_2(\zeta_{2,\vec{k}}) | u_{\vec{k}}^{(2)}, v_{-\vec{k}}^{(2)}(t) \rangle. \qquad (D.10) $$

To evaluate $\mathcal{N}(\mathbf{q}_{2,\vec{k}})$, we first use the factorisation (D.3) and compute

$$ \sum_{u,v=0}^{\infty} \langle u_{\vec{k}}^{(2)}, v_{-\vec{k}}^{(2)}(t) | \hat{D}_2(\zeta_{2,\vec{k}}) | u_{\vec{k}}^{(2)}, v_{-\vec{k}}^{(2)}(t) \rangle $$
$$ = \left[ \sum_{u=0}^{\infty} \langle u_{\vec{k}}^{(2)} | \hat{D}_{2,\vec{k}}(\zeta_{2,\vec{k}}) | u_{\vec{k}}^{(2)} \rangle \right] \cdot \left[ \sum_{v=0}^{\infty} \langle v_{-\vec{k}}^{(2)} | \hat{D}_{2,-\vec{k}}(\zeta_{2,-\vec{k}}) | v_{-\vec{k}}^{(2)} \rangle \right]. \quad (D.11) $$

Using the Baker-Campbell-Hausdorff formula, and recalling that $[\hat{a}_{i,\pm\vec{k}}, \hat{a}_{i,\pm\vec{k}}^{\dagger}] = 1$, the one-mode displacement operator (D.4) can be written as

$$ \hat{D}_{i, \pm \vec{k}}(\zeta_{i, \vec{k}}) = e^{-\frac{1}{2} |\gamma_{i, \pm \vec{k}}|^2} \cdot \exp(\gamma_{i, \pm \vec{k}} \hat{a}_{i, \pm \vec{k}}^{\dagger}) \cdot \exp(-\gamma_{i, \pm \vec{k}}^{*} \hat{a}_{i, \pm \vec{k}}), \quad (D.12) $$

where $\gamma_{i, \pm \vec{k}}$ is given in Eq. (D.5). This means that, when evaluating $\langle u_{\pm \vec{k}}^{(i)} | \hat{D}_{i, \pm \vec{k}}(\zeta_{i, \pm \vec{k}}) | u_{\pm \vec{k}}^{(i)} \rangle$, one first has to compute

$$ \exp(-\gamma_{i, \pm \vec{k}}^* \hat{a}_{i, \pm \vec{k}}) |u_{\pm \vec{k}}^{(i)}\rangle = \sum_{n=0}^{u} \frac{1}{n!} (-\gamma_{i, \pm \vec{k}}^*)^n \times \sqrt{u(u-1)\cdots(u-n+1)} |(u-n)_{\pm \vec{k}}^{(i)}\rangle \quad (D.13) $$

where we have simply Taylor expanded the exponential function, and then

$$ \exp(\gamma_{i, \pm \vec{k}} \hat{a}_{i, \pm \vec{k}}^{\dagger}) |(u-n)_{\pm \vec{k}}^{(i)}\rangle = \sum_{m=0}^{\infty} \frac{1}{n!} (\gamma_{i, \pm \vec{k}})^m \times \sqrt{(u-n+1)(u-n+2)\cdots(u-n+m)} \times |(u-n+m)_{\pm \vec{k}}^{(i)}\rangle. \quad (D.14) $$

This leads to

$$ \langle u_{\pm \vec{k}}^{(i)} | \hat{D}_{i, \pm \vec{k}}(\zeta_{i, \pm \vec{k}}) | u_{\pm \vec{k}}^{(i)} \rangle = e^{-\frac{1}{2} |\gamma_{i, \pm \vec{k}}|^2} \sum_{n=0}^{u} \frac{1}{n!} \left(-|\gamma_{i, \pm \vec{k}}|^2\right)^n \times \frac{u(u-1)\cdots(u-n+1)}{n!}. \quad (D.15) $$

The remaining sum over $n$ can be performed by means of the Laguerre polynomials $L_u(z)$, see Eq. (18.5.12) of Ref. [59],

$$ \langle u_{\pm \vec{k}}^{(i)} | \hat{D}_{i, \pm \vec{k}}(\zeta_{i, \pm \vec{k}}) | u_{\pm \vec{k}}^{(i)} \rangle = e^{-\frac{1}{2} |\gamma_{i, \pm \vec{k}}|^2} L_u \left( |\gamma_{i, \pm \vec{k}}|^2 \right). \quad (D.16) $$

Plugging this formula into Eq. (D.11), one obtains

$$ \begin{aligned} & \sum_{u,v=0}^{\infty} \langle u_{\vec{k}}^{(2)}, v_{-\vec{k}}^{(2)}(t) | \hat{D}_2(\zeta_{2,\vec{k}}) | u_{\vec{k}}^{(2)}, v_{-\vec{k}}^{(2)}(t) \rangle \\ &= e^{-\frac{1}{2} (|\gamma_{2,\vec{k}}|^2 + |\gamma_{2,-\vec{k}}|^2)} \sum_{u=0}^{\infty} L_u \left( |\gamma_{2,\vec{k}}|^2 \right) \\ & \quad \times \sum_{v=0}^{\infty} L_v \left( |\gamma_{2,-\vec{k}}|^2 \right). \end{aligned} \quad (D.17) $$

According to Eq. (D.10), we now need to evaluate the Fourier transform of the above expression with respect to $\zeta_{2,\vec{k}}$, in order to compute $N(q_{2,\vec{k}})$. To this end, we introduce the generating function of the Laguerre polynomials [see Eq. (18.12.13) of Ref. [59]],

$$ G_z(t) = \sum_{n=0}^{\infty} t^n L_n(z) = \frac{1}{1-t} \exp \left( -\frac{zt}{1-t} \right). \quad (D.18) $$

Evaluating this formula with $t = 1 - \epsilon$ in the limit $\epsilon \to 0$, one obtains

$$ \sum_{u=0}^{\infty} L_u \left( |\gamma_{2, \pm \vec{k}}|^2 \right) = \lim_{\epsilon_{\pm \vec{k}} \to 0} \left[ G_{|\gamma_{2, \pm \vec{k}}|^2} (1 - \epsilon_{\pm \vec{k}}) \right], \quad (D.19) $$

where we introduce $\epsilon_{\pm \vec{k}}$ for each sector, $\vec{k}$ and $-\vec{k}$. This allows one to write $N(q_{2,\vec{k}})$ as the product of two limits of the Fourier transform of a Gaussian, i.e.

$$ N(q_{2,\vec{k}}) = \prod_{s=+/-} \lim_{\epsilon_{s\vec{k}} \to 0} \left[ \frac{1}{\epsilon_{s\vec{k}}} \int_{\mathbb{R}^2} d^2 \xi_{2,s\vec{k}} \exp \left( i q_{2,s\vec{k}} \cdot \xi_{2,s\vec{k}} \right) \times \exp \left( -\frac{2 - \epsilon_{s\vec{k}}}{4\epsilon_{s\vec{k}}} \xi_{2,s\vec{k}}^2 \right) \right], \quad (D.20) $$

where we have used the fact that $|\gamma_{2,\pm \vec{k}}|^2 = (\xi_{2,\pm \vec{k}}^2 + \kappa_{2,\pm \vec{k}}^2)/2 = \xi_{2,\pm \vec{k}}^2/2$. The Fourier transform of a Gaussian is also a Gaussian and we have

$$ N(q_{2,\vec{k}}) = \prod_{s=+/-} \lim_{\epsilon_{s\vec{k}} \to 0} \left[ \left( \frac{4\pi}{2 - \epsilon_{s\vec{k}}} \right) \exp \left( -\frac{\epsilon_{s\vec{k}}}{2 - \epsilon_{s\vec{k}}} q_{2,s\vec{k}}^2 \right) \right] = (2\pi)^2. \quad (D.21) $$

We conclude that $N(q_{2,\vec{k}})$ simply corresponds to a global phase-space volume, with no dependence on the phase-space location.

The reduced Wigner function is finally given by Eq. (D.9), which leads to

$$ W_{\vec{k}, \text{red}}(\mathbf{q}_1, \vec{k}, t) = \int_{\mathbb{R}^4} d^4 \mathbf{q}_2, \vec{k} \ W_{\vec{k}}(\mathbf{q}_1, \vec{k}, \mathbf{q}_2, \vec{k}, t), \quad (D.22) $$
i.e. it is simply obtained by marginalising the full Wigner function over the environmental degrees of freedom in phase-space. This matches Eq. (6.12) given in the main text. Let us stress that the above result is fully generic and does not assume any specific shape for the Wigner function. It can be easily generalised to an arbitrary number of degrees of freedom both in the system and in the environment sectors, since the Wigner–Weyl transform is generated by a kernel that can be written in a fully factorisable form. Non-factorisability due to entanglement is entirely contained in the Wigner function. Therefore, the reduced Wigner function is always obtained from the full Wigner function by simply integrating over the phase-space variables describing the environmental sector. In a sense, marginalisation in phase space is the Wigner–Weyl representation of the partial trace.

## References

1. S.M. Barnett, P.M. Radmore, Methods in theoretical quantum optics (Clarendon Press Publication, Oxford, 1997)

2. C.M. Caves, B.L. Schumaker, New formalism for two-photon quantum optics. 1. Quadrature phases and squeezed states. Phys. Rev. A **31**, 3068 (1985)

3. B.L. Schumaker, C.M. Caves, New formalism for two-photon quantum optics. 2. Mathematical foundation and compact notation. Phys. Rev. A **31**, 3093 (1985)

4. V.V. Dodonov, ‘Nonclassical’ states in quantum optics: a ‘squeezed’ review of the first 75 years. J. Opt. B **4**, R1 (2002)

5. R. Schnabel, Squeezed states of light and their applications in laser interferometers. Phys. Rep. **684**, 1 (2017). arXiv:1611.03986

6. L. Grishchuk, Y. Sidorov, Squeezed quantum states of relic gravitons and primordial density fluctuations. Phys. Rev. D **42**, 3413 (1990)

7. L. Grishchuk, H. Haus, K. Bergman, Generation of squeezed radiation from vacuum in the cosmos and the laboratory. Phys. Rev. D **46**, 1440 (1992)

8. A. Albrecht, P. Ferreira, M. Joyce, T. Prokopec, Inflation and squeezed quantum states. Phys. Rev. D **50**, 4807 (1994). arXiv:astro-ph/9303001

9. D. Polarski, A.A. Starobinsky, Semiclassicality and decoherence of cosmological perturbations. Class. Quantum Gravity **13**, 377 (1996). arXiv:gr-qc/9504030

10. J. Lesgourgues, D. Polarski, A.A. Starobinsky, Quantum to classical transition of cosmological perturbations for nonvacuum initial states. Nucl. Phys. B **497**, 479 (1997). arXiv:gr-qc/9611019

11. C. Kiefer, D. Polarski, A.A. Starobinsky, Quantum to classical transition for fluctuations in the early universe. Int. J. Mod. Phys. D **7**, 455 (1998). arXiv:gr-qc/9802003

12. J. Martin, V. Vennin, Quantum discord of cosmic inflation: can we show that CMB anisotropies are of quantum-mechanical origin? Phys. Rev. D **93**, 023505 (2016). arXiv:1510.04038

13. J. Grain, V. Vennin, Squeezing formalism and canonical transformations in cosmology. JCAP **2002**, 022 (2020). arXiv:1910.01916

14. PLANCK collaboration, Planck 2018 results. I. Overview and the cosmological legacy of Planck. Astron. Astrophys. **641**, A1 (2020). arXiv:1807.06205

15. N. Turok, String driven inflation. Phys. Rev. Lett. **60**, 549 (1988)

16. T. Damour, A. Vilenkin, String theory and inflation. Phys. Rev. D **53**, 2981 (1996). arXiv:hep-th/9503149

17. S. Kachru, R. Kallosh, A.D. Linde, J.M. Maldacena, L.P. McAllister, S.P. Trivedi, Towards inflation in string theory. JCAP **10**, 013 (2003). arXiv:hep-th/0308055

18. A. Krause, E. Pajer, Chasing brane inflation in string-theory. JCAP **07**, 023 (2008). arXiv:0705.4682

19. D. Baumann, L. McAllister, *Inflation and string theory, Cambridge monographs on mathematical physics* (Cambridge University Press, Cambridge, 2015). https://doi.org/10.1017/CBO9781316105733. arXiv:1404.2601

20. W.H. Zurek, Pointer basis of quantum apparatus: into what mixture does the wave packet collapse? Phys. Rev. D **24**, 1516 (1981)

21. W.H. Zurek, Environment-induced superselection rules. Phys. Rev. D **26**, 1862 (1982)

22. E. Joos, H. Zeh, The Emergence of classical properties through interaction with the environment. Z. Phys. B **59**, 223 (1985)

23. J. Martin, V. Vennin, Observational constraints on quantum decoherence during inflation. JCAP **05**, 063 (2018). arXiv:1801.09949

24. E. Oudot, P. Sekatski, F. Fröwis, N. Gisin, N. Sangouard, Two-mode squeezed states as Schrödinger cat-like states. J Opt Soc Am B Opt Phys **32**, 2190 (2015). arXiv:1410.8421

25. J. Martin, V. Vennin, Leggett–Garg inequalities for squeezed states. Phys. Rev. A **94**, 052135 (2016). arXiv:1611.01785

26. S. Choudhury, S. Panda, Quantum entanglement in de Sitter space from stringy axion: an analysis using $\alpha$ vacua. Nucl. Phys. B **943**, 114606 (2019). arXiv:1712.08299

27. K. Ando, V. Vennin, Bipartite temporal Bell inequalities for two-mode squeezed states. Phys. Rev. A **102**, 052213 (2020). arXiv:2007.00458

28. S. Kanno, J. Soda, J. Tokuda, Indirect detection of gravitons through quantum entanglement. arXiv:2103.17053

29. Y.S. Kim, M.E. Noz, Dirac matrices and Feynman's rest of the universe. Symmetry **4**, 626 (2012). arXiv:1210.6251

30. K. Hasebe, $Sp(4; \mathbb{R})$ squeezing for Bloch four-hyperboloid via the non-compact Hopf map. J. Phys. A **53**, 055303 (2020). arXiv:1904.12259

31. A. Garcia-Chung, Symplectic group in polymer quantum mechanics. Phys. Rev. D **101**, 106004 (2020). arXiv:2003.00388

32. A. Garcia-Chung, Squeeze operator: a classical view. arXiv:2003.04257

33. G. Chacón-Acosta, A. García-Chung, The relation between the symplectic group $Sp(4, \mathbb{R})$ and its Lie algebra: its application in polymer quantum mechanics. arXiv:2102.12049

34. H. Goldstein, C. Poole, J. Safko, *Classical mechanics* (Addison Wesley, Boston, 2002)

35. C. Bloch, A. Messiah, The canonical form of an antisymmetric tensor and its application to the theory of superconductivity. Nucl. Phys. **39**, 95 (1962)

36. R. Puri, T. Asakura, K. Brenner, T. Hansch, F. Krausz, H. Weber et al., *Mathematical methods of quantum optics, Physics and astronomy online library* (Springer, Berlin, 2001)

37. A. Perelomov, *Generalized coherent states and their applications, modern methods of plant analysis* (Springer, Berlin, 1986)

38. S. Barnett, P. Radmore, Methods in theoretical quantum optics, Oxford series in optical and imaging sciences (Clarendon Press, Oxford, 2002)

39. D. Truax, Baker–Campbell–Hausdorff relations and unitarity of SU(2) and SU(1,1) squeeze operators. Phys. Rev. D **31**, 1988 (1985)

40. R. Simon, E.C.G. Sudarshan, N. Mukunda, Gaussian–Wigner distributions in quantum mechanics and optics. Phys. Rev. A **36**, 3868 (1987)

41. R. Simon, E. Sudarshan, N. Mukunda, Gaussian Wigner distributions: a complete characterization. Phys. Lett. A **124**, 223 (1987)

42. B. Arvind, N Mukunda Dutta, R. Simon, The real symplectic groups in quantum mechanics and optics. Pramana **45**, 471 (1995). arXiv:quant-ph/9509002
43. J. Martin, Inflationary perturbations: the cosmological Schwinger effect. Lect. Notes Phys. **738**, 193 (2008). [arXiv:0704.3540](https://arxiv.org/abs/0704.3540)

44. A. Matacz, The emergence of classical behavior in the quantum fluctuations of a scalar field in an expanding universe. Class. Quantum Gravity **10**, 509 (1993)

45. M.B. Einhorn, F. Larsen, Squeezed states in the de Sitter vacuum. Phys. Rev. D **68**, 064002 (2003). [arXiv:hep-th/0305056](https://arxiv.org/abs/0305056)

46. W.B. Case, Wigner functions and Weyl transforms for pedestrians. Am. J. Phys. **76**, 937 (2008). [https://doi.org/10.1119/1.2957889](https://doi.org/10.1119/1.2957889)

47. H. Weyl, Quantenmechanik und Gruppentheorie. Z. Angew. Phys. **46**, 1 (1927)

48. H.J. Groenewold, On the principles of elementary quantum mechanics. Physica **12**, 405 (1946)

49. J.E. Moyal, M.S. Bartlett, Quantum mechanics as a statistical theory. Proc. Camb. Philos. Soc. **45**, 99 (1949)

50. W.B. Case, Wigner functions and Weyl transforms for pedestrians. Am. J. Phys. **76**, 937 (2008)

51. U. Seyfarth, A.B. Klimov, H. de Guise, G. Leuchs, L.L. Sanchez-Soto, Wigner function for SU(1,1). Quantum **4**, 317 (2020). [arXiv:1911.11703](https://arxiv.org/abs/1911.11703)

52. A. Serafini, F. Illuminati, S. De Siena, Von Neumann entropy, mutual information and total correlations of Gaussian states. J. Phys. B **37**, L21 (2004). [arXiv:quant-ph/0307073](https://arxiv.org/abs/quant-ph/0307073)

53. T.F. Demarie, Pedagogical introduction to the entropy of entanglement for Gaussian states (2012). [arXiv:1209.2748](https://arXiv.org/abs/1209.2748)

54. J. Eisert, M. Cramer, M.B. Plenio, Area laws for the entanglement entropy—a review. Rev. Mod. Phys. **82**, 277 (2010). [arXiv:0808.3773](https://arXiv.org/abs/0808.3773)

55. J. Martin, V. Vennin, Non Gaussianities from quantum decoherence during inflation. JCAP **06**, 037 (2018). [arXiv:1805.05609](https://arXiv.org/abs/1805.05609)

56. L. Pinol, Multifield inflation beyond $N_{field} = 2$: non-Gaussianities and single-field effective theory. [arXiv:2011.05930](https://arXiv.org/abs/2011.05930)

57. J. Martin, V. Vennin, Bell inequalities for continuous-variable systems in generic squeezed states. Phys. Rev. A **93**, 062117 (2016). [arXiv:1605.02944](https://arXiv.org/abs/1605.02944)

58. J. Martin, V. Vennin, Obstructions to bell CMB experiments. Phys. Rev. D **96**, 063501 (2017). [arXiv:1706.05001](https://arXiv.org/abs/1706.05001)

59. NIST Digital Library of Mathematical Functions, Release 1.1.1 of 2021-03-15. [http://dlmf.nist.gov/](http://dlmf.nist.gov/)