Objective trajectories in hybrid classical-quantum dynamics

Consistent dynamics which couples classical and quantum degrees of freedom exists, provided it is stochastic. This dynamics is linear in the hybrid state, completely positive and trace preserving. One application of this is to study the back-reaction of quantum ﬁelds on space-time which does not suﬀer from the pathologies of the semi-classical equations. Here we introduce several toy models in which to study hybrid classical-quantum evolution, including a qubit coupled to a particle in a potential, and a quantum harmonic oscillator coupled to a classical one. We present an unravelling approach to calculate the dynamics, and provide code to numerically simulate it. Unlike the purely quantum case, the trajectories (or histories) of this unravelling can be unique, conditioned on the classical degrees of freedom for discrete realisations of the dynamics, when diﬀerent jumps in the classical degrees of freedom are accompanied by the action of unique operators on the quantum system. As a result, the “measurement postulate” of quantum theory is not needed; quantum systems become classical because they interact with a fundamentally classical ﬁeld.


Introduction
We are often interested in the dynamics of composite systems where one system can be considered classical while the other must be treated quantum mechanically.In quantum thermodynamics and quantum chemistry, we often have small molecules interacting with a large thermal reservoir which can be treated classically.In measurement theory the quantum system interacts with a macroscopic device which can be considered classical.In gravity, macroscopic objects such as evaporating black holes radiate thermally, and we imagine there is a regime where we can treat the black hole space-time classically even though the radiation must still be described by quantum field theory.
There is a lot of debate in the literature on whether one can consistently couple quantum systems and classical ones, and many ways of doing so which although useful in some regimes, are pathological [10,21,20,29,17,67,63,62,64,5,45].Many proposals for such dynamics [12,26,3,42,41,61,65,66,11,27] are not completely positive or haven't been shown to be, and since both the density matrix of a quantum system and the phase-space density of a classical system are positive functions or matrices, such maps lead to negative probabilities, and are at best an approximation to the true dynamics.Other proposals, such as semi-classical gravity, are non-linear and thus don't respect the statistical nature of the density matrix [55,35].Another approach, inspired by quantum measurement and control [71], sources classical degrees of freedom such as the Newtonian potential via feedback and measurement of quantum matter [39,40,68,69].
If on the other hand, one doesn't wish to introduce auxiliary processes such as measurement or collapse, then one can still construct consistent evolution laws which couple quantum and classical degrees of freedom.Two such evolution laws were introduced in [9,23].These dynamics are completely positive, trace preserving, and preserve the split between classical and quantum degrees of freedom.These dynamics are special cases of the master equation (4), derived in [51,53].These evolution laws have been used to study the collisionless Boltzmann equation [2] and even proposed as a fundamental theory coupling quantum matter to classical Newtonian gravity [24,60] and quantum fields to classical space-time in the context of General Relativity [51].Because we don't have a quantum theory of gravity, the prospect of understanding how classical space-time and quantum fields interact is exciting, whether the dynamics are fundamental or merely effective.In addition to the potential to better understand black hole evaporation, an effective theory opens up the possibility of gaining new insights into cosmology, since the quantum nature of vacuum fluctuations are important in structure formation in an expanding universe.
We will not discuss gravity in the present work, nor address the question of what degrees of freedom can or should be considered classical.Rather the aim here is to introduce basic techniques, and study some simple systems in order to build up an intuition for how classical systems couple to quantum ones.In particular, we shall see that such coupling not only leads to decoherence of the quantum system [38,59,51,68], but also a "collapse" of the wavefunction, in the sense that the quantum system jumps to a pure state which can be uniquely determined from the classical degrees of freedom.In this sense, we will see that the state of the system can be considered to have some objectivity, unlike the purely quantum case where a decomposition of the system's density matrix into pure states (or evolution into a history) is not unique [43,28].We also study another feature of hybrid dynamics, namely a trade-off between decoherence and diffusion.Quantum systems which have long coherence times necessarily have high diffusion in the classical degrees of freedom.Here, we shall exhibit systems which have this feature.That there is necessarily a tradeoff between decoherence of the quantum system and diffusion of the classical one, will be shown in [54].
We shall here consider a qubit and a quantum harmonic oscillator that are coupled to classical degrees of freedom that are described by a phase-space manifold M, and for simplicity in the following we take M ≡ R 2n , the space of position and momentum z = q, p for n particles.The quantum degrees of freedom are described by a Hilbert space H, which might be finite or infinite dimensional.Given the Hilbert space, we define the set of positive semi-definite operators as S(H).The object describing the state of this composite system at a given time is the map : M → S(H) such that, when the quantum degrees of freedom are traced out, returns a valid probability distribution over the phase-space.We will call this the hybrid-density, and in this paper we are using natural units, setting = 1.The hybrid density is a probability density over a phase space, taking values in the collective variable z = (q, p).As such, it need not be subnormalized for each z but it must be a normalized distribution once integrated over all classical configurations.Mathematically, we have that the distribution prob(z) = Tr [ (z)] is such that prob(z) ≥ 0 for all z ∈ M, and it is normalized M dz prob(z) = 1.From the above property it is easy to show that, if the classical degrees of freedom are traced out, the resulting state ρ = M dz (z) is a valid quantum state (i.e., a positive semi-definite operator with unit trace).The simplest such system is the hybrid qubit, whose state can be written as (q, p) = u 0 (q, p) c(q, p) c (q, p) u 1 (q, p) If we integrate over phase space, we obtain the density matrix of the qubit, while if we take the trace of the density matrix, we are left with the phase-space density of a classical particle.
A master equation which governs the dynamics of this hybrid state must posses several necessary properties; it is linear in the state of the composite classical-quantum system, and it is completely positive and trace preserving, so that it's action preserves the statistical interpretation of (q, p).A valid quantum-classical state is transformed to another valid quantum-classical state, with the master equation preserving the separation between classical and quantum degrees of freedom.Two generators of this kind of dynamics have previously been proposed.In [8], a master equation with diagonal coupling to Lindblad operators was considered, which reduces to the GKSL or Lindblad equation under suitable assumptions.This dynamics was shown to be the most general master equation in the case of discrete classical variables and bounded Lindblad operators by embedding the classical system in Hilbert space [60].However, for continuous classical degrees of freedom, one generally requires an uncountable number of Lindblad operators.Díosi, instead, considered a master equation exhibiting a diffusion term, which generates a continuous dynamics for the hybrid system [23].These two CQ master equations are special cases of Equation (4) [51]. 1t and other hybrid dynamics have also been proposed as a mechanism which leads to fundamental decoherence of quantum states [38,68,69,60,51] In this paper, we are interested in understanding the main features of the discontinuous dynamics generated by the discrete master equation, given in Eq. (9).In order to study these properties, we generalise a well-known tools from open quantum systems, known as the unravelling of the master equation [7,18,30,34,14].In open quantum systems [1,13], one studies the evolution of a quantum system coupled to an external environment.The most general Markovian evolution of the quantum system given some reasonable assumptions is the GKSL or Lindblad equation [44,36], whose action over the set of quantum states forms a semi-group [19].The dynamics of a quantum state ρ are given by with H the Hamiltonian of the system, {•, •} + is the anti-commutator and L α , L β ∈ B(H) are Lindblad operators.Such dynamics arise, for example, when a quantum system is weakly coupled to a very large environment.
The unravelling is a technique that allows us to numerically simulate the evolution of a quantum system generated by the GKSL equation [7,18,30,16,46,47].Instead of considering the evolution of the full density matrix describing the system, in the unravelling one consider a single pure state, which is evolved according to a stochastic dynamics.Multiple paths in the Hilbert space are thus generated by the stochastic dynamics, and by averaging over the paths we recover the evolution of the quantum system.This procedure is computationally advantageous compared to evolving the full density matrix, in particular when the dimension of the system is large.Moreover, the unravelling technique provides a different perspective on the evolution of open quantum systems; indeed, the evolution of the system can be understood as generated by a continuous, deterministic dynamics (when the system does not interact with the environment), scattered by discrete jumps of the wave-function (when system and environment interact), which occur stochastically in time [30].
In the quantum case, there is generally not a unique unravelling for the master equation, so one cannot regard a particular trajectory as actually happening.Indeed, the decomposition of the dynamics in terms of Lindblad operators L α and a Hamiltonian H in Eq. (1), is not in general unique.One can rewrite the L α in terms of some other basis of operators.This presents a problem [43,28] for interpretations of quantum theory which take the trajectories as giving a theory of microscopic state reduction [70], and likewise for the decohering histories approach [49,37,31,25].In contrast, by extending the tools of unravelling to the master equation which generates the dynamics of classical-quantum systems, we shall find that, under the assumption that a unique classical shift is associated with a different Lindblad operator L α , trajectories are unique when one conditions on the classical degrees of freedom.This imbues the unravelled trajectories with an ontological significance and allows one to interpret the state of the system as "collapsing" to a particular state.This collapse of the wave-function is here generated by the interaction between classical and quantum degrees of freedom as in [38,68,69,60,51], rather than the introduction of an ad-hoc field as in spontaneous collapse models [32,57,33,6].Objective trajectories in unitary quantum theory have previously been studied in [72,56,73,15], where the environment plays this role of the classical system in the limit where we take it's size to infinity.Classical systems can be in definite states, and we find that a classical-quantum interaction causes the quantum system to inherit this property.This provides a new resolution of the measurement problem of quantum theory -quantum systems become classical because they interact with another system which is itself intrinsically classical (e.g.space-time as in [51]).
A simple example illustrates this point.Consider the purely quantum case of a qubit in an equal superposition of basis states |0 and |1 .If the quantum state ρ undergoes pure decoherence via the Linblad equation then it's state at any time is given by While we might be inclined to describe the state of the system as evolving continuously from |+ := (|0 + |1 ) / √ 2 to the equal mixture I/2, others may describe the system's evolution in terms of its state starting in |+ and then at some random time, suddenly collapsing to the |0 or |1 state at a rate given by λ.Still others may describe it as suddenly collapsing to any other two orthogonal states.There is however no physical meaning to these statements, since there is no way to perform an experiment which would distinguish between these different descriptions.In quantum theory, the density matrix completely determines the system, and any decomposition into an ensemble of states is arbitrary and has no physical meaning in this case.In contrast to this, the evolution law we will consider allows for the sudden change of the quantum state to be accompanied by a change in a classical degree of freedom.Since a classical degree of freedom can be monitored without disturbing the system, one can perform measurements on the quantum system conditioned on when the classical system changes, to verify that the quantum system has undergone a sudden jump and to what state, unambiguously.This example is discussed in Subsection 3.1.
However, one needn't make such an interpretational commitment.The unravelling can merely be considered to be a calculations tool which allows us to better simulate the dynamics of some simple hybrid systems and we present an algorithm and code for this purpose.We study both the hybrid qubit and hybrid oscillator numerically and analytically, and explore a relation between the diffusion of the classical degrees of freedom in phase space and the rate of decoherence of the quantum degrees of freedom.Interestingly, this relation provides a way to experimentally test whether the master equation can be used to construct a "post-" quantum theory of gravity [51].We also study the energy of the composite system, defined in terms of the sum of the classical and quantum Hamiltonian operators.Because this composite energy is not the generator of time-translations, there is no reason to expect it to be conserved [4] -Noether's theorem doesn't apply if the dynamics is not unitary.Nonetheless, insofar as one can have a classical-quantum Hamiltonian, one might want to keep violations of its conservation small.There are various mechanisms one can employ for this purpose which we will briefly mention in the discussion.In the gravitational setting, the more natural question to consider is whether the constraints are preserved in time, a discussion which we reserve for the future [52].
The paper is structured as follow.In Sec. 2 we briefly review the hybrid master equation describing the evolution of classical-quantum systems, and we extend the unravelling technique to this setting.In Sec. 3 we study two different toy models where the classical degrees of freedom interact with a finite-dimensional quantum system (a qubit).We use the unravelling technique to simulate the hybrid dynamics, and to identify its main features.In Sec. 4 we consider the case in which the quantum system is infinite-dimensional (a harmonic oscillator), and we repeat our analysis in this setting, with a special focus on the decoherence of this system.We conclude in Sec. 5, and suggest various ways these models could naturally be extended, for example, by having the quantum dynamics always dependent on phase space degrees of freedom or by adding in friction terms which would control the diffusion of the classical degrees of freedom in phase space.

The classical-quantum master equation and its unravelling
In this section, we show that the most general master equation governing the dynamics of a classical-quantum (CQ) system can be unravelled, allowing us to efficiently simulate the evolution of different finite-dimensional hybrid systems.Recasting the dynamics under the unravelling formalism provides a useful perspective for understanding the evolution of CQ systems.Indeed, from this point of view the hybrid system evolves continuously in its classical and quantum degrees of freedom, and its evolution is interrupted at random times by jumps in both the classical phase-space and the Hilbert space.It is worth noting that the jumps are here due to the interactions between the classical and the quantum system, and they are not, as in the case of open quantum systems, due to the presence of an external environment.Since one can always consider the state of a hybrid classicalquantum system to be the restriction of a pure quantum state on part of an enlarged Hilbert-space, the evolution of Eq. (4) might in fact be the result of unitary dynamics on the entire system.However, strong constrains would be placed on this global dynamics due to the form of the reduced one, and it is doubtful that such dynamics can be made completely positive and trace preserving.Dynamics which is completely positive on the reduced state of a restricted class of system and environment states, needn't be completely positive on the full state of the system when we include the environment [51].
Another difference between the unravelling for hybrid dynamics and the one used in open system dynamics concerns the objectiveness of the jumps.In the case of open quantum systems, a given dynamics can be obtained from an entire family of jumps [13].There are usually many ways in which one can decompose the dynamics in terms of jumps (Lindblad operators) in Hilbert space, and so one cannot think of the unravelling as representing actual trajectories that the system makes.In contrast, for some classes of hybrid dynamics, each classical jump in phase-space uniquely determines which Lindblad operator was applied to the system.One can then imagine continuously monitoring the classical system through time without disturbing the system, enabling one to uniquely determine the sequence of Lindblad operators, and thus the evolution of the quantum state, which remains pure, conditioned on the classical degrees of freedom.We can thus think of the classical trajectories as being ones which objectively occur.
Before we can extend the unravelling technique to the case of classical-quantum dynamics, we need to formally introduce the class of hybrid systems we consider, and the form of the master equation describing their evolution.The most general Markovian CQ dynamics in the case of bounded operators L α is given by where z = (q 1 , p 1 , q 2 , p 2 , . ..) ∈ R 2n is the vector of phase space coordinates for n systems, and (z, t) is a classical-quantum state at time t.The Hamiltonian H(z), which appears in the commutator with the CQ state, controls the unitary evolution of the quantum system, and in principle depends on the classical degrees of freedom.For each α and β, the rate W αβ (z|z − ∆) is non-negative, and it governs the transition of the classical degrees of freedom from z − ∆ to z, as well as the jump of the quantum state due to the map L α • L † β .For the master equation to preserve the norm of the quantum state (once the classical degrees of freedom have been traced out), we need to require that ( The hybrid master equation ( 4) is completely positive (CP) over the classical and quantum degrees of freedom, and trace-preserving (TP) in the sense that the normalisation condition dzTr [ (z)] = 1 is preserved.Furthermore, its structure is such that the separation between classical and quantum degrees of freedom is preserved i.e. (z) is always a un-normalized density matrix.Finally, the master equation is clearly linear in the CQ state.This ensures that the theory is operationally non-signalling: if one considers two space-like separated systems then if the dynamics is local, local operations on one of the systems does not alter the reduced density matrix of the other systems, just like in standard quantum theory.All these properties ensure that a CQ state is mapped into another CQ state by the dynamics generated by Eq. (4).
In the previous paragraph we mentioned how Eq. ( 4) is the most general Markovian hybrid which is completely positive.It is important to notice that the dynamics is Markovian on the whole hybrid system, that is, on both the classical and quantum degrees of freedom.However, when we restrict our focus to either the classical or quantum degrees of freedom, this is not the case anymore.Indeed, the classical (quantum) system can act as a memory for the state of the quantum (classical) one, so that the reduced dynamics is not Markovian.
In Ref. [53], we have since shown that the master equation (4) gives rise to two different classes of hybrid dynamics.On the one hand, one can obtain a Fokker-Plank type of equation, where the classical-quantum degrees of freedom evolve continuously.This dynamics has been studied in [23] in the case where the dynamics is Markovian even when restricted to either the classical or quantum.One can also obtain dynamics which, as we shall see here, has finite jumps in the classical-quantum degrees of freedom.As we stressed in the introduction, in this paper we focus on the latter kind of dynamics, and explore the properties associated with a discontinuous hybrid dynamics.
Although Eq. (4) gives the general form of the dynamics, we are interested in evolution which in some limit reproduces classical Hamiltonin dynamics on the entire hybrid system.This was done in [51] by expanding W αβ in terms of moments of ∆ and demanding that the 1st moment reduce to the Poisson bracket of the state with some Hamiltonian.The particular form we will explore here, is given by taking matrix W αβ (z|z − ∆) defined through the following equation: where τ αβ > 0 can be understood as a rate, {•, •} is the Poisson bracket, and the functions h αβ (z) are associated with the CQ interaction Hamiltonian The Poisson bracket in the exponential acts on the CQ state as a linear operator: The resulting master equation is The dynamics in this equation are completely determined by the choice of the free parameters τ αβ , and the choice of the Lindblad operators L α and the functions h αβ (z).
We see that to 0'th order in the τ αβ , this master equation reproduces Hamiltonian dynamics in the appropriate limit.Namely, if we expand the exponential in the right hand side of the equation for small values of the parameters τ αβ 's, we obtain where the parameters τ αβ 's play the role of the relaxation rates.The first term of Eq. (10) is the dynamics of the quantum system, which can in principle depend on the classical degrees of freedom.In this way, it represents the evolution of the quantum system as controlled by the classical one.The terms on the first line of Eq. (10) which include a summation over α, β weighted by 1 τ αβ define purely Lindbladian dynamics and in the limit that τ αβ → 0, it can lead to rapid decoherence of the quantum system.The final term is the interaction Hamiltonian dynamics we desire, namely under the trace, it is Tr [{H, }].
We can think of it as the counter-part to the commutator in the sense that it gives the back-action of the quantum system on the classical one, while the commutator reflects the influence of the classical system on the quantum one.Indeed, in the classical limit, we expect −i[H(z), ] to act as the Poisson bracket with respect to the degrees of freedom of the quantum system, while in this limit one can verify that {h αβ , L α L † β } are the additional terms of the Poisson bracket with respect to the degrees of freedom of the classical system.We can take purely classical dynamics to correspond to identity Lindblad operators which we take to be L α=0 = I.If we wish this classical dynamics to be purely deterministic, then we can take τ 00 → 0, because in this limit, the only terms which remain in the expansion of the exponential are the pure Lindbladian term which acts only on the quantum system and the deterministic term in the dynamics given by {H C , } with the classical Hamiltonian defined as Higher order terms can lead to diffusion in phase space.For example, at first order in τ , Eq. (10) will acquire terms of the form where for ease of presentation we have taken all τ αβ = τ and used the Einstein summation convention to suppress the summation over α, β on the right hand side of Equation (12).Such terms lead to diffusion in phase-space, with D αβ p,ij (z) being the diffusion term usually found in the Fokker-Planck equation, or arising from Langevin dynamics.When W αβ (z|z− ∆) is deterministic in q, but diffusive in p, we shall refer to it as dynamics of the Langevin type.It is possible to also have terms appearing in the expansion which give rise to friction.We discuss such master equations in greater detail in Sec. 5.If the dynamics is of the Langevin type, then one can verify that the equation of motion for q is unchanged by the diffusive terms appearing in Eq. (10).Inverting this expression gives p(q, q).On the other hand, master equations which contain terms D αβ qp,ij and D αβ q,ij (z) in their expansion, we shall call q-diffusive and we will discuss such master equations in more detail elsewhere.Care should be taken in any expansion of W αβ (z|z − ∆) -although the full dynamics is completely positive, if we truncate our expansion to any finite order, it is not.Note that because the jumps are in classical momenta, and not position, the dynamics is local and there is no signalling.However, if we considered evolution operators which jump in position, we would encounter some violation of faster-than-light signalling.To address this problem, we have to go to a generalization of the classical-quantum theory, which involves Lorentz symmetry and quantum field theory, which we study in [50].

Unravelling of the CQ dynamics
In the quantum open system formalism, the unravelling of a master equation takes the master equation governing the evolution of the density matrix, and divides it up into trajectories -i.e.probabilities assigned to an ensemble of histories, where each history is the evolution of a pure quantum states which is described as undergoing a continuous evolution punctuated by stochastic jumps, where the system is mapped from one state to another as a result of a Lindblad operator being applied to it.Due to the stochastic nature of the unravelling, the initial state is mapped into different final states by the same evolution.If we record each evolution in time (or trajectory) and we average over them, we re-obtain the evolution law for the density matrix, as given via the master equation (in the limit of infinite trajectories and infinitesimal time steps).As noted in the Introduction, this formalism is merely one of many possible descriptions, and doesn't have physical meaning.If we consider only the quantum system, there is no measurement which can distinguish between continuous evolution of the density matrix in the Hilbert space, and one which is described in terms of jumps, nor is there physical meaning to what basis is used to describe the jumps.
On the other hand, while the unravelling of the CQ master equation proceeds in analogy with the quantum one, it can have physical content.In this section we provide an algorithm to unravel the dynamics given by Eq. (4), where in addition we consider a fully classical term, given by the Poisson bracket between a classical Hamiltonian H C (z) and the hybrid state, Specifically, we will provide the tools to generate different classical-quantum trajectories, and we will show that the average of these trajectories coincides with the evolution of the CQ state under the above master equation.
In the following, we consider the case in which the matrix of rates W αβ (z|z − ∆) are diagonal in α, β.If this matrix is not diagonal in the same basis at each point, we can introduce unitary operators u(z) that diagonalize it, and use the new Lindblad operators in the unravelling procedure, L γ (z) = u γ,α (z) L β .In this way, the same Lindblad operator acts on the left and the right of the CQ state in Eq. (4), and we can use the same procedure explained in the following, with the difference that now the Lindblad operators additionally depends on the phase space coordinate z.
The central object of this study is a pure CQ state, that is, where |φ(t) is a quantum state independent of the classical degrees of freedom.This object represent a classical-quantum system, centred in the point z(t) of the phase space, and described by the quantum state |φ(t) .The first step in the unravelling procedure is the discretisation of the dynamics of the CQ state; we fix a time interval δt, and the evolution of the CQ state is given for multiples of this interval.It is worth noting that the unravelling evolution is exact at first order in this parameter, and coincides with the one given by the master equation in the limit of δt → 0. We now outline the procedure for generating each trajectory, and we will then show that when we average over trajectories, we obtain the CQ master equation.For a given point z in phase space, at each time step we either update the state with a continuous evolution, or with a jump in the classical and quantum degrees of freedom.When the evolution is continuous, we apply to the quantum part of the CQ state a unitary evolution given by the (non-Hermitian) effective Hamiltonian which is composed by the Hamiltonian H(z), which appears in the commutator of Eq. (4), and by the operator appearing in the anti-commutator of the same equation.The classical part of the CQ state is instead subjected to a translation in phase space generated by the classical Hamiltonian H C (z), which appears in the Poisson bracket in the master equation.
In order to evolve the CQ state with a jump, we first need to choose at random a value of the tuple (α, ∆) from the (un-normalised) distribution W α (z|z − ∆).The α parameter specifies which Lindblad operator L α is performed over the quantum state, while the ∆ parameter quantifies the shift performed over the classical degrees of freedom.The unravelling algorithm is then defined by the following set of instructions, that specify how to update a pure CQ state at time t under continuous or jumping evolution, with p α,∆ (z, t).
(16) At each time step, the hybrid pure state can be updated either continuously, first line in the right hand side of the above equation, or with a jump, second line.The type of evolution is chosen at random; the probability of updating the state with a continuous evolution is given by p 0 (z, t), while the probability of evolving it with a jump is given by p α,∆ (z, t), where α specifies the quantum jump and ∆ specifies the classical one.In the continuous evolution, the quantum degrees of freedom are evolved by the unitary operator U (z) = e −iδt H eff (z) , that we have expanded at first order in δt, while the continuous shift in the classical degrees of freedom is generated by the operator e δt {H C (z),•} .The gradient vector is the symplectic matrix.By iterating the procedure in Eq. ( 16),we can produce a single trajectory.Different iterations provide different trajectories, since at every time step the CQ state is evolved either continuously or with a jump.The pure CQ states produced by the unravelling technique are normalized.In the following, we provide the normalization coefficient for both the continuous and jumping evolution, and the probability that either evolution occurs.For simplicity, in the rest of the section we suppress the time t dependence.The normalization coefficient are given as where we have used the form of the CQ pure state, given in Eq. ( 14), and we are only considering terms in the first order of the time step δt.In the next section we explicitly derive these coefficients.At a given time step, the probability of evolving the pure CQ state continuously or with a jump is given by, These probabilities are chosen such that the unravelling procedure, given the appropriate averaging, reproduces the dynamics of the master equation, which is discussed in the next subsection.It is easy to see that the above distribution is a valid probability distribution over the phase space.Clearly, each probability distribution is non-negative, and the distribution is normalised, where the last equality follows from Eq. ( 5).
Before concluding the section, let us comment on an interesting feature of the unravelling of the hybrid master equation, which is not shared by its quantum counterpart.If a unique shift ∆ in the classical degrees of freedom is associated with each Lindblad operator L α , then by monitoring the classical degrees of freedom we have complete knowledge on the jumps that occur in the quantum part of the hybrid state.Furthermore, the quantum evolution conditioned on the classical degrees of freedom is unique, since only a specific sequence of jumps could produce the sequence of shifts in the classical part.The hybrid state can be said to actually make the transitions given in Eq. ( 16).This is in contrast with the unravelling for open quantum systems, where the lack of a classical system which can be monitored makes it impossible to single out different trajectories in the Hilbert space.
Clearly, the above holds as long as the assumption that a unique shift ∆ in the classical degrees of freedom is associated with each Lindblad operator L α holds.In this paper, we consider toy models where we demand the assumption to hold.However, it is interesting to understand when the assumption holds in a physical setting, when there is an interaction Hamiltonian H(z) as in Eq. (7).From Eq. (9) we find that the classical jumps are given by the Hamiltonian components h αβ (z).If these components are different for each value of α and β, then each shift is uniquely assigned to a different Lindblad operator.Notice that additional freedom is given by the choice of the τ αβ 's, which contribute to the shifts.Then, in any physical situation where classical-quantum fields interact with different coupling strengths, one should expect the above assumption to hold.
We have since shown in [? ] that the continuous master equation has a unique unravelling if it saturates an inequality we call the decoherence-diffusion trade-off.

From unravelling to the CQ master equation
We can now show that the update rule presented in Eq. ( 16) reproduces, at first order in δt, the same dynamics of the hybrid master equation.First, let us derive the normalization coefficients shown in Eqs.(18).To compute the normalization of the continuously evolved state, we need to Taylor expand it and truncate the expansion at first order in δt.It is straightforward to show that The normalization coefficient is then given by where the second line follows from the Taylor expansion in Eq. ( 21), from the form of the effective Hamiltonian H eff (z) and from Leibniz's rule of the Poisson bracket, while the last line follows from the fact that | (z, t) is normalized.The normalization of the CQ pure state after a jump is simply We can now derive the master equation (4) (equipped with the additional fully-classical term) from the update rule given in Eq. ( 16).To so so, we express the state | (z, t + δt) as a density operator, and we re-write it as a mixture over the possible evolutions the system undergoes.Let us first consider the continuous update for the CQ density operator, which we refer to as ρ 0 (z, t + δt), where we truncate at first order in δt, where we have used the definition of H eff (z), see Eq. (15), and [•, •] is the commutator while {•, •} + is the anti-commutator.When the evolution is given by an α-jump in the quantum degrees of freedom, and by a ∆-jump in the classical ones, we find the following density operator We can now express the overall evolution of the state by weighting the different updates with the correct probabilities, given in Eq. (19).Thus, we find that and if we rearrange the above equation by moving | (z, t) (z, t)| from the right to the left hand side, we divide by δt, and we send δt → 0, we obtain As a last step, we need to move from a CQ pure state | (z, t) to a (possibly mixed) CQ state (z, t).To do so, we perform an average over the classical-quantum configurations of the initial state, that is, we average over a phase space distribution P ps of initial points z(t = 0) = z0 , as well as a distribution P q over the initial quantum states |φ(t = 0) = |φ 0 .Thus, by applying this average over the CQ pure state | (z, t) , we obtain If we now perform the average in Eq. ( 25), we obtain the CQ master equation.

Main features of hybrid dynamics in qubit toy models
In this section we present a few toy models to illustrate the main features of the hybrid dynamics generated by the master equation (4).The simplest hybrid dynamics we can consider is that of a spin half particle in a classical potential, which in our examples is taken to be linear.For each of the toy models we derive, both analytically and numerically, the dynamics of the corresponding hybrid system.To numerically evolve the CQ state we use the unravelling procedure presented in the previous section.We tailor the unravelling to the specific master equation used for the toy models, see Eq. (28).This unravelling procedure is shown in Appendix B.
The main features of the hybrid dynamics, that are highlighted in our toy models, are i) the presence of stochastic collapses of the quantum degrees of freedom, ii) a trade-off between the classical diffusion in phase space and the quantum decoherence, and iii) the fact that energy is not conserved by the dynamics.More specifically, we find that the interaction between classical and quantum degrees of freedom generates a sudden collapse in the latter.This collapse occurs when the CQ state is subjected to a jump, that is, when the classical degrees of freedom are shifted in phase space and a Lindblad operator acts over the quantum state.If the Lindblad operators are associated with unique shifts in the classical degrees of freedom, we find that the quantum dynamics can be unambiguously recorded by monitoring the classical one, see for example Fig. 1.Furthermore, we find that the rate of decoherence of the quantum degrees of freedom is linked to the rate of diffusion of the classical ones; the faster the quantum state decoheres, the slower the classical degrees of freedom spread in phase space.This is not a specific feature of the toy models we study, but rather a general property of the hybrid dynamics, as shown in Ref. [54].Finally, we find that the average energy of the toy models studied is not conserved by the hybrid dynamics.This should not come as a surprise, since the Hamiltonian operator appearing in Eq. (4) is not the generator of the dynamics.A detailed study of the conservation laws and symmetries present in the hybrid dynamics is performed in Ref. [52].
Let us introduce the toy models we study in this section, namely, a single spin half particle in a classical potential.The position and momentum of the particle are taken to be the classical degrees of freedom z = (q, p), while the spin of the particle (s = 1  2 ) is the quantum one.To simplify the notation, we take the transition matrix in the hybrid master equation (9) to be diagonal in α, β (if the transition matrix is not diagonal, we have shown in the previous section that we can always diagonalise it).For a spin-1 2 particle, the most general quantum Hamiltonian can be written as λ • σ where λ is a vector in R 4 and σ = {I, σ x , σ y , σ z } is the vector of Pauli matrices.By having this depend on the classical degrees of freedom, we cause the classical and quantum system to interact.Let us take the interaction term of the CQ-Hamiltonian to be where the coupling constant B(q) depends on the position of the particle q, and the operator is diagonal in the basis {|0 , |1 }.By comparing the above interaction Hamiltonian with the general form of Eq. ( 7), we can easily see that we are left some freedom in the choice of the Lindblad operators {L α } α =0 and of the functions h α (z).By exploiting this freedom, we will study two different CQ models describing a particle in a linear potential.In both cases, we choose a classical Hamiltonian H C (q, p) = p 2 2m , that is, the one for a free particle of mass m.One could imagine a purely quantum evolution determined by a Hamiltonian which is not dependent on classical degrees of freedom, but in analogy with the gravitational case we will not consider it here.Finally, we demand the rates τ α =0 's to be all equal to τ > 0, so that no non-trivial Lindblad operator acts more often than the others.
The CQ master equation (9) for the case we are considering is then given by The first term on the right is the evolution of the quantum degrees of freedom and it depends on the position of the classical particle.The next term is the purely classical evolution.If we take the rate τ 0 to be a finite positive value, this classical evolution is stochastic, while if it tends to 0 it is determinstic.In the former case, we will see that the equation can be solved analytically.When τ 0 → 0, the purely classical evolution is given by the Poisson brackets between the classical Hamiltonian and the CQ state, {H C (z), (z, t) }.In this second case, we analyse numerically the dynamics generated by the master equation, using the unravelling method shown in App.2.1.The final set of terms, gives the backreaction of the qubit on the classical degrees of freedom.

Qubit evolution with diagonal Lindblad operators in a linear potential
The first toy model we consider is a spin half particle interacting with a linear potential through Lindblad operators that are diagonal in the interaction Hamiltonian eigenbasis {|0 , |1 }.We make the following choice for the operators, Once the Lindblad operators are defined, they fix uniquely the functions h α (q, p).Indeed, in order to identify the interaction Hamiltonian H I (q, p) in Eq. ( 27) with the expression in Eq. (7), we need to define these functions as It is worth noting that, for a given interaction Hamiltonian, the above choice of Lindblad operators and h α (q, p) functions is not unique.In the next section, for instance, we consider non-diagonal operators that nonetheless return the same interaction Hamiltonian, but generate a completely different dynamics.Furthermore, one could chose the same Lindblad operators as we did above, but with a different normalization.To re-obtain the correct Hamiltonian H I (q, p) we then need to scale the functions h α (q, p) consistently.Notice that the normalization of the Lindblad operators influences the rate at which the operator is applied, while the function h α (q, p) regulates the amplitude of the shifts in the classical degrees of freedom.This freedom in the Lindblad operators and the functions h α (q, p) is linked to the fact that the interaction Hamiltonian is not the generator of the hybrid dynamics of Eq. (4), and this is highlighted by the fact that this operator is in fact not conserved by the evolution.
In this example, we consider a potential linearly dependent in position, so that B(q) = qB, where B has a constant value.We can express the CQ state (q, p, t) as where u i (q, p, t) is the population of the quantum state |i (for i = 0, 1), and c(q, p, t) is the coherence.We now re-write the CQ master equation given in Eq. (28), in terms of the populations and the coherence of the state (q, p, t).In the equations below, we take τ 0 → 0, so as to obtain a system of non-local differential equations.When we solve the equations analytically, however, we will require τ 0 to be a positive (albeit small) constant.The dynamics of the populations and coherence is thus given by It is worth noting that the equation for the coherence can be solve analytically, and the solution has the form for any function c(q − p m t) in C 1 , the set of all continuously differentiable functions.We thus see that the coherence term decays exponentially fast.We will see that this is a result of the quantum system collapsing to the 0 or 1 state while making a momentum jump.

Analytical and numerical evolution of populations
In order to study the evolution of the populations analytically, we can re-express Eq. (32a) as a stochastic equation, see App. A. To do so, we do not send to 0 the rate τ 0 in Eq. ( 28), and we obtain the following evolution for the populations, The solution of this equation is given in the appendix for general Hamiltonians, see Eq. (94).For the specific case we have are considering, the solution is −0.009 −0.003 0.003 Figure 1: The evolution of the population u 0 (q, p, t) and u 1 (q, p, t) in phase space, under the CQ qubit dynamics with diagonal Lindblad operators.The initial CQ state is δ(q)δ(p) |+ , where The interaction with the field makes the quantum state stochastically collapse into either |0 or |1 at a rate given by ω τ .The same collapse shifts the classical momentum of the particle by ±Bτ , depending on the quantum state it projects onto.By monitoring the classical degrees of freedom, we can establish the state of the quantum system.Indeed, a particle in the state |0 has negative momentum, whereas a particle in |1 has positive momentum.The evolution of the populations is shown for t = 0.02 s, 0.06 s and 0.1 s.The state is evolved for t = 0.1 s, with time steps of δt = 10 −4 s, and the jump rate is τ = 10 −2 s.The constant B = 1 J • s • m −1 , the mass m = 1 kg, and the frequencies are where k is the number of position shifts, n is the number of momentum shifts, and both P 0 (k) and P 1 (n) are Poisson distributions with mean value λ 0 = t τ 0 and λ 1 = t τ , respectively.The function u i (n, k) is given by where S n,k is a proper subset of the set of all permutations of the shift operators, with k momentum shifts and n position ones.Each element π creates a different combination of the n + k shifts.
The model under consideration can additionally be solved numerically, using the unravelling method presented in the previous section.In particular, we can use Eq.(28) to better understand the evolution of the CQ state.In this equation, both classical and quantum degrees of freedom evolve either continuously, or with a sudden jump.The latter evolution is the most interesting, since during a jump one of the Lindblad operators is applied to the quantum state, and this state is projected either in |0 or in |1 , so that the coherence of the state is destroyed, and the classical momentum receives a kick proportional to Bτ , whose direction depends on the quantum state.From this point of view, this hybrid model mimics some of the features of a standard Stern-Gerlach experiment, namely the fact that when the magnetic field measures the quantum state (and collapses it into |0 or |1 ), it kicks the particle (that is, increase the classical momentum) either upward or downward, depending on the state the spin is collapsed into.In Fig. 1 we provide the numerical solution for the evolution of the populations, obtained using the unravelling technique of Sec.(2.1).This is the example mentioned in the Introduction and contrasted with the purely quantum decoherence of Equation (2).For the Lindblad equation, and the state of the system initially in the |+ state, there is no meaning to the statement that the quantum state starts suddenly collapses to |0 or |1 , while here, one can continuously monitor the classical momentum without disturbing the quantum degree of freedom.Conditioned on finding a jump in momentum, one can then measure the qubit to verify that it is now in the |0 or |1 state, depending on the direction of the jump, while any previous measurement of the qubit in the ± basis, will yield the initial |+ state.0.000 0.025 0.050 0.075 0.100 t (s) 0.000 0.002 0.004 q (m) q fit σ q fit q analytic σ q analytic 0.000 0.025 0.050 0.075 0.100 t The mean value and standard deviation of the population u 0 (q, p, t) in phase space.Left.The mean value q and standard deviation σ q of the marginal distribution u 0 (q, t) (where we traced out over the momentum degree of freedom).The solid lines are the values obtained from the numerical simulation, whereas the dashed line are obtained from the analytical solution.Centre.The mean value p and standard deviation σ p of the marginal distribution u 0 (p, t).Right.The coherence |c| 2 of the quantum state located in q = 0 m and p = 0 kg • m • s −1 , as a function of time.It is worth noting that the numeric solution is obtained with the unravelling method, in which the position is updated during the continuous evolution step, see Eq. (118).Instead, the analytic solution is obtained for the case in which the position evolves through jumps, see Eq. (34).In order to match the two evolution, we have to ask the stochastic rate τ 0 to be of the order of the infinitesimal time steps δt used in the simulation.The mean value and variance of position and momentum are therefore computed using Eqs.(37) and (38) for a value of τ 0 = δt.

Diffusion in phase-space and trade-off with decoherence
We can also study the diffusion in phase-space of the populations of the CQ state, both numerically and analytically.Let us take the initial state of the CQ system to be at the origin of phase-space (position q(t = 0) = 0 and momentum p(t = 0) = 0), and in the quantum state In the unravelling picture when a quantum jump occurs, the state jumps to either |0 or |1 due to the application of the Lindblad operators of Eq. (29) to the state.After the jump, the quantum state cannot change anymore, but additional jumps do increase the momentum of the particle.Since the initial condition is symmetric in |0 and |1 , as the evolution of the populations is, we can focus on the population u 0 (q, p, t).The kind of decoherence affecting the quantum degrees of freedom in this model can be understood in terms of a leakage of information about the quantum degrees of freedom into the classical degrees of freedom.Indeed, each Lindblad operator L α is here associated to a distinct classical jump in phase-space, which allows an observer monitoring the classical degrees of freedom to know exactly in which quantum state the system is in and when the transition occurred.Thus, the possibility of monitoring the quantum system using its position in phase space removes any coherence in the basis {|0 , |1 }.Furthermore, by conditioning over the classical degrees of freedom, the quantum state remains pure and does not become mixed.This is different to the standard decoherence found in the usual Lindblad equation.There, the decomposition of the dynamics into Lindblad operators is not unique, and there is no physical meaning to the jumps -rather, the density matrix of the system evolves continuously in time with the off-diagonal matrix elements slowly decaying with time.In contrast, here, the quantum jump is accompanied by a discontinuous jump in p.The quantum jump thus has physical meaning, since an observer who is monitoring the classical degree of freedom can verify the quantum jump.If they observe a sudden increase (decrease) in momentum, they will expect that the quantum state is now in the |1 (|0 ) state, and can verify this by measuring the quantum state.The crucial ingredient here is that each classical jump corresponds to a single Lindblad operator being applied to the quantum state.
We will now see that the quantum system's coherence time is related to the diffusion in phase space.In order to obtain the diffusion constant from the numeric results, we compute the marginal distribution for position and momentum as a function of time, and we fit them with a Gaussian distribution to extrapolate mean value and variance.These quantities can also be computed using the solution in Eq. ( 35), as we show in App.A.2.1.From the analytic solution we find that the mean value of position and momentum is given by Bω 0 m t 2 + q (0), (37a) which is the expected solution for the position and momentum of a particle in a linear potential given by V (q) = qBω 0 .We notice that the expectation values of position and momentum are independent of the parameter τ 0 , which can therefore be sent to 0. The variance of position and momentum can be obtained too, as we show in the Appendix A, and are given by The variation in momentum is due to the fact that the number of momentum jumps will be normally distributed and as with a random walk, increase like √ t, while the variation in position is a consequence of the momentum having a variation.We can additionally link the variance in momentum to the diffusion term appearing in the expansion of the hybrid master equation, see Eq. (12).To do so, we need to expand the exponential operator in Eq. (28) to second order, obtaining the following diffusion term, D (u 0 (q, p, t)) = 1 τ τ h α=1 (q, p), τ h α=1 (q, p), u 0 (q, p, t) = (Bω 0 ) 2 τ ∂ 2 ∂ p 2 u 0 (q, p, t).
(39) It is then easy to see that the variance in momentum arises from the diffusion coefficient D α=1 p (q, p) = (Bω 0 ) 2 τ .This is consistent with the fact that the expansion of the master equation at first order in τ has the form of a Fokker-Planck equation with the above diffusion term.
In Fig. 2 we compare the numerical results with the ones obtained analytically, and we show that indeed the numeric simulation accurately describes the evolution of the CQ state.It is worth noting that, while the decoherence of the quantum system is inversely proportional to τ as seen from Eq. (33), the diffusion in phase space is directly proportional to it.Thus, we can have two opposite situations; for τ 1, the state of the quantum system quickly collapses and does not diffuses much in phase space, following an almost Liouvillian deterministic evolution.When τ 1, instead, the quantum system slowly decoheres, but the classical degrees of freedom quickly diffuse in phase space.We thus have a trade-off between the decoherence rate and the diffusion rate.This turns out to be a standard feature of hybrid dynamics and can be understood in terms of the moments of the Kramers-Moyal expansion of the CQ master equation, as proven in [54].

Energy conservation
We now turn to the question of energy conservation in this model of hybrid dynamics.We will see that there is a violation of energy conservation at a rate proportional to τ , although it can be made arbitrarily small.The average energy of the system can be computed using the total Hamiltonian H(q, p) = H C (q, p) + H I (q, p), and the solution obtained in Eq. (35) for the populations.Indeed, the average energy of the CQ system is given by Let us consider, for simplicity, that the initial state of the CQ system is (q, p, 0) = δ(q)δ(p) |0 0|, so that we only need to consider the energy contribution of the population u 0 .Furthermore, given the specific initial state of the system, we can express the ensemble u 0 (n, k), shown in Eq. (36), as u 0 (n, k) = dqdp P (q, p | n, k) u 0 (q, p), where u 0 (q, p) = δ(q − q)δ(p − p).The evolution of the population u 0 can therefore be expressed as   u 0 (q, p) = dqdp P (q, p) u 0 (q, p), ( 41) and by replacing it in Eq. (40) we find that the average energy is where the mean value and variance of the marginal probabilities P (q) and P (p) have been computed in the previous section, in Eqs.(37) and (38) respectively.Using these results, we can show that the energy associated with the classical Hamiltonian is given by H , while the quantum coupling contributes to the energy as H I (t) = − (Bω 0 ) 2 2m t 2 (notice that in the previous section we have considered the absolute value of the position for convenience).As a result, we find that the average energy of the system is not conserved, and instead grows linearly in time from the initial value of H (0) = 0, see Fig. 3.In the limit when τ tends to 0, one recovers the classical evolution given by the Liouville equation, and energy is conserved, but otherwise the energy increases at a rate proportional to the jump distance τ .
As we stressed at the beginning of the section, the fact that the average energy is not conserved does not come as a surprise.Indeed, we are taking the energy to be given by the operator H(z) even though it is not the generator of the symmetry of time-translation.In this model however, the failure of the system to conserve energy is not due to changes in the quantum system -the interaction of the qubit with the classical system merely causes it to decohere in a basis which commutes with the total CQ Hamiltonian.Instead, the term appearing in Eq. ( 43) is due to dispersion in the momentum of the classical system.This can be prevented by adding in a friction term, as discussed in Sec. 5.

Qubit evolution with non-diagonal Lindblad operators in linear potential
We now consider the same hybrid system, with the difference that now the Lindblad operators are non-diagonal in the interaction Hamiltonian eigenbasis.We will see that the dynamics of this CQ system is qualitatively and quantitatively different from the one in 0.000 0.025 0.050 0.075 0.100 t (s) 0.000 0.002 0.004 The solid lines are the numerical results, while the dashed ones are the analytical ones.The energy contribution given by the kinetic energy, E C = H C , grows faster than the one given by the potential energy, E I = H I .As a result, the total energy of the system, E tot = H , increases linearly as a function of time.
the previous section.The form of the Lindblad operators is In order to obtain the interaction Hamiltonian H I (q, p) of Eq. ( 27) from these operators, the functions h α (q, p) are h α=1 (q, p) = ω 0 qB, (45a) where again the particle is in a linear potential, i.e., B(q) = qB, with B = const.If we express the CQ state (q, p, t) as in Eq. (31), and we take the rate τ 0 → 0, then the master equation given in Eq. (28) can be expressed as a system of non-local differential equations.
In particular, the evolution of the populations is given by of the following differential equations, where ⊕ is addition modulo 2. Notice that the above equations have the same form as the ones obtained in the previous section, see Eq. (32a), but now each population depends non-locally on the other.The evolution for the coherence is unchanged from the one in Eq. (32b).

Analytical and numerical evolution of populations
To solve the model analytically we can express Eq. ( 46), which describes the evolution of the populations, as a stochastic equation, see App. A. To do so, we fix the rate τ 0 in Eq. (28) to be a (small) positive constant, so that the equation for the populations can be expressed as If τ 0 → 0, we come back to the original equation for the populations shown in Eq. (46).In order to solve this equation we can make use of the tools described in appendix, see ) is a coherent superposition of the eigenstates of the interaction Hamiltonian H I .Notice that, in the rest of the section, we instead focus on the case where the initial quantum state is |0 , which simplifies the analytical study of the dynamics; the qualitative behaviour of the phase-space evolution, however, is the same as the one shown in this figure.This interaction makes the quantum degrees of freedom stochastically collapse from |0 to |1 and vice versa, at a rate given by τ .When the quantum state is mapped as |0 → |1 , the classical momentum is increased by Bτ , while the transformation |1 → |0 is accompanied by a change in momentum equal to −Bτ .The simulation depicted here uses time steps of δt = 10 −4 s, and a jump rate τ = 10 −2 s.The interaction constant is B = 1 J • s • m −1 , the mass is m = 1 kg, and the frequencies are Sec. A.3, and the solution we obtain is, where P 0 and P 1 are Poisson distributions with mean value t τ 0 and t τ , respectively.For the specific case under consideration, we can express the state in terms of shift operators ,j π e − p m τ 0 ∂q e Bω i τ ∂p e Bω i⊕1 τ ∂p j u i (q, p, t = 0), (49a) e Bω i τ ∂p e Bω i⊕1 τ ∂p j e Bω i τ ∂p u i⊕1 (q, p, t = 0), where in the first equation above, S ,j is a proper subset of the set of all permutations of the shift operators.Each element π ∈ S (1) ,j creates a different combinations of the + 2j shift operators for position and momentum, while preserving the relative order of the shift operators for the momentum.The same applies to the set S (2) ,j in the second equation, with the difference that in this case there are 2j + 1 shift for the momentum operator.
To better understand the above solution, let us consider the case in which the initial CQ state is δ(q)δ(p) |0 0|, so that at t = 0 only the level |0 is populated.Additionally, we fix ω 0 = −ω 1 = ω > 0, so that opposite shifts in momentum cancel each others, see Fig. 4. At time t we find that both quantum levels are populated, and in partiular we have e Bωτ ∂p e −Bωτ ∂p j δ(q)δ(p), (50a) By counting the number of shift operators in momentum, we can see that the population u 0 spreads in position q while keeping the momentum p fixed at 0. On the other hand, the population u 1 has an odd number of momentum shifts, and therefore its momentum is fixed at −Bωτ .We can also use the unravelling technique to obtain a numerical solution for the CQ dynamics of Eq. (46).In this equation, we see that the CQ state can evolve either continuously, or through a jump in both the quantum and classical degrees of freedom.When a jump occur, the quantum state is projected in either |0 or |1 , and this change is accompanied by a positive or negative shift in momentum, respectively.Furthermore, due to the form of the Lindblad operators, we have that a CQ state whose quantum degree of freedom is described by |0 can only jump when it is hit by the operator L 0 , and vice versa for |1 .In our simulation, we require these jumps to be associated with an opposite change in momentum.In this case, an initial CQ state with well-defined momentum (for instance the one used in Fig. 4) can spread to at most 3 different values of momentum.

Diffusion in phase-space
We can additionally study the spreading, due to the interaction between classical and quantum degrees of freedom, of the position of the particle as time passes by.The numerical simulation shows that at later times, the populations u 0 (q, p, t) and u 1 (q, p, t) divide into two Gaussian distributions, one with zero momentum, the other with non-zero momentum, see Fig. 4. Here, we consider the same scenario as in the previous section, where the initial CQ state is (q, p) = δ(q)δ(p) |0 0|.As we noticed before, this state evolves into two different ensembles, one associated with the quantum level |0 , and the other with the quantum level |1 .In particular, the population u 0 spreads in position while keeping the momentum fixed at p = 0, while the population u 1 spreads in position with a fixed momentum p = −ωBτ .In Fig. 5, we show the average position and standard deviation of the two populations as a function of time.We obtain these values by fitting the numerical data, and by the following analytical considerations.
Let us consider the evolution in phase space of the initial state under consideration.At each time step, this state has a probability P jump = δt τ of jumping in momentum, and a probability 1 − P jump of jumping in position.As we noticed before, given the initial state under consideration, the momentum of the CQ state at time t can be either 0 or p min = −Bωτ .When the momentum is non-zero, a position jump modifies the position of the state by ∆q = p min δt = −Bωτ δt.We can then compute the average change in position between two jumps, when the momentum is p min .The (normalised) probability that n 0.000 0.025 0.050 0.075 0.100 t (s) 0.0000 0.0002 0.0004 q (m) u 0 (q, p = 0.0 kg • m • s −1 , t) q fit σ q fit q analytic σ q analytic 0.000 0.025 0.050 0.075 0.100 t (s) 0.0000 0.0002 0.0004 q fit σ q fit q analytic σ q analytic Figure 5: The mean value q (or better, its absolute value) and standard deviation σ q of populations u 0 and u 1 .For the numerical simulation we choose the frequencies , the mass m = 1 kg, and a jump rate τ = 10 −2 s.The time step used in the simulation is δt = 10 −4 s.The solid lines are obtained by fitting the numerical data with a normal distribution, while the dashed lines are obtained using analytical considerations discussed in the main text.It is interesting to notice that, while σ q is roughly the same for the two populations, the mean position is shifted.This is because the distribution of u 1 is the first to be populated, since only one jump in momentum is needed when starting from the initial quantum state |0 0|.On the other side, the distribution of u 0 is populated (for values of q = 0) only after two jumps in momentum are preformed.Since the rate of jump is given by τ , we find that the mean position of the two populations is actually shifted (in time) by this parameter.
time steps separate two jumps is given by Using the above distribution, we can compute the average distance travelled by the state in between two momentum jumps, which is given by After n time steps, the average number of jumps is N jump = nP jump .Furthermore, the position of the state can only change after an odd number of jumps, since only in that case the momentum is non-zero.As a result, we find that the average position of the state at time t is where we define t = n δt.It is worth noting that the above result only applies to the average position of the population u 1 .In this case, since the initial quantum state is |0 , a single jump in momentum is required for the CQ state to populate |1 .However, the CQ state starts to populate |0 (with a position q = 0), only after two jumps in momentum are performed.Since the rate of jumping is given by the parameter τ , we have that the average position of the population u 0 is roughly delayed by this amount of time, and it is therefore equal to q (t − τ ), see Fig. 5.
To compute the variance of the position for the populations, we notice that the number of jumps N jump is distributed according to a binomial distribution, and its variance is given by σ 2 jump = nP jump (1 − P jump ).Then, keeping the other terms constant in Eq. (53), we have As we noticed before, the CQ state starts in position q = 0 with momentum p = 0, and roughly starts moving only after t = τ , that is, the time interval before a momentum jumps occurs (on average).As a result, the standard deviation of the populations needs to be modified so as to account for this delay, and we find that for both u 0 and u 1 the standard deviation is given by σ q (t − τ ), at least for t τ .As a final remark we notice that also in this model, like in the previous one, the average energy, as computed by averaging the operator H(z) = H C (z) + H I (z) over the classical and quantum degrees of freedom, is not conserved.Perhaps surprisingly, though, the energy increases from the initial time until it saturates at a positive value, see Fig. 6.This is due to the form of the Lindblad operators considered, which flip a |0 state into a |1 and vice versa, and thus create, after an initial transient, an equilibrium between four possible states in the classical phasespace, highlighted in the last panel of Fig. 4. In this equilibrium, we find that the kinetic energy receives a contribution by those states associated with non-zero momentum, while the potential energy equilibrates since the contribution from the term u 0 is proportional to that of u 1 , but with opposite sign (see the first panel of Fig. 5, and the above discussion on the expected value of the position).0.000 0.025 0.050 0.075 0.100 t (s) 0.0000 0.0001 0.0002 .As in the first model in the previous section, the total energy of the system is not conserved.However, in this case the total energy E tot = H C + H I increases during the initial part of the evolution, until it reaches the maximum value 1 2 τ 2 .

Quantum harmonic oscillator in a classical potential
In this section, we will consider the dynamics of coupling two degrees of freedom, one quantum and one classical.It is a particle whose internal degree of freedom is a quantum harmonic oscillator, coupled to its classical position and momentum.The model we focus on here demonstrates some basic features of the classical-quantum dynamics, though it describes a somewhat artificial situation where raising the energy level of the harmonic oscillator kicks the momentum towards the left, and lowering the energy kicks the momentum towards right.One quickly sees that energy will not be conserved.
The simpler model we explore in this section has one important and universal feature: starting from a state that is a delta function in momentum and position, and a superposition of two states in the harmonic oscillator, as expected, the numerical simulations show the initially pure state going through a process of decoherence, with the parts of the state associated to different energy levels becoming more and more distinguishable in the phase space as time goes on.We numerically simulated the dynamics and analytically found the decoherence rate, for the large n approximation (where n marks the energy levels of the harmonic oscillator).
We will begin with a simple model with interaction Hamiltonian, depending on both lowering and raising operators, which are quantum operators acting on the internal quantum degree of freedom, and the classical position q and momentum p of the particle carrying it: This simple model is analogous to the one used for the spin-1 2 particle.We consider a model with two Lindblad operators, one of which is proportional to the creation operator and the other is proportional to the annihilation operator: As a result, the functions h α (q, p) are defined as In the following, we focus on the case in which B(q) = qB, where B is a constant, and for simplicity we set the coupling constant ω 0 = −ω 1 = 1 s −1 .Notice that, for this choice of constants, the interaction Hamiltonian reduces to H I (q, p) = qB I, and the commutator in Eq. (28) drops out.We will discuss more general models, including two coupled oscillators in Sec. 5.
The resulting CQ master equation, where we additionally account for the free-particle classical Hamiltonian H C = p 2 2m , is given by where τ is the rate of jump in both the classical and quantum degrees of freedom.The constants γ ↑ , γ ↓ determine the rate of damping (γ ↓ ) of the quantum oscillator compared with the rate of pumping (γ ↑ ).For simplicity we henceforth take γ ↑ = γ ↓ = 1 although from a physical point of view, one might want to have the damping term be larger, to drive the quantum oscillator to its ground state.In a field theory, taking the damping term to be larger than the anti-damping term, can result in the theory violating causality [60], a state of affairs which the authors of [60] tune down by allowing some degree of energy non-conservation.In this model, when the energy of the particle is increased with the raising operator a † , its momentum is lowered, and vice versa, when the energy is lowered with the annihilation operator a, its momentum is increased.It is worth noting that the momentum is here one dimensional; thus, increasing the momentum does not necessarily mean increasing also its absolute value, and therefore, for the free-particle Hamiltonian H C , we cannot conserve the energy in this model.An example of the dynamics generated by this master equation is shown in Fig. 7.
Because the CQ master equation introduces a new contribution to decoherence, coming from the interaction between the classical and quantum degrees of freedom, we now aim to study the decoherence rate of the harmonic oscillator.In order to analytically find this rate, we have to compute the time evolution of the CQ state; to simplify the analysis, however, we will neglect the contribution of the Poisson bracket with the classical Hamiltonian H C , and only focus on the jumping dynamics, which is here responsible for the decoherence of the system.First, we write out the CQ state (q, p, t) in terms of the Fock basis, We can additionally perform the Fourier transform of the above state with respect to the momentum p.Thus, we get the following equation This step simplifies our analysis, since in this way we remove from Eq. ( 58) the nonlocality in momentum.However, at the end of our analysis we will have to perform an inverse-Fourier transform to recover the state in the momentum domain.
If we sandwich the above equation with two energy eigenstates n| • |m , we obtain the following equation for the elements u n,m (q, x, t) of the CQ state, The fact that the matrix element (n, m) interacts only with itself, and its nearest neighbours (on a line parallel to the matrix diagonal) (n − 1, m − 1) and (n + 1, m + 1) will prove to be convenient.Indeed, this implies that if we are interested in the problem where the initial state of the CQ density matrix is a superposition of two energy eigenstates, then we only need to solve three differential vector equations, where the vector takes values along the three diagonal lines involved; the main diagonal, among whose elements we find (n 1 , n 1 ) and (n 2 , n 2 ), the second line parallel to the diagonal marked by the matrix element (n 1 , n 2 ), and the third one marked by (n 2 , n 1 ).

Solving the differential equations for the CQ harmonic oscillator
We now attempt to solve the differential equation with the previously mentioned initial conditions, where we know that the three diagonals along which we have to solve are the ones marked with the matrix elements (n 2 , n 1 ), (n 1 , n 1 ) and (n 1 , n 2 ).For example, we can focus on the differential equation of the third diagonal, Assuming that n 1 < n 2 , we can now rearrange the functions u n 1 +k,n 2 +k (q, x, t), where the { e k } k are unit vectors, and we suppressed the time and phase space dependence for simplicity.The goal is then to find the solution of the equation where M is a tridiagonal matrix, whose entries are determined by Eq. (62).The size of this matrix is infinite, but we can obtain an approximate solution to the above equation by considering only part of the elements in the vector v.In particular, we require the index k to take values in an interval − N 2 , N 2 , for a given integer N .This approximation is acceptable for short enough times, when the coherence u n 1 ,n 2 has not yet spread too much.By imposing these lower and upper bounds to the k index, we make the system effectively finite dimensional, and therefore the finite-dimensional matrix M can be diagonalised.In the following, we refer to the i-th eigenvector of this matrix as w i , with eigenvalue λ i .
The equation regulating the dynamics of the eigenvectors w i is then given by and the solution is given by w i (t) = e λ i t w i (0).The same procedure can be followed for the main diagonal, marked with the matrix elements (n 1 , n 1 ) and (n 2 , n 2 ), so as to obtain the evolution for the populations.The diagonalisation of the matrix M can be performed numerically, so as to obtain the decay rates for the coherences of the (Fourier transformed) CQ state.However, obtaining the evolution of the actual CQ state ρ(q, p, t) would require us to Fourier transform back to the momentum domain, which can be computational expensive, given that such transformation should be performed on every element u n,m (q, x, t) of the CQ state.For this reason, we preferred to make use of the unravelling technique to obtain the evolution of the CQ harmonic oscillator, and to compute the decoherence rate.

Decoherence rate in the large n approximation
We now compute an analytical expression for the decoherence rate in the case in which the Fock numbers n 1 and n 2 are large compared to the range of the k index.In this situation, we can approximate the matrix elements along the diagonals of M to be equal, so as to obtain a Toeplitz matrix.Concretely, if the k index takes value between − N 2 , N 2 , where N n 1 , n 2 is an even number, then the matrix M can be approximated as where we define e iφ = e iBτ x for convenience.We can now find the eigenvalues λ m of the above Toeplitz matrix, which are given by [48] The eigenvectors of the matrix { ω m } N m=1 are instead given by where B m,r = 2 N +1 e −irφ sin m r π N +1 is the unitary matrix diagonalizing M.Under the dynamics described by Eq. (64), the eigenvectors evolve as ω m (t) = e λm t ω m (0).Thus, by re-writing the vector v, defined in Eq. (63), in terms of the eigenvectors of M, we can compute the approximate time evolution of the coherence of the state, which is given by r,m e λm t B m, , = 1, . . ., N.
(69) We can now specialise the above approximate solution to the initial condition we introduced in Eq. (61).In this case, the sole contribution to the coherence at time t = 0 is u n 1 ,n 2 (q, x, 0) = 1  2 δ(q)δ(p) = 1 2 δ(q) dxe ipx , and its time evolution is given by m π e where the dominant rate of decoherence is given by Γ = n 1 +n 2 τ , with corrections of the order √ n 1 n 2 .The above analytical result can be compared to the numerical one obtained with the unravelling technique, see Fig. 8.When the difference between n 1 and n 2 is big, the coherence reaches lower values than in the case in which n 1 and n 2 are close together; this behaviour of the coherence is analogous to the case in which the dynamics is fully quantum, and the system decoheres due to the coupling with an external bath, rather than to the interaction with the classical degrees of freedom.
Notice that, in order to obtain the decoherence as a function of phase space variables, we have to perform the inverse Fourier transform on the coherence given in Eq. (70).However, since the time evolution of u n 1 ,n 2 is constant in the x-space, we find that in pspace the solution picks up a delta function δ(p), and the decoherence rate is not modified.The presence of a delta function for the momentum is consistent with the dynamics we are considering.Indeed, from Eq. (58) it is easy to see that the coherence u n 1 ,n 2 is non-zero for p = 0 only, while for any k = 0 the coherence u n 1 +k,n 2 +k is non-zero for p = −kBτ .Likewise, the position dependence is described by δ(q), since the master equation we are considering does not include the classical Hamiltonian H C = p 2 2m .It is worth noting that the evolution of the coherence induced by the CQ dynamics described in Eq. ( 58) is analogous to the evolution induced by a fully quantum master equation where the phase space degrees of freedom are neglected.The decoherence can however be interpreted in two different ways; for the CQ dynamics, this is due to the interaction with the classical degrees of freedom, and the fact that (quantum) information is leaking in the classical phase space.For the fully quantum dynamics, instead, the decoherence is due to the interaction with an external environment.The coherence evolution for these two settings can still be distinguished if we introduce the classical Hamiltonian H C in the CQ dynamics, since this term provides a non-trivial dependence on the position q in phase space.As a final remark, we notice that coherence in this model is "moving" from the initial superposition between |n 1 and |n 2 to superposition of nearby Fock states.This is shown in the left and centre plots of Fig. 8; the coherence of the initial state |n 1 + |n 2 monotonically decreases, while the coherence of the superposition |n 1 ± 1 + |n 2 ± 1 increases at early times and then start decreasing as well.

Uniqueness of unravelling
We have here considered general dynamics which consistently couples a quantum system with classical degrees of freedom.The combined system is described by a hybrid-density which gives the probability of the system to be at any point in phase space and if so, what the density matrix of the quantum system is.We introduced an unravelling approach to solve the equations of motion.In this approach, both the classical and quantum systems evolve along a particular trajectory with some probability, and when we average over these trajectories, we recover the dynamics for the hybrid-density.Unlike the purely quantum case, the trajectories considered here can lead to objectively certain state transformations, in which the quantum system transitions to a particular pure state, conditioned on the classical degrees of freedom.This happens when the Lindblad operators L α are uniquely determined by the classical transition W α (z|z − ∆) rate, so that the observation of a particular classical transition unambiguously informs us about the quantum transition.Since the classical system has a definite trajectory, the trajectory of the quantum system is also definite.This is contrast to the purely quantum case, where the decomposition of a density matrix into an ensemble of pure states is not generally unique.A particularly simple and illustrative example is the master equation where, if the ∆ α are different for each α, observing a sudden jump in momentum by an amount ∆ α implies that the quantum state has changed from |ψ to a suitably normalized L α |ψ .The classical system encodes the trajectory of the quantum system.Furthermore, for the class of hybrid master equations we considered, there is a very real sense in which the quantum system jumps from being in a superposition of several states, to being in one of those states.Ironically, this is in contrast to spontaneous collapse models [32,57,33] where collapse is meant to be caused by an unobservable field.In these theories, as in ordinary quantum theory, the density matrix of the quantum system fully describes the system, and so, one cannot distinguish between decoherence (which we think of as the gradual decay of off-diagonal elements of the density matrix in a chosen basis), and the case where the wavefunction has a probability per unit time of being projected to a particular states, such that on average, off-diagonal matrix elements of the density matrix decay with time.In both cases, the density matrix of the quantum system is the same.In contrast, in the theory studied here, by monitoring the classical degree of freedom one can determine that a sudden change has occurred.
The system which is treated classically could be merely a large quantum system that is effectively a classical system such as an environment or a measuring apparatus.But it could also be fundamentally classical, for example, space-time, or it could be a system whose classical nature is somehow assumed, such as the experimenter herself.This provides a potential solution to the measurement or "reality" problem of quantum theory.Namely, quantum systems in superposition collapse to a particular state in that superposition because they interact with a system which is classical.Which set of possible states they can collapse to (sometimes called the pointer basis) is determined by the interaction.
Turning towards more pedestrian phenomena, we have found that the toy models we studied exhibited a trade-off between the amount of classical diffusion that occurred in phase space and the amount of quantum decoherence or noise that occurred in the quantum system.To minimise the disturbance on the quantum state, there needs to be a lot of diffusion in the classical system.In Ref. [54] we will formally prove that this is indeed the case.

Discussion on energy non-conservation
We also found as in [4], that because the Hamiltonian was not the generator of timetranslations, it typically is not conserved.We saw that this effect does not have to be large, and there are a number of ways to keep it under control so that we can have stable solutions.Violations of energy conservation present a possible experimental signature of fundamental hybrid dynamics, but also feature in spontaneous collapse models, and gravitationally induced decoherence [22,58].A full study on the issue of energy conservation is beyond the scope of the current work, but in the following we present a brief discussion of the problem, at leading order.
It is instructive to recall the analogous issue in classical mechanics.If one has diffusion as in the Fokker-Planck equation, then the magnitude of the momentum will gradually drift and increase, resulting in energy increase.The introduction of a friction term prevents this, by causing the system to settle down to an equilibrium state.The Fokker-Planck equation for the phase space density ρ(q, p) is of the form ∂ ρ(q, p) ∂t = {H C (q, p), ρ} + ∂ 2 D(q, p)ρ(q, p) ∂p 2 + ∂ ∂p η(q, p)ρ(q, p) , (72) with the friction coefficient η := γ ∂H C ∂p .This is just the theory of Brownian motion as developed by Einstein and Smoluchowski, with a thermal stationary distribution given by where Z is the partition function and the inverse temperature is By adding a friction term, we can stop the diffusion in the classical degrees of freedom.
For the type of hybrid model presented in Sec.3.1, with Lindblad operators that are a basis of orthogonal projectors L α = |α α| onto energy eigenstates of H Q , and h α = ωBq (to take a simple example) we may add a friction term such as to the right hand side of the evolution law of Eq. (28).The friction coefficient γ can be chosen to depend on q.To first order in τ , the diffusion term is given by Eq. (39), and so from Eq. (74), we expect that to first order, the classical system will equilibriate to a thermal state at inverse temperature β = γ(q)/(Bω α ) 2 τ .Because the Lindblad operators project onto eigenstates of H Q , the quantum system also equilibrates by decohering into one of it's energy eigenstates.
The case where the Lindblad operations cause noise is more interesting, since then one has diffusion not only in phase space but also in the Hilbert space.What's more, the amount of diffusion in phase space can depend on the quantum state; for example, for a quantum harmonic oscillator the diffusion can grow with n, the energy level of the oscillator.This means that the diffusion term can be arbitrarily large, and overcome any classical friction term.We would thus like the friction term to also increase with n so that the classical system can reach a steady state.To explore this, let us consider a more sophisticated model than Eq.(72), where a classical harmonic oscillator with Hamiltonian H C = ω c (p 2 + q 2 )/2 is coupled to a quantum harmonic oscillators with free Hamiltonian H Q = ω Q a † a; the Lindblad operators are now L α = a, the annihilation operator.A natural coupling between the two oscillators would be to take the interaction Hamiltonian to be ), but it is hard to make the hybrid dynamics of such a coupling completely positive.An alternative, is to take the quantum part to be positive definite (while leaving the classical part arbitrary), e.g.H I = B(q)Q 2 , which is still a local coupling.Let us be slightly more general, and consider coupling to Q 2 + κP 2 with P = i √ 2 (a † − a) and κ ≥ 0. It will prove convenient to use [a, a † ] = 1 to rewrite this as where we have dropped the purely classical terms from the interaction.A master equation for such a theory is then ∂ (q, p) ∂t = {H c (q, p), (q, p)} + F( (q, p)) − iω Q a † a, (q, p) − iB(q) Q 2 + kP 2 , (q, p, t) where we are still to specify the transition rates W ↑ and W ↓ , associated to the pumping and dumping of the quantum oscillator, respectively, and the classical friction term F. The last term in the first line is the influence of the classical system on the harmonic oscillator, while the second and third lines is the influence of the quantum system on the classical one.In the following, we will require the expansion of W ↑ and W ↓ to contain a term proportional to τ B (q) ∂ (q,p) ∂p , so as to balance the energy increase of the classical system due to the interaction with the quantum one.
On the other hand, at leading order in ∆ the second and third lines give rise to a Lindblad equation with a damping and pumping term, and so when W ↑ (z) ≥ W ↓ (z), the state of the harmonic oscillator will continue to increase in energy.The quantum analogue of adding a friction term is to take the damping term to be strictly larger than the pumping term, and if the classical system were to remain at rest, the quantum oscillator will thermalise to an inverse temperature proportional to log W ↓ (z)/W ↑ (z), at which point one can verify that D( ) = 0. We can thus have the system reach equilibrium by taking log W ↓ (z)/W ↑ (z) to be strictly positive and adding a purely classical friction term F( ).
If on the other hand, we take B(q) and the rates W ↓ (z) and W ↑ (z) such that the classical system equilibriates in a region of phase space with sufficient probability of having W ↓ (z) ≤ W ↑ (z), then energy can be pumped into the harmonic oscillator indefinitely.In such a case, we would like to control this, so that it happens slowly.To explore this, let us consider the specific realisations, where the differential operator ∂ p (X•) acts over an hybrid operator ρ(q, p) as ∂ p (X ρ(q, p)).
The two realisations are positive as long as h ↑ and h ↓ is everywhere positive, and the master equation will preserve normalisation with W ↑ (q, p) = h ↑ (q)/τ and W ↓ (q, p) = h ↓ (q)/τ .Let us now choose X ↓ = ∆ ↓ (q) + η ↓ , X ↑ = ∆ ↑ (q) + η ↑ , where η ↑ and η ↓ are friction terms as in Eq. (72).Then, at zeroth order in τ we have the Lindblad equation (78) with temperature proportional to log h ↓ (q)/h ↑ (q).However, even if we don't have h ↓ (q) ≥ h ↑ (q), the rate at which energy is pumped into the oscillator is which is easily calculated from the master equation.This term can potentially be made small in the phase space region where the classical system equilibriates.The next order terms give, where we take the last two terms to be of roughly the same order as the previous diffusion terms.
The κ = 1 case is a natural one to consider, since it is reminiscent of the master equation of general relativity coupled to a scalar field [51].The classical oscillator couples to the energy of the quantum one, as the classical gravitational field couples to the scalar field expanded in terms of momentum modes.Likewise, analogy with the gravitational case suggests taking the pure commutator term of Eq. (77) since gravity is always coupling to the energy.In this case, we can take so that the first two terms in Eq. (81) will give the Poisson bracket {B(q) 1 2 Q 2 + 1 2 P 2 , } in the classical limit.When we put this all together, the master equation (77) expands as ∂ (q, p) ∂t = {H c (q, p), } + F( ) − iB(q) ω a † a, (q, p, t) In analogy with the field theory case, we would like to take h ↓ (q) = h ↑ (q) = h(q), since in that context, the theory is local with the consequence that the oscillator will not equilibrate and will instead increase in energy [60].This can be slowed down by any desired amount by taking h(q)/τ to be small.However, anticipating the discussion in Ref. [54], slowing the energy increase down results in a greater amount of classical diffusion.In particular, the requirement that Eq. (82) be satisfied, which in this case is encodes this trade-off.The right hand side governs the rate of evolution of p and is set by the dynamics, and the left hand side then requires that a realisation of this form have a trade-off between h(q) and ∆(q).Let us consider the simplest example of B(q) = Gq, where G is a positive constant, so that the classical system acts as if it is subject to a constant force F ≈ ωGn when the quantum system is in the Fock state |n , with n large.The rate at which energy is pumped into the quantum harmonic oscillator is δH Q = ωh(q)/τ , as shown in Eq. (80).Recalling the Fokker-Planck equation, the effective diffusion coefficient D is the term in front of 1 2 ∂ 2 ∂p 2 in Eq. (83), going as D ≈ F τ ∆(q)/2.Then, Eq. (84) gives δH Q D ≈ 1 4n ωF 2 so that a small amount of energy being pumped into the oscillator requires large diffusion in the classical system in relation to the total back-reaction exerted by the quantum system on the classical one.As n increases the relative trade-off, when the force F is held fixed, becomes less pronounced.
There appears to be significant freedom in the choice of h(q) and ∆(q), the main requirement other than Eq.(84) being that h(q) be positive.One can even take h(q) = |B (q)| and ∆(q) proportional to sign(B (q)).If on the other hand, ∆(q) is large, we can then get energy being pumped into the classical system, since large amounts of diffusion in momentum will increase the energy of the classical system.This can be controlled by increasing the friction term, to keep the momentum low.However, since the diffusion term goes like D ∝ nωGτ ∆(q) and thus increases with n, the friction term would also need to scale with n or risk getting overwhelmed by the diffusion.For this reason, the purely classical friction term F( ) will not be sufficient, while the hybrid friction terms in Eq. (83), should allow us to compensate the increase.This is seen explicitly by tracing out the oscillator, resulting in dynamics for the classical phase space density ρ n conditioned on the the oscillator being in state n where we have taken η ↑ = η ↓ = η for simplicity.To this order, this is just the Fokker-Planck equation, and if η = γ ∂Hc ∂p , then we expect the effective diffusion coefficient to be D ≈ τ ∆(q)B (q)ωn, and the effective friction coefficient to be γ ≈ nh(q)γ(q).Eq. (74) would then give the effective inverse "temperature" that the classical system is held at, although here, since it could depend on q, it is not a temperature but still defines the equilibrium state β = h(g)γ(q)/τ ∆(q)B (q)ω, which is now independent of n.Note that this system could act as if the classical system is being held at one temperature while the oscillator is held at a different temperature.Were we to have access to the environment that each system appears immersed in, we would presumably see heat flow from one environment to the other.Understanding how this would work in detail, either analytically or numerically would be an interesting research direction to pursue.One would like to better understand the various trade-offs involved in systems like this, and what realisations lead to equilibrium or steady states, or states which only change slowly.sical degrees of freedom.The classical system is here described by the ensemble ρ(q, p, t), and its evolution is given by ∂ρ(q, p, t) ∂t = 1 τ 0 e τ 0 {H C (q,p),•} − 1 ρ(q, p, t), (87) where the finite parameter τ 0 > 0 can be understood as the rate of jumping in phase space, due to the effect of the shift operator e τ {H C (q,p),•} .It is easy to show that the solution of Eq. (87) is given by where P (k) is the Poisson distribution which carries the time dependence of the solution, and the ensemble ρ(kτ 0 ) is obtained from the initial ensemble ρ(q, p, t = 0) by applying the operator k times, ρ(kτ 0 ) = e kτ 0 {H 0 (q,p),•} ρ(q, p, 0).
It is worth noting that each ensemble ρ(kτ 0 ) is the solution of the Liouville equation at time t = kτ 0 , so that we can interpret it as the state of the system after k jumps in phase space.Thus, the solution of Eq. (87) can be understood as a mixture of different ensembles, where each of these ensemble are the solution of Eq. (91) at discrete times kτ 0 , for k ∈ N, and the probability of making k jumps is given by the Poisson distribution whose mean value is t τ 0 .We can additionally use the simple master equation here introduced to better understand the role of the diffusion term discussed in Sec. 2, and more generally the role of the infinite number of corrections we obtain when the hybrid master equation is expanded with respect to τ 0 .Let us first notice that Eq. (87) can be re-express as where ρ(q, p, t + τ 0 ) is the solution of the Liouville equation at time t + τ 0 , as we show in Eq. (90).Thus, the solution of the master equation is given by an ensemble whose instantaneous derivative at time t depends on the solution of the Liouville equation at time t + τ 0 , and more in general on the phase-space trajectory imposed by the Hamiltonian H C (q, p) through the Liouville equation, see Fig. 9.As a result, the instantaneous derivative of the solution of the master equation cannot be given by the sole Poisson bracket of the Hamiltonian, but needs to be suitably corrected.These corrections are reflected in the infinite expansion of the hybrid dynamics, see Eq. (10) in the main text.

A.2 Solution to the stochastic CQ dynamics with diagonal Lindblad operators
We can now provide a solution for Eq.(28) in a simple case, which nevertheless turns out to be useful for obtaining an analytical solution in Sec.3.1.Specifically, we study the case The phase-space of a classical system with Hamiltonian H C (q, p).The blue line represent the continuous trajectory of the ensemble ρ(q, p, t), which initially is described by ρ(q, p, 0) = δ(q − q)δ(p − p).This trajectory is obtained through the Liouville equation, see Eq. (91).Each point on the trajectory corresponds to the position of the ensemble at discrete times t = kτ 0 , k ∈ N. The solution of Eq. ( 87) is a mixture of the ensembles at these discrete times, weighted by a Poisson distribution with mean value t τ0 .This equation can be understood as a constraint on the instantaneous derivative of the ensemble, as shown in Eq. (92).In particular, the derivative at time t needs to be equal to ρ(q,p,t+τ0)−ρ(q,p,t) τ0 , for a finite value of τ > 0.
in which a single population u(q, p, t) of the CQ state is considered, and only two shifts are considered, one associated with the classical Hamiltonian H C (q, p), and one associated with the quantum coupling H I (q, p).The equation we consider is then ∂u(q, p, t) ∂t = 1 τ 0 e τ 0 {H C (q,p),•} u(q, p, t) − u(q, p, t) + 1 τ 1 e τ 1 {H I (q,p),•} u(q, p, t) − u(q, p, t) , (93) where the action of the two Hamiltonians can either commute or not.
If the action of the two Hamiltonians commutes, {H I (q, p), H C (q, p)} = 0, then the solution of the above equation can be easily obtained, and it is given by where k is the number of jumps associated with H C , and n is the number of jumps associated with H I .Both P 0 (k) and P 1 (n) are Poission distributions, defined as and the population u(n, k) is obtained from the initial one u(q, p, t) by applying k classical jumps and n quantum jumps, that is, It is easy to check that this is the solution of Eq. (93).
When the action of the two Hamiltonians does not commute, i.e., {H I (q, p), H C (q, p)} = 0, the solution of the stochastic CQ equation is still in the form of Eq. (94), with P 0 (k) and P 1 (n) Poisson distributions.However, the population u(n, k) is now different, since the order of the shifts is now important, and we write it as   e τ 0 {H C (q,p),•} . . .e τ 0 {H C (q,p),•} k e τ 1 {H I (q,p),•} . . .e τ 1 {H I (q,p),•} n   u(q, p, 0), (97) where π i is a permutation on the n + k shifts which create a different combination, and i = 1, . . ., n+k k .
A.2.1 Solution to the stochastic CQ dynamics for a qubit in linear potential Let us now use the results obtained in the previous section to study the analytical solution of the qubit in linear potential, where the Lindblad operators are diagonal, see Sec. 3.1.It is easy to see that the differential equations for the populations, Eq. (32a), are almost in the same form of the one we have studied in the previous section.Indeed, one can obtain these equations from the stochastic one in Eq. (93) by setting the classical Hamiltonian H C = p 2 2 m , the quantum coupling H I = q Bω, and sending τ 0 → 0. Notice that the action of the two Hamiltonians does not commute, since their Poisson bracket is not zero, and the two shifts act as follow on the population, e τ 0 {H C (q,p),•} u(q, p) = u(q − p m τ 0 , p), (98a) e τ 1 {H I (q,p),•} u(q, p) = u(q, p + Bωτ 1 ).
It is clear that shifts in momentum affect the subsequent shifts in position, since the latter depend on the value of the momentum.We can define the jump units for position and momentum as ∆p = Bωτ 1 , and ∆q = Bωτ 1 m τ 0 .We are now interested in expressing the population u(k, n) of Eq. (97) in terms of its position and momentum, rather than in terms of its classical and quantum jumps.We can write it as u(k, n) = dq dp P (q, p|k, n)u(q, p) (99) where the conditional probability distribution P (q, p|k, n) can be divided into two distributions, one for the position and one for the momentum, P (q, p|k, n) = P (q|k, n) P (p|n).( 100) By looking at the effect of the quantum jump on the state, see Eq. (98b), it is easy to show that the distribution of the momentum, conditioned on the number of quantum shifts n, is given by The conditional probability distribution for the position, instead, is less straightforward to obtain.
In order to get the mean value and variance of the distribution P (q|k, n), we express the problem in a different way.For a fixed value of classical and quantum jumps k and n, we consider all possible combinations of jumps, or histories.We represent one such history with a n + k-bit string x, where 0 is a position jump, and 1 is a momentum jump.The final position of the particle for a given history can then be computed using the following formula, where x is the -th bit of the string x.Eq. (102) provides a link between a history x and the final position q, and can therefore be used to get the conditional distribution we are interested in.In particular, we can make the simplifying assumption that each x is an independent random variable, x = 1 with p = n n+k , 0 with 1 − p = k n+k . (103) The mean value and variance of x is the same for all , and they are given by x = n n+k and σ 2 = nk (n+k) 2 , respectively.If we now use Eq.(102), together with the fact that each x is independent, we find that where in the second equation we have used the fact that, given two random variables X, Y and two real numbers a, b, the variance of a random variable Z = aX + bY is given by σ 2 Z = a 2 σ 2 X + b 2 σ 2 Y .We can now use the above results for estimating the mean value and variance of the marginal distributions for position and momentum.These distributions are given by, respectively, P (q) = ∞ n,k=0 P (q|k, n) P 0 (k) P 1 (n), (105a) The mean value and variance for momentum are given by the Poisson distribution P 1 (n), since P (p|n) = δ(p − n∆p), and we get p = ∆p t τ 1 = Bωt, (106a) To compute the mean value and variance of position, instead, we can make use of the law of total expectation and the law of total variance.In the limit of n, k 1, we have that σ 2 q k,n ≈ 1 3 nk(n + k), and after straightforward calculations we find q = 1 2 The results of Eqs.(106) and (107) are then compared, in Fig. 2, to the values obtained by numerical simulation.

A.3 Solution to the stochastic CQ dynamics with non-diagonal Lindblad operators
We now consider the case in which the Lindblad operators of Eq. ( 28) are non-diagonal, and therefore the evolution of each population depends on the others.We consider the case of a two-level quantum system, since this is the situation studied in the main text, see Sec. 3.2.In particular, we study the evolution of the populations of a CQ system with classical Hamiltonian H C (q, p) and interaction Hamiltonian H I (q, p) (whose action we assume does not commute), and with Lindbald operators connecting the two levels (see Eqs. (44) for an example).The stochastic equation for the population u i , where i = 0, 1, is given by ∂u i (q, p, t) ∂t = 1 τ e τ 0 {H C (q,p),•} u i (q, p, t) − u i (q, p, t) Q (q,p),• u i⊕1 (q, p, t) − u i (q, p, t) , (108 where the interaction Hamiltonian additionally depends on the level is acting on. In order to solve this equation, we re-write the time derivative as (u i (t + δt) − u i (t)) 1 δt , and we express Eq.(108) as where we suppressed the dependence on the phase-space variables for compactness.If we now group the terms in the above equation into the following operators, acting on u i and u i⊕1 respectively, we can express the evolution of the populations as u 0 (t + δt) = A u 0 (t) + B 0 u 1 (t), (110) By considering a sequence of times t = {k δt} k∈N , we can recursively solve the above equations, and we obtain where the coefficients (when k is even) are given by In the above equations, π m ∈ S k are permutations of the operators A's and B i 's such that each combination obtained is different from the others, and the relative order between the B i 's is not modified.
If we now replace the operators A and B i , and we send k → ∞ and δt → 0 while asking for k δt = t, we find that the solution of Eq. (108) is P 0 ( ) P 1 (2j + 1) u i⊕1 ( , j ), (114) where both P 0 and P 1 are Poisson distributions with mean value t τ 0 and t τ 1 respectively, and  u i (q, p, t = 0), (115a) In the above solution, the probability distributions P 0 and P 1 carry the time dependence, while the functions u i ( , j) and u i⊕1 ( , j) contain the phase-space information.In order to find a more explicit solution, one needs to consider specific Hamiltonian operators H C (q, p) and H I (q, p), as we do in the main text in Sec.3.2.

B Unravelling code for CQ dynamics
The unravelling code is tailored to solve a specific model of CQ dynamics, whose master equation is given in Eq. (28) for a rate of classical jump τ 0 → 0, that is, ∂ (z, t) ∂t = − i [H I (z), (z, t)] + {H C (z), (z, t)} where H I (z) = α h α (z) L † α L α is the interaction Hamiltonian between the classical and quantum degrees of freedom.For simplicity, we consider the case in which z = (q, p), that is, we consider a single system in phase space.The unravelling of this master equation can be obtained following the steps shown in Sec.2.1 of the main text.However, here we will modify the updating rule of Eq. ( 16) to include a continuous evolution of the classical degrees of freedom together with that of the quantum degrees of freedom.
As stressed in the main text, the unravelling code evolves the CQ pure state, see Eq. ( 14), which we represent with the tuple (|φ , q, p, t).The first element of the tuple, |φ ∈ C d is the pure state of a qudit, while q ∈ R is the position and p ∈ R is the momentum of the system, and t ∈ R + is time.To update the state, we first sample from a uniform distribution over the interval [0, 1], and check if the outcome p is lower or equal that If this is the case, that is, p ≤ p 0 , we apply the continuous update to the CQ state, and we obtain the following tuple, where N 0 is the normalisation of the updated quantum state, and H eff = H I (z)− i 2τ α L † α L α .If p > p 0 , instead, we randomly choose an α-jump according to the distribution, In this case, we evolve the state using the following update rule, where N α is the normalisation of the new quantum state.
It is worth noting that the updating rule of Eq. (118) includes the continuous evolution of the classical degrees of freedom, as opposed to the case shown in the main text, Eq. ( 16).This addition in the updating rule is due to presence, in the above master equation, of the Poisson bracket between the classical Hamiltonian H C and the CQ state.In fact, it is easy to see that the first order in δt of the phase space shift operator coincides with the action of the Poisson bracket, since for any (analytic) function f over the phase space we have Using the above equation together with the fact that the unravelling provides a first order approximation of the evolution given by the master equation, one can show that the updating rule in this section correctly reproduces the dynamics of Eq. (116).
The code for the CQ unravelling can be found on GitHub.

Figure 2 :
Figure2: The mean value and standard deviation of the population u 0 (q, p, t) in phase space.Left.The mean value q and standard deviation σ q of the marginal distribution u 0 (q, t) (where we traced out over the momentum degree of freedom).The solid lines are the values obtained from the numerical simulation, whereas the dashed line are obtained from the analytical solution.Centre.The mean value p and standard deviation σ p of the marginal distribution u 0 (p, t).Right.The coherence |c| 2 of the quantum state located in q = 0 m and p = 0 kg • m • s −1 , as a function of time.It is worth noting that the numeric solution is obtained with the unravelling method, in which the position is updated during the continuous evolution step, see Eq. (118).Instead, the analytic solution is obtained for the case in which the position evolves through jumps, see Eq.(34).In order to match the two evolution, we have to ask the stochastic rate τ 0 to be of the order of the infinitesimal time steps δt used in the simulation.The mean value and variance of position and momentum are therefore computed using Eqs.(37) and (38) for a value of τ 0 = δt.

Figure 3 :
Figure 3: The different contributions to the total energy of the CQ dynamics, as a function of time.The solid lines are the numerical results, while the dashed ones are the analytical ones.The energy contribution given by the kinetic energy, E C = H C , grows faster than the one given by the potential energy, E I = H I .As a result, the total energy of the system, E tot = H , increases linearly as a function of time.

Figure 4 : 1 √ 2 (
Figure4: The evolution of the population u 0 (blue) and u 1 (green) under the CQ dynamics with non-diagonal Lindblad operators.The evolution of the populations is shown for t = 0 s, 0.04 s and 0.1 s.In this figure, we set the initial CQ state to be δ(q)δ(p) |+ , where|+ = 1 √ 2 (|0 + |1) is a coherent superposition of the eigenstates of the interaction Hamiltonian H I .Notice that, in the rest of the section, we instead focus on the case where the initial quantum state is |0 , which simplifies the analytical study of the dynamics; the qualitative behaviour of the phase-space evolution, however, is the same as the one shown in this figure.This interaction makes the quantum degrees of freedom stochastically collapse from |0 to |1 and vice versa, at a rate given by τ .When the quantum state is mapped as |0 → |1 , the classical momentum is increased by Bτ , while the transformation |1 → |0 is accompanied by a change in momentum equal to −Bτ .The simulation depicted here uses time steps of δt = 10 −4 s, and a jump rate τ = 10 −2 s.The interaction constant is B = 1 J • s • m −1 , the mass is m = 1 kg, and the frequencies are ω 0 = −ω 1 = ω = 1 s −1 .

Figure 6 :
Figure 6: The energy of the CQ state as a function of time, when the initial state is δ(q)δ(p) |0 0|.The different contributions to the energy are shown, such as the kinetic energy E C = H C and the potential energies E (i) I = H (i) I

Figure 7 :
Figure7: The evolution of the populations u n (q, p, t) of the harmonic oscillator in phase space,.The system is initially in a superposition over the state |n 1 = 12 and |n 2 = 16 , centered in the origin of the phase space.The evolution of the populations is shown for t = 0.01 s, 0.03 s and 0.5 s.The state is updated every δt = 2.5 • 10 −5 s, and the jump rate of the model is τ = 0.1 s.The constant B = 1 J • s • m −1 , and the mass of the harmonic oscillator is m = 1 kg.Top panels.The evolution of the population for n = n 1 is shown.Initially, this population is peaked around the origin of phase space, and it starts to diffuse in position due to the action of the rising and lowering operators.At later time, a second packet can be seen, with momentum p = 0.4 kg • m • s −1 .This packet originated from the population at n = n 2 , and reaches n = n 1 thanks to the repeated action of lowering operators.Central panels.The evolution of the population for n = n1+n2 2 is shown.The symmetric packets with negative and positive momentum (p = ±0.2kg • m • s −1 ) originated from n = n 1 and n = n 2 , respectively.Bottom panels.The evolution of the population for n = n 2 is shown.The general behaviour in phase space is analogous to that of the population for n = n 1 , but with a second packet with negative momentum.

|||Figure 8 :
Figure8: The evolution of the coherence for two different choices of initial state for the CQ system.In both cases, the system is evolved up to the time t = 2.5 • 10 −2 s, with time steps of δt = 2.5 • 10 −6 s, and a jumping rate τ = 2.5 • 10 −1 s.The interaction constant is B = 1 J • s • m −1 , and the mass is m = 1 kg.Left.The initial quantum state is an even superposition between |n 1 = 15 and |n 1 = 30 .The evolution of the quantum coherence |u m1,m2 | is shown for different values of m 1 and m 2 , and the classical degrees of freedom have been traced out.The numerical evolution is represented with the solid line, and the dashed line is the analytical solution given in Eq. (69), for N = 10.Centre.The initial quantum state is a even superposition between |n 1 = 22 and |n 1 = 23 , and the evolution of the coherence is shown for those m 1 and m 2 close to n 1 and n 2 .Again, the solid lines are obtained through the numerical simulation, using the unravelling code, and the dashed ones are the (approximate) analytical solutions.Right The evolution of the coherence |u n1,n2 | for the two choices of initial state.When the difference between n 1 and n 2 is big, the coherence reaches lower values than in the case in which n 1 and n 2 are close together; this behaviour of the coherence is analogous to the case in which the dynamics is fully quantum, and the system decoheres due to the coupling with an external bath, rather than to the interaction with the classical degrees of freedom.

Figure 9 :
Figure9: The phase-space of a classical system with Hamiltonian H C (q, p).The blue line represent the continuous trajectory of the ensemble ρ(q, p, t), which initially is described by ρ(q, p, 0) = δ(q − q)δ(p − p).This trajectory is obtained through the Liouville equation, see Eq. (91).Each point on the trajectory corresponds to the position of the ensemble at discrete times t = kτ 0 , k ∈ N. The solution of Eq. (87) is a mixture of the ensembles at these discrete times, weighted by a Poisson distribution with mean value t τ0 .This equation can be understood as a constraint on the instantaneous derivative of the ensemble, as shown in Eq. (92).In particular, the derivative at time t needs to be equal to