Combining Floquet and Lyapunov techniques for time-dependent problems in optomechanics and electromechanics

Cavity optomechanics and electromechanics form an established field of research investigating the interactions between electromagnetic fields and the motion of quantum mechanical resonators. In many applications, linearised form of the interaction is used, which allows for the system dynamics to be fully described using a Lyapunov equation for the covariance matrix of the Wigner function. This approach, however, is problematic in situations where the Hamiltonian becomes time dependent as is the case for systems driven at multiple frequencies simultaneously. This scenario is highly relevant as it leads to dissipative preparation of mechanical states or backaction-evading measurements of mechanical motion. The time-dependent dynamics can be solved with Floquet techniques whose application is, nevertheless, not straightforward. Here, we describe a general method for combining the Lyapunov approach with Floquet techniques that enables us to transform the initial time-dependent problem into a time-independent one, albeit in a larger Hilbert space. We show how the lengthy process of applying the Floquet formalism to the original equations of motion and deriving a Lyapunov equation from their time-independent form can be simplified with the use of properly defined Fourier components of the drift matrix of the original time-dependent system. We then use our formalism to comprehensively analyse dissipative generation of mechanical squeezing beyond the rotating wave approximation. Our method is applicable to various problems with multitone driving schemes in cavity optomechanics, electromechanics, and related disciplines.


I. INTRODUCTION
In the past years, cavity optomechanics and electromechanics [1] became a mature field. Starting with first experiments demonstrating ground-state cooling of mechanical motion [2,3], electromagnetic fields are now commonly used to manipulate and measure quantum states of mechanical resonators. Examples of nonclassical mechanical states include squeezed states, prepared using reservoir engineering [4] or parametric effects [5], as well as Fock states and their superpositions [6][7][8]. Preparation of these states presents a crucial first step towards creating more complex mechanical states such as truly macroscopic mechanical superpositions [9][10][11] or towards generation of entanglement in optomechanical and electromechanical systems [12][13][14].
In many applications, the system of interest is described by a Hamiltonian that is quadratic in the canonical operators of radiation and mechanical motion. As the unitary dynamics are accompanied by single-excitation loss and gain and-in case of backaction-evading measurements-homodyne detection (nonunitary processes which are linear in the canonical operators), the phase-space description of the quantum state retains a Gaussian character [15]. We can therefore characterise the state by the first and second statistical moments of the quasiprobability density-the mean vector and covariance matrix of the canonical operators-instead of the density operator. This presents a great advantage since we can formulate equations of motion for the mean and covariance matrix [16,17] which can be solved more efficiently than the master equation that describes the evolution of the density operator. Particularly the covariance matrix is important for investigating Gaussian quantum features of the resulting quantum state such as quadrature squeezing or entanglement [15]; in the steady state, its form can be obtained by solving an algebraic equation of the Lyapunov (for unconditional dynamics) or Riccati (for conditional evolution) type [16,17].
A crucial obstacle for a more widespread application of these techniques is the explicit time dependence of the driving electromagnetic fields. Dissipative preparation of mechanical states [4,[9][10][11][12][13][18][19][20][21][22][23][24][25] and tomographic backactionevading measurements of mechanical motion [26][27][28][29][30][31][32][33][34] rely on driving the system with multiple fields at different frequencies while parametric squeezing requires modulation of the optical spring [5,[35][36][37]; both of these approaches result in time-dependent optomechanical Hamiltonians. The steady-state Lyapunov equation can then only be applied under the rotating wave approximation (RWA) which neglects fast oscillating terms in the interaction and only keeps those that are resonant. This approximation can work well as long as both the cavity damping rate and optomechanical coupling strength are much smaller than the mechanical frequency. When one of these quantities becomes comparable to the mechanical frequency, the RWA breaks down and the fast oscillating terms can no longer be neglected. In this case, Floquet techniques can be used to find the steady state by expanding the system operators into their Fourier components [38][39][40][41][42]. Their application is, however, often cumbersome and lacks a simple, general approach.
In this article, we show how the Lyapunov and Floquet techniques can be combined to obtain a straightforward approach to finding the steady state of optomechanical and electromechanical systems with periodic Hamiltonians. We derive a general method which can be directly applied to the Lyapunov equation of a time-dependent quantum system to obtain a time-independent Lyapunov equation in the larger Floquet space. Importantly, the frequency components of the relevant quantities can be directly read off from the time-dependent equations of motion for the canonical operators; this feature makes our strategy particularly easy to implement for a wide variety of critical problems in optomechanics and electromechanics.
After presenting the general method, we apply it to the problem of dissipative mechanical squeezing. We consider both the usual squeezing obtained by cooling a mechanical Bogoliubov mode [4] as well as the recently proposed combination of dissipative and parametric squeezing possible for levitated particles [36] and analyse the resulting squeezing in both cases beyond the usual RWA. This approximation is invalidated not only by the nonzero sideband ratio (the ratio between the cavity decay rate and the mechanical frequency) but also by the ultrastrong coupling regime which is being reached in an evergrowing number of optomechanical and electromechanical systems [43][44][45][46][47][48][49][50][51][52]. A detailed systematic analysis of the effects that the counterrotating terms will thus have important consequences not only for immediate applications in state-of-the-art optomechanical and electromechanical systems but also for advanced state preparation techniques, possibly going beyond the Gaussian regime.

A. Lyapunov equation
The dynamics of an open quantum system can be described by the master equation where 2ρ † is the Lindblad superoperator with collapse operator [53]. The full dynamics can, in principle, be solved from this equation but this is not always feasible in practice. Optomechanical systems are one such example; the large thermal occupation of the mechanical bath requires prohibitively large Hilbert-space dimensions in the Fock-state basis to be considered for numerical simulations.
Full description of the dynamics without solving the master equation (1) is possible for Gaussian systems, which are characterised by Gaussian quasiprobability distributions in phase space and as such can be fully described by the first and second statistical moments of these distributions. A system is Gaussian when its HamiltonianĤ is quadratic and the collapse operators n linear in canonical operators; then, we can solve for the covariance matrix describing the state of the system instead of the density matrix. (Note that the first moments remain zero when the Hamiltonian contains no linear terms and the collapse operators no constant terms.) The covariance matrix Γ obeys the Lyapunov equatioṅ where Γ ij = ⟨r irj +r jri ⟩ − 2⟨r i ⟩⟨r j ⟩, (3a) with the canonical commutation relation where σ is the one-mode symplectic matrix. The drift (A) and diffusion (N) matrices can be found from the master equation [17] or from the Langevin equationsṙ whereξ = (q 1,in ,p 1,in , . . . ,q N,in ,p N,in ) ⊺ is the vector of noise operators. The diffusion matrix can be found from the correlation properties of the noise via The Lyapunov equation (2) is a powerful tool in optomechanics where it can be used to evaluate the dynamics of the system or to find its steady state from the algebraic form of the Lyapunov equation The process of solving Eq. (2) or (7) is straightforward as long as the system operators are not explicitly time dependent. When the Hamiltonian or coupling to the environment becomes time dependent, the drift matrix (and possibly also the diffusion matrix) becomes time dependent as well. While the differential Lyapunov equation (2) can still be solved in this case, a straightforward solution of the algebraic Lyapunov equation (7) is not possible.

B. The Floquet-Lyapunov method
A common approach to problems with periodic time-dependence is based on the Floquet formalism, which allows us to transform a set of periodic linear differential equations into a larger set of linear differential equations without explicit time dependence [54]. Here, we will show how to use this method to obtain a time-independent Lyapunov equation for explicitly time-dependent dynamics. For simplicity, we will assume that only the Hamiltonian is τperiodic,Ĥ(t + τ ) =Ĥ(t), leading to a periodic drift matrix A(t + τ ) = A(t), while the coupling to the environment (and the diffusion matrix N) remains explicitly time independent. We will only present a general recipe for obtaining a time-independent Lyapunov equation in this section; the mathematical justification of these steps is presented in Appendix A.
We start by expressing the operators and the drift matrix in terms of their Fourier components at frequency ω = 2π τ ,r The frequency componentsr (0) ,r The individual frequency components of the noise operators have the same correlation properties as the initial noise, (m)⊺ l ⟩ = 0 for n ≠ m or k ≠ l. Note also that in Eqs. (9), only the frequency components with ±m ± n > 0 should be included when evaluating the sums; effectively, we haver where we introduced the drift matrix with I N denoting the N -mode (i.e., 2N -dimensional) identity matrix andξ F = (ξ (0) ,ξ s , . . .) ⊺ being the vector of noise operators in terms of their frequency components. In these expressions, we use the subscript F to signify that the quantities are defined in the Floquet space. We can now formulate the Lyapunov equation in the Floquet spacė with the time-independent drift matrix given by Eq. (11) and diffusion matrix N F = diag(N, N, ...). The covariance matrix can be expressed in the block form where the blocks are defined as the covariances between the corresponding frequency components of the quadrature operators; we thus have The steady state can now be calculated from Eq. (12) by settingΓ F = 0, which gives a linear, time-independent algebraic equation.
The system now becomes effectively infinite dimensional and must be suitably truncated for numerical simulations. This truncation has to ensure that the solution (contained in the zeroth Brillouin zone, i.e., the covariance block Γ (0) ) converges. As we shall show in the following examples, only a small number of Brillouin zones is typically sufficient, making our approach more feasible than direct solution of the master equation or finding the long-time limit of the differential Lyapunov equation with time-dependent drift.

C. Illustration: optomechanical cooling
To illustrate our formalism, we now consider optomechanical sideband cooling and show how the Floquet-Lyapunov techniques can be used to obtain the correct steady-state mechanical occupation in the rotating frame. For an optical cavity (with annihilation operatorĉ) coupled to a mechanical oscillator (annihilation operatorb), the linearised optomechanical Hamiltonian (in the lab frame of the mechanical motion) takes the form [1] where ∆ = ω c − ω L is the detuning between the cavity frequency ω c and driving frequency ω L , ω m is the mechanical frequency, and g is the coupling constant. This is a time-independent Hamiltonian and the dynamics can be solved directly by using the Lyapunov equation.
Introducing the canonical operatorsq 1 = (ĉ +ĉ † ) 2 for the cavity field andq 2 ,p 2 for the mechanics (defined similarly), we can describe the dynamics using the Langevin equationṡ Here, κ and γ are the damping rates for the cavity field and the mechanical mode and the input noise operators of the optical and mechanical baths have the correlations withn denoting the average occupation of the thermal mechanical bath; all other correlations are zero. From these equations, we obtain the Lyapunov equatioṅ The optimal cooling performance is then achieved for driving on the lower mechanical sideband, ∆ = ω m , and with a sideband-resolved system, κ ≪ ω m [55,56]. This is due to the passive beam-splitter part of the optomechanical interaction being resonantly enhanced while the active two-mode-squeezing part (responsible for depositing energy into the mechanical mode) is far off resonant. The steady-state mechanical occupation, obtained by settingΓ = 0 in Eq. (18), can be found as n f = 1 4 (Γ 33 + Γ 44 − 2) and is plotted Fig. 1 as the solid blue line. To demonstrate the Floquet formalism, we start from the Hamiltonian (15) driven on the red sideband, ∆ = ω m , and move to the interaction picture with respect toĤ 0 = ω m (ĉ †ĉ +b †b ). Thus, the Hamiltonian becomeŝ By applying the RWA (valid for high mechanical frequencies, ω m ≫ κ, g), we are left with the beam-splitter Hamilto-nianĤ This results in the drift matrix while the diffusion matrix N remains the same as before. The RWA clearly demonstrates the mechanism of coolingmechanical excitations are swapped to the cavity field from which they leak out through the cavity mirrors-but does not properly estimate the final mechanical occupation since it does not include the effective temperature of the optical bath induced by the two-mode squeezing part of the optomechanical interaction. This is visible in regimes where the cavity damping or the optomechanical coupling start approaching the mechanical frequency, invalidating the assumptions necessary for the RWA (see the black dotted line in Fig. 1, calculated in the same way as for the full model).
The full rotating-frame Hamiltonian (19) is, however, periodic with frequency 2ω m so we can include the effect of the counterrotating terms in our simulations by applying the Floquet-Lyapunov approach. We begin by forming the equations of motion for the canonical operators, From these equations, we can readily read off the frequency components of the Floquet-space drift matrix, where A 0 = diag(−κ, −κ, −γ, −γ) is the part of the drift matrix associated with damping, Keeping the first three Brillouin zones (i.e., defining the Floquet vectorr F = (r (0) ,r (N, N, N).
The steady-state mechanical occupation (dashed yellow line in Fig. 1) now shows perfect agreement with the full model.

III. DISSIPATIVE SQUEEZING BEYOND THE ROTATING WAVE APPROXIMATION
In the previous example of optomechanical cooling, the use of the Floquet method is unnecessary as the system can be represented in a frame where the Hamiltonian is time independent. In the rest of this paper, we focus on more complicated driving schemes for which the optomechanical system remains time-dependent regardless of the reference frame. This is generally the case when the system is driven with multiple fields at different frequencies or, equivalently, using a single drive with a modulated amplitude. We focus on the simplest nontrivial example-one mechanical mode and one cavity field subject to a two-tone drive. We perform a detailed analysis of the mechanical squeezing that can be achieved in this system originally proposed by Kronwald et al. [4] in section III A. In section III B, we then analyse a recently proposed scheme for generating steady-state mechanical squeezing in levitated systems by a combination of parametric and dissipative squeezing [36].
A. Mechanical squeezing with a two-tone drive Following Ref. [4], we start from the full optomechanical Hamiltonian under two-tone driving, where g 0 is the single-photon coupling strength and η ± are the amplitudes of the drives at frequencies ω ± = ω c ± ω m . Under this driving, the cavity field acquires a large classical amplitude α + e −iω+t + α − e −iω−t ; introducing the linearised coupling rates g ± = g 0 α ± (assuming, for simplicity, α ± ∈ R) and moving to the interaction picture with respect tô H 0 = ω −ĉ †ĉ + ω mb †b , we obtain the interaction Hamiltonian To follow the Floquet approach, we first formulate the equations of motioṅ Moving to the Floquet basis, we obtain the drift matrix where where we introduced the matrix where Γ (0) is a 4 × 4 matrix. The steady-state squeezed and antisqueezed mechanical variances are The variances simulated by the Floquet-Lyapunov method and under the RWA are shown in Fig. 2. For small coupling (g ± ≪ ω m ) and a sideband resolved system (κ ≪ ω m ), there is a good agreement between the two approaches. As the coupling or the cavity decay is increased, particularly the squeezed variance is affected by the presence of the counterrotating terms. This result is not surprising-as the counterrotating terms give rise to a phase-independent backaction limit, they add a fixed amount of noise to both mechanical quadratures. As the squeezed variance is smaller, the same amount of added noise corresponds to a larger relative increase of the variance than for the antisqueezed quadrature. Similar to the case of mechanical cooling in Fig. 1, the Floquet solution becomes unstable for strong coupling, behaviour not visible with the RWA.

B. Parametric and dissipative squeezing for a levitated particle
In the previous example, only one frequency component (oscillating at ω = 2ω m ) was present in the system. Now, we turn our attention to a system where also the next component is relevant-a levitated particle squeezed by a combination of parametric and dissipative squeezing [36]. The potential for the particle's centre-of-mass motion is defined by the laser beam used for levitation; its scattering into an empty cavity mode provides the optomechanical interaction [46,51,57]. To achieve strong squeezing, the optical tweezer amplitude is modulated at twice the mechanical frequency, E tw (t) = E 0 [1 + α cos(2ω m t)], where α ∈ (0, 1) is the modulation depth, modulating both the mechanical potential and the optomechanical coupling rate. Since the potential is proportional to the tweezer intensity (i.e., the square of the amplitude) and the optomechanical coupling to the tweezer amplitude, the system dynamics are characterised by the Hamiltonian [36] The parametric squeezing is provided by the modulation of the trapping frequency according to the second term in Eq. (35). To see how the dissipative squeezing arises, we set the detuning to the red mechanical sideband, ∆ = ω m , and move to the rotating frame with respect to the free oscillations,Ĥ 0 = ω mĉ †ĉ + 1 2 ω m (q 2 2 +p 2 2 ); under the RWA, we then obtain the interaction-picture Hamiltonian Introducing the mechanical Bogoliubov modeβ = (2b + αb † ) √ 4 − α 2 (which satisfies [β,β † ] = 1), we can rewrite this Hamiltonian asĤ where λ eff = λ (4 − α 2 ) 8. From this expression, we see that the Bogoliubov mode undergoes slow oscillations at frequency ω m α 2 ≪ ω m (the first term on the right-hand side) while being squeezed parametrically (the second term in the Hamiltonian); dissipative squeezing is provided by the optomechanical interaction in the last term. Remarkably, these two squeezing techniques can combine and provide stronger squeezing than either strategy alone [36].
To describe the generated squeezing beyond the RWA, we write the interaction Hamiltonian (35) including the counterrotating terms, .
From this Hamiltonian, we derive the equations of motioṅ from which we obtain the drift matrix where Again, the diffusion matrix is the same as before and the RWA corresponds to working with the zeroth Brillouin zone only.
We compare the results of the Floquet-Lyapunov approach and the RWA in Fig. 3. Again, the squeezed quadrature is more sensitive to the validity of the RWA than the antisqueezed quadrature; generally, this effect seems more pronounced now owing to the larger number of counterrotating terms present in the system. Surprisingly, there exist regimes for the levitated particle where the antisqueezed variance is smaller with the counterrotating terms than under the RWA (see the red lines in Fig. 3(f)). This result points to a nontrivial role that the counterrotating FIG. 3. Parametric and dissipative mechanical squeezing for a levitated particle as a function of (a)-(c) the coupling strength g and (d)-(f) modulation depth α. The lines have the same meaning as in Fig. 2 and we again plot the ratio of the variances (top), the squeezed variance (middle) and the anti-squeezed variance (bottom). The parameters used for the simulations are γ ωm = 10 −9 andn = 2 × 10 7 .
terms play in this situation, which might be harnessed by optimising the modulation phase in the tweezer amplitude, E tw (t) = E 0 [1 + α cos(2ω m t + φ)] (we considered φ = 0 for simplicity here).
The dynamical stability of the system is now more complicated than in the case of dissipative squeezing via a two-tone drive. Similar to the previous case, the system becomes unstable for large coupling, which is not captured by the RWA calculation. Additionally, the system is unstable (both under RWA and with the Floquet-Lyapunov approach) also for large modulation depth, where the parametric squeezing effect becomes too strong and for weak optomechanical coupling. We can understand the latter result by noting that, for g → 0, only the parametric squeezing effect persists; for the modulation depth considered in Fig. 3, this would correspond to more than 3 dB of squeezing, making the system unstable.

IV. CONCLUSIONS
In this manuscript, we presented an efficient method for analysing dynamics of Gaussian systems with periodic Hamiltonians. Our approach is based on deriving equations of motion for the canonical operators in the Floquet space from which a Lyapunov equation for the covariance matrix of the system's Wigner function can be formulated. In this way, we arrive at a time-independent algebraic equation for the steady state of the system starting from a time-dependent one; the price we pay for this-increased size of the Hilbert space-is not prohibitive for numerical simulations. We expect the technique to find wide use in cavity optomechanics and electromechanics where timedependent coupling rates are commonly used for dissipative state preparation. In these systems, large thermal phonon occupations of the mechanical reservoirs preclude direct simulations of the master equation in the Fock basis.
We demonstrated the use of our techniques on dissipative generation of mechanical squeezing in conventional dispersive and novel levitated optomechanical systems and performed a detailed analysis of the role counterrotating terms play in these systems. We showed the intuitive result that the squeezed quadrature is more sensitive to the validity of the RWA than the antisqueezed quadrature owing to its overall lower amount of noise. Importantly, even though counterrotating terms have stronger effect in the recent proposal of combined parametric and dissipative squeezing for levitated particles [36], sub-vacuum squeezing is still possible with state-of-the-art systems; our simulations show that about 2 dB of squeezing (with about 11 dB of antisqueezing) can be generated with parameters similar to Ref. [52]. Finally, we note that the presence of counterrotating terms can also affect the dynamical stability of the system (particularly when approaching the ultrastrong coupling regime) which cannot be seen from simple analysis under the RWA.
Our approach can be readily adapted for various other state preparation schemes that rely on multitone driving, such as generation of Gaussian entanglement between two mechanical resonators [12,13]. Our technique allows a straightforward analysis of dissipative state preparation beyond the RWA and a direct evaluation of approaches to mitigate the undesired effects of the counterrotating terms. One such possible technique is using cavity fields with squeezed input noise. This strategy has already been used to suppress the heating associated with counterrotating terms in sideband cooling [58,59]; it would be interesting to see whether its benefits persist also in more complex dissipative dynamics such as generation of one-and two-mode mechanical squeezing.
Our strategy could also be used to analyse conditional dynamics of optomechanical systems under backactionevading measurements [31,34]. In these schemes, driving on both mechanical sidebands is used to perform a quantum nondemolition measurement of one of the mechanical quadratures; similar to dissipative state preparation, the effects of counterrotating terms are often neglected under the RWA. Since the measurement at the output is homodyne (projecting the system onto a Gaussian state), the system can still be fully characterised using the covariance matrix formalism; the important distinction is that the dynamics of the conditional state obey an algebraic equation of the Riccati type [16]. The Riccati equation is closely related to the Lyapunov equation-it operates with the same drift and diffusion matrices (which account for the unconditional part of the evolution) with additional, nonlinear terms stemming from the weak measurement. Since the measurement is typically time-independent in backaction-evading measurements, extension of our formalism to Riccati-type equations should be straightforward.
All in all, dissipative preparation of nonclassical mechanical states and their tomography using quantum nondemolition measurements are important techniques for modern optomechanical and electromechanical systems. As we have discussed, they are often based on the use of periodic Hamiltonians in which resonant parts give rise to the desired dynamics and fast oscillating terms are neglected under the RWA. The RWA can be invalidated by two factors: nonzero sideband ratio (between the cavity damping and mechanical frequency), which stems from practical limitations in realistic experimental systems, and ultrastrong coupling, associated with improved experimental design. With our techniques, the dynamics of these systems can be efficiently analysed including the effects of the fast oscillating terms, allowing the new physics associated with these processes to be uncovered.
Let us consider a time-dependent Hamiltonian that is τ -periodic,Ĥ(t) =Ĥ(t + τ ), and quadratic in the canonical operators. Representing all modes with creation (â † i ) and annihilation (â i ) operators, we can write the time development of their expectation values where we defined the vector y = (α 1 , α * 1 , . . . , α N , α * N ) ⊺ . Since the linear system (A2) is τ -periodic with the corresponding frequency ω = 2π τ , we can decompose y and B(t) in their frequency components according to Defining the vector of frequency components y (n) = (α where I N is the N -mode (and hence 2N -dimensional) identity matrix. With the frequency decomposition of B(t) and the matrix Q(t), we can evaluate Eq. (A1); a straightforward calculation reveals that individual frequency components obey the differential equationẏ Transformation from the frequency components of the creation and annihilation operators (defined in terms of the complex exponential e inωt ) to the sine and cosine components of the quadratures (as introduced in Eqs. (8)) is possible by introducing the definitions Note that these quantities are real since, owing to the definitions of the frequency components in Eq. (A3), we have . Moreover, the normalisation of the operators is chosen such that the basis change is unitary; this choice also fixes the prefactor in front of the sum in the Fourier transform Indeed, introducing the vectors r (0) = (q s , . . .) ⊺ , we can write where we introduced the unitary matrix We can now transform the Langevin equation in the Floquet space,ẏ F = B F y F , intȯ where A F is given by Eq. (11) in the main text. So far we assumed that the system is homogeneous, which is valid only for the classical expectation values α i = ⟨â i ⟩, q i = ⟨q i ⟩, p i = ⟨p i ⟩. To include quantum fluctuations, we append the original linear system (A2) with a vector of noise operatorsξ y = (â 1,in ,â † 1,in , . . . ,â N,in ,â † N,in ) ⊺ ,ẏ = B(t)ŷ +ξ y , where we now work with the vector of operatorsŷ = (â 1 ,â † 1 , . . . ,â N ,â † N ) ⊺ instead of their expectation values y = ⟨ŷ⟩.
where k ∈ {c, s} andξ is the input quadrature noise of the original system, and N F = diag(. . . , N, N, . . .). The Floquet-space covariance matrix Γ F then obeys the Lyapunov equatioṅ limit agrees perfectly with the next level of the approximation (i.e., adding the two components for the frequency n = 2). For the simulations in Sec. III A, we therefore use five Brillouin zones (frequency components n ≤ 2) as shown explicitly in the drift matrix (29).
With the joint dissipative and parametric generation for a levitated particle, more Brillouin zones have to be considered as the interaction Hamiltonian includes higher frequency components; here, agreement is reached one step later in the approximation between the cases of n ≤ 2 and n ≤ 3 as shown in Fig. 4(b). For the simulations in Sec. III B, we use the first seven Brillouin zones (frequencies n ≤ 3), corresponding to the drift matrix (39).