A gravitationally induced decoherence model using Ashtekar variables

Decoherence models play an important role in the semiclassical sector of quantum models in various areas of physics from condensed matter, quantum optics and quantum information [1–5] to cosmology [6–10], for a mathematical perspective see [11]. These models are formulated in the framework of open quantum systems [12–14] and allow for instance to address the question of the quantum-to-classical transition [15–17]. Because such models involve an interaction between the quantum system and its environment the choice of the latter is a crucial ingredient for any decoherence model. A very common choice is to couple the system to a bath of harmonic oscillators [12, 18, 19]. In this work we want to analyse gravitationally induced decoherence for which the environment is taken to be gravity. This provides a platform to formulate and analyse decoherence models based on existing models for quantum gravity. Compared to the situation in closed quantum models where the semiclassical sector is usually described by means of semiclassical states, see for instance [20–23], the semiclassical sector in the context of quantum gravity can be investigated from a different angle within decoherence models with a gravitational environment.

There already exist several models for gravitationally induced decoherence in the literature [24–38], see also the reviews in [39, 40] and the review [41] for further models in which quantum matter is coupled to gravity. Some of them have for instance been used to study the influence of decoherence on standard model matter such as neutrinos and their oscillations, see for instance [27–34]. All models have in common that the system's dynamics is encoded in a so-called master equation that is an effective first order differential equation for the density matrix of the system, where this density matrix is obtained once the environmental degrees of freedom in the total dynamics have been traced out. One can group the existing models in two classes, that we will denote as bottom-up [25, 26, 28, 30, 33] and top-down approaches [35–38, 42–44], respectively. In the first approach one takes as a starting point a given master equation, as for instance the well known Lindblad equation [45, 46], and then constrains existing free parameters phenomenologically for a given model under consideration. In the top-down approach one starts with a given classical model, performs a quantisation of the total system and afterwards derives the master equation within that quantum model. If one aims at formulating a decoherence model based on a specific quantum gravity approach then the top-down approach is advantageous, since an interesting question in this context is whether there exist characteristic features in the decoherence model that will distinguish one quantum gravity model from another. Furthermore, in the derivation of the Lindblad equation, generally several assumptions such as the first and second Markov approximation are applied and one needs to check carefully under which circumstances these can be applied in gravitationally induced decoherence models. Compared to decoherence models derived in a quantum mechanical context, choosing quantum gravity as the environment complicates the situation in two aspects. On the one hand we have to deal with a quantum field theory where also the coupling between the system and the environment is more complex and on the other hand we consider a gauge theory and both facts need to be addressed appropriately in these models.

In this article we are interested in deriving a gravitationally induced decoherence model inspired by loop quantum gravity, see [47–49] for textbooks. Most works on the semiclassical physics of loop quantum gravity rather focus on closed quantum systems with semiclassical states [50–53] or coarse graining in spin foam models [54, 55]. There exist a few exceptions see for instance [56] where a decoherence toy model for surface states in the context of quantum black holes is introduced or the work in [57] for a coupling of a spin foam model to an environment using the adiabatic elimination method [58] to obtain an effective Lindblad equation.

As a first step in this direction we consider a top-down approach in which we couple a scalar field to linearised gravity as it has been done in [35–37]. In contrast to [35–37] where ADM variables have been used we construct the model using Ashtekar variables [59, 60]. The latter allow in principle a quantisation using techniques from loop quantum gravity along the lines of [61, 62]. However, as a first step due to the complexity of the model and to be able to better compare our results to the existing literature, we apply a Fock quantisation to the linearised gravitational sector.

Already at the classical level the model considered in this work requires further generalisations to the already existing results. The models in [35–37] perturb the gravitational degrees of freedom only and we consider the coupling of the scalar field in the framework of a Post-Minkowski approximation scheme that allows to consider a setup where only the gravitational degrees of freedom are perturbed and still obtain a consistent formulation of the involved constraints in perturbation theory. Because we work with Ashtekar variables the model involves an additional Gauß constraint that needs to be taken into account for the gauge invariance. The route we follow here to deal with the gauge freedom is that we consider a reduced phase space quantisation. For this purpose we first construct a canonical transformation that separates the gauge and physical degrees of freedom into two sets. The latter are obtained by constructing Dirac observables in the relational formalism applying techniques from [63–68] where suitable reference fields need to be chosen. In our case, these will be so-called geometrical clocks that only depend on the gravitational degrees of freedom. These clocks differ from the ones chosen in [69], where complex Ashtekar variables have been used with the corresponding gauge fixing, whereas the Dirac observables used here have a natural relation to the gauge-fixed quantities different than the one presented in [61], which considers the vacuum case only. Geometrical clocks have also been applied in different contexts, see for instance [69–72]. Working in the relational formalism (physical) time can be defined in a relational way and circumvents some of the issues related to fluctuating light cones mentioned in [73]. In order to construct the above mentioned canonical transformation we need to generalise the observable map and introduce a kind of dual observable map where the role of constraints and reference fields is interchanged. When applied to the constraints themselves it reproduces the results for the perturbative formulation of weak Abelianisation in [69] but being more general in the sense that it can be applied to any phase space function, a property crucial for our application. As we work in the linearised theory when computing the algebra of Dirac observables that provides the starting point for the reduced phase space quantisation, it is necessary to consider Dirac observables up to second order in general to formulate the model at the gauge invariant level. In former works [36, 38] often specific gauge fixings were chosen or Dirac observables were constructed for the matter sector only [35], where to our understanding these Dirac observables are not left invariant under the spatial diffeomorphism constraint. In [37] Dirac quantisation was applied and the obtained results are consistent with the reduced phase space quantisation procedure we apply here. We also consider a second choice of geometrical clocks introduced in [74] which was carried over to Ashtekar variables in [69] and discuss its differences and similarities to the geometrical clocks chosen in this work. There exist also former work on decoherence effects caused by quantum clocks [75–77], whereas we use classical (perfect) clocks in the language of [77]. Such decoherence effects would be present in addition to the decoherence effects that we will discuss in this work and it will be interesting to compare the size of the two types of decoherence effects. A detailed analysis, requires to work with fundamental quantum clocks and perform the construction of Dirac observables in the quantum theory for instance along the lines of [78]. However, to the knowledge of the authors this has so far only been considered in the context of quantum mechanical systems but not a full quantum field theory setup that we will need in the case considered in this work.

The physical Hamiltonian that is itself a Dirac observable and generates the dynamics of the linearised Dirac observables of the total system is then quantised using a Fock quantisation. The strategy we follow in the quantisation is that we choose to normal order the entire physical Hamiltonian operator. This differs from [35, 37, 38] where either individual operators in products of operators are normal ordered separately or no normal ordering is applied at all and such differences appear in the operator ordering in the self-interaction part of the Hamiltonian operator. For deriving a master equation in a top-down approach two prominent techniques exist, the influence functional approach [79] and the projection operator technique [80–83], where we use the latter to obtain a time convolution-less (TCL) master equation for the model considered in this work. This allows us to understand in detail how certain assumptions and approximations enter into the derivation of the master equations which will be useful also for later applications and generalisations of the model. We choose a Gibbs state as the initial state of the quantum gravitational environment and use a regularisation by means of a box of finite volume in order to implement this on Fock space, similar to the case in [84]. After computing the environmental trace we remove the regulator by taking the infinite volume limit. In order to circumvent to work with regulators one can use KMS states [85, 86], which is beyond the scope of this work. We express the final master equation in terms of thermal Wightman functions, which on the one hand is a general formalism other models can also be formulated in and reflects the structure of the involved correlation functions. On the other hand this gives an easy access to the special case where the environment is chosen to be in the vacuum state. For a better comparison with the already existing models in [35, 37, 38] we present three equivalent but different forms of the final master equation and discuss in detail the similarities and differences to [35, 37, 38].

The paper is structured as follows: After the introduction in section 1 in section 2 we present the classical formulation of the decoherence model considered in this article. Section 2.2 briefly reviews the existing results on linearised gravity in terms of Ashtekar variables and we consider the coupling of a scalar field to linearised gravity in a Post-Minkowski approximation scheme in section 2.3. Section 2.4 presents the construction of Dirac observables for the linearised model using geometrical clocks. Furthermore, we introduce a canonical transformation on the full phase space in which the geometrical clocks as well as the constraints are chosen as some of the new canonical coordinates, allowing to clearly separate physical and gauge degrees of freedom in the model. This decomposition, often also called Kuchař decomposition [87, 88] that was recently considered in [89] to discuss different choices of gauge in the context of linear cosmological perturbation theory. In the generic situation where a chosen set of set of clocks is not necessary abelian a priori, as it was for instance be assumed in [69] and no comments on the case of non-abelian gauge fixing conditions can be found in [89], this requires the application of the already existing observable map as well as its dual version introduced in this section. This opens the possibility for constructing such canonical transformations and hence performing a Kuchar decomposition in principle also in other models and beyond linear order. For each choice of geometrical clocks there exists a corresponding gauge in which the gauge invariant variables have a natural and simple interpretation. The advantage of working at the gauge invariant level is that even if we choose a different gauge we still have a gauge invariant formulation of the model, with an in general different interpretation of the gauge invariant variables but the equations are valid in any chosen gauge. If one is interested in always choosing the set of geometrical clocks for which the Dirac observables have a natural interpretation, then this is also possible in this framework. Since in the work presented here the choice of geometrical clocks is implemented by a canonical transformation, it can be inverted and thus afterwards one can choose a different gauge and hence also switch between different choices of geometrical clocks, which can again be carried over to choosing different gauges. We have discussed to examples for different set of geometrical clocks in appendix A.3.4.

Given the reduced phase space derived in section 2.4 and the physical Hamiltonian, we briefly discuss the Fock quantisation of the model in section 3 with a focus on the individual contributions in the physical Hamiltonian operator. The derivation of the final master equation is presented in section 4. We include a brief review of the projection operator technique as well as the influence functional approach in sections 4.1 and 4.2, respectively. This serves as a preparation for the application of the first technique to the model in this work, presented in section 4.3 where we obtain the explicit form of the final master equation. In order to be able to compare our results better with the already existing ones in the literature we present the final master equation in three equivalent forms. The first version can be found in (4.52) and is written in terms of a dissipation and noise kernel. Next in (4.58) an equivalent form of the final master equation is presented which is very close to the notation used in [35, 37, 38] and allows a direct comparison to their work. Furthermore, the third equivalent version of the master equation, found in (4.60), clearly separates the contributions from the zero and finite temperature cases. These equivalent version of the final master equation are further compared to already existing master equations for similar models in the literature. The main differences to the existing results are: First the model presented here is formulated in terms of Dirac observables with respect to chosen geometrical clocks or also known as dynamical reference systems. The formalism allows different choices of such sets of geometrical clocks and thus also an analysis how such choices affect a given model in future work. Furthermore, the final master equation we obtain is not automatically of Lindblad type. Often approximations such as the second Markov approximation and the rotating wave approximation are applied to arrive at a Lindbladian master equation. At least to the authors it is not clear from the outset that such techniques, which are mostly used in the quantum mechanical context, can be directly transferred to the context of quantum field theories and, according to the authors, require further careful and detailed analysis, as also shown, for example, by the work in [90]. The resulting final master equation in this work provides an good starting point for such an investigation. Finally, in section 5 we summarise our results and conclude. We present more detailed computations in the appendix that we refer to in the main text of this article.

As the starting point for deriving the decoherence model in this work we choose the following action involving gravity coupled to a Klein–Gordon scalar field in the ADM decomposition [91] expressed in terms of Ashtekar variables [59, 60, 92, 93]:

$$\begin{align}\begin{array}{l} S = \mathop \smallint \limits_ dt\mathop \smallint \limits_\sigma {d^3}x\left( {\frac{1}{{\kappa \beta }}{{\dot A}_a}^i(\vec x,t){E^a}_i(\vec x,t) + \dot \phi (\vec x,t)\pi (\vec x,t)} \right.\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left. { - \left[ {{{\mathrm{\Lambda }}^i}(\vec x,t){G_i}(\vec x,t) + {N^a}(\vec x,t){C_a}(\vec x,t) + N(\vec x,t)C(\vec x,t)} \right]} \right), \end{array}\end{align} \tag{ 2.1 }$$

where β is the Barbero-Immirzi parameter, $\kappa = \frac{8\pi G_N}{c^4}$ with Newton's gravitational constant G_N , c denotes the speed of light and we work with the mostly-plus signature of the metric. From here on, we will set $c = \hbar = 1$. The gravitational degrees of freedom are encoded in the Ashtekar variables consisting of an SU(2)-connection ${A}{_a{}^i}(\vec{x},t)$ and the canonically conjugate densitised triads ${E}{^a{}_i}(\vec{x},t)$. The matter sector includes the scalar field $\phi(\vec{x},t)$ as well as its canonically conjugate momentum $\pi(\vec{x},t)$, an overdot denotes a derivative with respect to the temporal coordinate. That general relativity is a fully constrained theory is reflected by the fact that apart from the symplectic potential in the action, there appears a sum of constraints $G_i(\vec{x},t)$, $C_a(\vec{x},t)$ and $C(\vec{x},t)$ only multiplied by different Lagrange multipliers $\Lambda^i(\vec{x},t)$, $N^a(\vec{x},t)$ and $N(\vec{x},t)$, where the latter two are usually referred to as the shift vector and the lapse function, respectively. In Ashtekar variables the spatial diffeomorphism constraint reads

$$\begin{align} {C}{_a}(\vec{x},t) = \frac{1}{\kappa\beta}{F}{_{ab}{}^i}(\vec{x},t)\, {E}{^b{}_i}(\vec{x},t) + \pi(\vec{x},t)\, \partial_a \phi(\vec{x},t)\,, \end{align} \tag{ 2.2 }$$

where ${F}{_{ab}{}^i}(\vec{x},t) : = \partial_{a} {A}{_b{}^i}(\vec{x},t)-\partial_{b} {A}{_a{}^i}(\vec{x},t) + {\epsilon{_{ijk}}} {A}{_a{}^j}(\vec{x},t)\, {A}{_b{}^k}(\vec{x},t)$ denotes the curvature related to the connection variable A. The Gauß constraint, which is absent if working in ADM variables, is given by

$$\begin{align} {G{{}_i}}(\vec{x},t) = \frac{1}{2\kappa\beta} \left({{\partial}{_a}} {E}{^a{}_i}(\vec{x},t) + {{\epsilon}{_{ij}{}^k}} {A}{_a{}^j}(\vec{x},t)\, {{E}{^a{}_k}}(\vec{x},t)\right) \end{align} \tag{ 2.3 }$$

with the completely antisymmetric tensor ε. Finally, the Hamiltonian constraint takes the form

$$\begin{align} C(\vec{x},t) = &\frac{1}{2\kappa}\left( {F}{_{ab}{}^i}(\vec{x},t) - \frac{\beta^2 +1}{\beta^2} {{\epsilon}{_{ilm}}} ({{A}{_a{}^l}}(\vec{x},t)-{{\Gamma}{_a{}^l}}(\vec{x},t)) ({A}{_b{}^m}(\vec{x},t)-{{\Gamma}{_b{}^m}}(\vec{x},t)) \right) \nonumber\\ & \times \frac{{{\epsilon}^{ijk}} {{E}{^a{}_j}}(\vec{x},t)\, {E}{^b{}_k}(\vec{x},t)}{\sqrt{\det({E}{^c{}_n})}} +\frac{\pi^2(\vec{x},t)}{2\sqrt{\det({E}{^c{}_n})}} + \frac{1}{2 \sqrt{\det({E}{^c{}_n})}} \nonumber\\ & \times \delta^{ij}{E}{^a{}_i}(\vec{x},t)\, {E}{^b{}_j}(\vec{x},t)\, ({{\partial}{{}_a}} \phi(\vec{x},t))({{\partial}{_b}} \phi(\vec{x},t)) + \frac{\sqrt{\det({E}{^c{}_n})}}{2}\; m^2\phi^2(\vec{x},t)\,. \end{align} \tag{ 2.4 }$$

Here, m denotes the mass of the scalar field and ${\Gamma}{_a{}^i}$ is the spin connection, considered a functional of the densitised triads, i.e.

$$\begin{align} {\Gamma}{_a{}^i}= -\frac{1}{2} {{\epsilon}^{ijk}} {E}{^b{}_k} \left({{E}{_a{}^j{}_{,b}}} - {{E}{_b{}^j{}_{,a}}} + {E}{^c{}_j} {E}{_a{}^l} {{E}{_c{}^l{}_{,b}}} + {{E}{_a{}^j}} \frac{\det{({E}{^c{}_l})_{,b}}}{\det{({E}{^c{}_l})}} - {{E}{_b{}^j}} \frac{\det{({E}{^c{}_l})_{,a}}}{2\det{({E}{^c{}_l})}} \right)\,, \end{align} \tag{ 2.5 }$$

where we abbreviated partial derivatives with a comma, that is ${{E}{_a{}^j{}_{,b}}} : = {{\partial}{_b}} {{E}{_a{}^j}}$ and suppressed the spatial and temporal arguments of the involved field variables. The elementary phase space variables satisfy the following Poisson algebra on the classical phase space:

$$\begin{align} \{{E}{^a{}_i}(\vec{x},t),{E}{^b{}_j}(\vec{y},t)\} = \{{A}{_a{}^i}(\vec{x},t),\,{A}{_b{}^j}(\vec{y},t) \} = 0\,, \quad \{{A}{_a{}^i}(\vec{x},t),\, {E}{^b{}_j} (\vec{y},t)\} = \beta\kappa \delta_a^b\delta^i_j \delta^3(\vec{x}-\vec{y})\nonumber\\ \end{align} \tag{ 2.6 }$$

$$\begin{align} \{\phi(\vec{x},t), \phi(\vec{y},t) \} = \{\pi(\vec{x},t), \pi(\vec{y},t) \} = 0\,, \qquad \{\phi(\vec{x},t), \pi(\vec{y},t) \} = \delta^{(3)}(\vec{x}-\vec{y})\,, \end{align} \tag{ 2.7 }$$

where all the remaining ones vanish. As the decoherence model presented in this article requires to linearise the above action around a flat Minkowski background in the course of the upcoming section, we choose asymptotically flat boundary conditions [94] also addressed in [95, 96] for the case of real Ashtekar variables. With the fall-off behaviour of the Ashtekar variables discussed in section 2.2 and the linearisation of the Lagrange multipliers chosen such that both are consistent with an asymptotically flat Universe, it turns out that the Gauß constraint and the gravitational part of the spatial diffeomorphism constraint do not cause any issues and need no further modification in accordance with [69]. However, the gravitational part of the Hamiltonian constraint as well as its variation is not well-defined in the asymptotically flat limit and contains divergences. To circumvent these, we introduce a suitable boundary term given by

$$\begin{align} T[N] = -2 \int_{\partial \sigma} {{dS}{{}_a}} N {A}{_b{}^k} \epsilon^{ijk} \frac{{E}{^a{}_i} {E}{^b{}_j} }{\sqrt{\det E}} \end{align} \tag{ 2.8 }$$

to the smeared Hamiltonian constraint, which makes it well-defined and functionally differentiable. Note that on the constraint hypersurface where the Gauß constraint vanishes $T[1]$ is the so-called ADM-energy [97]. Further we assume, as is usually done, that the initial data of the matter variables have compact support, yielding matter contributions that are well-defined in the asymptotically flat case. Given this, the canonical Hamiltonian that we will work with in the following and that is consistent with our chosen boundary conditions reads

$$\begin{align} \mathbf{H}_\textrm{can} &= C(N)+\vec{C}(\vec{N})+\vec{G}(\vec{\Lambda})+T[N] \nonumber\\ &= \int_\sigma d^3x \bigg[ \frac{1}{\kappa\beta} N^a {F}{_{ab}{}^j} {E}{^b{}_j} + N^a \pi {{\nabla}{{}_a}} \phi+\frac{1}{2\kappa\beta} \Lambda^i{{\partial}{_a}} {E}{^a{}_i} + \frac{1}{2\kappa\beta} \Lambda^i {{\epsilon}_{ikl}} {{A}{_a{}^k}} {{E}{^a{}_l}}\nonumber \\ &\quad - N\epsilon^{jkl} \frac{1}{\sqrt{\det E}} {{E}{^a{}_k}} {{E}{^b{}_l}}\frac{1}{2\kappa\beta^2} {{\epsilon}{_{jmn}}} \left({{A}{_a{}^m}} {{A}{_b{}^n}} + (\beta^2+1) {{\Gamma}{_a{}^m}} {{\Gamma}{_b{}^n}} - 2 {{\Gamma}{_a{}^m}}{{A}{_b{}^n}}\right) \nonumber\\ &\quad + \frac{1}{\sqrt{\det E}} \bigg(-\frac{1}{\kappa} \epsilon^{jkl} {{E}{^a{}_k}} {{E}{^b{}_l}} {A}{_b{}^j}{{\partial}{_a}} N + N \frac{\pi^2}{2} + \frac{N}{2} {E}{^a{}_i} {E}{^b{}_i} ({{\partial}{{}_a}} \phi)({{\partial}{_b}} \phi)\bigg) \nonumber\\ &\quad + N \frac{\sqrt{\det E}}{2} \;m^2\phi^2\bigg]\,, \end{align} \tag{ 2.9 }$$

where ${C}({N}) : = \int_\sigma d^3x\; C(\vec x,t) N(\vec x,t)$, $\vec{C}(\vec{N}) : = \int_\sigma d^3x\; C_a(\vec x,t) N^a(\vec x,t)$ and $\vec{G}(\vec{\Lambda}) : = \int_\sigma d^3x$ $G_i(\vec x,t) \Lambda^i(\vec x,t)$.

In the following two subsections we will consider a linearisation of the gravity-matter system. For this purpose we first briefly review the linearisation of the vacuum case and afterwards discuss how the scalar field can be coupled to linearised gravity in the framework of a Post-Minkowski approximation scheme.

A quantisation of the full theory of general relativity is a highly complex task. As outlined in the introduction, one of the motivations for the work presented in this article comes from scenarios that involve weak couplings between the matter system and the gravitational field. Therefore, we will consider a linearisation of general relativity around a Minkowski spacetime and then consider a perturbation for the gravitational degrees of freedom using κ as an expansion parameter. Here we will focus on vacuum gravity using Ashtekar variables first, analogous to the work in for instance [61] or [98]. We denote the linear perturbations with a prefix δ. The Ashtekar variables as well as the Lagrange multipliers then become

$$\begin{align} {E}{^a{}_i} &= \overline{E}^{\,a}_i + \kappa\, \delta{E}{^a{}_i}=\delta^a_i + \kappa\, \delta{E}{^a{}_i}\qquad \delta E : O(r^{-1}) \text{even} \end{align} \tag{ 2.10 }$$

$$\begin{align} {A}{_a{}^i} &=\overline{A}^i_a + \kappa\, \delta{A}{_a{}^i}= 0 + \kappa\, \delta{A}{_a{}^i}\qquad\;\;\; \delta A : O(r^{-2}) \text{odd} \end{align} \tag{ 2.11 }$$

$$\begin{align} N &=\overline{N} + \kappa\, \delta N=1 + \kappa\, \delta N \qquad\qquad\;\;\;\; \delta N : O(r^{-1}) \end{align} \tag{ 2.12 }$$

$$\begin{align} N^a & = \overline{N}^{\,a} + \kappa\, \delta N^a=0 + \kappa\, \delta N^a \qquad\quad\, \delta N^a : O(r^{-1}) \end{align} \tag{ 2.13 }$$

$$\begin{align} \Lambda^i &= \overline{\Lambda}^i + \kappa\, \delta\Lambda^i=0 + \kappa\, \delta\Lambda^i\qquad\quad\;\; \delta \Lambda : O(r^{-2}) \text{even}\,, \end{align} \tag{ 2.14 }$$

where the fall-off behaviour in the limit of spatial infinity of the first order perturbations has been chosen in accordance with the asymptotically flat boundary condition [94]. The split into background and perturbed quantities is chosen such that the background corresponds to a Minkowski spacetime that remains unchanged by the action of the canonical Hamiltonian of the system, i.e.

$$\begin{align} \overline{\{{E}{^a{}_i},\mathbf{H}_\textrm{can}\}} = 0 \end{align} \tag{ 2.15 }$$

$$\begin{align} \overline{\{{A}{_a{}^i},\mathbf{H}_\textrm{can}\}} = 0\,, \end{align} \tag{ 2.16 }$$

where the overline stands for evaluation with respect to the background Minkowski spacetime. In the case of vacuum gravity (after computation of the Poisson brackets) this amounts to $\overline{E}^{\,a}_i = \delta^a_i$, $\overline{A}_a^i = 0$, $\overline{N} = 1$ and $\overline{N}^{\,a} = 0$. As one can easily compute, the spin connection vanishes in the background:

$$\begin{align} {\Gamma}{_a{}^i} &= \overline{\Gamma}^i_a +\kappa\, \delta{\Gamma}{_a{}^i} \approx 0+\kappa\, \delta{\Gamma}{_a{}^i} \nonumber\\ &=- \frac{\kappa\, }{2} \epsilon^{ijk} \delta_k^b \left[- \delta_a^l \delta_c^j {{\partial}{_b}} (\delta{E}{^c{}_l}) + \delta_c^j \delta_b^l{{\partial}{_a}} (\delta{E}{^c{}_l}) - \delta_{ac}{{\partial}{_b}} (\delta{E}{^c{}_j}) + \delta_a^j \delta_c^l {{\partial}{_b}} (\delta{E}{^c{}_l}) \right]\,. \end{align} \tag{ 2.17 }$$

The Poisson algebra of the linearised gravitational variables can be inherited from (2.6):

$$\begin{align} & \{\delta {E}{^a{}_i}(\vec{x},t),\delta {E}{^b{}_j}(\vec{y},t)\} = \{\delta {A}{_a{}^i}(\vec{x},t),\delta {A}{_b{}^j}(\vec{y},t) \} = 0\,,\nonumber\\ & \{\delta {A}{_a{}^i}(\vec{x},t), \delta {E}{^b{}_j} (\vec{y},t)\} = \frac{\beta}{\kappa} \delta_a^b\delta^i_j \delta^3(\vec{x}-\vec{y})\,. \end{align} \tag{ 2.18 }$$

In the next subsection we will address how the matter contributions can be included consistently in the linearised framework.

In the context of general relativistic perturbation theory, one usually chooses a background solution and then considers perturbations of the gravitational and matter sector around it assuming that the perturbations are small compared to the chosen background quantities. In this work we linearise around a flat Minkowski spacetime, which is a vacuum solution of general relativity. Hence, any considered perturbation in the matter sector will not be small compared to vanishing matter degrees of freedom in the vacuum case. However, we can formulate a model that involves the coupling between matter and linearised gravity in a Post-Minkowski approximation scheme. For the model presented in this article we will need to consider the zeroth and first order in the Post-Minkowski formalism only and this will be the underlying guiding principle of how the matter sector will be included into the Hamiltonian formulation and in particular into the constraints. We will briefly sketch the main steps that lead to the linearised action that we take as our starting point. The Post-Minkowski formalism is based on the Landau–Lifshitz formulation of general relativity, see for instance the book [99] for an introduction to the subject, and one starts by introducing the so-called gothic metric defined as

$$\begin{align} \mathfrak{g}^{\mu\nu} := \sqrt{-\det(g_{\rho\sigma})}\, g^{\mu\nu}, \end{align} \tag{ 2.19 }$$

by means of which we can introduce the following tensor density

$$\begin{align} H^{\mu\nu\rho\sigma} := \mathfrak{g}^{\mu\rho}\mathfrak{g}^{\nu\sigma}-\mathfrak{g}^{\mu\sigma}\mathfrak{g}^{\nu\rho}\,, \end{align} \tag{ 2.20 }$$

where $H^{\mu\nu\rho\sigma}$ carries the same symmetries as the Riemann tensor. Following the presentation in [99] one can use (2.20) to rewrite Einstein's equations $G_{\mu\nu} = \kappa T_{\mu\nu}$ in the following form:

$$\begin{align} \partial_\nu\partial_\sigma H^{\mu\nu\rho\sigma} = 2\kappa (-\det(g)) \left[ T^{\mu\rho} + t^{\mu\rho}_{LL}\right]\,. \end{align} \tag{ 2.21 }$$

Here, $T^{\mu\rho}$ denotes the energy-momentum tensor of the matter degrees of freedom and $t^{\mu\rho}_{LL}$ is the so-called Landau–Lifshitz pseudotensor which consists of a sum of (contractions of) terms $\partial \mathfrak{g}\, \partial \mathfrak{g}$ and does not transform as a tensor under coordinate transformations. This quantity corresponds to the distribution of energy of the gravitational field in spacetime. If one aims at formulating the corresponding linearised theory, as a first step, one applies a partial gauge fixing on the gothic metric by the harmonic coordinate condition

$$\begin{align} \partial_\mu \mathfrak{g}^{\mu\nu}=0\,. \end{align} \tag{ 2.22 }$$

For later convenience it is useful to introduce the following quantity also known as the metric potentials $\mathfrak{h}^{\mu\nu}$ that have the following relation to the gothic metric $\mathfrak{g}^{\mu\nu}$:

$$\begin{align} \mathfrak{g}^{\mu\nu}=\eta^{\mu\nu}-\mathfrak{h}^{\mu\nu}\,. \end{align} \tag{ 2.23 }$$

The harmonic coordinate condition then carries over to $\partial_\mu \mathfrak{h}^{\mu\nu} = 0$. In this gauge, Einstein's equations in (2.21) have the form, see e.g. equations (6.51) and (6.52) in [99]:

$$\begin{align} \Box \mathfrak{h}^{\mu\nu} = -2\kappa (-\det(g)) \left[ T^{\mu\nu}[\Phi,g] + t^{\mu\nu}_{LL}[\mathfrak{h}] + t_H^{\mu\nu}[\mathfrak{h}]\right], \end{align} \tag{ 2.24 }$$

with the d'Alembertian $\Box$ in flat spacetime and

$$\begin{align} t^{\mu\nu}_H[\mathfrak{h}] := \frac{1}{2\kappa (-\det(g))} \left( \partial_\rho \mathfrak{h}^{\mu\sigma} \; \partial_\sigma \mathfrak{h}^{\nu\rho} - \mathfrak{h}^{\rho\sigma} \; \partial_{\rho} \partial_\sigma \mathfrak{h}^{\mu\nu}\right) \,. \end{align} \tag{ 2.25 }$$

Among the individual contributions in (2.24) the energy-momentum tensor $T^{\mu\nu}$ depends on the matter variables here denoted by Φ and on the metric g, while the other two contributions depend on the modified gothic metric $\mathfrak{h}$ only ¹ . The idea of the Post-Minkowski approximation scheme is to construct an iterative solution for (2.24). For this purpose one chooses that in lowest order the gothic metric agrees with the Minkowski metric and for higher orders considers an expansion of $\mathfrak{h}^{\mu\nu}$ in terms of powers of κ according to

$$\begin{align} \mathfrak{g}^{\mu\nu}=\eta^{\mu\nu}+ \sum\limits_{n=1}^\infty \kappa^n \, \mathfrak{h}^{\mu\nu}_{(n)}\,. \end{align} \tag{ 2.26 }$$

As in this work we linearise our entire system around a flat Minkowski background and are therefore only interested in weak gravity, for the model in this article we can consider truncations of this perturbative expansion ² .

An iterative solution of Einstein's equations is constructed by the following procedure. Let us assume one has a solution for the potentials $\mathfrak{h}$ of order κ^k . Then one is able to compute the contributions of the quantities on the right hand side of (2.24) with it. Due to the prefactor κ, the right hand side will then be of order $\kappa^{k+1}$, so the entire equation will yield a solution for $\mathfrak{h}_{(k+1)}$. As can be seen from (2.26), the perturbative ansatz in the Post-Minkowski approximation scheme is constructed in a way that in the zeroth order where $\mathfrak{h} = 0$ one has $\mathfrak{g} = \eta$, from which one directly obtains $g = \eta$. With this, one can construct the right hand side of (2.24) using this zeroth order solution. As t_LL and t_H both contain $\mathfrak{h}$ in each summand, they vanish for $\mathfrak{h} = 0$ and the only non-trivial contribution comes from the energy-momentum tensor $T^{\mu\nu}$ on the right hand side. As in zeroth order $g = \eta$, it simplifies to the energy momentum tensor $T^{\mu\nu}$ on a flat Minkowski spacetime. Given the right hand side of (2.24), the left hand side amounts to the first order modified gothic metric, that is

$$\begin{align} \kappa\Box \mathfrak{h}^{\mu\nu}_{(1)} = -2\kappa T^{\mu\nu}[\Phi,\eta]\,. \end{align} \tag{ 2.27 }$$

Finally, one needs to relate $\mathfrak{h}^{\mu\nu}_{(1)}$ to the (inverse) metric perturbations $\delta h^{\mu\nu}$. Since in this work we are interested in the linearised theory only, we can use $g^{\mu\nu} = \eta^{\mu\nu}-\kappa \delta h^{\mu\nu}+O(\kappa^2)$ and obtain

$$\begin{align} \mathfrak{h}_{(1)}^{\mu\nu} = \delta h^{\mu\nu} - \frac{1}{2} \eta^{\mu\nu} \, \eta_{\rho\sigma} \delta h^{\rho\sigma}\,. \end{align} \tag{ 2.28 }$$

The last equations allows us to rewrite the linearised equations in a first order Minkowski approximation as

$$\begin{align} -\frac{\kappa}{2}\left( \Box (\delta h^{\mu\nu}) - \frac{1}{2} \eta^{\mu\nu} \, \eta_{\rho\sigma}\Box (\delta h^{\rho\sigma}) \right)= \kappa T^{\mu\nu}[\Phi,\eta] \,. \end{align} \tag{ 2.29 }$$

From (2.29) we can directly read off how the interaction term in the action needs to look like, namely

$$\begin{align} S_{I} = \frac{\kappa}{2} \int_{M} d^4x \; \delta h_{\mu\nu} T^{\mu\nu}[\Phi,\eta] \,. \end{align} \tag{ 2.30 }$$

The overall factor $\frac{\kappa}{2}$ takes into account that we start with an overall factor of $\frac{1}{2\kappa}$ and then consider linear perturbations of the metric $g_{\mu\nu} = \eta_{\mu\nu}+\kappa\delta h_{\mu\nu}$ yielding a factor κ² for second-order perturbations in the vacuum case. Note that, as expected, the same interaction contribution to the action can be obtained by minimally coupling matter to gravity and then linearising the metric degrees of freedom around a Minkowski spacetime [35–37, 100]. Let us denote the matter Lagrangian density by ${\mathcal L}_\Phi$, then we have

$$\begin{align} \int_M d^4x \; \sqrt{-\det(g)}\, \mathcal{L}_\Phi &\approx \int_M d^4x \; \mathcal{L}_\Phi \rvert_{g=\eta} + \int_M d^4x \; (-\kappa \delta h^{\mu\nu})\left(\frac{\delta (\sqrt{-\det(g)} \, \mathcal{L}_\Phi)}{\delta g^{\mu\nu}}\right) \Big|_{g=\eta} \nonumber\\ &= S^\eta_{\Phi} +\frac{\kappa}{2} \int_M d^4x\;\delta h^{\mu\nu} T_{\mu\nu}[\Phi,\eta] \,, \end{align} \tag{ 2.31 }$$

where $S^\eta_{\Phi}$ denotes the matter action on flat spacetime. For the model discussed in this work we choose the matter Lagrangian of a Klein–Gordon scalar field

$$\begin{align} {\mathcal L}_\phi=-\frac{1}{2}g^{\mu\nu}D_\mu\phi D_\nu \phi-\frac{m^2}{2}\phi^2, \end{align} \tag{ 2.32 }$$

where D_µ denotes the torsion-free, covariant metric-compatible derivative with respect to g. Performing an ADM decomposition as we did for the vacuum case we end up with

$$\begin{align} S_\phi = \int dt \int_\sigma d^3x \left( \dot{\phi} \pi - \left[ N^a C_a^\phi + N C^\phi \right] \right) \end{align} \tag{ 2.33 }$$

with

$$\begin{align} C_a^\phi := \pi \nabla_a \phi, \end{align} \tag{ 2.34 }$$

$$\begin{align} C^\phi := \frac{\pi^2}{2\sqrt{\det q}} + \frac{\sqrt{\det q}}{2} q^{ab} (\nabla_a \phi)(\nabla_b \phi) + \sqrt{\det q}\; \frac{m^2}{2}\phi^2\,, \end{align} \tag{ 2.35 }$$

where π denotes the canonically conjugate momentum associated to φ. If we linearise the model including the scalar field as discussed above and take into account all contributions up to first order in κ, we obtain the following constraints in the linearised theory:

$$\begin{align} \delta G_i &= \frac{\kappa}{2\beta} \left({{\partial}{_a}} ({{\delta E}{^a{}_i}}) + {{\epsilon}{_{ik}{}^l}} {\delta}_l^a\, {{\delta A}{_a{}^k}} \right), \end{align} \tag{ 2.36 }$$

$$\begin{align} \delta C_a &= \frac{\kappa}{\beta} \delta_j^b \left({{\partial}{_a}} ({{\delta A}{_b{}^j}}) - {{\partial}{_b}} ({{\delta A}{_a{}^j}}) \right) + \kappa \underbrace{\pi\,{{\partial}{_a}}\phi}_{=: p_a(\phi,\pi)}, \end{align} \tag{ 2.37 }$$

$$\begin{align} \delta C &= \kappa {{\epsilon}{_j{}^{kl}}} \delta_k^a \delta_l^b \,{{\partial}{{}_a}} ({{\delta A}{_b{}^j}}) + \kappa \underbrace{\frac{1}{2} \left[ \pi^2 +{{\partial}{_a}} \phi \,\partial^a \phi + m^2 \phi^2 \right]}_{=:\epsilon(\phi,\pi)} \,. \end{align} \tag{ 2.38 }$$

These constraints are also consistent with the 00 and 0a components of the linearised Einstein's equations that one obtains in the Post-Minkowski approximation scheme. Now these include gravitational- as well as matter contributions, where we introduced the momentum density $p_a(\phi,\pi)$ and the energy density $\epsilon(\phi,\pi)$ of the scalar field. Note that these constraints are abelian up to first order in κ. The corresponding background constraints of the zero order Post-Minkowski approximation scheme all vanish trivially, since $\overline{C} = \overline{C}^\textrm{geo} = 0, \overline{C}_a = \overline{C}^\textrm{geo}_a = 0$ and $\overline{G}_j = 0$. The action of the linearised theory is given by

$$\begin{align} S_\textrm{lin} = \int dt \Bigg(&\int_\sigma d^3x\left[\frac{\kappa}{\beta} {{\delta \dot{A}}{_a{}^i}} {{\delta E}{^a{}_i}} +\pi \dot{\phi} \right]\nonumber\\ &-\underbrace{\int_\sigma d^3x\bigg[ \epsilon(\phi,\pi) + \delta N \delta C + \delta N^a \delta C_a + \delta \Lambda^i \delta G_i + \delta \mathcal{H}_I +\frac{1}{\kappa} \delta^2 C^{geo} \bigg]}_{= \delta \mathbf{H}_\textrm{can}}\Bigg) \end{align} \tag{ 2.39 }$$

with the interaction part $\delta \mathcal{H}_I: = \frac{\kappa}{2}\delta h_{ab}T^{ab}$ with $\delta h_{ab}$ understood as a function of the densitised triads and the second order of the geometrical Hamiltonian constraint $\delta^2 C^{geo}$. Note that the overall $\frac{1}{\kappa}$ in front of the geometrical term has been partly cancelled by the κ involved in the linearised quantities for $N, N^a$ and $\Lambda^i$, while the constraints and the interaction term are still linear in κ. In the above form of the linearised action it was already used that in the linearised framework with asymptotically flat boundary conditions the boundary term precisely cancels the first order of the (geometrical) Hamiltonian constraint in accordance with [69], i.e. $T = -\overline{N}\delta C^{geo} = -\delta C^{geo}$ where we used that $\overline{N} = 1$. Hence, neither the linearised Hamiltonian constraint with background lapse nor the boundary term do appear in the action and canonical Hamiltonian anymore.

The linearised total canonical Hamiltonian, which follows from (2.9), is given by:

$$\begin{align} \delta \mathbf{H}_\textrm{can} =\int_\sigma d^3x \bigg[ &\epsilon(\phi,\pi) + \kappa \, \delta N^a\, \pi\, \partial_a \phi +\kappa\, \delta N\, \epsilon(\phi,\pi) +\kappa\, ({{\partial}{{}_a}}\phi) \,({{\partial}{{}^b}}\phi)\, {{\delta E}{^a{}_i}}\, \delta^i_b \nonumber\\ &+\frac{\kappa}{2} m^2 \phi^2\, {{\delta E}{^a{}_i}}\,\delta^i_a - \frac{\kappa}{2} {{\delta E}{^a{}_i}}\,\delta_a^i\, \epsilon(\phi,\pi)+ \frac{\kappa}{\beta} \delta N^a\, \delta_j^b \left({{\partial}{{}_a}} ({{\delta A}{_b{}^j}}) - {{\partial}{_b}} ({{\delta A}{_a{}^j}})\right) \nonumber\\ &+ \frac{\kappa}{2\beta} \delta \Lambda^i\,{{\partial}{_a}} ({{\delta E}{^a{}_i}}) + \frac{\kappa}{2\beta} \delta \Lambda^i\, {{\epsilon}{_{ik}{}^l}}\, \delta^a_l\, {{\delta A}{_a{}^k}} - \kappa\, {{\epsilon}{_j{}^{kl}}}\, \delta_k^a\, \delta_l^b\, {{\delta A}{_b{}^j}}\,{{\partial}{_a}} (\delta N)\nonumber\\ &- \frac{\kappa}{2\beta^2} \epsilon^{jkl} {{\epsilon}{_{jmn}}}\, \delta^a_k\, \delta^b_l \left( {{\delta A}{_a{}^m}} {{\delta A}{_b{}^n}}+ (\beta^2+1) {{\delta \Gamma}{_a{}^m}} {{\delta \Gamma}{_b{}^n}} -2 {{\delta \Gamma}{_a{}^m}} {{\delta A}{_b{}^n}} \right) \bigg]. \end{align} \tag{ 2.40 }$$

Let us briefly comment on the individual contributions in different orders of κ. The κ⁰-order encodes the equations of motion for the uncoupled system and environment. In our case these are the Klein–Gordon equation and the equations of motion for the background connection and densitised triad variables. Since the connection vanishes in the background and the triad is just given by $\overline{E}^{\,a}_i = \delta^a_i$ the corresponding equations of motion trivially vanish and this is again consistent with the vanishing of the background constraints that generate this trivial dynamics for the gravitational degrees of freedom. Note that the energy density of the scalar field $\epsilon(\phi,\pi)$ in κ⁰-order is not part of the background constraints but just contributes to the non-vanishing part of the Hamiltonian in this order. Without the boundary term a further term that contributes to the κ⁰ is the linearised Hamiltonian given by $\overline{N}\delta C^\textrm{geo} = \overline{N}(\delta C-\kappa\epsilon)$. On the linearised phase space this term generates the background equations. If we compute the Poisson brackets of $\delta{A}{_a{}^i},\delta {E}{^a{}_i}$ with $\overline{N}\delta C^\textrm{geo}$ they both vanish for $\overline{N} = 1$ demonstrating again that the background equations are trivially fulfilled. In linear order in κ in the covariant case we obtain a part that is quadratic in the perturbations of all metric components. Here, this corresponds to the terms that are either quadratic in the linearised Ashtekar variables or involve them linearly together with the linearised lapse and shift variables. These terms together in $\delta \mathbf{H}_\textrm{can}$ will generate the left hand side of all linearised Einstein's equations including the 00 and 0a component. The right hand side of the linearised Einstein's equations will be obtained from all contributions that involve gravitational as well as matter variables in the κ¹-order. These will be $\delta\mathcal{H}_I$ together with the parts where the matter variables occur in combination with the linearised lapse and shift.

As can easily be seen, most of the twenty phase space variables of the linearised theory,

$$\begin{align} (\phi,\pi) \hspace{0.4in} (\delta {A}{_a{}^i},\delta {E}{^b{}_j})\,, \end{align} \tag{ 2.41 }$$

are not observable, as they transform non-trivially under the linearised constraints. In the next section 2.4 we use the relational formalism to construct a set of independent gauge invariant quantities.

As mentioned above in a gauge theory we can either apply a specific gauge fixing or formulate the model in terms of gauge invariant quantities. In case we work at the level of perturbation theory, as done here, one often takes the approach that gauge invariance is guaranteed up to possible corrections that are of higher order in perturbation theory than one is truncating at, we will follow the same strategy here. To construct gauge invariant variables on the linearised phase space we consider the relational formalism [63–65] together with an observable map that maps the elementary variables of the linearised phase space to their corresponding gauge invariant quantities [66–68, 101]. The formalism is based on a choice of reference fields, one for each constraint in the system, and then the observable map returns for a given tensor field of certain rank on phase space its corresponding gauge invariant extension. Explicitly, it can be written as a power series in the reference fields weighted by contributions that involve nested Poisson brackets of the tensor fields and the constraints. The physical interpretation of these gauge invariant quantities, also known as Dirac observables in the framework of general relativity, is that they give the value of the tensor field at those values where the reference fields take specific values. For general relativity and its spatial diffeomorphism and Hamiltonian constraints this particularly means that we can formulate the dynamics of a given tensor field on phase space with respect to the reference fields and these provide a notion of physical spatial and temporal coordinates. The observable map introduced in [66–68] has the property that when applied to the chosen reference fields it maps them to phase space independent quantities and thus such a map can describe neither a canonical transformation on the full phase space one starts with nor a transformation that preserves the number of elementary phase space variables. For this reason we will slightly modify the map from [66–68] in two aspects. First we modify its action on reference fields and second we further introduce also the kind of dual version of this map to treat the choice of the reference fields and the constraints more on an equal footing, allowing to use the modified map to construct a canonical transformation for the variables $(\phi,\pi), (\delta {A}{_a{}^i},\delta {E}{^b{}_j})$. For the purpose of this work we can restrict our discussion to the linearised phase space. The transformation that we aim at constructing will map the set of elementary variables $(\phi,\pi), (\delta {A}{_a{}^i},\delta {E}{^b{}_j})$ to a new set of canonical variables containing the Dirac observables associated with $(\phi,\pi)$, seven chosen reference fields that are canonically conjugate to the seven linearised constraints as well as two canonical pairs of further Dirac observables in the gravitational sector. Using this kind of canonical transformation for the elementary phase space variables, the physical phase space can be more easily accessed because the constraints as well as the reference fields are among the new canonical variables.

As we will discuss below for a modification of the observable maps in [66–68, 101] in our work here, we add a so called dual observable map and the combination of both allow us to construct the canonical transformation on the full phase space. Let us briefly introduce the observable map from [66–68, 101] as well as the dual one and then apply them to the model. We consider a system with a set of first class constraints $\{C_I\}$ and elementary phase space variables $(q^A,p_A)$. Then we choose a set of reference variables $\{T^I\}$ that satisfy $\det(\{T^I,C_J\})\not = 0$ where we abbreviate $\{T^I,C_J\}: = M^{I}_{\,\,J}$. We can define an equivalent set of constraints given by

$$\begin{align} C^{\prime}_I=\sum\limits_{J}(M^{-1})^J_{\,\,I}C_J. \end{align} \tag{ 2.42 }$$

Then by construction we have that the constraints $C_I^{\prime}$ are weakly canonically conjugate to the T^I 's. Given this, one further chooses a set of functions $\{\tau^I\}$ that can depend on the spatial and temporal coordinates by means of which one can construct an observable map which maps a function f on the phase space to its corresponding observable denoted by ${\mathcal O}_{f,\{T\}}$

$$\begin{align} {\mathcal O}_{f,\{T\}}(\tau^I)&=\left[\exp( \xi^I\{C_I,\cdot\})\cdot f\right]\Big|_{\xi^I:=T^I-\tau^I}\\ &=f+(T^I-\tau^I)\{C_I,f\}+\frac{1}{2!}(T^I-\tau^I)(T^J-\tau^J)\{C_J,\{C_I,f\}\}+\ldots\, ,\nonumber \end{align} \tag{ 2.43 }$$

where the label $\{T\}$ refers to the chosen set of reference variables. The observable ${\mathcal O}_{f,\{T\}}(\tau^I)$ returns the value of f at those values where the reference variables T^I take the values τ^I . As can be shown [66, 68] we have $\{C_I,O_{f,T}\}\approx 0$ where ≈ denotes weak equivalence, that is on the constraint hypersurface defined by the C_I 's and thus the $O_{f,\{T\}}$ are weak Dirac observables. In case the constraints C_I are already canonically conjugate to the T^I 's then the $O_{f,\{T\}}$ are even strong Dirac observables and weak equalities are replaced by equality signs. Let us introduce ${\mathcal G}^I: = T^I-\tau^I$ that we can also understand as a choice of coordinate gauge fixing condition if we require ${\mathcal G}^I\approx 0$. Then (2.43) can be understood as a power series in the ${\mathcal G}^I$'s with nested Poisson brackets involving the constraints C_I that are in this case equal to the $C^{\prime}_I$'s. In case we have $\{T^I,C_J\} = \delta^I_J$ because the τ^I do not depend on phase space variables we also obtain $\{{\mathcal G}^I,C_J\} = \delta^I_J$. Thus, the gauge fixing conditions and the constraints build canonically conjugate pairs. If this is the case, we can construct a dual version of the observable map in (2.43) where the role of the gauge fixing conditions and constraints are interchanged given by

$$\begin{align} {\mathcal O}^\textrm{dual}_{f,\{C\}} &=\left[\exp( -\xi_I\{{\mathcal G}^I,\cdot\})\cdot f\right]\Big|_{\xi_I:=C_I}\\ &=f-C_I\{{\mathcal G}^I,f\}+\frac{1}{2!}C_I C_J\{{\mathcal G}^J,\{f, {\mathcal G}^I \}\}\pm\cdots\nonumber \end{align} \tag{ 2.44 }$$

where the label $\{C\}$ denotes that we have interchanged the role of gauge fixing conditions and the constraints for the dual map. This is also the reason for the additional factor of $(-1)^n$ compared to the observable map in (2.43). Note that this dual map does not depend on any functions τ^I since the constraint hypersurface is defined by $C_I\approx 0$. Similar to ${\mathcal O}_{f,\{T\}}(\tau^I)$ which can be understood as a family of gauge invariant extension of f parameterised by τ^I for ${\mathcal G}^I\not = 0$, the quantity ${\mathcal O}^\textrm{dual}_{f,\{C\}}$ is an extension from the $C_I\approx0$ constraint hypersurface with the property that it commutes with all the ${\mathcal G}^I$ by construction.

Our strategy is to construct a canonical transformation to a new set of elementary variables. For this purpose, we apply a combination of the two observable maps yielding quantities that commute by construction with all constraints C_I and all ${\mathcal G}_I$'s. Hence, we can choose the set $({\mathcal G}^I,C_I)$ as new canonical variables together with the number of independent Dirac observables obtained from the observable and its dual map associated with the set $(q^I,p_I)$. How many independent such Dirac observables exist depends on the number of constraints.

Note that the introduction of the equivalent set of constraints in (2.42) is also called weak Abelianisation [66] and ensures that the corresponding Hamiltonian vector fields of the C_I 's weakly commute. As a consequence, the order in which we apply the constraints C_I in the nested Poisson bracket is irrelevant. Now the situation relevant for us is that the constraints satisfy $\{T^I,C_J\} = \delta^I_J$ and thus the C_I 's can be expressed linearly in the momenta conjugate to the T^I 's. Using that the set $\{C_I\}$ is first class, this is sufficient to show that the constraints are abelian and hence their corresponding Hamiltonian vector fields commute and the order how we apply the C_I also does not matter in our case. However, if we consider the entire set $({\mathcal G}^I,C_I)$ then the ${\mathcal G}^I$ and C_I do not commute, not even weakly and neither is the entire set first class and in general also the subset of the ${\mathcal G}^I$ is not abelian. This has the effect that the final observable that we obtain by applying the observable map in combination with its dual in general depends on the order in which we apply the two maps as well as the order in which the ${\mathcal G}^I$ occur in the nested Poisson brackets. This causes no problem for the model considered in this article, it just means that there exist different choices of possible coordinate transformations on phase space, but important for us is rather that we can choose one among those.

In the following we will generalise the observable map to the field theoretic setup and perturbations theory which can easily be done. In the context of perturbation theory it is sufficient to require gauge invariance or some specific form of the Poisson algebra up to corrections that are higher than the order that is considered in perturbation theory. This means we can truncate the power series for the map and its dual at some order in accordance with the desired order in perturbation theory that we consider for the linearised model. For instance if we want to compute Poisson brackets of observables up to linear order then we need to perturb both observables up to second order and collect all terms that contribute up to linear order to the final result. Such a perturbative approach for constructing observables has for instance also been used in [69–71]. Alternatively, if available we can also take the result of the corresponding Poisson bracket in full general relativity and perturb it up to linear order. But since many quantities we work with will be known at the perturbative level only the latter option is often not possible.

The first step in the relational formalism consist of choosing suitable reference fields for the given set of constraints. In our case these will be the Hamiltonian, spatial diffeomorphism and Gauß constraint. If we are interested in results up to corrections of second order in the perturbations inside the Poisson brackets we need to consider perturbations up to second order of these constraints. The properties our chosen reference fields should have are:

• (i)  
  Each chosen reference field consists of linear perturbations of the elementary gravitational variables and its derivatives only. These are also known as linearised geometrical clocks.
• (ii)  
  Each chosen reference field is in lowest order canonically conjugated to one of the constraints.
• (iii)  
  Each reference field commutes in lowest order with all constraints except the one that it is in lowest order canonically conjugated to.
• (iv)  
  All reference fields mutually commute.

Let us introduce the notation

$$\begin{align*} \delta{\mathcal C}_I(\vec{x},t):=(\delta C(\vec{x},t), \delta C_a(\vec{x},t), \delta G_i(\vec{x},t)), \quad \delta^2{\mathcal C}_I(\vec{x},t):=(\delta^2 C(\vec{x},t), \delta^2 C_a(\vec{x},t), \delta^2 G_i(\vec{x},t)) \end{align*}$$

for the set of linearised constraints $\delta{\mathcal C}_I(\vec{x},t)$ and $\delta^2{\mathcal C}_I(\vec{x},t)$ for the set of second-order perturbations. For each of the individual constraints we have to choose one reference field and we introduce the following notation ${\mathcal G}^I(\vec{x},t)$ for this set with

$$\begin{align} & {\mathcal G}^I(\vec{x},t):=(\delta{\mathcal G}(\vec{x},t),\delta{\mathcal G}^a(\vec{x},t),\delta{\mathcal G}^j(\vec{x},t))\nonumber \\ &\delta{\mathcal G}(\vec{x},t):=\delta T(\vec{x},t)-\delta\tau(\vec{x},t), \quad \delta {\mathcal G}^a(\vec{x},t):=\delta T^a(\vec{x},t)-\delta\sigma^a(\vec{x},t),\nonumber \\ &\delta {\mathcal G}^j(\vec{x},t):=\delta\Xi^j(\vec{x},t)-\delta\xi^j(\vec{x},t), \end{align} \tag{ 2.45 }$$

where $\delta T, \delta T^a, \delta \Xi^j$ denote the individual reference fields for the Hamiltonian, spatial diffeomorphism and Gauß constraint respectively. We denote background quantities with a bar and assume that in the background the gauge fixing conditions as well as the constraints are satisfied, that is

$$\begin{align*} &\overline{\mathcal G}(\vec{x},t):=\overline{T}(\vec{x},t)-\overline{\tau}(\vec{x},t)=0, \overline{C}=0,\quad \overline{\mathcal G}^a(\vec{x},t):=\overline{T}^a(\vec{x},t)-\overline{\sigma}^a(\vec{x},t)=0, \overline{C}_a=0,\nonumber \\ &\overline{\mathcal G}^j(\vec{x},t):=\overline{T}^j(\vec{x},t)-\overline{\xi}^j(\vec{x},t)=0, \overline{G}_j=0. \end{align*}$$

With the assumptions (i)-(iv) above we know that $\delta^2{\mathcal G}^I = 0$ and further we have

$$\begin{align} \{{\mathcal G}^I(x),{\mathcal C}_J(y)\} :&=\{\delta{\mathcal G}^I(x),\delta {\mathcal C}_J(y)\} +\{\delta{\mathcal G}^I(x),\delta^2 {\mathcal C}_J(y)\} +O(\delta^2,\kappa^2)\nonumber \\ &=\frac{1}{\kappa}\delta^I_J\delta^{(3)}(x,y) +\{\delta{\mathcal G}^I(x),\delta^2 {\mathcal C}_J(y)\} +O(\delta^2,\kappa^2), \end{align} \tag{ 2.46 }$$

where $O(\delta^2,\kappa^2)$ means that we neglect all terms that are second-order in the perturbations and/or of order κ² and the factor $\frac{1}{\kappa}$ has been chosen because it is also involved in the Poisson bracket of the elementary gravitational variables. This means there could be terms being of order κⁿ with $n\leqslant 1$ that do not contribute because they involve second or higher orders of δ. That those terms can be present is caused by an asymmetry in δ and κ due to the fact that we only perturb the gravitational degrees of freedom but not the matter ones and further that the individual terms in the action involve different powers of κ from the beginning. We realise that the linearised constraints $\delta{\mathcal C}_I$ are canonically conjugated to the reference fields ${\mathcal G}^I$ but in general there can be non-vanishing contribution in linear order coming from the Poisson bracket $\{\delta{\mathcal G}^I(x),\delta^2 {\mathcal C}_J(y)\}$. To ensure that we have a vanishing contribution in linear order we will use the dual observable map as discussed below. For this purpose we need to adapt the observable map and its dual to field theory and perturbation theory and for both maps we need the observable map up to second order. Taking this into account the observable formula up to second order reads

$$\begin{align} {\mathcal O}_{f,\{T\}}(\delta\tau,\delta\sigma^a,\delta\xi^j)& = \delta f +\delta^2f +\kappa\int\limits_{\sigma}d^3y\;\delta{\mathcal G}^I(y)\{\delta {\mathcal C}_I(y),\delta f\}\nonumber\\ & \quad + \kappa\int\limits_{\sigma}d^3y\;\delta^2{\mathcal G}^I(y)\{\delta {\mathcal C}_I(y),\delta f\}\nonumber\\ & \quad + \kappa\int\limits_{\sigma}d^3y\;\delta{\mathcal G}^I(y)\Big(\{\delta^2 {\mathcal C}_I(y),\delta f\}+\{\delta {\mathcal C}_I(y),\delta^2 f\}\Big)\nonumber \\ & \quad + \frac{\kappa^2}{2}\int\limits_{\sigma}d^3y\;\delta{\mathcal G}^I(y)\int\limits_{\sigma}d^3z\;\delta{\mathcal G}^J(z)\Big(\{\delta {\mathcal C}_J(z),\{\delta^2 {\mathcal C}_I(y),\delta f\}\}\nonumber\\ & \quad +\{\delta {\mathcal C}_J(z),\{\delta{\mathcal C}_I(y),\delta^2f\}\}\Big)+O(\delta^3,\kappa^3), \end{align} \tag{ 2.47 }$$

where we used that $\overline{\mathcal G} = \overline{\mathcal G}^a = \overline{\mathcal G}^j = 0$ and allowed possible non-vanishing second-order perturbations of ${\mathcal G}^I$ that might be present if one wants to relax the assumptions (i)–(iv) from above. Because we consider perturbations around flat spacetime, in the background it holds that $\overline{C} = 0, \overline{C}_a = 0$ as well as $\overline{G}_i = 0$, all trivially vanish and hence we also have $\overline{\mathcal G} = 0, \overline{\mathcal G}^a = 0$ and $\overline{\mathcal G}^j = 0$. Hence, for the background there is no gauge freedom we have to deal with and thus no corresponding observables to construct. Note that if we construct the observable for elementary phase space variables then $\delta^2f = 0$ and the observable formula above simplifies. This construction of observables order by order in κ plays also a role when we consider the linearised Hamiltonian $\delta\textbf{H}_\textrm{can}$ which, as can be seen from (2.40), has contributions in κ⁰ and linear order in κ. The transformation behaviour under the linearised constraints of the matter variables can be found in the appendix A.1 in (A.7) and (A.11) and the results are an expression linear in κ. This again demonstrates that in the limit where κ is sent to zero, which corresponds to the situation that we consider a scalar field on Minkowski only with no coupling to linearised gravity, the elementary variables $\phi,\pi$ are suitable observables. Once we consider the coupling with linearised gravity perturbations, this is no longer the case and we need the gauge invariant version of these variables. For the $\delta\textbf{H}_\textrm{can}$ this means that in the κ⁰ term we can still work with the original $\phi,\pi$ whereas for the linear order in κ we also need to construct Dirac observables for the matter sector, see for instance also the discussion in [100] in the context of the covariant theory.

Note that as before, the linearised phase space in the Post-Minkowski approximation scheme involves the linear perturbations of the gravitational degrees of freedom as well as the variables $(\phi,\pi)$ for the matter sector. In a similar way we obtain the linearised dual observable map given by

$$\begin{align} {\mathcal O}^\textrm{dual}_{f,\{C\}}& = \delta f +\delta^2f -\kappa\int\limits_{\sigma}d^3y\;\delta{\mathcal C}_I(y)\{\delta {\mathcal G}^I(y),\delta f\}\nonumber\\ & \quad - \kappa\int\limits_{\sigma}d^3y\;\delta^2{\mathcal C}_I(y)\{\delta {\mathcal G}^I(y),\delta f\}\nonumber\\ & \quad - \kappa\int\limits_{\sigma}d^3y\;\delta{\mathcal C}_I(y)\Big(\{\delta^2 {\mathcal G}^I(y),\delta f\}+\{\delta {\mathcal G}^I(y),\delta^2 f\}\Big)\nonumber\\ & \quad + \frac{\kappa^2}{2}\int\limits_{\sigma}d^3y\;\delta{\mathcal C}_I(y)\int\limits_{\sigma}d^3z\;\delta{\mathcal C}_J(z)\Big(\{\delta {\mathcal G}^J(z),\{\delta^2 {\mathcal G}^I(y),\delta f\}\}\nonumber\\ & \quad + \{\delta {\mathcal G}^J(z),\{\delta{\mathcal G}^I(y),\delta^2f\}\}\Big)+O(\delta^3,\kappa^3). \end{align} \tag{ 2.48 }$$

The alternating signs compared to the observable map in (2.43) are needed since the order of how the constraints and the reference fields enter into the Poisson bracket is switched. Note that if we choose for instance ${\mathcal G}^I$ such that $\delta^2{\mathcal G}^I = 0$ and not the non-linear ${\mathcal C}_I$ but only the linearised parts as the clocks then we drop all terms that involve $\delta^2{\mathcal G}^I,\delta^2{\mathcal C}_I$ in the dual observable map (2.48) and this will be exactly the application we need later.

As we will work with so-called geometrical clocks, in a first step we will choose the set of reference fields among the linearised elementary gravitational degrees of freedom. As a consequence, we will work with geometrical clocks for which $\delta^2{\mathcal G}^I = 0$ and hence they satisfy assumption (i). To also fulfil (ii) we further choose ${\mathcal G}^I$ such that

$$\begin{align*} \{\delta{\mathcal G}^I(x), \delta C_J(y)\} = \frac{1}{\kappa}\delta^I_J\delta^{(3)}(x,y), \end{align*}$$

where the factor $\frac{1}{\kappa}$ is chosen because it is also involved in Poisson bracket of the elementary geometrical variables. These requirements still allow several choices for suitable ${\mathcal G}^I$ and the specific choice for ${\mathcal G}^I$ taken here is motivated by the fact that we can relate the set of gauge invariant variables to a gauge fixing often chosen in the context of linearised gravity as will be discussed at the end of this subsection. For the reference fields $\delta T(\vec{x},t), \delta\Xi^j(\vec{x},t)$, for the Hamiltonian constraint and the Gauß constraint respectively we choose

$$\begin{align} \delta T(\vec{x},t) &:= -\frac{1}{\kappa\beta} \left[\frac{1}{2} \delta_c^i {{\epsilon}{_a{}^{cb}}} \partial_b \left({{\delta E}{^a{}_i}} * G^\Delta \right)(\vec x,t) + \delta^a_i \left({{\delta A}{_a{}^i}} * G^\Delta\right)(\vec{x},t)\right]\,, \end{align} \tag{ 2.49 }$$

$$\begin{align} \kern-3.5pt \delta \Xi^i(\vec{x},t)&:= \frac{2}{\kappa} \partial^a \left( {{\delta A}{_a{}^i}} * G^\Delta \right)(\vec{x},t), \end{align} \tag{ 2.50 }$$

with the three dimensional, spatial convolution

$$\begin{align} \left({{\delta E}{^a{}_i}} * G^\Delta \right)(\vec x,t) := \int d^3y \; {{\delta E}{^a{}_i}}(\vec y,t) \, G^\Delta(\vec{x}-\vec y)\,, \end{align} \tag{ 2.51 }$$

the derivative of the convolution acting with respect to $\vec x$, hence only on $G^\Delta(\vec x -\vec y)$, and the Green's function of the Laplacian

$$\begin{align} G^\Delta(\vec x-\vec y) := \int \frac{d^3k}{(2\pi)^3} \frac{1}{||\vec{k}||^2} e^{i\vec{k}(\vec{x}-\vec{y})}\,, \end{align} \tag{ 2.52 }$$

such that $\Delta G^\Delta(\vec x-\vec y) = -\delta(\vec x -\vec y)$. Thus the reference fields are, as in general expected, non-local quantities because they are local in momentum space. As discussed in the appendix A.1 these reference fields satisfy

$$\begin{align*} \{\delta T(\vec{x},t),\delta C(\vec{y},t)\} &= \frac{1}{\kappa}\delta(\vec{x}-\vec{y}) ,\quad \{\delta \Xi^i(\vec{x},t),\delta G_j(\vec{y},t)\} = \frac{1}{\kappa}\delta_j^i \delta(\vec{x}-\vec{y}), \quad \{\delta T(\vec{x},t), \delta \Xi^i(\vec{y},t)\} =0\,. \end{align*}$$

Furthermore, the reference fields $\delta T,\delta \Xi^j$ commute with all remaining linearised constraints and hence they satisfy (iii) and (iv) as well. A more detailed discussion can be found in the appendix A.1. Reference fields that are canonically conjugate to the components of the spatial diffeomorphism constraint and commute with the linearised Hamiltonian as well as the linearised Gauß constraint are given by

$$\begin{align} \delta\widetilde{T}^a(\vec{x},t):= \frac{1}{\kappa} \left( \delta^a_b \delta^i_c \partial^c -\frac{1}{2} \delta^i_b \partial^a + \delta^{ac}\delta_c^i \partial_b \right) \left(\delta {E}{^b{}_i} * G^\Delta\right)(\vec{x},t). \end{align} \tag{ 2.53 }$$

Although the reference fields $\delta\widetilde{T}^a(\vec{x},t)$ mutually commute, they do not have vanishing Poisson brackets with all the remaining reference fields and violate assumption (iv). To obtain reference fields for the spatial diffeomorphism constraint we can employ the dual observables map for $\delta\widetilde{T}^a(\vec{x},t)$ yielding a quantity that by construction commutes with $\delta{\mathcal G}, \delta{\mathcal G}^a, \delta{\mathcal G}^j$ and hence with all reference fields. As the clock for the dual observable map we choose the geometric part of the linearised Hamiltonian constraint, that is $\delta C-\kappa\epsilon$ as well as the linearised Gauß constraint. Note that $\delta C-\kappa\epsilon\not\approx 0$ and as a consequence $\delta\widetilde{T}^a\not\approx\delta T^a$ but this causes no problems because two different choices of set of clocks need not necessarily to be weakly equivalent. In order that the observable formula can be applied we just need $\{\delta T^I,\delta{\mathcal C}_J\} = \{\delta {\mathcal G}^I,\delta{\mathcal C}_J\} = \delta^I_J$ which is satisfied also if we use the geometric contribution of $\delta C$ only. For $\delta\widetilde{T}^a(\vec{x},t)$ the part corresponding to the spatial diffeomorphism constraint will not contribute in the dual observable map. Since we use the linearised constraint as the clock and not ${\mathcal C}_I$ we can drop all terms in involving $\delta^2C_I$. Using this together with $\delta^2{\mathcal G}^I = 0$ the dual observable map in (2.44) simplifies in this case to

$$\begin{align*} &{\mathcal O}^\textrm{dual}_{\delta\widetilde{T}^a,\{\delta C\}}=:\delta T^a = \delta \widetilde{T}^a -\kappa\int\limits_{\sigma}d^3y\;(\delta C-\kappa\epsilon)(y)\{\delta T(y),\delta \widetilde{T}^a\} -\kappa\int\limits_{\sigma}d^3y\;\delta G_j(y)\{\delta \Xi^j(y),\delta \widetilde{T}^a\}. \end{align*}$$

The explicit result we obtain has the form

$$\begin{align} \delta T^a(\vec{x},t) :=& \frac{1}{\kappa} \left( \delta^a_b \delta^i_c \partial^c -\frac{1}{2} \delta^i_b \partial^a + \delta^{ac}\delta_c^i \partial_b \right) \left(\delta {E}{^b{}_i} * G^\Delta\right)(\vec{x},t) \nonumber\\ &+\frac{2\beta}{\kappa^2}\left[ \frac{1}{2} \delta^i_b \partial^a \partial^b \left( \delta G_i * G^{\Delta\Delta}\right)(\vec x,t) - \delta^{ab}\delta_b^i \left( \delta G_i * G^\Delta\right)(\vec x,t) \right]\nonumber\\ & + \frac{1}{2\kappa^2} \partial^a \left[\left(\delta C-\kappa\epsilon\right) * G^{\Delta\Delta} \right](\vec{x},t)\,, \end{align} \tag{ 2.54 }$$

where $G^{\Delta\Delta}$ denotes the Green's function of the squared Laplacian, that is

$$\begin{align} G^{\Delta\Delta}(\vec x-\vec y) = \int \frac{d^3k}{(2\pi)^3} \frac{1}{||\vec k||^4} e^{i\vec{k}(\vec x-\vec y)} = \int d^3z \; G^\Delta(\vec x -\vec z) \, G^\Delta(\vec y-\vec z)\,. \end{align} \tag{ 2.55 }$$

Note that in [69] also geometrical clocks were used based on the ADM clocks introduced in [74, 102]. This set of clocks is not equivalent to ours and when requiring that all clocks vanish this corresponds to a different gauge fixing, as is discussed at the end of appendix B.1.

Now we are in the situation, that the seven linearised constraints as well as the seven reference fields form in lowest order two abelian subalgebras and further obey the following Poisson algebra:

$$\begin{align} \{\delta T(\vec{x},t),\delta C(\vec{y},t)\} = \frac{1}{\kappa}\delta(\vec{x}-\vec{y})\quad \{\delta T(\vec{x},t),\delta C_a(\vec{y},t)\} = 0 \quad \{\delta T(\vec{x},t),\delta G_i(\vec{y},t)\} = 0 \end{align} \tag{ 2.56 }$$

$$\begin{align} \{\delta T^a(\vec{x},t),\delta C(\vec{y},t)\} = 0\quad \{\delta T^a(\vec{x},t),\delta C_b(\vec{y},t)\} = \frac{1}{\kappa}\delta_b^a \delta(\vec{x}-\vec{y})\quad \{\delta T^a(\vec{x},t),\delta G_i(\vec{y},t)\} = 0\nonumber\\ \end{align} \tag{ 2.57 }$$

$$\begin{align} \{\delta \Xi^i(\vec{x},t),\delta C(\vec{y},t)\} = 0\quad \{\delta \Xi^i(\vec{x},t),\delta C_a(\vec{y},t)\} = 0\quad \{\delta\Xi^i(\vec{x},t),\delta G_j(\vec{y},t)\} = \frac{1}{\kappa}\delta_j^i \delta(\vec{x}-\vec{y})\,.\nonumber\\ \end{align} \tag{ 2.58 }$$

Since $\delta{\mathcal G},\delta{\mathcal G}^c,\delta{\mathcal G}^j$ differ from $\delta T,\delta T^c,\delta\Xi^j$ by some phase space independent quantity only, we can replace the reference fields by the corresponding $\delta{\mathcal G}^I$s and the Poisson algebra above will not change. We have defined the Poisson algebra with a factor $\frac{1}{\kappa}$ because the gravitational degrees of freedom involve a similar factor. The algebra reference fields have been defined such that the algebra of them and the linearised constraints no longer involves the Barbero-Immirzi parameter on the right hand side for the reason that κ is the parameter labelling the order of perturbations. Going back to the Poisson algebra of ${\mathcal G}^I$ and ${\mathcal C}_I$ in (2.46) so far we have chosen reference fields that satisfy this algebra. However, the in general non-vanishing contribution in linear order will prevent us from chosen ${\mathcal G}^I$ and ${\mathcal C}_J$ as new canonical coordinates. To achieve this we will apply the dual observable map to ${\mathcal C}_I$ and obtain an equivalent set of constraints ${\mathcal C}_I^{\prime}$ satisfying

$$\begin{align} \{{\mathcal G}^I(x),{\mathcal C}^{\prime}_J(y)\}=\{\delta{\mathcal G}^I(x),\delta {\mathcal C}^{\prime}_J(y)+\delta^2 {\mathcal C}^{\prime}_J(y)\}=\frac{1}{\kappa}\delta^I_J\delta^{(3)}(x,y) +O(\delta^2,\kappa^2). \end{align} \tag{ 2.59 }$$

With this algebra given we can then choose $\delta{\mathcal G}^I$ as new configuration variables and $\delta{\mathcal C}^{\prime}_I$ as new momentum variables. Because in lowest order we have already the correct algebra relations we only need a modification in the linear order. To accomplish this we apply the dual observable map to all ${\mathcal C}_I$ and drop the linear order terms. The latter would just abelianise the lowest order and cancel the ${\mathcal C}_I$ contribution in the observable something we do not want. The second order that we keep will then modify the linear order in a way that we obtain the algebra shown in (2.59). Since we want an equivalent set of constraints here we are forced to choose as clocks the linearised constraints $\delta{\mathcal C}_I$ involving also the matter contributions in $\delta C$ and $\delta C_a$. Taking all this into account the equivalent constraints ${\mathcal C}_I^{\prime}$ are given by

$$\begin{align} {\mathcal C}_I^{\prime}:=& {\mathcal O}^\textrm{dual}_{{\mathcal C}_I,\{\delta C\}}=\delta{\mathcal C}_I+\delta^2{\mathcal C}_I -\kappa\int\limits_{\sigma}d^3y\;\delta{\mathcal C}_I(y)\{\delta {\mathcal G}^I(y),\delta^2 {\mathcal C}_I\}\nonumber \\ & +\frac{\kappa^2}{2}\int\limits_{\sigma}d^3y\;\delta{\mathcal C}_J(y)\int\limits_{\sigma}d^3z \;\delta{\mathcal C}_K(z)\{\delta {\mathcal G}_J(z),\{\delta{\mathcal G}_K(y),\delta^2{\mathcal C}_I\}\} +O(\delta^3,\kappa^3), \end{align} \tag{ 2.60 }$$

where we used again that $\delta^2{\mathcal G}^I = 0$ in our case. Up to linear order we have $\delta{\mathcal C}_I^{\prime} = \delta{\mathcal C}_I$ showing that as expected we have modified the constraint only in second order. This modifications become not relevant when we work with the constraints directly in the linearised theory since they are of second order but need to be considered when we compute Poisson brackets and want to consider the result up to linear order. For a set of mutually commuting reference fields, that we consider here, in a more general context it was proven in [69] that the constraints ${\mathcal C}_I^{\prime}$ are abelian up to corrections of second order and applying this result here we have

$$\begin{align*} \{{\mathcal C}_I^{\prime}(x), {\mathcal C}_J^{\prime}(y)\}=\{\delta{\mathcal C}_I^{\prime}(x)+\delta^2{\mathcal C}_I^{\prime}(x), \delta{\mathcal C}_J^{\prime}(y)+\delta^2{\mathcal C}_J^{\prime}(y)\}=0+O(\delta^2,\kappa^2) \end{align*}$$

and thus in the order of perturbation theory we consider here, we can treat them as abelian constraints. Because we can just consider higher orders in the dual observable map we can always ensure that the constraints are abelian up to a chosen order in perturbation theory and if the entire sum of the dual observable map converges even in the full unperturbed theory. Note that a similar result has already been discussed in [69] in the context of applying weak abelianisation order by order in perturbation theory and afterwards modifying the constraints by terms that involve higher order powers of the linearised constraints. Here we can rediscover that case in the framework of the dual observable map which seems to us to be slightly more general because it cannot only be applied to the constraints as this corresponds to the case discussed in [69] but to any phase space function as we did for instance for the clock of the spatial diffeomorphism constraint before. Further, the strategy discussed in [69] separates weak Abelianisation and the addition of terms higher order in the linearised constraints where this happens all at once here using the dual observable map from second order on.

Given the abelian constraints ${\mathcal C}^{\prime}_I(x)$ now we can use them in the observable map in (2.43) and construct gauge invariant quantities that commute by construction with the constraints ${\mathcal C}_I^{\prime}(x)$ up to $O(\delta^2,\kappa^2)$ corrections. We want to apply the observable map to the elementary phase space variables of the matter and gravitational degrees of freedom and use $\delta{\mathcal G}^I$ as our clocks. Therefore, we can drop all terms involving $\delta^2f$ and $\delta^2{\mathcal G}^I$ and the observable map reduces to

$$\begin{align} (\delta f)^{GI}(\delta\tau,\delta\sigma^a,\delta\xi^j)&:={\mathcal O}_{\delta f,\{\delta T\}}(\delta\tau,\delta\sigma^a,\delta\xi^j)\\ &={\mathcal O}^{(1)}_{\delta f,\{\delta T\}}(\delta\tau,\delta\sigma^a,\delta\xi^j)+{\mathcal O}^{(2)}_{\delta f,\{\delta T\}}(\delta\tau,\delta\sigma^a,\delta\xi^j)+O(\delta^3,\kappa^3),\nonumber \end{align} \tag{ 2.61 }$$

where the label GI means gauge invariant with

$$\begin{align} {\mathcal O}^{(1)}_{\delta f,\{\delta T\}}(\delta\tau,\delta\sigma^a,\delta\xi^j)&:=\delta f +\kappa\int\limits_{\sigma}d^3y\;\delta{\mathcal G}^I(y)\{\delta {\mathcal C}^{\prime}_I(y),\delta f\}\nonumber\\ {\mathcal O}^{(2)}_{\delta f,\{\delta T\}}(\delta\tau,\delta\sigma^a,\delta\xi^j)&:= \kappa\int\limits_{\sigma}d^3y\;\delta{\mathcal G}^I(y)\{\delta^2 {\mathcal C}^{\prime}_I(y),\delta f\} \nonumber \\ & \quad +\frac{\kappa^2}{2}\int\limits_{\sigma}d^3y\;\delta{\mathcal G}^I(y)\int\limits_{\sigma}d^3z \;\delta{\mathcal G}^J(z)\Big(\{\delta {\mathcal C}^{\prime}_J(z),\{\delta^2 {\mathcal C}^{\prime}_I(y),\delta f\}\} \Big)\nonumber\\ \end{align} \tag{ 2.62 }$$

where ${\mathcal O}^{(1)}_{\delta f,\{\delta T\}}$ and ${\mathcal O}^{(2)}_{\delta f,\{\delta T\}}$ denote all contributions to the observable in linear and second order respectively. Since in this work we perturb the geometrical degrees of freedom only for the matter sector we just replace $\delta f$ by $f = \phi,\pi$ in the formula above. That this works so easily in the matter sector might not be obvious in the first place for the following reason: If we consider in ${\mathcal O}^{(1)}_{\delta f,\{\delta T\}}$ the second term then for $\delta f$ chosen from geometric degrees of freedom this will be of zeroth order and hence trivially commute with the linearised or higher order constraints. In contrast, for the matter variables the corresponding term still involves the matter variables linearly and hence does not commute with the linearised constraints or higher order ones and thus yields an additional contribution compared to the geometric case. However, since the Poisson brackets of the gravitational degrees of freedom involve a factor $\frac{1}{\kappa}$, whereas the Poisson brackets of the matter variables do not, these additional contributions come with an extra factor of κ. In the linearised theory we require for Dirac observables to commute with the constraints up to order δ and κ and hence, these kind of contributions will only occur in higher orders in κ. A similar argument applies to these kind of additional contributions in ${\mathcal O}^{(2)}_{\delta f,\{\delta T\}}$ and therefore we can indeed apply the observable formula also to the matter variables. To obtain the explicit form of these Dirac observables we apply the observable map in (2.62) to $(\phi,\pi,{{\delta A}{_a{}^i}},{{\delta E}{^a{}_i}})$ and obtain

$$\begin{align} \phi^{GI}(\delta\sigma^c,\delta\tau,\delta\xi^j) &= \phi(\vec{x},t) - \kappa^2 (\delta T(\vec{x},t) - \delta\tau(\vec{x},t)) \pi(\vec{x},t)\nonumber\\ & \quad - \kappa^2 (\delta T^c(\vec{x},t) - \delta\sigma^c(\vec{x},t)) \partial_c \phi(\vec{x},t) \nonumber\\ & \quad +{\mathcal O}^{(2)}_{\phi,\{\delta T\}}(\delta\tau,\delta\sigma^a,\delta\xi^j)+O(\delta^3,\kappa^3), \end{align} \tag{ 2.63 }$$

$$\begin{align} \pi^{GI}(\delta\sigma^c,\delta\tau,\delta\xi^j) &= \pi(\vec{x},t) - \kappa^2 \partial_a [(\delta T(\vec{x},t) - \delta\tau(\vec{x},t)) \partial^a\phi(\vec{x},t)] \nonumber\\ & \quad - \kappa^2 \partial_c[(\delta T^c(\vec{x},t) - \delta\sigma^c(\vec{x},t)) \pi(\vec{x},t)] + \kappa^2 (\delta T(\vec{x},t) - \delta\tau(\vec{x},t))\, m^2 \phi(\vec{x},t)\nonumber \\ & \quad +{\mathcal O}^{(2)}_{\pi,\{\delta T\}}(\delta\tau,\delta\sigma^a,\delta\xi^j)+O(\delta^3,\kappa^3), \end{align} \tag{ 2.64 }$$

$$\begin{align} ({{\delta E}{^a{}_i}})^{GI}(\delta\sigma^a,\delta\tau,\delta\xi^j) &= {{\delta E}{^a{}_i}}(\vec{x},t) - \kappa(\delta_i^a \partial_c -\delta_c^a \delta_i^b \partial_b)(\delta T^c(\vec{x},t) - \delta\sigma^c(\vec{x},t))\nonumber\\ & \quad +\kappa\beta {{\epsilon}{_c{}^{ab}}} \delta_i^c \partial_b(\delta T(\vec{x},t) - \delta \tau(\vec{x},t))\nonumber -\frac{\kappa}{2} {{\epsilon}{_{ij}{}^k}} \delta_k^a (\delta \Xi^j(\vec{x},t)-\delta\xi^j(\vec{x},t))\nonumber \\ & \quad +{\mathcal O}^{(2)}_{{{\delta E}{^a{}_i}},\{\delta T\}}(\delta\tau,\delta\sigma^a,\delta\xi^j)+O(\delta^3,\kappa^2), \end{align} \tag{ 2.65 }$$

$$\begin{align} ({{\delta A}{_a{}^i}})^{GI}(\delta\sigma^a,\delta\tau,\delta\xi^j)&={{\delta A}{_a{}^i}}(\vec{x},t) +\frac{\kappa}{2} \partial_a (\delta\Xi^i(\vec{x},t)-\delta\xi^i(\vec{x},t))\nonumber\\ & \quad +{\mathcal O}^{(2)}_{{{\delta A}{_a{}^i}},\{\delta T\}}(\delta\tau,\delta\sigma^a,\delta\xi^j)+O(\delta^3,\kappa^2), \end{align} \tag{ 2.66 }$$

where we displayed only the linear order in explicit form and $(\delta\tau(\vec{x},t), \delta\sigma^c(\vec{x},t), \delta\xi^j(\vec{x},t))$ denote the spacetime functions corresponding to the reference fields $(\delta T, \delta T^c, \delta\Xi^j)$, see $\delta{\mathcal G},\delta{\mathcal G}^c$ and $\delta{\mathcal G}^j$ in (2.45). The so constructed observables are Dirac observables for all values of the functions $(\delta\tau(\vec{x},t), \delta\sigma^c(\vec{x},t), \delta\xi^j(\vec{x},t))$ and thus each can be understood as a family of Dirac observables parameterised by $\delta\tau(\vec{x},t), \delta\sigma^c(\vec{x},t)$ and $\delta\xi^j(\vec{x},t)$ respectively. Now let us choose specific spacetime functions for $\delta\tau(\vec{x},t), \delta\sigma^c(\vec{x},t)$ and $\delta\xi^j(\vec{x},t)$. A specific choice for $\delta\tau(\vec{x},t), \delta\sigma^c(\vec{x},t)$ and $\delta\xi^j(\vec{x},t)$ can be chosen as follows: For the Gauß constraint we choose vanishing parameters $\delta\xi^j(\vec{x},t)$. The remaining parameters we associate with the temporal and spatial coordinates respectively. Taking into account that the full non-linear parameters τ and σ^c with $\tau = \overline{\tau}+\kappa \delta\tau + O(\kappa^2)$ and $\sigma^c = \overline{\sigma}^c+\kappa\delta\sigma^c + O(\kappa^2)$ should be set to t and $x_\sigma^c$ respectively and that the clocks vanish in the background as well as for higher orders than the linear one, we only need to choose the linearised parameters δτ and $\delta \sigma^c$ yielding

$$\begin{align} \delta\tau(\vec{x},t):= \frac{t}{\kappa},\quad \delta\sigma^c(\vec{x},t):=\frac{x^c_{\sigma}}{\kappa},\quad \delta\xi^j(\vec{x},t):=0, \end{align} \tag{ 2.67 }$$

with $x^c_\sigma$ being the unique solution where $\delta T^c(\vec{x},t) = \frac{1}{\kappa}\sigma^c$ with $\sigma^c = \textrm{const}$. We reinsert these choices into $\delta{\mathcal G},\delta{\mathcal G}^c$ and $\delta{\mathcal G}^j$, respectively, in (2.45). Then we can understand the Dirac observables as functions of $(\vec{x}_\sigma,t)$.

We have that $(\phi^{GI},\pi^{GI})$ are already the independent Dirac observables in the matter sector. However, the 18 observables $({{\delta A}{_a{}^i}})^{GI}, ({{\delta E}{^a{}_i}})^{GI}$ are not all independent, but only four of them are in total. For completeness we have presented the explicit forms of $({{\delta A}{_a{}^i}})^{GI}$ and $({{\delta E}{^a{}_i}})^{GI}$ in Fourier space in (A.45) and (A.46) respectively. We realise that δτ and $\delta\sigma^c$ enter explicitly into these formulae. However, for the choice discussed here where δτ is linearly in t and $\delta\sigma^c$ is constant the contributions both vanish because they come with spatial derivatives acting on $\delta t$ and $\delta\sigma^c$. In case we choose $\delta\sigma^c$ linearly in x^c then $({{\delta E}{^a{}_i}})^{GI}$ has an additional term involving a Kronecker delta that also survives when we express $({{\delta E}{^a{}_i}})^{GI}$ in terms of the independent gauge invariant degrees of freedom similar to what was observed in [74]. As a consequence, considering $({{\delta E}{^a{}_i}})^{GI}$ as an isolated quantity its fall-off behaviour gets modified and does no longer satisfy the requirement in (2.10). This, however causes no issue because we need this quantity as an intermediate step only to reinsert it into the physical Hamiltonian. The final result of the physical Hamiltonian has a suitable fall-off behaviour because these critical terms are combined with scalar field contributions for which we assumed, as usually done, that the initial data has compact support, see also our discussion below (2.76). The four independent degrees of freedom in the gravitational sector are encoded in the symmetric ³ transverse-traceless part of the the variables $({{\delta A}{_a{}^i}})^{GI}, ({{\delta E}{^a{}_i}})^{GI}$ and therefore the 4 independent Dirac observables in the gravitational sector are given by

$$\begin{align} {\delta{\mathcal A}}{_a{}^i}(\vec{x}_\sigma,t)&:= {{P}{_a{}^i{}_j{}^b}}({{\delta{A}}{_b{}^j}})^{GI}(\vec{x}_\sigma,t) \quad {\delta{\mathcal E}}{^a{}_i}(\vec{x}_\sigma,t):={{P}{^a{}_i{}^{\,j}{}_b}}({{\delta{E}}{^b{}_j}})^{GI}(\vec{x}_\sigma,t), \end{align} \tag{ 2.68 }$$

where ${{P}{_a{}^i{}_j{}^b}}$ and ${{P}{^a{}_i{}^{\,j}{}_b}}$ denote the projector on the transverse-traceless part in position space. Its explicit form is given by

$$\begin{align} {{P}{_a{}^i{}_j{}^b}} X(\vec x_\sigma,t)= \frac{1}{2}&\Bigg\{ \left[ \delta_{aj} \delta^{bi} + \delta_a^b \delta_j^i - \delta_a^i \delta_j^b \right] X(\vec x_\sigma,t) + \partial_a \partial^i \partial_j \partial^b \left( X * G^{\Delta \Delta}\right)(\vec x_\sigma,t) \nonumber\\ &\quad +\left[ \delta_{aj} \partial^i \partial^b + \delta^{ib} \partial_a \partial_j + \delta_a^b \partial^i \partial_j + \delta^i_j \partial_a \partial^b - \delta_a^i \partial_j \partial^b - \delta_j^b \partial_a \partial^i \right] \left(X * G^\Delta\right) (\vec x_\sigma,t)\Bigg\}, \nonumber\\ \end{align} \tag{ 2.69 }$$

where $\partial_i : = \delta_i^c \partial_c$, $\partial^i : = \delta^i_c \partial^c$, $\delta_{aj} : = \delta_{ac} \delta^c_j$, and

$$\begin{align} {{P}{^a{}_i{}^{\,j}{}_b}} X(\vec x_\sigma,t)= \delta^{ac}\delta_{ik}\delta^{jm}\delta_{bd}{{P}{_c{}^k{}_m{}^d}} X(\vec x_\sigma,t)\,. \end{align} \tag{ 2.70 }$$

As expected these projectors are non-local in position space but local in momentum space, see the appendix A.2 for more details. A comparison between the local projector in position space also often used in the literature in the context of the transverse-traceless gauge has been analysed in [103].

From now on we will drop the label σ at $\vec{x}_\sigma$ and by an abuse of notation denote it by $\vec{x}$ again and always keep in mind that these x^a coordinates are related to values that the reference field T^a takes. By construction the Dirac observables Poisson commute with all constraints ${\mathcal C}_I^{\prime}$. Since the set ${\mathcal C}_I^{\prime}$ is weakly equivalent to the set ${\mathcal C}_I$ the Dirac observables are weak Dirac observables with respect to ${\mathcal C}_I$ where the weak equivalence refers to the constraint hypersurface defined by the linearised constraints $\delta{\mathcal C}_I$ that, as discussed above, agree with ${\delta C}^{\prime}_I$. Furthermore, as discussed in the appendix A all Dirac observables truncated at linear order, that is neglecting contributions from ${\mathcal O}^{(2)}_{\delta f,\{\delta T\}}$, Poisson commute with all linearised constraints. What still remains to be discussed is the algebra of the independent Dirac observables. For this purpose we will consider the relation between the algebra of Dirac observables and the observable associated to the corresponding Dirac bracket. As proven in theorem 2.2 in [101], see also [66] for an alternative proof in general we have

$$\begin{align*} \{O_{f,\{T\}},O_{g,\{T\}} \}\approx O_{\{f,g\}^*,\{T\}} \end{align*}$$

where $\{f,g\}^*$ denotes the Dirac bracket associated with the total set $({\mathcal G}^I,{\mathcal C}_J)$, that, when requiring ${\mathcal G}^I$ to be gauge fixing conditions is a second class system of constraints for which a Dirac bracket can be constructed. If we have such a relation also at the perturbative level then it provides an efficient way to compute the observable algebra because if we are interested in the algebra up to linear order then we can just expand the observable on the right hand side up to linear order whereas we need to consider the observables up to second order on the left hand side. In the following we will not consider the most general case here but just discuss the necessary result we need for the model derived in this work.

As proven in the appendix under the assumption that we have linearised clocks, that is ${\mathcal G}^I = \delta{\mathcal G}^I$ with $\delta^k{\mathcal G}^I = 0$ for all $k\geqslant 2$ and we consider observables of quantities $f = \delta f$ for which $\delta^k f = 0$ with $k\geqslant 2$ then one can show that

$$\begin{align} \left\{O^{(1)}_{\delta f,\{\delta T\}}+O^{(2)}_{\delta f,\{\delta T\}}\, ,\, O^{(1)}_{\delta g,\{\delta T\}}+O^{(2)}_{\delta g,\{\delta T\}}\right \}=\{\delta f,\delta g\}^*+ O^{(1)}_{\{\delta f,\delta g\}^*,\{\delta T\}} +O(\delta^2,\kappa^2), \end{align} \tag{ 2.71 }$$

where the Dirac bracket $\{\delta f,\delta g\}^*$ corresponds to the set $({\mathcal G}^I = \delta{\mathcal G}^I,{\mathcal C}_I^{\prime} = {\mathcal C}^{^{\prime} (1)}_I+{\mathcal C}^{^{\prime} (2)}_I)$. Again as above because we do not perturb the matter degrees of freedom we can apply the same result here if we replace $\delta f$ by f again with $f = \phi,\pi$. As discussed in appendix A.3.2 the case interesting for us is when $f,g$ are the elementary phase space variables of the linearised theory. Given this and using that for the reference fields we have ${\mathcal G}^I = \delta{\mathcal G}^I$ as well as $\{{\mathcal G}^I(x),{\mathcal C}^{\prime}_J(y)\} = \frac{1}{\kappa}\delta^I_J\delta^{(3)}(x,y)+O(\delta^2,\kappa^2)$ the explicit form of the Dirac bracket up to linear order reads

$$\begin{align} \{\delta f(x),\delta g(y)\}^*:= &\{\delta f(x), \delta g(y)\}-\kappa\int\limits_\sigma d^3z\{\delta f(x), \delta{\mathcal C}^{\prime}_L(z)\}\delta^L_M\{\delta{\mathcal G}^M(z),\delta g(y)\}\nonumber \\ & -\kappa\int\limits_\sigma d^3z\{\delta f(x), \delta^2{\mathcal C}^{\prime}_L(z)\}\delta^L_M\{\delta{\mathcal G}^M(z),\delta g(y)\}\nonumber \\ & +\kappa\int\limits_\sigma d^3z\{\delta g(x), \delta{\mathcal C}^{\prime}_L(z)\}\delta^L_M\{\delta{\mathcal G}^M(z),\delta f(y)\}\nonumber \\ & +\kappa\int\limits_\sigma d^3z\{\delta g(x), \delta^2{\mathcal C}^{\prime}_L(z)\}\delta^L_M\{\delta{\mathcal G}^M(z),\delta f(y)\}+O(\delta^2,\kappa^2).\nonumber \end{align}$$

Now for $\delta f$ and $\delta g$ that have the property that they commute with all $\delta{\mathcal G}^I$ the Dirac bracket reduces to the ordinary Poisson bracket. If $\delta f,\delta g$ are further elementary phase space variables that satisfy standard canonical Poisson brackets then their Poisson bracket is phase space independent and as a consequence, the Dirac observable of the Dirac bracket has only a zeroth order contributions that agrees with the original Poisson bracket. As shown in appendix A.3.1 the variables $(\phi,\pi), ({\delta{\mathcal A}}{_a{}^i},{\delta{\mathcal E}}{^a{}_i})$ all commute with all $\delta {T}^I$ which is equivalent that they commute with all $\delta{\mathcal G}^I$. Thus, for these set of variables the Dirac bracket agrees with the original Poisson bracket. Consequently, we can immediately conclude that the non-vanishing Poisson brackets of the set of Dirac observables $\phi^{GI},\pi^{GI},{\delta{\mathcal A}}{_a{}^i},{\delta{\mathcal E}}{^a{}_i}$ read

$$\begin{align} \{\phi^{GI}(\vec{x},t),\pi^{GI}(\vec{y},t)\}&=\delta(\vec{x}-\vec{y})+O(\delta^2,\kappa^2), \end{align} \tag{ 2.72 }$$

$$\begin{align} \{{\delta{\mathcal A}}{_a{}^i}(\vec{x},t),{{\delta{\mathcal E}}{^b{}_j}}(\vec{y},t)\}&=\frac{\beta}{\kappa} {{P}{_a{}^i{}_j{}^b}} \delta(\vec{x}-\vec{y})+O(\delta^2,\kappa^2).\, \end{align} \tag{ 2.73 }$$

All remaining ones vanish up to corrections of order $O(\delta^2,\kappa^2)$. The matter variables satisfy the same Poisson algebra as their gauge variant counter parts and for the geometric gauge invariant degrees of freedom the Poisson bracket involves on the right hand side as well the expected projector on the transverse and traceless part. Considering this together with the Poisson algebra of the ${\mathcal G}^I$ and ${\mathcal C}_I^{\prime}$ in (2.59) we realise we can perform a canonical transformation from the set of variables

$$\begin{align*} (\phi,\pi), \quad ({\delta{A}}{_a{}^i},{\delta{E}}{^a{}_i}) \end{align*}$$

to the new set of variables

$$\begin{align*} (\delta{\mathcal G},C^{\prime}),\quad (\delta {\mathcal G}^a,C^{\prime}_a), \quad (\delta {\mathcal G}^j,G^{\prime}_j), \quad (\phi^{GI},\pi^{GI}), \quad ({{\delta{\mathcal A}}{_a{}^i}},{{\delta{\mathcal E}}{^a{}_i}}) \end{align*}$$

with

$$\begin{align} \delta{\mathcal G}=\delta T -t,\quad \delta{\mathcal G}^a=\delta T^a-x^a_\sigma,\quad \delta {\mathcal G}^j=\delta\Xi^j. \end{align} \tag{ 2.74 }$$

As discussed in the seminal work in [74], one can reformulate the corresponding forms of the geometric clocks so that they no longer explicitly involve the coordinates. This is possible here because the way the geometric clocks are defined is non-local due to the Green's function involved in the position space representation. The equivalent local conditions are obtained by applying the Laplacian, which has the consequence that neither t nor x^a are explicitly involved any more. It is also these local quantities that are involved in the perturbed quantities, and therefore the use of such geometric clocks does not contradict the requirement that the perturbations be small. The reason is that these coordinate conditions only set certain components of the metric perturbations and their conjugate moments to zero, see also the more detailed discussion in [74] and references therein.

Because the definition of ${\mathcal C}_I^{\prime}$ as can be seen from (2.60) involves second-order contributions to implement this canonical transformation on the linearised phase space we restrict ${\mathcal C}_I^{\prime}$ to its linear part and using that in our case $\delta{\mathcal C}^{\prime}_I = \delta{\mathcal C}_I$ the new set of canonical variables on the linearised phase space is given by

$$\begin{align*} (\delta{\mathcal G},\delta C),\quad (\delta {\mathcal G}^a,\delta C_a), \quad (\delta {\mathcal G}^j,\delta G_j), \quad (\phi^{GI},\pi^{GI}), \quad ({{\delta{\mathcal A}}{_a{}^i}},{{\delta{\mathcal E}}{^a{}_i}})\,, \end{align*}$$

the explicit transformation can be found in appendix A.4. An advantage of this new set of elementary phase space variables is that it allows to clearly separate the physical degrees of freedom from the remaining gauge degrees of freedom at the level of the full phase space. The physical phase space corresponds to the subspace involving the subalgebra of all Dirac observables $(\phi^{GI},\pi^{GI}), ({{\delta{\mathcal A}}{_a{}^i}},{{\delta{\mathcal E}}{^a{}_i}})$ including the expected six physical degrees of freedom. The price to pay in this context is that, as far as the position space is considered, we need to work with highly non-local clocks that become local in momentum space. This is not unexpected here considering the fact that the projector to the transverse-traceless part of the gravitational perturbations is local in momentum but non-local in position space ⁴ .

Finally, let us briefly sketch the strategy we followed here that can also be applied to a different choice of reference fields. Suppose we are interested in perturbation theory up to linear order.

• First choose a set of linearised reference fields, that is $\delta^2{\mathcal G}^I = 0$ such that they build canonical pairs with the linearised constraints, that is $\{\delta\widetilde{\mathcal G}^I(x),\delta{\mathcal C}_J(y)\} = \delta^I_J\delta^{(3)}(x,y)$.
• In case not all reference fields and hence $\delta\widetilde{\mathcal G}^I(x)$ mutually commute then apply the dual observable map with the linearised constraints as the clocks up to second order ⁵ to those that do not to obtain ${\mathcal G}^I$.
• Next apply the dual observable map up to second order neglecting the linear order to all constraints ${\mathcal C}_I$ to get the constraints ${\mathcal C}_I^{\prime}$ that are abelian and canonically conjugate to the reference fields ${\mathcal G}^I$ up to corrections of second order.
• Define the observable map by means of the constraints ${\mathcal C}_I^{\prime}$ using as reference fields ${\mathcal G}^I$.
• Choose next to $(\delta{\mathcal G}^I,\delta{\mathcal C}^{\prime}_I)$ further independent phase space variables on the linearised phase space denoted by $\delta q^I,\delta p_J$.
• Apply the observable observable map to $(\delta q^I,\delta p_I)$ to obtain the physical gauge invariant degrees of freedom.
• Compute the algebra of the Dirac observables of $(\delta q^I,\delta p_I)$ using (2.71).
• If we further apply the dual observable map to these Dirac observables and they still satisfy the same Poisson algebra as their gauge variant counterparts we can also say that given the set of variables $(\delta{\mathcal G}^I,\delta{\mathcal C}^{\prime}_I,\delta q^I,\delta p_I)$ then applying the observable map and its dual to the entire set with the modification that we exclude the linear order in the observable map and its dual for the set $(\delta{\mathcal G}^I,\delta{\mathcal C}^{\prime}_I)$ defines a canonical transformation on the entire phase space.

In our case the final application of the dual observable map on the variables $(\delta q^I,\delta p_J)$ corresponding to $\phi^\textrm{GI},\pi^\textrm{GI},\delta{\mathcal A}^i_a,\delta{\mathcal E}^a_i$ acts trivially since the original gauge variant quantities already commute with our chosen geometrical clocks. We expect that one can extend this strategy without any problems to higher orders in perturbation theory and for the aspect of constructing abelian constraints this has be done in [69]. How and if also the further steps can be generalised to higher order goes beyond the scope of this work and will be discussed more in detail elsewhere.

The work in [35] also aims at constructing Dirac observables in the framework of the ADM formalism. In appendix A.3.4 we compare their construction to the one presented here. To our understanding the explicit form of the observable for the scalar field momentum presented in [35] is not invariant under transformations generated by the linearised spatial diffeomorphism constraint.

We will close this subsection with discussing the dynamics of our constructed Dirac observables. This is the last ingredient missing in order to take the gauge invariant formulation as the starting point for the quantisation in the next section. As we know from full general relativity, the dynamics of Dirac observables cannot be generated by the canonical Hamiltonian because by construction, all Dirac observables commute with the constraints and $\mathbf{H}_\textrm{can}$ is just a linear combination of the smeared constraints. The generator of the dynamics of the Dirac observables is the so-called physical Hamiltonian [101]. If we rewrite the Hamiltonian constraint linearly in the clock momentum of the temporal reference field, that is $C = P_T+h$, then the physical Hamiltonian corresponds to the Dirac observable associated with the phase space function h. Here we need to adapt this to the framework of perturbations theory and consider the fact that in addition we have an interaction Hamiltonian as well as a non-vanishing Hamiltonian in κ⁰ order. The perturbed Hamiltonian constraint up to second order is given by $C = \delta C+\delta^2C^\textrm{geo}+O(\delta^2,\kappa^2)$. The reference field $\delta T$ is canonically conjugate to the linearised Hamiltonian constraint $\delta C$ and thus we can identify $\delta C$ with the momentum variable conjugate to $\delta T$ that we denote by $\delta P_T$, then we have $C = \delta P_T+\delta^2C^\textrm{geo}+O(\delta^2,\kappa^2)$. ⁶ The Dirac observable associated with $\delta^2C^\textrm{geo}$ corresponds to that part of the physical Hamiltonian that is generating the dynamics of the purely gravitational physical degrees of freedom. The Dirac observable corresponding to the interaction Hamiltonian ${\mathcal H}_I = \frac{\kappa}{2}\delta h_{ab}T^{ab}$ with $\delta h_{ab} = -\delta_b^i\delta_{ac} \delta{{E}{^c{}_i}}-\delta_a^i\delta_{bc} \delta{{E}{^c{}_i}}+\delta_{ab} \delta_c^i {{\delta E}{^c{}_i}}$ encodes the interaction between the matter and gravitational physical degrees of freedom. As we will see below when expressed in terms of the independent Dirac observables it decomposes into a term where the physical triad variables $\delta{\mathcal E}{}{^a{}_i}$ couple to the physical matter variables as well as a self-interaction term that involves a coupling between the physical matter variables only. Finally, in κ⁰-order the physical Hamiltonian involves the matter Hamiltonian for the scalar field on Minkowski spacetime. Summarising the above discussion we denote the physical Hamiltonian of the model considered in this work by $\delta\textbf{H}$ that takes the following form

$$\begin{align} \delta\textbf{H} ={\mathcal H}_{\phi} + O_{{\mathcal H}_\textrm{geo},\{T\}}+O_{{\mathcal H}_I,\{T\}}\quad\textrm{with}\quad {\mathcal H}_\textrm{geo}:=\int\limits_\sigma d^3x\; \delta^2C^\textrm{geo}(x),\quad {\mathcal H}_\phi:=\int\limits_\sigma d^3x\; \epsilon(\vec{x},t).\nonumber\\ \end{align} \tag{ 2.75 }$$

As can be seen in (2.76) it can be written entirely in terms of the Dirac observables $(\phi^{GI},\pi^{GI})$, $({{\delta{\mathcal A}}{_a{}^i}},{{\delta{\mathcal E}}{^a{}_i}})$ and therefore is by construction gauge invariant up to corrections of order $O(\delta^2,\kappa^2)$. Here, as discussed above, we insert the Dirac observables in all terms being linearly in κ but keep $\phi,\pi$ in the κ⁰ contribution because at zeroth order $\phi,\pi$ are suitable observables. Because in second-order perturbation theory the constraints ${\mathcal C}^{\prime}_I$ are only weakly equivalent to the original constraints ${\mathcal C}_I$ and the observable of the interaction Hamiltonian ${\mathcal H}_I$ is linearly in the perturbations we need to consider the Dirac observable of ${{\delta E}{^a{}_i}}$ up to second order in ${\mathcal H}_I$ when computing Poisson brackets with the constraints. Hence, as far as the set of original constraints ${\mathcal C}_I$ are considered $\delta \textbf{H}$ is a weak Dirac observable.

Explicitly we obtain, see appendix A.5:

$$\begin{align} \delta\mathbf{H} = &\int_{\mathbb{R}^3} d^3x \;\epsilon(\vec{x},t) +\kappa \int_{\mathbb{R}^3}d^3k\; \frac{1}{2\beta^2}\, \sum_{r\in\{\pm\}}\Big( \underset{\smile}{\delta A}^r(\vec{k},t) \,\underset{\smile}{\delta A}^r (-\vec{k},t)+ 2r\,||\vec{k}|| \underset{\smile}{\delta E}^r(\vec{k},t) \, \underset{\smile}{\delta A}^r(-\vec{k},t) \nonumber\\ & + (\beta^2+1) ||\vec{k}||^2 \underset{\smile}{\delta E}^r(\vec{k},t)\, \underset{\smile}{\delta E}^r(-\vec{k},t)\Big) +\kappa \int_{\mathbb{R}^3} d^3k \;\delta_{b}^i \underset{\smile}{\varphi} {}_a^b(\vec{k},t) \, \underset{\smile}{\delta \mathcal{E}}{}{^a{}_i}(-\vec{k},t)\nonumber\\ &{+\kappa \int d^3k \;ik_b \kappa\underset{\smile}{\delta \sigma}^{^{\prime} a} \left[ \underset{\smile}{\varphi}{}_a^b - \delta_a^b \left( \underset{\smile}{\varphi}{}_c^c +m^2 \underset{\smile}{\widetilde{V}}-\underset{\smile}{\epsilon} \right)\right]+\kappa \int d^3k \left[ i \kappa ||\vec k|| \underset{\smile}{\delta \tau} \hat{k}^b \underset{\smile}{p}^{^{\prime}}{}_b \right]}\nonumber\\ &-\kappa\int_{\mathbb{R}^3} \frac{d^3k}{||\vec{k}||^2} \Bigg\{- \underset{\smile}{p}{{}{_c}}(\vec{k},t) \,\underset{\smile}{p}{}{_d}(-\vec{k},t) \left[ \delta^{cd}-\frac{1}{4} \hat{k}^c \hat{k}^d \right] + \underset{\smile}{\epsilon}(\vec{k},t)\,\underset{\smile}{\epsilon}(-\vec{k},t) \left[\frac{1}{4}+\frac{3}{2}\right] \nonumber\\ &\quad\qquad\qquad\qquad -\frac{3}{2} m^2\,\underset{\smile}{\widetilde{V}}(\vec{k},t) \,\underset{\smile}{\epsilon}(-\vec{k},t)-\underset{\smile}{\varphi} {}_a^a(\vec{k},t) \,\underset{\smile}{\epsilon}(-\vec{k},t) \Bigg\}\,, \end{align} \tag{ 2.76 }$$

with

$$\begin{align} \epsilon(\vec{x},t):=\epsilon(\phi^{GI}(\vec{x},t),\pi^{GI}(\vec{x},t)), \; p_a(\vec{x},t):=p_a(\phi^{GI}(\vec{x},t),\pi^{GI}(\vec{x},t)), \end{align} \tag{ 2.77 }$$

$$\begin{align} \widetilde{V}(\vec{x},t) := \phi^{GI}(\vec{x},t)^2, \, \varphi_a^b(\vec{x},t) := (\partial_a \phi^{GI}(\vec{x},t))(\partial^b \phi^{GI}(\vec{x},t))\,, \end{align} \tag{ 2.78 }$$

where as before we considered terms up to linear order in κ only. Note that an underbow denotes the spatial Fourier transform of the corresponding quantity, i.e.

$$\begin{align*} \underset{\smile}{f}(\vec{k},t):=\frac{1}{(2\pi)^{\frac{3}{2}}}\int_{\mathbb{R}^3}d^3x\; e^{-i\vec{k}\vec{x}} f(\vec{x},t), \qquad \hat{k}^a := \frac{1}{||\vec{k}||} k^a. \end{align*}$$

Note that for the κ⁰ order we also use the operator valued distributions of the physical matter variables because the difference to the gauge variant variables is of higher order than κ⁰ and is neglected to obtain a model that is consistent with the Post-Minkowski approximation order by order in κ. As can be seen from the explicit form of the physical Hamiltonian in (2.76) for our choice of δτ and $\delta\sigma^c$ in (2.67), δτ being linear in t does not contribute to the physical Hamiltonian and $\delta\sigma^c$ is also absent if we chose it to be constant and enters via Kronecker delta in case we choose it to be linearly in x^c in exact agreement to what was found in [74] for the ADM variables. This additional Kronecker delta contribution for the second choice combines with terms including the scalar field and its derivatives respectively and has thus a suitable fall-off behaviour.

The different contributions of the Hamiltonian in the individual lines can be interpreted as follows: The first line corresponds to the Hamiltonian of a scalar field on a Minkowski background. The next two lines encode the Hamiltonian of a gravitational field in vacuum, the third integral denotes the interaction between the gravitational degrees of freedom and the matter field and finally the remaining contributions encode the gravitational self-interaction of the scalar field. Even though the latter does not contain any gravitational degrees of freedom, it only appears due to the coupling between matter and gravity and results from the coupling of the gauge degrees of freedom in the gravitational sector to the scalar field. Once these gauge degrees of freedom are expressed in terms of the independent physical degrees of freedom we end up with this result. The latter contribution was neglected in [69] by using the argument that only the transverse-traceless components need to be considered. While in the vacuum case one can use a gauge fixing that sets all but the transverse-traceless gravitational degrees of freedom to zero, this is no longer a valid choice for a gauge-fixing if we couple the scalar field, because the constraints involve additional contributions from the scalar field and suitable gauge fixing will result in such a contribution of gravitational self-interaction of the scalar field. Hence, to our understanding, it is not justified to neglect this contributions in the order of perturbation theory considered in [69]. This finalises the discussion on the classical model that will be the starting point for the Fock quantisation in the next section. Having identified the physical degrees of freedom in the model, now we are in the situation to separate the total system into a system and an environmental part, as usually done in the framework of decoherence models. We choose the following separation:

$$\begin{align*} \textrm{system:}\,\, \phi^{GI},\, \pi^{GI} \quad\quad \textrm{environment:}\,\, {{\delta{\mathcal A}}{_a{}^i}},\, {\delta {\mathcal E}}{^a{}_i} \end{align*}$$

and hence the gauge invariant matter sector becomes the system and the physical degrees freedom in the gravitational sector provide the environment.

Finally let us comment on the relation of the fully gauge invariant framework to the gauge fixing discussed in the appendix B.1. As discussed for instance in [104] for certain choices of the ${\mathcal G}^I$ one can relate a model formulated in terms of fully gauge invariant quantities to a corresponding gauge-fixed model. The observables constructed here also fall into this class of models, which means that the gauge-fixed Hamiltonian and the gauge invariant physical Hamiltonian formally agree if we replace all gauge invariant quantities by their gauge-fixed counter parts. Practically, this can be achieved by setting ${\mathcal G} = {\mathcal G}^a = {\mathcal G}^j = 0$, that is strongly equal to zero and then we can formally identify the Dirac observables with the gauge-fixed quantities involved in the gauge fixing discussed in the appendix B.1. For this specific gauge fixing, the constructed observables have a standard interpretation. The Hamiltonian of the gauge fixed theory can be obtained by inserting all constraints, gauge fixing conditions as well as the Lagrange multipliers that ensure the stability of the gauge fixing into $\delta \mathbf{H}_\textrm{can}$. As can be seen from the explicit form of $\delta \mathbf{H}_{can}$ in (A.106) the Hamiltonian of the gauge fixed theory agrees formally with $\delta\mathbf{H}$ under the identification of the Dirac observables with their corresponding gauge-fixed quantities. In order to compare our results with the existing literature we discuss in appendix B.2 the geometrical clocks used in [69] that are denoted as ADM clocks there because they were introduced in the seminal paper [74]. Also this set of geometrical clocks encodes the physical gravitational degrees of freedom in the symmetric transverse-traceless linearised connection and triad variables. The main difference we see among the two sets of geometrical clocks is how the reference field associated with the linearised Gauß constraint is chosen. In our work we choose a Lorentz-like condition if we set the corresponding gauge fixing condition $\delta{\mathcal G}^j = 0$, whereas in [69] the resulting $\delta {\mathcal G}^j = 0$ is only equivalent to a Lorentz-like condition in the vacuum case but not if we couple a scalar field. Then the resulting reference fields for the linearised spatial diffeomorphism and Hamiltonian constraint also differ because for both sets of clocks we require that the all geometrical clocks mutually commute. The physical Hamiltonian for the ADM clocks is also presented in (B.2) and it turns out that compared to the physical Hamiltonian in (2.76) it only differs in the term describing the self-interaction of the scalar field. There we find that it differs only by numerical factors in all but one term and it involves one additional term where in our case only a trace enters, but for the ADM clocks an additional contraction with the Fourier basis is present, see (B.18) in appendix B.2.

However, working at the gauge invariant level and not choosing one specific gauge fixing allows us also to choose a different gauge fixing than the one discussed in the appendix B.1 and then the observables and their dynamics discussed here are still valid because everything was formulated in a gauge invariant manner. The only difference for other gauge fixing choices is that the interpretation of the gauge invariant Dirac observables and their relation to the gauge-fixed quantities is modified because then, the gauge fixing does not necessarily correspond to setting all ${\mathcal G}^I$ equal to zero. Hence in this sense, if we consider a class of models where a relation to one gauge fixing is possible, then there is always a convenient choice for the reference fields in order to make such a relation as simple as possible. This was exactly our motivation for choosing the reference fields for this model the way we did, apart from the additional requirements (i) to (iii) listed below (2.44).

The starting point in this top-down approach of the derivation of every quantum master equation is the microscopic dynamics in the form of the Liouville-von Neumann equation, regardless of the exact approach that is taken thereafter. With an interaction Hamiltonian operator $\tilde{H}(t)$ and the density operator of the total system in the interaction picture given by $\tilde{\rho}(t)$, where we denote operators in the interaction picture with a tilde, the evolution equation initially amounts to:

$$\begin{align} \frac{\partial}{\partial t} \tilde{\rho}(t) = -i \alpha [\tilde{H}(t), \tilde{\rho}(t)] =: \alpha \mathcal{L}(t) \tilde{\rho}(t), \end{align} \tag{ 4.1 }$$

where α denotes the dimensionless ⁷ coupling strength between the system and the environment and $\mathcal{L}(t)$ is the so-called Liouville superoperator. The goal is to eliminate the explicit evolution of the many environmental degrees of freedom in order to arrive at an effective equation for the relevant system's degrees of freedom, generally of the following form:

$$\begin{align} \frac{\partial}{\partial t} \tilde{\rho}_S(t) &= -i \alpha [\tilde{H}_{LS}(t), \tilde{\rho}_S(t)] + \tilde{\mathcal{D}}\big( \tilde{\rho}_S(t) \big) \qquad \qquad \qquad \textrm{(interaction picture)} \end{align} \tag{ 4.2 }$$

$$\begin{align} \frac{\partial}{\partial t} \rho_S(t) = -i \alpha [H_S(t) + H_{LS}(t), \rho_S(t)] + \mathcal{D}\big( \rho_S(t) \big) \qquad \; \textrm{(Schr}\ddot{\mathrm{o}}\textrm{dinger picture)}. \end{align} \tag{ 4.3 }$$

In general, there can be a Lamb shift-type contribution $H_{LS}(t)$ in addition to the unitary dynamics generated by $H_S(t)$ in the Schrödinger picture equation, alongside the dissipator $\mathcal{D}\big( \tilde{\rho}_S(t) \big)$, the latter encodes both decoherence and dissipation. Note that the dissipator can in general retain some time dependence even after the transition to the Schrödinger picture, this critically depends on the model and the applied approximations. At this point there are various possibilities to proceed, the easiest and most frequently used path is towards the so-called Lindblad equation. We would like to mention that the Lindblad equation is a very specific form of a master equation in the sense that it originates from a bounded generator of a quantum-dynamical semigroup [45, 46]. Naturally, the Lindblad equation is completely positive and trace-preserving, a feature not all master equations exhibit [13] without further adjustments as we will see an the end of this section. We will start with a fully analytical and in principle exact approach towards the Nakajima-Zwanzig equation, which will provide the basis for the TCL master equation we will later use for the model considered in this work. This relies fundamentally on the results in [80–83]. Hereby we will closely follow the pedagogical outline given in [13] and [12]. The density operator $\tilde{\rho}(t)$ in (4.1) is trace class by construction, hence is contained in the set $\mathcal{B}(\mathcal{H})$ of linear, bounded operators. We proceed to define a pair of projection superoperators called $\mathcal{P}$ and $\mathcal{Q}$:

$$\begin{align} \mathcal{P}: \mathcal{B}(\mathcal{H}) \to \mathcal{B}(\mathcal{H}), \qquad \mathcal{P}(\rho) &:= tr_{\mathcal{E}}(\rho) \otimes \rho_{\mathcal{E}}, \end{align} \tag{ 4.4 }$$

$$\begin{align} \mathcal{Q}: \mathcal{B}(\mathcal{H}) \to \mathcal{B}(\mathcal{H}), \qquad \mathcal{Q}(\rho) &:= \rho - \mathcal{P}(\rho), \end{align} \tag{ 4.5 }$$

where $\rho_{\mathcal{E}}$ is a stationary environmental state that we will specify later and $tr_{\mathcal{E}}$ defines the environmental trace, i.e. the trace over all degrees of freedom that do not belong to the system of interest. The fact that these operators satisfy $\mathcal{P}^2 = \mathcal{P}$, $\mathcal{Q}^2 = \mathcal{Q}$, $\mathcal{P}+\mathcal{Q} = \unicode{x1D7D9}$ and $\mathcal{P}\mathcal{Q} = \mathcal{Q}\mathcal{P} = 0$ and thus are indeed projectors is directly evident from the definitions above. The operators $\mathcal{P}$ and $\mathcal{Q}$ signify the relevant and irrelevant part of the combined dynamics, respectively, in the sense that everything projected with $\mathcal{P}$ exclusively contains the dynamics of the system degrees of freedom, with only an effective influence coming from the traced-out environment. These operators do furthermore commute with the partial time-derivative and hence can be directly applied to the Liouville-von Neumann equation (4.1) to obtain:

$$\begin{align} \mathcal{P} \frac{\partial}{\partial t} \tilde{\rho}(t) &= \frac{\partial}{\partial t} \mathcal{P} \tilde{\rho}(t) = \alpha \mathcal{P} \mathcal{L}(t) \tilde{\rho}(t) = \alpha \big( \mathcal{P} \mathcal{L}(t) \mathcal{P} \tilde{\rho}(t) + \mathcal{P} \mathcal{L}(t) \mathcal{Q} \tilde{\rho}(t) \big), \end{align} \tag{ 4.6 }$$

$$\begin{align} \mathcal{Q} \frac{\partial}{\partial t} \tilde{\rho}(t) &= \frac{\partial}{\partial t} \mathcal{Q} \tilde{\rho}(t) = \alpha \mathcal{Q} \mathcal{L}(t) \tilde{\rho}(t) = \alpha \big( \mathcal{Q} \mathcal{L}(t) \mathcal{P} \tilde{\rho}(t) + \mathcal{Q} \mathcal{L}(t) \mathcal{Q} \tilde{\rho}(t) \big). \end{align} \tag{ 4.7 }$$

Hereby we used that the projection operators are linear. In the last step of each line we used that the two projectors combine to unity. In the end we would like to find a solution for the relevant part $\mathcal{P} \tilde{\rho}(t)$ of the density operator, the projectors do however not commute with the Liouvillian $\mathcal{L}(t)$. As a first step to eliminate the appearance of $\mathcal{Q} \tilde{\rho}(t)$ in the first and hence for us interesting and relevant equation (4.6), we write down a formal solution of equation (4.7):

$$\begin{align} \mathcal{Q} \tilde{\rho}(t) = \mathcal{G}(t,t_0) \mathcal{Q} \tilde{\rho}(t_0) + \alpha \int^t_{t_0} ds \, \mathcal{G}(t,s) \mathcal{Q} \mathcal{L}(s) \mathcal{P} \tilde{\rho}(s), \quad \mathcal{G}(t,s) := \mathcal{T}_{\leftarrow} \exp \Big[ \int^t_s d\tau \, \mathcal{Q} \mathcal{L}(\tau) \Big], \end{align} \tag{ 4.8 }$$

where $\mathcal{T}_{\leftarrow}$ stands for time-ordering with the largest argument to the very left. Applying a partial derivative with respect to t yields equation (4.7) while the initial condition for $t = t_0$ is satisfied due to $\mathcal{G}(t_0,t_0) = \unicode{x1D7D9}$ because of a vanishing integral, so the stated solution is unique. Inserting (4.8) into the relevant equation (4.6) gives:

$$\begin{align} \frac{\partial}{\partial t} \mathcal{P} \tilde{\rho}(t) = \alpha \mathcal{P} \mathcal{L}(t) \mathcal{P} \tilde{\rho}(t) + \alpha \mathcal{P} \mathcal{L}(t) \mathcal{G}(t,t_0) \mathcal{Q} \tilde{\rho}(t_0) + \alpha^2 \int^t_{t_0} ds \, \mathcal{P} \mathcal{L}(t) \mathcal{G}(t,s) \mathcal{Q} \mathcal{L}(s) \mathcal{P} \tilde{\rho}(s). \end{align} \tag{ 4.9 }$$

This is the Nakajima-Zwanzig equation, albeit admittedly not less complicated than the attempt of solving the initial problem, mostly due to the appearance of the integral kernel

$$\begin{align} \mathcal{K}(t,s) := \mathcal{P} \mathcal{L}(t) \mathcal{G}(t,s) \mathcal{Q} \mathcal{L}(s) \mathcal{P}. \end{align} \tag{ 4.10 }$$

It is worth noting that this manifestly non-Markovian equation (4.9) is in principle exact, there was no need for approximations whatsoever until this point. If we assume that the interaction Hamiltonian density is of the general form $\sum_\alpha S_\alpha \otimes E_\alpha$ with a discrete summation index α (omitting model-specific time- or mode dependencies for the sake of brevity) and E_α consisting of linear contributions of annihilation- and creation-type variables, the first contribution in equation (4.9) vanishes. This is simply due to the fact that odd numbers of the environmental interaction term do not contribute in the environmental trace, given a suitable initial state, specifically a Gibbs state. This observation is contained in the very definition (4.4) of the relevant projector $\mathcal{P}$, which turns out to be useful even beyond the mere simplification of equation (4.9) as we will see later in this section. For this broad class of models, a method called cumulant expansion turns out to be a powerful tool for a step-by-step approximation of a time-convolutionless approach to non-Markovian dynamics, first developed in [105] and [106], respectively. Up to now the main complexity of the Nakajima-Zwanzig equation remains, the convolution integral makes it nearly impossible to be solved analytically. To this extent, we introduce the two propagators $F(t,s)$ and $G(t,s)$ such that $G(t,s)\tilde{\rho}(s) = \tilde{\rho}(t)$ and $F(t,s) \circ G(t,s) = \unicode{x1D7D9}$ with $s \leqslant t \in \mathbb{R}$. Note the distinction $\mathcal{G}(t,s) \neq G(t,s)$ in the previously used notation, the latter being the propagator of the total system and not merely the one of the irrelevant projection. Hence:

$$\begin{align} \frac{\partial}{\partial t} \Big( F(t,s) \circ G(t,s) \tilde{\rho}(s) \Big) = \Big( \frac{\partial}{\partial t} F(t,s) \Big) \tilde{\rho}(t) + \alpha F(t,s) \mathcal{L}(t) \tilde{\rho}(t) = 0. \end{align} \tag{ 4.11 }$$

This can be easily solved formally for the inverse propagator:

$$\begin{align} \frac{\partial}{\partial t} F(t,s) = - \alpha F(t,s) \mathcal{L}(t) \quad \Longleftrightarrow \quad F(t,s) = \mathcal{T}_{\rightarrow} \exp \Big[- \alpha \int^t_s d\tau \, \mathcal{L}(\tau) \Big], \end{align} \tag{ 4.12 }$$

where $\mathcal{T}_{\rightarrow}$ denotes anti-time-ordering, that is with the largest argument to the very right. As the next step, we insert $\tilde{\rho}(s) = F(t,s) \big( \mathcal{P} + \mathcal{Q} \big) \tilde{\rho}(t)$ into the solution of equation (4.7) and obtain:

$$\begin{align} \mathcal{Q} \tilde{\rho}(t) &= \mathcal{G}(t,t_0) \mathcal{Q} \tilde{\rho}(t_0) + \alpha \int_{t_0}^t ds \, \mathcal{G}(t,s) \mathcal{Q} \mathcal{L}(s) \mathcal{P} F(t,s)(\mathcal{P} + \mathcal{Q}) \tilde{\rho}(t)\nonumber \\ &=: \mathcal{G}(t,t_0) \mathcal{Q} \tilde{\rho}(t_0) + \Sigma(t,t_0)(\mathcal{P} + \mathcal{Q}) \tilde{\rho}(t), \end{align} \tag{ 4.13 }$$

where we collected the convolution part into a separate, now central entity $\Sigma(t,t_0)$:

$$\begin{align} \Sigma(t,t_0) := \alpha \int_{t_0}^t ds \, \mathcal{G}(t,s) \mathcal{Q} \mathcal{L}(s) \mathcal{P} F(t,s). \end{align} \tag{ 4.14 }$$

Now our aim is to solve (4.13) for $\mathcal{Q} \tilde{\rho}(t)$ by collecting all terms on one side:

$$\begin{align} \big( \unicode{x1D7D9} - \Sigma(t,t_0) \big) \mathcal{Q} \tilde{\rho}(t) = \mathcal{G}(t,t_0) \mathcal{Q} \tilde{\rho}(t_0) + \Sigma(t,t_0) \mathcal{P} \tilde{\rho}(t). \end{align} \tag{ 4.15 }$$

The key point is the treatment of $\big( \unicode{x1D7D9} - \Sigma(t,t_0) \big)$ or rather the inverse of it. Evidently, the initial conditions $\Sigma(t_0, t_0) = 0$ and $\Sigma(t,t_0)|_{\alpha = 0} = 0$ are easy to check, our approach involves the representation of $\big( \unicode{x1D7D9} - \Sigma(t,t_0) \big)^{-1}$ as the limit of a geometric series of operators, that is:

$$\begin{align} \big( \unicode{x1D7D9} - \Sigma(t,t_0) \big)^{-1} = \sum\limits_{n=0}^\infty \big( \Sigma(t,t_0)\big)^n. \end{align} \tag{ 4.16 }$$

It is indeed the case that there exists some neighbourhood around t₀ where the expression in (4.16) is invertible. These steps require some additional mathematical justification, for a more detailed discussion on this point we refer the reader to appendix C in the appendix. Hence, under the assumption that we restrict to an interval $[t_0,t]$ where invertibility is guaranteed:

$$\begin{align} \mathcal{Q} \tilde{\rho}(t) = \big( \unicode{x1D7D9} - \Sigma(t,t_0) \big)^{-1} \mathcal{G}(t,t_0) \mathcal{Q} \tilde{\rho}(t_0) + \big( \unicode{x1D7D9} - \Sigma(t,t_0) \big)^{-1} \Sigma(t,t_0) \mathcal{P} \tilde{\rho}(t). \end{align} \tag{ 4.17 }$$

Finally inserting this into the Nakajima-Zwanzig equation (4.9) for $\mathcal{P}\tilde{\rho}(t)$ yields:

$$\begin{align} \frac{\partial}{\partial t} \mathcal{P} \tilde{\rho}(t) &= \alpha \mathcal{P} \mathcal{L}(t) \Big( \unicode{x1D7D9} + \big( \unicode{x1D7D9} - \Sigma(t,t_0) \big)^{-1} \Sigma(t,t_0) \Big) \mathcal{P} \tilde{\rho}(t) + \alpha \mathcal{P} \mathcal{L}(t) \big( \unicode{x1D7D9} - \Sigma(t,t_0) \big)^{-1} \mathcal{G}(t,t_0) \mathcal{Q} \tilde{\rho}(t_0) \nonumber \\[0.4em] &=: K(t,t_0) \mathcal{P} \tilde{\rho}(t) + \mathcal{I}(t,t_0) \mathcal{Q} \tilde{\rho}(t_0). \end{align} \tag{ 4.18 }$$

This is the so-called TCL master equation. Note that this master equation, although time-local, still conveys non-Markovian processes in contributions with higher than second order in α. In principle it is even exact in the neighbourhood where $\big( \unicode{x1D7D9} - \Sigma(t,t_0) \big)^{-1}$ is invertible, considering terms of all orders. We will not further investigate the inhomogeneity $\mathcal{I}(t,t_0) \mathcal{Q} \tilde{\rho}(t_0)$ at this point since we are focused on separable initial conditions for now, for which this expression vanishes. To this extent, consider the explicit form of the TCL generator $K(t,t_0)$ in terms of the geometric series expansion of the inverse of $\big( \unicode{x1D7D9} - \Sigma(t,t_0) \big)$, neglecting $\mathcal{P} \mathcal{L}(t) \mathcal{P}$ based on the assumption that the environmental part of the interaction Hamiltonian density is linearly comprised of creation- and annihilation-type operators:

$$\begin{align} K(t,t_0) &= \alpha \mathcal{P} \mathcal{L}(t) \big( \unicode{x1D7D9} - \Sigma(t,t_0) \big)^{-1} \Sigma(t,t_0) \mathcal{P} = \alpha \mathcal{P} \mathcal{L}(t) \sum\limits_{n=1}^\infty \big( \Sigma(t,t_0)\big)^n \mathcal{P} \nonumber \\ &= \alpha \mathcal{P} \mathcal{L}(t) \sum\limits_{n=1}^\infty \left( \sum\limits_{m=1}^\infty \alpha^m \Sigma_m(t,t_0) \right)^n \mathcal{P}, \end{align} \tag{ 4.19 }$$

where we noticed that it is apparent that $\Sigma(t,t_0)$, given its explicit form in equation (4.14), can be expanded in powers of α as well. Since we are interested in a weak coupling scenario, it is sufficient to consider contributions up to second order in α, that is only the lowest-order terms that occur in the series expansion of $\Sigma(t,t_0)$. In our particular case this entails only the second-order contributions since $\mathcal{P} \mathcal{L}(t) \mathcal{P} = 0$. This fact amounts to the Born approximation on the one hand, by virtue of the assumption of separable initial conditions with a stationary thermal environment and to the (first) Markov approximation on the other hand. The latter fact becomes evident once we realise that by eliminating the convolution integral with the density operator altogether without correcting for it with higher-order terms, we end up with an equation that has simply replaced $\tilde{\rho}(s)$ with $\tilde{\rho}(t)$. The second Markov approximation is commonly depicted as the limit $t_0 \to -\infty$, based on the assumption that the characteristic time scales of system and environment are vastly different, inducing a peakedness of the correlation functions. The latter arise naturally as functions of the form $\mathrm{tr}_\varepsilon \big(E_\alpha(t) E_\beta(s) \rho_\varepsilon \big)$ in the explicit computation of the environmental trace with $E_\alpha(t)$ meaning the environmental part of the interaction Hamiltonian. This peakedness property means that an extension of the integration interval has a negligible effect on the effective dynamics of the system, which however has to be checked explicitly for any given model. In the specific decoherence model we consider, we do not take this limit yet, all of the previous assumptions are nevertheless justified based on the weak coupling of matter to gravity. The only relevant contribution is the one where $\mathcal{G}(t,s)$ and $F(t,s)$ do not contribute, since everything beyond the identity contains the coupling constant.

$$\begin{align} K(t,t_0) = \alpha^2 \mathcal{P} \mathcal{L}(t) \Sigma_1(t,t_0) \mathcal{P} =\alpha^2 \mathcal{P} \mathcal{L}(t) \int_{t_0}^t ds \, \mathcal{Q} \mathcal{L}(s) \mathcal{P} = \alpha^2 \mathcal{P} \int_{t_0}^t ds \, \mathcal{L}(t) \mathcal{L}(s) \mathcal{P}. \end{align} \tag{ 4.20 }$$

Note that higher-order contributions in the coupling constant get increasingly difficult to evaluate explicitly, based on the complexity of (4.19), see for example [107] for the fourth-order TCL generator. The cumulant expansion procedure [105, 106] provides a systematic framework in which the individual contributions in $K(t,t_0)$ can be evaluated, especially for models in which all terms involving $\mathcal{P} \mathcal{L}(t) \mathcal{P}$ vanish due to the explicit form of the interaction. Now we can explicitly write down the desired master equation we will use for the model considered in this article later:

$$\begin{align} \frac{\partial}{\partial t} \mathcal{P} \tilde{\rho}(t) &= \alpha^2 \mathcal{P} \int_{t_0}^t ds \, \mathcal{L}(t) \mathcal{L}(s) \tilde{\rho}(t) \nonumber \\ &= -\alpha^2 \int_{t_0}^t ds \, tr_{\mathcal{E}} \left\{\Big[\tilde{H}(t), \Big[\tilde{H}(s), tr_\mathcal{E} \left\{\tilde{\rho}(t)\right\} \otimes \tilde{\rho}_\mathcal{E} \Big] \Big] \right\} \otimes \tilde{\rho}_\mathcal{E}. \end{align} \tag{ 4.21 }$$

This second-order TCL master equation could have been obtained in a straightforward fashion: First integrate the von Neumann equation (4.1) and subsequently insert the formal solution into the original equation, akin to the second-order truncation of the Dyson series. Secondly, assume separable initial conditions and replace $\tilde{\rho}_S(s) \to \tilde{\rho}_S(t)$, eliminating the convolution integral. It is however not clear from this perspective how non-Markovian corrections or initial conditions with nonvanishing entanglement entropy can be dealt with. The presented projection operator technique in conjunction with the TCL formalism answers these questions order by order and simultaneously provides a differential (as opposed to an integro-differential) equation for $\tilde{\rho}(t)$. An alternative but largely equivalent route towards a time-local master equation is presented in the next section. Lastly we would like to add a few sentences regarding the relation between the TCL master equation (4.21) and the frequently used Lindblad equation [45, 46]. Trace preservation and the notion of positivity are of paramount importance when it comes to solutions of the master equation. However, not every master equation is guaranteed to admit solutions that adhere to these principles. The Lindblad form ensures that these criteria are automatically met:

$$\begin{align} \frac{\partial}{\partial t} \rho_S(t) = -i \alpha [H_S(t) + H_{LS}(t), \rho_S(t)] + \sum\limits_{\alpha \beta} \lambda_{\alpha \beta} \Big( L_\alpha \rho_S(t) L^\dagger_\beta - \frac{1}{2} \{L^\dagger_\beta L_\alpha, \rho_S(t) \} \Big), \end{align} \tag{ 4.22 }$$

with a finite set of traceless Lindblad operators L_α and positive, semi-definite coefficients $\lambda_{\alpha \beta}$. This equation of so-called first standard form can be diagonalised by a unitary transformation of the operators L_α , after which it is commonly called the second standard form. It is immediately clear that not every equation of the form (4.21) exhibits Lindblad form, firstly because of the remaining time dependence that has to be removed by hand based on the properties of the environmental correlation functions. Secondly, the operator structure also differs in a way that positivity is not automatically ensured, even if the derivation has been performed in a top-down manner [108]. There is no universal pathway that remedies these issues, for certain models it turns out to be feasible to expand the dissipator to obtain Lindblad form [12], whereas it can be sufficient to neglect certain contributions in the course of a rotating-wave approximation [109].

A different albeit mostly equivalent approach of formulating a master equation is the influence phase functional approach, which has the same starting point as the previously introduced projection operator formalism. During the course of this brief review we will follow the idea first presented in [79] while we also follow the pedagogical presentation in [12]. Regarding the notation in terms of propagators and thermal Wightman functions, we work with the same notation as for instance used in [110] or [10]. As before, we start at the von Neumann equation in the interaction picture:

$$\begin{align} \frac{\partial}{\partial t} \tilde{\rho}(t) = -i \alpha [\tilde{H}(t), \tilde{\rho}(t)] =: \alpha \mathcal{L}(t) \tilde{\rho}(t), \quad \tilde{\rho}_S(t) = tr_{\mathcal{E}} \left\{\mathcal{T}_\leftarrow \exp \left\{\alpha \int_{t_0}^t \int_{\mathbb{R}^3} d\tau \, d^3x \, \ell(x) \right\} \tilde{\rho}(t_0) \right\},\nonumber\\ \end{align} \tag{ 4.23 }$$

where $\ell(x) = \ell(\tau, \vec{x})$ is the superoperator density associated to $\mathcal{L}(\tau)$ and $\mathcal{T}_\leftarrow : = \mathcal{T}^{(S)}_\leftarrow \mathcal{T}^{(\varepsilon)}_\leftarrow$ corresponds to chronological time-ordering of system an environmental operators, respectively, with larger time arguments to the left of the expression. Contrary to the Nakajima-Zwanzig formalism, the strategy is to directly eliminate the environmental time-ordering at the level of the iterated Dyson series (4.23). Note that system's operators can be considered as commuting under the time-ordering $\mathcal{T}^{(S)}_\leftarrow$. We assume an interaction Hamiltonian density of the general form $\mathcal{H}(x) = \sum\nolimits_\alpha S_\alpha(x) \otimes E_\alpha(x)$, the indices $\alpha, \beta$ can be understood as (a mixture of either co- or contravariant spacetime- or internal) multi-indices depending on the specific model. The goal of this section is to reformulate equation (4.23) into a more tractable form similar to

$$\begin{align} \tilde{\rho}_S(t) = \exp \{i \Phi_t[S^+, S^-] \} \tilde{\rho}_S(t_0), \end{align} \tag{ 4.24 }$$

where $\Phi_t[S^+, S^-]$ is called the influence phase functional and the $S^\pm$ are superoperators associated to the system of interest. A convenient decomposition of the time ordering in equation (4.23) can be straightforwardly shown by induction for almost commuting environmental operators, i.e. the case that $[E_\alpha(x), E_\beta(x^{^{\prime}})] \propto ic_{\alpha \beta}(x, x^{^{\prime}}) \unicode{x1D7D9}_\varepsilon$ with $c_{\alpha \beta}(x, x^{^{\prime}}) \in \mathbb{R}$. This is done by introducing an ansatz where the time-ordering amounts to a phase factor times the ordinary, non-time-ordered exponential of the Liouvillian, leading to a differential equation which can be directly integrated. A very similar procedure has been performed at the level of the time-evolution operator in [111] for the spin-boson model. This can be carried over almost immediately by recalling the adjoint action of these operators on the density operator and the fact that the $S_\alpha(t)$ can be ordered arbitrarily inside the time-ordering $\mathcal{T}^{(S)}_\leftarrow$. This is the setup this work will rely on, however in the more general case of strictly non-commuting environmental operators, the proof still holds but relies on the application of the finite-temperature Wick theorem [12, 112]. The result is:

$$\begin{align} \tilde{\rho}_S(t) &= \mathcal{T}^{(S)}_\leftarrow tr_{\mathcal{E}} \left\{\mathcal{T}^{(\varepsilon)}_\leftarrow \exp \left\{\alpha \int_{t_0}^t \, \int_{\mathbb{R}^3} d\tau \, d^3x \, \ell(x) \right\} \tilde{\rho}(t_0) \right\} \nonumber \\ &= \mathcal{T}^{(S)}_\leftarrow \exp \left\{\frac{\alpha^2}{2} \int_{t_0}^t \int_{\mathbb{R}^3} d\tau d^3x \int_{t_0}^\tau \int_{\mathbb{R}^3} d\tau^{^{\prime}} d^3x^{^{\prime}} \, \left[ \ell(x), \ell(x^{^{\prime}}) \right] \right\} tr_{\mathcal{E}} \nonumber \\ & \quad \times \left\{\exp \left\{\alpha \int_{t_0}^t \, \int_{\mathbb{R}^3} d\tau \, d^3x \, \ell(x) \right\} \tilde{\rho}(t_0) \right\} \nonumber\\ &= \mathcal{T}^{(S)}_\leftarrow \exp \left\{\frac{\alpha^2}{2} \int_{t_0}^t \int_{\mathbb{R}^3} d\tau d^3x \int_{t_0}^\tau \int_{\mathbb{R}^3} d\tau^{^{\prime}} d^3x^{^{\prime}} \, \left[ \ell(x), \ell(x^{^{\prime}}) \right] \right\} \nonumber\\ & \quad \times \exp \left\{\frac{\alpha^2}{2} \int_{t_0}^t \int_{\mathbb{R}^3} d\tau d^3x \int_{t_0}^t \int_{\mathbb{R}^3} d\tau^{^{\prime}} d^3x^{^{\prime}} \, \left\langle \ell(x) \ell(x^{^{\prime}}) \right\rangle_\varepsilon \right\} \tilde{\rho}_S(t_0), \end{align} \tag{ 4.25 }$$

with the environmental expectation value $\big\langle \ell(x) \ell(x^{^{\prime}}) \big\rangle_\varepsilon : = tr_\varepsilon \big\{\ell(x) \ell(x^{^{\prime}}) \rho_\varepsilon \big\}$ under the assumption that the initial state $\tilde{\rho}(t_0) = \tilde{\rho}_S(t_0) \otimes \tilde{\rho}_\varepsilon(t_0)$ is uncorrelated and hence separable. Note that the integrals in the first exponent are nested, i.e. the commutator of the Liouvillians is explicitly time-ordered whereas the second exponent holds integrals with identical limits. Furthermore we would like to point out that the last line of (4.25) only acts upon $\tilde{\rho}_S(t_0)$. This is due to the fact that for an uncorrelated thermal initial state, the action of the exponentiated Liouvillian can be directly computed in terms of a cumulant expansion with a vanishing first order, terminating after the second order, assuming that $\tilde{\rho}_\varepsilon$ describes a Gaussian state with respect to the environmental interaction variables. Hence it is possible to absorb the expectation value $\big\langle \ell(x) \ell(x^{^{\prime}}) \big\rangle_\varepsilon$ into the exponential. Use of the Jacobi identity and the above definitions lets us further simplify (4.25):

$$\begin{align*} \left[ \ell(x), \ell(x^{^{\prime}}) \right] \tilde{\rho} &= - \left[ \left[\mathcal{H}(x), \mathcal{H}(x^{^{\prime}}) \right], \tilde{\rho} \right] = - \sum\limits_\alpha \sum\limits_\beta \left[ \left[S_\alpha(x) \otimes E_\alpha(x), S_\beta(x^{^{\prime}}) \otimes E_\beta(x^{^{\prime}}) \right], \tilde{\rho}_S \otimes \tilde{\rho}_\varepsilon \right] \\ &\cong - \sum\limits_\alpha \sum\limits_\beta \left[ E_\alpha(x), E_\beta(x^{^{\prime}}) \right] \left[ S_\alpha(x) S_\beta(x^{^{\prime}}) \otimes \unicode{x1D7D9}_\varepsilon, \tilde{\rho}_S \otimes \tilde{\rho}_\varepsilon \right]. \end{align*}$$

Note that the last equality holds under consideration of the overall system time-ordering to the very left of expression (4.25) and under the assumption that the commutator $\big[ E_\alpha(x), E_\beta(x^{^{\prime}}) \big]$ is purely a phase factor, as it is the case for the model in this work. We have indicated that a certain equality only holds under additional time-ordering $\mathcal{T}^{(S)}_\leftarrow$ with $\cong$. This justifies acting upon $\tilde{\rho}_S$ instead of $\tilde{\rho}_S \otimes \tilde{\rho}_\varepsilon$ in (4.25), where the environmental part has been completely reduced to the expectation value in the second exponential. The action of these operator densities on the density operator is to be understood in a formal sense, all quantities are spatially integrated once we insert these identities back into the phase functional. In order to represent the action of the exponentiated commutators, let us introduce the following notation:

$$\begin{align*} S_\alpha^+(x) \tilde{\rho}_S := S_\alpha(x) \tilde{\rho}_S, \qquad S_\alpha^-(x) \tilde{\rho}_S := \tilde{\rho}_S S_\alpha(x), \end{align*}$$

where the commutator of Liouvillians $\big[ \ell(x), \ell(x^{^{\prime}}) \big]$ immediately assumes the form

$$\begin{align*} \big[ \ell(x), \ell(x^{^{\prime}}) \big] = - \sum\limits_\alpha \sum\limits_\beta \Big( \big[ E_\alpha(x), E_\beta(x^{^{\prime}}) \big] S_\alpha^+(x) S_\beta^+(x^{^{\prime}}) + \big[ E_\alpha(x), E_\beta(x^{^{\prime}}) \big] S_\alpha^-(x) S_\beta^-(x^{^{\prime}}) \Big). \end{align*}$$

This concludes the first part of the derivation of the influence phase functional based on (4.25):

$$\begin{align} \tilde{\rho}_S(t) &= \mathcal{T}^{(S)}_\leftarrow \exp \left\{- \frac{\alpha^2}{2} \int_{t_0}^t \int_{\mathbb{R}^3} d\tau d^3x \int_{t_0}^t \int_{\mathbb{R}^3} d\tau^{^{\prime}} d^3x^{^{\prime}} \, \theta(\tau - \tau^{^{\prime}}) \sum\limits_\alpha \sum\limits_\beta \left[ E_\alpha(x), E_\beta(x^{^{\prime}}) \right] S_\alpha^+(x) S_\beta^+(x^{^{\prime}}) \right\} \nonumber \\ & \quad \times \exp \left\{\frac{\alpha^2}{2} \int_{t_0}^t \int_{\mathbb{R}^3} d\tau d^3x \int_{t_0}^t \int_{\mathbb{R}^3} d\tau^{^{\prime}} d^3x^{^{\prime}} \, \theta(\tau - \tau^{^{\prime}}) \sum\limits_\alpha \sum\limits_\beta \left[ E_\alpha(x), E_\beta(x^{^{\prime}}) \right] S_\alpha^-(x) S_\beta^-(x^{^{\prime}}) \right\} \nonumber \\ & \quad \times \exp \left\{\frac{\alpha^2}{2} \int_{t_0}^t \int_{\mathbb{R}^3} d\tau d^3x \int_{t_0}^t \int_{\mathbb{R}^3} d\tau^{^{\prime}} d^3x^{^{\prime}} \, \left\langle \ell(x) \ell(x^{^{\prime}}) \right\rangle_\varepsilon \right\} \tilde{\rho}_S(t_0). \end{align} \tag{ 4.26 }$$

Next we have a closer look at the thermal expectation value $\big\langle \ell(x) \ell(x^{^{\prime}}) \big\rangle_\varepsilon$ and apply its definition:

$$\begin{align*} & \int_{t_0}^t \int_{\mathbb{R}^3} d\tau d^3x \int_{t_0}^t \int_{\mathbb{R}^3} d\tau^{^{\prime}} d^3x^{^{\prime}} \, \big\langle \ell(x) \ell(x^{^{\prime}}) \big\rangle_\varepsilon \, \tilde{\rho}_S \nonumber \\ & \quad = \int_{t_0}^t \int_{\mathbb{R}^3} d\tau d^3x \int_{t_0}^t \int_{\mathbb{R}^3} d\tau^{^{\prime}} d^3x^{^{\prime}} \, tr_\varepsilon\Big(- \big[\mathcal{H}(x), \big[\mathcal{H}(x), \tilde{\rho}_S \otimes \tilde{\rho}_\varepsilon \big] \big] \Big)\nonumber \\ & \quad \cong \int_{t_0}^t \int_{\mathbb{R}^3} d\tau d^3x \int_{t_0}^t \int_{\mathbb{R}^3} d\tau^{^{\prime}} d^3x^{^{\prime}} \, \sum\limits_\alpha \sum\limits_\beta \Big( - \big\langle E_\beta(x^{^{\prime}}) E_\alpha(x) \big\rangle_\varepsilon S_\alpha^+(x) S_\beta^+(x^{^{\prime}})\nonumber \\ & \qquad - \big\langle E_\alpha(x) E_\beta(x^{^{\prime}}) \big\rangle_\varepsilon S_\alpha^-(x) S_\beta^-(x^{^{\prime}}) \!+\! \big\langle E_\beta(x^{^{\prime}}) E_\alpha(x) \big\rangle_\varepsilon S_\alpha^+(x) S_\beta^-(x^{^{\prime}})\! +\! \big\langle E_\alpha(x) E_\beta(x^{^{\prime}}) \big\rangle_\varepsilon S_\alpha^-(x) S_\beta^+(x^{^{\prime}}) \Big), \end{align*}$$

where it is important that the last equality again only holds under the time-ordering still present in (4.26). This is because we exchanged the operators $S_\alpha^+(x), S_\beta^+(x^{^{\prime}})$ and switched the integration variables and summation indices respectively. Afterwards we ended up with what seem to be switched arguments in the correlation function $\big\langle E_\beta(x^{^{\prime}}) E_\alpha(x) \big\rangle_\varepsilon$, likewise for the contribution from $S_\alpha^-(x) S_\beta^+(x^{^{\prime}})$. This gives us the possibility to group the individual terms in the influence functional (4.25) such that we can express it in terms of thermal Wightman functions.

$$\begin{align} \tilde{\rho}_S(t) &= \mathcal{T}^{(S)}_\leftarrow \exp \bigg\{\frac{\alpha^2}{2} \int_{t_0}^t \int_{\mathbb{R}^3} d\tau d^3x \int_{t_0}^t \int_{\mathbb{R}^3} d\tau^{^{\prime}} d^3x^{^{\prime}} \sum\limits_\alpha \sum\limits_\beta \nonumber \\ &\quad\times \bigg[ - \Big( \theta(\tau - \tau^{^{\prime}}) \big[ E_\alpha(x), E_\beta(x^{^{\prime}}) \big] + \big\langle E_\beta(x^{^{\prime}}) E_\alpha(x) \big\rangle_\varepsilon \Big) S_\alpha^+(x) S_\beta^+(x^{^{\prime}}) \nonumber \\ &\quad+ \Big( \theta(\tau - \tau^{^{\prime}}) \big[ E_\alpha(x), E_\beta(x^{^{\prime}}) \big] - \big\langle E_\alpha(x) E_\beta(x^{^{\prime}}) \big\rangle_\varepsilon \Big) S_\alpha^-(x) S_\beta^-(x^{^{\prime}}) \nonumber \\ &\quad+ \big\langle E_\beta(x^{^{\prime}}) E_\alpha(x) \big\rangle_\varepsilon S_\alpha^+(x) S_\beta^-(x^{^{\prime}}) + \big\langle E_\alpha(x) E_\beta(x^{^{\prime}}) \big\rangle_\varepsilon S_\alpha^-(x) S_\beta^+(x^{^{\prime}})\bigg] \bigg\} \tilde{\rho}_S(t_0). \end{align} \tag{ 4.27 }$$

Upon closer inspection, these combinations of correlation functions amount to the propagators:

$$\begin{align} G^{++}_{\alpha \beta}(x-x^{^{\prime}}) &:= \big\langle \mathcal{T}^{(\varepsilon)}_\leftarrow \big( E_\alpha(x) E_\beta(x^{^{\prime}}) \big) \big\rangle_\varepsilon = \theta(\tau - \tau^{^{\prime}}) \big[ E_\alpha(x), E_\beta(x^{^{\prime}}) \big] + \big\langle E_\beta(x^{^{\prime}}) E_\alpha(x) \big\rangle_\varepsilon \nonumber \\[0.3em] &= \big\langle E_\alpha(x) E_\beta(x^{^{\prime}}) \big\rangle_\varepsilon \theta(\tau - \tau^{^{\prime}}) + \big\langle E_\beta(x^{^{\prime}}) E_\alpha(x) \big\rangle_\varepsilon \theta(\tau^{^{\prime}} - \tau) \nonumber \\[0.3em] &=: G_{\alpha \beta}^{>}(x-x^{^{\prime}}) \theta(\tau - \tau^{^{\prime}}) + G_{\alpha \beta}^{\lt}(x-x^{^{\prime}}) \theta(\tau^{^{\prime}} - \tau), \end{align} \tag{ 4.28 }$$

and similarly although anti-time-ordered for the other contribution:

$$\begin{align} -G^{-}_{\alpha \beta}(x-x^{^{\prime}}) &:= - \big\langle \mathcal{T}^{(\varepsilon)}_\rightarrow \big( E_\alpha(x) E_\beta(x^{^{\prime}}) \big) \big\rangle_\varepsilon = \theta(\tau - \tau^{^{\prime}}) \big[ E_\alpha(x), E_\beta(x^{^{\prime}}) \big] - \big\langle E_\alpha(x) E_\beta(x^{^{\prime}}) \big\rangle_\varepsilon \nonumber \\[0.3em] &= - \big\langle E_\alpha(x) E_\beta(x^{^{\prime}}) \big\rangle_\varepsilon \theta(\tau^{^{\prime}} - \tau) - \big\langle E_\beta(x^{^{\prime}}) E_\alpha(x) \big\rangle_\varepsilon \theta(\tau - \tau^{^{\prime}}) \nonumber \\[0.3em] &=: - G_{\alpha \beta}^{>}(x-x^{^{\prime}}) \theta(\tau^{^{\prime}} - \tau) - G_{\alpha \beta}^{\lt}(x-x^{^{\prime}}) \theta(\tau - \tau^{^{\prime}}). \end{align} \tag{ 4.29 }$$

Inserting $G^{++}_{\alpha \beta}(x-x^{^{\prime}})$ and $G^{-}_{\alpha \beta}(x-x^{^{\prime}})$ into equation (4.27) yields

$$\begin{align*} \tilde{\rho}_S(t) &= \mathcal{T}^{(S)}_\leftarrow \exp \bigg\{\frac{\alpha^2}{2} \int_{t_0}^t \int_{\mathbb{R}^3} d\tau d^3x \int_{t_0}^t \int_{\mathbb{R}^3} d\tau^{^{\prime}} d^3x^{^{\prime}} \sum\limits_\alpha \sum\limits_\beta \\ &\quad\times \bigg[ - \Big( G_{\alpha \beta}^{>}(x-x^{^{\prime}}) \theta(\tau - \tau^{^{\prime}}) + G_{\alpha \beta}^{\lt}(x-x^{^{\prime}}) \theta(\tau^{^{\prime}} - \tau) \Big) S_\alpha^+(x) S_\beta^+(x^{^{\prime}}) \\ &\quad - \Big( G_{\alpha \beta}^{>}(x-x^{^{\prime}}) \theta(\tau^{^{\prime}} - \tau) + G_{\alpha \beta}^{\lt}(x-x^{^{\prime}}) \theta(\tau - \tau^{^{\prime}}) \Big) S_\alpha^-(x) S_\beta^-(x^{^{\prime}}) \\ &\quad+ G_{\alpha \beta}^{\lt}(x-x^{^{\prime}}) S_\alpha^+(x) S_\beta^-(x^{^{\prime}}) + G_{\alpha \beta}^{>}(x-x^{^{\prime}}) S_\alpha^-(x) S_\beta^+(x^{^{\prime}})\bigg] \bigg\} \tilde{\rho}_S(t_0). \\ \end{align*}$$

Now we would like to eliminate the Heaviside functions $\theta(\tau-\tau^{^{\prime}})$ from this expression in order to be left purely with the thermal Wightman functions $G_{\alpha \beta}^{\lessgtr}(x-x^{^{\prime}})$. This is achieved by expanding the mixed contributions with $1 = \theta(\tau - \tau^{^{\prime}}) + \theta(\tau^{^{\prime}} - \tau)$ and renaming the integration variables in the individual terms while changing the integration domain into a nested integral [110]:

$$\begin{align} \tilde{\rho}_S(t) &= \mathcal{T}^{(S)}_\leftarrow \exp \bigg\{\alpha^2 \int_{t_0}^t \int_{\mathbb{R}^3} d\tau d^3x \int_{t_0}^\tau \int_{\mathbb{R}^3} d\tau^{^{\prime}} d^3x^{^{\prime}} \sum\limits_\alpha \sum\limits_\beta \nonumber \\ &\quad\times \bigg[ - G_{\alpha \beta}^{>}(x-x^{^{\prime}}) S_\alpha^+(x) S_\beta^+(x^{^{\prime}}) - G_{\alpha \beta}^{\lt}(x-x^{^{\prime}}) S_\alpha^-(x) S_\beta^-(x^{^{\prime}}) \nonumber \\ &\quad+ G_{\alpha \beta}^{\lt}(x-x^{^{\prime}}) S_\alpha^+(x) S_\beta^-(x^{^{\prime}}) + G_{\alpha \beta}^{>}(x-x^{^{\prime}}) S_\alpha^-(x) S_\beta^+(x^{^{\prime}})\bigg] \bigg\} \tilde{\rho}_S(t_0).\nonumber\\ \end{align} \tag{ 4.30 }$$

This concludes the derivation of the influence phase functional, while strictly speaking $\Phi_t[S^+, S^-]$ is both a functional in the $S^{\pm}$ operators as well as a superoperator that acts on the space of density operators associated to the system. Explicitly, we obtain

$$\begin{align} i\Phi_t[S^+, S^-] &:= \alpha^2 \int_{t_0}^t \int_{\mathbb{R}^3} d\tau d^3x \int_{t_0}^\tau \int_{\mathbb{R}^3} d\tau^{^{\prime}} d^3x^{^{\prime}} \sum\limits_\alpha \sum\limits_\beta \nonumber \\ &\quad\,\times \bigg[ -G_{\alpha \beta}^{>}(x-x^{^{\prime}}) S_\alpha^+(x) S_\beta^+(x^{^{\prime}}) - G_{\alpha \beta}^{\lt}(x-x^{^{\prime}}) S_\alpha^-(x) S_\beta^-(x^{^{\prime}}) \nonumber \\ &\quad\, + G_{\alpha \beta}^{\lt}(x-x^{^{\prime}}) S_\alpha^+(x) S_\beta^-(x^{^{\prime}}) + G_{\alpha \beta}^{>}(x-x^{^{\prime}}) S_\alpha^-(x) S_\beta^+(x^{^{\prime}})\bigg].\nonumber\\ \end{align} \tag{ 4.31 }$$

The influence phase functional is often depicted in a slightly different form, including functions called the dissipation- and noise kernel based on their interpretation in well-established decoherence models. The transition from the expression in (4.31) is rather straightforward. To this extent, let us first introduce $D_{\alpha \beta}(x-x^{^{\prime}})$ and $D_{\alpha \beta}^1(x-x^{^{\prime}})$ based on the thermal Wightman functions:

$$\begin{align} D_{\alpha \beta}(x-x^{^{\prime}}) &:= \big\langle \big[E_\alpha(x), E_\beta(x^{^{\prime}}) \big] \big\rangle_\varepsilon = i \big(G_{\alpha \beta}^{>}(x-x^{^{\prime}}) - G_{\alpha \beta}^{\lt}(x-x^{^{\prime}})\big), \end{align} \tag{ 4.32 }$$

$$\begin{align} D_{\alpha \beta}^1(x-x^{^{\prime}}) &:= \big\langle \big\{E_\alpha(x), E_\beta(x^{^{\prime}}) \big\} \big\rangle_\varepsilon = G_{\alpha \beta}^{>}(x-x^{^{\prime}}) + G_{\alpha \beta}^{\lt}(x-x^{^{\prime}}). \end{align} \tag{ 4.33 }$$

These relations can be readily inverted:

$$\begin{align*} G_{\alpha \beta}^{>}(x-x^{^{\prime}}) &= \frac{1}{2} \big( D_{\alpha \beta}^1(x-x^{^{\prime}}) - iD_{\alpha \beta}(x-x^{^{\prime}}) \big), \quad\nonumber\\ G_{\alpha \beta}^{\lt}(x-x^{^{\prime}}) &= \frac{1}{2} \big( D_{\alpha \beta}^1(x-x^{^{\prime}}) + iD_{\alpha \beta}(x-x^{^{\prime}}) \big). \end{align*}$$

In terms of the noise- and dissipation kernel, the influence phase functional is given by:

$$\begin{align} i\Phi_t[S^+, S^-] &:= \alpha^2 \int_{t_0}^t \int_{\mathbb{R}^3} d\tau d^3x \int_{t_0}^\tau \int_{\mathbb{R}^3} d\tau^{^{\prime}} d^3x^{^{\prime}} \sum\limits_\alpha \sum\limits_\beta \nonumber \\ &\quad\,\times \bigg[ \frac{i}{2}D_{\alpha \beta}(x-x^{^{\prime}}) \Big( S_\alpha^+(x) S_\beta^+(x^{^{\prime}}) - S_\alpha^-(x) S_\beta^-(x^{^{\prime}}) + S_\alpha^+(x) S_\beta^-(x^{^{\prime}}) - S_\alpha^-(x) S_\beta^+(x^{^{\prime}}) \Big) \nonumber \\ &\quad\,- \frac{1}{2}D_{\alpha \beta}^1(x-x^{^{\prime}}) \Big( S_\alpha^+(x) S_\beta^+(x^{^{\prime}}) + S_\alpha^-(x) S_\beta^-(x^{^{\prime}}) - S_\alpha^+(x) S_\beta^-(x^{^{\prime}}) - S_\alpha^-(x) S_\beta^+(x^{^{\prime}}) \Big) \bigg].\nonumber \\ \end{align} \tag{ 4.34 }$$

We realise that these particular combinations of $S^{\pm}$-superoperators can be cast into a familiar form in terms of combinations of commutators and anti-commutators

$$\begin{align*} S_\alpha^c(x)(\cdot) := [S_\alpha(x), \cdot] = \big( S_\alpha^+(x) - S_\alpha^-(x) \big) (\cdot), \quad S_\alpha^a(x)(\cdot) := \{S_\alpha(x), \cdot \} = \big( S_\alpha^+(x) + S_\alpha^-(x) \big) (\cdot). \end{align*}$$

Finally, the most convenient (and most abundant) version of the influence phase functional reads:

$$\begin{align} i\Phi_t[S^+, S^-] &:= \alpha^2 \int_{t_0}^t \int_{\mathbb{R}^3} d\tau d^3x \int_{t_0}^\tau \int_{\mathbb{R}^3} d\tau^{^{\prime}} d^3x^{^{\prime}} \sum\limits_\alpha \sum\limits_\beta \nonumber \\ &\quad\,\times \bigg[ \frac{i}{2}D_{\alpha \beta}(x-x^{^{\prime}}) S_\alpha^c(x) S_\beta^a(x^{^{\prime}}) - \frac{1}{2}D_{\alpha \beta}^1(x-x^{^{\prime}}) S_\alpha^c(x) S_\beta^c(x^{^{\prime}}) \bigg]. \end{align} \tag{ 4.35 }$$

At this stage it is evident that the projection operator formalism and the influence phase functional approach are equivalent formulations. Strictly speaking, this is only the case for Gaussian environmental initial states [12]. Up to second order this equality can be easily verified by deriving equation (4.31) with respect to time and consequently applying all the commutators and anti-commutators to the density operator to the right. Once we recall the definition of the kernels, this precisely amounts to the Born-Markov master equation prior to the second Markov approximation regarding the upper integral limit. A similar procedure can be implemented order by order [107].

In this section we apply the projection operator method reviewed in section 4.1 to the model in this work. Assuming factorising initial conditions and a Gibbs state for the environment, we can take the TCL master equation (4.21) to second order in the coupling strength $\alpha = \sqrt{\kappa}$ and evaluate it for the model under our consideration described by the interaction Hamiltonian operator shown in (3.17):

$$\begin{align} \frac{\partial}{\partial t} \widetilde{\rho}_S(t) = -i\kappa \left[ \widetilde{U}(t),\widetilde{\rho}_S(t) \right]- \kappa \int_{0}^t ds\; tr_{\mathcal{E}} \left\{\left[ \widetilde{H}_{TI}(t), \left[ \widetilde{H}_{TI}(s), \widetilde{\rho}_S(t) \otimes \rho_{\mathcal{E}} \right] \right] \right\}\,. \end{align} \tag{ 4.36 }$$

Now to obtain the final master equations we have to evaluate the trace over the environmental degrees of freedom in this expression and thus to obtain the correlation functions. As a first step we take into account that the second term on the right hand side of (4.36) can be written in terms of thermal Wightman functions, using the expression of $\widetilde{H}_{TI}(t)$ in position space given by

$$\begin{align} \widetilde{H}_{TI} = \sqrt{\kappa}\int_{\mathbb{R}^3} d^3x \; \delta^i_b\; \varphi^b_a(\vec{x},t)\; \delta \mathcal{E}{}{^a{}_i}(\vec{x},t), \end{align} \tag{ 4.37 }$$

where $\varphi^b_a(\vec{x},t)$ denotes the three dimensional Fourier transforms of $\underset{\smile}{\varphi}{}^b_a(\vec{k},t)$ which was defined in (3.16). Following the procedure outlined in [110], where in our case the environmental part of the interaction is linear in the environmental fields and we use a Gibbs state for the environment, hence all one-point correlation functions vanish, we can define the thermal Wightman functions:

$$\begin{align} G^{>}{}{^a{}_i{}^b{}_j}(\vec{x}-\vec{y},t-s) &:= \big\langle \delta \mathcal{E}{}{^a{}_i}(\vec{x},t) \: \delta \mathcal{E}{}{^b{}_j}(\vec{y},s) \big\rangle_\varepsilon = {{P}{^a{}_i{}^b{}_j}} G^{>}(\vec{x}-\vec{y},t-s) \end{align} \tag{ 4.38 }$$

$$\begin{align} G^{\lt}{}{^a{}_i{}^b{}_j}(\vec{x}-\vec{y},t-s) &:= \big\langle \delta \mathcal{E}{}{^b{}_j}(\vec{y},s) \: \delta \mathcal{E}{}{^a{}_i}(\vec{x},t) \big\rangle_\varepsilon = {{P}{^a{}_i{}^b{}_j}} G^{\lt}(\vec{x}-\vec{y},t-s) =G^{>}{}{^a{}_i{}^b{}_j}(\vec{y}-\vec{x},s-t)\,,\nonumber\\ \end{align} \tag{ 4.39 }$$

where $\big\langle \cdot \big\rangle_\varepsilon$ denotes the expectation value with respect to a thermal Gibbs state and ${{P}{^a{}_i{}^b{}_j}}$ is the transverse-traceless projector given in (2.69). The detailed derivations of the identities in (4.38) and (4.39) as well as the derivations for the results up to (4.46) can be found in appendix D. Note that this definition coincides with the one in (4.28) adapted to the model we consider here, hence the α and β indices from (4.28) here correspond each to the indices carried by $\delta \mathcal{E}$, i.e. to a spatial index a as well as an internal index i. The spectral representations of the thermal Wightman functions in terms of a spectral density denoted by $\rho(k^0,||\vec{k}||)$ read:

$$\begin{align} G^{>}(\vec{x}-\vec{y},t-s) &= \int_{\mathbb{R}^4} \frac{d^4k}{(2\pi)^{\frac{4}{2}}} \rho^{>}(k^0,||\vec{k}||) e^{-ik^0 (t-s) +i\vec{k}(\vec{x}-\vec{y})} \end{align} \tag{ 4.40 }$$

$$\begin{align} G^{\lt}(\vec{x}-\vec{y},t-s) &= \int_{\mathbb{R}^4} \frac{d^4k}{(2\pi)^{\frac{4}{2}}} \rho^{\lt}(k^0,||\vec{k}||) e^{-ik^0 (t-s) +i\vec{k}(\vec{x}-\vec{y})} \end{align} \tag{ 4.41 }$$

with

$$\begin{align} \rho^{>}(k^0,||\vec{k}||) := (1+ N(k^0)) \rho(k^0,||\vec{k}||) \ \ \ \rho^{\lt}(k^0,||\vec{k}||) := N(k^0) \rho(k^0,||\vec{k}||) \end{align} \tag{ 4.42 }$$

and the spectral density

$$\begin{align} \rho(k^0,||\vec{k}||) := \sqrt{\frac{\pi}{2}} \frac{1}{\Omega_k} \left[ \delta(k^0-\Omega_k) - \delta(k^0+\Omega_k) \right]\,. \end{align} \tag{ 4.43 }$$

From this decomposition, it is immediately possible to see the additional effect caused by the finite temperature in the environment, evident by the presence of the Bose–Einstein distribution

$$\begin{align} N(k^0) := \frac{1}{e^{\beta k^0}-1} = \frac{\coth\left(\frac{\beta k^0}{2}\right)}{2}-\frac{1}{2}\,, \end{align} \tag{ 4.44 }$$

which vanishes for zero temperature parameter $\Theta = \beta^{-1} = 0$ of the Gibbs state. As usual the two thermal Wightman functions can be combined to build the thermal Feynman propagator:

$$\begin{align} G^{(F)}{}{^a{}_i{}^b{}_j}(\vec{x}-\vec{y},t-s) &= \big\langle \mathcal{T}_{\leftarrow} \; \delta \mathcal{E}{}{^a{}_i}(\vec{x},t) \: \delta \mathcal{E}{}{^b{}_j}(\vec{y},s) \big\rangle_\varepsilon \nonumber\\ &= {{P}{^a{}_i{}^b{}_j}}G^{>}(\vec{x}-\vec{y},t-s) \theta(t-s) + {{P}{^a{}_i{}^b{}_j}}G^{\lt}(\vec{x}-\vec{y},t-s) \theta(s-t) \,, \end{align} \tag{ 4.45 }$$

where $\theta(t)$ denotes the Heaviside step function which is one for non-negative arguments and zero otherwise. A four-dimensional Fourier transform yields

$$\begin{align} \underset{\smile}{G}^{(F)}{}{^a{}_i{}^b{}_j}(p) = \underset{\smile}{P}{}{^a{}_i{}^b{}_j}(\vec{p}) \left[ \frac{-i}{p^\mu p_\mu-i\epsilon} + 2\pi N(\Omega_p) \delta(p^\mu p_\mu)\right]\,, \end{align} \tag{ 4.46 }$$

where the first part is the ordinary Feynman propagator one obtains at a vanishing temperature parameter Θ, that is when the Gibbs state merges into a vacuum state, and a second part the thermal contribution that obviously vanishes for $\Theta = 0$. Note that we obtain a decomposition into the ordinary Feynman propagator and a thermal correction, this is caused by the fact that we use a mode expansion that involves a splitting into positive and negative frequency modes. As a consequence, normal ordered expectation values with respect to thermal states are in general non-vanishing as for instance $\lt b_k^\dagger b_k\gt_\epsilon = N(\Omega_k) \neq 0$. Vanishing expectation values of normal ordered products are usually an important property in the proof of Wick's theorem. Given that we are working with a Gaussian state for the environment and an interaction Hamiltonian being linear in the environmental fields, we can calculate the relevant expectation values in the model we consider directly without the use of Wick's theorem. However, for different and more complicated models the application of Wick's theorem might be of advantage and for this purpose it is often convenient to consider a different splitting than in the zero temperature case such that the expectation values of normal ordered products with respect to thermal states also vanish. A detailed discussion on possible splittings in this context as well as the relation between different choices of splittings is presented in [112]. Now we can write down the second term on the right side of in the master equation (4.36) in terms of the thermal Wightman functions:

$$\begin{align} &\kappa \int_{0}^t ds\; tr_{\mathcal{E}} \left\{\left[ \widetilde{H}_{TI}(t), \left[ \widetilde{H}_{TI}(s), \widetilde{\rho}_S(t) \otimes \rho_{\mathcal{E}} \right] \right] \right\} \nonumber\\[4pt] &\quad = -\kappa \int_0^t ds\; \int_{\mathbb{R}^3} d^3x \int_{\mathbb{R}^3} d^3y \sum\limits_{r\in\{\pm\}} \Bigg\{\left[ J_r(\vec{x},t) J_r(\vec{y},s) \widetilde{\rho}_S(t)- J_r(\vec{y},s)\widetilde{\rho}_S(t) J_r(\vec{x},t) \right] \nonumber\\[4pt] &\qquad \times G^>(\vec{x}-\vec{y},t-s) +\left[ \widetilde{\rho}_S(t) J_r(\vec{y},s) J_r(\vec{x},t)- J_r(\vec{x},t)\widetilde{\rho}_S(t) J_r(\vec{y},s) \right] G^{\lt}(\vec{x}-\vec{y},t-s) \Bigg\}\,.\nonumber\\ \end{align} \tag{ 4.47 }$$

Since we consider a Gibbs state on Fock space which is not well defined, we would need to work with KMS states [85, 86] or alternatively regularise the system within a finite volume, where we will follow the latter in this work. For this purpose we put the system into a box of volume $V = L^3$ allowing us to explicitly show the above identities and evaluate thermal two-point functions. This kind of regularisation leads to the discreteness of modes that belong to the set $\mathbb{K}$ now and to a replacement of the operator-valued distributions by operators:

$$\begin{align} \int_{\mathbb{R}^3} \frac{d^3k}{(2\pi)^3} \longrightarrow \frac{1}{V} \sum\limits_{\vec{k}\in\mathbb{K}} \hspace{1in} b_k^r \longrightarrow \sqrt{\frac{V}{(2\pi)^3}}\; b_{\vec{k}}^r\,, \end{align} \tag{ 4.48 }$$

where a rescaling of the operators was introduced in order to keep the exponential $e^{-\beta H_{\mathcal{E}}}$ dimensionless in the regularised model:

$$\begin{align} e^{-\beta H_{\mathcal{E}}} = \exp\left\{-\beta \sum\limits_{r\in\{+,-\}} \sum\limits_{\vec{k}\in\mathbb{K}} \Omega_{\vec{k}} \left[(b_{\vec{k}}^r)^\dagger\; b_{\vec{k}}^r\right] \right\}\,. \end{align} \tag{ 4.49 }$$

The Gibbs state for the environment is then of the form

$$\begin{align} \rho_{\mathcal{E}} = Z_{\mathcal{E}}^{-1} e^{-\beta H_{\mathcal{E}}} \end{align} \tag{ 4.50 }$$

with the partition sum given by

$$\begin{align} Z_{\mathcal{E}} := tr_{\mathcal{E}} \left\{e^{-\beta H_{\mathcal{E}}} \right\} = \left[ \prod_{j\in\mathbb{N}} \frac{1}{1-e^{-\beta \Omega_{\vec{k}_j}}} \right]^2\,, \end{align} \tag{ 4.51 }$$

where the $j\in\mathbb{N}$ label the $\vec{k}_j\in\mathbb{K}$. Moreover, for the calculation we had to assume $\vec{k} = 0 \not\in \mathbb{K}$, which corresponds to $\Omega_{\vec{k}} = 0$ and is the usual infrared divergence. More details on computing this partition sum can be found in appendix D. Removing the regulator by taking the limit $L\to\infty$ in the Schrödinger picture we obtain the following TCL master equation for the system density operator $\rho_S(t)$:

$$\begin{align} \frac{\partial}{\partial t} \rho_S(t) &= -i \left[H_S+\kappa \: {U}, \rho_S(t)\right]-\frac{\kappa}{2}\int_{0}^t ds \sum\limits_{r} \int_{\mathbb{R}^3} d^3k \; \bigg\{i D(\vec{k},t-s)\nonumber\\ &\quad\times \left[{J}_{r}(\vec k),\left\{J_{r}(-\vec k,s-t), {\rho}_S(t)\right\}\right] + D_1(\vec{k},t-s) \left[{J}_{r}(\vec k),\left[ J_{r}(-\vec k,s-t),{\rho}_S(t)\right]\right] \bigg\}\nonumber\\ \end{align} \tag{ 4.52 }$$

with the matter system's Hamiltonian operator

$$\begin{align} H_S := \int_{\mathbb{R}^3} d^3k \; \omega_k \, a_k^\dagger \, a_k \end{align} \tag{ 4.53 }$$

as well as the operator U for the gravitational self-interaction of the scalar field that was defined in (3.12), the operator-valued distribution $J_r(\vec{k},t)$ defined in (3.18) and $J_r(\vec{k}) = J_r(\vec{k},0)$. The quantities $D(\vec{k},t-s)$ and $D_1(\vec{k},t-s)$ are the three-dimensional Fourier transforms of the dissipation and noise kernel. They can be related to the commutator's and anticommutator's correlation functions and also to the thermal Wightman functions $G^\gt(\vec{x}-\vec{y},t-s)$ and $G^\lt(\vec{x}-\vec{y},t-s)$, see also (4.32) and (4.33):

$$\begin{align} \left[\delta \mathcal{E}{}{^a{}_i}(\vec{x},t), \delta \mathcal{E}{}{^b{}_j}(\vec{y},s)\right] &= {{P}{^a{}_i{}^b{}_j}} [G^>(\vec{x}-\vec{y},t-s) - G^{\lt}(\vec{x}-\vec{y},t-s)]\nonumber\\ &=i {{P}{^a{}_i{}^b{}_j}} \int_{\mathbb{R}^3} \frac{d^3k}{(2\pi)^{\frac{3}{2}}} e^{i\vec{k}(\vec{x}-\vec{y})}\; D(\vec{k},t-s) \end{align} \tag{ 4.54 }$$

$$\begin{align} tr_{\mathcal{E}}\left(\left\{\delta \mathcal{E}{}{^a{}_i}(\vec{x},t), \delta \mathcal{E}{}{^b{}_j}(\vec{y},s)\right\}\right) &= {{P}{^a{}_i{}^b{}_j}} [G^>(\vec{x}-\vec{y},t-s) + G^{\lt}(\vec{x}-\vec{y},t-s)]\nonumber\\ &={{P}{^a{}_i{}^b{}_j}} \int_{\mathbb{R}^3} \frac{d^3k}{(2\pi)^{\frac{3}{2}}} e^{i\vec{k}(\vec{x}-\vec{y})}\; D_1(\vec{k},t-s)\,, \end{align} \tag{ 4.55 }$$

where we used that for the model considered here the commutator in the first equation is proportional to the identity operator and thus can be pulled out of the trace. Hence, we get

$$\begin{align} D(\vec{k},t-s) &:= -\frac{\sin(\Omega_k(t-s))}{\Omega_k} \end{align} \tag{ 4.56 }$$

$$\begin{align} D_1(\vec{k},t-s) &:= [2 N(\Omega_k)+1]\,\frac{\cos(\Omega_k(t-s))}{\Omega_k}\,. \end{align} \tag{ 4.57 }$$

To compare our results better with those already existing in the literature like for instance the ones in [35, 37, 38] it is of advantage to rewrite the master equation in (4.52) in two equivalent ways. The first alternative and equivalent form is given by

$$\begin{align}\begin{array}{l} \frac{\partial }{{\partial t}}{\rho _S}(t) = - i\left[ {{H_S} + \kappa U,{\rho _S}(t)} \right]\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + \frac{\kappa }{4}\int {{d^3}} k\sum\limits_r {\left\{ {\frac{{\coth \left( {\frac{{\beta {\Omega _k}}}{2}} \right)}}{{{\Omega _k}}}\left[ {{J_r}(\vec k),\left[ {{\rho _S}(t),{\textrm{ }}\tilde J_r^\dagger (\vec k,t)} \right]} \right]} \right.} + {\textrm{h}}.{\textrm{c}}.\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + \frac{1}{{{\Omega _k}}}\left. {\left[ {{J_r}(\vec k),\left\{ {{\rho _S}(t),{\textrm{ }}\tilde J_r^\dagger (\vec k,t)} \right\}} \right] + {\textrm{h}}.{\textrm{c}}.} \right\}, \end{array}\end{align} \tag{ 4.58 }$$

with

$$\begin{align} \widetilde{J}_r(\vec{k},t) := \int_0^t ds \; e^{-i\Omega_k(t-s)} J_r(\vec{k},s-t) \end{align} \tag{ 4.59 }$$

and where h.c. denotes the hermitian conjugate. This form of master equation is similar to the master equations derived in [35, 38]. While the first reference investigates also a model where a scalar field is coupled to linearised gravity, the second one replaces the scalar field by photons. Let us briefly discuss the similarities and difference of these models and the one considered here: At the classical level one of the main difference is that we formulate the model in terms of Ashtekar variables whereas the models in [35, 38] are based on ADM variables. Thus, we have to deal with an additional Gauß constraint in the model that also needs to be taken into account when gauge invariance is considered. While we presented a way to work with a gauge invariant formulation of the classical theory by means of constructing suitable Dirac observables and also showed that one can construct a canonical transformation on the original phase space that provides a natural separation of the physical and gauge sector of the phase space by the transformation in appendix A.4, in the papers [35, 38] a specific gauge fixing is used to eliminate the gauge degrees of freedom. Furthermore, we consider a standard boundary term [94] that ensures that the action in terms of Ashtekar variables as well as its variation is well defined and we expect that the corresponding ADM boundary term should be included in the formulation of the models in [35, 38] as well. Note that in our case the boundary terms cancels the term in $\mathbf{\delta H}_{can}$ that is linear in the linearised geometric contribution to the Hamiltonian constraint and the remaining matter part combines with the $\frac{1}{\kappa}$ in front of the action to the energy density of the scalar field in κ⁰ order. In contrast in [35, 38] no perturbations of linear order have been considered in the canonical Hamiltonian but the κ⁰ order is present. Therefore, effects of the missing boundary term cannot be seen in the final master equation when comparing the two results.

For the quantised model let us comment on the work in [35, 38] separately, starting with the latter. One difference in the Hamiltonian operators in our result and the one in [35] is the form of the self-interaction part because we choose a different normal ordering compared to [35]. We chose to normal order the entire self-interaction part of the Hamiltonian operator, whereas in [35] only the individual contributions $\hat{\underset{\smile}{p}}{}_a$, $\hat{\underset{\smile}{\epsilon}}$, $\hat{\underset{\smile}{\widetilde{V}}}$ and $\underset{\smile}{\hat{\varphi}}{}^a_b$ were normal ordered. Albeit the current operator valued distribution J_r defined in (3.11) is the same as in [35], it appears in their Hamiltonian (3.10) with a different factor of $\sqrt{2}$ compared to our coupling. This difference has no physical effect because it is absorbed at the level of the master equation for the following reason: The different factor of $\sqrt{2}$ arises because the Hamiltonian for the pure gravitational part contains an additional factor of $\frac{1}{2}$ when expressed in Ashtekar variables and working with $\delta \mathcal{E}$ and $\delta \mathcal{C}$ compared to the one in ADM variables in [35]. This leads to an additional factor $\frac{1}{\sqrt{2}}$ in our Fock quantised gravitational variables compared to the ADM variables. In addition in terms of ADM variables the interaction term in the Hamiltonian density reads $-\frac{\kappa}{2} \delta h_{ab}^{TT}\, T^{ab}$ while in Ashtekar variables it is given by $+\kappa \delta \mathcal{E}{}{^a{}_i} \delta^{bi}\: T_{ab}$. This difference can be explained by the relation between ADM and Ashtekar variables. In case of ADM variables the gravitational physical degrees of freedom are given by the transverse-traceless components $\delta h^{TT}_{ab} : = P_{ab}^{cd} \delta h_{cd}$ and the interaction reads $\delta h^{TT}_{ab} T^{ab}$ which is then quantised in [35] and where the coupling to the gravitational gauge degrees of freedom to the energy momentum tensor is contained in the self-interaction part. In contrast in our case, writing $\delta h_{ab}$ in terms of the perturbed cotriad, that is $\delta h_{ab}(E) = -\delta_b^i\delta_{ac} \delta{{E}{^c{}_i}}-\delta_a^i\delta_{bc} \delta{{E}{^c{}_i}}+\delta_{ab} \delta_c^i {{\delta E}{^c{}_i}}$ and then using symmetry of the energy-momentum tensor, we obtain $\delta h_{ab}(E) T^{ab} \rightarrow [-2\delta_a^i \delta_{bc} +\delta_{ab}\delta_c^i]{{\delta E}{^c{}_i}} T^{ab}$, where the gravitational physical degrees of freedom only enter in the first term and we obtain in the linearised theory just $-2\delta {\mathcal{E}}{^a{}_i} \delta_b^i T^b_a$. This difference already present in the classical Hamiltonian has no effect on the final equations of motion because with our convention the Poisson bracket between $\delta {\mathcal{A}}{_a{}^i},\delta {\mathcal{E}}{^a{}_i}$ involves an additional factor of $\frac{1}{2}$ compared to the Poisson brackets used in [35, 38], cancelling the additional factor of 2. The minus sign is also cancelled because in our case the momentum variables couple to the energy momentum tensor whereas in [35, 38] the coupling is via configuration variables. This kind of cancellation carries of course over to the commutator at the quantum level.

A further difference is that the final master equation presented in [35] is of Lindblad form (4.22). This is a significant difference because as we will show in the following, the master equation in (4.52) we obtain is not of Lindblad type yet and it requires further assumptions and approximations respectively to be of this form. For the purpose of discussing this we rewrite the master equation in (4.52) in another equivalent form

$$\begin{align} \frac{\partial}{\partial t} \rho_S(t) & = -i [H_S+\kappa\: U, \rho_S(t)] -\kappa \sum\limits_{r} \int_{\mathbb{R}^3} \frac{d^3k}{2 \Omega_k} \;\nonumber\\ & \quad \times \Big\{[J_r^\dagger(\vec{k}),\widetilde{J}_r(\vec{k},t)\rho_S(t)] + N(\Omega_k) [J_r^\dagger(\vec{k}), [\widetilde{J}_r(\vec{k},t),\rho_S(t)]]+ \textrm{h.c.} \Big\}\,. \end{align} \tag{ 4.60 }$$

In this form, one can clearly separate the dissipator into a part associated with zero temperature (terms without $N(\Omega_k)$) and a second part corresponding to finite temperature (terms with $N(\Omega_k)$). All three forms of the master equation derived in this paper are time-convolutionless. However, they cannot (yet) be brought into Lindblad form, as there are still time-dependent quantities apart from the reduced density operator itself present in the dissipator. Nevertheless, we can bring the master equation into a form which is similar to the first standard form introduced in (4.22) and which admits the splitting of the dissipator into a Lamb-shift Hamiltonian and a remaining contribution. To this extent, we start from the formulation in (4.60). First, we note that the J_r -operators can be split in the following way into pairs of creation and / or annihilation operator-valued distributions:

$$\begin{align} J_r(\vec{k}) &= -\int \frac{d^3p}{(2\pi)^\frac{3}{2}} \sum\limits_{a=1}^4 j_r^a(\vec{k},\vec{p}) \end{align} \tag{ 4.61 }$$

$$\begin{align} \widetilde{J}_r(\vec{k},t) &= -\int \frac{d^3p}{(2\pi)^\frac{3}{2}} \sum\limits_{a=1}^4 j_r^a(\vec{k},\vec{p}) \; f(\Omega_k + \omega_a(\vec{k},\vec{p});t) \end{align} \tag{ 4.62 }$$

with

$$\begin{align} j_r^1(\vec{k},\vec{p})&:= a_p^\dagger a_{k+p} \frac{1}{2\sqrt{\omega_p \omega_{k+p}}} \left[p_ap^b [P^{-r}(\vec{k})]^a_b \right]\quad\;\;\; \omega_1(\vec{k},\vec{p}) := \omega_p - \omega_{k+p} \end{align} \tag{ 4.63 }$$

$$\begin{align} j_r^2(\vec{k},\vec{p})&:= a_{-p-k}^\dagger a_{-p} \frac{1}{2\sqrt{\omega_p \omega_{k+p}}} \left[p_ap^b [P^{-r}(\vec{k})]^a_b \right] \;\; \omega_2(\vec{k},\vec{p}) := \omega_{k+p} - \omega_{p} \end{align} \tag{ 4.64 }$$

$$\begin{align} j_r^3(\vec{k},\vec{p})&:= a_{-p} a_{k+p} \frac{1}{2\sqrt{\omega_p \omega_{k+p}}} \left[p_ap^b [P^{-r}(\vec{k})]^a_b \right]\quad \omega_3(\vec{k},\vec{p}) := -\omega_p - \omega_{k+p} \end{align} \tag{ 4.65 }$$

$$\begin{align} j_r^4(\vec{k},\vec{p})&:= a_p^\dagger a_{-k-p}^\dagger \frac{1}{2\sqrt{\omega_p \omega_{k+p}}} \left[p_ap^b [P^{-r}(\vec{k})]^a_b \right]\quad \omega_4(\vec{k},\vec{p}) := \omega_p + \omega_{k+p} \end{align} \tag{ 4.66 }$$

and

$$\begin{align} f(\omega;t) := \int_0^t ds\; e^{-i\omega(t-s)} = \frac{i}{\omega} (e^{-i\omega t} -1)\,. \end{align} \tag{ 4.67 }$$

Using this expansion in (4.60), the dissipator can be rewritten as

$$\begin{align} \mathcal{D}[\rho_S]= \frac{\kappa}{2}\sum\limits_{r;a,b} \int_{\mathbb{R}^3} \frac{d^3k\: d^3p\: d^3l}{(2\pi)^\frac{6}{2}} \Bigg\{\Delta_{ab}(\vec{p},\vec{l};\vec{k},t) \left[j_r^b(\vec{k},\vec{l})\rho_S(t), j^a_r(\vec{k},\vec{p})^\dagger\right]+h.c.\Bigg\} \end{align} \tag{ 4.68 }$$

with

$$\begin{align} \Delta_{ab}(\vec{p},\vec{l};\vec{k},t) & := \frac{1}{\Omega_k}\left[(N(\Omega_k)+1) \:f(\Omega_k+\omega_b(\vec{k},\vec{l});t) +N(\Omega_k)\;f(-\Omega_k+\omega_b(\vec{k},\vec{l});t) \right]\nonumber\\ & = 2 \int_0^t ds\; G^>(\vec{k},t-s) e^{-i\omega_b(\vec{k},\vec{l})(t-s)} \,, \end{align} \tag{ 4.69 }$$

where $G^\gt(\vec{k},t-s)$ denotes the three dimensional Fourier transform of the Wightman function $G^\gt(\vec{x}-\vec{y},t-s)$. Combining the new coefficient functions in the following form

$$\begin{align} S_{ab}(\vec{p},\vec{l};\vec{k},t) & := \frac{1}{2i} \left( \Delta_{ab}(\vec{p},\vec{l};\vec{k},t) - \Delta^*_{ba}(\vec{l},\vec{p};\vec{k},t) \right)\nonumber\\[6pt] & = \frac{1}{2\Omega_k} \Bigg[ [N(\Omega_k)+1] \Bigg\{\frac{e^{-i(\Omega_k + \omega_b(\vec k,\vec l))t}-1}{\Omega_k + \omega_b(\vec k,\vec l)} + \frac{e^{i(\Omega_k + \omega_a(\vec k,\vec p))t}-1}{\Omega_k + \omega_a(\vec k,\vec p)}\Bigg\}\nonumber\\[6pt] & \qquad\quad \ \ - N(\Omega_k) \Bigg\{\frac{e^{i(\Omega_k - \omega_b(\vec k,\vec l))t}-1}{\Omega_k - \omega_b(\vec k,\vec l)} + \frac{e^{-i(\Omega_k - \omega_a(\vec k,\vec p))t}-1}{\Omega_k - \omega_a(\vec k,\vec p)}\Bigg\}\Bigg] \end{align} \tag{ 4.70 }$$

$$\begin{align} R_{ab}(\vec{p},\vec{l};\vec{k},t) & := \Delta_{ab}(\vec{p},\vec{l};\vec{k},t) + \Delta^*_{ba}(\vec{l},\vec{p};\vec{k},t)\nonumber\\[6pt] & = \frac{i}{\Omega_k} \Bigg[ [N(\Omega_k)+1] \Bigg\{\frac{e^{-i(\Omega_k + \omega_b(\vec k,\vec l))t}-1}{\Omega_k + \omega_b(\vec k,\vec l)} - \frac{e^{i(\Omega_k + \omega_a(\vec k,\vec p))t}-1}{\Omega_k + \omega_a(\vec k,\vec p)}\Bigg\}\nonumber\\[6pt] & \qquad\quad - N(\Omega_k) \Bigg\{\frac{e^{i(\Omega_k - \omega_b(\vec k,\vec l))t}-1}{\Omega_k - \omega_b(\vec k,\vec l)} - \frac{e^{-i(\Omega_k - \omega_a(\vec k,\vec p))t}-1}{\Omega_k - \omega_a(\vec k,\vec p)}\Bigg\}\Bigg]\,, \end{align} \tag{ 4.71 }$$

we can split the dissipator into two parts:

$$\begin{align} \mathcal{D}[\rho_S] = -i\kappa[H_{LS},\rho_S(t)] + \mathcal{D}_\textrm{first}[\rho_S] \end{align} \tag{ 4.72 }$$

with the Lamb-shift Hamiltonian operator

$$\begin{align} H_{LS} := \frac{1}{2} \int \frac{d^3k\, d^3p\, d^3l}{(2\pi)^\frac{6}{2}} \sum\limits_{r;a,b} \; S_{ab}(\vec{p},\vec{l};\vec{k},t) \; j_r^a(\vec{k},\vec{p})^\dagger\, j_r^b(\vec{k},\vec{l}) \end{align} \tag{ 4.73 }$$

and a new dissipator term being in a form similar to the first standard form:

$$\begin{align} \mathcal{D}_\textrm{first}[\rho_S] & := \frac{\kappa}{2} \int \frac{d^3k\, d^3p\, d^3l}{(2\pi)^\frac{6}{2}} \sum\limits_{r;a,b} \; {R}_{ab}(\vec{p},\vec{l};\vec{k},t)\nonumber\\ & \quad \times \left( j_r^b(\vec{k},\vec{l}) \rho_S(t)j^a_r(\vec{k},\vec{p})^\dagger-\frac{1}{2} \left\{j^a_r(\vec{k},\vec{p})^\dagger j_r^b(\vec{k},\vec{l}), \rho_S(t)\right\} \right)\,, \end{align} \tag{ 4.74 }$$

where the label'first' refers to the dissipator in first standard form. The difference to the first standard form shown in (4.22) is that the coefficient function R_ab is still time-dependent and thus the corresponding master equation is not of Lindblad type. Often a Lindblad form can be obtained for a given master equation by applying the second Markov approximation, that is by formally sending the upper limit of the integral involved in (4.69) or directly t in (4.70) and (4.71) to infinity leading to time-independent coefficients. Note that this is a priori problematic here as the complex exponentials in (4.70) and (4.71) do not have a well-defined limit for $t\rightarrow \infty$. This can be considered in the context of distributions along the lines of the Fourier transform of the step functions, see for instance [113] and this will be analysed together with the remaining integration over the modes in detail in future work in [114]. If we ignore such subtleties for the moment and simply send t formally to infinity then the master equation from [35] coincides structurally with the one in (4.58) apart from the explicit form of their last term, which is is given by $\{J_r,[\rho_S,\widetilde{J}_r^\dagger]\} + h.c.$. In contrast in our derivation we actually end up with $[J_r,\{\rho_S,\widetilde{J}_r^\dagger\}] +h.c.$, which yields different terms and from our calculation we do not see a way how their result can be reproduced. In [35] no detailed derivation of this result is presented but the authors cite [12], where however also such terms are not involved in the master equations.

As already mentioned, the master equation in [38] involves a photon field instead of a scalar field and thus the projectors involved are those that project on the physical degrees of freedom of the photon field. Most of the differences have already been pointed out above because the procedure in [38] is quite similar to the one followed in [35]. Compared to [35], more details on the derivation of the Lindblad equation are presented in [38] and it is derived by means of a Born-Markov and a (weak) rotating wave approximation. The form of their master equations after the Born-Markov approximation is structurally the same as the one in (4.58), also the commutator-anticommutator structure agrees.

Note that the master equation in (4.60) looks, apart from the expected differences due to their usage of ADM variables, similar to the one derived in [37], where a Dirac quantisation was carried out and the physical degrees of freedom were identified by imposing the gauge conditions in the quantum theory and specialising to the transverse-traceless degrees of freedom. In a second step, they then derived a master equation using the influence functional approach for a general bosonic field. Their final master equation is a TCL master equation that, similar to the result derived in this paper, does also not exhibit Lindblad form.

The physical investigation of the final master equation is complex, as the equation is in general not completely positive, in contrast to an equation of Lindblad-type. Apart from applying a second Markov approximation in order to arrive at second standard form, it is required to diagonalise the coefficient functions with respect to the labels $(a,\vec{p})$ and $(b,\vec{l})$. In this context an application of a rotating wave approximation is beneficial as has for instance been used in [38]. For the model considered in our work this needs a further detailed analysis because applying the above-proposed additional approximations to cast it into a completely positive form in general one looses some features of the dynamics of the system. Hence, we postpone this kind of detailed analysis into our companion paper [114] where we also will discuss the special case of a one-particle master equations.

Defining the smeared linearised Gauß, Hamiltonian and spatial diffeomorphism constraints as

$$\begin{align} \delta G[\delta \Lambda](t) &:= \int_\sigma d^3y\; \delta \Lambda^i(\vec{x},t)\, \delta G_i(\vec{x},t) \end{align} \tag{ A.1 }$$

$$\begin{align} \delta C[\delta N](t) &:= \int_\sigma d^3x \; \delta N(\vec{x},t) \,\delta C(\vec{x},t) \end{align} \tag{ A.2 }$$

$$\begin{align} \overrightarrow{\delta C}[\overrightarrow{\delta N}](t) &:= \int_\sigma d^3x\; \delta N^a(\vec{x},t)\, \delta C_a(\vec{x},t) \end{align} \tag{ A.3 }$$

one can evaluate the gauge transformations they infer on the phase space variables. As the Gauß constraint does not contain matter degrees of freedom, it leaves the matter fields unchanged and just modifies the geometrical degrees of freedom:

$$\begin{align} \{\phi(\vec{x},t),\delta G[\delta \Lambda](t)\} =\{\pi(\vec{x},t),\delta G[\delta \Lambda](t)\} &=0 \end{align} \tag{ A.4 }$$

$$\begin{align} \{{{\delta E}{^a{}_i}}(\vec{x},t) ,\delta G[\delta \Lambda](t)\} &= \frac{1}{2} \epsilon^{ijk} {\delta \Lambda}{_j}(\vec{x},t)\, \delta_k^a \end{align} \tag{ A.5 }$$

$$\begin{align} \{{{\delta A}{_a{}^i}}(\vec{x},t) , \delta G[\delta \Lambda](t)\} &=- \frac{1}{2}{{\partial}{_a}} (\delta \Lambda^i(\vec{x},t)) \,. \end{align} \tag{ A.6 }$$

For the Hamiltonian constraint we find the following gauge transformations:

$$\begin{align} \{\phi(\vec{x},t),\delta C[ \delta N](t) \} &= \kappa \delta N(\vec{x},t) \, \pi(\vec{x},t) \end{align} \tag{ A.7 }$$

$$\begin{align} \{\pi(\vec{x},t),\delta C[ \delta N](t) \} &= \kappa \partial_a \left[\delta N(\vec{x},t)\, \partial^a\phi(\vec{x},t)\right] -\kappa \delta N(\vec{x},t)\, m^2 \,\phi(\vec{x},t) \end{align} \tag{ A.8 }$$

$$\begin{align} \{{{\delta E}{^a{}_i}}(\vec{x},t) ,\delta C[\delta N](t)\} &= -\beta\, \epsilon^{ijk}\, \delta_j^a\, \delta_k^b\, {{\partial}{_b}} (\delta N(\vec{x},t)) \end{align} \tag{ A.9 }$$

$$\begin{align} \{{{\delta A}{_a{}^i}}(\vec{x},t) ,\delta C[\delta N](t)\} &=0\,. \end{align} \tag{ A.10 }$$

Finally, the gauge transformations induced by the spatial diffeomorphism constraint are

$$\begin{align} \{\phi(\vec{x},t),\overrightarrow{\delta C}[\overrightarrow{\delta N}](t)\} &=\kappa \delta N^a(\vec{x},t)\, \partial_a \phi(\vec{x},t) \end{align} \tag{ A.11 }$$

$$\begin{align} \{\pi(\vec{x},t),\overrightarrow{\delta C}[\overrightarrow{\delta N}](t) \} &= \kappa \partial_a [\delta N^a(\vec{x},t)\, \pi(\vec{x},t)] \end{align} \tag{ A.12 }$$

$$\begin{align}\{{{\delta E}{^a{}_i}}(\vec{x},t),\overrightarrow{\delta C}[\overrightarrow{\delta N}](t) \} &= -\left( \delta_i^b\, {{\partial}{_b}} (\delta N^a(\vec{x},t)) - \delta_i^a\, {{\partial}{_b}} (\delta N^b(\vec{x},t)) \right) \end{align} \tag{ A.13 }$$

$$\begin{align} \{{{\delta A}{_a{}^i}}(\vec{x},t),\overrightarrow{\delta C}[\overrightarrow{\delta N}](t) \} &=0\,. \end{align} \tag{ A.14 }$$

A.1.1. Choice of the reference field for the linearised Gauß constraint.

As the choice of reference fields can be understood as gauge fixing constraints in a corresponding gauge fixed theory, we choose a reference field for the Gauß constraint that implements a Lorentz-like gauge condition analogous to [61] for the connection perturbation, i.e.

$$\begin{align} \partial^a (\delta {A}{_a{}^i}(\vec{x},t))=0\,. \end{align} \tag{ A.15 }$$

Evaluation of the Poisson bracket between this condition and the Gauß constraint (2.36) using the linearised Poisson bracket (2.18) yields

$$\begin{align} \{\partial^a (\delta {A}{_a{}^i}(\vec{x},t)), \delta G_j(\vec{y},t)\} =- \frac{1}{2} \delta_j^i \Delta_{\vec{x}} \delta^3(\vec{x}-\vec{y})\,, \end{align} \tag{ A.16 }$$

where $\Delta_{\vec{x}}$ denotes the Laplacian with respect to $\vec{x}$. We will drop the subscript if only one coordinate is involved. In order to have the commutation relation

$$\begin{align} \{\delta\Xi^i(\vec{x},t), \delta G_j(\vec{y},t)\} = \frac{1}{\kappa}\delta_j^i \delta^3(\vec{x}-\vec{y})\,, \end{align} \tag{ A.17 }$$

we define the reference fields for the Gauß constraint as

$$\begin{align} \delta\Xi^i(\vec{x},t) := \frac{2}{\kappa} \partial^a \left(\delta {A}{_a{}^i} * G^\Delta\right)(\vec{x},t) = \frac{2}{\kappa} \int d^3y \; \partial_{\vec{x}}^a G^\Delta(\vec x -\vec y) \, {{\delta A}{_a{}^i}}(\vec y) \end{align} \tag{ A.18 }$$

with $\partial_{\vec{x}}^a : = \frac{\partial}{\partial x_a}$. The abbreviation $G^\Delta(\vec x-\vec y)$ denotes the Green's function of the Laplacian,

$$\begin{align} G^\Delta(\vec x -\vec y) = \int \frac{d^3k}{(2\pi)^3} \frac{1}{||\vec k||^2} e^{i\vec{k}(\vec{x}-\vec{y})}\,. \end{align} \tag{ A.19 }$$

Due to the fact that the connection remains invariant under the Hamiltonian and spatial diffeomorphism constraint, this reference field also remains invariant under transformations induced by these two constraints.

A.1.2. Choice of the reference field for the linearised Hamiltonian constraint.

The Hamiltonian constraint leaves the connection invariant and only transforms the part ${{\epsilon}{_a{}^{cb}}} \delta_c^i \partial_b \delta E{}{^a{}_i}(\vec{x},t)$ of the densitised triad in the following manner:

$$\begin{align} \{{{\epsilon}{_a{}^{cb}}}\delta_c^i \partial_b \delta E{}{^a{}_i}(\vec{x},t), \delta C(\vec{y},t)\} = 2\beta \Delta_{\vec{x}} \delta^3(\vec{x}-\vec{y})\,. \end{align} \tag{ A.20 }$$

This suggests to define a reference field

$$\begin{align} \delta\widetilde{T}(\vec{x},t) := -\frac{1}{2\beta\kappa} {{\epsilon}{_a{}^{cb}}}\delta_c^i \partial_b \left(\delta E{}{^a{}_i} * G^\Delta \right)(\vec{x},t) \,. \end{align} \tag{ A.21 }$$

However, in this form the reference field is not invariant under gauge transformations induced by the Gauß constraint:

$$\begin{align} \{\delta\widetilde{T}(\vec{x},t), \delta G_j(\vec{y},t)\} =-\frac{1}{2\beta\kappa } \delta_j^a\partial_a^{\vec{x}} G^\Delta(\vec{x}-\vec{y})\,, \end{align} \tag{ A.22 }$$

where $\partial_a^{\vec{x}}$ denotes the partial derivative with respect to x^a . To cure this, we seek to subtract some combination of geometrical phase space variables that transforms precisely the same way as $\delta\widetilde{T}$ under the Gauß constraint and remains invariant under the Hamiltonian constraint. The latter is true for any combination of the connection variables, and it turns out that the trace of the connection solves the problem. Hence a good choice for a reference field corresponding to the Hamiltonian constraint is

$$\begin{align} \delta T(\vec{x},t) := -\frac{1}{\kappa\beta} \left[ \frac{1}{2} {{\epsilon}{_a{}^{cb}}} \delta_c^i\partial_b \left({{\delta E}{^a{}_i}} * G^\Delta\right)(\vec{x},t) + \delta^a_i \left({{\delta A}{_a{}^i}} * G^\Delta\right)(\vec{x},t)\right]\,, \end{align} \tag{ A.23 }$$

which also commutes with the spatial diffeomorphism constraint. Additionally, it commutes with the reference field for the Gauß constraint:

$$\begin{align} \{\delta T(\vec{x},t), \Xi^i(\vec{y},t)\} =0\,. \end{align} \tag{ A.24 }$$

A.1.3. Choice of the reference field for the linearised spatial diffeomorphism constraint.

In a last step we have to find a suitable reference field for the spatial diffeomorphism constraint. The list of requirements for this field is motivated by the application in the main part of this work and we will construct it in several steps. In the end, the final reference field $\delta T^a(\vec{x},t)$ should

• consists of linearised elementary gravitational degrees of freedom only
• fulfill $\{\delta T^a(\vec{x},t), \delta C_b(\vec{y},t)\} = \frac{1}{\kappa} \delta_b^a \delta^3(\vec{x}-\vec{y})$,
• commute with the remaining linearised constraints and
• commute with the remaining reference fields.

For the first step it is helpful to realise from (A.13) that only the trace and one of the longitudinal parts of the densitised triad will be modified by this constraint in the following manner:

$$\begin{align} \{\delta^i_c\partial^c \delta {E}{^a{}_i}(\vec{x},t) , \delta C_b(\vec{y},t) \}&= (\partial^a_{\vec{x}} \partial_b^{\vec{x}} - \delta_b^a \Delta_{\vec{x}} ) \delta^3(\vec{x}-\vec{y}) \end{align} \tag{ A.25 }$$

$$\begin{align} \{\delta_a^i \delta {E}{^a{}_i}(\vec{x},t) , \delta C_b(\vec{y},t) \}&=2 \partial_b^{\vec{x}} \delta^3(\vec{x}-\vec{y})\,. \end{align} \tag{ A.26 }$$

A suitable combination for a reference field that is canonically conjugated to the linearised spatial diffeomorphism constraint is therefore

$$\begin{align} \delta \widetilde{\widetilde{T}} {}^a(\vec{x},t):=\frac{1}{\kappa} \left( \delta^a_b \delta_c^i \partial^c -\frac{1}{2} \delta^i_b \partial^a \right) \left(\delta {E}{^b{}_i} * G^\Delta \right)(\vec{x},t)\,. \end{align} \tag{ A.27 }$$

Unfortunately, this combination is not invariant under the linearised Gauß constraint:

$$\begin{align} \{\delta \widetilde{\widetilde{T}} {}^a(\vec{x},t), \delta G_j(\vec{y},t)\} = -\frac{1}{2\kappa} {{\epsilon}{_{cb}{}^a}}\delta^c_j \partial^b_{\vec{x}} G^\Delta(\vec{x}-\vec{y})\,. \end{align} \tag{ A.28 }$$

To cure this, we can add a suitable form of

$$\begin{align} \partial_a \delta {E}{^a{}_i}(\vec{x},t)\,, \end{align} \tag{ A.29 }$$

which is just a scalar density and therefore invariant under the linearised spatial diffeomorphism constraint. It turns out that

$$\begin{align} \delta \widetilde{T}^a(\vec{x},t):= \frac{1}{\kappa} \left( \delta^a_b \delta^i_c \partial^c -\frac{1}{2} \delta^i_b \partial^a + \delta^{ac}\delta_c^i \partial_b \right) \left(\delta {E}{^b{}_i} * G^\Delta\right)(\vec{x},t) \end{align} \tag{ A.30 }$$

is invariant under the linearised Gauß constraint and still canonically conjugated to the linearised spatial diffeomorphism constraint. A quick calculation also shows that it is invariant under the gauge transformations generated by the linearised Hamiltonian constraint. The last bullet point which we required, namely vanishing Poisson bracket with the other reference fields, is not yet fulfilled. To construct an extension of (A.30) that commutes with the remaining reference fields in the lowest order and maintains the additional requirements listed in the beginning, we use the linearised dual observable map introduced in (2.48). The advantage is that we already know that the constraints are abelian up to first order in κ, so this extension will have exactly the same behaviour as $\delta\widetilde{T}^a(\vec{x},t)$ in the Poisson brackets with the linearised constraints. As the clock in the dual observable map we choose the geometric contribution of the linearised Hamiltonian constraint given by $\delta C-\kappa\epsilon$ because we want the final clock to depend on the geometrical degrees of freedom only. Since we have $\{\delta T(x),(\delta C-\kappa\epsilon)(y)\} = \frac{1}{\kappa}\delta(x,y)$ this is a suitable choice. The result of the application of the dual linearised observable map on $\delta\widetilde{T}{}^a$ yields a suitable reference field corresponding to the spatial diffeomorphism constraint:

$$\begin{align} \delta T^a(\vec{x},t) & := \frac{1}{\kappa} \left( \delta^a_b \delta^i_c \partial^c -\frac{1}{2} \delta^i_b \partial^a + \delta^{ac}\delta_c^i \partial_b \right) \left(\delta {E}{^b{}_i} * G^\Delta\right)(\vec{x},t) \nonumber\\[4pt] & \quad +\frac{2\beta}{\kappa^2}\left[ \frac{1}{2} \delta^i_b \partial^a \partial^b \left( \delta G_i * G^{\Delta\Delta}\right)(\vec x,t) - \delta^{ab}\delta_b^i \left( \delta G_i * G^\Delta\right)(\vec x,t) \right] \nonumber\\[4pt] & \quad + \frac{1}{2\kappa^2} \partial^a \left[\left(\delta C-\kappa\epsilon\right) * G^{\Delta\Delta} \right](\vec{x},t)\,, \end{align} \tag{ A.31 }$$

where $G^{\Delta\Delta}$ denotes the Green's function of the squared Laplacian, that is

$$\begin{align} G^{\Delta\Delta}(\vec x-\vec y) = \int \frac{d^3k}{(2\pi)^3} \frac{1}{||\vec k||^4} e^{i\vec{k}(\vec x-\vec y)} = \int d^3z \; G^\Delta(\vec x -\vec z) \, G^\Delta(\vec y-\vec z)\,. \end{align} \tag{ A.32 }$$

This reference field is defined on the entire phase space. On the constraint hypersurface, the additional constraint terms drop and we are left again with $\delta\widetilde{T}^a(\vec{x},t)$, which still commutes with the constraints in the desired way. This concludes the construction of the reference fields for the model under consideration.

When discussing the physical degrees of freedom and also in the context of gauge fixing and the relation between the fully gauge invariant and gauge fixed model, it turns out to be helpful to work in Fourier space. For these calculations we choose a certain complex basis in Fourier space, which is the same as in [61] consisting of

$$\begin{align} \hat{k}^a := \frac{k^a}{||\vec{k}||} \hspace{0.3in} m^a(\vec{k}) \hspace{0.3in} \overline{m}^a(\vec{k})\,, \end{align} \tag{ A.33 }$$

where $\vec{k}$ is the momentum and we understand the plane perpendicular to it as a complex plane spanned by the two unit vectors $m^a(\vec{k})$ and $\overline{m}^a(\vec{k})$, where the bar denotes complex conjugation. This basis is orthonormal in the following manner:

$$\begin{align} \hat{k}^a \hat{k}_a = m^a(\vec{k}) \overline{m}_a(\vec{k}) = \overline{m}^a(\vec{k}) m_a(\vec{k}) &= 1 \end{align} \tag{ A.34 }$$

$$\begin{align} \hat{k}^a m_a(\vec{k}) = \hat{k}^a \overline{m}_a(\vec{k}) = m^a(\vec{k}) m_a(\vec{k}) = \overline{m}^a(\vec{k}) \overline{m}_a(\vec{k}) &= 0\,. \end{align} \tag{ A.35 }$$

The expansion of the metric in this basis is

$$\begin{align} q_{ab}(\vec{k}) = \hat{k}_a \hat{k}_b + m_a(\vec{k})\, \overline{m}_b(\vec{k}) + \overline{m}_a(\vec{k})\, m_b(\vec{k})\,, \end{align} \tag{ A.36 }$$

the soldering form reads

$$\begin{align} \delta_a^i(\vec{k}) = \hat{k}_a \hat{k}^i + m_a(\vec{k})\, \overline{m}^i(\vec{k}) + \overline{m}_a(\vec{k})\, m^i(\vec{k})\,, \end{align} \tag{ A.37 }$$

and we fix the orientation of the basis vectors in the complex plane relative to each other by imposing

$$\begin{align} \epsilon_{abc}\, \hat{k}^a\, m^b(\vec{k})\, \overline{m}^c(\vec{k}) = -i\,. \end{align} \tag{ A.38 }$$

We additionally impose (see [120] section 5.3):

$$\begin{align} m^a(-\vec{k})= \overline{m}^a(\vec{k})\,. \end{align} \tag{ A.39 }$$

We expand the Fourier transform of the gravitational fields in this basis, following the notation in [61]:

$$\begin{align} \underset{\smile}{\delta E}{}{^a{}_i}(\vec{k}) & = \underset{\smile}{\delta E}^{+}(\vec{k}) m^a(\vec{k})m_i(\vec{k}) + \underset{\smile}{\delta E}^{-}(\vec{k}) \overline{m}^a(\vec{k}) \overline{m}_i(\vec{k}) + \underset{\smile}{\delta E}^1(\vec{k}) \hat{k}^a m_i(\vec{k}) + \underset{\smile}{\delta E}^{\overline{1}}(\vec{k}) \hat{k}^a \overline{m}_i(\vec{k}) \nonumber\\ & \quad + \underset{\smile}{\delta E}^2(\vec{k}) m^a(\vec{k}) \hat{k}_i + \underset{\smile}{\delta E}^{\overline{2}}(\vec{k}) \overline{m}^a(\vec{k}) \hat{k}_i + \underset{\smile}{\delta E}^3(\vec{k}) \hat{k}^a \hat{k}_i \nonumber\\ & \quad + \underset{\smile}{\delta E}^4(\vec{k}) m^a(\vec{k}) \overline{m}_i(\vec{k}) + \underset{\smile}{\delta E}^5(\vec{k}) \overline{m}^a(\vec{k}) m_i(\vec{k}) \end{align} \tag{ A.40 }$$

and

$$\begin{align} \underset{\smile}{\delta A}{}{_a{}^i}(\vec{k}) & = \underset{\smile}{\delta A}^{+}(\vec{k}) m_a(\vec{k})m^i(\vec{k}) + \underset{\smile}{\delta A}^{-}(\vec{k}) \overline{m}_a(\vec{k}) \overline{m}^i(\vec{k}) + \underset{\smile}{\delta A}^1(\vec{k}) \hat{k}_a m^i(\vec{k}) + \underset{\smile}{\delta A}^{\overline{1}}(\vec{k}) \hat{k}_a \overline{m}^i(\vec{k}) \nonumber\\[4pt] & \quad + \underset{\smile}{\delta A}^2(\vec{k}) m_a(\vec{k}) \hat{k}^i + \underset{\smile}{\delta A}^{\overline{2}}(\vec{k}) \overline{m}_a(\vec{k}) \hat{k}^i + \underset{\smile}{\delta A}^3(\vec{k}) \hat{k}_a \hat{k}^i\nonumber\\[4pt] & \quad + \underset{\smile}{\delta A}^4(\vec{k}) m_a(\vec{k}) \overline{m}^i(\vec{k}) + \underset{\smile}{\delta A}^5(\vec{k}) \overline{m}_a(\vec{k}) m^i(\vec{k})\,. \end{align} \tag{ A.41 }$$

From this expansion it is evident that $\underset{\smile}{\delta A}^\pm$ and $\underset{\smile}{\delta E}^\pm$ encode the transverse traceless degrees of freedom with corresponding projectors

$$\begin{align} \underset{\smile}{P}{}{^a{}_i{}^j{}_b}(\vec{k}) := \overline{m}^a(\vec{k}) \overline{m}_i(\vec{k}) m^j(\vec{k}) m_b(\vec{k}) + m^a(\vec{k}) m_i(\vec{k}) \overline{m}^j(\vec{k}) \overline{m}_b(\vec{k}) \end{align} \tag{ A.42 }$$

for $\underset{\smile}{\delta E}$ and $\underset{\smile}{P}{}{_a{}^i{}_j{}^b}(\vec{k}) : = \delta_{ac} \delta^{il} \delta_{jm} \delta^{bd}\underset{\smile}{P}{}{^c{}_l{}^m{}_d}(\vec{k})$ for $\underset{\smile}{\delta A}$. Their position space version was given in (2.69) and involves highly-nonlocal terms. Hence we set

$$\begin{align} \underset{\smile}{\delta \mathcal E}{}{^a{}_i}(\vec k) &:= \underset{\smile}{\delta E}^{+}(\vec{k}) m^a(\vec{k})m_i(\vec{k}) + \underset{\smile}{\delta E}^{-}(\vec{k}) \overline{m}^a(\vec{k}) \overline{m}_i(\vec{k}) \end{align} \tag{ A.43 }$$

$$\begin{align} \underset{\smile}{\delta \mathcal A}{}{_a{}^i}(\vec k) &:= \underset{\smile}{\delta A}^{+}(\vec{k}) m_a(\vec{k})m^i(\vec{k}) + \underset{\smile}{\delta A}^{-}(\vec{k}) \overline{m}_a(\vec{k}) \overline{m}^i(\vec{k})\,. \end{align} \tag{ A.44 }$$

The Dirac observables constructed in (2.65) and (2.66) can be expressed in terms of this basis up to linear order and read, where we drop the $\vec k$ dependency for better readability and already set $\delta \xi^j = 0$:

$$\begin{align} (\underset{\smile}{\delta E}{}{^a{}_i})^{GI} & = \underset{\smile}{\delta \mathcal E}{}{^a{}_i}+ \hat{k}^a m_i \left(\underset{\smile}{\delta E}^1 +\frac{1}{||\vec k||} \underset{\smile}{\delta A}^1\right) + \hat{k}^a \overline{m}_i \left(\underset{\smile}{\delta E}^{\overline{1}} -\frac{1}{||\vec k||} \underset{\smile}{\delta A}^{\overline{1}}\right) -m^a \hat{k}_i \frac{1}{||\vec k||} \underset{\smile}{\delta A}^2\nonumber\\[4pt] & \quad +\overline{m}^a \hat{k}_i \frac{1}{||\vec k||} \underset{\smile}{\delta A}^{\overline{2}} + \hat{k}^a\hat{k}_i \underset{\smile}{\delta E}^3 - m^a \overline{m}_i \frac{2}{||\vec k||} \underset{\smile}{\delta A}^4 + \overline{m}^a m_i \frac{2}{||\vec k||} \underset{\smile}{\delta A}^5\nonumber\\[4pt] & \quad +i\kappa\left(\delta_i^a k_c - \delta^a_c \delta_i^b k_b \right) \underset{\smile}{\delta \sigma}^c- i\beta\kappa {{\epsilon}{_c{}^{ab}}}\delta_i^c k_b \underset{\smile}{\delta \tau} +O(\delta^2,\kappa) \end{align} \tag{ A.45 }$$

$$\begin{align} (\underset{\smile}{\delta A}{}{_a{}^i})^{GI}=& \underset{\smile}{\delta \mathcal A}{}{_a{}^i}+\underset{\smile}{\delta A}^2 m_a \hat{k}^i + \underset{\smile}{\delta A}^{\overline{2}} \overline{m}_a \hat{k}^i + \underset{\smile}{\delta A}^4 m_a \overline{m}^i + \underset{\smile}{\delta A}^5 \overline{m}_a m^i+O(\delta^2,\kappa)\,. \end{align} \tag{ A.46 }$$

In this appendix we present further details on the algebra of Dirac observables used in the main text.

A.3.1. Original matter and geometrical variables commute with clocks.

A.3.2. Computing the algebra of the Dirac observables by means of using the Dirac observable of the corresponding Dirac bracket.

In this appendix we will show how the algebra of the Dirac observables is related to the Dirac observable of the Dirac bracket. Here we will discuss the special case and order in perturbation theory that is relevant for the model under consideration in this work and we refer the reader for the proof in the general case to [67, 101]. In the application where we need this relation we will assume that we consider linearised geometrical clocks $\delta T^I(x)$ with $\delta^2T^I(x) = 0$ that only depend on the gravitational degrees of freedom. We then introduce $\delta{\mathcal G}^I(x): = \delta T^I-\tau^I(x)$ that mutually commute with a set of constraints ${\mathcal C}^{^{\prime}}_J(x)$ up to corrections of order δ² that is

$$\begin{align*} M^I_J(x,y):=\{\delta{\mathcal G}^I(x),\delta{\mathcal C}^{^{\prime}}_J(y)\}=\frac{1}{\kappa}\delta^I_J\delta^{(3)}(x,y)+O(\delta^2). \end{align*}$$

Given the definition of $M^I_J(x,y)$ above its inverse is simply given by

$$\begin{align*} (M^{-1})^I_J(x,y):=\kappa\delta^I_J\delta^{(3)}(x,y)+O(\delta^2,\kappa^2) \end{align*}$$

with

$$\begin{align*} \int d^3z M^I_J(x,z)(M^{-1})^J_K(z,y) = \delta^I_K\delta^{(3)}(x,y). \end{align*}$$

Note that we chose to denote the constraints with a prime here because this is the notation used in the main text. This is exactly the setup needed to analyse the algebra of Dirac observables for the model in this work. Since we only perturb the gravitational degrees of freedom but not the matter ones, it is convenient to discuss these two types separately. We will start with the geometric sector and afterwards discuss in detail how the results can be carried over to the algebra of the matter variables. In the framework of field theory the Dirac bracket $\{\cdot,\cdot\}^*$ reads

$$\begin{align} \{f(x), g(y)\}^* & = \{f(x), g(y)\} -\int d^3z^{\prime}\int d^3z^{^{\prime\prime}}\{f(x), {\mathcal C}^{^{\prime}}_L(z^{\prime})\}(M^{-1})^L_M(z^{\prime}, z^{^{\prime\prime}})\{{\mathcal G}^M(z^{^{\prime\prime}}),g(y)\} \nonumber \\[6pt] & \quad +\int d^3z^{\prime}\int d^3z^{^{\prime\prime}}\{g(y), {\mathcal C}^{^{\prime}}_L(z^{\prime})\}(M^{-1})^L_M(z^{\prime}, z^{^{\prime\prime}})\{{\mathcal G}^M(z^{^{\prime\prime}}),f(x)\} \nonumber \\ & =\{f(x), g(y)\}-\int d^3z\{f(x), {\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{T^M(z),g(y)\} \nonumber \\[6pt] & \quad +\int d^3z\{g(y), {\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{T^M(z),f(x)\}, \end{align} \tag{ A.47 }$$

where we used the explicit form of $(M^{-1})^J_K(z,y)$ together with the fact that inside the Poisson bracket we can replace ${\mathcal G}^M$ by T^M because they differ only by a phase space independent function. Quantities involving gravitational variables are perturbed in the model in this article. Let us denote their linear and second-order perturbations by $\delta f$ and $\delta^2 f$ respectively. Then we consider the perturbation of the Dirac bracket up to linear order that is given by

$$\begin{align} \{f(x), g(y)\}^* & = \{\delta f(x),\delta g(y)\}-\kappa\int d^3z\{\delta f(x), \delta {\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta g(y)\}\nonumber\\[4pt] & \quad +\kappa\int d^3z\{\delta g(y), \delta {\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta f(x)\} \nonumber\\[4pt] & \quad \times \{\delta f(x),\delta^2 g(y)\}+ \{\delta^2 f(x),\delta g(y)\} \nonumber \\[4pt] & \quad -\kappa\int d^3z\{\delta^2 f(x), \delta {\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta g(y)\} \nonumber\\[4pt] & \quad -\kappa\int d^3z\{\delta f(x), \delta^2{\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta g(y)\} \nonumber \\[4pt] & \quad -\kappa\int d^3z\{\delta f(x), \delta{\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta^2 g(y)\} \nonumber \\[4pt] & \quad +\kappa\int d^3z\{\delta^2 g(y), \delta {\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta f(x)\} \nonumber \\[4pt] & \quad +\kappa\int d^3z\{\delta g(y), \delta^2 {\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta f(x)\} \nonumber \\[4pt] & \quad +\kappa\int d^3z\{\delta g(y), \delta {\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta^2 f(x)\} +O(\delta^2,\kappa^2). \end{align} \tag{ A.48 }$$

Now for elementary phase space variables we have $\delta^2 f = 0$ and $\delta^2 g = 0$ then the Dirac bracket simplifies to

$$\begin{align} \{f(x), g(y)\}^* & = \{\delta f(x),\delta g(y)\} -\kappa\int d^3z\{\delta f(x), \delta {\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta g(y)\} \\[4pt] &\quad +\kappa\int d^3z\{\delta g(y), \delta {\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta f(x)\} \nonumber \\[4pt] &\quad -\kappa\int d^3z\{\delta f(x), \delta^2{\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta g(y)\} \nonumber \end{align} \tag{ A.49 }$$

$$\begin{align} &\quad +\kappa\int d^3z\{\delta g(y), \delta^2 {\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta f(x)\} +O(\delta^2,\kappa^2). \end{align} \tag{ A.50 }$$

The Dirac observable associated with $\delta f$ is up to second order given by

$$\begin{align} O_{\delta f,\{T\}} & = \delta f -\kappa\int d^3z \delta {\mathcal G}^{K}(z)\{\delta f,\delta{\mathcal C}^{^{\prime}}_K(z)\} -\kappa\int d^3z \delta {\mathcal G}^{K}(z)\{\delta f,\delta^2{\mathcal C}^{^{\prime}}_K(z)\} \nonumber \\ & \quad -\frac{\kappa^2}{2}\int d^3z\int d^3z^{^{\prime}} \delta {\mathcal G}^{K}(z) \delta {\mathcal G}^{L}(z^{^{\prime}})\{\{\delta f,\delta^2{\mathcal C}^{^{\prime}}(z)\},\delta{\mathcal C}^{^{\prime}}_M(z^{^{\prime}})\}+O(\delta^3,\kappa^2) \nonumber \\ & =: O^{(1)}_{\delta f,\{T\}}+ O^{(2)}_{\delta f,\{T\}}+O(\delta^3,\kappa^2), \end{align} \tag{ A.51 }$$

where, as in the main text, we include all terms of order δ and δ² in $O^{(1)}_{\delta f,\{T\}}$ and $O^{(2)}_{\delta f,\{T\}}$ respectively. Our aim is to show the following: For linearised quantities $\delta f,\delta g$ with $\delta^2f = 0$ and $\delta^2g = 0$ with a set of constraints ${\mathcal C}^{^{\prime}}_I$ and reference fields ${\mathcal G}^I$ that satisfy $\{{\mathcal G}^I(x),{\mathcal C}^{^{\prime}}_J(y)\} = \delta^I_J\delta^{(3)}(x,y)+O(\delta^2,\kappa^2)$ we have

$$\begin{align} \{O_{\delta f,\{T\}}, O_{\delta g,\{T\}}\} = O_{\{\delta f,\delta g\}^*,\{T\}} +O(\delta^2,\kappa^2), \end{align} \tag{ A.52 }$$

where the crucial property for our later application is that we have a strong equality here and not a weak one that is present in the general case [67, 101]. For the purpose of showing the equality in (A.53) we introduce the following notation:

$$\begin{align} \{f,g\}^*\big|_{\delta^0} & := \{\delta f(x),\delta g(y)\} -\kappa\int d^3z\{\delta f(x), \delta {\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta g(y)\}\nonumber\\[6pt] & \quad +\kappa\int d^3z\{\delta g(y), \delta {\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta f(x)\} \nonumber \\[6pt] \{f,g\}^*\big|_{\delta^1} & :=-\kappa\int d^3z\{\delta f(x), \delta^2{\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta g(y)\} \nonumber \\[6pt] & \quad +\kappa\int d^3z\{\delta g(y), \delta^2 {\mathcal C}^{^{\prime}}_L(z)\}\delta^L_M\{\delta T^M(z),\delta f(x)\}\, , \end{align} \tag{ A.53 }$$

where we again used that $\delta^2f = \delta^2g = 0$ and that we can replace $\delta{\mathcal G}^M$ by $\delta T^M$ inside the Poisson brackets. Given this notation we can write $O_{\{\delta f,\delta g\}^*,\{T\}}$ up to linear order as

$$\begin{align} O_{\{\delta f,\delta g\}^*,\{T\}} & = \{f,g\}^*\big|_{\delta^0}+\{f,g\}^*\big|_{\delta^1} -\kappa\int d^3z \delta{\mathcal G}^K(z)\{\{f,g\}^*\big|_{\delta^1}, \delta{\mathcal C}^{^{\prime}}_K(z)\} + O(\delta^2,\kappa^2)\nonumber \\[6pt] & =\{f,g\}^*\big|_{\delta^0}+\{f,g\}^*\big|_{\delta^1} -\kappa^2\int d^3z\int d^3z^{^{\prime}} \delta{\mathcal G}^L(z)\nonumber \\[6pt] & \quad\times \Big(-\delta^K_M\{\delta T^M(z^{^{\prime}}),\delta g\}\{\{\delta f, \delta^2{\mathcal C}^{^{\prime}}_K(z^{^{\prime}}) \},\delta {\mathcal C}^{^{\prime}}_L(z)\} \nonumber\\[6pt] & \qquad\quad -\delta^K_M\{\delta T^M(z^{^{\prime}}),\delta f\}\{\{\delta g, \delta^2{\mathcal C}^{^{\prime}}_K(z^{^{\prime}}) \},\delta {\mathcal C}^{^{\prime}}_L(z)\}\Big)+O(\delta^2,\kappa^2). \end{align} \tag{ A.54 }$$

Hence, we need to show that $\{O^{(1)}_{\delta f,\{T\}}+ O^{(2)}_{\delta f,\{T\}},O^{(1)}_{\delta g,\{T\}}+ O^{(2)}_{\delta g,\{T\}}\}$ agrees with the right hand side of (A.54). To confirm this we insert the Dirac observables of $\delta f$ and $\delta g$ up to second order in (A.51) into the Poisson bracket and collect all terms in the individual orders. In order δ⁰ we obtain

$$\begin{align} & \delta^0 : \{\delta f(x), \delta g(y)\} -\kappa\int d^3z \{\delta f(x),\delta T^L(z)\}\{\delta g(y), \delta {\mathcal C}^{^{\prime}}_L(z)\}\nonumber\\[6pt] & \qquad -\kappa\int d^3z \{\delta T^L(z),\delta g(y)\} \{\delta f(x), \delta {\mathcal C}^{^{\prime}}_L(z)\} \nonumber \\[6pt] & \quad = \{f(x),g(y)\}^*\big|_{\delta^0}, \end{align} \tag{ A.55 }$$

where we used the antisymmetry of the Poisson bracket as well as again used that we can replace $\delta{\mathcal G}^M$ by $\delta T^M$ inside a given Poisson bracket. For the linear order in δ the result is given by

$$\begin{align} \delta^1 : & \underbrace{-\kappa\int d^3z{\mathcal G}^M(z)\{\delta f(x),\{\delta g(y), \delta^2{\mathcal C}^{^{\prime}}_M(z)\}\}}_{=:A_1} \underbrace{-\kappa\int d^3z \{\delta f(x),\delta T^M(z)\}\{\delta g(y),\delta^2{\mathcal C}^{^{\prime}}_M(z)\}}_{=:B_1}\nonumber \\[4pt] &\underbrace{-\kappa\int d^3z {\mathcal G}^K(z)\{\{\delta f(x), \delta^2{\mathcal C}^{^{\prime}}_K(z)\},\delta g(y)\}}_{=:A_2} \underbrace{-\kappa\int d^3z \{\delta T^K(z),\delta g(y)\}\{\delta f(x),\delta^2{\mathcal C}^{^{\prime}}_K(z)\}}_{=:B_2}\nonumber \\[4pt] & \underbrace{-\kappa^2\int d^3z \int d^3z^{^{\prime}}\{\delta f(x), \delta T^{M}(z)\}\delta{\mathcal G}^N(z^{^{\prime}})\{\{\delta g(y) ,\delta^2{\mathcal C}^{^{\prime}}_M(z)\}, \delta {\mathcal C}^{^{\prime}}_N(z^{^{\prime}})\}}_{:=C_1}\nonumber \\[4pt] &\underbrace{-\kappa^2\int d^3z \int d^3z^{^{\prime}}\{\delta T^{K}(z),\delta g(y)\}\delta{\mathcal G}^L(z^{^{\prime}})\{\{\delta f(x) ,\delta^2{\mathcal C}^{^{\prime}}_K(z)\}, \delta {\mathcal C}^{^{\prime}}_L(z^{^{\prime}})\}}_{:=C_2}, \end{align} \tag{ A.56 }$$

where we again used that we can replace $\delta{\mathcal G}^M$ by $\delta T^M$ inside a given Poisson bracket.

Now let us discuss the individual contributions separately. First we show that $A_1+A_2 = 0$. Using the Jacobi identity for A₁ we obtain

$$\begin{align} A_1 & = -\kappa\int d^3z{\mathcal G}^M(z)\{\delta f(x),\{\delta g(y), \delta^2{\mathcal C}^{^{\prime}}_M(z)\}\}\\[3pt] & = \kappa\int d^3z {\mathcal G}^M(z)\Big( \{\delta g(y),\{\delta^2{\mathcal C}^{^{\prime}}_M(z), \delta f(x)\}\} + \underbrace{\{\delta^2{\mathcal C}^{^{\prime}}_M(z),\underbrace{\{\delta f(x),\delta g(y)\}}_{\sim \delta^0}\} }_{=0}\Big)\nonumber \end{align} \tag{ A.57 }$$

$$\begin{align} & =\kappa\int d^3z{\mathcal G}^M(z) \{\delta g(y),\{\delta^2{\mathcal C}^{^{\prime}}_M(z), \delta f(x)\}\} \nonumber \\[3pt] & = \kappa\int d^3z {\mathcal G}^M(z) \{\{\delta f(x), \delta^2{\mathcal C}^{^{\prime}}_M(z)\},\delta g(y)\}\}\nonumber \\[3pt] & = -A_2. \end{align} \tag{ A.58 }$$

Comparing $B_1+B_2$ with $\{f,g\}^*\big|_{\delta^1}$ in (A.53) and using the antisymmetry of the Poisson bracket we realise that we have

$$\begin{align} B_1+B_2 = \{f,g\}^*\big|_{\delta^1}. \end{align} \tag{ A.59 }$$

Comparing $C_1+C_2$ with the last two terms in (A.54) and applying a suitable relabelling of the integration variables it turns out that these two expressions are exactly identical. Thus, collecting all intermediate results we have shown

$$\begin{align*} O_{\{\delta f,\delta g\}^*,\{T\}} = A_1+A_2+B_1+B_2+C_1+C_2 +O(\delta^2, \kappa^2) =\{O_{\delta f,\{T\}}, O_{\delta g,\{T\}}\} +O(\delta^2, \kappa^2), \end{align*}$$

which is exactly the equality in (A.52) we wanted to prove.

Now let us discuss the case of elementary variables from the matter sector. Because we use geometrical clocks if we consider for $\delta f$ and $\delta g$ both matter variables, then since these will commute with all geometrical clocks we immediately have $\{\delta f, \delta g\}^* = \{\delta f,\delta g\}$ so that we can work with the Poisson bracket instead of Dirac bracket in (A.52) and as discussed in the main text the observable algebra drastically simplifies. The other case is if we consider the algebra of observables of one geometric and one matter quantity and without loss of the generality let us assume that $\delta f$ contains geometric and $\delta g$ matter variables. Then in general the Dirac bracket can differ from the Poisson bracket by those terms where $\delta f$ is involved in a Poisson bracket with the geometric clocks. In this case the proof just presented carries over if we just replace $\delta g$ by g for the following reason: The Poisson bracket $\{g(y), \delta {\mathcal C}^{\prime}_L(z)\}$ for g(y) a function of elementary matter variables is of order δ⁰ and hence also contributes in this case to the δ⁰-order of the Dirac bracket. The second relevant Poisson bracket $\{g(y), \delta^2{\mathcal C}^{\prime}_L(z)\}$ is as before of order δ¹ because in the order of perturbation theory we consider, the contributions to this Poisson bracket come from terms that involve the gravitational perturbation linearly and the matter variables in quadratic order so that the final result will still be of order δ¹. Therefore, we can also apply the equality shown in (A.52) to matter variables by simply replacing $\delta f$ by f and $\delta g$ by g if we restrict $f,g$ to be elementary phase space variables from the set $\phi(x),\pi(y)$. Further because the Poisson brackets of the matter variables do not involve a factor $\frac{1}{\kappa}$ in contrast to the gravitational variables additional factors κ can arise when the Dirac bracket with matter variables is considered. The application we need in the main text is even a more special case so we will not discuss these terms more in detail here since they will not be needed in any further computation. The main motivation for presenting this proof here in detail was to understand under which assumptions regarding perturbation theory we can ensure that the Dirac bracket agrees with the Poisson bracket and how, if this is not the case, the explicit form of the Dirac observable of the Dirac bracket looks like in perturbation theory.

A.3.3. Poisson algebra of linearised Dirac observables.

In this appendix we will compute the Poisson algebra of the linearised Dirac observables explicitly and show that the zeroth order contributions are given by the standard canonical Poisson brackets, a result that we use in the main text of this work. In the notation of the main text for all observables $O_{f,\{T\}} = O^{(1)}_{f,\{T\}}+O^{(2)}_{f,\{T\}}+O(\delta^3,\kappa^3)$ we will neglect the contributions coming from $O^{(2)}_{f,\{T\}}+O(\delta^3,\kappa^3)$. This means that the explicit computations presented in this appendix here do only consider the zeroth order result of the corresponding Poisson algebra. This causes no problem since we also know this as well as the linear contribution to the Poisson algebra already form the relation to the Dirac observable of the corresponding Dirac bracket as discussed in the main text. The computation presented here should be rather understood as double checking our results for the explicit form of the Dirac observables that would yield incorrect results for the algebra in zeroth order if they were not computed correctly. For the gauge invariant geometrical degrees of freedom, one can use the following way to obtain their algebra:

$$\begin{align} \{{{\delta A}{_a{}^i}}(\vec x,t),{\delta E}{^b{}_j}(\vec y,t)\} = \frac{\beta}{\kappa} \delta_j^i\delta_a^b \delta(\vec x-\vec y) \implies \{\underset{\smile}{\delta A}{}{_a{}^i}(\vec k,t),\underset{\smile}{\delta E}{}{^b{}_j}(\vec p,t)\} = \frac{\beta}{\kappa} \delta_j^i\delta_a^b \delta(\vec k+\vec p) \end{align} \tag{ A.60 }$$

and hence for $r,u\in\{\pm\}$:

$$\begin{align} \{\underset{\smile}{\delta A}^r(\vec k,t) ,\underset{\smile}{\delta E}^u(\vec p,t) \} & = \overline{m}^a(r\vec{k})\overline{m}_i(r\vec{k})\overline{m}_b(u\vec{p})\overline{m}^j(u\vec{p}) \; \{\underset{\smile}{\delta A}{}{_a{}^i}(\vec k,t), \underset{\smile}{\delta E}{}{^b{}_j}(\vec p,t) \}\nonumber\\ & = \frac{\beta}{\kappa} \overline{m}^a(r\vec{k})\overline{m}_i(r\vec{k}){m}_a(u\vec{k}){m}^b(u\vec{k}) \delta(\vec{k}+\vec{p}) \nonumber\\ & = \delta^{ru} \frac{\beta}{\kappa} \delta(\vec{k}+\vec{p})\,. \end{align} \tag{ A.61 }$$

Analogously one can show that $\{\underset{\smile}{\delta A}^r(\vec k,t) ,\underset{\smile}{\delta A}^u(\vec p,t) \} = \{\underset{\smile}{\delta E}^r(\vec k,t) ,\underset{\smile}{\delta E}^u(\vec p,t) \} = 0$. From this and the fact that $O^{(1)}_{\delta \mathcal{A}{}{_a{}^i},\{\delta T\}}(\vec x,t) = O^{(1)}_{P\delta A{}{_a{}^i},\{\delta T\}}(\vec x,t) = PO^{(1)}_{\delta A{}{_a{}^i},\{\delta T\}}(\vec x,t) = P{A}{_a{}^i}(\vec x,t)$ and similarly for $\delta E$, follows that

$$\begin{align} \{O^{(1)}_{\delta \mathcal{A}{}{_a{}^i},\{\delta T\}}(\vec x,t),O^{(1)}_{\delta \mathcal{E}{}{^b{}_j},\{\delta T\}}(\vec y,t)\} &= \int \frac{d^3k d^3p}{(2\pi)^3} e^{i\vec k\vec x + i\vec p \vec y} \underset{\smile}{P}{}{_a{}^i{}_l{}^c}(\vec k) \underset{\smile}{P}{}{^b{}_j{}^m{}_d}(\vec p) \{\underset{\smile}{\delta A}{}{_c{}^l}(\vec k,t),\underset{\smile}{\delta E}{}{^d{}_m}(\vec p,t)\}\nonumber\\[6pt] &= \frac{\beta}{\kappa}\int \frac{d^3k}{(2\pi)^3} e^{i\vec k(\vec x -\vec y)} \underset{\smile}{P}{}{_a{}^i{}_l{}^c}(\vec k) \underset{\smile}{P}{}{^b{}_j{}^l{}_c}(-\vec k)\nonumber\\[6pt] &=\frac{\beta}{\kappa}\int \frac{d^3k}{(2\pi)^3} e^{i\vec k(\vec x -\vec y)} \underset{\smile}{P}{}{_a{}^i{}_j{}^b}(\vec k)\nonumber\\[6pt] &= \frac{\beta}{\kappa} {{P}{_a{}^i{}_j{}^b}}\delta(\vec x-\vec y)\,, \end{align} \tag{ A.62 }$$

where we used in the fourth line that $\underset{\smile}{P}{}{^b{}_j{}^l{}_c}(-\vec k) = \underset{\smile}{P}{}{^b{}_j{}^l{}_c}(\vec k) = \underset{\smile}{P}{}{^l{}_c{}^b{}_j}(\vec k)$, as well as

$$\begin{align} \{O^{(1)}_{\delta \mathcal{A}{}{_a{}^i},\{\delta T\}}(\vec x,t),O^{(1)}_{\delta \mathcal{A}{}{_b{}^j},\{\delta T\}}(\vec y,t)\} &= 0 \end{align} \tag{ A.63 }$$

$$\begin{align} \{O^{(1)}_{\delta \mathcal{E}{}{^a{}_i},\{\delta T\}}(\vec x,t),O^{(1)}_{\delta \mathcal{E}{}{^b{}_j},\{\delta T\}}(\vec y,t)\} &= 0\,. \end{align} \tag{ A.64 }$$

For the matter variables we obtain, dropping terms in $O(\delta^2,\kappa^2)$ and higher:

$$\begin{align} \{O^{(1)}_{\phi,\{\delta T\}}(\vec x,t), O^{(1)}_{\phi,\{\delta T\}}(\vec y,t) \} &= -\kappa^2 (\delta \mathcal{G}^c(\vec x,t)) \delta^3(\vec x-\vec y) +\kappa^2 (\delta\mathcal{G}^c(\vec y,t)) \delta^3(\vec x-\vec y)=0\nonumber\\ \end{align} \tag{ A.65 }$$

$$\begin{align} \{O^{(1)}_{\pi,\{\delta T\}}(\vec x,t), O^{(1)}_{\pi,\{\delta T\}}(\vec y,t)\} &= \kappa^2 ((\partial_a^y \:\delta \mathcal{G}(\vec{y},t) + \partial_a^x \:\delta \mathcal{G}(\vec x,t)) \partial_y^a \delta^3(\vec x-\vec y) \nonumber\\[6pt] & \quad + (\delta \mathcal{G}(\vec y,t)-\delta \mathcal{G}(\vec x,t)) \Delta \delta^3(\vec x-\vec y)] \end{align} \tag{ A.66 }$$

$$\begin{align} \{O^{(1)}_{\phi,\{\delta T\}}(\vec x,t),O^{(1)}_{\pi,\{\delta T\}}(\vec y,t)\} & = \delta^3(\vec x-\vec y) -\kappa^2 \partial_a^y[\delta\mathcal{G}^a(\vec y,t) \delta^3(\vec x-\vec y)] \nonumber\\[6pt] & \quad -\kappa^2 \;\delta \mathcal{G}^a(\vec x,t) \partial_a^x\delta^3(\vec x-\vec y)\,. \end{align} \tag{ A.67 }$$

Due to the presence of derivatives acting on the delta distributions, it is not immediately evident that all additional terms (which are of $O(\kappa)$ due to the presence of a factor κ⁻¹ in the reference fields) vanish. To see this, it is convenient to consider the smeared version of the linearised observables, then it turns out that:

$$\begin{align} \int d^3x \int d^3y \;f(\vec x) g(\vec y) \{O^{(1)}_{\pi,\{\delta T\}}(\vec x,t), O^{(1)}_{\pi,\{\delta T\}}(\vec y,t)\} &= 0 \end{align} \tag{ A.68 }$$

$$\begin{align} \int d^3x \int d^3y \;f(\vec x) g(\vec y) \{O^{(1)}_{\phi,\{\delta T\}}(\vec x,t),O^{(1)}_{\pi,\{\delta T\}}(\vec y,t)\} &= \int d^3x \; f(\vec x) g(\vec x)\,. \end{align} \tag{ A.69 }$$

Hence we indeed end up with the desired algebra also for the matter observables:

$$\begin{align} \{O^{(1)}_{\phi,\{\delta T\}}(\vec x,t), O^{(1)}_{\phi,\{\delta T\}}(\vec y,t) \} & =0 \nonumber\\[6pt] \{O^{(1)}_{\pi,\{\delta T\}}(\vec x,t), O^{(1)}_{\pi,\{\delta T\}}(\vec y,t) \} & =0 \nonumber \\[6pt] \{O^{(1)}_{\phi,\{\delta T\}}(\vec x,t), O^{(1)}_{\pi,\{\delta T\}}(\vec y,t) \} & =\delta^3(\vec x-\vec y)\,. \end{align} \tag{ A.70 }$$

As $O^{(1)}_{\delta \mathcal{A}{}{_a{}^i},\{\delta T\}}$ and $O^{(1)}_{\delta \mathcal{E}{}{^a{}_i},\{\delta T\}}$ commute with the clocks and do not depend on the original matter fields, all Poisson brackets among the $(O^{(1)}_{\delta \mathcal{A},\{\delta T\}},O^{(1)}_{\delta \mathcal{E},\{\delta T\}})$ vanish up to $O(\delta^2,\kappa^2)$ and the $(O^{(1)}_{\phi,\{\delta T\}},O^{(1)}_{\pi,\{\delta T\}})$ vanish up to $O(\kappa^2)$.

A.3.4. Comparison with the Dirac observables used by Anastopoulos and Hu.

The model in [35] also considers gauge invariant quantities for the matter variables in the ADM framework and we would like to compare their quantities to the constructed Dirac observables in this work. There, an extension of the basic gauge variant matter variables φ and π is proposed that is required to commute with the linearised constraints. These observables take the form

$$\begin{align*} \widetilde{\phi}(\vec x,t) = \phi(\vec x-\vec q,t-\tau) \quad \widetilde{\pi}(\vec x,t) = \pi(\vec x-\vec q,t-\tau) -\partial_a q^a\,, \end{align*}$$

where $\vec q$ and τ, similarly to the geometrical clocks used in our work, consist of purely geometrical quantities in ADM variables. For these as well as the matter fields in [35] the following transformation behaviour under the Hamiltonian and spatial diffeomorphism constraint respectively is used ⁸ :

$$\begin{align*} \begin{array}{r@{}ll} \text{Under the Hamilton constraint}: & \phi \rightarrow \phi + \lambda \frac{\delta H_0}{\delta \pi} & \pi \rightarrow \pi -\lambda \frac{\delta H_0}{\delta \phi}\\ & \vec{q} \rightarrow \vec{q} & \kern1.7pt\tau \rightarrow \tau + \lambda\\ \text{Under the spatial diffeomorphism constraint}: & \phi \rightarrow \phi + \kappa \lambda^a \partial_a\phi & \pi \rightarrow \pi +\kappa \partial_a (\lambda^a \pi)\\ & \kern-1.5pt q^i \rightarrow q^i+\lambda^i & \kern1.7pt\tau \rightarrow \tau\,, \end{array} \end{align*}$$

where $H_0 = \int d^3x \; \epsilon(\vec x,t)$. Starting with the transformations under the Hamilton constraint, which is just $\epsilon(\vec x,t)$ in the matter part, we obtain

$$\begin{align*} \{\phi(\vec x), \int d^3y \; \lambda(\vec y) \epsilon(\vec y)\} = \lambda(\vec x) \pi(\vec x) = \lambda(\vec x) \frac{\delta H_0}{\delta \pi}(\vec x) \end{align*}$$

which hence indeed confirms their proposed transformation of the matter field, its momentum transforms differently reflecting the fact that it is a scalar density of weight one:

$$\begin{align*} \left\{\pi(\vec x), \int d^3y \; \lambda(\vec y) \epsilon(\vec y)\right\} & = \partial_a(\lambda(\vec x) \partial^a \phi(\vec x)) - \lambda(\vec x) V^{^{\prime}}(\phi(\vec x)) = -\lambda(\vec x) \frac{\delta H_0}{\delta\phi}(\vec x) \nonumber\\ & \quad + (\partial_a\lambda(\vec x)) \; \partial^a \phi(\vec x), \end{align*}$$

where we assumed that the scalar field obeys a Klein–Gordon equation and that the energy density possesses an arbitrary potential $V(\phi)$. The transformation behaviour used in [35] under the Hamiltonian constraint only agrees with the above expression if λ is constant that in general cannot be assumed and hence the second term coming from the Poisson bracket needs to be included in the transformation of $\pi(\vec x)$ and to our understanding has been omitted in [35].

Given the correct general transformation behaviours, we want to check whether the Dirac observables used in [35] $\tilde{\phi}$ and $\widetilde{\pi}$ are indeed gauge invariant. First it is worth noticing that due to the fact that we are working in a linearised theory, to our understanding the quantity $\widetilde{\phi}(\vec x,t)$ has to be understood as a Taylor expansion in the modifications truncated after linear order. In general this way of constructing observables looks like the canonical version of the strategy followed in [100]. To show this, we will add an additional κ in front of the modifications and also include it into the constraints. In [35] this was done for the spatial diffeomorphism constraint only, as its action on the matter fields indeed yields a correction term of order κ. As shown in the discussion on the Post-Minkowski approximation scheme in 2.3, for consistency such a factor is also required in the Hamiltonian constraint. Therefore, we obtain

$$\begin{align*} \widetilde{\phi}(\vec x,t)& = \phi(\vec x - \kappa \vec q(\vec x,t), t-\kappa\tau(\vec x,t)) = \phi(\vec x,t) - \kappa q^a(\vec x,t) \partial_a\phi(\vec x,t) -\kappa \tau(\vec x,t) \dot{\phi}(\vec x,t)\\ \widetilde{\pi}(\vec x,t)& = \pi(\vec x - \kappa \vec q(\vec x,t), t-\kappa\tau(\vec x,t))-\partial_a q^a(\vec x,t) = \pi(\vec x,t) - \kappa q^a \partial_a\pi(\vec x,t) -\kappa \tau \dot{\pi}(\vec x,t)\\ & \quad - \partial_a q^a(\vec x,t)\,. \end{align*}$$

Now under the Hamilton constraint $\widetilde{\phi}$ transforms as:

$$\begin{align*} \widetilde{\phi}(\vec x,t) \rightarrow\widetilde{\phi}(\vec x,t) + \kappa(\lambda\pi)(\vec x,t) - \kappa (\lambda\dot{\phi})(\vec x,t) = \widetilde{\phi}(\vec x,t) \,, \end{align*}$$

where we neglected terms of higher order in κ and see that indeed $\widetilde{\phi}$ remains invariant under the transformations induced by the Hamilton constraint. For $\widetilde{\pi}$ one obtains:

$$\begin{align*} \widetilde{\pi}(\vec x,t) \rightarrow \widetilde{\pi}(\vec x,t)+ \kappa(\lambda \dot{\pi})(\vec x,t) - \kappa (\lambda \dot{\pi})(\vec x,t) =\widetilde{\pi}(\vec x,t)\,. \end{align*}$$

Thus, working with the correct transformation behaviour of π under the Hamiltonian constraint, also $\widetilde{\pi}$ remains invariant. For the spatial diffeomorphism constraint a quick calculation shows that its action on the matter fields is indeed given by the terms stated above and in [35]. Under its action we obtain:

$$\begin{align*} \widetilde{\phi} &\rightarrow\widetilde{\phi} +\kappa \lambda^a \partial_a \phi -\kappa \lambda^a \partial_a \phi = \widetilde{\phi}\\ \widetilde{\pi} &\rightarrow\widetilde{\pi} -\kappa \lambda^a \partial_a\pi -\partial_a\lambda^a \neq \widetilde{\pi}\, \end{align*}$$

which demonstrates that the $\tilde{\pi}$ proposed in [35] is no (linearised) Dirac observable with respect to the spatial diffeomorphism constraint. Thus, actual linearised observables for $\phi(\vec x,t),\pi(\vec x,t)$ in [35] would be given by

$$\begin{align*} \widetilde{\phi}(\vec x,t)& = \phi(\vec x,t) - \kappa q^a(\vec x,t) \partial_a\phi(\vec x,t) -\kappa \tau(\vec x,t) \dot{\phi}(\vec x,t)\\ \widetilde{\pi}(\vec x,t)& = \pi(\vec x,t) - \kappa q^a \partial_a(q^a\pi)(\vec x,t) -\kappa \tau \dot{\pi}(\vec x,t)\,, \end{align*}$$

where it is crucial that the partial derivative in the second term in $\widetilde{\pi}$ acts on both $q^a(\vec x,t)$ and $\pi(\vec x,t)$. Compared to our framework, the basic idea is similar to extend the matter fields with geometrical variables to obtain gauge invariant observables. Hence, the quantities q^a and τ play the role of geometrical clocks in our setup. In our approach we in addition require that q^a and τ mutually commute which requires to use an explicit choice of geometrical clocks that satisfies this additional property.

In this appendix the aim is to find an expression of the original phase space variables in terms of the new ones which consist of the physical degrees of freedom, the constraints and the clocks. For this, we first transform the clocks and constraints into momentum space using the basis introduced above in appendix A.2 and end up with the following expressions:

$$\begin{align} \underset{\smile}{\delta \mathcal{G}} &= \frac{1}{2\beta \kappa||\vec{k}||} (\underset{\smile}{\delta E}^5 - \underset{\smile}{\delta E}^4) - \frac{1}{\beta \kappa||\vec{k}||^2} (\underset{\smile}{\delta A}^3+ \underset{\smile}{\delta A}^4+\underset{\smile}{\delta A}^5)-\underset{\smile}{\delta \tau} \end{align} \tag{ A.71 }$$

$$\begin{align} \underset{\smile}{\delta C} &= \kappa ||\vec{k}|| (\underset{\smile}{\delta A}^4 - \underset{\smile}{\delta A}^5 ) +\kappa \underset{\smile}{\epsilon} \end{align} \tag{ A.72 }$$

$$\begin{align} m_a \underset{\smile}{\delta \mathcal{G}}^a &= \frac{i}{\kappa||\vec{k}||} (\underset{\smile}{\delta E}^{\overline{1}} + \underset{\smile}{\delta E}^{\overline{2}} ) - \frac{2\beta}{\kappa^2 ||\vec{k}||^2} m^i \underset{\smile}{\delta G}{}_i - m_a \underset{\smile}{\delta \sigma}^a \end{align} \tag{ A.73 }$$

$$\begin{align} \overline{m}_a \underset{\smile}{\delta \mathcal{G}}^a &= \frac{i}{\kappa||\vec{k}||} (\underset{\smile}{\delta E}^{1} + \underset{\smile}{\delta E}^{2} ) - \frac{2\beta}{\kappa^2||\vec{k}||^2} \overline{m}^i \underset{\smile}{\delta G}{}_i- \overline{m}_a \underset{\smile}{\delta \sigma}^a \end{align} \tag{ A.74 }$$

$$\begin{align} \hat{k}_a \underset{\smile}{\delta \mathcal{G}}^a &= \frac{i}{2\kappa||\vec{k}||} (3\underset{\smile}{\delta E}^{3} - \underset{\smile}{\delta E}^{4} - \underset{\smile}{\delta E}^5 ) - \frac{3\beta}{\kappa^2 ||\vec{k}||^2} \hat{k}^i \underset{\smile}{\delta G}{}_i +\frac{i}{2\kappa ||\vec{k}||^3} \left(\frac{1}{\kappa}\underset{\smile}{\delta C}-\underset{\smile}{\epsilon}\right)- \hat{k}_a \underset{\smile}{\delta \sigma}^a\nonumber\\ \end{align} \tag{ A.75 }$$

$$\begin{align} m^a \underset{\smile}{\delta C}{}_a &= -\frac{i\kappa||\vec{k}||}{\beta} \underset{\smile}{\delta A}^{\overline{2}} + \kappa m^a \underset{\smile}{p}{}_a \end{align} \tag{ A.76 }$$

$$\begin{align} \overline{m}^a \underset{\smile}{\delta C}{}_a &= -\frac{i\kappa||\vec{k}||}{\beta} \underset{\smile}{\delta A}^{2} + \kappa \overline{m}^a \underset{\smile}{p}{}_a \end{align} \tag{ A.77 }$$

$$\begin{align} \hat{k}^a \underset{\smile}{\delta C}{}_a &= \frac{i\kappa||\vec{k}||}{\beta} (\underset{\smile}{\delta A}^{4} + \underset{\smile}{\delta A}^5 ) + \kappa \hat{k}^a \underset{\smile}{p}{}_a \end{align} \tag{ A.78 }$$

$$\begin{align} m^i \underset{\smile}{\delta G}{}_i &= \frac{\kappa i}{2\beta} (||\vec{k}|| \underset{\smile}{\delta E}^{\overline{1}} - \underset{\smile}{\delta A}^{\overline{1}} + \underset{\smile}{\delta A}^{\overline{2}} ) \end{align} \tag{ A.79 }$$

$$\begin{align} \overline{m}^i \underset{\smile}{\delta G}{}_i &= \frac{\kappa i}{2\beta} (||\vec{k}|| \underset{\smile}{\delta E}^{1} - \underset{\smile}{\delta A}^{2} + \underset{\smile}{\delta A}^{1} ) \end{align} \tag{ A.80 }$$

$$\begin{align} \hat{k}^i \underset{\smile}{\delta G}{}_i &= \frac{\kappa i}{2\beta} (||\vec{k}|| \underset{\smile}{\delta E}^{3} - \underset{\smile}{\delta A}^{5} + \underset{\smile}{\delta A}^{4} ) \end{align} \tag{ A.81 }$$

$$\begin{align} m_i \underset{\smile}{\delta \mathcal{G}}^i &= \frac{2i}{\kappa||\vec{k}||} \underset{\smile}{\delta A}^{\overline{1}} - m_i \underset{\smile}{\delta \xi}^i \end{align} \tag{ A.82 }$$

$$\begin{align} \overline{m}_i \underset{\smile}{\delta \mathcal{G}}^i &= \frac{2i}{\kappa||\vec{k}||} \underset{\smile}{\delta A}^{1} -\overline{m}_i \underset{\smile}{\delta \xi}^i \end{align} \tag{ A.83 }$$

$$\begin{align} \hat{k}_i \underset{\smile}{\delta \mathcal{G}}^i &= \frac{2i}{\kappa||\vec{k}||} \underset{\smile}{\delta A}^{3}-\hat{k}_i \underset{\smile}{\delta \xi}^i\,. \end{align} \tag{ A.84 }$$

In total, these are 14 equations that can be solved uniquely for the 14 old variables. The solution then is

$$\begin{align} \underset{\smile}{\delta A}^1 &= -\frac{i\kappa||\vec{k}||}{2} \overline{m}_i \left[\underset{\smile}{\delta \mathcal{G}}^i+ \underset{\smile}{\delta \xi}^i\right] \end{align} \tag{ A.85 }$$

$$\begin{align} \underset{\smile}{\delta A}^{\overline{1}} &= -\frac{i\kappa||\vec{k}||}{2} m_i \left[\underset{\smile}{\delta\mathcal{G}}^i+ \underset{\smile}{\delta \xi}^i\right] \end{align} \tag{ A.86 }$$

$$\begin{align} \underset{\smile}{\delta A}^2 &= \frac{i\beta}{||\vec{k}||} \overline{m}^a \left(\frac{1}{\kappa}\underset{\smile}{\delta C}{}_a - \underset{\smile}{p}{}_a \right) \end{align} \tag{ A.87 }$$

$$\begin{align} \underset{\smile}{\delta A}^{\overline{2}} &= \frac{i\beta}{||\vec{k}||} m^a \left(\frac{1}{\kappa}\underset{\smile}{\delta C}{}_a - \underset{\smile}{p}{}_a \right) \end{align} \tag{ A.88 }$$

$$\begin{align} \underset{\smile}{\delta A}^3 &= -\frac{i\kappa||\vec{k}||}{2} \hat{k}_i \left[\underset{\smile}{\delta\mathcal{G}}^i+ \underset{\smile}{\delta \xi}^i\right] \end{align} \tag{ A.89 }$$

$$\begin{align} \underset{\smile}{\delta A}^4 &=- \frac{i\beta}{2||\vec{k}||} \hat{k}^a \left( \frac{1}{\kappa} \underset{\smile}{\delta C}{}_a - \underset{\smile}{p}{}_a \right) + \frac{1}{2||\vec{k}||} \left(\frac{1}{\kappa} \underset{\smile}{\delta C} - \underset{\smile}{\epsilon} \right) \end{align} \tag{ A.90 }$$

$$\begin{align} \underset{\smile}{\delta A}^5 &=- \frac{i\beta}{2||\vec{k}||} \hat{k}^a \left( \frac{1}{\kappa} \underset{\smile}{\delta C}{}_a - \underset{\smile}{p}{}_a \right) - \frac{1}{2||\vec{k}||} \left(\frac{1}{\kappa} \underset{\smile}{\delta C} - \underset{\smile}{\epsilon} \right) \end{align} \tag{ A.91 }$$

$$\begin{align} \underset{\smile}{\delta E}^1 &= -\frac{2i\beta}{\kappa||\vec{k}||} \overline{m}^i \underset{\smile}{\delta G}{}_i + \frac{i\beta}{||\vec{k}||^2} \overline{m}^a\left(\frac{1}{\kappa} \underset{\smile}{\delta C}{}_a- \underset{\smile}{p}{}_a\right) +\frac{i\kappa}{2} \overline{m}^i \left[\underset{\smile}{\delta\mathcal{G}}{}_i + \underset{\smile}{\delta \xi}{}_i\right] \end{align} \tag{ A.92 }$$

$$\begin{align} \underset{\smile}{\delta E}^{\overline{1}} &= -\frac{2i\beta}{\kappa||\vec{k}||} m^i \underset{\smile}{\delta G}{}_i - \frac{i\beta}{||\vec{k}||^2} m^a \left(\frac{1}{\kappa} \underset{\smile}{\delta C}{}_a- \underset{\smile}{p}{}_a\right) -\frac{i\kappa}{2} m^i \left[\underset{\smile}{\delta\mathcal{G}}{}_i+ \underset{\smile}{\delta \xi}{}_i\right] \end{align} \tag{ A.93 }$$

$$\begin{align} \underset{\smile}{\delta E}^2 &= -i\kappa||\vec{k}|| \overline{m}_a \left[\underset{\smile}{\delta \mathcal{G}}^a+ \underset{\smile}{\delta \sigma}^a\right] -\frac{i\kappa}{2} \overline{m}_i \left[\underset{\smile}{\delta\mathcal{G}}^i+ \underset{\smile}{\delta \xi}^i\right] -\frac{\beta i}{||\vec k||^2} \overline{m}_a \left(\frac{1}{\kappa} \underset{\smile}{\delta C}^a - \underset{\smile}{p}{}^a\right) \end{align} \tag{ A.94 }$$

$$\begin{align} \underset{\smile}{\delta E}^{\overline{2}} &= -i\kappa||\vec{k}|| m_a\left[ \underset{\smile}{\delta \mathcal{G}}^a+ \underset{\smile}{\delta \sigma}^a\right] +\frac{i\kappa}{2} m_i \left[\underset{\smile}{\delta\mathcal{G}}^i+ \underset{\smile}{\delta \xi}^i\right] +\frac{\beta i}{||\vec k||^2} m_a \left(\frac{1}{\kappa} \underset{\smile}{\delta C}{}_a- \underset{\smile}{p}{}^a\right) \end{align} \tag{ A.95 }$$

$$\begin{align} \underset{\smile}{\delta E}^3 &=-\frac{2i\beta}{\kappa ||\vec{k}||} \hat{k}^i \underset{\smile}{\delta G}{}_i- \frac{1}{||\vec{k}||^2} \left( \frac{1}{\kappa} \underset{\smile}{\delta C} - \underset{\smile}{\epsilon}\right) \end{align} \tag{ A.96 }$$

$$\begin{align} \underset{\smile}{\delta E}^4 & = -\kappa\beta ||\vec{k}||\left[ \underset{\smile}{\delta \mathcal{G}}+ \underset{\smile}{\delta \tau}\right] + i\kappa ||\vec{k}|| \hat{k}^a \left[\underset{\smile}{\delta \mathcal{G}}{}_a+ \underset{\smile}{\delta \sigma}{}_a\right] -\frac{1}{||\vec{k}||^2} \left(\frac{1}{\kappa}\underset{\smile}{\delta C}-\underset{\smile}{\epsilon}\right) +\frac{i\kappa}{2} \hat{k}^i\left[ \underset{\smile}{\delta \mathcal{G}}{}_i+ \underset{\smile}{\delta \xi}^i\right] \nonumber\\ & \quad + \frac{i \beta}{||\vec{k}||^2} \hat{k}^a \left( \frac{1}{\kappa} \underset{\smile}{\delta C}{}_a - \underset{\smile}{p}{}_a\right) \end{align} \tag{ A.97 }$$

$$\begin{align} \underset{\smile}{\delta E}^5 & = \kappa\beta ||\vec{k}|| \left[\underset{\smile}{\delta \mathcal{G}}+ \underset{\smile}{\delta \tau}\right] + i\kappa ||\vec{k}|| \hat{k}^a\left[ \underset{\smile}{\delta \mathcal{G}}{}_a+ \underset{\smile}{\delta \sigma}{}_a\right] -\frac{1}{||\vec{k}||^2}\left(\frac{1}{\kappa} \underset{\smile}{\delta C}-\underset{\smile}{\epsilon}\right) -\frac{i\kappa}{2} \hat{k}^i\left[ \underset{\smile}{\delta \mathcal{G}}{}_i+ \underset{\smile}{\delta \xi}^i\right] \nonumber\\ & \quad - \frac{i \beta}{||\vec{k}||^2} \hat{k}^a \left( \frac{1}{\kappa} \underset{\smile}{\delta C}{}_a - \underset{\smile}{p}{}_a\right) \,. \end{align} \tag{ A.98 }$$

All the appearing matter degrees of freedom must be substituted by their gauge invariant extensions plus the corresponding correction terms, what can be read off from (2.63) and (2.64), from which we obtain up to first order in κ:

$$\begin{align} \phi(\vec{x},t) &= \phi^{GI}(\delta\sigma^c,\delta\tau,\delta\xi^j) + \kappa^2 (\delta T(\vec{x},t) - \delta\tau(\vec{x},t)) \pi(\vec{x},t) + \kappa^2 (\delta T^c(\vec{x},t) - \delta\sigma^c(\vec{x},t)) \partial_c \phi(\vec{x},t)\nonumber\\ \end{align} \tag{ A.99 }$$

$$\begin{align} \pi(\vec{x},t) & =\pi^{GI}(\delta\sigma^c,\delta\tau,\delta\xi^j) + \kappa^2 \partial_a [(\delta T(\vec{x},t) - \delta\tau(\vec{x},t)) \partial^a\phi(\vec{x},t)] \nonumber\\ & \quad + \kappa^2 \partial_c[(\delta T^c(\vec{x},t) - \delta\sigma^c(\vec{x},t)) \pi(\vec{x},t)] - \kappa^2 (\delta T(\vec{x},t) - \delta\tau(\vec{x},t))\, m^2 \phi(\vec{x},t)\,. \end{align} \tag{ A.100 }$$

Note that when considering terms where any of the gravitational variables $\underset{\smile}{\delta A}$ or $\underset{\smile}{\delta E}$ appear with a prefactor κ, then we can drop all the correction terms of order κ in (A.99) and (A.100) as they would then lead to terms of order κ².

Using the transformation derived in the previous subsection, it is possible to express the canonical Hamiltonian $\delta{\mathbf{H}}_\textrm{can}$ in (2.40) in terms of the new phase space variables consisting of the physical gauge invariant degrees of freedom, the constraints and the clocks. In a first step we rewrite the Hamiltonian in the following way:

$$\begin{align} \delta{\mathbf{H}}_\textrm{can} =\int_\sigma d^3x \bigg[ &\epsilon(\phi,\pi) + \delta N^a\, \delta C_a + \delta N\, \delta C + \delta \Lambda^i\, \delta G_i \nonumber\\[6pt] & +\kappa\, ({{\partial}{{}_a}}\phi) \,({{\partial}{{}^b}}\phi)\, {{\delta E}{^a{}_i}}\, \delta^i_b +\frac{\kappa}{2} m^2 \phi^2\, {{\delta E}{^a{}_i}}\,\delta^i_a - \frac{\kappa}{2} {{\delta E}{^a{}_i}}\,\delta_a^i\, \epsilon(\phi,\pi)+ \nonumber\\[6pt] &- \frac{\kappa}{2\beta^2} \epsilon^{jkl} {{\epsilon}{_{jmn}}}\, \delta^a_k\, \delta^b_l \left( {{\delta A}{_a{}^m}} {{\delta A}{_b{}^n}}+ (\beta^2+1) {{\delta \Gamma}{_a{}^m}} {{\delta \Gamma}{_b{}^n}} -2 {{\delta \Gamma}{_a{}^m}} {{\delta A}{_b{}^n}} \right) \bigg],\nonumber\\ \end{align} \tag{ A.101 }$$

where in the first line there is the background matter Hamiltonian and the constraints, in the second line the interaction part and in the last line the second order of the Hamilton constraint. The first line contains the background Hamiltonian of the scalar field, which does not need to be transformed as there are no background constraints and hence it is already gauge invariant, see discussion above. In order to obtain the physical Hamiltonian as described in (2.75) in terms of the new variables, one first has to determine the observables (A.45) and (A.46) in terms of the new variables, where we already set $\delta \xi^j = 0$ as explained in the main text:

$$\begin{align} (\underset{\smile}{\delta E}{}{^a{}_i})^{GI} & = \underset{\smile}{\delta \mathcal E}{}{^a{}_i}-\frac{2i\beta}{\kappa||\vec k||} \hat{k}^a \underset{\smile}{\delta G}{}_i + \frac{\beta}{\kappa||\vec k||^2} {\epsilon}{^{ba}{}_i} \left[\underset{\smile}{\delta C}{}_b-\kappa\underset{\smile}{p}{}_b\right]-\frac{1}{\kappa||\vec k||^2} \delta^a_i \left[ \underset{\smile}{\delta C}-\kappa\underset{\smile}{\epsilon}\right] \nonumber\\[6pt] & \quad -i\kappa||\vec k|| \underset{\smile}{\delta \sigma}^b \left( \delta^a_b \hat{k}_i - \delta^a_i \hat{k}_b \right) -\kappa ||\vec k||\beta \underset{\smile}{\delta \tau} i\hat{k}_c {\epsilon}{_i{}^{ac}{}} +O(\delta^3,\kappa^2) \end{align} \tag{ A.102 }$$

$$\begin{align} (\underset{\smile}{\delta A}{}{_a{}^i})^{GI} & = \underset{\smile}{\delta \mathcal A}{}{_a{}^i}+\frac{i\beta}{2\kappa||\vec k||} (2\hat{k}_i\delta^{ab} -\hat{k}^b \delta^a_i -\hat{k}^i\hat{k}_a\hat{k}^b) \left[ \underset{\smile}{\delta C}{}_b - \kappa\underset{\smile}{p}{}_b\right] -\frac{1}{2\kappa||\vec k||} i{\epsilon}{_{ba}{}^i} \hat{k}^b \left[ \underset{\smile}{\delta C}-\kappa\underset{\smile}{\epsilon}\right]\nonumber\\ & \quad +O(\delta^3,\kappa^2) \,. \end{align} \tag{ A.103 }$$

The interaction part in (2.75), $O_{{\mathcal H}_I,\{T\}}$, then assumes the following form:

$$\begin{align} O_{{\mathcal H}_I,\{T\}} & = \kappa\int d^3k \; \underset{\smile}{\varphi}{}_a^b \; \underset{\smile}{\delta \mathcal{E}}^{^{\prime}}{}{^a{}_i}\delta_b^i +\kappa \int d^3k ik_b \underset{\smile}{\delta \sigma}^{^{\prime} a} \left[ \underset{\smile}{\varphi}{}_a^b - \delta_a^b \left( \underset{\smile}{\varphi}{}_c^c +m^2 \underset{\smile}{\widetilde{V}}-\underset{\smile}{\epsilon} \right)\right] \nonumber\\[9pt] & \quad +\int \frac{d^3k}{||\vec k||^2} \;\left(- \left[ \underset{\smile}{\varphi}{}_a^a + \frac{3}{2} m^2\underset{\smile}{\widetilde{V}} -\frac{3}{2}\underset{\smile}{\epsilon}\right] \left(\underset{\smile}{\delta C^{^{\prime}}}-\kappa\underset{\smile}{\epsilon}^{^{\prime}}\right) +2i\beta k^a \underset{\smile}{\delta G}^{^{\prime}}{}_i\right. \nonumber\\[9pt] & \qquad\qquad\quad\qquad \times \left.\left[ \underset{\smile}{\varphi}{}_a^b \delta^i_b +\left(\frac{1}{2}m^2 \underset{\smile}{\widetilde{V}}-\frac{1}{2}\underset{\smile}{\epsilon}\right) \delta_a^i\right] \right) + O(\delta^3,\kappa^2)\,, \end{align} \tag{ A.104 }$$

where $\epsilon = \epsilon(\phi^{GI},\pi^{GI})$ and $p_b = p_b(\phi^{GI},\pi^{GI})$. We used here a notation where no argument means dependency on $(\vec k,t)$ and a prime after a quantity means dependency on $(-\vec k,t)$. The last line in (A.101) yields $O_{{\mathcal H}_\textrm{geo},\{T\}}$:

$$\begin{align} O_{{\mathcal H}_\textrm{geo},\{T\}}& = \frac{\kappa}{2\beta^2}\int d^3k \sum_{r\in\{\pm\}}\left[\underset{\smile}{\delta A}^r \underset{\smile}{\delta A}^{^{\prime} r} +||\vec{k}||^2 (\beta^2+1) \underset{\smile}{\delta E}^r \underset{\smile}{\delta E}^{^{\prime} r} + 2r||\vec{k}|| \underset{\smile}{\delta E}^r \underset{\smile}{\delta A}^{^{\prime} r} \right] \nonumber\\[9pt] & \quad -\kappa\int \frac{d^3k}{||\vec{k}||^2} \bigg\{-\left( \frac{1}{\kappa} \underset{\smile}{\delta C}{}_a-\underset{\smile}{p}{}_a\right)\left( \frac{1}{\kappa} \underset{\smile}{\delta C}^{^{\prime}}{}_b-\underset{\smile}{p}^{^{\prime}}{}_b\right) \left(\delta^{ab}-\frac{1}{4}\hat{k}{}^a\hat{k}{}^b\right)\nonumber\\[9pt] & \quad\qquad\qquad\quad\quad +\frac{1}{4} \left( \frac{1}{\kappa} \underset{\smile}{\delta C} -\underset{\smile}{\epsilon}\right) \left( \frac{1}{\kappa} \underset{\smile}{\delta C}^{^{\prime}} -\underset{\smile}{\epsilon}^{^{\prime}}\right)\Bigg\}\nonumber\\[9pt] & \quad - \kappa\int d^3k \; \underset{\smile}{\delta G}{}_i \left[ \frac{\beta^2+1}{\kappa^2} \hat{k}^i\hat{k}^j\ \underset{\smile}{\delta G}^{^{\prime}}{}_j + \frac{i\beta}{\kappa||\vec{k}||} \hat{k}^i \left(\frac{1}{\kappa} \underset{\smile}{\delta C}^{^{\prime}}-\underset{\smile}{\epsilon}^{^{\prime}} \right)\right. \nonumber \\[9pt] & \qquad\qquad\qquad\qquad \left. - \frac{2i}{\kappa||\vec{k}||} \epsilon_{abi} \hat{k}^a \left( \frac{1}{\kappa} \underset{\smile}{\delta C}^{^{\prime}}{}_b-\underset{\smile}{p}^{^{\prime}}{}_b\right)\right]\nonumber \\[9pt] & \quad -\kappa \int d^3k \left[ i \kappa ||\vec k|| \underset{\smile}{\delta \tau} \hat{k}^b \left(\frac{1}{\kappa} \underset{\smile}{\delta C}^{^{\prime}}{}_b-\underset{\smile}{p}^{^{\prime}}{}_b \right) \right] +O(\delta^3,\kappa^2)\,. \end{align} \tag{ A.105 }$$

As discussed in the main text in section 2.4 in (2.75) the physical Hamiltonian $\delta \mathbf{H}$ consists in zeroth order of the Hamiltonian of the scalar field on Minkowski spacetime plus the Dirac observable associated with $\delta^2C^\textrm{geo}$ as well as the Dirac observable corresponding to the interaction Hamiltonian ${\mathcal H}_I$. In the following equations these contributions have been marked in blue:

$$\begin{align} \delta \mathbf{H}_\textrm{can} & = {\int_{\mathbb{R}^3} d^3x \;\epsilon(\vec{x},t)} {+\kappa \int_{\mathbb{R}^3}d^3k\; \frac{1}{2\beta^2}\sum_{r\in\{\pm\}}\Big( \underset{\smile}{\delta A}^r(\vec{k},t) \,\underset{\smile}{\delta A}^r (-\vec{k},t)+ 2r\,||\vec{k}|| \underset{\smile}{\delta E}^r(\vec{k},t) \, \underset{\smile}{\delta A}^r(-\vec{k},t) }\nonumber\\ & \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad {+ (\beta^2+1) ||\vec{k}||^2 \underset{\smile}{\delta E}^r(\vec{k},t)\, \underset{\smile}{\delta E}^r(-\vec{k},t)\Big) }\nonumber\\[-3pt] & \quad {+\kappa \int_{\mathbb{R}^3} d^3k \;\delta_{b}^i \underset{\smile}{\varphi} {}_a^b(\vec{k},t) \, \underset{\smile}{\delta \mathcal{E}}{}{^a{}_i}(-\vec{k},t)}\nonumber\\ & \quad {+\kappa \int d^3k ik_b \kappa\underset{\smile}{\delta \sigma}^{\prime{a}} \left[ \underset{\smile}{\varphi}{}_a^b - \delta_a^b \left( \underset{\smile}{\varphi}{}_c^c +m^2 \underset{\smile}{\widetilde{V}}-\underset{\smile}{\epsilon} \right)\right]+\kappa \int d^3k \left[i \kappa ||\vec k|| \underset{\smile}{\delta \tau} \hat{k}^b \underset{\smile}{p}^{^{\prime}}{}_b \right]}\nonumber\\ & \quad {-\kappa\int_{\mathbb{R}^3} \frac{d^3k}{||\vec{k}||^2} \Bigg( -\underset{\smile}{p}{{}{_c}}(\vec{k},t) \,\underset{\smile}{p}{}{_d}(-\vec{k},t) \left[ \delta^{cd}- \frac{1}{4} \hat{k}^c \hat{k}^d \right] + \underset{\smile}{\epsilon}(\vec{k},t)\,\underset{\smile}{\epsilon}(-\vec{k},t) \left[\frac{1}{4}+\frac{3}{2}\right]} \nonumber\\ & \qquad\qquad\qquad\quad{-\frac{3}{2} m^2\,\underset{\smile}{\widetilde{V}}(\vec{k},t) \,\underset{\smile}{\epsilon}(-\vec{k},t)-\underset{\smile}{\varphi} {}_a^a(\vec{k},t) \,\underset{\smile}{\epsilon}(-\vec{k},t) \Bigg)}\nonumber\\ & \quad + \int {d^3k} \; \underset{\smile}{\delta G}{}_i \Bigg( \underset{\smile}{\delta \Lambda}{}^{^{\prime} i} - 2i\beta k^a \frac{1}{||\vec k||^2}\left[ \underset{\smile}{\varphi}^{^{\prime}}{}_a^b \delta^i_b +\frac{1}{2}\left(m^2 \widetilde{V}^{^{\prime}}-\epsilon^{^{\prime}}\right) \delta_a^i\right] -\frac{\beta^2+1}{\kappa^2} \hat{k}^i\hat{k}^j\ \underset{\smile}{\delta G}^{^{\prime}}{}_j \nonumber\\ & \qquad\qquad\qquad\quad - \frac{i\beta}{\kappa||\vec{k}||} \hat{k}^i \left(\frac{1}{\kappa} \underset{\smile}{\delta C}^{^{\prime}}-\underset{\smile}{\epsilon}^{^{\prime}} \right) + \frac{2i}{||\vec{k}||} \epsilon_{bia} \hat{k}^a \left( \frac{1}{\kappa} \underset{\smile}{\delta C}^{^{\prime}}{}_b-\underset{\smile}{p}^{^{\prime}}{}_b\right) \Bigg)\nonumber\\ & \quad + \int d^3k \; \underset{\smile}{\delta C}{}_a \Bigg(\underset{\smile}{\delta N}^{^{\prime} a} +\frac{1}{||\vec k||^2} \left( \frac{1}{\kappa} \underset{\smile}{\delta C}^{^{\prime}}{}_b -2 \underset{\smile}{p}^{^{\prime}}{}_b \right) \left(\delta^{ab}-\frac{1}{4}\hat{k}{}^a\hat{k}{}^b\right)\Bigg) \nonumber\\ & \quad + \int d^3k \; \underset{\smile}{\delta C} \Bigg( \underset{\smile}{\delta N}^{^{\prime}}- \frac{1}{||\vec k||^2}\left[ \underset{\smile}{\varphi}^{^{\prime}}{}_a^a + \frac{3}{2} m^2\underset{\smile}{\widetilde{V}}^{^{\prime}} -\frac{3}{2}\underset{\smile}{\epsilon}^{^{\prime}}\right] -\frac{1}{4||\vec k||^2} \left( \frac{1}{\kappa} \underset{\smile}{\delta C}^{^{\prime}} -2\underset{\smile}{\epsilon}^{^{\prime}} \right) \Bigg)\nonumber\\ & \quad -\kappa \int d^3k \left[ i ||\vec k|| \underset{\smile}{\delta \tau} \hat{k}^b \underset{\smile}{\delta C}^{^{\prime}}{}_b\right] + O(\delta^3,\kappa^2)\, \end{align} \tag{ A.106 }$$

with

$$\begin{align} \delta \mathbf{H} & = {\int_{\mathbb{R}^3} d^3x \;\epsilon(\vec{x},t)} {+\kappa \int_{\mathbb{R}^3}d^3k\; \frac{1}{2\beta^2}\, \sum_{r\in\{\pm\}}\Big( \underset{\smile}{\delta A}^r(\vec{k},t) \,\underset{\smile}{\delta A}^r (-\vec{k},t)+ 2r\,||\vec{k}|| \underset{\smile}{\delta E}^r(\vec{k},t) \, \underset{\smile}{\delta A}^r(-\vec{k},t) }\nonumber\\ & \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \ {+ (\beta^2+1) ||\vec{k}||^2 \underset{\smile}{\delta E}^r(\vec{k},t)\, \underset{\smile}{\delta E}^r(-\vec{k},t)\Big) }\nonumber\\[-3pt] & \quad {+\kappa \int_{\mathbb{R}^3} d^3k \;\delta_{b}^i \underset{\smile}{\varphi} {}_a^b(\vec{k},t) \, \underset{\smile}{\delta \mathcal{E}}{}{^a{}_i}(-\vec{k},t)}\nonumber\\ & \quad {+\kappa \int d^3k\; ik_b \kappa\underset{\smile}{\delta \sigma}^{^{\prime} a} \left[ \underset{\smile}{\varphi}{}_a^b - \delta_a^b \left( \underset{\smile}{\varphi}{}_c^c +m^2 \underset{\smile}{\widetilde{V}}-\underset{\smile}{\epsilon} \right)\right] +\kappa \int d^3k \left[i \kappa ||\vec k|| \underset{\smile}{\delta \tau} \hat{k}^b \underset{\smile}{p}^{^{\prime}}{}_b \right]}\nonumber\\ & \quad {-\kappa\int_{\mathbb{R}^3} \frac{d^3k}{||\vec{k}||^2} \Bigg( -\underset{\smile}{p}{{}{_c}}(\vec{k},t) \,\underset{\smile}{p}{}{_d}(-\vec{k},t) \left[ \delta^{cd}-\frac{1}{4} \hat{k}^c \hat{k}^d \right] + \underset{\smile}{\epsilon}(\vec{k},t)\,\underset{\smile}{\epsilon}(-\vec{k},t) \left[\frac{1}{4}+\frac{3}{2}\right]} \nonumber\\ & \qquad\qquad\qquad\quad \ {-\frac{3}{2} m^2\,\underset{\smile}{\widetilde{V}}(\vec{k},t) \,\underset{\smile}{\epsilon}(-\vec{k},t)-\underset{\smile}{\varphi} {}_a^a(\vec{k},t) \,\underset{\smile}{\epsilon}(-\vec{k},t) \Bigg)} + O(\delta^3,\kappa^2). \end{align} \tag{ A.107 }$$

As expected the physical Hamiltonian differs from $\mathbf{\delta H_{can}}$ only by terms that involve the linearised constraints at least linearly and $\delta \mathbf{H}$ agrees exactly with the expression for the physical Hamiltonian shown in (2.75) in the main text. As discussed in the main text in section 2.4 below (2.67) as well as (2.76) the fall-off behaviour of the physical Hamiltonian is suitable for the choices of δτ and $\delta\sigma^c$ in (2.67) where we also discuss the fall-off behaviour of the gauge invariant variables for these choices.

In this appendix we present a possible gauge fixing that naturally arises from the choice of the clocks introduced in section 2.4 in the context of constructing Dirac observables. The induced gauge fixing is obtained by setting the gauge fixing conditions ${\mathcal G}^I$ shown in (2.74) in section 2.4 equal to zero. Evaluating the constraints in momentum space, using the basis introduced in appendix A.2, yields the following seven relations among the elementary phase space variables at every point $(\vec k,t)$:

$$\begin{align} ||\vec{k}|| \underset{\smile}{\delta E}^{\overline{1}} &= \underset{\smile}{\delta A}^{\overline{1}} - \underset{\smile}{\delta A}^{\overline{2}} \qquad\qquad\qquad i ||\vec{k}|| (\underset{\smile}{\delta A}^4 + \underset{\smile}{\delta A}^5) = - \beta \hat{k}^a \underset{\smile}{p}{}{_a{}} \end{align} \tag{ B.1 }$$

$$\begin{align} ||\vec{k}|| \underset{\smile}{\delta E}^{1} &= \underset{\smile}{\delta A}^{2} - \underset{\smile}{\delta A}^{1} \qquad\qquad\qquad i ||\vec{k}|| \underset{\smile}{\delta A}^{\overline{2}} = \beta m^a(\vec{k}) \underset{\smile}{p}{}{_a{}} \end{align} \tag{ B.2 }$$

$$\begin{align} ||\vec{k}|| \underset{\smile}{\delta E}^{3} &= \underset{\smile}{\delta A}^{5} - \underset{\smile}{\delta A}^{4} \qquad\qquad\qquad i ||\vec{k}|| \underset{\smile}{\delta A}^{2} = \beta \overline{m}^a(\vec{k}) \underset{\smile}{p}{}{_a{}} \end{align} \tag{ B.3 }$$

$$\begin{align} ||\vec{k}|| (\underset{\smile}{\delta A}^4 - \underset{\smile}{\delta A}^5) &= -\underset{\smile}{\epsilon}\,, \end{align} \tag{ B.4 }$$

where $\underset{\smile}{\epsilon}$ and $\underset{\smile}{p}{}_a$ denote the three-dimensional Fourier transforms of the energy and momentum density of the scalar field. The gauge fixing conditions, obtained by requiring ${\mathcal G}^I$ to vanish and in addition choosing $\tau(t,\vec{x}) = t, \sigma^a = x^a_\sigma = \textrm{const}$ and $\xi^j = 0$, hence $\delta \tau = \frac{t}{\kappa}$, $\delta\sigma^a = \frac{x_\sigma^a}{\kappa} = \textrm{const}$, give seven additional conditions in momentum space for every $\vec{k}$ and t:

$$\begin{align} \underset{\smile}{\delta A}^3 &=0 \qquad\qquad\qquad\qquad\qquad\qquad\quad\kern57.4pt \underset{\smile} {\delta E}^2 = -\underset{\smile}{\delta E}^1 \end{align} \tag{ B.5 }$$

$$\begin{align} \underset{\smile}{\delta A}^1 &= 0 \qquad\qquad\qquad\qquad\qquad\qquad\quad\kern57.4pt \underset{\smile} {\delta E}^{\overline{2}} = -\underset{\smile}{\delta E}^{\overline{1}} \end{align} \tag{ B.6 }$$

$$\begin{align} \underset{\smile}{\delta A}^{\overline{1}} &= 0 \qquad\qquad\qquad\qquad\qquad\qquad\quad 3 \underset{\smile}{\delta E}^3 - \underset{\smile}{\delta E}^4 - \underset{\smile}{\delta E}^5 = \frac{1}{||\vec k||^2} \underset{\smile}{\epsilon} \end{align} \tag{ B.7 }$$

$$\begin{align} \underset{\smile}{\delta E}^5- \underset{\smile}{\delta E}^4 &= \frac{2}{||\vec{k}||} \left( \underset{\smile}{\delta A}^3 + \underset{\smile}{\delta A}^4 + \underset{\smile}{\delta A}^5 \right)\,. \end{align} \tag{ B.8 }$$

Note that δτ or δσ were chosen to be constant in position and hence vanish in all appearing combinations in the gauge fixing conditions. Substituting these results into (A.40) and (A.41) yields an expression for $\underset{\smile}{\delta E}$ and $\underset{\smile}{\delta A}$ in terms of the physical degrees of freedom:

$$\begin{align} \underset{\smile}{\delta E}{}{^a{}_i} & = \underset{\smile}{\delta E}^{+} m^a m_i + \underset{\smile}{\delta E}^{-} \overline{m}^a \overline{m}_i + \frac{\underset{\smile}{\epsilon}}{||\vec{k}||^2}\delta^a_i+ \frac{\beta i}{||\vec{k}||^2} \underset{\smile}{p}{{}{_c}} \nonumber\\ & \quad \times \left[ m^c \left(\hat{k}^a \overline{m}_i - \overline{m}^a \hat{k}_i\right) - \overline{m}^c \left( \hat{k}^a m_i- m^a \hat{k}_i \right) - \hat{k}^c \left(m^a \overline{m}_i - \overline{m}^a m_i \right) \right] \nonumber\\ & = \underset{\smile}{\delta E}^{+} m^a m_i + \underset{\smile}{\delta E}^{-} \overline{m}^a \overline{m}_i + \frac{\underset{\smile}{\epsilon}}{||\vec{k}||^2}\delta^a_i -\frac{\beta}{||\vec k||^2} \underset{\smile}{p}{}_c {\epsilon}{^{ca}{}_i} \end{align} \tag{ B.9 }$$

$$\begin{align} \underset{\smile}{\delta A}{}{_a{}^i} & = \underset{\smile}{\delta A}^{+} m_a m^i + \underset{\smile}{\delta A}^{-} \overline{m}_a \overline{m}^i - \frac{\beta i}{||\vec{k}||} \underset{\smile}{p}{{}{_c}} \left[ \overline{m}^c m_a \hat{k}^i + m^c \overline{m}_a \hat{k}^i - \frac{1}{2} \hat{k}^c \left( m_a \overline{m}^i + \overline{m}_a m^i \right) \right] \nonumber\\ & \quad + \frac{\underset{\smile}{\epsilon}}{2 ||\vec{k}||} \left(\overline{m}_a m^i - m_a \overline{m}^i \right)\nonumber\\ & = \underset{\smile}{\delta A}^{+} m_a m^i + \underset{\smile}{\delta A}^{-} \overline{m}_a \overline{m}^i - \frac{\beta i}{||\vec{k}||} \underset{\smile}{p}{{}{_c}} \left[ \hat{k}^i \delta_a^c -\frac{1}{2} \hat{k}^c \delta_a^i -\frac{1}{2}\hat{k}^c\hat{k}_a\hat{k}^i\right] + \frac{i\underset{\smile}{\epsilon}}{2||\vec k||} {\epsilon}{_{ba}{}^i}\hat{k}^b\,. \end{align} \tag{ B.10 }$$

Note that one can also obtain setting the constraints equal to zero in the expressions shown in (A.102) and (A.103). The independent physical degrees of freedom are encoded in the matter fields $(\phi,\pi)$ as well as the symmetric transverse-traceless parts of the connection and densitised triads that we denote as

$$\begin{align} \underset{\smile}{\delta \mathcal{A}}{}{_a{}^i}(\vec{k},t) &:= \underset{\smile}{\delta A}^{+}(\vec{k},t) \, m_a(\vec{k}) \, m^i(\vec{k}) + \underset{\smile}{\delta A}^{-}(\vec{k},t) \, \overline{m}_a(\vec{k}) \, \overline{m}^i(\vec{k}) = \underset{\smile}{P}{}{_a{}^i{}_j{}^b}(\vec{k}) \underset{\smile}{\delta A}{}{_b{}^j}(\vec{k},t) \end{align} \tag{ B.11 }$$

$$\begin{align} \underset{\smile}{\delta \mathcal{E}}{}{^a{}_i}(\vec{k},t) &:= \underset{\smile}{\delta E}^{+}(\vec{k},t) \, m^a(\vec{k}) \, m_i(\vec{k}) + \underset{\smile}{\delta E}^{-}(\vec{k},t) \, \overline{m}^a(\vec{k}) \, \overline{m}_i(\vec{k}) = \underset{\smile}{P}{}{^a{}_i{}^j{}_b}(\vec{k}) \underset{\smile}{\delta E}{}{^b{}_j}(\vec{k},t) \end{align} \tag{ B.12 }$$

with the projectors $\underset{\smile}{P}{}{^a{}_i{}^j{}_b}(\vec{k})$ and $\underset{\smile}{P}{}{_a{}^i{}_j{}^b}(\vec{k})$ onto the symmetric transverse-traceless parts that were defined in (A.42). As expected, with the components of $\underset{\smile}{\delta \mathcal{A}}{}{_a{}^i},\underset{\smile}{\delta \mathcal{E}}{}{^a{}_i}(\vec{k},t)$ we end up with the six physical components given by

$$\begin{align} \underset{\smile}{\delta E}^{+}, \underset{\smile}{\delta E}^{-}, \underset{\smile}{\delta A}^{+}, \underset{\smile}{\delta A}^{-}, \phi, \pi\,. \end{align} \tag{ B.13 }$$

Without the matter fields, one gets the well-known four phase space field degrees of freedom of gravitational waves in vacuum, corresponding to two field degrees of freedom in the Lagrangian framework, and the real scalar matter field leads to two additional phase space field degrees of freedom. Reinserting (B.9) and (B.10) into the Hamiltonian (2.40) yields the total Hamiltonian of the linearised theory on the reduced phase space, (2.76). For this, we defined the transverse-traceless projectors

$$\begin{align} [P^\pm(\vec{k})]_a^b &:= \overline{m}_a(\pm\vec{k})\, \overline{m}^b(\pm\vec{k})\,, \end{align} \tag{ B.14 }$$

such that $[P^\pm(\vec{k})]_a^b \delta_b^i \, \underset{\smile}{\delta E}{}{^a{}_i} = \underset{\smile}{\delta E}^\pm$, $[P^\pm(\vec{k})]_b^a \, \delta^b_i \underset{\smile}{\delta A}{}{_a{}^i} = \underset{\smile}{\delta A}^\pm$ and $\underset{\smile}{P}{}{^a{}_i{}^j{}_b}(\vec{k}) = \sum\limits_{r\in\{\pm\}} [P^r(\vec{k})]^a_i$ $[P^{-r}(\vec{k})]^j_b$.

At this point, we would like to make a short comparison to the clocks introduced in [69] based on seminal work in [74, 102]. Likewise to [69] we denotes this set of clocks as ADM clocks. These clocks, similarly to the ones we use in this work, also do not contain the symmetric transverse-traceless gravitational degrees of freedom $\delta A^\pm$ and $\delta E^\pm$. Hence, these four phase space degrees of freedom are again identified as the physical degrees of freedom in the gravitational sector, as in our case. However, the gauge fixing induced by setting the corresponding gauge fixing conditions $\delta{\mathcal G}^I = \delta T^I-\tau^I$ used in [69] to zero is different from the one used by in section 2.4 in this work here, hence this set of clocks is not equivalent to the one discussed in section 2.4. This can be easily seen by considering the clock $^GT^i$ associated with the Gauß constraint from [69] in momentum space. Requiring that the components $m_i(\vec k) ^G\underset{\smile}{T}^i(\vec k)$ and $\overline{m}_i(\vec k) ^G\underset{\smile}{T}^i(\vec k)$ vanish yields $\underset{\smile}{\delta E}^1 = \underset{\smile}{\delta E}^{\overline{1}} = 0$.

In contrast the gauge fixing obtained by setting the constraints as well as gauge fixing conditions $\delta{\mathcal G}^I$ used in section 2.4 equal to zero, we can read off from the transformation introduced in appendix A.4 that $\underset{\smile}{\delta E}^1 = -\frac{i\beta}{||\vec{k}||^2} \overline{m}^a \underset{\smile}{p}{}_a$ and $\underset{\smile}{\delta E}^{\overline{1}} = \frac{i\beta}{||\vec{k}||^2} m^a \underset{\smile}{p}{}_a$. This gauge fixing condition would be equivalent to the Lorenz-like condition for the Gaußconstraint chosen in this work in section 2.4 in 2.49 only for the vacuum case, as can be seen in [61]. ⁹ Setting the gauge fixing conditions ${\mathcal G}^I_\textrm{ADM}$ induced by the ADM clocks used in [69] equal to zero imposes the following set of conditions for the phase space variables:

$$\begin{align} \underset{\smile}{\delta E}^1 = \underset{\smile}{\delta E}^{\overline{1}} = \underset{\smile}{\delta E}^2 =\underset{\smile}{\delta E}^{\overline{2}} = 2\underset{\smile}{\delta A}^3+\underset{\smile}{\delta A}^4 +\underset{\smile}{\delta A}^5 = \underset{\smile}{\delta E}^5-\underset{\smile}{\delta E}^4 = \underset{\smile}{\delta E}^3-\underset{\smile}{\delta E}^4-\underset{\smile}{\delta E}^5=0\,, \end{align} \tag{ B.15 }$$

which implies on the constraint hypersurface that

$$\begin{align} {\underset{\smile}{\delta E}}{}{^a{}_i} &= {\underset{\smile}{\delta \mathcal{E}}}{}{^a{}_i} + \frac{\underset{\smile}{\epsilon}}{2||\vec{k}||^2} (\delta^a_i + \hat{k}^a\hat{k}_i) \end{align} \tag{ B.16 }$$

$$\begin{align} {\underset{\smile}{\delta A}}{}{_a{}^i} &= {\underset{\smile}{\delta \mathcal{A}}}{}{_a{}^i}-\frac{i\beta}{||\vec k||} \underset{\smile}{p}{}_b \left( \hat{k}^i \delta^b_a + \hat{k}_a \delta^{ib} -\frac{1}{2} \hat{k}^b \delta_a^i -\hat{k}_a\hat{k}^i \hat{k}^b\right) +\frac{i\underset{\smile}{\epsilon}}{2||\vec k||} \hat{k}^b {\epsilon}{_{ba}{}^i}. \end{align} \tag{ B.17 }$$

In this particular gauge, the physical Hamiltonian becomes

$$\begin{align} \delta \mathbf{H} & = \int_{\mathbb{R}^3} d^3x \;\epsilon(\vec{x},t) +\kappa \int_{\mathbb{R}^3}d^3k\; \frac{1}{2\beta^2}\, \sum_{r\in\{\pm\}}\Big( \underset{\smile}{\delta A}^r(\vec{k},t) \,\underset{\smile}{\delta A}^r (-\vec{k},t)+ 2r\,||\vec{k}|| \underset{\smile}{\delta E}^r(\vec{k},t) \, \underset{\smile}{\delta A}^r(-\vec{k},t) \nonumber\\[6pt] & \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \ \ + (\beta^2+1) ||\vec{k}||^2 \underset{\smile}{\delta E}^r(\vec{k},t)\, \underset{\smile}{\delta E}^r(-\vec{k},t)\Big) \nonumber\\[6pt] & \quad +\kappa \int_{\mathbb{R}^3} d^3k \;\delta_{b}^i \underset{\smile}{\varphi} {}_a^b(\vec{k},t) \, \underset{\smile}{\delta \mathcal{E}}{}{^a{}_i}(-\vec{k},t)\nonumber\\[6pt] & \quad -\kappa\int_{\mathbb{R}^3} \frac{d^3k}{||\vec{k}||^2} \Bigg( -\underset{\smile}{p}{{}{_c}}(\vec{k},t) \,\underset{\smile}{p}{}{_d}(-\vec{k},t) \left[ \delta^{cd}- \frac{3}{4}\hat{k}^c \hat{k}^d \right] + \underset{\smile}{\epsilon}(\vec{k},t)\,\underset{\smile}{\epsilon}(-\vec{k},t) \left[\frac{1}{4}+1\right] \nonumber\\[6pt] & \qquad\qquad\qquad\quad \ \ - m^2\,\underset{\smile}{\widetilde{V}}(\vec{k},t) \,\underset{\smile}{\epsilon}(-\vec{k},t)-\underset{\smile}{\varphi} {}_a^b(\vec{k},t) \,\underset{\smile}{\epsilon}(-\vec{k},t) \left[ \frac{1}{2}\delta^a_b +\frac{1}{2}\hat{k}^a\hat{k}_b \right] \Bigg) + O(\delta^3,\kappa^2)\,.\nonumber\\ \end{align} \tag{ B.18 }$$

As can be readily seen, it differs from the Hamiltonian in the gauge used in this work in (A.107) only in the term encoding self-interaction of the scalar field. The reason is that by choosing the clocks, different choices for physical temporal and spatial coordinates were established and these are determined by the choice of different sets of clocks that on the gauge fixed surface are equal to δτ and $\delta\sigma^a$ respectively. As discussed in section 2.4 by an abuse of notation we understand the $\vec x$ and t arguments of the fields to be the associated values of the diffeomorphism and Hamiltonian clock, δτ and $\delta \sigma^a$ respectively. By inspection of the clocks, it turns out that the ADM clocks used in this section, corresponding to the parameters $x^a_{ADM}$ and t_ADM are related to the parameters used in this works, t and x^a , by

$$\begin{align} t_{ADM} &= \frac{1}{2}t \end{align} \tag{ B.19 }$$

$$\begin{align} x^a_{ADM} &= x^a+ \frac{1}{2}\partial^a (\epsilon * G^{\Delta\Delta}), \end{align} \tag{ B.20 }$$

where we neglected factors of κ and β, as the notation regarding these factors in [69] partly differs from the one used in this work.

We realise that in the case of vacuum gravity the choices for physical temporal and spatial coordinates agree but differ for non-vanishing momentum and energy density of the scalar field

Compared to the work in [35], where for the ADM constraints the same gauge fixing as in this section was used, their gauge fixed Hamiltonian appears with different prefactors in the self-interaction part compared to the one in (B.18). A reason might be, that we were not able to reproduce the expression for V in their equation (6) in [35] which consists of the second order terms of the expansion of $\sqrt{q_{ab}}\: ^{(3)}R$, where q_ab denotes the spatial ADM metric and $^{(3)}R$ the three-dimensional Ricci scalar. However, in contrast the result for this expansion given by the same authors in [38] could be reproduced by our computations.