Multi-Qubit Joint Measurements in Circuit QED: Stochastic Master Equation Analysis

We derive a family of stochastic master equations describing homodyne measurement of multi-qubit diagonal observables in circuit quantum electrodynamics. In the regime where qubit decay can be neglected, our approach replaces the polaron-like transformation of previous work, which required a lengthy calculation for the physically interesting case of three qubits and two resonator modes. The technique introduced here makes this calculation straightforward and manifestly correct. Using this technique, we are able to show that registers larger than one qubit evolve under a non-Markovian master equation. We perform numerical simulations of the three-qubit, two-mode case from previous work, obtaining an average post-measurement state fidelity of $\sim 94\%$, limited by measurement-induced decoherence and dephasing.


Introduction
Circuit QED provides a promising avenue for the realization of quantum algorithms, with recent experiments showing increases in both coherence time and precision of control [1][2][3]. Quantum algorithms are thought to require error correction as a prerequisite [4], and quantum error correction requires non-demolition measurement of joint operators, most often Pauli operators of low weight [5][6][7][8]. This can be accomplished using an ancilla register which is prepared in a specific state, interacts with the encoded state, and is then measured (possibly destructively) [9][10][11][12]. In circuit QED, ancilla measurement has been accomplished by coupling the qubit to photons passing through a resonator, and observing the accrued phase using homodyne detection [13].
Recent work has begun to consider direct joint measurements in circuit QED, in which all qubits in the support of the measured operator are coupled to one or more internal resonator modes, using homodyne detection to observe an output mode, requiring no ancilla qubit. Difficulty in calculating the reduced qubit dynamics has restricted previous analysis of direct measurement schemes to systems containing two [14][15][16] or three [17,18] qubits. In this paper, we simplify this calculation, deriving reduced qubit dynamics for an arbitrary number of qubits and resonator modes. We then use the resulting stochastic master equation to extend the analysis of the three-qubit, two-mode scheme presented in [17,18].
The rest of this paper is organized as follows. We write the multi-qubit, multi-mode Lindbladian in Section 2, and incorporate it into a stochastic master equation corresponding to homodyne measurement of the output mode. In Sections 3 and 4 we determine the reduced equations of motion for the resonator and register states, respectively. Using these equations, we proceed to simulate multi-qubit measurement dynamics in section 5. We discuss what can be done to increase post-measurement state fidelity and conclude in section 6.

Parameters
We begin with a set of qubits Q, called the register, and a set of internal resonator modes C. For convenience, we denote as B the set of 0/1 assignments to the register (comprising n = |Q| bits). We consider a scenario in which the coupling between the internal modes and the input/output mode is described by the relation (Recall that the discrepancy in units is explained by comparing the waveguide commutation relation, a out (ω), a † out (ω ′ ) = δ(ω − ω ′ ), and the cavity commutation relation a k , a † k ′ = δ k,k ′ .) A single measurement tone is used, and we describe the dynamics in a frame rotating at the carrier frequency of this tone. We consider a model Hamiltonian for the coupled-qubit system that describes the essential aspects of the dispersive-coupling regime, that is, when the qubit frequencies are far detuned from the cavity frequencies: l∈Q Ω l 2 + k∈C χ k,l σ z,l + l∈Q,k∈C χ k,l σ z,l a † k a k + ǫ(t) k∈C √ κ k a k + a † k .
Here, = 1, ∆ k ω k −ω is the difference between the kth resonator mode frequency and the measurement tone frequency ω, a k a † k is the lowering (raising) operator on the kth resonator mode, χ k,l is the coupling frequency between the kth resonator mode and the lth qubit, Ω l + k∈C χ k,l is the Lamb-shifted qubit frequency (with Ω l being the bare frequency), and ǫ(t) is the time-dependent measurement tone amplitude.
We note that if we use a Schrieffer-Wolff analysis to derive the dispersive-coupling Hamiltonian from an underlying multi-qubit Jaynes-Cummings model [19], additional terms appear, which describe qubit-qubit and resonator-resonator couplings [20,21]. In Appendix A we give a full derivation of this expression. We give arguments for why these terms can be neglected (within a rotating-wave approximation) or incorporated into parity measurement schemes with straightforward modifications of the analysis described below. It is also interesting to note that, with a more general starting point provided by circuit Hamiltonians [22], couplings more general than the Jaynes-Cummings form appear, and some qubit-qubit and mode-mode coupling terms can be arranged to cancel, as explored in [22]. We proceed with the model Hamiltonian Eq. (2), as it permits a full and clear exploration of all the issues connected with parity measurement, without the inessential complicating features introduced by the additional coupling terms.
We assume that decoherence can be described using a Lindblad master equation, The noise sources we include in the model, as in previous work [13][14][15][16]18] are resonator photon loss, intrinsic dephasing/amplitude damping of the qubit states, and the Purcell effect [19]:ρ where κ k is the photon loss rate in the kth resonator mode, γ z,l is the intrinsic dephasing rate for the lth qubit, γ −,l is the intrinsic amplitude damping rate for the lth qubit, and λ k,l is an effective Purcell factor for the kth mode interacting with the lth qubit. Since, as emphasized by the derivation of Appendix A, the factors λ k,l can have either sign, it is perfectly possible to arrange these factors so that the k sum in the final term of Eq. (3) is zero for every l. In other words, an effective Purcell filter [23] can be created by taking advantage of flexibility provided by the multi-mode structure. Given the ongoing advances in qubit coherence, we believe it is also reasonable to ignore intrinsic qubit damping, i.e., we can set γ −,l = 0. Thus, from this point onward, we will ignore qubit damping effects (but we will retain qubit dephasing terms).
To model the evolution of the state ρ and measurement record j under homodyne measurement, we use the stochastic master equation [24]: Here, η ∈ [0, 1] is the quantum efficiency of the homodyne measurement, φ is the homodyne phase (which we set to 0, corresponding to measurement of the real part of the operator a out ) and dW is a Wiener increment (a normal variate with mean 0 and variance dt) [25].
As proved in Appendix B, a family of solutions to Eq. (3) or Eq. (4) can be expressed using pointer states [26,27]: Here, |α k,j (t) is a coherent state (an eigenstate of the lowering operator a k with eigenvalue α k,j (t)) corresponding to the bitstring j. This simplifies the numerical solution of the (deterministic or stochastic) master equation in the event that the initial state is coherent (the vacuum is such a state, with α k,j = 0 ∀k, j), by requiring only a fixed number of coherent state trajectories to be calculated, rather than the full, infinite-dimensional resonator state.
In the following sections, we simplify the equations of motion further, by deriving the equations of motion for the coherent state amplitudes {α k,j }, and incorporating these into the equation of motion for ρ Q , the register reduced state.

Resonator Equations of Motion
Equations of motion for the resonator mode lowering operator a k can be derived using input-output theory [28]: The corresponding equation forα k,j (t) can be found as in [18]: where j l is the value of the lth bit of j. The amplitude of the output of the system is (see Eq. (1)): In the following subsection, we place this system of linear, first-order ordinary differential equations (ODEs) into a canonical form which can be used for further analysis.

State-Space Representation
Systems of first-order linear ODEs with time-invariant coefficients can be represented using a vector x, called the state, a vector u(t), called the input, and a vector y(t), called the output. They can be written in a standard form: This is the popular state-space representation of linear time-invariant (LTI) systems [29,30].
For the evolution described in Eq. (10), and fixing the register state to a specific bitstring j, these matrices can be written explicitly: The response to a given input can be calculated in the Laplace domain using the transfer matrix G(s): In the following subsection, we use these matrices to tailor the values of the resonator parameters in order to perform a specific measurement.

Steady States
In general, the resonator system will exhibit a different response to a given input ǫ(t) for every distinct state of the register |j , j ∈ B. This is not useful for measurements of joint degrees of freedom, which should not distinguish between given subsets of B. Consider, as a first example, measurements of the Hamming weight h(j) (the number of qubits in the ground state). Given that the system in Eq. (13) only depends on j through the term l∈Q (−1) j l χ k,l in the A matrix, systems with identical Hamming weights have identical (A, B, C, D) matrices if all χ k,l are equal to a single constant χ: To perform a parity measurement, it is furthermore necessary for the responses to be close for all even h and all odd h, with the even and odd responses being different, to ensure distinguishability. To determine whether this is possible for the three-qubit case, we analyze the difference between real transfer matrices corresponding to h = 0 and h = 2: where a(s), b(s), and d(s) are functions of {κ k }, {∆ k } too lengthy to include here. Seeing that this matrix is full-rank, we determine that there is no complex driving function ǫ(t) which results in exactly equal outputs at all time.
To obtain an approximate three-qubit parity measurement for large measurement times, it is sufficient to set the steady-state responses equal for h = 0 and h = 2 (equality of the h = 1 and h = 3 responses follows by symmetry). In order to accomplish this, we note that: A can be expressed as the sum of a diagonal matrix and a rank-one matrix. This allows us to invert it analytically, using the Sherman-Morrison formula [31]: The output corresponding to the steady state of the resonator system is: To match steady states in the |Q| = 3 parity measurement from [18], we impose the following condition on the detunings of the two resonator modes (∆ 0 , ∆ 1 ): In order to constrain κ, we maximize the difference in Re(α out ) between the h = 0 and the h = 1 cases. This is achieved when κ 0 = κ 1 = 2χ.

Register Equations of Motion
The most economical equation for simulating the desired homodyne measurement is a reduced master equation acting only on the qubit register. Obtaining and solving such an equation has been the strategy of previous works [13][14][15][16]18]. Here, we provide a simplified derivation of the reduced master equation. We begin by expressing the register state using a partial trace acting on the state ansatz given in Eq. (7): In order to determine the dynamics of this reduced state, we use the integral representation of the partial trace over the harmonic oscillator [32]: Being a dummy variable, the vector of amplitudes α in the integral above does not depend on time or on the Wiener increment, resulting in a simple expression for the reduced master equation: We evaluate the deterministic and stochastic terms separately in the following subsections.

Deterministic
Before calculating the partial trace of the deterministic Lindbladian L, we note that, for an operator R supported only on the resonator Hilbert spaces, the partial trace acting on a commutator or dissipator annihilates it: Also, for any operator Q supported only on the register, the partial trace acting on a commutator or dissipator is expressed in terms of inner products of coherent states, similar to the partial trace of ρ in Eq. (23): The only partial trace which does not benefit from these simplifications is To write it succinctly, we introduce modified density matrices where O k ′ is an operator on the kth resonator space, and o k ′ is a scalar: noting that O k ′ and o k ′ may also depend on the register basis state |j , or on other variables.
The commutator and its partial trace can now be easily expressed: Permuting terms, we write the integral in Eq. (33) as a matrix element of the state of the k ′ th resonator mode: We use the resolution of polynomials in the raising and lowering operators [32, Section 2.9.3]: m,n c m,n a m a †n = 1 π α∈C c m,n α mᾱn |α α|d 2 α The integral in Eq. (34) reduces to: We substitute back in: Here, (A • B) i,j A i,j B i,j defines the elementwise (or Hadamard) matrix product. We write the unconditional master equation by transforming the qubit-only Lindbladian, and adding on a term which accounts for the qubit-resonator coupling (see Eq. 3, 26 -29, and 37): where H Q = 1 2 l∈Q Ω l + k∈C χ k,l σ z,l . In order to simulate evolution under this master equation, it is convenient to eliminate fast-rotating terms by expressing the master equation in a frame rotating with H Q . This has the effect of eliminating H Q from the Lindbladian:

Markovianity
It is interesting to note that the coupling Lindbladian in Eq. (37), though it generates a completely-positive trace-preserving map, is non-Markovian. This is not surprising, since the Markov approximation is the result of a weak-coupling assumption, and fast quantum measurement requires strong coupling. This has been confirmed in the case of a single-qubit measurement. The coupling Lindbladian, though it can be written in explicit Lindblad form, has a decay rate associated with the dephasing operator which is not necessarily positive [13]. In this section, we prove non-Markovianity of the coupling Lindbladian in the general case, and we examine the consequences of this property of the Lindbladian by numerical simulation.
To show that the coupling Lindbladian is non-Markovian, we note that the action of a Markovian Lindbladian on a density matrix in an N -dimensional Hilbert space can be expressed as [33,34]: where the coefficients ai,j Fi Fj form a Hermitian, positive-semidefinite matrix, and the operators F i form an orthogonal basis under the Hilbert-Schmidt inner product, with F 0 =1/ √ N . In the remainder of this section, we will derive a minimal set of operators F i for our coupling Lindbladian, and show that the resulting coefficient matrix is non-positive, precluding a true Lindblad representation.
In order to express the Hadamard product from Eq. (37) in the form given by Eq. (41), we use a dyadic product formula (see [35] for further applications of Hadamard and dyadic products to Lindbladians). The term of interest is given by: The matrixP k given in Eq. (38) is a dyadic product, so we can rewrite the elementwise product as a conjugation by diagonal matrices: whereα k i,j = α k,i δ i,j , and we have used Eq. (5) from [35]. We can now express the action of the Lindbladian in terms of a non-orthonormal set of operators: whereσ z,k l∈Q χ k,l σ z,l . We now show the non-positivity of the coupling Lindbladian in the case of a single mode (dropping the index k). To put Eq. (44) in the form of Eq. (41), we decompose the operatorsα andασ z in terms of1, F 1 and F 2 : Substituting into Eq. (44), we can derive the elements of the coefficient matrix: Therefore, Diagonalizing, we see that the eigenvalues are x ± √ x 2 + 1; the matrix has one negative eigenvalue, since x is real. This, in turn, implies that the coupling Lindbladian is always non-Markovian for |C| = 1, as long as |Q| > 1.
Note, however, that if only one qubit is present, the operators F 0 , F 1 and F 2 derived above are always linearly dependent. Therefore, the above argument is inapplicable in the one-qubit case. However, we observe that the coefficient of the dissipator term of the single-qubit pseudo-Lindblad equation can be negative in some time intervals during transient evolution, showing that the one-qubit evolution also has non-Markovian features [13]. We expect that non-Markovianity will be the general case for multi-mode coupling Lindbladians, since the sum of multiple coefficient matrices with negative eigenvalues is not necessarily positive.
To see the effect of this non-Markovianity on the performance of the measurement, we introduce two parity eigenstates: and simulate their evolution under the deterministic master equation (Eq. (40)), subject to a piecewise-quadratic input pulse, detailed in Fig. (1). Markovian dynamics produce a trace distance 1 /2 tr between the time-dependent states corresponding to initial states |ψ ± which is monotonically decreasing [36]. Plotting this quantity in Fig. (2), we see that, as the pulse is turned off, there is a clear increase, indicating non-Markovian behaviour.
Since the trace distance is a measure of state distinguishability, this increase implies  Increase of trace distance between ideal post-measurement states as a result of pulse turn-off, in the three-qubit parity measurement described in [17,18]. Pulse is shown in Fig. (1). System parameters are, in units of χ: ∆ 0 = 3, ∆ 1 = −3, κ 0 = κ 1 = 2. No intrinsic sources of decoherence have been included.
that the ability of this measurement to produce high-accuracy post-measurement states is higher than an otherwise-identical Markovian measurement would allow.

Stochastic
The stochastic master equation governing homodyne measurement differs from the deterministic master equation only in the term M a out e −iφ (ρ) dW , given in Eq. (6). Upon tracing out the resonator, the stochastic term in the master equation becomes: Since the third term contains a full trace, we calculate it first: This is also the expected value of the photocurrent in Eq. (5). We calculate the other two terms simultaneously, using the notation of Eq. (31) for concision: We see that tr C M a out e −iφ (ρ) = M[c Q ](ρ Q ). Since c Q is diagonal, it commutes with the Hamiltonian, and is identical in the rotating frame, yielding the following stochastic master equation: In the following section, we numerically integrate the stochastic master equation, using the resulting state fidelity as a measure of performance.

Simulation
In order to assess the accuracy with which joint measurements can be made directly, we focus on the three-qubit parity measurement from [17,18]. We set the resonator parameters according to Eq. (22) (with both κ 0 and κ 1 set to 2χ).

Methods
We simulate the evolution of ρ Q over the interval [t, t + dt] in two steps: we first determine the time-dependent amplitudes {α k,j (t)} through numerical integration, using a 4 th /5 th -order adaptive Runge-Kutta stepper, using the pulses detailed in Fig. (1), producing the state-dependent response shown in Fig. (3). We then use these values of {α k,j (t)} to formulate the time-dependent reduced master equation (Eq. (55)), and use an order-1.5 stochastic Runge-Kutta method [25] to integrate it. We repeat this for 10 5 uniformly-spaced timesteps on the interval [0, τ ], where τ is the total measurement time (taken to be 13.5 /χ throughout).
In order to verify the correctness of the above simulation, we calculate the minimum/maximum eigenvalues, traces, deviations from hermiticity ρ − ρ † ∞ and purities tr ρ 2 for a typical trajectory. In order for these quantities to be meaningful when plotted for a single trajectory, the algorithm used has to be strong (in the terminology of stochastic differential equations [25]). Deviations from hermiticity on the order 10 −15 are typical, as are deviations from unit trace on the order 10 −13 , likely caused by numerical rounding error. For the simulations discussed in this paper, the remaining checks are satisfied to within ∼ 10 −49 .

Figure of Merit
To calculate the performance of the continuous measurement, we evaluate its ability to produce one of two definite-parity entangled states, projected from an initial state |+ + + , and to produce a measurement record which correctly identifies the parity of the final state. To this end, we calculate the quantum state fidelity: where |ψ ± is one of the parity eigenstates from Eq. (51). This fidelity is calculated for a post-selected ensemble determined by the signal s(τ ), which is a weighted integral of the measurement record: If this signal is positive (negative), we infer that the parity of the post-measurement state is even (odd). The choice of filter function f (t) has a significant impact on the performance of the measurement. Though it has no effect on the post-measurement states themselves, it can increase or decrease the probability of an incorrect assignment of states to ensembles, affecting the fidelity indirectly. We compare (in Fig. (4)) the distributions of s(τ ) corresponding to the uniform filter (f (t) = 1) and a matched filter (a filter function proportional to the measurement record, see [37,38]), shown in Fig. (5).
This matched filter is derived from a simplified model of the dynamics, in which we assume that the state immediately collapses to a uniform mixture of the computational basis states from one of the parity eigenspaces. The filter is then the expected value of the nominal photocurrent in the even-parity subspace, normalized  The approximate matched filter given by Eq. (58), for use in obtaining a measurement signal s(τ ) given a photocurrent time trace j(t). The filter function can be thought of as assigning an importance to the photocurrent emitted at time t; this filter function assigns high importance during the steady-state phase of the measurement, and low importance at other times.
to have a mean value of 1: where j + (t) = tr c Q + c † Q Π + (see Eq. (53)). Here, we have defined Π ± to be projectors onto the even/odd parity subspaces of B, and we have taken advantage of the fact that tr c Q + c † Q Π − = −tr c Q + c † Q Π + . In the following section, we use the post-selected state fidelity to examine the performance of a nominal measurement.

Results & Discussion
We simulate the action of the measurement in the presence of decoherence from the measurement itself and qubit dephasing noise. Selecting χ as the natural scale for frequency, we take γ z to be χ /300, as in [18]. If χ ∼ 1 × 2π MHz, this would correspond to a dephasing time of ∼ 50µs.
A histogram of integrated photocurrents is given in Fig. (4). The separation between the two Gaussian peaks therein is visibly greater when using the approximate matched filter from Eq. (58), indicating that the probability of incorrect assignment is decreased. The resulting state fidelities are given in Fig.  (6). To understand the features of the fidelity distribution, we plot (in Fig. (7 the expected fidelity for a simplified model of measurement in which the state immediately collapses to one of the parity eigenstates, and is then subject to both intrinsic and measurement-induced decoherence. In order to compare this average state fidelity with known performance thresholds for error-correcting architectures [39], we calculate the post-measurement state fidelity obtained by performing a circuit-based measurement [12]. This measurement uses the circuit in Fig. (8).
We assign to each operation in the circuit (state preparation, memory, CNOT and single-qubit measurement) a failure probability p. Failure corresponds to: • unintentional preparation of the |1 state on the ancilla qubit, (|100 + |010 + |001 ) without the coupling Lindbladian (solid black line) coincides exactly with the case in which the coupling Lindbladian is present (black circles). The initial state 1 2 (|100 + |100 + |100 + |111 ) produces a markedly different effect, its fidelity decaying much more quickly when the coupling lindbladian is present (solid blue line) than when it is absent (dashed black line).

•
• • |0 Figure 8 A measurement circuit for the operator ZZZ, similar to the four-qubit circuit used in [12].
• insertion of a random one-qubit Pauli operator (X, Y , or Z) after a memory operation with probability 1 /3, Plotting this fidelity in Fig. (9), we see that a state fidelity of 94% corresponds to an error rate of 1.4−1.5%, and that an output state fidelity of 96−98% would match the performance of the gate-based model close to the error-correction threshold. This indicates that, given the parameters for existing state-of-the-art transmon/cavity systems, multi-qubit measurements of quality near that required for fault-tolerant error correction could be performed, but that further improvements in control design and/or hardware parameters would be needed for threshold error rates to be attained. State fidelity F ± = ψ ± |ρ|ψ ± for the gate-based measurement acting on the initial state |+ + + . Inset: small p regime, in which F ± ∼ 1 − 4p.

Conclusions and Future Work
The formalism presented in this article facilitates the design of a class of quantum measurement devices, with applications in fault-tolerant quantum computing architectures and remote entanglement preparation. Using a pulse with few free parameters, and an approximate matched filter, it is possible to limit transient-induced decoherence, achieving a high state fidelity, comparable to error models studied in the fault tolerance literature operating slightly above the fault-tolerance threshold. Further advances in the design of quantum hardware, such as decreases in the dephasing rate, will permit higher-fidelity implementations of this protocol. In addition, there are several purely theoretical avenues to be explored, which can inform the feasibility of this idea.
In future work, we will attempt to eliminate or minimize the decoherent portion of the coupling Lindbladian, by selecting a control pulse optimally. Such a control pulse can have a quadrature-phase component, unlike the pulse used in this manuscript. Decreasing this decoherence will permit stronger driving, which in turn enables shorter measurement times, reducing the effective strength of intrinsic decoherence. We will also attempt to minimize or correct unwanted rotations, which occur as a result of the ac-Stark shift and the non-zero imaginary part of c Q . It is known that non-linear filtering of the output photocurrent and multi-qubit gates can be used for this purpose [15], but the performance of a scheme involving efficient filters and single-qubit gates has yet to be examined.

Remarks
Numerical simulations were performed using the libraries homodyne sim and sde solve, available on Github. determines the value of the "dielectric constant" of the linear resonator system, whose frequency can be probed by the phase shift of scattered coherent radiation. From this point of view, this mode-mode coupling does not alter the fundamental strategy, or the basic formalism, of the parity-measurement setup.
However, this term does produce some new complications. The response is now not only determined by the χ k,l parameters, but depends separately on the g and ∆ parameters. The response is not automatically determined by only the Hamming weight of the state of the qubits; but the difference in the response of states with the same Hamming weight can be made weak if the overall detuning is large compared with the differences of qubit frequencies. If the κ parameters of the two modes are made significantly different, then, according to Eqs. (22), the bare resonator frequencies can be made significantly different; this in turn will make the influence of the mode-mode coupling term very small, in accordance with another RWA.
We will report in a forthcoming publication [41] on the re-optimization of the parity measurement with all the effects of this Hamiltonian Eq. (61) taken into account. It is especially worthwhile to do this, in light of the fact that a similar analysis of the dissipative terms shows that we can design a structure which both does a good parity measurement and has built-in protection against Purcell decay. The relevant decay term can be written [41]: The inner sum over the cavity modes gives the opportunity for cancellation of the coefficient of each σ − l term, since the g factors can be set to have either sign. Our preliminary work indicates that this indeed can be done, in a way that is consistent with satisfying the conditions for the parity measurement. We will provide the details of this analysis in an upcoming publication.

Appendix B Pointer State Solution
Here we sketch the proof that the pointer states in Eq. (7) constitute a family of solutions of the multi-qubit multi-mode master equation (Eq. 3), when no qubit damping is present γ −,l = λ k,l = 0. First, one establishes by direct substitution that the "diagonal states" are a solution to Eq. (3), where the functions α k,i (t) are the solutions to Eq. (10). This is established by observing that the entire right-hand side of the equation commutes with σ z,l . The form of the remaining terms (involving the off-diagonal qubit operator |i j|) does not follow from any simple linearity argument. However, we can seek an extended solution having these new terms, HereÔ i,j (t) is an undetermined operator on the Hilbert space of the resonators, which is a function of the register states |i and |j , as well as being an undertermined function of time.Ô is fixed by positivity; it can be shown (details in [41]) that ρ ≥ 0 forcesÔ i,j (t) to be of the form This confirms Eq. (7), with the added information that the diagonal coefficients are independent of time, which is a property of the equation of motion that we subsequently derive for ρ Q (Eq. (39)).