We investigate the oneparticle sector for the fieldtheoretical model of gravitationally induced decoherence for a scalar field in 1 with a special focus on the renormalisation of the oneparticle master equation. In contrast to existing models in the literature where the renormalisation is usually performed after the Markov and rotating wave approximation and often only for certain limits such as the non or ultrarelativistic limit here we apply the renormalisation directly after the oneparticle projection. With this strategy we show that UVdivergent contributions in the oneparticle master equation can be identified with the vacuum contributions in the selfenergy of the scalar field in the effective quantum field theory and depending on the chosen oneparticle projection method its vacuum bubbles while the additional thermal contributions in the selfenergy are all UVfinite. To obtain the renormalised oneparticle master equation we use an onshell renormalisation procedure of the underlying effective QFT. We then apply the Markov and rotating wave approximation specifying a condition under which the Markov approximation can be applied in the case of the ultrarelativistic limit. We compare our results with those available in the literature. This includes an analysis of two different kinds of oneparticle projections a comparison of the application and effects of renormalisation of quantum mechanical and field theoretical models the nonrelativistic and ultrarelativistic limits of the renormalised oneparticle master equations and a comparison with a quantum mechanical toy model for gravitationally induced decoherence in the context of neutrino oscillations.
['High Energy Physics - Theory', 'General Relativity and Quantum Cosmology']
Gravitationally induced decoherence of a scalar field investigating the oneparticle sector and its interplay with renormalisation
{'Gravitationally induced decoherence of a scalar field: investigating the one-particle sector and its interplay with renormalisation': 'Max Joseph Fahn ∗ and Kristina Giesel †  \nInstitute for Quantum Gravity, Theoretical Physics III, Department of Physics, FAU Erlangen-Nurnberg, Staudtstr. 7, 91058 Erlangen, Germany.  \nWe investigate the one-particle sector for the field-theoretical model of gravitationally induced decoherence for a scalar field in [1] with a special focus on the renormalisation of the one-particle master equation. In contrast to existing models in the literature, where the renormalisation is usually performed after the Markov and rotating wave approximation and often only for certain limits such as the non- or ultra-relativistic limit, here we apply the renormalisation directly after the one-particle projection. With this strategy, we show that UV-divergent contributions in the one-particle master equation can be identified with the vacuum contributions in the self-energy of the scalar field in the effective quantum field theory and depending on the chosen one-particle projection method, its vacuum bubbles, while the additional thermal contributions in the self-energy are all UV-finite. To obtain the renormalised one-particle master equation, we use an on-shell renormalisation procedure of the underlying effective QFT. We then apply the Markov and rotating wave approximation, specifying a condition under which the Markov approximation can be applied in the case of the ultra-relativistic limit. We compare our results with those available in the literature. This includes an analysis of two different kinds of one-particle projections, a comparison of the application and effects of renormalisation of quantum mechanical and field theoretical models, the non-relativistic and ultra-relativistic limits of the renormalised one-particle master equations, and a comparison with a quantum mechanical toy model for gravitationally induced decoherence in the context of neutrino oscillations.', 'CONTENTS': 'References  \n74', 'I. INTRODUCTION': "Since we expect closed and thus isolated quantum systems to be an idealisation, there is increasing interest in various areas of physics in the investigation of open quantum systems. These are usually modelled by a total system consisting of two parts: the system under consideration and the environment as well as the interaction between the two. By tracing out the degrees of freedom of the environment, the effective dynamics of the system under consideration are obtained, which, in contrast to the closed system, is still affected by the interaction with the environment [2-4]. The effective dynamics of the system under consideration is often formulated in the form of a so-called master equation, of which the Lindblad equation is a special and prominent example [5, 6]. An important question in the context of open quantum systems is, what to choose as the environment and how to model the interaction with the system. Since we assume that all standard matter is at least coupled to gravity, interest in open quantum systems with a gravitational environment has increased in recent years, see for instance [7-13] and the reviews in [14, 15] and references therein. This has been investigated in the context of gravitationally induced decoherence due to a quantum gravitational environment in the framework of linearised gravity in several works in the literature. For instance in [1, 8, 9] a scalar field is considered as the matter field, while in [12, 16] a photon field is considered and in [11] a model for a generic bosonic matter field is investigated. All these works have in common that they start from the underlying field theoretical model and then derive a corresponding master equation for the matter system under consideration.  \nIn addition, there are also phenomenological models, see for instance [7, 10, 17-20], that do not necessarily derive a master equation from an underlying action, but rather set up physically motivated ansatze for the form of the dissipator, which encodes the effective influence of the environment in the master equation. The latter are often formulated in the context of quantum mechanics and thus for finitely many degrees of freedom, where a given master equation is often less difficult to work with or even to solve. Furthermore, many of these models already used the Lindblad equation as a starting point, so that the physical properties of the environment are often less accessible compared to the underlying field theoretical models. An interesting question is therefore, on the one hand, how the underlying field theoretical model can be linked to a given quantum-mechanical phenomenological model and, on the other hand, what additional insights the underlying field-theoretical model can provide that are no longer accessible to us in the quantum-mechanical model. In this work we want to address both question in the context of models for gravitationally induced decoherence.  \nTo answer these questions, the one-particle sector of field theory needs to be investigated to establish a link to microscopic quantum mechanical models. Furthermore, in order to investigate the connection to existing phenomenological models for gravitationally induced decoherence, one needs to understand more precisely how certain methods such as renormalisation and specific approximations such as the Markov or rotating wave approximation, which are often performed to finally arrive at a Lindblad-type master equation, affect the field-theoretical model respectively its one-particle sector. Some of these questions have already been discussed and answered in the works in [8, 9, 12]. There, the one-particle sector was derived from a field-theoretical model and quantum mechanical master equations were derived for the non-relativistic and ultra-relativistic cases. The new aspect we aim at investigating in this work is the role of renormalisation in this context and its interplay with the further approximations, such as the Markov and rotation wave approximation, that one needs to apply in the derivation of the final master equation.  \nTo the best of the authors' knowledge, renormalisation for field-theoretical models of gravitationally induced decoherence has been performed in the existing literature after applying the Markov  \nand rotational wave approximation, often at the level of the corresponding one-particle sector in certain limits such as the ultra-relativistic limit [12] or the non-relativistic limit [9]. In contrast in this work we will perform the renormalisation at the level of the effective field theory before we apply any of the above mentioned approximations or limits. This strategy allows us to obtain a more detailed understanding about the UV-divergent contributions in the one-particle master equations. We will follow the methods introduced in [21], where one scalar field was coupled to a second scalar field as an environment and extend those techniques to the case of a gravitational environment. These methods allow us to identify individual contributions in the one-particle master equations with specific Feynman diagrams of the underlying effective field-theoretical model. Since the starting point of the model in this paper is the canonical formulation of a scalar field coupled to linearised gravity, in a first step we introduce non-covariant Feynman rules adapted to the canonical model, following [22, 23], where this was introduced for the case of QED. Interestingly, the connection between the covariant and non-covariant Feynman rules can be used to show that the divergent contributions in the one-particle master equation are involved in the self-energy of the scalar field in the effective field theory. The self-energy can be decomposed into a vacuum and a thermal contribution, the latter vanishing when we consider the zero temperature limit of thermal gravitational waves in the environment. The vacuum contributions can then be renormalised using a standard procedure, and the renormalised one-particle master equation can be obtained.  \nEquipped with this result, we can use it and investigate what kind of effect a Markov and rotating wave approximation have and compare it to the existing literature, where these two approximations are mostly done before a renormalisation is preformed. At least for the case of the ultra-relativistic limit, we are able to provide a condition under which the Markov approximation can be applied in the model considered here. For the general case, this is quite a challenging task, since the integrals involved, which have to be analysed for the environmental correlation functions, are quite complicated. In the case of the rotating wave approximation in the existing literature the pre- and post-trace application (see for instance [12, 24-26] for analyses and applications) is discussed and we apply the latter in this work and determine the final one-particle master equation where both the Markov and rotating wave approximation have been applied. This is used together with the intermediate results before the individual approximations in some applications to compare with the existing results in the literature, in particular the work in [9, 11, 12] for field-theoretical models and the quantum mechanical model in [27] in the context of gravitationally induced decoherence in neutrino oscillations, which is based on a quantum mechanical toy model for gravitationally induced decoherence from [28]. For the comparison with the quantum mechanical model, we are especially interested in the extent to which we can relate the application and effect of renormalisation in the field-theoretical model and in the quantum mechanical model and how we can thereby gain new insights into the differences, and similarities respectively as well as the physical properties of these models.  \nThe paper is structured as follows: after the introduction in section I, we briefly review the field-theoretical model from [1] in section II whose one-particle projection is derived in section III, where we consider two types of projections, a non-extended and an extended one, which differ in whether each individual operator in the master equation preserves the one-particle sector or whether only the combinations of operators that enter the final master equation must do so. The renormalisaton of the one-particle master equation is discussed in section IV. First we identify those contributions in the one-particle master equations that are UV-divergent in subsection IV A. Subsection IV B and IV C introduce the non-covariant and covariant Feynman rules and discuss their relation as well as show how the UV-divergent contributions can be identified with the vacuum part of the scalar field's self-energy. The renormalisation of the self-energy is discussed in  \nsubsection IV D and in subsection IV E we use these results to determine the final form of the renormalised one-particle master equation. Afterwards in section V we apply the Markov and rotating wave approximation that are separately discussed in subsection V A and V B respectively. As possible applications in section VI we discuss in subsection VI A the evolution of the populations of the one-particle master equation before and after renormalisation as well as after in addition the Markov approximation has been applied, see subsections VI A 1, VI A 2 and VI A 3 and compare our results with the ones in [11]. To compare more in detail to the existing results in [8, 9, 12] we consider the non-relativistic and ultra-relativistic limit in subsections VI B and VI C. The comparison with quantum mechanical model from [27] that considers gravitationally induced decoherence in the context of neutrino oscillations can be found in subsection VI D. Finally we summarise and conclude in section VII. In addition, details of the calculations required to obtain the results of each section are provided in the appendix to make the article self-contained.", 'II. REVIEW OF THE UNDERLYING FIELD-THEORETICAL MASTER EQUATION': "As we aim at investigating the one-particle sector of model considered in [1] in this work and compare the results obtained here to results in the existing literature such as in [8, 9, 12], we briefly review the main results from [1] that are taken here as a starting point for the further analysis. In [1] a step-by-step derivation of a second order, time-convolutionless master equation for gravitationally induced decoherence of a scalar field is presented with the aim to compute the effective dynamics of the scalar field evolving in an environment consisting of thermal gravitational waves. The starting point in [1] is the classical action of general relativity coupled to a scalar field, where the mostly plus signature is used for the metric. To be able to apply canonical quantisation later on, general relativity is then formulated in the Hamiltonian (ADM) framework ([29]) formulated in terms of Ashtekar-Barbero variables ([30-33]), which encodes the gravitational degrees of freedom in terms of an SU(2)-connection A i a ( ⃗x, t ) and its canonically conjugate momenta denoted as densitised triads E a i ( ⃗x, t ). Here, the first index a is a spatial one and the second one i an SU(2)-Lie algebra index. The metric can be determined by the densitised triads only up to a rotation and this yields to an additional so-called Gauß constraint next to the Hamiltonian and spatial diffeomorphism constraint already present in the ADM formulation.  \nThe reason that [1] works with these elementary variables is that they form the elementary canonical variables in the context Loop Quantum Gravity (LQG, see [34, 35] for an introduction) and hence formulating the classical model in terms of them makes the application of a quantisation using LQG techniques possible, as discussed in [1]. After the inclusion of a boundary term for an asymptotic flat spacetime ([36]), the gravitational sector of the system is linearised around a Minkowski spacetime, where κ = 8 πG N c 4 plays the role of the perturbation parameter and G N denotes Newton's constant. The matter contribution is included in terms of a post-Minkowski approximation [37]. As the linearised model still includes gauge symmetries, suitable Dirac observables are then constructed in perturbation theory using the relational formalism ([38-44]) to extract the physical degrees of freedom. For gravity they correspond to the linearised symmetric transverse traceless components of the connection and the triad fields. These constitute of two degrees of freedom each, equivalent to the two polarisations of gravitational waves and their conjugate momenta which are the linearised excitations one obtains when using the metric and its conjugate momenta as the elementary variable. In addition one can construct Dirac observables for the scalar field and its canonically conjugate momentum for the matter sector.  \nAfter a Fock quantisation, a time-convolutionless master equation ([2]) is derived by tracing out the gravitational environment, which is described by a Gibbs state with temperature 1 parameter  \nΘ, and in this way treating the scalar field as an open quantum system. For this, the projection operator techniques (see [2, 45, 46]) are employed, for which one splits the Hilbert space into a relevant and an irrelevant part. In the model considered here the scalar field is the relevant and the gravitational environment the irrelevant part. Then the strategy is to solve the Liouville-von Neumann equation of the relevant part perturbatively using an approximated, truncated solution for the irrelevant part. The final master equation in [1] is time-convolutionless and truncated after second order in √ κ . It reads  \n∂ ∂t ρ S ( t ) = -i [ H S + κ U + κ H LS , ρ S ( t )] + D [ ρ S ( t )] . (2.1)  \nHere, ρ S ( t ) denotes the density matrix describing the state of the scalar field at physical 2 time t . Furthermore, H S is the Hamiltonian of the free scalar field,  \nH S := ∫ R 3 d 3 k ω k a † k a k , (2.2)  \nwhere a ( † ) k denote the annihilation/creation operator valued distributions for momentum ⃗ k corresponding to the scalar field, that obey the standard commutation relations [ a k , a † l ] = δ 3 ( ⃗ k, ⃗ l ), and ω k = √ ⃗ k 2 + m 2 its dispersion relation. The other terms in (2.1) arise due to the coupling to linearised gravity and represent a self-interaction term U and a Lamb-shift-like term H LS (given in (A.I.2)) as well as a dissipator term D (given in (A.I.3)). The first one arises when the Hamiltonian is expressed in terms of the independent Dirac observables mentioned above and the Lamb-shiftlike term H LS contributes to the unitary evolution of the scalar field, whereas the dissipator leads to to non-unitary evolution encoding effects like decoherence. A similar form for a master equation can be obtained by using the Born approximation instead of the projection operator technique, see section 4 in [1] for more details.  \nWhile this master equation predicts the effective dynamics of the scalar field under the influence of a gravitational environment, we still need to deal with infinitely many degrees of freedom which makes investigating the solution of the master equation very challenging. For this reason in this work the master equation for a single scalar particle that follows from the master equation in (2.1) is derived.  \nIn general, this master equation, as well as the full field theoretical one from [1], might not be completely positive and usually to check this is a challenging task for the field theoretical models. To obtain a completely positive master equation that allows a probability interpretation of the density matrix for all times larger than the initial time, a common procedure is to invoke several approximations. In this work, we also discuss two frequently applied approximations of this kind in detail, which are the Markov and the rotating wave approximation, and investigate their applicability as well as their interplay with the renormalisation of the model.", 'III. PROJECTION OF THE MASTER EQUATION TO THE ONE-PARTICLE CASE': 'In this section, we discuss how such a master equation for a single scalar particle can be obtained starting with equation (2.1) in the field theory context. When working with field theoretical master equations such a one-particle projection is commonly applied to investigate some features of the master equation, see for instance the works in [9, 11, 12, 21]. There exist however different methods  \non how to perform this projection in detail. In [9, 11] the procedure is carried out such that the final one-particle master equation is probability conserving, which requires to neglect some terms that would otherwise be present in a direct projection. In contrast, in [21] a different strategy based on Thermo Field Dynamics (TFD), which is a formulation of Quantum field theory at finite temperature (see [47, 48] and for an introduction [49]), is employed, in which these terms still contribute to the one-particle master equation. Here we follow the method used in [9, 11], but keep all possible terms and investigate their influence on the one-particle master equation. It will turn out that after applying an on-shell renormalisation and Markov approximation, they will not play any role for decoherence, but will remove the remaining contribution of the Lamb-shift-like term to the unitary evolution after the rotating wave approximation has been applied.  \nTo obtain the one-particle projection of the master equation in (2.1), we replace the density matrix with the corresponding density matrix for a single particle in momentum representation  \nρ 1 ( t ) = ∫ R 3 d 3 u ∫ R 3 d 3 v ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v , (3.1)  \nin the master equation and neglect all contributions that project out of the single particle space. In this formulation, ρ ( ⃗u, ⃗v, t ) is the (quantum mechanical) density matrix of a single particle in momentum representation.  \nIn the following we will discuss the corresponding individual contributions in (2.1) separately and further will discuss the assumptions used in the model considered here as well as compare them to the existing literature:  \nThe first term of the master equation, representing the evolution of a free scalar particle, can be computed immediately to yield  \n-i [ H S , ρ 1 ( t )] = -i ∫ R 3 d 3 u ∫ R 3 d 3 v ( ω u -ω v ) ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v . (3.2)  \nThe contribution of the second term depends on the structure of the form of the operator U . The detailed expression is given in [1] in equation (3.12), for the analysis here it is however enough to know that U consists of different combinations of always four creation and/or annihilation operators for the scalar field, i.e.  \nU = ∫ d 2 k 1 ∫ d 2 k 2 ∫ d 2 k 3 ∫ d 2 k 4 [ ζ 1 a k 1 a k 2 a k 3 a k 4 + ζ 2 a † k 1 a k 2 a k 3 a k 4 + ζ 3 a k 1 a † k 2 a k 3 a k 4 + ... + ζ N a † k 1 a † k 2 a † k 3 a † k 4 ] , (3.3)  \nwith the coefficient distributions ζ i = ζ i ( ⃗ k 1 , ⃗ k 2 , ⃗ k 3 , ⃗ k 4 ) that contain delta distributions that relate some of the momenta. For more details, we refer to the definition in equation (3.12) in [1]. When applying normal ordering to this operator, as it is done in [1], then it will not contribute after the one-particle projection: in the summands of U , where the number of creation operators is not equal to the number of annihilation operators, the resulting terms would project out of the one-particle space. In the other summands, there are exactly two creation and two annihilation operators which, when normal ordered, annihilate any one-particle state. In [9, 12] the normal ordering of U is applied differently: in their work the four annihilation and/or creation operators are normal ordered pairwise 3 . In that case contributions of the form : a † k 1 a k 2 : : a † k 3 a k 4 : preserve the one-particle space and thus still contribute after the one-particle projection. To distinguish  \nthese two types of operator orderings, we denote the first one, where U is normal ordered, total normal ordering, and the second one partial normal ordering. In this work we consider a totally normal ordered Hamiltonian as in [1].  \nThe third term, the Lamb-shift-like Hamiltonian, as well as the fourth contribution, the dissipator, both contain the same building blocks. To evaluate them, it is sufficient to consider the following three combinations:  \nj A r ( ⃗ k, ⃗p ) † j B r ( ⃗ k, ⃗ l ) ρ 1 ( t ) ︸ ︷︷ ︸ ( I ) and ρ 1 ( t ) j A r ( ⃗ k, ⃗p ) † j B r ( ⃗ k, ⃗ l ) ︸ ︷︷ ︸ ( II ) and j B r ( ⃗ k, ⃗ l ) ρ 1 ( t ) j A r ( ⃗ k, ⃗p ) † ︸ ︷︷ ︸ ( III ) , (3.4)  \nwhere the j A r ( ⃗ k, ⃗p ) denote individual and different normal-ordered current operators labeled by A ∈ { 1 , 2 , 3 , 4 } that carry a polarisation label r ∈ {±} and two momentum arguments. These current operators are defined in detail starting in equation (A.I.8) and are of the form  \nj 1 r ( ⃗ k, ⃗p ) ∝ a † p a k + p j 2 r ( ⃗ k, ⃗p ) ∝ a † -p -k a -p (3.5)  \nj 3 r ( ⃗ k, ⃗p ) ∝ a -p a k + p j 4 r ( ⃗ k, ⃗p ) ∝ a † p a † -k -p . (3.6)  \nHence they consist of two creation and/or annihilation operators with different momentum labels. At this point arises the question whether we want to enforce trace preservation in the one-particle master equation, which corresponds to probability conservation. In [9, 11] this is done, which results in the exclusion of specific terms from the one-particle master equation. These terms can be identified from the general form of the master equation in (2.1) as we will discuss now: it is evident that when applying the trace the commutator vanishes and one is left with  \n∂ ∂t tr { ρ S ( t ) } = tr {D [ ρ S ( t )] } . (3.7)  \nInserting the definition of the dissipator given in [1] in equation (4.74) then yields  \n∂ ∂t tr { ρ S ( t ) } = κ 2 ∫ d 3 k d 3 p d 3 l (2 π ) 3 ∑ r ∈{±} 4 ∑ A,B =1 R AB ( ⃗ p, ⃗ l ; ⃗ k, t ) · tr { j B r ( ⃗ k, ⃗ l ) ρ S ( t ) j A r ( ⃗ k, ⃗p ) † -1 2 { j A r ( ⃗ k, ⃗p ) † j B r ( ⃗ k, ⃗ l ) , ρ S ( t ) } } , (3.8)  \nwhere the R AB are time-dependent coefficients. When the current operators j A r are individually projected onto the one-particle space, due to the cyclicity of the trace all terms in the difference of the two traces are exactly canceled and one obtains a preserved trace of the density matrix, hence probability conservation. This is the approach used for instance in [9, 11]. Another option is to apply the one-particle projection in such a way that each entire term in the master equation has to preserve the one-particle space. This is for instance done in [21], where two scalar fields are considered, one as the system and the other one as the environment. In this case there will remain terms in the one-particle projection of the product of two current operators in the last term of (3.8) that have no counterpart in the first term of (3.8) and thus will not cancel in the difference of the two traces.  \nTo keep our analysis as generic as possible, we will include these terms in this work and investigate their effect in the one-particle master equation and denote this one-particle projection the extended one-particle projection. To take them into account in our further calculations, we will introduce a  \nfactor δ P in these contributions to be able to switch between the extended one-particle projection ( δ P = 1) and the non-extended one ( δ P = 0).  \nThe detailed derivation of the one-particle projection of the master equation following these methods introduced here can be found in appendix A.I. The additional terms that are present in the extended one-particle projection correspond to physical situations in QFT in which in the intermediate steps two particles are created and annihilated afterwards. This also includes the case where the original particle is left invariant and a vacuum bubble is created. The latter case thus requires a renormalisation, which is also carried out in appendix A.I. Following these projection methods, the one-particle master equation for the density matrix in momentum representation is then given in (A.I.47) and reads  \n∂ ∂t ρ ( ⃗u, ⃗v, t ) = -iρ ( ⃗u, ⃗v, t ) ( ω u -ω v ) -κ 2 ∫ d 3 k (2 π ) 3 { P u ( ⃗ k ) ω u -k ω u [ C ( ⃗u, ⃗ k, t ) + δ P C P ( ⃗u, ⃗ k, t ) ] + P v ( ⃗ k ) ω v -k ω v [ C ∗ ( ⃗v, ⃗ k, t ) + δ P C ∗ P ( ⃗v, ⃗ k, t ) ] } ρ ( ⃗u, ⃗v, t ) + κ 2 ∫ d 3 k (2 π ) 3 P ijln ( ⃗ k ) u i u j v l v n √ ω u + k ω u ω v + k ω v { C ( ⃗u + ⃗ k, ⃗ k, t ) + C ∗ ( ⃗v + ⃗ k, ⃗ k, t ) } ρ ( ⃗u + ⃗ k,⃗v + ⃗ k, t ) (3.9)  \nwith P u ( ⃗ k ) := P ijln ( ⃗ k ) u i u j u l u n which contains the symmetric transverse traceless (STT-)projector  \nP ijln ( ⃗ k ) = 1 2 [ P il ( ⃗ k ) P jn ( ⃗ k ) + P in ( ⃗ k ) P jl ( ⃗ k ) -P ij ( ⃗ k ) P ln ( ⃗ k )] , (3.10)  \nthat in turn consists of combinations of the transverse projectors  \nP ij ( ⃗ k ) = δ ij -k i k j ⃗ k 2 . (3.11)  \nThe presence of this projector is a consequence of the chosen Dirac observables, and thus the physical degrees of freedom of the linearised gravitational field. The coefficients in (3.9) are defined as  \nC ( ⃗u, ⃗ k, t ) = ∫ t -t 0 0 dτ Ω k { [ N ( k ) + 1] e -i (Ω k + ω u -k -ω u ) τ + N ( k ) e i (Ω k -ω u -k + ω u ) τ } (3.12)  \nC P ( ⃗u, ⃗ k, t ) = ∫ t -t 0 0 dτ Ω k { [ N ( k ) + 1] e -i (Ω k + ω u -k + ω u ) τ + N ( k ) e i (Ω k -ω u -k -ω u ) τ } , (3.13)  \nwhere N ( k ) = 1 e β Ω k -1 is the Bose-Einstein distribution of the gravitational waves in the environment with frequencies Ω k := √ ⃗ k 2 and β = ( k B Θ) -1 , where k B is the Boltzmann constant and Θ the temperature parameter of the Gibbs state that characterises the environment of thermal gravitational waves. The term in the first line of the master equation in (3.9) represents the standard unitary evolution of the free scalar particle. The remaining terms describe the influence of the environment and encode in general different physical processes like energy shifts, dissipation and decoherence. While the expressions in lines two and three only depend on the state ρ ( ⃗u, ⃗v, t ) considered at time t , the last line links this state to other states ρ ( ⃗u + ⃗ k,⃗v + ⃗ k, t ) at time t . This master equation however still has some problematic contributions that possess UV-divergences and hence needs to be UV-renormalised, as will be discussed in the next section.', 'IV. RENORMALISATION OF THE TCL ONE-PARTICLE MASTER EQUATION': 'Upon investigation of the individual contributions in (3.9) it becomes evident that some terms exhibit divergences as will be discussed in detail below. This raises the question of at what stage of the derivation of the master equation the renormalisation procedure should be carried out. In the literature, there are different strategies how to deal with these divergent contributions. For gravitationally induced decoherence, they have often not been computed in detail due to the reason that they are expected to not modify decoherence but only influence the unitary evolution, see for instance the discussion in [8, 11]. In [9, 12] the renormalisation of these contributions has been performed in the end after a Markov and rotating wave approximation have already been applied.  \nIn this work we choose the strategy to renormalise the master equation first before applying further approximations or deriving physical implications. It will turn out that effects predicted with a non-renormalised master equation might get modified or even vanish when working with the renormalised version instead. An example of this kind is also discussed in [27], where a quantum mechanical toy model for gravitationally induced decoherence based on the model in [28] is applied in the context of neutrino oscillations. In that case the necessary renormalisation is very trivial compared to the model considered in this work. There, the renormalisation causes the contributions of the Lamb-shift Hamiltonian to cancel exactly. Consequently, all physical implications involving contributions of the Lamb-shift Hamiltonian, as discussed for example in [50], would be absent when working with the renormalised model presented in [27].  \nIn order to carry out the renormalisation, we will first identify the diverging terms. As we will discuss below in more detail, these are in particular the terms in the second and third line of the master equation (3.9) that will also be present in the case where the temperature parameter vanishes, that is for Θ = 0, in which the thermal state merges into a vacuum state. They are of the form ∫ d 3 k 1 | ⃗ k | 3 and thus yield a logarithmic UV-divergence. Once these contributions are identified, we express them in the form of Feynman diagrams of the underlying effective QFT. For this purpose, we follow the strategy in [21], where a master equation for a scalar field with an environment consisting of another scalar field is presented. Here the treatment is extended so that the linearised gravitational field can be included as an environment. We will proceed in five steps: first in subsection IV A we will identify the divergent contributions in the master equation and then present the corresponding Feynman rules following from the underlying effective QFT based on the non-covariant formulation in subsection IV B. Afterwards in subsection IV C we provide a set of equivalent, covariant Feynman rules in terms of which we perform the renormalisation of the divergent contributions in subsection IV D. Finally we discuss the resulting renormalised one-particle master equation in subsection IV E.', 'A. Identification of the UV divergences in the one-particle master equation': 'Starting from the master equation in (3.9), we want to investigate which terms on the righthand side are UV-divergent with respect to the ∫ R 3 d 3 k integration. As the projector P ijln ( ⃗ k ) is independent of the absolute value of ⃗ k , it does not influence the UV behaviour. Then, one can identify four different types of contributions in the integrands after performing the τ -integration and introducing the following sign factors σ, σ 1 , σ 2 ∈ {± 1 } :  \n- (a) 1 ω u -k Ω k 1 Ω k + ω u -k + σω u : for large | ⃗ k | , i.e. for | ⃗ k | >> | ⃗u | , m this term becomes 1 | ⃗ k | 3 and thus leads to a logarithmic UV-divergence under the integral.\n- (b) 1 ω u -k Ω k 1 Ω k + ω u -k + σω u ρ ( ⃗u + ⃗ k,⃗v + ⃗ k, t ): assuming that ρ ( ⃗x, ⃗y, t ) is a proper, normalisable density matrix in position space for which the Fourier transform exists leads to the requirement that ρ ( ⃗u, ⃗v, t ) has to decrease rapidly for large ⃗u, ⃗v . Therefore this expression is UV-finite.\n- (c) N ( k ) 1 ω u -k Ω k 1 Ω k + σ 1 ω u -k + σ 2 ω u : a series expansion of the denominator of N ( k ) yields N ( k ) = 1 β | ⃗ k | + O ( | ⃗ k | 2 ) . Hence, this term tends to zero for large | ⃗ k | and also the combination x n e x -1 decreases rapidly for x →∞ , thus this kind of contribution is UV-finite.\n- (d) N ( k ) 1 ω u -k Ω k 1 Ω k + σ 1 ω u -k + σ 2 ω u ρ ( ⃗u + ⃗ k,⃗v + ⃗ k, t ): this contribution is a combination of cases (b) and (c) and also UV-finite.  \nFrom this analysis follows that the expressions involving N ( k ), that would be absent in the vacuum case and are thus denoted as thermal contributions in the following, are all UV-finite. Some of the vacuum contributions, these are the ones that do not involve N ( k ), lead to UV-divergences which we want to cure by a renormalisation. To achieve this, in the next section we show in a first step that these terms correspond exactly to the self-energy diagrams for the scalar particle in the form of Feynman diagrams.', 'B. Non-covariant Feynman rules and self-energy': "In this section we present the Feynman rules in non-covariant form corresponding to the effective quantum field theory containing a scalar field coupled to linearised gravity in [1], where the latter is considered as the environment, which is the basis for the master equation in (3.9). Then, we rewrite these rules in the next section in a covariant form to be able to follow the strategy of [21], where a suitable renormalisation for a master equation for a scalar fields with a second scalar field as the environment is discussed. Here, we slightly extend these methods in order to apply them to the case where the linearised gravitational field is treated as the environment. The Feynman rules can be constructed from [1]:  \n- · The scalar field has the standard propagator , which follows from its quantised mode expansion in (3.4) and (3.5) in [1], which we denote by a solid line and which reads in the mostly plus signature convention:  \n= -i k 2 + m 2 -iϵ . (4.1)  \nHere, k 2 = -( k 0 ) 2 + ⃗ k 2 .  \n- · The propagator of the triad field was derived in [1] in equation (4.46) and is denoted in terms of a curved line:  \n= 1 κ P abcd ( ⃗ k ) [ -i k 2 -iϵ +2 πN ( k ) δ ( k 2 ) ] , (4.2)  \nwhere the first summand is the vacuum and the second one the thermal part. The existing tensor structure manifests in the form of the tensor structure of the STT-projector defined  \nabove in (3.10). When contracted with a quantity that is symmetric in ( ab ) as well as in ( cd ), like it is the case for the interaction vertex introduced below, the STT-projector reduces to  \nP abcd ( ⃗ k ) = δ ac δ bd -1 2 δ ab δ cd + 1 2 ⃗ k 4 k a k b k c k d -1 ⃗ k 2 ( 2 δ ac k b k d -1 2 δ ab k c k d -1 2 δ cd k a k b )] . (4.3)  \n- · The coupling between the scalar field and linearised gravity is encoded in the interaction part of the total action in [1] that is given as a reformulation of equation (3.17) in that work by:  \nS int = -∫ dt ∫ d 3 x H int ( ⃗x, t ) = -∫ dt ∫ d 3 k ⃗ k 2 κ { ⃗ k 2 δ E ab ( ⃗ k, t ) T ab ( -⃗ k, t ) -1 2 T 00 ( ⃗ k, t ) T a a ( -⃗ k, t ) -1 4 T 00 ( ⃗ k, t ) T 00 ( -⃗ k, t ) + T 0 a ( ⃗ k, t ) T 0 b ( -⃗ k, t ) ( δ ab -k a k b 4 ⃗ k 2 )} , (4.4)  \nwhere T µν ( ⃗ k, t ) denotes the Fourier transform of the scalar field's energy momentum tensor 4 . The interaction vertex between the scalar and the triad field can then be read off and is related to T ab . Due to the fact that T ab depends on derivatives of the scalar field, the expression for the triad-scalar-field vertex is different depending on the direction of the momenta involved in the diagrams. Considering the Fourier transform of T ab ( ⃗x, t ), where the the scalar fields can be factorised, we find the expression ˜ T ab ( p, q ) ϕ ( p ) ϕ ( q ) with  \n˜ T ab ( p, q ) := 1 2 δ ab [ -p 0 q 0 + ⃗ p · ⃗q -m 2 ] -1 2 [ p a q b + p b q a ] . (4.5)  \nHence, the triad-scalar-field-vertex is given by  \np q = p q = iκ ˜ T ab ( σ p p, σ q q ) ,  \nwhere σ p is + if particle p is incoming and -if it is outgoing.  \n- · The remaining terms in the interaction part of the action (4.4) give rise to an additional second order vertex which cannot be split into first order vertices due to the lack of a suitable intermediate particle in this effective field theory. They have the form  \n-κ ∫ d 3 x ∫ d 3 y -1 4 T 00 ( ⃗x ) T 00 ( ⃗y ) + T 0 a ( ⃗x ) T 0 b ( ⃗y ) δ ab -1 4 ( T 00 ( ⃗x ) T a a ( ⃗y ) + T a a ( ⃗x ) T 00 ( ⃗y )) 4 π || ⃗x -⃗y || -κ ∫ d 3 x ∫ d 3 y ∫ d 3 z -1 4 ∂ a ⃗x T 0 a ( ⃗x ) ∂ b ⃗y T 0 b ( ⃗y ) (4 π ) 2 || ⃗x -⃗z |||| ⃗y -⃗z || . (4.6)  \nT 00 = T 00 = 1 2 [ π 2 +( ∂ a ϕ ) ( ∂ a ϕ ) + m 2 ϕ 2 ] = ϵ ( ϕ, π ) , T 0 a = T a 0 = -δ ab T 0 b = -δ ab π ∂ b ϕ = -δ ab p b ( ϕ, π ) , T ab = δ ac δ bd T cd = δ ab [ π 2 -( ∂ c ϕ ) ( ∂ c ϕ ) -m 2 ϕ 2 ] +( ∂ a ϕ ) ( ∂ b ϕ ) .  \n2  \nThe corresponding symmetrised Feynman rule reads:  \nk p q u v = -iκ ⃗ k 2 NI ( p, q, u, v )  \nwith  \nNI ( p, q, u, v ) := [ -1 4 ˜ T 00 ( p, q ) ˜ T 00 ( -u, -v ) + ˜ T 0 a ( p, q ) ˜ T 0 b ( -u, -v ) ( δ ab -k a k b 4 ⃗ k 2 ) -1 4 ˜ T 00 ( p, q ) ˜ T a a ( -u, -v ) -1 4 ˜ T a a ( p, q ) ˜ T 00 ( -u, -v ) ] , (4.7)  \nwith ⃗ k defined using momentum conservation as ⃗ k = ⃗ p + ⃗q = ⃗u + ⃗v . Note that a similar vertex does also appear in QED when quantised in Coulomb gauge and there it represents the Coulomb interaction, see for instance [22, 23].  \n- · As external lines we only have the scalar field in the cases we are interested in here. This follows the standard case for a quantised scalar field, the detailed expressions are however not required for the following discussions.  \nEquipped with these Feynman rules, we will now show in the subsequent section that the divergent contributions in the one-particle master equation can be identified exactly with the contribution of the self-energy diagram constructed with the above Feynman rules. In the model considered here this corresponds to the following Feynman diagram:  \nu u -k u k (4.8)  \nWith the Feynman rules introduced above, the amplitude represented by this diagram and denoted by Π( u 2 ) has the form  \nΠ( u 2 ) := ∫ R 4 d 4 k (2 π ) 4 1 κ P abcd ( ⃗ k ) [ -i k 2 -iϵ +2 πN ( k ) δ ( k 2 ) ] [ iκ ˜ T ab ( u, -( u -k )) ] · -i ( u -k ) 2 + m 2 -iϵ [ iκ ˜ T cd ( u -k, -u ) ] = κu a u b u c u d ∫ R 4 d 4 k (2 π ) 4 P abcd ( ⃗ k ) [ -i k 2 -iϵ +2 πN ( k ) δ ( k 2 ) ] 1 ( u -k ) 2 + m 2 -iϵ = κ ∫ R 3 d 3 k 2(2 π ) 3 P u ( ⃗ k ) { (Ω k + ω u -k ) Ω k ω u -k (Ω k + ω u -k + u 0 -iϵ )(Ω k + ω u -k -u 0 -iϵ ) -N (Ω k ) [ 1 Ω k ( u 0 -Ω k + ω u -k -iϵ )( u 0 -Ω k -ω u -k + iϵ ) + 1 Ω k ( u 0 +Ω k + ω u -k -iϵ )( u 0 +Ω k -ω u -k + iϵ ) ]} =Π vac ( u 2 ) + Π Θ ( u 2 ) (4.9)  \nIn the first step the definition of ˜ T ab in (4.5) was used and in the second step the k 0 -integration was performed, where for the vacuum part the residue theorem was applied. In the last step, we have defined the vacuum and thermal contribution to the self-energy as  \nΠ vac ( u 2 ) := κ ∫ R 3 d 3 k 2(2 π ) 3 P u ( ⃗ k ) { (Ω k + ω u -k ) Ω k ω u -k (Ω k + ω u -k + u 0 -iϵ )(Ω k + ω u -k -u 0 -iϵ )  \nΠ Θ ( u 2 ) := -κ ∫ R 3 d 3 k 2(2 π ) 3 P u ( ⃗ k ) N (Ω k ) 1 Ω k ( u 0 -Ω k + ω u -k -iϵ )( u 0 -Ω k -ω u -k + iϵ  \n(4.10) [ ) + 1 Ω k ( u 0 +Ω k + ω u -k -iϵ )( u 0 +Ω k -ω u -k + iϵ ) ]} . (4.11)  \nIf we now want to identify contributions in the one-particle master equation with the self-energy, the following subtlety results: a key difference between the master equation in (3.9) and standard quantum field theory is that the latter is constructed for the limit t 0 →-∞ , t →∞ when evaluating scattering amplitudes. To take this into account, we apply the method presented in [21] for situations where there is a finite temporal interval. In this way, we can transform the self-energy diagram into the second line of the right-hand side of the master equation in (3.9):  \nΞ( ω u , ⃗u, t 0 , t ) := ∫ t t 0 dτ ∫ R du 0 Π( u 2 ) cos[( u 0 -ω u )( t -τ )] = ∫ t -t 0 0 dτ ∫ R du 0 Π( u 2 ) cos[( u 0 -ω u ) τ ] . (4.12)  \nThe standard QFT-limit can be recovered, in which t → ∞ and t 0 → -∞ , and using this the integral over the temporal interval can be rewritten as a δ -distribution as ∫ ∞ 0 dτ cos[( u 0 -ω u ) τ ] = 1 2 ∫ ∞ -∞ dτe -i ( u 0 -ω u ) τ = πδ ( u 0 -ω u ), that will set the external momentum u on-shell.  \nAfter evaluating the u 0 integration we obtain for finite times t and t 0 :  \nΞ( ω u , ⃗u, t 0 , t ) = κ 2 ∫ R 3 d 3 k (2 π ) 3 P u ( ⃗ k ) πi 2 ω u -k [ C ( ⃗u, ⃗ k, t -t 0 ) + C P ( ⃗u, ⃗ k, t -t 0 ) ] , (4.13)  \nwhich indeed, multiplied by a factor 2 i πω u , can be identified with the first term on the right hand side of the second line of the one-particle master equation in (3.9) in the extended one-particle projection. Given this results, it is now also easy to discuss the case of the non-extended one-particle projection: the master equation for this case, i.e. for δ P = 0, can just be obtained by replacing the cosine in (4.12) by 1 2 e -i ( u 0 -ω u ) τ . We find that in the QFT-limit the difference between the extended and non-extended one-particle projection manifests itself in a factor of 2.  \nTo obtain the second term in the second line of (3.9), we can follow the same steps and just have to replace ⃗u by ⃗v and take the complex conjugate.  \nWith the results in this section we have shown that the UV-divergent terms in the one-particle master equation correspond to the self-energies of the scalar particle. What remains is to discuss the renormalisation of this self-energy. In order to be able to apply the standard procedure for renormalisation in this case, however, the corresponding covariant Feynman rules of the model considered here must first be derived.", 'C. Covariant Feynman rules': "To be able to employ the standard renormalisation technique for the loop associated with the scalar particle's self-energy, we introduce in this section the covariant Feynman rules corresponding to the effective QFT under consideration here. For this, we follow [22, 23], where the procedure is outlined for QED.  \nIn a first step we will demonstrate that specific sums of non-covariant Feynman diagrams add up to the corresponding covariant Feynman diagram. For this purpose, we consider the sum of the second order vertex in (4.6) with a second order combination of the non-covariant scalar field-triad vertex, shown below (4.5). As will be derived below, the second order vertex in (4.6) is precisely that term which restores covariance if we work with a fully covariant triad propagator and a covariant vertex.  \nWe will restrict our discussions mainly to a Coulomb-scattering type of diagram here which is sufficient for our later applications. At the end of this section we will also briefly discuss the diagram associated with the scalar particle's self-energy. In the case of the Coulomb-scattering type diagram, the above mentioned equivalence in terms of Feynman diagrams reads  \n<!-- image -->  \nwhere the curly line corresponds to the covariant triad propagator.  \nNext, we will present the covariant Feynman rules, then specialise them to the case of the Coulombscattering type diagram to show the above equivalence. The corresponding covariant Feynman rules for the propagators and vertices discussed in the last section are as follows:  \n- · The scalar propagator remains unchanged  \n= -i k 2 + m 2 -iϵ . (4.14)  \n- · The covariant triad propagator becomes  \n= 1 κ -i k 2 -iϵ P µνρσ (4.15)  \nwith  \nP µνρσ := 1 2 [ η µρ η νσ + η µσ η νρ -η µν η ρσ ] . (4.16)  \nIn the context of a linearised gravitational environment there is no multi-triad vertex and therefore in the effective QFT considered in this work the triad propagator always couples only to the scalar field-triad-vertex, see also below. The latter is symmetric in ( µν ) as well as in ( ρσ ). This allows us to slightly simplify the projector P µνρσ whenever it occurs in combination with scalar field-triad-vertices and express it as  \nP µνρσ := η µρ η νσ -1 2 η µν η ρσ . (4.17)  \n- · The covariant vertex is given by  \np q =  \np q = -iκ 2 ( -2 ˜ T µν ( σ p p, σ q q ) + η µν ˜ T ρ ρ ( σ p p, σ q q )) ,  \nwhere  \n˜ T µν ( p, q ) = 1 2 η µν ( p ρ q ρ -m 2 ) -1 2 ( p µ q ν + p ν q µ ) . (4.18)  \nWhenever this is combined with a triad propagator, the second term of the vertex contribution vanishes because we have  \nη µν ˜ T α α [ η µρ η νσ + η µσ η νρ -η µν η ρσ ] = 2 η ρσ ˜ T α α -2 η ρσ ˜ T α α = 0 , (4.19)  \nwhere we used that η µν η µν = 4. Hence, for processes like Coulomb scattering we can replace the expression for the vertex with  \niκ ˜ T µν ( σ p p, σ q q ) . (4.20)  \n- · Since the second order vertex was used to obtain a covariant propagator and vertex, there is no analogue of the second order vertex in the covariant case.\n- · The external lines for the scalar field remain unmodified and the ones for the triad field are not important for this work.  \nThese covariant Feynman rules are in accordance with the ones 5 presented in [51, 52], where also a scalar field is coupled to a linearised gravitational field, except for the usual differences caused by the choice of different signatures for the metric, as they use the mostly minus signature. Using momentum conservation, i.e. k = p + q = u + v , yields  \nk µ ˜ T µ 0 ( p, q ) = k 0 ˜ T 00 ( p, q ) + k a ˜ T a 0 ( p, q ) = 1 2 [ q 0 ( p 2 0 -⃗ p 2 -m 2 ) + p 0 ( q 2 0 -⃗q 2 -m 2 )] , (4.21)  \nk µ ˜ T µa ( p, q ) = k 0 ˜ T 0 a ( p, q ) + k b ˜ T ba ( p, q ) = 1 2 [ q a ( p 2 0 -⃗ p 2 -m 2 ) + p a ( q 2 0 -⃗q 2 -m 2 )] . (4.22)  \nIf the scalar field is on-shell, which we assume for a moment for the Coulomb-scattering diagram, then the right hand side of both expressions vanishes. With this, one can directly show the equivalence of using the covariant set of Feynman rules for the Coulomb scattering diagrams discussed above of this. We do present this in appendix A.II.  \nNote that in perturbation series, the diagram containing the two vertices is of second order in the expansion and hence obtains an additional factor 1 2 compared to the second order vertex diagram.  \nAs a next step we would like to show a similar equivalence between the non-covariant and covariant  \nFeynman rules for the self-energy diagram, namely the following equivalence in terms of Feynman diagrams  \n+ =  \nNote that the diverging term actually only contains the first of the two Feynman diagrams on the left hand side. However, the second term is the self-energy contribution which vanishes in the one-particle projection of the master equation, so we can add this diagram as in the one-particle master equation its contribution vanishes for normal ordering.  \nThe equivalence for the Feynman diagrams on both sides of this equation is however much more difficult to prove compared to the Coulomb scattering diagram, which is also the case in QED, since the momentum inside the loop is not on-shell, which prevents a similar calculations as done for the Coulomb scattering tree level graph. Given that, to prove this equivalence goes beyond the scope of this work here and we refer here to the fact that the covariant set of Feynman rules can also be derived from the same underlying action using a different approach and gauge, which are then used for instance in [51, 52]. Hence, independently of the derivation, we expect that they describe the same physics. Based on this, it is now possible to specify the expression corresponding to the scalar particle self-energy diagram in covariant form and renormalise it, which will be discussed in detail in the next section. As mentioned at the beginning of this section, such a replacement of non-covariant by covariant Feynman rules along the lines presented here is also employed in QED when quantising in Coulomb gauge, as for instance in [22, 23].", 'D. UV-renormalisation of the self-energy of the scalar particle': 'In terms of the Feynman rules introduced in the previous section, the self-energy diagram for the vacuum propagator, which was defined in (4.10), can be expressed as  \nΠ vac ( u 2 ) = ∫ d 4 k (2 π ) 4 [ iκ ˜ T µν ( u, -( u -k )) ] [ iκ ˜ T ρσ ( u -k, -u ) ] 1 κ -i k 2 -iϵ · 1 2 [ η µρ η νσ + η µσ η νρ -η µν η ρσ ] -i ( u -k ) 2 + m 2 -iϵ = κ 2 ∫ d 4 k u 2 k 2 +2 m 2 uk -2 m 4 [( k + u ) 2 -iϵ ][ k 2 + m 2 -iϵ ] . (4.23)  \nAs the thermal part Π Θ ( u 2 ), defined in (4.11), is not divergent, as it has been discussed in subsection IVA, we only consider the vacuum part here. For the renormalisation we follow the strategy in [53]. Using dimensional regularisation with d = 4 -ϵ , the STT-projector is slightly modified in d dimensions and reads (see e.g. [52]):  \nP µνρσ := 1 2 [ η µρ η νσ + η µσ η νρ -2 2 -ϵ η µν η ρσ ] . (4.24)  \nDue to this, the expression for the self-energy diagram slightly changes and becomes (for ϵ ≪ 1):  \nΠ vac ( u 2 ) = κ 2 µ ϵ ∫ d d k u 2 k 2 +2 m 2 uk -2 m 4 ( 1 + ϵ 4 ) [( k + u ) 2 -iϵ ][ k 2 + m 2 -iϵ ] , (4.25)  \nwhich coincides with the expression derived in [52]. Here, we rescaled κ → κµ ϵ to keep the dimension of κ for any value of d . As later we will encounter also IR-divergences, we introduce a small artificial triad mass λ in the triad propagator that becomes  \n1 κ -i k 2 + λ 2 -iϵ P µνρσ . (4.26)  \nWith this, the self-energy diagram reads  \nΠ vac ( u 2 ) = κ 2 µ ϵ ∫ d d k u 2 k 2 +2 m 2 uk -2 m 4 ( 1 + ϵ 4 ) [( k + u ) 2 + λ 2 -iϵ ][ k 2 + m 2 -iϵ ] , (4.27)  \nThis can then be evaluated using the standard methods for dimensional regularisation (see appendix A.III). The result is that the divergent part can be isolated such that one obtains  \nΠ vac ( u 2 ) = -2 π 2 κm 2 ϵ ( m 2 + u 2 ) + Π reg vac ( u 2 ) (4.28)  \nwith the finite part Π reg vac ( u 2 ). The infinite part then has to be renormalised by introducing a suitable counter term. As the finite part of this counter term can in principle be chosen arbitrarily, Π reg vac ( u 2 ) can still change. In our case, we choose the finite part of the counter term according to the on-shell renormalisation procedure. This then yields for the final renormalised loop Π R vac ( u 2 ):  \nΠ R vac ( u 2 ) := Π reg vac ( u 2 ) -Π reg vac ( -m 2 ) -( u 2 + m 2 ) ∂ ∂u 2 Π reg vac ( -m 2 ) . (4.29)  \nThis specific form is determined by the on-shell renormalisation scheme that sets the pole of the scalar propagator to m 2 and also fixes its residue according to the following two conditions:  \nΠ R vac ( u 2 = -m 2 ) ! = 0 (4.30)  \n∂ ∂u 2 Π R vac ( u 2 = -m 2 ) ! = 0 . (4.31)  \nIt can readily be seen that the definition in (4.29) satisfies these two conditions. Note that in [52], they apply a similar procedure without fixing the residue of the pole and therefore also not including an artificial triad mass, because for their purposes it is sufficient to fix the pole of the propagator.  \nThe consideration above suggests that we have to include the following counter term:  \nδ Π( u 2 ) = 2 π 2 κm 2 ϵ ( m 2 + u 2 ) -Π reg vac ( -m 2 ) -( u 2 + m 2 ) ∂ ∂u 2 Π reg vac ( -m 2 ) (4.32)  \nsuch that  \nΠ( u 2 ) + δ Π( u 2 ) = Π R vac ( u 2 ) + Π Θ ( u 2 ) , (4.33)  \nwhere Π Θ ( u 2 ) denotes the finite thermal contribution to the loop defined in (4.11). From (A.III.19) we have for λ → 0:  \nΠ reg vac ( -m 2 ) = 0 , (4.34)  \nthus  \nδ Π( u 2 ) = [ 2 π 2 κm 2 ϵ -∂ ∂u 2 Π reg vac ( -m 2 ) ] ( m 2 + u 2 ) , (4.35)  \nwhere the expression in the square brackets only depends on m 2 . In order to implement a suitable counter term, we introduce a renormalised mass m R by m 2 = m 2 R + m 2 R δ m , where m denotes the bare mass we have used so far and δ m a mass counterterm, as well as a renormalised wave function  \nφ R = 1 √ Z 2 φ . Then the renormalised scalar field propagator (a Greens function containing twice φ ) reads up to the one-loop contribution:  \niG (1) ( u 2 ) = 1 Z 2 -i u 2 + m 2 = -i u 2 + m 2 R + -i u 2 + m 2 R [ -i ( u 2 δ 2 + m 2 R ( δ 2 + δ m ))+Π( u 2 )] -i u 2 + m 2 R + O ( κ 2 ) , (4.36)  \nwhere we expanded Z 2 = 1 + δ 2 . From this follows that  \n-i ( u 2 δ 2 + m 2 R ( δ 2 + δ m )) ! = δ Π( u 2 ) (4.37)  \nand it becomes evident that only the wave function has to be renormalised in the following manner:  \nδ m = 0 (4.38)  \nδ 2 = i [ 2 π 2 κm 2 R ϵ -∂ ∂u 2 Π reg ( -m 2 R ) ] .  \nDue to the counter-term, We have to replace in the old set of Feynman rules m by m R and obtain the following additional Feynman rule of order κ :  \n= -iδ 2 ( u 2 + m 2 R ) . (4.39)  \nTo simplify notation, we continue to use m , in particular as we have seen that m R = m . The result is therefore an additional counterterm in the Lagrangian which leads in renormalised perturbation theory to an additional interaction of order κ that we have to include when evaluating the loop. Then the former diverging term Π( u 2 ) becomes finite and only Π R vac ( u 2 ) is left. This yields a modification of the right hand side of the master equation. As calculated in appendix A.III.2, the contribution of the renormalised vacuum loop terms to the master equation vanishes:  \nΞ R ( ω u , ⃗u, t 0 , t ) = ∫ t t 0 dτ ∫ R du 0 Π R vac ( u 2 ) cos[( u 0 -ω u )( t -τ )] = 0 , (4.40)  \nΞ R ( ω v , ⃗v, t 0 , t ) = ∫ t t 0 dτ ∫ R dv 0 Π R vac ( v 2 ) cos[( v 0 -ω v )( t -τ )] = 0 . (4.41)  \nTherefore neither the dimensional constant µ , nor the artificial triad mass λ play a role in the physical predictions made with the master equation.', 'E. Renormalised one-particle master equation': "With the renormalisation carried out in the previous subsections, a first renormalised version of the one-particle master equation (3.9), where only the former diverging terms are modified, reads:  \n∂ ∂t ρ ( ⃗u, ⃗v, t ) = -iρ ( ⃗u, ⃗v, t ) ( ω u -ω v ) -κ 2 ∫ d 3 k (2 π ) 3 { P u ( ⃗ k ) ω u -k ω u [ C R ( ⃗u, ⃗ k, t ) + δ P C R P ( ⃗u, ⃗ k, t ) ] + P v ( ⃗ k ) ω v -k ω v [( C R ( ⃗v, ⃗ k, t ) ) ∗ + δ P ( C R P ( ⃗v, ⃗ k, t ) ) ∗ ] } ρ ( ⃗u, ⃗v, t ) + κ 2 ∫ d 3 k (2 π ) 3 P ijln ( ⃗ k ) u i u j v l v n √ ω u + k ω u ω v + k ω v { C ( ⃗u + ⃗ k, ⃗ k, t ) + C ∗ ( ⃗v + ⃗ k, ⃗ k, t ) } ρ ( ⃗u + ⃗ k,⃗v + ⃗ k, t ) (4.42)  \nwith  \nand  \nC ( ⃗u, ⃗ k, t ) = ∫ t -t 0 0 dτ Ω k { [ N ( k ) + 1] e -i (Ω k + ω u -k -ω u ) τ + N ( k ) e i (Ω k -ω u -k + ω u ) τ } . (4.45)  \nAt the level of the operator equation, the renormalisation removed the Θ-independent terms from the terms in the second and third line of (4.42), hence leaving us with the following dissipator:  \nD [ ρ S ] = -κ 2 ∑ r ∈{ + , -} 4 ∑ a,b =1 ∫ R 3 d 3 k d 3 p d 3 l (2 π ) 3 1 Ω k { -( j b r ( ⃗ k, ⃗ l ) ρ S ( t ) j a r ( ⃗ k, ⃗p ) † ) f (Ω k + ω b ( ⃗ k, ⃗ l )) + h.c. + N ( k ) [ j a r ( ⃗ k, ⃗p ) † , [ j b r ( ⃗ k, ⃗ l ) , ρ S ( t ) ]] f (Ω k + ω b ( ⃗ k, ⃗ l )) + h.c. } . (4.46)  \nIf working with the non-extended projection δ P = 0, then there was probability conservation before the renormalisation, i.e. ∫ d 3 u ∂ ∂t ρ ( ⃗u, ⃗u, t ) = 0. Now, due to the vacuum term in (4.45) this probability conservation is destroyed. As the renormalisation is a purely technical procedure that should not change the physics, in particular not basic principles as probability conservation, we also replace C ( ⃗u, ⃗ k, t ) by C R ( ⃗u, ⃗ k, t ) in the last line of the master equation. Another reason for this is that the term in the last line of the master equation is based on the same QFT as the terms in the second and third line, hence they should be renormalised in the same way 6 . The final renormalised one-particle master equation is thus  \n∂ ∂t ρ ( ⃗u, ⃗v, t ) = -iρ ( ⃗u, ⃗v, t ) ( ω u -ω v ) -κ 2 ∫ d 3 k (2 π ) 3 { P u ( ⃗ k ) ω u -k ω u [ C R ( ⃗u, ⃗ k, t ) + δ P C R P ( ⃗u, ⃗ k, t ) ] + P v ( ⃗ k ) ω v -k ω v [( C R ( ⃗v, ⃗ k, t ) ) ∗ + δ P ( C R P ( ⃗v, ⃗ k, t ) ) ∗ ] } ρ ( ⃗u, ⃗v, t ) + κ 2 ∫ d 3 k (2 π ) 3 P ijln ( ⃗ k ) u i u j v l v n √ ω u + k ω u ω v + k ω v { C R ( ⃗u + ⃗ k, ⃗ k, t ) + ( C R ( ⃗v + ⃗ k, ⃗ k, t ) ) ∗ } ρ ( ⃗u + ⃗ k,⃗v + ⃗ k, t ) (4.47)  \nand the dissipator at operator level  \nD [ ρ S ] = -κ 2 ∑ r ∈{ + , -} 4 ∑ a,b =1 ∫ R 3 d 3 k d 3 p d 3 l (2 π ) 3 N ( k ) Ω k { [ j a r ( ⃗ k, ⃗p ) † , [ j b r ( ⃗ k, ⃗ l ) , ρ S ( t ) ]] f (Ω k + ω b ( ⃗ k, ⃗ l )) + h.c. } . (4.48)  \nC R ( ⃗u, ⃗ k, t ) = 2 ∫ t -t 0 0 dτ Ω k N ( k ) cos[Ω k τ ] e -i ( ω u -k -ω u ) τ (4.43)  \nC R P ( ⃗u, ⃗ k, t ) = 2 ∫ t -t 0 0 dτ Ω k N ( k ) cos[Ω k τ ] e -i ( ω u -k + ω u ) τ (4.44)  \nCompared to the expression before renormalisation given in [1] in equation (4.60), the vacuum contribution in the dissipator vanishes.  \nIn [27] a quantum mechanical model based on the model in [28] is considered where a system is coupled to an environment of harmonic oscillators. The bath of harmonic oscillators mimics the thermal gravitational waves and the model serves as toy model for gravitationally induced decoherence. In [27] the system was then specified to neutrinos in order to investigate gravitationally induced decoherence in the context of neutrino oscillations. A more detailed discussion on the model from [27] can be found in section VI D.  \nThe quantum mechanical master equation in that work consists of a dissipator term and a Lambshift, where the latter contained divergences and is finally removed by a renormalisation in [27]. We are now interested to discuss the similarities and differences of that quantum mechanical renormalisation and the renormalisation of the one-particle master equation presented in this work. In the quantum mechanical case, the interaction Hamiltonian has the following form:  \nˆ H I = ˆ H S ⊗ N ∑ i =1 g i ˆ q i (4.49)  \nwith the Hamiltonian ˆ H S , describing the neutrino propagation, coupling constants g i and the position operators of the harmonic oscillators in the environment ˆ q i . This provides a toy model for the coupling of the energy momentum tensor and the metric perturbations in the field theoretical model. The form of the coupling in the quantum mechanical model implies that the coefficients Λ of the Lamb-shift and Γ of the dissipator in the final Lindblad equation only depend on the environment and are defined in the following way, see [27] equations (4) and (5):  \nΛ( t -t 0 ) := ∫ t -t 0 0 dτ ∫ ∞ 0 dω J ( ω ; Ω , η ) sin( ωτ ) (4.50)  \nΓ( t -t 0 ) := 2 ∫ t -t 0 0 dτ ∫ ∞ 0 dω J ( ω ; Ω , η ) cos( ωτ ) coth ( β ℏ ω 2 ) , (4.51)  \nwhere a spectral density J ( ω ; Ω , η ) was used which depends on the frequencies of the harmonic oscillators denoted by ω , an effective coupling parameter η and a UV-cutoff frequency Ω. From this form it becomes evident that the Lamb-shift term is independent of the temperature parameter β = 1 k B Θ , where k B denotes Boltzmann's constant. Thus the Lamb-shift contribution only encodes vacuum effects, while the prefactor Γ of the dissipator depends on Θ and yields a non-vanishing contribution for Θ = 0. The renormalisation in this model applied in [27] then removes the Lamb-shift contribution completely, as it depends on the unphysical cutoff frequency Ω, while the prefactor of the dissipator is not altered, as here the dependency on Ω vanishes after the Markov approximation.  \nNext, let us discuss to what extent it is possible to connect the renormalisation and its effects of the one-particle master equation presented above with the renormalisation applied in the quantum mechanical toy model. For this purpose first we discuss the two forms of the original full field theoretical master equation derived in [1]. The first form is given in that work in equation (4.52):  \n∂ ∂t ρ S ( t ) = -i [ H S + κ U, ρ S ( t )] -κ 2 ∫ t 0 ds ∑ r ∫ R 3 d 3 k { iD ( ⃗ k, t -s ) [ J r ( ⃗ k ) , { J r ( -⃗ k, s -t ) , ρ S ( t ) }] + D 1 ( ⃗ k, t -s ) [ J r ( ⃗ k ) , [ J r ( -⃗ k, s -t ) , ρ S ( t ) ]] } , (4.52)  \nwhere D and D 1 are two coefficients that arise from combinations of the environmental correlation functions similar to Λ and Γ in the quantum mechanical model and they read:  \nD ( ⃗ k, t -s ) := -sin(Ω k ( t -s )) Ω k (4.53)  \nD 1 ( ⃗ k, t -s ) := coth ( β Ω k 2 ) cos(Ω k ( t -s )) Ω k . (4.54)  \nThe operators J r ( ⃗ k, t ) and J r ( ⃗ k ) := J r ( ⃗ k, 0) were defined in [1] in equation (3.18) and contain a combination of two creation and/or annihilation operators of the scalar field along with their time evolution. As a first difference to the master equation in [27] it turns out that the term proportional to D in (4.52) cannot be written as a simple commutator, as it is the case with the Lamb-shift contribution in the quantum mechanical toy model. If the system's operator J were to commute with the system Hamiltonian, then this would be possible, and this is the case in the quantum mechanical toy model in [27]. Then this would imply that the Lamb-shift is independent of the temperature parameter Θ and therefore a pure vacuum effect.  \nFor the field-theoretical model, a similar form where one has a Lamb-shift contribution and a dissipator is the one given at the beginning of Appendix A.I.1. Here, the coefficients of the Lamb-shift term are S ab and the prefactors of the dissipator are R ab . From their definitions in (A.I.5) and (A.I.6) one can see that in general they have a different form as Λ and Γ in (4.50) and (4.51) above and the Lamb-shift includes vacuum as well as thermal contributions. A similar result is obtained in [2], where quantum electrodynamics is discussed from the point of view of open quantum systems with the standard interaction Hamiltonian of QED. There, the resulting Lamb-shift Hamiltonian is therefore split into a vacuum part, denoted as Lamb-shift and a thermal part, denoted as Stark-shift.  \nAs discussed above, we would expect that if the system's operator in the interaction Hamiltonian commutes with the system Hamiltonian, that we can then recover the form of the quantum mechanical model in [27]. Indeed, if we had [ J, H S ] = 0, then the phases e ± iω a ( ⃗ k,⃗p ) t and e ± iω b ( ⃗ k, ⃗ l ) t coming from the time evolution of the J operators would vanish in the definitions of S ab and R ab in (A.I.5) and (A.I.6). This allows the remaining terms to be combined into a form similar to Λ and Γ. In particular the thermal contribution of the Lamb-shift would vanish, as it is the case in the quantum mechanical toy model.  \nLet us now compare the renormalisations of the two models: in the quantum mechanical model in [27], the effect of the renormalisation is to remove the Lamb-shift Hamiltonian which only consisted of a vacuum part. The renormalisation applied in the present work removes the vacuum parts in the Lamb-shift Hamiltonian and the dissipator. The thermal part of the Lamb-shift Hamiltonian however remains. From the discussion of the open QED model from [2], one would expect a similar result in a quantum mechanical model where a thermal contribution in the Lamb-shift is present. The dissipator of the quantum mechanical model [27] is left unmodified by the renormalisation, in particular the vacuum contribution is present there. This is in contrast to the procedure here, where the renormalisation removes all vacuum terms, also the ones from the dissipator. In the quantum mechanical toy model, these contributions are however removed at a later stage when the Markov approximation is applied and hence also not present in the final Lindblad equation.  \nThis concludes the discussion on the renormalisation of the one-particle master equation. In the next section, we discuss how one can apply specific physical approximations to draw physical implications from the renormalised one-article master equation.", 'V. APPLICATION OF THE MARKOV AND ROTATING WAVE APPROXIMATIONS TO TRANSFORM THE TCL MASTER EQUATION INTO LINDBLAD FORM': 'The renormalised TCL one-particle master equation (4.47) describes the evolution of a single scalar particle in an environment filled with thermal gravitational waves. Since this master equation is not in Lindblad form, we cannot directly conclude that it is completely positive and provides physically meaningful implications based on positive probabilities for all chosen time intervals. For such models, one usually has to investigate case by case whether further assumptions such as the Markov and rotating wave approximation are justified that are usually used to obtain a master equation in Lindblad form. It is often possible to understand from the involved time scales in the system and environment of the open quantum model in which scenarios these approximations are a good choice, see for example [2] for a discussion in quantum optics. For models with finitely many degrees of freedom, there are also results that suggest time scales which allow to judge when the Markov approximation can be applied that are completely determined by the properties of the environment, such as its spectral density as well as the coupling constant, which encodes the strength of the coupling to the system in the interaction Hamiltonian [54].  \nThe derivation of master equations in the context of field-theoretical models with gravity as an environment is less well explored in the literature in comparison and has been presented in the context of gravitationally induced decoherence recently for instance in [1, 8, 9, 11, 12, 16]. While the works in [1, 8, 11, 16] focus on the derivation of a TCL master equation, in [9, 12] a Lindblad equation is used, for which further approximations are employed, among these the Markov approximation and the rotating wave approximation.  \nCompared to the above-mentioned open quantum mechanical models for gravitationally induced decoherence, a detailed analysis of the applicability of such approximations is much more challenging and beyond the scope of this article. An important difference to the present work is that in [9, 12] the approximations are applied on the non-renormalised one-particle master equation. Given the results of the last section, we can instead perform the Markov and rotating wave approximations for the renormalised one-particle master equation and investigate whether applying these approximations before or after renormalisation leads to differences in the final one-particle master equation, considering both the extended and non-extended one-particle projection.  \nIn this section we consider both approximations separately, in subsection V A we discuss the Markov approximation and in subsection V B the rotating wave approximation. In addition, for the case of an ultra-relativistic limit, we also specify some conditions when the Markov approximation can be used for the model considered here.', 'A. Markov approximation': 'The Markov approximation consists in the assumption that the correlation functions of the environment are strongly peaked around the initial time and decay rapidly. If this is given, the integral ∫ t -t 0 0 dτ over these environmental correlation functions has the main contribution from around their peak. Thus the error obtained when shifting the initial time t 0 → -∞ and therefore the upper integration limit t -t 0 →∞ is negligibly small. As a consequence, the parameters involved in the dissipator of the final Lindblad equation will no longer depend on the temporal coordinate.  \nAs discussed above for a field theoretical model, even in the single particle sector, to develop generic criteria for which the Markov approximation can be applied that can easily be checked for a given model, is difficult. For instance the methods developed in [54] strongly rely on the fact that the model is formulated in a quantum mechanical context. Motivated by the physical applications  \nin section VI D to ultra-relativistic particles, in particular neutrinos and their oscillations, as a first step, we investigate this special case more in detail in this context and present a condition under which the Markov approximation can be applied for the model under consideration. The details are discussed in appendix A.IV.1 by analysing the individual parts of the master equation and where we show that a suitable condition for the applicability of the Markov approximation in the ultra-relativistic case is the requirement that  \nu, v ≫ 1 cβ , (5.1)  \nwhere u := | ⃗u | , v := | ⃗v | and c denotes the speed of light. The reason why β and thus the temperature parameter Θ are involved here is because we use a Gibbs state to trace out the environmental and thus gravitational degrees of freedom. The specific time scales then yield that for the application in VI D the correction terms to the Markov approximation are negligible. If another than the ultra-relativistic case is considered, the above condition could be violated, and therefore a more comprehensive analysis is needed to understand when and under what conditions the Markov approximation can be applied, which we envisage for future work.  \nBefore applying the Markov approximation to the renormalised one-particle master equation, we briefly discuss the main steps that are involved: the first step is to perform the τ -integration using the identity  \n∫ ∞ 0 dτ e -iωτ = ∫ ∞ -∞ dτ Θ( τ ) e -iωτ = 2 π ∫ ∞ -∞ dx Θ( x ) e -2 πiωx = πδ ( ω ) -PV ( i ω ) , (5.2)  \nwhere PV denotes the Cauchy principal value, and secondly the evaluation of the ⃗ k -integration which simplifies due to the first step. With this, the affected terms in the master equation can then be split into two classes: one class that consists of contributions involving the delta distribution that will yield a real contribution to the master equation and hence lead to decoherence. Another class that contains terms including the principal value, which result in an imaginary contribution that affects the unitary evolution. The detailed computation for the contributions leading to decoherence can be found in appendix A.IV.2. Those contributions that involve the principal value  \nare evaluated in appendix A.IV.3. The final result, given in (A.IV.43), takes the form:  \n∂ ∂t ρ ( ⃗u, ⃗v, t ) = -iρ ( ⃗u, ⃗v, t ) ( ω u -ω v ) + κ 16 πβ { -10 3 ( u 2 + v 2 ) + 2( ω 2 u + ω 2 v ) -2 ω u u m 4 arctanh ( u ω u ) -2 ω v v m 4 arctanh ( v ω v ) -ω u ω v ( v 2 -3 ( ⃗u · ⃗v ) 2 u 2 )[ 2 3 -m 2 u 2 + m 4 ω u u 3 arctanh ( u ω u )] -ω v ω u ( u 2 -3 ( ⃗u · ⃗v ) 2 v 2 )[ 2 3 -m 2 v 2 + m 4 ω v v 3 arctanh ( v ω v )] } ρ ( ⃗u, ⃗v, t ) + iκ 2(2 π ) 3 1 √ ω u ω v ∫ d 3 k P ijln ( ⃗ k ) u i u j v l v n √ ω u + k ω v + k N ( k ) Ω k ρ ( ⃗u + ⃗ k,⃗v + ⃗ k, t ) { PV ( 1 Ω k + ω u -ω u + k ) -PV ( 1 Ω k -ω u + ω u + k ) -PV ( 1 Ω k + ω v -ω v + k ) +PV ( 1 Ω k -ω v + ω v + k ) } -(1 -δ P ) iκ 2(2 π ) 2 lim ϵ → 0 [ u 4 ω u ∫ ∞ ϵ dk ∫ π 0 dθ sin 5 ( θ ) kN ( k ) 1 -ω u ω u -k k 2 -( ω u -k -ω u ) 2 -v 4 ω v ∫ ∞ ϵ dk ∫ π 0 dθ sin 5 ( θ ) kN ( k ) 1 -ω v ω v -k k 2 -( ω v -k -ω v ) 2 ] ρ ( ⃗u, ⃗v, t ) .  \n(5.3)  \nThe contributions in lines two to four arose from the delta distributions and are real, so they cause decoherence in the evolution of the scalar particle. The remaining terms are imaginary and therefore contribute to the unitary evolution of the density matrix. When working with the extended one-particle projection, then the expressions in the last two lines vanish. The real part in lines two to four remains unaffected by the rotating wave approximation that will be applied in the next subsection, hence it already possesses its final form.  \nIn the existing literature, the master equation in the one-particle picture is usually directly specified or derived for the non-relativistic (see e.g. [9] for scalar particles) or the ultra-relativistic case (see e.g. [12] for photons). In these cases, the dissipator has a simpler form and there are no arctanhterms present as it is the case here in the general master equation, where neither the non- nor the ultra-relativistic limits have been applied yet.  \nTo further compare with the existing literature, in section VI D we will consider the non- and ultra-relativistic limit of this master equation above and show that the arctanh does not appear in either limit. Thus, it indeed leads to the results obtained in the literature.  \nAnother difference to similar work in [9, 12] is that the master equation in the present work has already been renormalised, i.e. all vacuum contributions that are independent of the temperature parameter Θ have been removed in these terms, which arise due to the gravitational influence in (5.3), which in particular contains all vacuum fluctuations of the gravitational field. This can be seen by setting the temperature parameter equal to zero, as then all terms including the gravitational influence vanish.', 'B. Rotating wave approximation': 'After having applied the Markov approximation, the rotating wave approximation is usually a next step in order to cast the master equation into a completely positive Lindblad form. The physical idea behind the approximation is to take into account that detectors only have a finite resolution and cannot resolve arbitrarily fast oscillations, but only measure a coarse-grained result. In the literature, there exist different ways to apply the rotating wave approximation. One possibility is, following the nomenclature in [24], the pre-trace RWA, where the approximation is applied at the level of the interaction Hamiltonian by dropping counter-rotating terms. This is often employed e.g. in quantum optics and leads to the Jaynes-Cummings model, see [55, 56], which is nowadays extensively studied for instance in quantum technology, see [57]. This pre-trace RWA, which is also applicable in closed quantum systems, has been studied from several angles yielding different results in the last years among other things on its higher order corrections and a renormalisation of the resulting series (see [26]) as well as also on the bounds of its applicability (see [25]). From the analysis in [24] it follows that in open quantum systems the second version of the RWA, the post-trace rotating wave approximation which is applied at the level of the master equation after tracing out the environment, yields dynamics which are expected to be closer to the true system dynamics. This analysis in [24] is carried out for quantum mechanical models and we expect that more work is required to extend it to the full field theoretical case. Nevertheless, we take this discussion as a motivation to apply in this work the post-trace RWA, which was also employed in similar analyses, for instance in [12]. This post-trace RWA is implemented by considering the master equation in the interaction picture and then dropping terms that rotate very fast compared to the other ones. The detailed implementation and computation can be found in appendix A.V, here we only state the result:  \n∂ ∂t ρ ( ⃗u, ⃗v, t ) = -iρ ( ⃗u, ⃗v, t ) ( ω u -ω v ) + κ 16 πβ { -10 3 ( u 2 + v 2 ) + 2( ω 2 u + ω 2 v ) -2 ω u u m 4 arctanh ( u ω u ) -2 ω v v m 4 arctanh ( v ω v ) -ω u ω v ( v 2 -3 ( ⃗u · ⃗v ) 2 u 2 )[ 2 3 -m 2 u 2 + m 4 ω u u 3 arctanh ( u ω u )] -ω v ω u ( u 2 -3 ( ⃗u · ⃗v ) 2 v 2 )[ 2 3 -m 2 v 2 + m 4 ω v v 3 arctanh ( v ω v )] } ρ ( ⃗u, ⃗v, t ) -(1 -δ P ) iκ 2(2 π ) 2 lim ϵ → 0 [ u 4 ω u ∫ ∞ ϵ dk ∫ π 0 dθ sin 5 ( θ ) kN ( k ) 1 -ω u ω u -k k 2 -( ω u -k -ω u ) 2 -v 4 ω v ∫ ∞ ϵ dk ∫ π 0 dθ sin 5 ( θ ) kN ( k ) 1 -ω v ω v -k k 2 -( ω v -k -ω v ) 2 ] ρ ( ⃗u, ⃗v, t ) . (5.4)  \nIt can be seen when comparing this result to (5.3) that the effect of the rotating wave approximation is to remove the remaining part of the Lamb-shift in the extended projection. In the non-extended projection, there still survives one term of the Lamb-shift which corresponds to the last two lines in (5.4). Apart from that, the rotating wave approximation causes no further modifications on the master equation. This is due to the fact that all other terms that would be removed by the approximation were already dropped when performing the one-particle projection of the master equation. The general dissipator at the operator level can however be written in Lindblad form  \nafter the RWA, see (A.V.45) in the appendix:  \nD δ [ ρ S ] = κ ∑ r ∈{ + , -} ∫ R 3 d 3 k (2 π ) 2 δ ( k ) N ( k ) Ω k ( L r ( ⃗ k ) ρL r ( ⃗ k ) † -1 2 { ρ, L r ( ⃗ k ) † L r ( ⃗ k ) } ) (5.5)  \nwith Lindblad operators  \nL r ( ⃗ k ) := ∫ R 3 d 3 p 1 √ 1 -( k + p cos( θ p )) 2 ω 2 k + p J 2 r ( ⃗ k, ⃗p ) , (5.6)  \nwhere θ p denotes the angle between ⃗ k and ⃗ p and J 2 r ( ⃗ k, ⃗p ) = 2 j 1 r ( ⃗ k, ⃗p ) with the latter being defined in (A.I.8). With this, we have derived the final form of the renormalised one-particle master equation after Markov and rotating wave approximation. In the next section, we discuss some applications and investigate some features of the master equation at different intermediate stages before, during and after the applied approximations.', 'VI. APPLICATIONS OF THE ONE-PARTICLE MASTER EQUATION': 'In this section, we discuss some applications of the one-particle master equation derived in the previous sections. We start with analysing the evolution of the populations of the one-particle density matrix with a special focus on the interplay between the renormalisation and Markov and rotating wave approximations in subsection VI A and compare the results to [11] where the evolution of the populations of the non-renormalised TCL master equation is derived. Next we discuss the non-relativistic limit of the one-particle master equation in subsection VI B and compare the results to the ones in [9]. Furthermore, we investigate the ultra-relativistic limit in section VI C, compare it to [12], and discuss the relation to the quantum mechanical model for gravitationally induced decoherence in neutrino oscillations in [27] in section VI D. This further allows to connect to phenomenological models that investigate the influence of gravity on neutrino oscillations, like for instance in [7, 10, 17, 58].', 'A. Evolution of the populations of the one-particle master equation': 'We start by analysing the dynamics of the populations, that is the diagonal elements, in momentum representation predicted by the master equation at different stages in the derivation of the final Lindblad equation. We have chosen this application because it is an example that allows us to discuss and compare the implications that arise depending on the stage of the calculation at which the renormalisation procedure is performed.  \nTo investigate the evolution of the populations, we take the different versions of the master equation and compute it for ρ ( ⃗ k, t ) := ρ ( ⃗ k, ⃗ k, t ) before and after the renormalisation as well as after the Markov approximation. As the rotating wave approximation only affects the Lamb-shift Hamiltonian, the dynamics of the populations will not get modified after its application.', '1. Before renormalisation': "The dynamics of the populations in the one-particle master equation (3.9) before renormalisation and further approximations can be obtained by evaluating the master equation for ρ ( ⃗ k, t ) := ρ ( ⃗ k, ⃗ k, t ). In this case, we have no contribution from the unitary dynamics and in the dissipator all  \nimaginary parts will vanish 7 and one obtains a dissipator that is purely real. In this subsection we adapt the notation to the one used in [11] in order to better facilitate the comparison with their results:  \n˙ ρ ( ⃗ k, t ) = -κ ∫ d 3 k ' (2 π ) 3 P k ( ⃗ k ' -⃗ k ) Ω k ' -k ω k ' ω k · · {[ [ N ( k ' -k ) + 1] sin( χ ( t -t 0 )) χ + N ( k ' -k ) sin( χ ' ( t -t 0 )) χ ' + δ P [ N ( k ' -k ) + 1] sin( η ( t -t 0 )) η + δ P N ( k ' -k ) sin( η ' ( t -t 0 )) η ' ] ρ ( ⃗ k, t ) -[ [ N ( k ' -k ) + 1] sin( χ ' ( t -t 0 )) χ ' + N ( k ' -k ) sin( χ ( t -t 0 )) χ ] ρ ( ⃗ k ' , t ) } (6.1)  \nwith ˙ ρ ( ⃗ k, t ) = ∂ t ρ ( ⃗ k, t ) and χ := Ω k ' -k -ω k + ω k ' , χ ' := Ω k ' -k + ω k -ω k ' , η := Ω k ' -k + ω k + ω k ' and η ' := Ω k ' -k -ω k -ω k ' as well as P k ( ⃗ k ' -⃗ k ) = P ijln ( ⃗ k ' -⃗ k ) k i k j k l k n . Using that P ijln ( ⃗ k -⃗ k ' )( ⃗ k -⃗ k ' ) i = 0 as P ijln ( ⃗ k -⃗ k ' ) projects onto the symmetric transverse traceless part and therefore removes the longitudinal part ∝ ⃗ k -⃗ k ' , which can be seen from the definition in (3.10), we can use  \nP ijln ( ⃗ k -⃗ k ' ) k i = P ijln ( ⃗ k -⃗ k ' ) k ' i (6.2)  \nand hence rewrite  \nP k ( ⃗ k ' -⃗ k ) = P ijln ( ⃗ k ' -⃗ k ) k i k j k ' l k ' n . (6.3)  \nFrom equation (6.1) one can also once more see the implication of the chosen projection, i.e. whether δ P = 0 or δ P = 1, on the probability conservation, which was discussed below equation (3.8). When working with the non-extended projection δ P = 0, then we have  \n∫ R 3 d 3 k ˙ ρ ( ⃗ k, t ) = 0 (6.4)  \ndue to symmetry and thus probability in the scalar particle's subsystem is conserved. If working with the extended one-particle projection δ P = 1 instead, it can be seen in equation (6.1) that the terms containing η and η ' lack a symmetric counterpart to be cancelled and hence in that case probability conservation is not given any more, as it was also discussed below equation (3.8).  \nIn [11] the dynamics of the population for a master equation of a photon coupled to linearised gravity are discussed. We obtain an agreement with their result if we specialise to a massless scalar particle and choose as the initial time t 0 = 0. In addition we need to consider the non-extended one particle projection (i.e. δ P = 0), in order to adapt to their chosen normal ordering as well as choose the temperature parameter Θ to be zero. The latter corresponds to a vacuum state of the gravitational waves environment. Inserting these assumptions in the evolution of the populations this equation becomes  \n˙ ρ ( ⃗ k, t ) = -κ ∫ d 3 k ' (2 π ) 3 P k ( ⃗ k ' -⃗ k ) Ω k ' -k ω k ' ω k [ sin( χt ) χ ρ ( ⃗ k, t ) -sin( χ ' t ) χ ' ρ ( ⃗ k ' , t ) ] , (6.5)  \nwhich has a very similar form as the one in [11] for a photon. The only difference arises due to the fact that for the photons in [11] the polarisation vectors couple to the symmetric transverse traceless projector while here for the scalar particles, as they do not carry any polarisation, This role is taken over by the momentum, which is the only direction-dependent quantity that scalar particles possess.", '2. After renormalisation': "As discussed in section IV E, the effect of the renormalisation was that the vacuum part in the one-particle master equations, these are the contributions not involving N ( ⃗ k ), vanishes. At the practical level this can be implemented by replacing everywhere N ( k ' -k ) + 1 → N ( k ' -k ). Then the dynamics of the populations becomes  \n˙ ρ ( ⃗ k, t ) = -κ ∫ d 3 k ' (2 π ) 3 P k ( ⃗ k ' -⃗ k ) Ω k ' -k ω k ' ω k · N ( k ' -k ) · {[ sin( χ ( t -t 0 )) χ + sin( χ ' ( t -t 0 )) χ ' + δ P sin( η ( t -t 0 )) η + δ P sin( η ' ( t -t 0 )) η ' ] ρ ( ⃗ k, t ) -[ sin( χ ' ( t -t 0 )) χ ' + sin( χ ( t -t 0 )) χ ] ρ ( ⃗ k ' , t ) } . (6.6)  \nWe realise that now all terms depend on N ( k ' -k ). As a consequence, the entire evolution of the populations is trivial, that is ˙ ρ ( ⃗ k, t ) vanishes, if we consider the specific case of a vanishing temperature parameter Θ = 0 yielding directly N ( k ' -k ) = 0 for all k, k ' . The comparison to the non-renormalised master equation shows that the physical properties of the two one-particle master equations are quite different as far as the dynamics of the populations is concerned. For this reason the discussions and physical implications drawn in [11] based on the dynamics of the populations in the non-renormalised equation (6.5) are problematic, as the evolution of the diagonal terms vanishes after renormalisation in the zero temperature limit.", '3. After the Markov approximation': 'Due to the fact that for the diagonal elements the coefficients always enter in the form C + C ∗ , only real terms in the one-particle master equation after the second Markov approximation in (5.3) remain:  \n∂ ∂t ρ ( ⃗ k, t ) = ˙ ρ ( ⃗ k, t ) = κ 16 πβ { -20 3 k 2 +4 ω 2 k -4 ω k k m 4 arctanh ( k ω k ) +4 k 2 [ 2 3 -m 2 k 2 + m 4 ω k k 3 arctanh ( k ω k )] } ρ ( ⃗ k, t ) =0 . (6.7)  \nThis means that the Markov approximation removes the dynamics of the populations also in the case of non-vanishing temperature and independently of the extended projection δ P . It therefore  \nalso restores probability conservation, as it removes all terms from the extended projection in the dissipator. The rotating wave approximation only affects the imaginary parts of the master equation, thus it does not change the evolution of the populations any more.  \nThis result that the dynamics of the populations vanishes is also obtained in [9] for a non-relativistic one-particle master equation that was renormalised after the application of Markov and rotating wave approximation, and for the one of a photon after renormalisation and application of the same two approximations in [12].', 'B. Non-relativistic limit': 'In the following, we apply the renormalised one-particle master equation after Markov and rotating wave approximation (5.4) to non-relativistic particles in order to compare the decoherence with the one derived in [9].  \nIn the non-relativistic limit we have u 2 m 2 ≪ 1 and v 2 m 2 ≪ 1 and due to this the one-particle master equation simplifies. In this case we can expand the arctanh as  \n-2 m 2 m u 1 √ 1 + u 2 m 2 arctanh   u m 1 √ 1 + u 2 m 2   = -2 m 2 m u [ u m -2 3 u 3 m 3 + 8 15 u 5 m 5 + O ( u 6 m 6 )] . (6.8)  \nGiven this we find for the contribution from lines two to four in (5.4) which is the part leading to decoherence the following expression:  \nκ 16 πβ { -16 15 u 4 + v 4 m 2 -16 15 u 2 v 2 m 2 (1 -3 cos 2 ( γ ) } = -κ 5 πβm 2 [ 1 3 ( u 4 + v 4 ) + u 2 v 2 ( 1 3 -cos 2 ( γ ) )] , (6.9)  \nwhere γ is defined as the angle between ⃗u and ⃗v , i.e. ⃗u · ⃗v = uv cos( γ ). We work with the extended projection δ P = 1 here, as a consequence there is no Lamb-shift contribution left and the master equation becomes  \n∂ ∂t ρ ( ⃗u, ⃗v, t ) = -iρ ( ⃗u, ⃗v, t )( ω u -ω v ) -κ 5 πβm 2 [ 1 3 ( ⃗u 4 + ⃗v 4 ) + ⃗u 2 ⃗v 2 ( 1 3 -cos 2 ( γ ) )] ρ ( ⃗u, ⃗v, t ) = -iρ ( ⃗u, ⃗v, t ) ( u 2 2 m -v 2 2 m ) -κ 5 πβm 2 [ 1 3 ( ⃗u 4 + ⃗v 4 ) + ⃗u 2 ⃗v 2 ( 1 3 -cos 2 ( γ ) )] ρ ( ⃗u, ⃗v, t ) , (6.10)  \nwhere we used in the last step ω u = √ ⃗u 2 + m 2 = m √ u 2 m 2 +1 ≈ m + u 2 2 m and likewise for ω v . The master equation (60) in [9], where a Lindblad equation is used after Markov and rotating wave approximation to also describe a scalar field coupled to a linearised gravitational field, reads in momentum representation 8  \n∂ ∂t ρ ( ⃗u, ⃗v, t ) = -iρ ( ⃗u, ⃗v, t ) ( u 2 2 m R -v 2 2 m R ) -2 κ 3 βm 2 R [ 1 3 ( ⃗u 4 + ⃗v 4 ) -1 3 ⃗u 2 ⃗v 2 ( 1 + cos 2 ( γ ) ) ] ρ ( ⃗u, ⃗v, t ) . (6.11)  \nWhile in that work, they use the same underlying physical system, one of the differences is that there a gauge fixing is used while in this work the elementary physical variables were identified in [1]  \nby choosing geometrical clocks with respect to which suitable Dirac observables were constructed. Additionally, the Hamiltonian used in [9] is not completely normal ordered, while here we worked with a completely normal ordered one. A more detailed discussion of these two points can be found in [1]. Furthermore, the renormalisation is carried out in a different manner: in this work it is done before the Markov and rotating wave approximation are applied. In contrast in [9] the final master equation is renormalised after the application of these two approximations and after going into the non-relativistic limit. Their renormalisation procedure involves the introduction of a cutoff Λ ≪ m which is later absorbed in a redefinition of the renormalised mass m → m R , while here we found in equation (4.38) that only the wave function needs to be renormalised, see section IV D. In [9] compared to our result here, there are some additional unitary terms left due to using the non-extended one-particle projection. These contributions are proportional to the UV-cutoff Λ and to u 4 m 2 R , which is why they are dropped in [9] from the final master equation in the non-relativistic limit, even though in the limit Λ → ∞ they would diverge. As our results demonstrate, using the extended projection and a renormalisation before the application of the approximations hence removes the necessity to drop diverging terms by hand.  \nAdditional differences between (6.10) and (6.11) are the prefactor in front of the dissipator and the structure inside the square brackets. In these two points the results derived here do not agree with the results in [9]. Particularly regarding the last point, our result however agrees with a similar derivation for photons in [12] where more intermediate steps are provided and where the final structure in the square brackets is the same as in (6.10).', 'C. Ultra-relativistic limit': "In this subsection we apply the one-particle master equation to ultra-relativistic particles. Possible applications are one-particle master equations for photons as discussed in [12] as well as gravitationally induced decoherence in neutrino oscillations as for instance discussed in [27], where a quantum mechanical toy model was used.  \nIn the ultra-relativistic limit we have m 2 u 2 ≪ 1 as well as m 2 v 2 ≪ 1. Taking this into account, we neglect all terms of order O ( m 2 u 2 ) and O ( m 2 v 2 ) respectively and higher order contributions. Note, that this also includes terms involving arctanh function because  \nm 2 u 2 1 √ 1 + m 2 u 2 arctanh   1 √ 1 + m 2 u 2   = O ( m 2 u 2 ) .  \nThis leads to the the following simplification for the decoherence term in (5.4):  \nκ 16 πβ { -4 3 ( u 2 + v 2 ) -4 3 [ uv -3( ⃗u · ⃗v ) 2 ] } = -κ 4 πβ [ 1 3 ( u 2 + v 2 ) + uv ( 1 3 -cos 2 ( γ ) )] . (6.12)  \nThe remaining computation of the imaginary part in the dissipator can be found in (A.IV.40) in the appendix. Combining all results, the renormalised one-particle master equation in the ultra-  \nrelativistic limit can be written in the form  \n∂ ∂t ρ ( ⃗u, ⃗v, t ) = -iρ ( ⃗u, ⃗v, t ) ( u -v ) -κ 4 πβ [ 1 3 ( u 2 + v 2 ) + uv ( 1 3 -cos 2 ( γ ) )] ρ ( ⃗u, ⃗v, t ) -(1 -δ P ) iκu 4 ρ ( ⃗u, ⃗v, t ) 105(2 π ) 2 ω u { 4 [ π 4 15 β 4 u 4 -7 π 2 6 β 2 u 2 -6 ln(1 -e -βu ) βu +4 Li 2 ( e -βu ) β 2 u 2 -6 Li 3 ( e -βu ) β 3 u 3 -6 Li 4 ( e -βu ) β 4 u 4 ] +35 -35 ln ( e βu -1 ) βu -14 u ∫ ∞ u dk N ( k ) k 2 +3 u 3 ∫ ∞ u dk N ( k ) k 4 } +(1 -δ P ) iκv 4 ρ ( ⃗u, ⃗v, t ) 105(2 π ) 2 ω v { 4 [ π 4 15 β 4 v 4 -7 π 2 6 β 2 v 2 -6 ln(1 -e -βv ) βv +4 Li 2 ( e -βv ) β 2 v 2 -6 Li 3 ( e -βv ) β 3 v 3 -6 Li 4 ( e -βv ) β 4 v 4 ] +35 -35 ln ( e βv -1 ) βv -14 v ∫ ∞ v dk N ( k ) k 2 +3 v 3 ∫ ∞ v dk N ( k ) k 4 } , (6.13)  \nwhere γ denotes the angle between ⃗u and ⃗v and Li s ( x ) denotes the poly-logarithm function defined in (A.IV.41). In the extended projection, i.e. for δ P = 1, this becomes  \n∂ ∂t ρ ( ⃗u, ⃗v, t ) = -iρ ( ⃗u, ⃗v, t ) ( u -v ) -κ 4 πβ [ 1 3 ( u 2 + v 2 ) + uv ( 1 3 -cos 2 ( γ ) )] ρ ( ⃗u, ⃗v, t ) , (6.14)  \nwhich can be rewritten in terms of an operator equation as  \n∂ ∂t ˆ ρ ( t ) = -i [ ˆ H, ˆ ρ ( t )] + κ 8 πβ ( δ il δ jm -1 3 δ ij δ lm )[ ˆ p i ˆ p j ˆ p 0 , [ ˆ ρ ( t ) , ˆ p l ˆ p m ˆ p 0 ]] , (6.15)  \nwith ˆ p 0 := √ ˆ p n ˆ p n + ξ m 2 1 . In this definition, ξ m 2 is a small regulator that removes the eigenvalue zero from the spectrum of ˆ p 0 , as in that case the operator would not be invertible. For massive particles, this regulator corresponds to the mass squared m 2 which is still present in the ultrarelativistic limit, even though very small compared to the other summand. This is, up to a factor of 2, the same result for decoherence as derived in [12] for gravitationally induced decoherence of photons. This difference of a factor of 2 is already present when comparing the field-theoretical models of [12] and [1]. Note that in [12] the derivation and in particular the application of the approximations is performed without a prior renormalisation of the one-particle master equation, which is done in the end to get rid of the diverging Lamb-shift term. As expected from the analysis in this work, they find a logarithmic divergence in the end. The derivation of the master equation in [12] is very similar to the one in [9], hence we refer for a detailed comparison to the discussion in subsection VI B. The renormalisation in [12] is done after performing the approximations and the ultra-relativistic limit such that the detailed procedure depends on the cutoff frequency Λ and its relation to the photon frequency ω u (in our case the scalar particle's frequency). In the end in [12] the electric and magnetic fields as well as the coupling constant are renormalised.", 'D. Application to neutrino oscillations': "Finally, we want to discuss the relation of the results obtained in this work with the one presented in [27], where gravitationally induced decoherence in neutrino oscillations is investigated based on a quantum mechanical toy model [28] with neutrinos as system of interest and a collection of Harmonic oscillators to model the thermal gravitational waves environment. For this purpose we consider the decoherence of neutrinos predicted by the ultra-relativistic one-particle master equation derived in this work.  \nThe connection to the work in [27] is of interest to us in several respects: first due to the quantum mechanical nature of the the toy model used in [27], the coupling parameter encoding the strength for the coupling between system and environment cannot be determined from first principles. Instead a free parameter was introduced, that after introduction of a spectral density was denoted by 9 η , similar to what was done also in [28]. Second, the quantum mechanical model requires the choice of a spectral density to derive the final master equation together with an appropriate cut-off function that regulates the integral over the frequency domain. In [27], four different commonly used cut-off functions were considered and it was shown that the final result of the master equations cannot distinguish between the different choices. For the spectral density, the usual linear dependence on the frequency, which is widely used in the context of quantum optics, was considered. As the one-particle master equation in this work is derived from an underlying field theory model presented in [1] the situation is different here. As the matter, which is a scalar field in the present work, couples to linearised gravity, the coupling constant in the interaction Hamiltonian is naturally build into the model and given by κ = 8 πG N , where G N is Newton's constant. Furthermore, due to the field-theoretical character of the model, it is not necessary to introduce a spectral density by hand, since the interaction Hamiltonian contains an integral over all modes from the beginning. As a third aspect we want to compare the application of the Markov approximation in the quantum mechanical toy model and in the one-particle master equation derived in this work.  \nEven though the equation derived here is, taken strictly, only applicable to scalar particles, we still apply it to the case of neutrinos in order to discuss the relation with the results in [27]. That this can be done in this context is due to the reason that the quantum mechanical toy model investigated in [27] treats the neutrinos as plane waves and thus does not take the full spinorial nature of neutrinos into account.  \nWe assume that the two momenta ⃗u and ⃗v are approximately parallel to each other in order to have intersection probability and to be able to measure them in a neutrino detector. With this, the one-particle master equation in the ultra-relativistic limit becomes  \n∂ ∂t ρ ( ⃗u, ⃗v, t ) = -iρ ( ⃗u, ⃗v, t ) ( u -v ) -κ 12 πβ ( u -v ) 2 ρ ( ⃗u, ⃗v, t ) . (6.16)  \nAs for ultra-relativistic particles we have ˆ H ˆ ρ ≈ u ˆ ρ , we indeed obtain the same form for the master equation as in [27] in the effective mass basis 10 :  \n∂ ∂t ˆ ρ ( t ) = -i [ ˆ H, ˆ ρ ( t )] + κ 6 πβ ( ˆ H ˆ ρ ( t ) ˆ H -1 2 { ˆ H 2 , ˆ ρ ( t ) } ) . (6.17)  \nThe master equation above can be solved in the energy eigenbasis, where we denote the energy eigenvalues by E i . With respect to this basis we denote the individual elements of ˆ ρ by ρ ij whose  \nsolutions read  \nρ ij ( t ) = ρ ij (0) e -i ( E i -E j ) t -κ 12 πβ ( E i -E j ) 2 t . (6.18)  \nThis result agrees, up to the different prefactor of 2 in front of the decoherence term mentioned in subsection VI C, with the one obtained for the one-particle projection evaluated for motion in one dimension for a photon in [12].  \nWe will now discuss the comparison of the three aspects mentioned above. We start with the comparison of the coupling parameter. Such a comparison can be obtained by comparing the prefactors of the decoherence terms in the model from [27] and here. For the latter we have  \nκ 12 πβ = 2 G N k B Θ 3 c 4 → 2 G N k B Θ 3 ℏ 2 c 5 , (6.19)  \nwhere we restored the correct units in the last step. Comparing with the decoherence rate in [27], which is 4 η 2 k B Θ ℏ 3 , and introducing the Planck length ℓ P , we find that  \nη 2 = ℏ G N 6 c 5 = ℓ 2 P 6 c 2 ≈ 5 · 10 -88 s 2 . (6.20)  \nAn estimate for the value of η 2 inspired from field theory was already discussed in the appendix in [27]. The difference in the orders of magnitude compared to the analysis in the appendix in [27] arises due to the numerical prefactors that could not be determined precisely by the analogy analysis and estimate in [27]. Similar results of the coupling strength can be found in [9, 12]. As the one-particle master equation considers the case of a scalar field with a thermal gravitational background, more work is needed in order to develop more sophisticated models for neutrinos or fermions in general to derive a similar master equation for a fermionic system under consideration and for more general environments. In addition the model is based on linearising gravity around a flat Minkowski background, whereas it would be interesting to also consider decoherence models for longer propagation distances and consider master equations based on a model on a cosmological background as the presence of a scale factor could modify the decoherence effect, as analysed for instance in [59-61].  \nCompared to the quantum mechanical model [27], as mentioned above, here it was not necessary to introduce a spectral density, which is a kind of continuum limit for the frequencies of the oscillators in the environment that is typically used in similar quantum mechanical models to avoid Poincare recurrences, see for instance [2]. Furthermore, the cut-off function that needed to be used in the quantum mechanical model to regularise divergent integrals is not required for the one-particle master equation here. Instead the divergent contributions in the one-particle master equation could be linked to Feynman diagrams of a corresponding effective field theory giving a clearer physical interpretation than in the quantum mechanical toy model. With that given, the divergent contributions were treated using a standard renormalisation procedure known from quantum field theory that would be applied also in other situations in quantum field theory where such kind of diagrams play a role.  \nNot entirely unrelated to the latter paragraph is the discussion of the application of the Markov approximation in the model considered here and the quantum mechanical toy model in [27]. In the latter, due to its simplicity compared to the one-particle master equation considered here, it was explicitly shown that the environmental correlation functions are strongly peaked around the initial time and decay rapidly after the peak. Such environmental correlation functions depend on both  \nthe chosen spectral density and the chosen cut-off function. In this work, however, none of these choices are made, but the corresponding quantities are determined and set from the beginning when formulating the model. In section V A a condition was discussed under which the Markov approximation can be applied for the ultra-relativistic limit. Considering here the application to neutrinos we can discuss whether this condition is satisfied in this application and how it relates to the application of the Markov approximation in the quantum mechanical toy model in [27]. The condition applied in this work, equation (5.1), states that the Markov approximation is justified if 11  \nu, v ≫ k B Θ c . (6.21)  \nIn the case of ultra-relativistic neutrinos where we neglect their masses, this is equivalent to  \nE u , E v ≫ k B Θ , (6.22)  \nwhere E u and E v denote the neutrino energies. Typical neutrino energies investigated in [27] start at energies of 1 GeV ≈ 1 . 6 · 10 -10 J . Given the Boltzmann constant k B ≈ 1 . 4 · 10 -23 J K and the temperature parameter Θ of the thermal gravitational waves of around 1 K used in [27], is the condition for the applicability of the Markov approximation used in this work  \n1 . 6 · 10 -10 J ≫ 1 . 4 · 10 -23 J . (6.23)  \nBoth sides of the inequality still differ by more than ten orders of magnitude, so the approximation can also be used for neutrinos with lower energies or for higher values of the temperature parameters. Thus, we can conclude given the proof presented in this work, the application of the Markov approximation to the physical scenario used in [27] is not only justified at the level of the quantum mechanical model as shown in [27] but also if one derives that model from the one-particle master equation of the QFT model.  \nFinally, compared to the quantum mechanical toy model in [27], the more general one-particle master equation derived here also allows to consider the generalisation to decoherence models with wave packets, whereas in [27] plane waves were considered. This is due to the general form of the one particle density matrix defined in (3.1) as  \nρ 1 ( t ) = ∫ R 3 d 3 u ∫ R 3 d 3 v ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v . (6.24)  \nWhen choosing suitable initial conditions for ρ ( ⃗u, ⃗v, 0), one can model different descriptions for neutrinos like wave packets or plane waves, where the latter just correspond to delta distributions in this context.", 'VII. CONCLUSIONS': 'In this paper, we investigate the one-particle sector of the field-theoretic model in [1] for gravitationally induced decoherence for a scalar field coupled to linearised gravity. One of the main foci of this work is the renormalisation of the one-particle master equation, which, in contrast to the existing literature [9, 12], is performed before applying approximations such as the Markov or rotating wave approximation. As our results show this strategy provides a physical interpretation  \nof the UV-divergent contributions in the one-particle master equation as being the vacuum part of the self-energy of the scalar field. To demonstrate this explicitly several steps are necessary: first the one-particle projection, where we considered two different kinds in this work, the non-extended and extended one. The latter also includes those processes in which in an intermediate steps two particles are created and annihilated afterwards or when the initial particle is left invariant and a vacuum bubble is created. Compared to the non-extended projection the last situation requires and additional renormalisation. While the non-extended projection yields a one-particle master equation with probability conservation, the extended projection does not have that property before any renormalisation, approximations or limits are taken into account. Then for both one-particle projections we identify the UV divergent contributions in the one-particle master equation. It turns out that all thermal parts of the one-particle master equation are UV finite and only the vacuum contributions in those terms which do not include other states than the one under consideration, are the divergent ones. Since in open quantum mechanical models the Lamb-shift Hamiltonians often require renormalisation and their physical interpretation is given, we wanted to address the question of the physical interpretation of the UV-divergent contributions in the one-particle master equation determined from an underlying field-theoretic model.  \nFor this purpose, we used the methods introduced in [21] and applied them to an open quantum model with thermal gravitational waves as environment, instead of an environment consisting of a scalar field as in [21]. These methods allow us to identify contributions in the one-particle master equation with certain Feynman diagrams of the effective field theory for the scalar field. Since the model in [1] is based on a canonical quantisation of the master equation in a first step we consider the interaction part of the effective quantum field theory model in its canonical from and introduce the corresponding non-covariant Feynman rules along the lines of [22, 23], where similar methods are used in the framework of QED. Equipped with them, we consider the (non-covariant) self-energy diagram, which can be decomposed into a thermal and a vacuum part, and show that the latter can be identified with the UV-divergent contributions in the one-particle master equation. In order to apply a standard on-shell renormalisation procedure for the self-energy of the scalar field, we relate the non-covariant Feynmann rules to the covariant rules. Interestingly, a sum of two non-covariant Feynman diagrams, expressed by non-covariant propagators, can be combined into one diagram containing the corresponding covariant propagator. In the work in [21] the Feynman rules were directly available in covariant form and therefore the introduction of non-covariant Feynman rules was not necessary there. We then present the final renormalised one-particle master equation.  \nAs our results show, the renormalisation leads to a renormalisation of the wave function in the model considered in the present work and the effect of the renormalisation is that all vacuum contributions in the Lamb-shift terms as well as in the dissipator are no longer present in the renormalised one-particle master equation. Comparing our results to the one in [9], where the renormalisation is performed after a Markov and rotating wave approximation have been applied and the non-relativistic limit has been considered, there the renormalisation procedure consists of introducing a cutoff function that is later absorbed in a redefinition of the renormalised mass. In [12] the effective model for the photon is also renormalised after applying the Markov and rotational wave approximation and considering the ultra-relativistic limit, and there the electric and magnetic field as well as the coupling constant are renormalised. Our results also allow a comparison to the renormalisaiton performed in the quantum mechanical model in [27]. While there the Lamb-shift Hamiltonian, that consists in that model only of vacuum contributions, needs to be renormalised and is absent after renormalisation, in the renormalised one-particle master equation in this work, the thermal part of the Lamb-shift Hamiltonian remains. This agrees with the situation one has when QED is treated as an open quantum model [2], where only the vacuum  \ncontribution in the Lamb-shift terms is denoted as Lamb-shift and the thermal contribution as Stark shift. From the discussion of the open QED model from [2], one would expect a similar result in a quantum mechanical model where a thermal contribution in the Lamb-shift is present. We also see differences in how renormalisation affects the dissipator. While in the model presented here the renormalisation removes the vacuum contributions from the dissipator, in [27] the dissipator is not changed when the renormalisation is applied and therefore the vacuum contributions remain. However, these are removed as soon as the Markov approximation is applied and therefore do not contribute to the final Lindblad equation in [27].  \nGiven the renormalised one-particle master equation we then discuss the application of the Markov and rotation wave approximation. While a general analysis of the applicability of the Markov approximation is beyond the scope of this article, we can find a condition for the ultra-relativistic case under which the Markov approximation can be applied. For the applications to neutrino oscillations considered in [27], this condition is very mild and does not lead to severe restrictions. This fits with the fact that it could be explicitly shown for the quantum mechanical model in [27] that the Markov approximation can be applied. After the Markov approximation we also applied a post-trace rotatating wave approximation. Here we obtain a difference for the extended and non-extended one-particle projection. While for the extended projection, the rotating wave approximation has the effect that it removes the remaining part of the Lamb-shift term, in the case of the non-extended projection there is still one term left in the Lamb-shift contribution. The dissipator term is not affected by this approximation because those terms that would be potentially affected have been already removed by the one-particle projection.  \nAs the first application of the renormalised one-particle master equation we discuss the evolution of the population and compare our results to the one in [11]. We demonstrate that although there exists a non-trivial evolution of the populations before renormalistion, which is consistent with the results obtained in [11] for the effective model of a photon if we choose the non-extended one-particle projection and specialise to the case of a vanishing temperature parameter, after renormalisation the evolution of the populations, that is the diagonal elements of the effective density matrix, becomes trivial. This shows the relevance of a renormalisation when analysing physical effects with the master equation, as the non-trivial dynamics of the populations and hence the effect discussed in [11] is removed by the renormalisation applied here. Hence, we conclude that if renormalisation is taken into account, we have no physically interesting effect in the evolution of the populations. In addition from our results for the evolution of the populations we learn that after applying the Markov approximation to the renormalised one-particle master equation, the extended and non-extended one-particle projection yield the same result.  \nIn addition, we analyse the form of the one-particle master equation in the non- and ultrarelativistic limit. In these cases, the dissipator obtains a rather simple form which is similar to the ones in [9] for a non-relativistic scalar particle and in [12] for a photon. While the detailed structure of the dissipator in [9] differs from ours, [12] has the same one we have.  \nAnother application we present is the comparison with the model in [27], where a quantum mechanical toy model for gravitationally induced decoherence, based on the model in [28], is investigated in the context of neutrino oscillations. Considering the renormalised one-particle master equation after applying the Markov and rotating wave approximation and in the special case of the ultra-relativistic limit, we obtain the renormalised Lindblad equation derived from the underlying quantum mechanical microscopic model in [27]. This comparison allows to fix the coupling constant in the interaction Hamiltonian that was an open parameter in [27]. As already discussed in [27],  \nwhere a first estimate of this coupling constant was given as well as the discussions in [12, 62-65] the value is fairly tiny. Whether more sophisticated models for gravitationally induced decoherence will involve larger values of the coupling parameter is an interesting question to discuss in future research. In addition, in contrast to the quantum mechanical model [27], no spectral density, which can be understood as a kind of continuum limit of the frequencies of the harmonic oscillators in the environment of the model in [27], has to be introduced into the field-theoretical model in this work, but the corresponding integrals are automatically determined from the underlying action of the model. The results obtained here also allow an extension of the model in [27], where the plane wave approach for the neutrinos is used, to models including wave packets, which can be implemented by an appropriate choice of the one-particle density matrix. We plan to consider an analysis of such a model in future work.  \nThere are several directions for further extensions and generalisations of the results obtained in this article. First, one can consider field-theoretical models where an operator ordering is chosen for which the self-interaction term of the scalar field does not vanish, as in [9]. In this case, we do not expect the renormalisation procedure to change, since we expect the same identification of non-covariant and covariant Feynman diagrams to hold for the self-energy, with the only difference that the second diagram in the sum of non-covariant diagrams does not vanish, as is the case in the model considered here. Instead of introducing the covariant Feynman rules in order to apply the renormalisation, one can also study if and how the methods from [66] for one-loop renormalisation of QED in Coulomb gauge can be employed to perform the renormalisation directly based on the non-covariant Feynman rules. Furthermore, a more general analysis of the question under which criteria the Markov approximation can be applied would be interesting. In particular, whether it is also possible to formulate conditions that depend only on the properties of the environment, as it was done in [54] for open quantum models with finitely many degrees of freedom. Also a study of the applicability and effect of different versions of the Markov approximation for gravitationally induced decoherence models is an interesting point, as some of them directly yield a Lindblad equation and therefore remove the necessity to apply an additional rotating wave approximation. Such modified Markov approximations is for instance the one discussed in [67, 68]. There, in addition to the density matrix, the operators of the system in the Schrodinger picture, which appear in the double commutator in the master equation, are Taylor expanded around the final time. Another possibility is a generalisation to field theory of the approximation for quantum mechanical models in [69], where one obtains a completely positive master equation by replacing an arithmetic mean of the spectral density by a geometric one.  \nAlthough there are quantum field theoretical models for gravitationally induced decoherence for scalar fields [1, 8, 9], photons [12, 16] or generic bosons [11], a model that includes fermions is still missing in the literature. Such a model would also be interesting in connection with the application to neutrino oscillations. Another possible generalisation is to consider quantisations other than Fock quantisations for the one-particle model, see for example [70], where a quantisation inspired by loop quantum gravity was used in the framework of polymerised quantum mechanics to formulate an open scattering model.', 'ACKNOWLEDGMENTS': 'MJF and KG would like to thank Michael Kobler for the valuable discussions at an early stage of the project. MJF would like to thank Renata Ferrero and Roman Kemper for the fruitful discussions on some computational aspects of the project. MJF thanks the Heinrich-Boll foundation for financial support.', 'Appendix A.I ONE-PARTICLE PROJECTION OF THE MASTER EQUATION': 'In this section, we explicitly carry out the projection of the master equation (2.1) on the oneparticle sector. As discussed in the main text, we proceed by inserting the density matrix of a single particle  \nρ 1 ( t ) = ∫ d 3 u ∫ d 3 v ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v (A.I.1)  \ninto the master equation (2.1), where ⃗u and ⃗v denote labels of the momentum and then by neglecting all contributions that do not preserve the one-particle subspace. The projection of the Hamiltonian H S of the free scalar field and the gravitational self-interaction of the scalar field U have already been discussed in the main text. Here we focus on the remaining terms, the Lamb-shift-like Hamiltonian H LS and the dissipator D . First, we give the detailed expressions for these two contributions of the master equation in subsection A.I.1 and then project them separately on the one-particle space in subsection A.I.2.', 'A.I.1 Individual contributions to the master equation': "In [1] the original dissipator was decomposed into two parts by separating a purely complex contribution entering the master equation analoguously to H S and U in terms of a commutator with ρ S . This term was then defined to be the so-called Lamb-shift Hamiltonian that contributes to the unitary evolution of the density matrix and the remaining terms were identified as a new dissipator term D that has a form similar to the first standard form (see e.g. [2] for the different forms of dissipators). The Lamb-shift Hamiltonian is given in [1] as  \nH LS := 1 2 ∫ d 3 k d 3 p d 3 l (2 π ) 6 2 ∑ r ; a,b S ab ( ⃗ p, ⃗ l ; ⃗ k, t ) j a r ( ⃗ k, ⃗p ) † j b r ( ⃗ k, ⃗ l ) , (A.I.2)  \nand the dissipator as  \nD [ ρ S ] := κ 2 ∫ d 3 k d 3 p d 3 l (2 π ) 6 2 ∑ r ; a,b R ab ( ⃗ p, ⃗ l ; ⃗ k, t ) ( j b r ( ⃗ k, ⃗ l ) ρ S ( t ) j a r ( ⃗ k, ⃗p ) † -1 2 { j a r ( ⃗ k, ⃗p ) † j b r ( ⃗ k, ⃗ l ) , ρ S ( t ) } ) , (A.I.3)  \nwhere S ab and R ab are coefficient functions and the j a r are operators containing combinations of scalar field's creation and annihilation operators respectively. These quantities will be defined below. From the structure of H LS and D follows that three different kinds of terms have to be evaluated for the projection:  \nj a r ( ⃗ k, ⃗p ) † j b r ( ⃗ k, ⃗ l ) ρ 1 ( t ) ︸ ︷︷ ︸ ( I ) and ρ 1 ( t ) j a r ( ⃗ k, ⃗p ) † j b r ( ⃗ k, ⃗ l ) ︸ ︷︷ ︸ ( II ) and j b r ( ⃗ k, ⃗ l ) ρ 1 ( t ) j a r ( ⃗ k, ⃗p ) † ︸ ︷︷ ︸ ( III ) . (A.I.4)  \nThe coefficient functions S and R are defined as  \nS ab ( ⃗ p, ⃗ l ; ⃗ k, t ) := 1 2Ω k [ [ N (Ω k ) + 1] { e -i (Ω k + ω b ( ⃗ k, ⃗ l )) t -1 Ω k + ω b ( ⃗ k, ⃗ l ) + e i (Ω k + ω a ( ⃗ k,⃗p )) t -1 Ω k + ω a ( ⃗ k, ⃗p ) } -N (Ω k ) { e i (Ω k -ω b ( ⃗ k, ⃗ l )) t -1 Ω k -ω b ( ⃗ k, ⃗ l ) + e -i (Ω k -ω a ( ⃗ k,⃗p )) t -1 Ω k -ω a ( ⃗ k, ⃗p ) }] (A.I.5)  \nand  \nR ab ( ⃗ p, ⃗ l ; ⃗ k, t ) := i Ω k [ [ N (Ω k ) + 1] { e -i (Ω k + ω b ( ⃗ k, ⃗ l )) t -1 Ω k + ω b ( ⃗ k, ⃗ l ) -e i (Ω k + ω a ( ⃗ k,⃗p )) t -1 Ω k + ω a ( ⃗ k, ⃗p ) } -N (Ω k ) { e i (Ω k -ω b ( ⃗ k, ⃗ l )) t -1 Ω k -ω b ( ⃗ k, ⃗ l ) -e -i (Ω k -ω a ( ⃗ k,⃗p )) t -1 Ω k -ω a ( ⃗ k, ⃗p ) }] (A.I.6)  \nwith the Bose-Einstein distribution  \nN (Ω k ) =: N ( k ) = 1 e βk -1 , (A.I.7)  \nwhere we use the graviton frequency Ω k = √ ⃗ k 2 =: k and β = 1 k B Θ with the Boltzmann constant k B and a temperature parameter Θ that determines the Gibbs state of the gravitational waves environment. Furthermore we have, defining ω k := √ m 2 + ⃗ k 2 with the scalar field's mass m :  \nj 1 r ( ⃗ k, ⃗p ) := a † p a k + p 1 2 √ ω p ω k + p [ p a p b [ P -r ( ⃗ k )] a b ] ω 1 ( ⃗ k, ⃗p ) := ω p -ω k + p (A.I.8)  \nj 2 r ( ⃗ k, ⃗p ) := a † -p -k a -p 1 2 √ ω p ω k + p [ p a p b [ P -r ( ⃗ k )] a b ] ω 2 ( ⃗ k, ⃗p ) := ω k + p -ω p (A.I.9)  \nj 3 r ( ⃗ k, ⃗p ) := a -p a k + p 1 2 √ ω p ω k + p [ p a p b [ P -r ( ⃗ k )] a b ] ω 3 ( ⃗ k, ⃗p ) := -ω p -ω k + p (A.I.10)  \nj 4 r ( ⃗ k, ⃗p ) := a † p a † -k -p 1 2 √ ω p ω k + p [ p a p b [ P -r ( ⃗ k )] a b ] ω 4 ( ⃗ k, ⃗p ) := ω p + ω k + p . (A.I.11)  \nHere, a k and a † k are annihilation and creation operator valued distributions acting on the underlying bosonic Fock space for a scalar field in the standard way. [ P ± r ( ⃗ k )] a b are the projectors on the individual transverse modes defined via  \n[ P ± ( ⃗ k )] b a := m a ( ± ⃗ k ) m b ( ± ⃗ k ) , (A.I.12)  \nwhere { ⃗ k | ⃗ k | , m ( ⃗ k ) , m ( ⃗ k ) } form an orthonormal basis of R 3 , details can be found in [1]. In what follows, we only need the properties of orthonormality of the three basis elements and the fact that the symmetric transverse-traceless projector is  \nP abcd ( ⃗ k ) = ∑ r ∈{±} [ P r ( ⃗ k )] a i [ P -r ( ⃗ k )] c j δ ib δ cd = 1 2 [ P ac ( ⃗ k ) P bd ( ⃗ k ) + P ad ( ⃗ k ) P bc ( ⃗ k ) -P ab ( ⃗ k ) P cd ( ⃗ k )] .  \n(A.I.13)  \nwith the transverse projector  \nP ab ( ⃗ k ) = δ ab -k a k b ⃗ k 2 . (A.I.14)  \nNow we can proceed to project H LS and D on the one-particle space.", 'A.I.2 Computation of the one-particle projection': 'We start with the evaluation of the three kinds of terms (I)-(III) and consider all possible combinations ( a, b ) that give a one-particle state after application. To keep track of all different combination in the three cases, we do this using a table. Considering only the creation and annihilation operator valued distributions a ( † ) k in the j -operators, we find:  \n- (I) j a r ( ⃗ k, ⃗p ) † j b r ( ⃗ k, ⃗ l ) a † u | 0 ⟩ ⟨ 0 | a v\n- (1,1) a † k + p a p a † l a l + k a † u | 0 ⟩ ⟨ 0 | a v = δ 3 ( ⃗ l + ⃗ k -⃗u ) a † k + p a p a † l | 0 ⟩ ⟨ 0 | a v = δ 3 ( ⃗ p + ⃗ k -⃗u ) δ 3 ( ⃗ p -⃗ l ) a † u | 0 ⟩ ⟨ 0 | a v\n- (1,2) a † k + p a p a † -l -k a -l a † u | 0 ⟩ ⟨ 0 | a v = δ 3 ( ⃗ l + ⃗u ) δ 3 ( ⃗ p + ⃗ l + ⃗ k ) a † u | 0 ⟩ ⟨ 0 | a v\n- (2,1) a † -p a -p -k a † l a k + l a † u | 0 ⟩ ⟨ 0 | a v = δ 3 ( ⃗ k + ⃗ l -⃗u ) δ 3 ( ⃗ p + ⃗ l + ⃗ k ) a † u | 0 ⟩ ⟨ 0 | a v\n- (2,2) a † -p a -p -k a † -l -k a -l a † u | 0 ⟩ ⟨ 0 | a v = δ 3 ( ⃗ l + ⃗u ) δ 3 ( ⃗ p -⃗ l ) a † u | 0 ⟩ ⟨ 0 | a v\n- (4,4) δ P a -k -p a p a † l a † -k -l a † u | 0 ⟩ ⟨ 0 | a v = δ P ( a ) see below\n- (II) a † u | 0 ⟩ ⟨ 0 | a v j a r ( ⃗ k, ⃗p ) † j b r ( ⃗ k, ⃗ l )\n- (1,1) a † u | 0 ⟩ ⟨ 0 | a v a † k + p a p a † l a l + k = δ 3 ( ⃗ p + ⃗ k -⃗v ) δ 3 ( ⃗ p -⃗ l ) a † u | 0 ⟩ ⟨ 0 | a v\n- (1,2) a † u | 0 ⟩ ⟨ 0 | a v a † k + p a p a † -l -k a -l = δ 3 ( ⃗ k + ⃗ p -⃗v ) δ 3 ( ⃗ p + ⃗ l + ⃗ k ) a † u | 0 ⟩ ⟨ 0 | a v\n- (2,1) a † u | 0 ⟩ ⟨ 0 | a v a † -p a -p -k a † l a k + l = δ 3 ( ⃗ p + ⃗v ) δ 3 ( ⃗ p + ⃗ l + ⃗ k ) a † u | 0 ⟩ ⟨ 0 | a v\n- (2,2) a † u | 0 ⟩ ⟨ 0 | a v a † -p a -p -k a † -l -k a -l = δ 3 ( ⃗ p + ⃗v ) δ 3 ( ⃗ p -⃗ l ) a † u | 0 ⟩ ⟨ 0 | a v\n- (4,4) δ P a † u | 0 ⟩ ⟨ 0 | a v a -k -p a p a † l a † -k -l = δ P (b) see below', '(III) j b r ( ⃗ k, ⃗ l ) a † u | 0 ⟩ ⟨ 0 | a v j a r ( ⃗ k, ⃗p ) †': "- (1,1) a † l a l + k a † u | 0 ⟩ ⟨ 0 | a v a † k + p a p = δ 3 ( ⃗ l + ⃗ k -⃗u ) δ 3 ( ⃗ k + ⃗ p -⃗v ) a † u -k | 0 ⟩ ⟨ 0 | a v -k\n- (1,2) a † -l -k a -l a † u | 0 ⟩ ⟨ 0 | a v a † k + p a p = δ 3 ( ⃗ k + ⃗ p -⃗v ) δ 3 ( ⃗u + ⃗ l ) a † u -k | 0 ⟩ ⟨ 0 | a v -k\n- (2,1) a † l a k + l a † u | 0 ⟩ ⟨ 0 | a v a † -p a -p -k = δ 3 ( ⃗ p + ⃗v ) δ 3 ( ⃗ l + ⃗ k -⃗u ) a † u -k | 0 ⟩ ⟨ 0 | a v -k\n- (2,2) a † -l -k a -l a † u | 0 ⟩ ⟨ 0 | a v a † -p a -p -k = δ 3 ( ⃗ p + ⃗v ) δ 3 ( ⃗u + ⃗ l ) a † u -k | 0 ⟩ ⟨ 0 | a v -k  \nThe expressions for (a) and (b) are calculated separately because they include the vacuum bubbles mentioned in the main text that require a renormalisation. By applying the commutators and commuting all annihilation operators towards the vacuum state one obtains:  \n( a ) = a -k -p a p a † l a † -k -l a † u | 0 ⟩ ⟨ 0 | a v =[ δ 3 ( ⃗ p + ⃗ k + ⃗ l ) δ 3 ξ ( ⃗ p + ⃗ k + ⃗ l ) + δ 3 ( ⃗ p -⃗ l ) δ 3 ξ ( ⃗ p -⃗ l ) + δ 3 ( ⃗ p + ⃗ k + ⃗ l ) δ 3 ( ⃗u + ⃗ k + ⃗ p ) + δ 3 ( ⃗ p -⃗ l ) δ 3 ( ⃗u + ⃗ k + ⃗ p ) + δ 3 ( ⃗ p -⃗u ) δ 3 ( ⃗ p -⃗ l ) + δ 3 ( ⃗ p -⃗u ) δ 3 ( ⃗ k + ⃗ p + ⃗ l )] a † u | 0 ⟩ ⟨ 0 | a v , (A.I.15)  \n( b ) = a † u | 0 ⟩ ⟨ 0 | a v a -k -p a p a † l a † -k -l =[ δ 3 ( ⃗ p + ⃗ k + ⃗ l ) δ 3 ξ ( ⃗ p + ⃗ k + ⃗ l ) + δ 3 ( ⃗ p -⃗ l ) δ 3 ξ ( ⃗ p -⃗ l ) + δ 3 ( ⃗ p + ⃗ k + ⃗ l ) δ 3 ( ⃗v + ⃗ k + ⃗ l ) + δ 3 ( ⃗ p -⃗ l ) δ 3 ( ⃗v + ⃗ k + ⃗ l ) + δ 3 ( ⃗v -⃗ l ) δ 3 ( ⃗ p -⃗ l ) + δ 3 ( ⃗v -⃗ l ) δ 3 ( ⃗ k + ⃗ p + ⃗ l )] a † u | 0 ⟩ ⟨ 0 | a v . (A.I.16)  \nNote that under the map ( ⃗u, ⃗p, ⃗ l ) ↔ ( ⃗v, ⃗ l, ⃗p ) we can get from (a) to (b) and vice versa. An important remark here is that the first two terms in both expressions contain the square of a Dirac delta distribution. This is a problematic term since, as can be shown, the corresponding integral over this expression still diverges when for the individual delta distributions the regularised version is considered and the regulator is removed after the integration is performed. We will deal with this issue further below in this subsection where we renormalise the density matrix to handle the divergent contributions involved. For now, we replace one of the two delta distributions by a function including a regulator δ 3 ( ⃗ k ) → δ 3 ξ ( ⃗ k ), where the regulator is sent to zero after performing the corresponding integrations.  \nIn a next step, we evaluate the entire expressions appearing in the Lamb-shift Hamiltonian  \nand the dissipator. For each term, where the delta distributions resolve two integrals, that is for all terms but the diverging ones, we choose to resolve the integrals over the range of the variables ⃗ p and ⃗ l respectively. Then we are left with the integration over the variables ⃗ k , ⃗u and ⃗v as well as with the sum over the polarisation labels r . Considering for V ab which is either R ab or S ab , that were defined above in (A.I.6) and (A.I.5), the three expressions  \n(I): ∑ r ∫ d 3 p ∫ d 3 l V ab ( ⃗ p, ⃗ l ; ⃗ k, t ) j a r ( ⃗ k, ⃗p ) † j b r ( ⃗ k, ⃗ l ) a † u | 0 ⟩ ⟨ 0 | a v (A.I.17)  \n(II): ∑ r ∫ d 3 p ∫ d 3 l V ab ( ⃗ p, ⃗ l ; ⃗ k, t ) a † u | 0 ⟩ ⟨ 0 | a v j a r ( ⃗ k, ⃗p ) † j b r ( ⃗ k, ⃗ l ) (A.I.18)  \n(III): ∑ r ∫ d 3 p ∫ d 3 l V ab ( ⃗ p, ⃗ l ; ⃗ k, t ) j a r ( ⃗ k, ⃗p ) † a † u | 0 ⟩ ⟨ 0 | a v j b r ( ⃗ k, ⃗ l ) , (A.I.19)  \nthat appear in (A.I.2) and (A.I.3), we obtain:  \n- (I) ∑ r ∫ d 3 p ∫ d 3 l V ab ( ⃗ p, ⃗ l ; ⃗ k, t ) j a r ( ⃗ k, ⃗p ) † j b r ( ⃗ k, ⃗ l ) a † u | 0 ⟩ ⟨ 0 | a v\n- (1,1) V 11 ( ⃗u -⃗ k, ⃗u -⃗ k ; ⃗ k, t ) ρ ( ⃗u, ⃗v ) 1 4 ω u ω u -k P u ( ⃗ k ) a † u | 0 ⟩ ⟨ 0 | a v\n- (1,2) V 12 ( ⃗u -⃗ k, -⃗u ; ⃗ k, t ) ρ ( ⃗u, ⃗v ) 1 4 ω u ω u -k P u ( ⃗ k ) a † u | 0 ⟩ ⟨ 0 | a v\n- (2,1) V 21 ( -⃗u, ⃗u -⃗ k ; ⃗ k, t ) ρ ( ⃗u, ⃗v ) 1 4 ω u ω u -k P u ( ⃗ k ) a † u | 0 ⟩ ⟨ 0 | a v\n- (2,2) V 22 ( -⃗u, -⃗u ; ⃗ k, t ) ρ ( ⃗u, ⃗v ) 1 4 ω u ω u -k P u ( ⃗ k ) a † u | 0 ⟩ ⟨ 0 | a v\n- (4,4) δ P (a) (see below)  \n(II) ∑ r ∫ d 3 p ∫ d 3 l V ab ( ⃗ p, ⃗ l ; ⃗ k, t ) a † u | 0 ⟩ ⟨ 0 | a v j a r ( ⃗ k, ⃗p ) † j b r ( ⃗ k, ⃗ l )  \n- (1,1) V 11 ( ⃗v -⃗ k,⃗v -⃗ k ; ⃗ k, t ) ρ ( ⃗u, ⃗v ) 1 4 ω v ω v -k P v ( ⃗ k ) a † u | 0 ⟩ ⟨ 0 | a v\n- (1,2) V 12 ( ⃗v -⃗ k, -⃗v ; ⃗ k, t ) ρ ( ⃗u, ⃗v ) 1 4 ω v ω v -k P v ( ⃗ k ) a † u | 0 ⟩ ⟨ 0 | a v\n- (2,1) V 21 ( -⃗v, ⃗v -⃗ k ; ⃗ k, t ) ρ ( ⃗u, ⃗v ) 1 4 ω v ω v -k P v ( ⃗ k ) a † u | 0 ⟩ ⟨ 0 | a v\n- (2,2) V 22 ( -⃗v, -⃗v ; ⃗ k, t ) ρ ( ⃗u, ⃗v ) 1 4 ω v ω v -k P v ( ⃗ k ) a † u | 0 ⟩ ⟨ 0 | a v\n- (4,4) δ P (b) (see below)  \n(III) ∑ r ∫ d 3 p ∫ d 3 l V ab ( ⃗ p, ⃗ l ; ⃗ k, t ) j a r ( ⃗ k, ⃗p ) † a † u | 0 ⟩ ⟨ 0 | a v j b r ( ⃗ k, ⃗ l )  \n- (1,1) V 11 ( ⃗v -⃗ k, ⃗u -⃗ k ; ⃗ k, t ) ρ ( ⃗u, ⃗v ) 1 4 √ ω v ω v -k ω u ω u -k P u,v ( ⃗ k ) a † u -k | 0 ⟩ ⟨ 0 | a v -k\n- (1,2) V 12 ( ⃗v -⃗ k, -⃗u ; ⃗ k, t ) ρ ( ⃗u, ⃗v ) 1 4 √ ω v ω v -k ω u ω u -k P u,v ( ⃗ k ) a † u -k | 0 ⟩ ⟨ 0 | a v -k\n- (2,1) V 11 ( -⃗v, ⃗u -⃗ k ; ⃗ k, t ) ρ ( ⃗u, ⃗v ) 1 4 √ ω v ω v -k ω u ω u -k P u,v ( ⃗ k ) a † u -k | 0 ⟩ ⟨ 0 | a v -k  \n(2,2) V 22 ( -⃗v, -⃗u ; ⃗ k, t ) ρ ( ⃗u, ⃗v ) 1 4 √ ω v ω v -k ω u ω u -k P u,v ( ⃗ k ) a † u -k | 0 ⟩ ⟨ 0 | a v -k ,  \nwhere we defined 12  \nP u ( ⃗ k ) := ∑ r [ ⃗u · ⃗ m ( -r ⃗ k )] 2 [ ⃗u · ⃗ m ( r ⃗ k )] 2 = 2[ ⃗u · ⃗ m ( -⃗ k )] 2 [ ⃗u · ⃗ m ( ⃗ k )] 2 = P abcd ( ⃗ k ) u a u b u c u d (A.I.20)  \nP u,v ( ⃗ k ) := ∑ r [ ⃗u · ⃗ m ( -r ⃗ k )] 2 [ ⃗v · ⃗ m ( r ⃗ k )] 2 = P abcd ( ⃗ k ) u a u b v c v d . (A.I.21)  \nNext, we evaluate in more detail the coefficient functions S and R which were defined in (A.I.5) and (A.I.6). Firstly excluding all terms from the extended projection, i.e. the ones arising from the combination (4 , 4), it turns out that the coefficients appearing in ( I ) -( III ) are equal to each other in each group:  \nS 11 ( ⃗u -⃗ k, ⃗u -⃗ k ; ⃗ k, t ) = S 12 ( ⃗u -⃗ k, -⃗u ; ⃗ k, t ) = S 21 ( -⃗u, ⃗u -⃗ k ; ⃗ k, t ) = S 22 ( -⃗u, -⃗u ; ⃗ k, t ) =: S ( -) 1 P ( ⃗u, ⃗ k, t ) = 1 Ω k { N ( k ) + 1 Ω k + ω u -k -ω u (cos [(Ω k + ω u -k -ω u ) t ] -1) -N ( k ) Ω k -ω u -k + ω u (cos [(Ω k -ω u -k + ω u ) t ] -1) } , (A.I.22)  \nR 11 ( ⃗u -⃗ k, ⃗u -⃗ k ; ⃗ k, t ) = R 12 ( ⃗u -⃗ k, -⃗u ; ⃗ k, t ) = R 21 ( -⃗u, ⃗u -⃗ k ; ⃗ k, t ) = R 22 ( -⃗u, -⃗u ; ⃗ k, t ) =: 2 R ( -) 1 P ( ⃗u, ⃗ k, t ) = 2 Ω k { N ( k ) + 1 Ω k + ω u -k -ω u sin [(Ω k + ω u -k -ω u ) t ] + N ( k ) Ω k -ω u -k + ω u sin [(Ω k -ω u -k + ω u ) t ] } (A.I.23)  \nand  \nR 11 ( ⃗v -⃗ k, ⃗u -⃗ k ; ⃗ k, t ) = R 12 ( ⃗v -⃗ k, -⃗u ; ⃗ k, t ) = R 21 ( -⃗v, ⃗u -⃗ k ; ⃗ k, t ) = R 22 ( -⃗v, -⃗u ; ⃗ k, t ) =: 2 R (2) 1 P ( ⃗u, ⃗v, ⃗ k, t ) = i Ω k { [ N ( k ) + 1] ( e -i (Ω k + ω u -k -ω u ) t -1 Ω k + ω u -k -ω u -e i (Ω k + ω v -k -ω v ) t -1 Ω k + ω v -k -ω v ) -N ( k ) ( e i (Ω k -ω u -k + ω u ) t -1 Ω k -ω u -k + ω u -e -i (Ω k -ω v -k + ω v ) t -1 Ω k -ω v -k + ω v )} . (A.I.24)  \nWith these expressions, we can rewrite the equations for the Lamb-shift Hamiltonian and the dissipator that contain the contributions of ( a, b ) ∈ { 1 , 2 } :  \n-i [ H { 1 , 2 } LS , ρ 1 ( t )] = -iκ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P u ( ⃗ k ) ω u -k ω u S ( -) 1 P ( ⃗u, ⃗ k, t ) + iκ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P v ( ⃗ k ) ω v -k ω v S ( -) 1 P ( ⃗v, ⃗ k, t ) (A.I.25)  \nas well as  \nD { 1 , 2 } [ ρ 1 ( t )] = -κ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P u ( ⃗ k ) ω u -k ω u R ( -) 1 P ( ⃗u, ⃗ k, t ) -κ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P v ( ⃗ k ) ω v -k ω v R ( -) 1 P ( ⃗v, ⃗ k, t ) + κ ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u -k | 0 ⟩ ⟨ 0 | a v -k P u,v ( ⃗ k ) √ ω u -k ω u ω v -k ω v R (2) 1 P ( ⃗u, ⃗v, ⃗ k, t ) . (A.I.26)  \nIt remains to deal with the (4 , 4)-terms that arise when using the extended projection. The four summands from (a) and (b) without the terms containing the ξ -regulator yield  \n-i [ H { 4 , 4 } LS , ρ 1 ( t )] = -iκ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P u ( ⃗ k ) ω u + k ω u ˜ S (+) 1 P ( ⃗u, ⃗ k, t ) δ P + iκ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P v ( ⃗ k ) ω v + k ω v ˜ S (+) 1 P ( ⃗v, ⃗ k, t ) δ P (A.I.27)  \nwith  \n˜ S (+) 1 P ( ⃗u, ⃗ k, t ) := S 44 ( -⃗u -⃗ k, ⃗u ; ⃗ k, t ) = S 44 ( -⃗u -⃗ k, -⃗u -⃗ k ; ⃗ k, t ) = S 44 ( ⃗u, ⃗u ; ⃗ k, t ) = S 44 ( ⃗u, -⃗u -⃗ k ; ⃗ k, t ) = 1 Ω k { N ( k ) + 1 Ω k + ω u + k + ω u (cos [(Ω k + ω u + k + ω u ) t ] -1) -N ( k ) Ω k -ω u + k -ω u (cos [(Ω k -ω u + k -ω u ) t ] -1) } . (A.I.28)  \nAfter a substitution ⃗ k →-⃗ k in the integration we find:  \n-i [ H { 4 , 4 } LS , ρ 1 ( t )] = -iκ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P u ( ⃗ k ) ω u -k ω u S (+) 1 P ( ⃗u, ⃗ k, t ) δ P + iκ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P v ( ⃗ k ) ω v -k ω v S (+) 1 P ( ⃗v, ⃗ k, t ) δ P (A.I.29)  \nwhere  \nS (+) 1 P ( ⃗u, ⃗ k, t ) := ˜ S (+) 1 P ( ⃗u, -⃗ k, t ) . (A.I.30)  \nThe same substitution and the definition  \nR (+) 1 P ( ⃗u, ⃗ k, t ) := 1 Ω k { N ( k ) + 1 Ω k + ω u -k + ω u sin [(Ω k + ω u -k + ω u ) t ] + N ( k ) Ω k -ω u -k -ω u sin [(Ω k -ω u -k -ω u ) t ] } (A.I.31)  \nlead to:  \nD { 4 , 4 } [ ρ 1 ] = -κ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P u ( ⃗ k ) ω u -k ω u R (+) 1 P ( ⃗u, ⃗ k, t ) δ P -κ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P v ( ⃗ k ) ω v -k ω v R (+) 1 P ( ⃗v, ⃗ k, t ) δ P . (A.I.32)  \nDefining additionally  \nS (1) 1 P ( ⃗u, ⃗ k, t ) := S ( -) 1 P ( ⃗u, ⃗ k, t ) + S (+) 1 P ( ⃗u, ⃗ k, t ) δ P (A.I.33)  \nR (1) 1 P ( ⃗u, ⃗ k, t ) := R ( -) 1 P ( ⃗u, ⃗ k, t ) + R (+) 1 P ( ⃗u, ⃗ k, t ) δ P (A.I.34)  \ngives us the opportunity to rewrite  \n-i [ H LS , ρ 1 ( t )] = -iκ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P u ( ⃗ k ) ω u -k ω u S (1) 1 P ( ⃗u, ⃗ k, t ) + iκ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P v ( ⃗ k ) ω v -k ω v S (1) 1 P ( ⃗v, ⃗ k, t ) (A.I.35)  \nas well as  \nD [ ρ 1 ] = -κ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P u ( ⃗ k ) ω u -k ω u R (1) 1 P ( ⃗u, ⃗ k, t ) -κ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P v ( ⃗ k ) ω v -k ω v R (1) 1 P ( ⃗v, ⃗ k, t ) + κ ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u -k | 0 ⟩ ⟨ 0 | a v -k P u,v ( ⃗ k ) √ ω u -k ω u ω v -k ω v R (2) 1 P ( ⃗u, ⃗v, ⃗ k, t ) . (A.I.36)  \nThis can be summarised as  \n-i [ H LS ,ρ 1 ( t )] + D [ ρ 1 ] = = -κ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P u ( ⃗ k ) ω u -k ω u [ R (1) 1 P ( ⃗u, ⃗ k, t ) + iS (1) 1 P ( ⃗u, ⃗ k, t )] -κ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u | 0 ⟩ ⟨ 0 | a v P v ( ⃗ k ) ω v -k ω v [ R (1) 1 P ( ⃗v, ⃗ k, t ) -iS (1) 1 P ( ⃗u, ⃗ k, t )] + κ ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ⃗u, ⃗v, t ) a † u -k | 0 ⟩ ⟨ 0 | a v -k P u,v ( ⃗ k ) √ ω u -k ω u ω v -k ω v R (2) 1 P ( ⃗u, ⃗v, ⃗ k, t ) . (A.I.37)  \nFinally, the terms in the extended projection that contain the ξ -regulator, that is the expressions arising from the first two terms in (A.I.15) and (A.I.16), are analysed. These terms are equal for (a) and (b) and read:  \nδ P lim ξ → 0 ∫ d 3 p d 3 l d 3 k d 3 ud 3 v (2 π ) 3 2 a † u | 0 ⟩ ⟨ 0 | a v V 44 ( ⃗ p, ⃗ l ; ⃗ k, t ) P p,l ( ⃗ k ) 4 √ ω p ω k + p ω l ω k + l ρ ( ⃗u, ⃗v, t ) · { δ 3 ( ⃗ p + ⃗ k + ⃗ l ) δ 3 ξ ( ⃗ p + ⃗ k + ⃗ l ) + δ 3 ( ⃗ p -⃗ l ) δ 3 ξ ( ⃗ p -⃗ l ) } = δ P lim ξ → 0 ∫ d 3 p d 3 k (2 π ) 3 2 V 44 ( ⃗ p, ⃗p ; ⃗ k, t ) P p ( ⃗ k ) 2 ω p ω k + p δ 3 ξ ( ⃗ 0) ρ 1 ( t ) , (A.I.38)  \nwhere we used that V 44 ( ⃗ p, ⃗p ; ⃗ k, t ) = V 44 ( ⃗ p, ⃗ -k -p ; ⃗ k, t ). Due to the equality of these terms for (a) and (b), they drop out of the Lamb-shift Hamiltonian and are only left in the dissipator term:  \nD div [ ρ 1 ] = -δ P κ 4 lim ξ → 0 ∫ d 3 k d 3 p (2 π ) 3 2 R 44 ( ⃗ p, ⃗p ; ⃗ k, t ) P p ( ⃗ k ) 2 ω p ω k + p δ 3 ξ ( ⃗ 0) ρ 1 ( t ) =: Z ( t ) ρ 1 ( t ) . (A.I.39)  \nWritten in this form, it becomes evident that they do act as a multiplicative constant and do not modify the state ρ 1 . Therefore they are nothing but vacuum bubbles expressed in QFT language. With the definition of Z ( t ) above the entire master equation for the single particle is given by  \n∂ ∂t ρ 1 ( t ) = -i [ H S + H LS , ρ 1 ( t )] + D [ ρ 1 ] + Z ( t ) ρ 1 ( t ) , (A.I.40)  \nwhere from now on the terms absorbed in Z ( t ) are dropped from the definition of D [ ρ 1 ]. We can see that the diverging term Z ( t ) can be absorbed by a renormalisation of the density matrix, likewise to a renormalisation of the wave function known from QFT:  \nρ 1 ( t ) → ρ ( ren ) 1 ( t ) := exp (∫ t 0 dt ' Z ( t ' ) ) ρ 1 ( t ) . (A.I.41)  \nIn terms of the renormalised density matrix, the one-particle master equation then reads:  \n∂ ∂t ρ ( ren ) 1 ( t ) = -i ∫ d 3 ud 3 v ρ ( ren ) ( ⃗u, ⃗v, t ) | ⃗u ⟩ ⟨ ⃗v | ( ω u -ω v ) -κ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ren ) ( ⃗u, ⃗v, t ) | ⃗u ⟩ ⟨ ⃗v | P u ( ⃗ k ) ω u -k ω u [ R (1) 1 P ( ⃗u, ⃗ k, t ) + iS (1) 1 P ( ⃗u, ⃗ k, t )] -κ 2 ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ren ) ( ⃗u, ⃗v, t ) | ⃗u ⟩ ⟨ ⃗v | P v ( ⃗ k ) ω v -k ω v [ R (1) 1 P ( ⃗v, ⃗ k, t ) -iS (1) 1 P ( ⃗v, ⃗ k, t )] + κ ∫ d 3 k d 3 ud 3 v (2 π ) 3 ρ ( ren ) ( ⃗u, ⃗v, t ) | ⃗u -⃗ k ⟩ ⟨ ⃗v -⃗ k | P u,v ( ⃗ k ) √ ω u -k ω u ω v -k ω v R (2) 1 P ( ⃗u, ⃗v, ⃗ k, t ) . (A.I.42)  \nTo simplify the notation, we drop the label ( ren ) from the density matrix from this point on. The complex combinations R ± iS can be evaluated further and one finds  \nR (1) 1 P ( ⃗u, ⃗ k, t ) ± iS (1) 1 P ( ⃗u, ⃗ k, t ) = = ± i Ω k { N ( k ) + 1 Ω k + ω u -k -ω u ( e ∓ i (Ω k + ω u -k -ω u ) t -1 ) -N ( k ) Ω k -ω u -k + ω u ( e ± i (Ω k -ω u -k + ω u ) t -1 ) + δ P N ( k ) + 1 Ω k + ω u -k + ω u ( e ∓ i (Ω k + ω u -k + ω u ) t -1 ) -δ P N ( k ) Ω k -ω u -k -ω u ( e ± i (Ω k -ω u -k -ω u ) t -1 ) } = ∫ t 0 dτ Ω k { [ N ( k ) + 1] e ∓ i (Ω k + ω u -k -ω u ) τ + N ( k ) e ± i (Ω k -ω u -k + ω u ) τ + δ P [ N ( k ) + 1] e ∓ i (Ω k + ω u -k + ω u ) τ + δ P N ( k ) e ± i (Ω k -ω u -k -ω u ) τ } . (A.I.43)  \nFor R (2) we have:  \nR (2) 1 P ( ⃗u, ⃗v, ⃗ k, t ) = i 2Ω k { [ N ( k ) + 1] ( e -i (Ω k + ω u -k -ω u ) t -1 Ω k + ω u -k -ω u -e i (Ω k + ω v -k -ω v ) t -1 Ω k + ω v -k -ω v ) -N ( k ) ( e i (Ω k -ω u -k + ω u ) t -1 Ω k -ω u -k + ω u -e -i (Ω k -ω v -k + ω v ) t -1 Ω k -ω v -k + ω v )} = ∫ t 0 dτ 2Ω k { [ N ( k ) + 1] e -i (Ω k + ω u -k -ω u ) τ + N ( k ) e i (Ω k -ω u -k + ω u ) τ +[ N ( k ) + 1] e i (Ω k + ω v -k -ω v ) τ + N ( k ) e -i (Ω k -ω v -k + ω v ) τ } . (A.I.44)  \nDefining then  \nC ( ⃗u, ⃗ k, t ) = ∫ t -t 0 0 dτ Ω k { [ N ( k ) + 1] e -i (Ω k + ω u -k -ω u ) τ + N ( k ) e i (Ω k -ω u -k + ω u ) τ } (A.I.45)  \nC P ( ⃗u, ⃗ k, t ) = ∫ t -t 0 0 dτ Ω k { [ N ( k ) + 1] e -i (Ω k + ω u -k + ω u ) τ + N ( k ) e i (Ω k -ω u -k -ω u ) τ } , (A.I.46)  \nwhere we have restored the initial time 13 t 0 that was set to 0 in [1], the one-particle master equation in momentum representation has the form  \n∂ ∂t ρ ( ⃗u, ⃗v, t ) = -iρ ( ⃗u, ⃗v, t ) ( ω u -ω v ) -κ 2 ∫ d 3 k (2 π ) 3 { P u ( ⃗ k ) ω u -k ω u [ C ( ⃗u, ⃗ k, t ) + δ P C P ( ⃗u, ⃗ k, t ) ] + P v ( ⃗ k ) ω v -k ω v [ C ∗ ( ⃗v, ⃗ k, t ) + δ P C ∗ P ( ⃗v, ⃗ k, t ) ] } ρ ( ⃗u, ⃗v, t ) + κ 2 ∫ d 3 k (2 π ) 3 P u,v ( ⃗ k ) √ ω u + k ω u ω v + k ω v { C ( ⃗u + ⃗ k, ⃗ k, t ) + C ∗ ( ⃗v + ⃗ k, ⃗ k, t ) } ρ ( ⃗u + ⃗ k,⃗v + ⃗ k, t ) . (A.I.47)", 'Appendix A.II EQUIVALENCE OF THE NON-COVARIANT AND COVARIANT FEYNMAN RULES FOR THE COULOMB-LIKE SCATTERING DIAGRAM': 'Here we show that the non-covariant and the covariant Feynman rules introduced in sections IV B and IVC respectively yield the same result for the Coulomb-like scattering diagram  \n<!-- image -->  \nThe left hand side of the equation is expressed in terms of the non-covariant Feynman rules and the first diagram, which we label A , reads  \nA = 1 2 ( -iκ ) 2 ˜ T ab ( p, q ) ˜ T cd ( -u, -v ) -i κ 1 k 2 P abcd ( ⃗ k ) = iκ 2 k 2 ˜ T ab ( p, q ) ˜ T cd ( -u, -v ) [ ( δ ac -k a k c ⃗ k 2 )( δ bd -k b k d ⃗ k 2 ) -1 2 ( δ ab -k a k b ⃗ k 2 )( δ cd -k c k d ⃗ k 2 ) ] = iκ k 2 [ 1 2 ˜ T ab ( p, q ) ˜ T ab ( -u, -v ) -1 4 ˜ T a a ( p, q ) ˜ T b b ( -u, -v ) + 1 4 ⃗ k 4 k a k b k c k d ˜ T ab ( p, q ) ˜ T cd ( -u, -v ) -1 ⃗ k 2 { ˜ T ab ( p, q ) ˜ T a c ( -u, -v ) k b k c -1 4 ˜ T ab ( p, q ) ˜ T c c ( -u, -v ) k a k b -1 4 ˜ T a a ( p, q ) ˜ T cd ( -u, -v ) k c k d } ] , (A.II.1)  \nwhere the overall factor of 1 2 arises due to the fact that it is a diagram of second order in the expansion of the Dyson series. Next, we can make use of energy-momentum conservation, as we assumed the scalar particles to be on-shell, which reads, as introduced in the main text in (4.21) and (4.22) for k = p + q = u + v :  \n1  \n2  \nk µ ˜ T µ 0 ( p, q ) = k 0 ˜ T 00 ( p, q ) + k a ˜ T a 0 ( p, q ) = 2 [ q 0 ( p 2 0 -⃗ p 2 -m 2 ) + p 0 ( q 2 0 -⃗q 2 -m 2 )] = 0 (A.II.2) k µ ˜ T µa ( p, q ) = k 0 ˜ T 0 a ( p, q ) + k b ˜ T ba ( p, q ) = 1 [ q a ( p 2 0 -⃗ p 2 -m 2 ) + p a ( q 2 0 -⃗q 2 -m 2 )] = 0 . (A.II.3)  \nFrom these we find  \nk a k b T ab ( p, q ) = -k 0 k b T 0 b ( p, q ) = ( k 0 ) 2 T 00 ( p, q )  \n˜ ˜ ˜ (A.II.4) k a ˜ T ab ( p, q ) = -k 0 ˜ T 0 b ( p, q ) (A.II.5) k a ˜ T 0 a ( p, q ) = -k 0 ˜ T 00 ( p, q ) (A.II.6)  \nwhich leads to  \nA = iκ k 2 [ 1 2 ˜ T ab ( p, q ) ˜ T ab ( -u, -v ) -1 4 ˜ T a a ( p, q ) ˜ T b b ( -u, -v ) + ( k 0 ) 4 4 ⃗ k 4 ˜ T 00 ( p, q ) ˜ T 00 ( -u, -v ) -( k 0 ) 2 ⃗ k 2 { ˜ T a 0 ( p, q ) ˜ T a 0 ( -u, -v ) -1 4 ˜ T 00 ( p, q ) ˜ T c c ( -u, -v ) -1 4 ˜ T a a ( p, q ) ˜ T 00 ( -u, -v ) } ] . (A.II.7)  \nCombining this with the expression for the second term in the above Feynman diagram for Coulomb-like scattering, which we call B and which reads  \nB = -iκ ⃗ k 2 [ -1 4 ˜ T 00 ( p, q ) ˜ T 00 ( -u, -v ) + ˜ T 0 a ( p, q ) ˜ T 0 b ( -u, -v ) ( δ ab -k a k b 4 ⃗ k 2 ) -1 4 ˜ T 00 ( p, q ) ˜ T a a ( -u, -v ) -1 4 ˜ T a a ( p, q ) ˜ T 00 ( -u, -v ) ] = -iκ k 2 [ -k 2 4 ⃗ k 2 ˜ T 00 ( p, q ) ˜ T 00 ( -u, -v ) + k 2 ⃗ k 2 ˜ T 0 a ( p, q ) ˜ T a 0 ( -u, -v ) -k 2 ( k 0 ) 2 4 ⃗ k 4 ˜ T 00 ( p, q ) ˜ T 00 ( -u, -v ) -k 2 4 ⃗ k 2 ˜ T 00 ( p, q ) ˜ T a a ( -u, -v ) -k 2 4 ⃗ k 2 ˜ T a a ( p, q ) ˜ T 00 ( -u, -v ) ] , (A.II.8)  \none can obtain  \nA + B = iκ k 2 [( ( k 0 ) 4 + k 2 ⃗ k 2 + k 2 ( k 0 ) 2 4 ⃗ k 4 ) ˜ T 00 ( p, q ) ˜ T 00 ( -u, -v ) -( ( k 0 ) 2 + k 2 ⃗ k 2 ) ˜ T 0 a ( p, q ) ˜ T a 0 ( -u, -v ) + 1 2 ˜ T ab ( p, q ) ˜ T ab ( -u, -v ) -1 4 ˜ T a a ( p, q ) ˜ T b b ( -u, -v ) + ( ( k 0 ) 2 + k 2 4 ⃗ k 2 ) ˜ T 00 ( p, q ) ˜ T c c ( -u, -v ) + ( ( k 0 ) 2 + k 2 4 ⃗ k 2 ) ˜ T a a ( p, q ) ˜ T 00 ( -u, -v ) ] = iκ k 2 [ 1 4 ˜ T 00 ( p, q ) ˜ T 00 ( -u, -v ) -˜ T 0 a ( p, q ) ˜ T a 0 ( -u, -v ) + 1 2 ˜ T ab ( p, q ) ˜ T ab ( -u, -v ) -1 4 ˜ T a a ( p, q ) ˜ T b b ( -u, -v ) + 1 4 ˜ T 00 ( p, q ) ˜ T c c ( -u, -v ) + 1 4 ˜ T a a ( p, q ) ˜ T 00 ( -u, -v ) ] . (A.II.9)  \nOn the other hand, we obtain with the covariant Feynman rules for the right side of the Coulomblike scattering diagram above 14 , which we name C :  \nC  \n= 1 2 iκ k 2 ˜ T µν ( p, q ) ˜ T ρσ ( -u, -v ) ( η µρ η νσ -1 2 η µν η ρσ ) = iκ 2 k 2 [ ˜ T µν ( p, q ) ˜ T µν ( -u, -v ) -1 2 ˜ T µ µ ( p, q ) ˜ T ρ ρ ( -u, -v ) ] = iκ 4 k 2 [ 2 ˜ T 00 ( p, q ) ˜ T 00 ( -u, -v ) + 2 ˜ T ab ( p, q ) ˜ T ab ( -u, -v ) + 4 ˜ T 0 a ( p, q ) ˜ T 0 a ( -u, -v ) -˜ T 0 0 ( p, q ) ˜ T 0 0 ( -u, -v ) -˜ T 0 0 ( p, q ) ˜ T c c ( -u, -v ) -˜ T a a ( p, q ) ˜ T 0 0 ( -u, -v ) -˜ T a a ( p, q ) ˜ T b b ( -u, -v ) ] = iκ k 2 [ 1 4 ˜ T 00 ( p, q ) ˜ T 00 ( -u, -v ) ˜ T 0 a ( p, q ) ˜ T 0 a ( -u, -v ) + 1 2 ˜ T ab ( p, q ) ˜ T ab ( -u, -v ) --1 4 ˜ T a a ( p, q ) ˜ T b b ( -u, -v ) + 1 4 ˜ T 00 ( p, q ) ˜ T c c ( -u, -v ) + 1 4 ˜ T a a ( p, q ) ˜ T 00 ( -u, -v ) ] . (A.II.10)  \nNote that we use the mostly plus signature of the metric, hence pulling a temporal index results in a sign change. By comparing (A.II.9) and (A.II.10) we can see that they are identical, therefore we indeed have that the non-covariant and the covariant Feynman rules produce the same result for Coulomb-like scattering.  \nWhen considering a loop diagram as discussed at the end of section IV C, then the momentum inside the loop is not on-shell. Due to this, in a similar calculation as the one shown in this appendix, there will remain correction terms to the relations in (A.II.4) - (A.II.6) that will prevent one from directly seeing the equivalence. However, as we discuss in the main text, as both sets can be derived from the same underlying Lagrangian, we expect that they yield the same physics.', 'Appendix A.III DETAILED COMPUTATION OF THE UV-RENORMALISATION': 'In this appendix we present the evaluation of the vacuum self-energy diagram from equation (4.27),  \nΠ vac ( u 2 ) = κ 2 µ ϵ ∫ d d k u 2 k 2 +2 m 2 uk -2 m 4 ( 1 + ϵ 4 ) [( k + u ) 2 + λ 2 -iϵ ][ k 2 + m 2 -iϵ ] , (A.III.1)  \nwhich encodes the self-energy in vacuum of the scalar particle of mass m caused by an internal triad particle of small auxiliary mass λ . For the evaluation, we follow the strategy in [53] for the regularisation of the QED self-energy at one-loop level.', 'A.III.1 Computation of the Loop integral': 'In a first step we use the identity  \n1 x -iϵ = i ∫ ∞ 0 dz e -iz ( x -iϵ ) (A.III.2)  \nin order to rewrite equation (A.III.1) as  \nΠ vac ( u 2 ) = -κ 2 µ ϵ ∫ d d k [ u 2 k 2 +2 m 2 uk -2 m 4 ( 1 + ϵ 4 )] · ∫ ∞ 0 dz 1 ∫ ∞ 0 dz 2 e -iz 1 ( k 2 + m 2 -iϵ ) e -iz 2 (( k + u ) 2 + λ 2 -iϵ ) = -κ 2 µ ϵ ∫ ∞ 0 dz 1 ∫ ∞ 0 dz 2 ∫ d d k [ u 2 k 2 +2 m 2 uk -2 m 4 ( 1 + ϵ 4 )] · e -i ( z 1 + z 2 ) ( k + z 2 u z 1 + z 2 ) 2 -[ ϵ ( z 1 + z 2 )+ iλ 2 z 2 + im 2 z 1 + iu 2 z 2 -i z 2 2 u 2 z 1 + z 2 ] . (A.III.3)  \nSubstituting k → k -z 2 u z 1 + z 2 then yields  \nΠ vac ( u 2 ) = -κ 2 µ ϵ ∫ ∞ 0 dz 1 ∫ ∞ 0 dz 2 ∫ d d k [ u 2 k 2 -2 u 2 z 2 uk z 1 + z 2 + z 2 2 u 4 ( z 1 + z 2 ) 2 +2 m 2 uk -2 m 2 z 2 u 2 z 1 + z 2 -2 m 4 ( 1 + ϵ 4 )] · e -i ( z 1 + z 2 ) k 2 -[ ϵ ( z 1 + z 2 )+ i ( λ 2 + u 2 ) z 2 + im 2 z 1 -i z 2 2 u 2 z 1 + z 2 ] . (A.III.4)  \n=  \n-  \n=  \n-  \nϵ  \nµ  \nϵ  \nµ  \ndz  \nlim  \ndz  \nlim  \nξ  \n→  \n0  \nξ  \n→  \n0  \n-  \nβ  \ndβ β e  \nDue to symmetry, all terms linear in k vanish and only terms of two different kinds remain for which the k -integrations can be performed directly:  \n∫ d d k e -i ( z 1 + z 2 ) k 2 = √ -πi z 1 + z 2 d (A.III.5) ∫ d d k k 2 e -i ( z 1 + z 2 ) k 2 = i d d ( z 1 + z 2 ) ∫ d d k e -i ( z 1 + z 2 ) k 2 = d 2 π √ -πi z 1 + z 2 d +2 . (A.III.6)  \nEmploying these, one can rewrite the self-energy as:  \nΠ vac ( u 2 ) = -κ 2 µ ϵ ∫ ∞ 0 dz 1 ∫ ∞ 0 dz 2 e -[ ϵ ( z 1 + z 2 )+ i ( λ 2 + u 2 ) z 2 + im 2 z 1 -i z 2 2 u 2 z 1 + z 2 ] · { √ -πi z 1 + z 2 d [ z 2 2 u 4 ( z 1 + z 2 ) 2 -2 m 2 z 2 u 2 z 1 + z 2 -2 m 4 ( 1 + ϵ 4 ) ] + u 2 d 2 π √ -πi z 1 + z 2 d +2 } . (A.III.7)  \nTo continue, we use 15  \n1 = lim ξ → 0 ∫ ∞ ξ dβ β δ ( 1 -z 1 + z 2 β ) . (A.III.8)  \nBy substituting z 1 → z 1 β and z 2 → z 2 β we then obtain  \nΠ vac ( u 2 ) = -κ 2 µ ϵ ∫ ∞ 0 dz 1 ∫ ∞ 0 dz 2 lim ξ → 0 ∫ ∞ ξ dβ δ (1 -( z 1 + z 2 )) β e -β [ ϵ ( z 1 + z 2 )+ i ( λ 2 + u 2 ) z 2 + im 2 z 1 -i z 2 2 u 2 z 1 + z 2 ] · {√ -πi β ( z 1 + z 2 ) d [ z 2 2 u 4 ( z 1 + z 2 ) 2 -2 m 2 z 2 u 2 z 1 + z 2 -2 m 4 ( 1 + ϵ 4 ) ] }  \nκ  \n2  \nκ  \n2  \n∫  \n0  \n∫  \n0  \n1  \n1  \n∫  \nξ  \n∫  \nξ  \n∞  \n[  \nϵ  \n+  \ni  \n(  \nλ  \n2  \n+  \nu  \n2  \n)  \nz  \n+  \nim  \n2  \n(1  \n-  \nz  \n)  \n-  \niz  \n2  \nu  \n2  \n]  \n·    √ -πi β d [ z 2 u 4 -2 m 2 zu 2 -2 m 4 ( 1 + ϵ 4 )] + u 2 d 2 π √ -πi β d +2     \n∞  \n[  \ndβ e  \n-  \niβ  \n-  \niϵ  \n+(  \nλ  \n2  \n+  \nu  \n2  \n)  \nz  \n+  \nm  \n2  \n(1  \n-  \nz  \n)  \n-  \nz  \n2  \nu  \n2  \n]  \n· { -π 2 β ϵ 2 β √ -iπ -ϵ [ z 2 u 4 -2 m 2 zu 2 -2 m 4 ( 1 + ϵ 4 )]  \n+ iπ 2 4 -ϵ 2 u 2 √ -iπ -ϵ β ϵ 2 β 2 } , (A.III.9)  \n̸  \n+ u 2 d 2 π √ -πi β ( z 1 + z 2 ) d +2  \nwhere in the last step we used d = 4 -ϵ . Next, we perform the β -integration, that is  \nlim ξ → 0 ∫ ∞ ξ dβ e -iβa β ϵ 2 β n , (A.III.10)  \nwhere we already set ϵ = 0 in the exponential, with a = ( λ 2 + u 2 ) z + m 2 (1 -z ) -z 2 u 2 and n ∈ { 1 , 2 } . To obtain the result, we apply the residue theorem. For a > 0 the contour can be closed by a quarter circle from ∞ to -i ∞ + ξ and a line from -i ∞ + ξ to ξ . With this, the pole at β = 0 can be avoided. The integral then becomes  \nlim ξ → 0 ∫ ∞ ξ dβ e -iβa β ϵ 2 β n = -lim ξ → 0 ∫ ∞ 0 dt e -t ( t + iξa ) ϵ 2 -n ( -i ) ϵ 2 -n +1 a -ϵ 2 + n -1 , (A.III.11)  \nwhere we substituted t = βa . Expanding the term ( t + iξa ) ϵ 2 -n for small ξ yields  \n( t + iξa ) ϵ 2 -n = t ϵ 2 -n + iξa ( ϵ 2 -n ) t ϵ 2 -n -1 + O ( ξ 2 ) . (A.III.12)  \nWith the definition of the Gamma function  \nΓ( z ) = ∫ ∞ 0 dt t z -1 e -t (A.III.13)  \nit then follows that  \nlim ξ → 0 ∫ ∞ ξ dβ e -iβa β ϵ 2 β n = ( ia ) n -1 -ϵ 2 [ Γ ( ϵ 2 -n +1 ) + lim ξ → 0 O ( ξ ) ] = ( ia ) n -1 -ϵ 2 Γ ( ϵ 2 -n +1 ) . (A.III.14)  \nFor a < 0, the contour is closed by a quarter circle from ∞ to i ∞ + ξ and a line from i ∞ + ξ to ξ . The result turns out to be the same as for the case a > 0. This thus yields  \nΠ vac ( u 2 ) = -κ 2 µ ϵ ∫ 1 0 dz { -π 2 √ -iπ -ϵ [ z 2 u 4 -2 m 2 zu 2 -2 m 4 ( 1 + ϵ 4 )] i -ϵ 2 · [( λ 2 + u 2 ) z + m 2 (1 -z ) -z 2 u 2 ] -ϵ 2 Γ ( ϵ 2 ) + iπ 2 ( 2 -ϵ 2 ) u 2 √ -iπ -ϵ i 1 -ϵ 2 [( λ 2 + u 2 ) z + m 2 (1 -z ) -z 2 u 2 ] 1 -ϵ 2 Γ ( ϵ 2 -1 ) } . (A.III.15)  \nNext, we expand all terms in ϵ and then perform the z -integration which results in  \nΠ vac ( u 2 ) = -π 2 κ ϵ [ 2 m 2 ( m 2 + u 2 ) + λ 2 u 2 ] + π 2 κ 2 u 2 { u 2 [2 m 2 ( m 2 + u 2 )( γ E -2) + λ 2 ( u 2 ( γ E -1) -m 2 )] + m 2 ( m 2 +2 u 2 + λ 2 ) √ ( m 2 + u 2 ) 2 -2 λ 2 ( m 2 -u 2 ) + λ 4 · [ arctanh ( m 2 + u 2 -λ 2 √ ( m 2 + u 2 ) 2 -2 λ 2 ( m 2 -u 2 ) + λ 4 ) -arctanh ( m 2 -u 2 -λ 2 √ ( m 2 + u 2 ) 2 -2 λ 2 ( m 2 -u 2 ) + λ 4 )] +ln( π )[ u 4 λ 2 +2 m 2 u 2 ( m 2 + u 2 )] + ln ( λ m ) [ m 6 -3 m 2 u 2 λ 2 -m 2 λ 4 ] +ln( λ )[3 u 2 m 4 +2 u 4 ( m 2 + λ 2 )] + ln( m )[ u 2 m 2 (2 u 2 + m 2 )] -ln( µ 2 )[2 u 2 m 2 ( m 2 + u 2 ) + λ 2 u 4 ] } + O ( ϵ ) (A.III.16)  \nGiven this, it becomes evident that the pole in ϵ arises from the term -π 2 κ ϵ [ 2 m 2 ( m 2 + u 2 ) + λ 2 u 2 ] , which yields in the limit of vanishing graviton mass -π 2 κ ϵ [ 2 m 2 ( m 2 + u 2 ) ] . A suitable counter term should remove this divergence, which is discussed in the main text, and we are left with the regularised version  \nΠ reg vac ( u 2 ) = π 2 κ 2 u 2 { u 2 [2 m 2 ( m 2 + u 2 )( γ E -2) + λ 2 ( u 2 ( γ E -1) -m 2 )] + m 2 ( m 2 +2 u 2 + λ 2 ) √ ( m 2 + u 2 ) 2 -2 λ 2 ( m 2 -u 2 ) + λ 4 · [ arctanh ( m 2 + u 2 -λ 2 √ ( m 2 + u 2 ) 2 -2 λ 2 ( m 2 -u 2 ) + λ 4 ) -arctanh ( m 2 -u 2 -λ 2 √ ( m 2 + u 2 ) 2 -2 λ 2 ( m 2 -u 2 ) + λ 4 )] +ln( π )[ u 4 λ 2 +2 m 2 u 2 ( m 2 + u 2 )] + ln ( λ m ) [ m 6 -3 m 2 u 2 λ 2 -m 2 λ 4 ] +ln( λ )[3 u 2 m 4 +2 u 4 ( m 2 + λ 2 )] + ln( m )[ u 2 m 2 (2 u 2 + m 2 )] -ln( µ 2 )[2 u 2 m 2 ( m 2 + u 2 ) + λ 2 u 4 ] } . (A.III.17)  \nIn order to determine the (arbitrary) finite part of the counter term, we pick the on-shell renormalisation condition. This requires the pole in the propagator to be at u 2 = -m 2 , that is it imposes the condition Π( u 2 = -m 2 ) ! = 0, and the fixing of its residue such that ∂ ∂u 2 Π( u 2 = -m 2 ) ! = 0. Together, these two conditions imply the following formula for the renormalised Π R vac :  \nΠ R vac ( u 2 ) := Π reg vac ( u 2 ) -Π reg vac ( -m 2 ) -( u 2 + m 2 ) ∂ ∂u 2 Π reg vac ( -m 2 ) . (A.III.18)  \nApplied to (A.III.16), one obtains for the additional terms:  \nΠ reg vac ( -m 2 ) = -π 2 κ 2 λ { ( m 2 -λ 2 ) √ λ 2 -4 m 2 [ arctanh ( 2 m 2 -λ 2 λ √ λ 2 -4 m 2 ) +arctanh ( λ √ λ 2 -4 m 2 )] + λ 3 ln ( m λ ) + m 2 λ [ γ E +ln ( πλ 5 µ 2 m 3 )] } . (A.III.19)  \nas well as  \n∂ ∂u 2 Π reg vac ( -m 2 ) = π 2 κ 2 m 2 { 5 m 4 λ +2 m 2 λ 3 -λ 5 √ λ 2 -4 m 2 [ arctanh ( λ 2 -2 m 2 λ √ λ 2 -4 m 2 ) -arctanh ( λ √ λ 2 -4 m 2 )] +(2 γ E -3) m 4 + m 2 λ 2 ( γ E -2 + ln( π )) + ( m 4 -λ 4 ) ln ( m λ ) +2 m 2 λ 2 ln( λ ) + 2 m 4 ln( mπλ ) -m 2 (2 m 2 + λ 2 ) ln ( µ 2 ) } . (A.III.20)  \nIf we had not introduced the small triad mass λ , then this last expression would be divergent, see for instance [52, 53]. Hence we continue to work now with Π R vac ( u 2 ). This concludes the discussion of the renormalisation of the self-energy diagram.', 'A.III.2 Contribution to the master equation': 'In order to see the effect of the renormalisation on the master equation, one has to evaluate the following expression, as discussed at the end of section IV B:  \nΞ R ( ω u , ⃗u, t 0 , t ) = ∫ t t 0 dτ ∫ R du 0 Π R vac ( u 2 ) cos[( u 0 -ω u )( t -τ )] . (A.III.21)  \nSubstituting t -τ → τ yields  \nΞ R ( ω u , ⃗u, t 0 , t ) = ∫ t -t 0 0 dτ ∫ R du 0 Π R vac ( u 2 ) cos[( u 0 -ω u ) τ ] (A.III.22)  \nand due to symmetry it holds that  \n∫ R du 0 Π R vac ( u 2 ) cos[( u 0 -ω u ) τ ] = ∫ R du 0 Π R vac ( ⃗u 2 -u 2 0 ) [cos( u 0 τ ) cos( ω u τ ) + sin( u 0 τ ) sin( ω u τ )] = cos( ω u τ ) ∫ R du 0 Π R vac ( ⃗u 2 -u 2 0 ) cos( u 0 τ ) . (A.III.23)  \nTo solve the integrations, we first consider all terms that depend on u 0 in the form ( u 2 + m 2 ) = ( ω 2 u -u 2 0 ) with ω u = √ ⃗u 2 + m 2 .  \nTo evaluate this, we would like to use the distributional integration ∫ R du 0 cos( u 0 τ ) = πδ ( τ ). However, in order for this to be true, we would need to have a Schwartz function paired with the distribution under the τ integration and its integration domain should be R . To have this, we modify the cosine slightly and we will see when evaluating the integration that this modification does not affect the final value. We introduce the Schwartz function S c ( ω u τ ) which coincides on the interval [ ϵ, t -t 0 -ϵ ] with cos( ω u τ ), where ϵ ≪ 1. On the interval [ -ϵ, ϵ ] it is a smooth function constructed in such a way that S c (0) = 1 2 cos(0) = 1 2 , d 2 dτ 2 S c ( ω u τ ) | τ =0 = 1 2 d 2 dτ 2 cos( ω u τ ) | τ =0 = -ω 2 u 2  \nand S c ( ω u τ ) = 0 for τ < ϵ . For [ t -t 0 -ϵ, t -t 0 + ϵ ] it is also constructed as a smooth function from S c ( ω u ( t -t 0 -ϵ )) = cos( t -t 0 -ϵ ) to 0 for S c ( ω u τ ) for τ > t -t 0 + ϵ . The construction is in such a way, that S c ( ω u τ ) is a Schwartz function on R . Then we have:  \n∫ t -t 0 0 dτ cos( ω u τ ) ∫ R du 0 ( ω 2 u -u 2 0 ) cos( u 0 τ ) = ∫ R dτ S c ( ω u τ ) [ πδ ( τ ) ω 2 u + d 2 dτ 2 ∫ R du 0 cos( u 0 τ ) ] = ∫ R dτ S c ( ω u τ ) [ πδ ( τ ) ω 2 u + π d 2 dτ 2 δ ( τ ) ] = ∫ R dτ [ πδ ( τ ) ω 2 u S c ( ω u τ ) + πδ ( τ ) d 2 dτ 2 S c ( ω u τ ) ] = π [ ω 2 u 2 -ω 2 u 2 ] = 0 . (A.III.24)  \nUsing this 16 , what remains is  \nπ 2 κ 2 ∫ t -t 0 0 dτ cos( ω u τ ) ∫ R du 0 cos( u 0 τ ) { ( m 4 -λ 4 ) m 2 + u 2 u 2 ln ( λ m ) + m 2 m 2 +2 u 2 + λ 2 u 2 √ ( m 2 + u 2 ) 2 -2 λ 2 ( m 2 -u 2 ) + λ 4 · · [ arctanh ( m 2 + u 2 -λ 2 √ ( m 2 + u 2 ) 2 -2 λ 2 ( m 2 -u 2 ) + λ 4 ) -arctanh ( m 2 -u 2 -λ 2 √ ( m 2 + u 2 ) 2 -2 λ 2 ( m 2 -u 2 ) + λ 4 )] + λ ( m 2 -λ 2 ) √ λ 2 -4 m 2 [ arctanh ( 2 m 2 -λ 2 λ √ λ 2 -4 m 2 ) +arctanh ( λ √ λ 2 -4 m 2 ) ]} , (A.III.25)  \nwhich is independent of µ . Before explicitly evaluating the integrations, we simplify the integrands by taking the limit λ → 0 where possible, as λ was only introduced as artificial small graviton mass to be able to fix the residuum of the pole in the propagator. As arctanh( i 0) and arctanh( i ∞ ) are finite, the last two lines vanish. Also, lim λ → 0 λ 4 ln( λ ) = 0 and, when expanding the square root in the second line, we find that lim λ → 0 λ arctanh(1 + λ ) = 0 as arctanh( x ) = 1 2 ln(1 + x ) -1 2 ln(1 -x ).  \nThis leaves us with  \nπ 2 κ 2 ∫ t -t 0 0 dτ cos( ω u τ ) ∫ R du 0 cos( u 0 τ ) { m 4 m 2 + u 2 u 2 ln ( λ m ) + m 2 m 2 +2 u 2 u 2 ( m 2 + u 2 ) · · [ arctanh ( m 2 + u 2 -λ 2 √ ( m 2 + u 2 ) 2 -2 λ 2 ( m 2 -u 2 ) + λ 4 ) -arctanh ( m 2 -u 2 m 2 + u 2 ) ]} . (A.III.26)  \nFrom the above named relation arctanh( x ) = 1 2 ln(1 + x ) -1 2 ln(1 -x ) follows that  \narctanh ( x y ) = 1 2 ln ( x + y y -x ) . (A.III.27)  \nThis yields for the last line:  \narctanh ( m 2 -u 2 m 2 + u 2 ) = 1 2 ln ( m 2 u 2 ) . (A.III.28)  \nFor the line before, we expand the argument in second order for small λ and obtain  \narctanh ( m 2 + u 2 -λ 2 √ ( m 2 + u 2 ) 2 -2 λ 2 ( m 2 -u 2 ) + λ 4 ) ≈ arctanh ( 1 -2 λ 2 u 2 ( m 2 + u 2 ) 2 ) . (A.III.29)  \nHigher orders will not contribute in the final limit λ → 0. Expressing the arctanh again in terms of logarithms, we obtain  \narctanh ( 1 -2 λ 2 u 2 ( m 2 + u 2 ) 2 ) = 1 2 ln ( 2 -2 λ 2 u 2 ( m 2 + u 2 ) 2 ) -1 2 ln ( 2 λ 2 u 2 ( m 2 + u 2 ) 2 ) ≈ 1 2 ln (2) -1 2 ln ( 2 λ 2 u 2 ( m 2 + u 2 ) 2 ) , (A.III.30)  \nwhere we neglected terms of order λ 2 and higher. With these simplifications one then finds that  \nΞ R ( ω u , ⃗u, t 0 , t ) = π 2 κ 2 ∫ t -t 0 0 dτ cos( ω u τ ) ∫ R du 0 cos( u 0 τ ) { m 4 m 2 + u 2 u 2 ln ( λ m ) -m 2 2 m 2 +2 u 2 u 2 ( m 2 + u 2 ) · · [ ln ( 2 λ 2 u 2 ( m 2 + u 2 ) 2 ) +ln ( m 2 2 u 2 ) ]} = π 2 κ 2 ∫ t -t 0 0 dτ cos( ω u τ ) ∫ R du 0 cos( u 0 τ ) { m 4 m 2 + u 2 u 2 ln ( λ m ) -m 2 2 m 2 +2 u 2 u 2 ( m 2 + u 2 ) · · ln ( λ 2 m 2 ( m 2 + u 2 ) 2 ) } . (A.III.31)  \nThe next step is to solve the following two integrations:  \n(A) m 4 ln ( λ m )∫ R du 0 cos( u 0 τ ) ω 2 u -u 2 0 ⃗u 2 -u 2 0 (A.III.32)  \n(B) -m 2 2 ∫ R du 0 cos( u 0 τ ) m 2 +2 ⃗u 2 -2 u 2 0 ⃗u 2 -u 2 0 ( ω 2 u -u 2 0 ) ln ( λ 2 m 2 ( ω 2 u -u 2 0 ) 2 ) . (A.III.33)  \nWe start with (A) :  \n∫ R du 0 cos( u 0 τ ) ω 2 u -u 2 0 ⃗u 2 -u 2 0 = -∫ R du 0 e iu 0 τ ω 2 u -u 2 0 ( u 0 -| ⃗u | )( u 0 + | ⃗u | ) . (A.III.34)  \nApplying the residue theorem by closing the contour with a semi-circle in the upper half plane leaves us with the contributions of the poles that lie on the contour, hence contribute + 1 2 their residua, yielding  \n-∫ R du 0 e iu 0 τ ω 2 u -u 2 0 ( u 0 -| ⃗u | )( u 0 + | ⃗u | ) = -πi [ e i | ⃗u | τ m 2 2 | ⃗u | -e -i | ⃗u | τ m 2 2 | ⃗u | ] = πm 2 | ⃗u | sin( | ⃗u | τ ) . (A.III.35)  \nTherefor we have  \n(A) = πm 6 | ⃗u | sin( | ⃗u | τ ) ln ( λ m ) . (A.III.36)  \nWe proceed with (B) :  \n∫ R du 0 cos( u 0 τ ) m 2 +2 ⃗u 2 -2 u 2 0 ⃗u 2 -u 2 0 ( ω 2 u -u 2 0 ) ln ( λ 2 m 2 ( ω 2 u -u 2 0 ) 2 ) = -∫ R du 0 e iu 0 τ m 2 +2 ⃗u 2 -2 u 2 0 ( u 0 -| ⃗u | )( u 0 + | ⃗u | ) ( ω 2 u -u 2 0 ) ln ( λ 2 m 2 ( ω 2 u -u 2 0 ) 2 ) . (A.III.37)  \nWe apply the residue theorem once again. The integrand has singularities at u 0 = ±| ⃗u | . Working with the principal value logarithm, i.e. the complex logarithm with branch cut at the negative real axis, the logarithm here does not have any branch cuts given that ω u > 0, hence the equation λ 2 m 2 ( ω 2 u -u 2 0 ) 2 ! = -α with α ∈ R , α > 0 does not have any solution for u 0 ∈ C . Due to the prefactor ( ω 2 u -u 2 0 ), there is no singularity in u 0 = ± ω u . We pick a closed integration contour from -∞ to -| ⃗u | -ϵ , then go in a semi-circle clockwise around the pole to -| ⃗u | + ϵ , continue to | ⃗u | -ϵ , again go around the singularity in a semi-circle clockwise to | ⃗u | + ϵ , continue to + ∞ and close it with a semicircle in the upper half-plane. The closed contour does not contain any singularities, hence its contribution vanishes. Due to the exponential e iu 0 τ , also the semi-circle at infinite radius in the upper half-plane vanishes. Hence only the contributions of the singularities at u 0 = ±| ⃗u | remain. As we went for the closed contour around them clockwise, the singularity contributions have to be evaluated counter-clockwise and added to the closed contour. We start by investigating the one at  \nu 0 = + | ⃗u | and replace u 0 = | ⃗u | + ϵe iϕ :  \n-lim ϵ → 0 ∫ π 0 dϕ iϵe iϕ e iτϵe iϕ + iτ | ⃗u | m 2 -4 | ⃗u | ϵe iϕ -2 ϵ 2 e 2 iϕ ϵe iϕ (2 | ⃗u | + ϵe iϕ ) ( m 2 -2 | ⃗u | ϵe iϕ -ϵ 2 e 2 iϕ ) · ln ( λ 2 m 2 ( m 2 -2 | ⃗u | ϵe iϕ -ϵ 2 e 2 iϕ ) 2 ) = -i lim ϵ → 0 ∫ π 0 dϕ e iτ | ⃗u | m 2 -4 | ⃗u | ϵe iϕ -2 ϵ 2 e 2 iϕ 2 | ⃗u | + ϵe iϕ ( m 2 -2 | ⃗u | ϵe iϕ -ϵ 2 e 2 iϕ ) ln ( λ 2 m 2 ( m 2 -2 | ⃗u | ϵe iϕ -ϵ 2 e 2 iϕ ) 2 ) = -i lim ϵ → 0 ∫ π 0 dϕ e iτ | ⃗u | m 2 1 2 | ⃗u | m 2 ln ( λ 2 m 2 ( m 2 -2 | ⃗u | ϵe iϕ -ϵ 2 e 2 iϕ ) 2 ) = -i m 4 2 | ⃗u | e iτ | ⃗u | lim ϵ → 0 ∫ π 0 dϕ ln ( λ 2 m 2 ) = -i πm 4 | ⃗u | e iτ | ⃗u | ln ( λ m ) , (A.III.38)  \nwhere in the first step we expanded e iτϵe iϕ + iτ | ⃗u | = e iτ | ⃗u | (1 + O ( ϵ )) and neglected all but the zeroth order due to the limit, in the second step we expanded 1 2 | ⃗u | + ϵe iϕ = 1 2 | ⃗u | ( 1 -ϵ 2 | ⃗u | e iϕ + O ( ϵ ) ) and applied the limit to the terms depending on ϵ where possible. In the third step, we expanded  \nln ( λ 2 m 2 ( m 2 -2 | ⃗u | ϵe iϕ -ϵ 2 e 2 iϕ ) 2 ) = ln ( λ 2 m 2 ) + ϵ 4 | ⃗u | m 2 e iϕ + O ( ϵ 2 ) . (A.III.39)  \nFor the other singularity at u 0 = -| ⃗u | we find analogously for u 0 = -| ⃗u | + ϵe iϕ :  \nlim ϵ → 0 ∫ π dϕ iϵe iϕ e iτϵe iϕ -iτ | ⃗u | m 2 +4 | ⃗u | ϵe iϕ -2 ϵ 2 e 2 iϕ ϵe iϕ ( -2 | ⃗u | + ϵe iϕ ) ( m 2 +2 | ⃗u | ϵe iϕ -ϵ 2 e 2 iϕ  \n-0 ) · ln ( λ 2 m 2 ( m 2 +2 | ⃗u | ϵe iϕ -ϵ 2 e 2 iϕ ) 2 ) = + i lim ϵ → 0 ∫ π 0 dϕ e -iτ | ⃗u | m 2 1 2 | ⃗u | m 2 ln ( λ 2 m 2 ( m 2 +2 | ⃗u | ϵe iϕ -ϵ 2 e 2 iϕ ) 2 ) = i m 4 2 | ⃗u | e -iτ | ⃗u | lim ϵ → 0 ∫ π 0 dϕ ln ( λ 2 m 2 ) = i πm 4 | ⃗u | e -iτ | ⃗u | ln ( λ m ) . (A.III.40)  \nCombining these two yields  \n(B) = -πm 6 | ⃗u | sin( | ⃗u | τ ) ln ( λ m ) = -(A) . (A.III.41)  \nFrom this follows that we have  \nΞ R ( ω u , ⃗u, t 0 , t ) = 0 . (A.III.42)  \nDue to the way the renormalised quantities entered in the master equation, the result is now independent of the scale µ as well as of the artificial graviton mass λ , whose limit to zero can therefore be taken without problems.', 'Appendix A.IV APPLICATION OF THE MARKOV APPROXIMATION': 'In this appendix, we apply the Markov approximation to the renormalised one-particle master equation (4.47). Applying the formula given in (5.2), we obtain two classes of terms, one class that contains the δ -distributions and the other one containing the Cauchy principal value. Before we evaluate them in subsections A.IV.2 and A.IV.3, we first discuss the applicability of the Markov approximation for ultra-relativistic particles with focus on neutrinos in subsection A.IV.1.', 'A.IV.1 Applicability of the Markov approximation for ultra-relativistic particles': 'In general, the Markov approximation can be applied if the timescales τ B on which the correlation functions decay are much smaller than the timescales τ R on which the state of the system varies (see [2]). The identification of these timescales is however hard without solving the one-particle master equation before the application of the approximation. As the Markov approximation corresponds to sending t 0 → -∞ and hence ∫ t -t 0 0 dτ -→ ∫ ∞ 0 dτ , we will analyse the error one makes when extending the integration domain from t -t 0 to ∞ . If the integrand is strongly peaked around τ = 0, which is usually assumed when deriving Markovian master equations, then the error of the additional contribution should be negligible. For this, we analyse the different parts of the renormalised one-particle master equation in (4.47):  \n∂ ∂t ρ ( ⃗u, ⃗v, t ) = -iρ ( ⃗u, ⃗v, t ) ( ω u -ω v ) -κ 2 ∫ d 3 k (2 π ) 3 { P u ( ⃗ k ) ω u -k ω u [ C R ( ⃗u, ⃗ k, t ) + δ P C R P ( ⃗u, ⃗ k, t ) ] + P v ( ⃗ k ) ω v -k ω v [( C R ( ⃗v, ⃗ k, t ) ) ∗ + δ P ( C R P ( ⃗v, ⃗ k, t ) ) ∗ ] } ρ ( ⃗u, ⃗v, t ) + κ 2 ∫ d 3 k (2 π ) 3 P ijln ( ⃗ k ) u i u j v l v n √ ω u + k ω u ω v + k ω v { C R ( ⃗u + ⃗ k, ⃗ k, t ) + ( C R ( ⃗v + ⃗ k, ⃗ k, t ) ) ∗ } ρ ( ⃗u + ⃗ k,⃗v + ⃗ k, t ) (A.IV.1)  \nwith  \nC R ( ⃗u, ⃗ k, t ) = 2 ∫ t -t 0 0 dτ Ω k N ( k ) cos[Ω k τ ] e -i ( ω u -k -ω u ) τ (A.IV.2)  \nC R P ( ⃗u, ⃗ k, t ) = 2 ∫ t -t 0 0 dτ Ω k N ( k ) cos[Ω k τ ] e -i ( ω u -k + ω u ) τ . (A.IV.3)  \nWe start with the real part of C R in the second line, which will later lead to decoherence. The term we have to take into account is  \n-κ 2 ∫ d 3 k (2 π ) 3 P u ( ⃗ k ) ω u -k ω u ∫ t -t 0 0 dτ Ω k N ( k ) [cos((Ω k + ω u -k -ω u ) τ ) + cos((Ω k -ω u -k + ω u ) τ )] = -κ 2(2 π ) 2 ∫ ∞ 0 dk ∫ π 0 dθ u 4 2 sin 5 ( θ ) ω u -k ω u kN ( k ) · [ sin(( k + ω u -k -ω u )( t -t 0 )) k + ω u -k -ω u + sin(( k -ω u -k + ω u )( t -t 0 )) k -ω u -k + ω u ] , (A.IV.4)  \nwhere in the second line we went to spherical coordinates with ⃗ k z || ⃗u and performed the τ -integration. The main contribution from the integrand will come from small k , because for k = 0 the denominator tends to zero and N ( k ) decreases rapidly for large k . For the Markov approximation however the behaviour depending on ( t -t 0 ) is important. To extract this, we first substitute  \nµ := ω u -k = √ ω 2 u + k 2 -2 uk cos( θ ) (A.IV.5)  \nand assume, already adapting the scenario of section VI D where we will apply the one-particle master equation to ultra-relativistic particles, that ω u ≈ u , which then yields  \nµ (0) = | k -u | µ ( π ) = k + u dθ = dµ µ uk sin( θ ) (A.IV.6)  \nas well as  \nsin 2 ( θ ) = 1 -( u 2 + k 2 -µ 2 ) 2 4 u 2 k 2 . (A.IV.7)  \nUsing this and defining ∆ t := t -t 0 , equation (A.IV.4) reads  \n-u 2 κ 4(2 π ) 2 ∫ ∞ 0 dk N ( k ) ∫ k + u | k -u | dµ [ 1 -( u 2 + k 2 -µ 2 ) 2 4 u 2 k 2 ] 2 [ sin(( k + µ -ω u )∆ t ) k + µ -ω u + sin(( k -µ + ω u )∆ t ) k -µ + ω u ] . (A.IV.8)  \nNow we make some assumptions on the involved quantities motivated by the application to ultrarelativistic neutrinos to simplify the integration. In order to continue, we assume that  \nu ≫ 1 cβ . (A.IV.9)  \nTo resolve the absolute value in the µ -integration, we split the k -integration into two regions, one with k < u and another one with k > u :  \n∫ u 0 dk N ( k ) ∫ u + k u -k + ∫ ∞ u dk N ( k ) ∫ u + k k -u . (A.IV.10)  \nThe second integral is now negligible compared to the first one, as on the one hand side the dominant contribution of the integrand is around k = 0 where the root of the denominator lies, and on the other hand due to condition (A.IV.9) for k > u ≫ 1 cβ also N ( k ) will strongly damp the integrand in that region. Hence we only continue with the first integral:  \n-u 2 κ 4(2 π ) 2 ∫ u 0 dk N ( k ) ∫ k + u u -k dµ [ 1 -( u 2 + k 2 -µ 2 ) 2 4 u 2 k 2 ] 2 [ sin(( k + µ -ω u )∆ t ) k + µ -ω u + sin(( k -µ + ω u )∆ t ) k -µ + ω u ] . (A.IV.11)  \nThe µ -integration can be solved and one obtains  \n-u 2 κ 4(2 π ) 2 ∫ u 0 dk N ( k ) k 4 (∆ t ) 8 u 4 [ 16 k ∆ t (45 + (∆ t ) 2 ( -6 k 2 +( -27 + 4 k 2 (∆ t ) 2 ) u 2 +(∆ t ) 2 u 4 )) +4 k ∆ t (135 + 2(∆ t ) 2 ( -3 k 2 +( -36 + k 2 (∆ t ) 2 ) u 2 +(∆ t ) 2 u 4 )) cos(2 k ∆ t ) -6(105 + 2(∆ t ) 2 ( -15 k 2 +( -30 + 7 k 2 (∆ t ) 2 ) u 2 +(∆ t ) 2 u 4 )) sin[2 k ∆ t ] ] . (A.IV.12)  \nFIG. 1. Decay of the integral I (∆ t ) in (A.IV.12) for different values of u and β (see main text) in percent normalised such that I (10 -11 s ) = 1.  \n<!-- image -->  \nThis can be integrated numerically. For 17 the three different values of u = 2 . 4 · 10 22 1 s , u = 2 . 4 · 10 24 1 s , u = 2 . 4 · 10 26 1 s and for the two values of β = 10 -11 s, β = 10 -13 s one obtains that the result vanishes rapidly for increasing times ∆ t , see figure 1.  \nIt hence is visible that the error made when not only considering time until t -t 0 but until ∞ is negligible, given that its value drops rapidly when increasing t -t 0 . In section VI D typical propagation times of neutrinos through the Earth are used, which are of the order t -t 0 ≈ 10 -2 s .  \nIn the general case, in which we do not specify to ultra-relativistic neutrinos, the integration boundaries for the µ -integration in (A.IV.8) would be µ (0) = √ m 2 +( k -u ) 2 as well as µ ( π ) = √ m 2 +( k + u ) 2 . As the temperature parameter Θ characterising the environment is finite, one can always find a parameter ˜ k such that ˜ k ≫ 1 cβ and hence the k -integration can be approximated ∫ ∞ 0 dk ≈ ∫ ˜ k 0 dk , since N ( k ) damps the result strongly on the remaining interval. The general integration of this quantity goes beyond the scope of this work. To facilitate its computation, one can use the specific properties of the considered model similarly to the way it is applied to ultra-relativistic neutrinos in this work.  \nNext, we consider the real terms arising from the extended projection in line two of equation (A.IV.1). The same steps for the ultra-relativistic conditions as above yield  \n-u 2 κ 4(2 π ) 2 ∫ u 0 dk N ( k ) ∫ k + u u -k dµ [ 1 -( u 2 + k 2 -µ 2 ) 2 4 u 2 k 2 ] 2 [ sin(( k + µ + ω u )∆ t ) k + µ + ω u + sin(( k -µ -ω u )∆ t ) k -µ -ω u ] . (A.IV.13)  \nNote that in contrast to the previous case, now the denominator is never zero, hence the contribution of this term is less dominant than the one of the previous terms. Solving the µ -integration and plotting the term with the same parameters as above again shows that the result vanishes rapidly for increasing values of t -t 0 . The absolute value of the integral is however several orders of magnitude below the contribution in (A.IV.12), hence it is negligible compared to the above term.  \nFor the imaginary terms in the second line of (A.IV.1), one can perform the same analysis and also obtains a strong decay in t -t 0 for the terms. The analogous arguments (if v ≫ 1 cβ ) hold for the  \nterms in the third line of (A.IV.1). For the last line, also the wave function ρ depends on ⃗ k , but as the structure of the involved terms is very similar and the same main arguments can be carried over (roots of the denominator, same main quantities u , v , β and damping due to N ( k )), also the behaviour of this term depending on t -t 0 is similar as above. Hence the Markov approximation is justified here under the above named assumptions of ultra-relativistic particles that fulfil (A.IV.9).', 'A.IV.2 Evaluation of the delta terms': "The main step in the application of the Markov approximation consists in replacing  \n∫ t -t 0 0 dτ e -iωτ -→ ∫ ∞ 0 dτ e -iωτ = πδ ( ω ) -PV ( i ω ) , (A.IV.14)  \nas argued in detail in the main text. In this subsection, we investigate the terms containing the delta distribution and in the next subsection we evaluate the principal value terms. The delta distribution parts of the terms in lines two and three of (A.IV.1) become  \n-πκ 2 ∫ d 3 k (2 π ) 3 { P u ( ⃗ k ) ω u -k ω u N ( k ) Ω k [ δ (Ω k -ω u -k + ω u ) + δ (Ω k + ω u -k -ω u ) + δ P δ (Ω k -ω u -k -ω u ) + δ P δ (Ω k + ω u -k + ω u ) ] + P v ( ⃗ k ) ω v -k ω v N ( k ) Ω k [ δ (Ω k -ω v -k + ω v ) + δ (Ω k + ω v -k -ω v ) + δ P δ (Ω k -ω v -k -ω v ) + δ P δ (Ω k + ω v -k + ω v ) ]} . (A.IV.15)  \nTo evaluate the delta distributions, we first have to determine the zeroes of the arguments. As the part for v is exactly the same as the one for u , we focus on the latter one. The four equations to solve then read, expressing the ⃗ k -integration in spherical coordinates ( k, θ, ϕ ) and picking them such that ⃗u || ⃗ k z :  \n± k = √ ω 2 u + k 2 -2 uk cos( θ ) -ω u (A.IV.16)  \n± k = √ ω 2 u + k 2 -2 uk cos( θ ) + ω u , (A.IV.17)  \nwhere u := | ⃗u | . With the limits set by the spherical coordinates, i.e. k ∈ { 0 , ∞} and cos( θ ) ∈ {-1 , 1 } it is evident that the right hand side of the equation in the second line is always positive, hence δ (Ω k + ω u -k + ω u ) = 0. For the positive sign in the second line we find  \n( k -ω u ) = √ ω 2 u + k 2 -2 uk cos( θ ) ⇐⇒ kω u = uk cos( θ ) , (A.IV.18)  \nwhich is solved for k = 0. There is no other solution, as the remaining equation ω u = √ m 2 + u 2 = u cos( θ ) is never fulfilled for m> 0. Note however that k = 0 is here only a solution of the squared equation and not of the original one 18 in (A.IV.17), hence we also have δ ( -Ω k + ω u -k + ω u ) = 0.  \nThis means that here all additional terms arising from the extended projection vanish. For the two equations in the first line we have  \n( ± k + ω u ) = √ ω 2 u + k 2 -2 uk cos( θ ) ⇐⇒ ∓ kω u = uk cos( θ ) . (A.IV.19)  \nThis is solved by k = 0, which is also a solution to both non-squared equations. Apart from that, there is no other solution as again | ∓ ω u | = | ∓ √ m 2 + u 2 | > u while | u cos( θ ) | ≤ u . From this we obtain  \nδ ( ± Ω k -ω u -k + ω u ) = δ ( k ) ∣ ∣ ∣ ± 1 + u cos( θ ) ω u ∣ ∣ ∣ = δ ( k ) 1 ± u cos( θ ) ω u . (A.IV.20)  \nApplying this to the original expression in spherical coordinates, we find  \n-πκ 2 1 (2 π ) 2 ∫ ∞ 0 dk ∫ π 0 dθ k sin( θ ) { P u ( ⃗ k ) ω u -k ω u N ( k ) [ δ ( k ) 1 + u cos( θ ) ω u + δ ( k ) 1 -u cos( θ ) ω u ] + P v ( ⃗ k ) ω v -k ω v N ( k ) [ δ ( k ) 1 + v cos( θ ) ω v + δ ( k ) 1 -v cos( θ ) ω v ]} , (A.IV.21)  \nwhere in the second line we picked different spherical coordinates with ⃗v || ⃗ k z and defined v := | ⃗v | . Evaluation of the delta yields, using  \nP u ( k, θ, ϕ ) = 1 2 [ u 2 -( uk cos( θ )) 2 k 2 ] 2 = u 4 2 [ 1 -cos 2 ( θ ) ] 2 = u 4 2 sin 4 ( θ ) = P u ( θ, ϕ ) (A.IV.22)  \nand with l'Hospital's limit  \nlim k → 0 k N ( k ) = lim k → 0 k e βk -1 = lim k → 0 1 βe βk = 1 β , (A.IV.23)  \nthe following result 19 , where we get an additional factor of 1 2 as the point k = 0, where the delta distribution does not vanish, lies at the edge of the integration area:  \n-πκ 8 β 1 (2 π ) 2 ∫ π 0 dθ sin 5 ( θ ) { u 4 ω 2 u [ 1 1 + u cos( θ ) ω u + 1 1 -u cos( θ ) ω u ] + v 4 ω 2 v [ 1 1 + v cos( θ ) ω v + 1 1 -v cos( θ ) ω v ]} = πκ 4 β 1 (2 π ) 2 ∫ π 0 dθ sin 5 ( θ ) { u 2 cos 2 ( θ ) -ω 2 u u 2 + v 2 cos 2 ( θ ) -ω 2 v v 2 } = πκ 4 β 1 (2 π ) 2 [ -10 3 ( u 2 + v 2 ) + 2( ω 2 u + ω 2 v ) -2 ω u u m 4 arccoth ( ω u u ) -2 ω v v m 4 arccoth ( ω v v ) ] . (A.IV.24)  \nNext, we compute the last line of equation (A.IV.1). The terms containing the delta distributions read  \nπκ 2 ∫ d 3 k (2 π ) 3 P ijln ( ⃗ k ) u i u j v l v n √ ω u + k ω u ω v + k ω v N ( k ) Ω k · { δ (Ω k + ω u -ω u + k ) + δ (Ω k -ω u + ω u + k ) + δ (Ω k + ω v -ω v + k ) + δ (Ω k -ω v + ω v + k ) } ρ ( ⃗u + ⃗ k,⃗v + ⃗ k, t ) . (A.IV.25)  \nProceeding the same way as above it follows that  \nδ ( ± Ω k + ω u -ω u + k ) = δ ( k ) ∣ ∣ ∣ 1 ∓ u cos( α uk ) ω u ∣ ∣ ∣ = ω u δ ( k ) ω u ∓ u cos( α uk ) , (A.IV.26)  \nwhere α uk is the angle between ⃗u and ⃗ k . Due to the presence of the two directions ⃗u and ⃗v in the prefactor, it is not possible to pick ⃗ k z parallel to both ⃗u and ⃗v , as in general ⃗u || ⃗v does not hold. Given that P ijln ( ⃗ k ) = P ijln ( θ, ϕ ), the delta distribution can be applied and one obtains 20  \nπκ 2(2 π ) 3 β ∫ π 0 dθ ∫ 2 π 0 dϕ sin( θ ) P ijln ( θ, ϕ ) u i u j v l v n ω u ω v · { ω 2 u ω 2 u -u 2 cos 2 ( α uk ) + ω 2 v ω 2 v -v 2 cos 2 ( α vk ) } ρ ( ⃗u, ⃗v, t ) . (A.IV.27)  \nNow we choose in the first term the spherical coordinates such that ⃗u || ⃗ k z and in the second one ⃗v || ⃗ k z and obtain  \nπκ 2(2 π ) 3 β ∫ π 0 dθ ∫ 2 π 0 dϕ sin( θ ) P ijln ( θ, ϕ ) u i u j v l v n ω u ω v { ω 2 u ω 2 u -u 2 cos 2 ( θ ) + ω 2 v ω 2 v -v 2 cos 2 ( θ ) } ρ ( ⃗u, ⃗v, t ) . (A.IV.28)  \nThe the only quantity depending on ϕ is P ijln ( θ, ϕ ). The appearing contraction is, assuming ⃗u || ⃗ k z :  \nP ijln ( θ, ϕ ) u i u j v l v n = [ ⃗u · ⃗v -( ⃗u · ⃗ k )( ⃗v · ⃗ k ) ⃗ k 2 ] 2 -1 2 [ ⃗u 2 -( ⃗u · ⃗ k ) 2 ⃗ k 2 ][ ⃗v 2 -( ⃗v · ⃗ k ) 2 ⃗ k 2 ] = [ ⃗u · ⃗v -u cos( θ )( ⃗v · ˆ ⃗ k ) ] 2 -1 2 [ u 2 -u 2 cos 2 ( θ ) ] [ v 2 -( ⃗v · ˆ ⃗ k ) 2 ] = ( ⃗u · ⃗v ) 2 + u 2 cos 2 ( θ )( ⃗v · ˆ ⃗ k ) 2 -2 u cos( θ )( ⃗u · ⃗v )( ⃗v · ˆ ⃗ k ) -1 2 u 2 sin 2 ( θ ) [ v 2 -( ⃗v · ˆ ⃗ k ) 2 ] , (A.IV.29)  \nwhere the unit vector ˆ ⃗ k is defined as  \nˆ ⃗ k =    sin( θ ) cos( ϕ ) sin( θ ) sin( ϕ ) cos( θ )    . (A.IV.30)  \nFrom this follows that  \n∫ 2 π 0 dϕ P ijln ( θ,ϕ ) u i u j v l v n =2 π ( ⃗u · ⃗v ) 2 + u 2 cos 2 ( θ )[2 π cos 2 ( θ ) v 2 z + π sin 2 ( θ )( v 2 x + v 2 y )] -4 πuv z cos 2 ( θ )( ⃗u · ⃗v ) -πu 2 2 sin 2 ( θ )[2 v 2 -2 v 2 3 cos 2 ( θ ) -sin 2 ( θ )( v 2 1 + v 2 2 )] = π [ 2( ⃗u · ⃗v ) 2 -u 2 2 ( v 2 + v 2 z ) ] + π cos 2 ( θ ) [ u 2 ( v 2 + v 2 z ) -4 uv z ( ⃗u · ⃗v ) ] + π cos 4 ( θ ) [ u 2 2 (3 v 2 z -v 2 ) ] = -π u 2 2 ( v 2 -3 v 2 z ) + π cos 2 ( θ ) u 2 ( v 2 -3 v 2 z ) -π cos 4 ( θ ) u 2 2 ( v 2 -3 v 2 z ) = -π 2 u 2 ( v 2 -3 v 2 z ) [ 1 -cos 2 ( θ ) ] 2 = -π 2 u 2 ( v 2 -3 v 2 z ) sin 4 ( θ ) , (A.IV.31)  \nwhere we used in the second step that we chose the coordinate system such that ⃗u = u⃗e z , hence ⃗u · ⃗v = uv z . The θ -integration then has the form  \n-π 2 u 2 ( v 2 -3 v 2 z ) ∫ π 0 dθ sin 5 ( θ ) ω 2 u ω 2 u -u 2 cos 2 ( θ ) = -πω u ( v 2 -3 v 2 z ) [ 5 3 ω u -ω 3 u u 2 + m 4 u 3 arctanh ( u ω u ) ] . (A.IV.32)  \nDue to symmetry, we get the same result for the other term just with the replacement u ↔ v , as for the corresponding terms we can analogously choose ⃗v = v⃗e z . The contribution of the last line of equation (A.IV.1) is therefore 21 :  \nκ 16 πβ { ω u ω v ( v 2 -3 ( ⃗u · ⃗v ) 2 u 2 )[ 5 3 -( 1 + m 2 u 2 ) + m 4 ω u u 3 arctanh ( u ω u )] + ω v ω u ( u 2 -3 ( ⃗u · ⃗v ) 2 v 2 )[ 5 3 -( 1 + m 2 v 2 ) + m 4 ω v v 3 arctanh ( v ω v )] } . (A.IV.33)  \nCollecting all contributions from the delta terms yields  \nκ 16 πβ { -10 3 ( u 2 + v 2 ) + 2( ω 2 u + ω 2 v ) -2 ω u u m 4 arctanh ( u ω u ) -2 ω v v m 4 arctanh ( v ω v ) -ω u ω v ( v 2 -3 ( ⃗u · ⃗v ) 2 u 2 )[ 2 3 -m 2 u 2 + m 4 ω u u 3 arctanh ( u ω u )] -ω v ω u ( u 2 -3 ( ⃗u · ⃗v ) 2 v 2 )[ 2 3 -m 2 v 2 + m 4 ω v v 3 arctanh ( v ω v )] } . (A.IV.34)  \nThis is the final form of the real part of the dissipator that causes decoherence. The rotating wave approximation which is carried out as a next step leaves this part of the master equation invariant.", 'A.IV.3 Evaluation of the Cauchy principal value contributions': 'For the Markov approximation, it remains to compute the terms that contain the Cauchy principal value in (A.IV.1) after the approximation (A.IV.14). The terms in line two of (A.IV.1) read  \niκ 2 ∫ d 3 k (2 π ) 3 u 4 2 sin 4 ( α uk ) 1 ω u -k ω u N ( k ) Ω k [ PV ( 1 Ω k + ω u -k -ω u ) -PV ( 1 Ω k -ω u -k + ω u ) + δ P PV ( 1 Ω k + ω u -k + ω u ) -δ P PV ( 1 Ω k -ω u -k -ω u ) ] , (A.IV.35)  \nwhere α uk is the angle between ⃗u and ⃗ k . It can be seen that the term inside the principal value causes problems for ⃗ k → ⃗ 0. Thus we exclude a small region of radius ϵ around ⃗ k = ⃗ 0, perform the integration and take the limit ϵ → 0 in the end. We then obtain in spherical coordinates 22 with ⃗ k z || ⃗u :  \nlim ϵ → 0 iκu 4 4(2 π ) 2 ω u ∫ ∞ ϵ dk ∫ π 0 dθ sin 5 ( θ ) 1 ω u -k kN ( k ) [ 1 Ω k + ω u -k -ω u -1 Ω k -ω u -k + ω u + δ P 1 Ω k + ω u -k + ω u -δ P 1 Ω k -ω u -k -ω u ] . (A.IV.36)  \nFor δ P = 1 this can be simplified to  \nlim ϵ → 0 iκu 4 4(2 π ) 2 ω u ∫ ∞ ϵ dk ∫ π 0 dθ sin 5 ( θ ) 1 ω u -k kN ( k ) [ 2 u cos( θ ) ω u -k kω 2 u -ku 2 cos 2 ( θ ) ] = lim ϵ → 0 iκu 5 2(2 π ) 2 ω u [ ∫ ∞ ϵ dk N ( k ) ][ ∫ π 0 dθ cos( θ ) sin 5 ( θ ) ω 2 u -u 2 cos 2 ( θ ) ] ︸ ︷︷ ︸ =0 = lim ϵ → 0 0 = 0 . (A.IV.37)  \nWithout the principal value the k -integration would diverge:  \n∫ ∞ 0 dk N ( k ) = ∫ ∞ 0 dk 1 e βk -1 = iπ β -lim ϵ → 0 2 β arctanh ( 1 -2 e ϵβ ) . (A.IV.38)  \nWithout prior renormalisation, some terms arising due to the additional term present in the nonrenormalised coefficients C ( ⃗u, ⃗ k, t ) in (4.45) compared to the renormalised one would remain here and lead to logarithmic divergences, as expected from the discussion in section IV A.  \nFor the non-extended projection, i.e. for δ P = 0, the situation is more complicated. In that case we find, again using spherical coordinates and implementing the principal value by excluding a sphere of radius ϵ around the critical point ⃗ k = ⃗ 0:  \nlim ϵ → 0 iκu 4 4(2 π ) 2 ω u ∫ ∞ ϵ dk ∫ π 0 dθ sin 5 ( θ ) 1 ω u -k kN ( k ) [ 1 Ω k + ω u -k -ω u -1 Ω k -ω u -k + ω u ] = -lim ϵ → 0 iκu 4 2(2 π ) 2 ω u ∫ ∞ ϵ dk ∫ π 0 dθ sin 5 ( θ ) kN ( k ) 1 -ω u ω u -k k 2 -( ω u -k -ω u ) 2 . (A.IV.39)  \nThe θ -integration leads to a complicated result that can be simplified when considering e.g. the ultra-relativistic limit. Then it yields, where the limit ϵ → 0 can be taken also before the integration:  \n-iκu 4 105(2 π ) 2 ω u { -4 ∫ u 0 dk N ( k ) ( k 3 u 4 -7 k u 2 ) + ∫ ∞ u dk N ( k ) ( 35 1 u -14 u k 2 +3 u 3 k 4 ) } = -iκu 4 105(2 π ) 2 ω u { 4 [ π 4 15 β 4 u 4 -7 π 2 6 β 2 u 2 -6 ln(1 -e -βu ) βu +4 Li 2 ( e -βu ) β 2 u 2 -6 Li 3 ( e -βu ) β 3 u 3 -6 Li 4 ( e -βu ) β 4 u 4 ] +35 -35 ln ( e βu -1 ) βu -14 u ∫ ∞ u dk N ( k ) k 2 +3 u 3 ∫ ∞ u dk N ( k ) k 4 } . (A.IV.40)  \nHere, Li s ( x ) denotes the poly-logarithm function defined by  \nLi s ( x ) := ∞ ∑ n =1 x n n s . (A.IV.41)  \nThe remaining two integrations cannot be performed analytically, but they can be solved numerically given a specific temperature Θ and a value for u , and are finite as long as u > 0. For the terms in line three of (A.IV.1) we get the same results when replacing ⃗u → ⃗v and applying complex conjugation. For the terms in line four we obtain:  \niκ 2(2 π ) 3 1 √ ω u ω v ∫ d 3 k P ijln ( ⃗ k ) u i u j v l v n √ ω u + k ω v + k N ( k ) Ω k ρ ( ⃗u + ⃗ k,⃗v + ⃗ k, t ) { PV ( 1 Ω k + ω u -ω u + k ) -PV ( 1 Ω k -ω u + ω u + k ) -PV ( 1 Ω k + ω v -ω v + k ) +PV ( 1 Ω k -ω v + ω v + k ) } . (A.IV.42)  \nWithout further specification of ρ , this cannot be simplified further at this point.  \nIn summary, after the second Markov approximation we hence have  \n∂ ∂t ρ ( ⃗u, ⃗v, t ) = )  \n-iρ ( ⃗u, ⃗v, t ) ( ω u -ω v ) + κ 16 πβ { -10 3 ( u 2 + v 2 ) + 2( ω 2 u + ω 2 v ) -2 ω u u m 4 arctanh ( u ω u ) -2 ω v v m 4 arctanh ( v ω v -ω u ω v ( v 2 -3 ( ⃗u · ⃗v ) 2 u 2 )[ 2 3 -m 2 u 2 + m 4 ω u u 3 arctanh ( u ω u )] -ω v ω u ( u 2 -3 ( ⃗u · ⃗v ) 2 v 2 )[ 2 3 -m 2 v 2 + m 4 ω v v 3 arctanh ( v ω v )] } ρ ( ⃗u, ⃗v, t ) + iκ 2(2 π ) 3 1 √ ω u ω v ∫ d 3 k P ijln ( ⃗ k ) u i u j v l v n √ ω u + k ω v + k N ( k ) Ω k ρ ( ⃗u + ⃗ k,⃗v + ⃗ k, t ) { PV ( 1 Ω k + ω u -ω u + k ) -PV ( 1 Ω k -ω u + ω u + k ) -PV ( 1 Ω k + ω v -ω v + k ) +PV ( 1 Ω k -ω v + ω v + k ) } -(1 -δ P ) iκ 2(2 π ) 2 lim ϵ → 0 [ u 4 ω u ∫ ∞ ϵ dk ∫ π 0 dθ sin 5 ( θ ) kN ( k ) 1 -ω u ω u -k k 2 -( ω u -k -ω u ) 2 -v 4 ω v ∫ ∞ ϵ dk ∫ π 0 dθ sin 5 ( θ ) kN ( k ) 1 -ω v ω v -k k 2 -( ω v -k -ω v ) 2 ] ρ ( ⃗u, ⃗v, t ) , (A.IV.43)  \nwhere the last two lines vanish when working with the extended projection.', 'Appendix A.V APPLICATION OF THE ROTATING WAVE APPROXIMATION (RWA)': "In this appendix the detailed implementation of the rotating wave approximation for the renormalised Markovian master equation under consideration in this work is discussed. We proceed in the standard way by considering the master equation in interaction picture and removing all terms that oscillate fast (see e.g. [2, 12]. For this we start in the field theory and consider the full dissipator from (4.60) in [1]:  \nD [ ρ S ] = -κ 2 ∑ r ∈{ + , -} 4 ∑ a,b =1 ∫ R 3 d 3 k d 3 p d 3 l (2 π ) 3 1 Ω k { [ j a r ( ⃗ k, ⃗p ) † , j b r ( ⃗ k, ⃗ l ) ρ S ( t ) ] f (Ω k + ω b ( ⃗ k, ⃗ l )) + h.c. + N ( k ) [ j a r ( ⃗ k, ⃗p ) † , [ j b r ( ⃗ k, ⃗ l ) , ρ S ( t ) ]] f (Ω k + ω b ( ⃗ k, ⃗ l )) + h.c. } ,  \n(A.V.1)  \nwhere f ( ω ; t ) := ∫ t -t 0 0 dτ e -iωτ and we assume that the Markov approximation has already been applied, thus the former functions f ( ω ; t ) are now distributions f ( ω ) = πδ ( ω ) -iPV ( 1 ω ) independent of time. The different ω a and j r a were defined starting in (A.I.8). Taking into account that the renormalisation removed the terms that are independent of N ( k ), the dissipator becomes (see  \n(4.48)):  \nD [ ρ S ] = -κ 2 ∑ r ∈{ + , -} 4 ∑ a,b =1 ∫ R 3 d 3 k d 3 p d 3 l (2 π ) 3 N ( k ) Ω k { [ j a r ( ⃗ k, ⃗p ) † , [ j b r ( ⃗ k, ⃗ l ) , ρ S ( t ) ]] f (Ω k + ω b ( ⃗ k, ⃗ l )) + h.c. } . (A.V.2)  \nThe basis for the RWA is the dissipator in interaction picture, which is obtained by substituting  \nj a r ( ⃗ k, ⃗p ) -→ j a r ( ⃗ k, ⃗p ) e iω a ( ⃗ k,⃗p ) t . (A.V.3)  \nThus we get as time-dependent frequencies that cause oscillations terms of the form  \ne ± i [ ω a ( ⃗ k,⃗p ) -ω b ( ⃗ k, ⃗ l )] t . (A.V.4)  \nNext we apply the RWA, which means that we discard all the rapidly oscillating terms and only keep those where ω a ( ⃗ k, ⃗p ) = ω b ( ⃗ k, ⃗ l ) holds. For a = 3 and a = 4 this means that only b = 3 and b = 4 survive. However, from the definition of the ω a ( ⃗ k, ⃗p ) follows that ω 1 ( ⃗ k, ⃗p ) = ω 2 ( ⃗ k, -⃗ k -⃗ p ), so also terms of the form a = 1, b = 2 and vice versa will remain. To simplify this, we introduce J A r ( ⃗ k, ⃗p ) and ω A ( ⃗ k, ⃗p ) with A ∈ { 2 , 3 , 4 } , similar as in [12], such that  \nJ 2 r ( ⃗ k, ⃗p ) := 2 j 1 r ( ⃗ k, ⃗p ) ω 2 ( ⃗ k, ⃗p ) := ω 1 ( ⃗ k, ⃗p ) = ω 2 ( ⃗ k, -⃗ k -⃗ p ) = ω p -ω k + p (A.V.5)  \nJ 3 r ( ⃗ k, ⃗p ) := j 3 r ( ⃗ k, ⃗p ) ω 3 ( ⃗ k, ⃗p ) := ω 3 ( ⃗ k, ⃗p ) = -ω p -ω k + p (A.V.6)  \nJ 4 r ( ⃗ k, ⃗p ) := j 4 r ( ⃗ k, ⃗p ) ω 4 ( ⃗ k, ⃗p ) := ω 4 ( ⃗ k, ⃗p ) = ω p + ω k + p . (A.V.7)  \nMaking use of the fact that j 1 r ( ⃗ k, ⃗p ) = j 2 r ( ⃗ k, -⃗ k -⃗ p ), we can rewrite the dissipator in terms of a sum over capital letters A,B :  \nD [ ρ S ] = -κ 2 ∑ r ∈{ + , -} 4 ∑ A,B =2 ∫ R 3 d 3 k d 3 p d 3 l (2 π ) 3 N ( k ) Ω k { [ J A r ( ⃗ k, ⃗p ) † , [ J B r ( ⃗ k, ⃗ l ) , ρ S ( t ) ]] f (Ω k + ω B ( ⃗ k, ⃗ l )) + h.c. } . (A.V.8)  \nThe RWA-requirement to keep only the terms where ω A ( ⃗ k, ⃗p ) = ω B ( ⃗ k, ⃗ l ) then keeps the following summands  \nA = B = 2 A = B = 3 A = B = 4 (A.V.9)  \ncompletely, which were also exactly the same summands that survived the one-particle projection, while for the remaining six summands it yields the following conditions:  \nA = 2 , B = 3 : ω p -ω k + p + ω l + ω k + l = 0 (A.V.10)  \nA = 2 , B = 4 : ω p -ω k + p -ω l -ω k + l = 0 (A.V.11)  \nA = 3 , B = 4 : -ω p -ω k + p -ω l -ω k + l = 0 , (A.V.12)  \nand for the other three summands the same conditions with the role of ⃗ p and ⃗ l swapped. For mass m > 0, which implies ω k > 0, the last condition is never fulfilled. For m = 0 there is a solution, namely ⃗ k = ⃗ p = ⃗ l = 0. The first condition reads  \nω k + p = ω p + ω l + ω k + l . (A.V.13)  \nThis is never fulfilled, as we show in the following (where we define p := || ⃗ p || , l := || ⃗ l || , k := || ⃗ k || and assume 23 m> 0):  \nω k + p ≤ √ ( k + p ) 2 + m 2 < √ p 2 + m 2 + √ l 2 + m 2 + √ ( k -l ) 2 + m 2 ≤ ω p + ω l + ω k + l . (A.V.14)  \nThe inequality in the middle can be proven by considering the square of both sides (as each summand individually is positive, the direction of the inequality remains unaffected), which yields:  \nkp + kl < m 2 + l 2 + √ p 2 + m 2 √ ( k -l ) 2 + m 2 + √ p 2 + m 2 √ l 2 + m 2 + √ l 2 + m 2 √ ( k -l ) 2 + m 2 . (A.V.15)  \nWe can now estimate the right hand side downwards as  \nm 2 + l 2 + p | k -l | + pl + l | k -l | < m 2 + l 2 + √ p 2 + m 2 √ ( k -l ) 2 + m 2 + √ p 2 + m 2 √ l 2 + m 2 + √ l 2 + m 2 √ ( k -l ) 2 + m 2 , (A.V.16)  \nwhich yields  \nkp + kl < m 2 + l 2 + p | k -l | + pl + l | k -l | (A.V.17)  \nto be proven. Due to the absolute value, we consider two different cases: Let's first assume that k > l . Then the inequality reads  \nkp + kl < m 2 + kp + kl , (A.V.18)  \nwhich is true as long as m> 0. In the second case, i.e. for l ≥ k , we are left with  \nkp + kl < m 2 +2 l 2 +2 pl -pk -kl , (A.V.19)  \nor equivalently  \n0 < m 2 +2[ l 2 -kl ] + 2[ pl -pk ] . (A.V.20)  \nHowever, as we are considering the case l ≥ k , both brackets yield non-negative results and thus the inequality is fulfilled for m > 0. Hence (A.V.10) does also not have any solutions as long as m> 0. For m = 0 we have to solve the following equality:  \n| ⃗ k + ⃗ p | -p ! = | ⃗ k + ⃗ l | -l . (A.V.21)  \nHowever, we also have  \n| ⃗ k + ⃗ p | -p ≤ k ≤ | k -l | + l ≤ | ⃗ k + ⃗ l | -⃗ l , (A.V.22)  \nhence in order for equation (A.V.21) to hold, all ≤ signs must become equalities. For the first one this is the case if ⃗ k ∥ ⃗ p and for the last one if ⃗ k ∥ -⃗ l . For the one in the middle, we have to consider two cases on how to resolve | k -l | . If k ≥ l , then we can directly drop the absolute value and the middle ≤ becomes an equality. In case k < l we find | k -l | + l = 2 l -k > 2 k -k = k , so there is no equality. Hence for m = 0 the following solutions exist:  \n⃗ k ∥ ⃗ p ∧ ⃗ k ∥ -⃗ l ∧ k ≥ l . (A.V.23)  \nIt remains to investigate (A.V.11). Isolating ω p on one side and following the same argumentation as above (squaring the inequality and estimating downwards the rights hand side) we end up with  \nkp + kl < k 2 + m 2 + l 2 + l | k -p | + | k -p || k -l | + l | k -l | . (A.V.24)  \nHere we have to consider four different cases:  \n- · k ≥ p and k ≥ l : We obtain 2 kp < 2 k 2 + m 2 , which has no solution for m> 0 and for m = 0 and k = p we get solutions, thus k ≥ l, k = p, ⃗ k ∥ -⃗ l, ⃗ k ∥ -⃗ p .\n- · k ≥ p and k < l : We obtain 0 < m 2 +2 l 2 -2 pl which has no solution for m ≥ 0.\n- · k < p and k ≥ l : We obtain 0 < m 2 , which has only for m = 0 a solution, thus there k ≥ l, k < p, ⃗ k ∥ -⃗ l, ⃗ k ∥ -⃗ p .\n- · k < p and k < l : We obtain 0 < m 2 +2[( p -k )( l -k ) + l ( l -k )] and thus no solution as every bracket is positive.  \nSummarising, for equality (A.V.11) we get again no solution if m> 0 and for m = 0 we have  \n⃗ k ∥ -⃗ p ∧ ⃗ k ∥ -⃗ l ∧ p ≥ k ≥ l . (A.V.25)  \nSo in total, for a positive mass non of the non-diagonal terms survives the rotating wave approximation. If the mass is zero, the following non-diagonal terms survive:  \nA = 2 , B = 3 : ⃗ k ∥ ⃗ p ∧ ⃗ k ∥ -⃗ l ∧ k ≥ l (A.V.26) A = 2 , B = 4 : ⃗ k ∥ -⃗ p ∧ ⃗ k ∥ -⃗ l ∧ p ≥ k ≥ l (A.V.27) A = 3 , B = 4 : ⃗ k = ⃗ l = ⃗ p = 0 (A.V.28) A = 3 , B = 2 : ⃗ k ∥ ⃗ l ∧ ⃗ k ∥ -⃗ p ∧ k ≥ p (A.V.29) A = 4 , B = 2 : ⃗ k ∥ -⃗ l ∧ ⃗ k ∥ -⃗ p ∧ l ≥ k ≥ p (A.V.30) A = 4 , B = 3 : ⃗ k = ⃗ l = ⃗ p = 0 . (A.V.31)  \nIf we plug these special cases into the dissipator, which contains a projection of ⃗ l and also of ⃗ p onto the plane perpendicular to ⃗ k , all extra terms containing ⃗ k ∥ ± ⃗ l and ⃗ k ∥ ± ⃗ p vanish. Then only the special solution ⃗ k = ⃗ l = ⃗ p = 0 remains in all six cases (due to polar coordinates and thus also spherical coordinates being non-unique for zero radius, in that case still all directions are possible). However, as this is only one point regarding the radius integrations, it will vanish under the integral. Thus all the extra correction terms vanish 24 and we are left with the dissipator after the RWA in the form  \nD [ ρ S ] = -κ 2 ∑ r ∈{ + , -} 4 ∑ A =2 ∫ R 3 d 3 k d 3 p d 3 l (2 π ) 3 N ( k ) Ω k · { [ J A r ( ⃗ k, ⃗p ) † , [ J A r ( ⃗ k, ⃗ l ) , ρ S ( t ) ]] f (Ω k + ω A ( ⃗ k, ⃗ l )) + h.c. }∣ ∣ ∣ ∣ ∣ ω A ( ⃗ k,⃗p )= ω A ( ⃗ k, ⃗ l ) (A.V.32)  \n.  \nIn order to explicitly write the hermitian conjugate, we split f ( ω ) = f δ ( ω ) + f PV ( ω ), as the two parts behave differently under complex conjugation (the first part is real, the second one purely imaginary) and compute them in the following two subsections.", 'A.V.1 Computation of the delta terms in the RWA': 'The delta terms remains unaffected by the complex conjugation 25 , thus we get, using that the terms in the first line, which are independent of N ( k ), vanished due to the Markov approximation:  \nD δ [ ρ S ] = κ ∑ r ∈{ + , -} 4 ∑ A =2 ∫ R 3 d 3 k d 3 p d 3 l (2 π ) 3 f δ (Ω k + ω A ( ⃗ k, ⃗ l )) Ω k N ( k ) · { ( J A r ( ⃗ k, ⃗p ) † ρJ A r ( ⃗ k, ⃗ l ) -1 2 { ρ, J A r ( ⃗ k, ⃗ l ) J A r ( ⃗ k, ⃗p ) † } ) + ( J A r ( ⃗ k, ⃗ l ) ρJ A r ( ⃗ k, ⃗p ) † -1 2 { ρ, J A r ( ⃗ k, ⃗p ) † J A r ( ⃗ k, ⃗ l ) } ) }∣ ∣ ∣ ∣ ∣ ω A ( ⃗ k,⃗p )= ω A ( ⃗ k, ⃗ l ) . (A.V.33)  \nUsing the following equalities:  \nJ 2 r ( ⃗ k, ⃗p ) † = J 2 r ( -⃗ k, ⃗ k + ⃗ p ) ω 2 ( -⃗ k, ⃗p + ⃗ k ) = -ω 2 ( ⃗ k, ⃗p ) (A.V.34)  \nJ 3 r ( ⃗ k, ⃗p ) † = J 4 r ( -⃗ k, -⃗ p ) ω 3 ( -⃗ k, -⃗ p ) = -ω 4 ( ⃗ k, ⃗p ) (A.V.35)  \nJ 4 r ( ⃗ k, ⃗p ) † = J 3 r ( -⃗ k, -⃗ p ) ω 4 ( -⃗ k, -⃗ p ) = -ω 3 ( ⃗ k, ⃗p ) , (A.V.36)  \nwe can simplify the δ -part of the dissipator and obtain  \nD δ [ ρ S ] = κ ∑ r ∈{ + , -} 4 ∑ A =2 ∫ R 3 d 3 k d 3 p d 3 l (2 π ) 3 f δ (Ω k + ω A ( ⃗ k, ⃗ l )) + f δ (Ω k -ω A ( ⃗ k, ⃗ l )) Ω k N ( k ) · { ( J A r ( ⃗ k, ⃗ l ) ρJ A r ( ⃗ k, ⃗p ) † -1 2 { ρ, J A r ( ⃗ k, ⃗p ) † J A r ( ⃗ k, ⃗ l ) } ) }∣ ∣ ∣ ∣ ∣ ω A ( ⃗ k,⃗p )= ω A ( ⃗ k, ⃗ l ) . (A.V.37)  \nAdditionally, all terms involving the extended projection, i.e. all terms where A = 3 or A = 4, vanished due to the Markov approximation, hence the δ -part of the dissipator simplifies even further:  \nD δ [ ρ S ] = κ ∑ r ∈{ + , -} ∫ R 3 d 3 k d 3 p d 3 l (2 π ) 3 f δ (Ω k + ω 2 ( ⃗ k, ⃗ l )) + f δ (Ω k -ω 2 ( ⃗ k, ⃗ l )) Ω k N ( k ) · { ( J 2 r ( ⃗ k, ⃗ l ) ρJ 2 r ( ⃗ k, ⃗p ) † -1 2 { ρ, J 2 r ( ⃗ k, ⃗p ) † J 2 r ( ⃗ k, ⃗ l ) } ) }∣ ∣ ∣ ∣ ∣ ω 2 ( ⃗ k,⃗p )= ω 2 ( ⃗ k, ⃗ l ) (A.V.38)  \nAs shown above, f δ ∝ δ ( k ), hence the RWA condition reads  \nω 2 ( ⃗ k, ⃗p ) = ω p -ω k + p = 0 = ω l -ω k + l = ω 2 ( ⃗ k, ⃗ l ) (A.V.39)  \n.  \nand is therefore automatically fulfilled, thus it can be dropped. One can furthermore evaluate the f δ and obtains, using (A.IV.20):  \nf δ (Ω k + ω 2 ( ⃗ k, ⃗ l )) + f δ (Ω k -ω 2 ( ⃗ k, ⃗ l )) = πδ ( k ) [ 1 1 + l cos( θ l ) ω l + 1 1 -l cos( θ l ) ω l ] = δ ( k ) 2 π 1 -l 2 cos 2 ( θ l ) ω 2 l , (A.V.40)  \nwhere θ l denotes the angle between ⃗ k and ⃗ l and l = | ⃗ l | , k = | ⃗ k | . This yields  \nD δ [ ρ S ] = κ ∑ r ∈{ + , -} ∫ R 3 d 3 k d 3 p d 3 l (2 π ) 3 δ ( k ) Ω k 2 π 1 -l 2 cos 2 ( θ l ) ω 2 l N ( k ) · { ( J 2 r ( ⃗ k, ⃗ l ) ρJ 2 r ( ⃗ k, ⃗p ) † -1 2 { ρ, J 2 r ( ⃗ k, ⃗p ) † J 2 r ( ⃗ k, ⃗ l ) } ) } . (A.V.41)  \nThe rotating wave condition further implies  \n∂ ∂ | k | ω 2 ( ⃗ k, ⃗p ) = ∂ ∂ | k | ω 2 ( ⃗ k, ⃗ l ) ⇐⇒ k + p cos( θ p ) ω k + p = k + l cos( θ l ) ω k + l . (A.V.42)  \nWe therefore can use that  \nδ ( k ) 2 π 1 -l 2 cos 2 ( θ l ) ω 2 l = 2 πδ ( k ) 1 1 -( k + l cos( θ l )) 2 ω 2 k + l = 2 πδ ( k ) 1 √ 1 -( k + l cos( θ l )) 2 ω 2 k + l 1 √ 1 -( k + p cos( θ p )) 2 ω 2 k + p (A.V.43)  \nand, defining the Lindblad operators  \nL r ( ⃗ k ) := ∫ R 3 d 3 p 1 √ 1 -( k + p cos( θ p )) 2 ω 2 k + p J 2 r ( ⃗ k, ⃗p ) , (A.V.44)  \nwe can recast the dissipator in Lindblad form:  \nD δ [ ρ S ] = κ ∑ r ∈{ + , -} ∫ R 3 d 3 k (2 π ) 2 δ ( k ) N ( k ) Ω k ( L r ( ⃗ k ) ρL r ( ⃗ k ) † -1 2 { ρ, L r ( ⃗ k ) † L r ( ⃗ k ) } ) . (A.V.45)  \nAs the rotating wave approximation dropped the same terms as the single-particle projection and led to a condition on the frequencies that is already implemented in δ ( k ), which is present in every term of the δ -part of the dissipator after the Markov approximation, the RWA does not change the form of the δ -part of the dissipator compared to its form after the Markov approximation in (A.IV.34).', 'A.V.2 Computation of the Cauchy principal value terms in the RWA': "The contributions involving the Cauchy principal value, denoted as the PV-part of the dissipator (A.V.32) in the following reads after renormalisation, which removes the terms independent of  \nN ( k ):  \nD PV [ ρ S ]  \n= -κ 2 ∑ r ∈{ + , -} 4 ∑ A =2 ∫ R 3 d 3 k d 3 p d 3 l (2 π ) 3 N ( k ) Ω k · { [ J A r ( ⃗ k, ⃗p ) † , [ J A r ( ⃗ k, ⃗ l ) , ρ S ( t ) ]] f PV (Ω k + ω A ( ⃗ k, ⃗ l )) + h.c. }∣ ∣ ∣ ∣ ∣ ω A ( ⃗ k,⃗p )= ω A ( ⃗ k, ⃗ l ) . (A.V.46)  \nAs f ( ω ) is purely imaginary, it switches sign under the hermitian conjugation and we obtain:  \n.  \nD PV [ ρ S ] = -κ 2 ∑ r ∈{ + , -} 4 ∑ A =2 ∫ R 3 d 3 k d 3 p d 3 l (2 π ) 3 N ( k ) Ω k · { [[ J A r ( ⃗ k, ⃗p ) † , J A r ( ⃗ k, ⃗ l ) ] , ρ S ( t ) ] f PV (Ω k + ω A ( ⃗ k, ⃗ l )) }∣ ∣ ∣ ∣ ∣ ω A ( ⃗ k,⃗p )= ω A ( ⃗ k, ⃗ l ) (A.V.47)  \nWe can then rewrite this part of the dissipator as  \nD PV [ ρ S ] = -i κ 2 [ V LS , ρ ] (A.V.48)  \nwith  \nV LS = -∑ r ∈{ + , -} 4 ∑ A =2 ∫ R 3 d 3 k d 3 p d 3 l (2 π ) 3 N ( k ) Ω k PV ( 1 Ω k + ω A ( ⃗ k, ⃗ l ) ) [ J A r ( ⃗ k, ⃗p ) † , J A r ( ⃗ k, ⃗ l ) ] ∣ ∣ ∣ ∣ ∣ ω A ( ⃗ k,⃗p )= ω A ( ⃗ k, ⃗ l ) . (A.V.49)  \nNote that the rotating wave approximation hence removed the imaginary terms in the fifth to seventh line of the Markovian master equation in (5.3), while it did not change the other imaginary terms. These were vanishing when working with the extended projection, hence in that case we find V LS = 0. If the non-extended projection is used, only the terms for A = 2 are left and there the RWA condition is already implemented in the one-particle projection, as the case (1 , 1) includes δ ( ⃗ p -⃗ l ), see table (I) and (II) at the beginning of section A.I.2. Hence the RWA does not change anything in the remaining PV-terms. The final one-particle master equation then becomes  \n∂ ∂t ρ ( ⃗u, ⃗v, t ) =  \n-iρ ( ⃗u, ⃗v, t ) ( ω u -ω v ) + κ 16 πβ { -10 3 ( u 2 + v 2 ) + 2( ω 2 u + ω 2 v ) -2 ω u u m 4 arctanh ( u ω u ) -2 ω v v m 4 arctanh ( v ω v ) -ω u ω v ( v 2 -3 ( ⃗u · ⃗v ) 2 u 2 )[ 2 3 -m 2 u 2 + m 4 ω u u 3 arctanh ( u ω u )] -ω v ω u ( u 2 -3 ( ⃗u · ⃗v ) 2 v 2 )[ 2 3 -m 2 v 2 + m 4 ω v v 3 arctanh ( v ω v )] } ρ ( ⃗u, ⃗v, t ) -(1 -δ P ) iκ 2(2 π ) 2 lim ϵ → 0 [ u 4 ω u ∫ ∞ ϵ dk ∫ π 0 dθ sin 5 ( θ ) kN ( k ) 1 -ω u ω u -k k 2 -( ω u -k -ω u ) 2 -v 4 ω v ∫ ∞ ϵ dk ∫ π 0 dθ sin 5 ( θ ) kN ( k ) 1 -ω v ω v -k k 2 -( ω v -k -ω v ) 2 ] ρ ( ⃗u, ⃗v, t ) . (A.V.50)  \n- [1] M. J. Fahn, K. Giesel, and M. Kobler, Class. Quant. Grav. 40 , 094002 (2023), 2206.06397.\n- [2] H.-P. Breuer et al. , The theory of open quantum systems (Oxford University Press on Demand, 2002).\n- [3] K. Hornberger, Introduction to Decoherence Theory (Springer Berlin Heidelberg, 2009), p. 221-276.\n- [4] U. Weiss, Quantum Dissipative Systems (World Scientific, 2021).\n- [5] G. Lindblad, Commun. Math. Phys. 48 , 119 (1976).\n- [6] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, J. Math. Phys. 17 , 821 (1976).\n- [7] F. Benatti and R. Floreanini, JHEP 02 , 032 (2000), hep-ph/0002221.\n- [8] M. P. Blencowe, Phys. Rev. Lett. 111 , 021302 (2013), 1211.4751.\n- [9] C. Anastopoulos and B. L. Hu, Class. Quant. Grav. 30 , 165007 (2013), 1305.5231.\n- [10] M. M. Guzzo, P. C. de Holanda, and R. L. N. Oliveira, Nucl. Phys. B 908 , 408 (2016), 1408.0823.\n- [11] T. Oniga and C. H. T. Wang, Phys. Rev. D 93 , 044027 (2016), 1511.06678.\n- [12] M. Lagouvardos and C. Anastopoulos, Class. Quant. Grav. 38 , 115012 (2021), 2011.08270.\n- [13] V. D'Esposito and G. Gubitosi, (2023), 2306.14778.\n- [14] A. Bassi, A. Großardt, and H. Ulbricht, Class. Quant. Grav. 34 , 193002 (2017), 1706.05677.\n- [15] C. Anastopoulos and B.-L. Hu, AVS Quantum Sci. 4 , 015602 (2022), 2111.02462.\n- [16] M. J. Fahn, K. Giesel, and R. Kemper, In preparation (2024).\n- [17] E. Lisi, A. Marrone, and D. Montanino, Phys. Rev. Lett. 85 , 1166 (2000), hep-ph/0002053.\n- [18] A. Sakharov, N. Mavromatos, A. Meregaglia, A. Rubbia, and S. Sarkar, J. Phys. Conf. Ser. 171 , 012038 (2009), 0903.4985.\n- [19] J. a. A. B. Coelho and W. A. Mann, Phys. Rev. D 96 , 093009 (2017), 1708.05495.\n- [20] J. A. Carpio, E. Massoni, and A. M. Gago, Phys. Rev. D 100 , 015035 (2019), 1811.07923.\n- [21] C. Burrage, C. Kading, P. Millington, and J. Min'aˇr, Phys. Rev. D 100 , 076003 (2019), 1812.08760.\n- [22] S. Weinberg, The quantum theory of fields: Volume 1, foundations (Cambridge university press, 2005).\n- [23] D. Tong, Quantum field theory (lecture notes), 2007.\n- [24] C. Fleming, N. Cummings, C. Anastopoulos, and B.-L. Hu, Journal of Physics A: Mathematical and Theoretical 43 , 405304 (2010).\n- [25] D. Burgarth, P. Facchi, R. Hillier, and M. Ligab'o, Quantum 8 , 1262 (2024), 2301.02269.\n- [26] P. Wang, E. Hiltunen, and J. C. Schotland, (2023), 2311.02670.\n- [27] A. Domi et al. , (2024), 2403.03106.\n- [28] Q. Xu and M. P. Blencowe, New J. Phys. 24 , 113048 (2022), 2005.02554.\n- [29] S. Deser, R. Arnowitt, and C. W. Misner, J. Math. Phys. 1 , 434 (1960).\n- [30] A. Ashtekar, Phys. Rev. Lett. 57 , 2244 (1986).\n- [31] J. F. Barbero G., Phys. Rev. D 51 , 5507 (1995), gr-qc/9410014.\n- [32] G. Immirzi, Class. Quant. Grav. 14 , L177 (1997), gr-qc/9612030.\n- [33] A. Ashtekar, J. D. Romano, and R. S. Tate, Phys. Rev. D 40 , 2572 (1989).\n- [34] C. Rovelli, Quantum Gravity Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2004).\n- [35] T. Thiemann, Modern Canonical Quantum General Relativity Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2007).\n- [36] T. Thiemann, Class. Quant. Grav. 12 , 181 (1995), gr-qc/9910008.\n- [37] E. Poisson and C. M. Will, Gravity: Newtonian, Post-Newtonian, Relativistic (Cambridge University Press, 2014).\n- [38] C. Rovelli, Class. Quant. Grav. 8 , 297 (1991).\n- [39] C. Rovelli, Class. Quant. Grav. 8 , 317 (1991).\n- [40] C. Rovelli, Phys. Rev. D 65 , 124013 (2002), gr-qc/0110035.\n- [41] A. S. Vytheeswaran, Annals Phys. 236 , 297 (1994).\n- [42] B. Dittrich, Gen. Rel. Grav. 39 , 1891 (2007), gr-qc/0411013.\n- [43] B. Dittrich, Class. Quant. Grav. 23 , 6155 (2006), gr-qc/0507106.\n- [44] B. Dittrich and J. Tambornino, Class. Quant. Grav. 24 , 757 (2007), gr-qc/0610060.\n- [45] R. Zwanzig, J. Chem. Phys. 33 , 1338 (1960).\n- [46] S. Nakajima, Prog. Theor. Phys. 20 , 948 (1958).\n- [47] T. Arimitsu and H. Umezawa, Prog. Theor. Phys. 74 , 429 (1985).\n- [48] Y. Takahashi and H. Umezawa, Int. J. Mod. Phys. B 10 , 1755 (1996).\n- [49] F. C. Khanna, A. P. C. Malbouisson, J. M. C. Malbouisson, and A. R. Santana, Thermal quantum field theory - Algebraic aspects and applications (, 2009).\n- [50] F. Benatti and R. Floreanini, Phys. Rev. D 64 , 085015 (2001), hep-ph/0105303.\n- [51] J. F. Donoghue, Phys. Rev. D 50 , 3874 (1994), gr-qc/9405057.  \n- [52] D. Arteaga, R. Parentani, and E. Verdaguer, Phys. Rev. D 70 , 044019 (2004), gr-qc/0311065.\n- [53] B. Hatfield, Quantum Field Theory Of Point Particles And Strings (CRC Press, 2019).\n- [54] F. Nathan and M. S. Rudner, Physical Review B 102 (2020).\n- [55] E. T. Jaynes and F. W. Cummings, IEEE Proc. 51 , 89 (1963).\n- [56] B. W. Shore and P. L. Knight, Journal of Modern Optics 40 , 1195 (1993).\n- [57] J. Larson and T. Mavrogordatos, The Jaynes-Cummings model and its descendants: modern research directions (IoP Publishing, 2021).\n- [58] G. Balieiro Gomes, D. V. Forero, M. M. Guzzo, P. C. De Holanda, and R. L. N. Oliveira, Phys. Rev. D 100 , 055023 (2019), 1805.09818.\n- [59] D. Boriero, D. J. Schwarz, and H. Velten, Universe 5 , 203 (2019), 1704.06139.\n- [60] D. Pfenniger and V. M. G. Observatory, Astron. Astrophys. 456 , 45 (2006), astro-ph/0605354.\n- [61] A. E. Bernardini, EPL 103 , 30005 (2013), 1204.1504.\n- [62] G. J. Milburn, Phys. Rev. A 44 , 5401 (1991).\n- [63] G. J. Milburn, New J. Phys. 8 , 96 (2006), gr-qc/0308021.\n- [64] L. Diosi, Braz. J. Phys. 35 , 260 (2005), quant-ph/0412154.\n- [65] H. P. Breuer, E. Goklu, and C. Lammerzahl, Class. Quant. Grav. 26 , 105012 (2009), 0812.0420.\n- [66] G. S. Adkins, Phys. Rev. D 27 , 1814 (1983).\n- [67] G. Kaplanek and E. Tjoa, Phys. Rev. A 107 , 012208 (2023), 2207.13750.\n- [68] T. Colas, J. Grain, and V. Vennin, Eur. Phys. J. C 82 , 1085 (2022), 2209.01929.\n- [69] D. Davidovi'c, Quantum 4 , 326 (2020).\n- [70] K. Giesel and M. Kobler, Mathematics 10 , 4248 (2022), 2207.08749.\n- [71] H.-T. Cho and B.-L. Hu, Phys. Rev. D 105 , 086004 (2022), 2112.08174."}