Simulating quantum dynamical phenomena using classical oscillators: Landau-Zener-Stückelberg-Majorana interferometry, latching modulation, and motional averaging

A quantum system can be driven by either sinusoidal, rectangular, or noisy signals. In the literature, these regimes are referred to as Landau-Zener-Stückelberg-Majorana (LZSM) interferometry, latching modulation, and motional averaging, respectively. We demonstrate that these pronounced and interesting effects are also inherent in the dynamics of classical two-state systems. We discuss how such classical systems are realized using either mechanical, electrical, or optical resonators. In addition to the fundamental interest of such dynamical phenomena linking classical and quantum physics, we believe that these are attractive for the classical analogue simulation of quantum systems.

. Classical analogues of qubits. In (a) the two qubit eigenenergies E ± are shown to depend on the bias ε 0 and to display avoided-level crossing at ε 0 = 0 with a minimal distance Δ. The system is excited when the characteristic qubit frequency is about a multiple of the driving frequency, ΔE/ℏ = kω. Several possible classical realizations have been demonstrated: (b) two weakly coupled spring oscillators 12 , (c) a two-mode nano-beam 14 , (d) optomechanical system with two cantilevers 13 , (e) two coupled electrical resonators 5,46 , (f) two coupled polarization modes in an optical ring resonator 47 .
SCIENTIfIC REPORtS | (2018) 8:12218 | DOI: 10.1038/s41598-018-28993 -8 In the rest of this paper, we first present details of how the description of two classical oscillators can be reduced to the qubit equations of motion. We further consider situations when the mechanical-resonator system is driven by strong periodic or noise signals. We demonstrate the resulting interference fringes, which are remarkably similar to those in the quantum analogues.

Schrödinger-like classical equation of motion
To be specific, among diverse classical analogues of qubits, we consider mechanical layouts. Such oscillators can be realized as separate resonators [43][44][45] or as two modes of one resonator 10,11 . Among the advantages of such mechanical systems are that they operate at room temperature, have small size and high quality factors, and are only weakly coupled to the environment. Such system is described by the equations of motion where the displacement u i relates to the i-th oscillator (i = 1, 2, j ≠ i). For oscillators we assume equal effective masses m, equal damping rates γ, different spring constants ω = k m i i 2 , and weak coupling,  k k i c , see Fig. 1(b). It is important to stress that the other layouts in Fig. 1 are also described by the same equation (1). Therefore, what is written below is equally applicable to any classical two-state system. For example, for two coupled electrical resonators, Fig. 1(e), instead of displacement we have charges on respective capacitances, u i → q i , the role of masses is played by inductances, m → L, capacitances replace the spring constants, and k c → C c , and the damping is defined by the resistances, γ → R/L. Thus, the present scheme, inspired by ref. 5 , requires driving via capacitances. Alternatively, the scheme may be modified so that to be driven via inductances, as in ref. 46 . In addition, Fig. 1(d) presents two cantilevers, one of which is coupled to the optical cavity 13 ; and in Fig. 1(f) the two polarization modes of the light propagating in counter-clockwise (ccw) direction are shown to be tuned by the electro-optic modulators, EOM1 and EOM2, with the tuning parameter being the electric field inside the EOM1 47 .
To obtain the analogue of a driven qubit, following refs 9,12 , we assume 1,2 0 Here Δk(t) contains, in general, both dc and ac components. Then, introducing the interaction-shifted frequency (1) can be rewritten in the matrix form using the notation of the Pauli matrices σ x,z : Using the ansatz Instead of first directly solving these classical equations of motion for specific realizations (as e.g. in refs 5,10,16,46 ), we rather demonstrate the equivalence of these to the motion equation of a qubit, and then (in the next section) we will make use of the available solutions. We find this appropriate and pedagogical to first demonstrate the equivalence and then make use of the techniques and solutions available for qubit dynamics. Our approach differs only in details from the other authors' methods, though. For example, in ref. 12 the eigenfrequencies are found and only then the Bloch-like equation is obtained, while we do vice versa.
Equation (6) is simplified as follows. First, due to small dissipation, γ is neglected next to Ω 0 . Then, the slowly-varying envelope approximation 8,12 allows neglecting the second derivative (cf. ref. 43 ). This means that ψ n changes little during the time 2π/Ω 0 , i.e. its characteristic evolution frequency ω is much smaller than Ω 0 , ω Ω  0 . Then introducing new notations, from Eq. (6) we obtain In the absence of dissipation, γ = 0, equation (7a) formally coincides with the Schrödinger equation for a two-level system with the Hamiltonian H(t), applying ħ = 1.
Dissipation can be eliminated from the problem by the substitution ψ Alternatively, the "density matrix" can be introduced as ρ ψ ψ = , where ψ ψ ψ = ⁎ ⁎ : ( , ) 1 2 . Then for the derivative we obtain This coincides with the Bloch equation for a two-level system with the Hamiltonian H(t), assuming ħ = 1, and with equal relaxation rates, Introducing a convenient parametrization for the "Hamiltonian" and the "density matrix", 1 , x y z the "Bloch" Eq. (9) can be rewritten in the form of the Landau-Lifshitz-Bloch equation: where σ σ σ ε The diagonal components of the density matrix ρ and the Z-component of the Bloch vector X define the occupation of the respective states, while the off-diagonal components of ρ and X and Y describe the coherence.
We now see the analogy between the classical system and a qubit, described by the Bloch equation. This allows one to expect very similar dynamical phenomena. This was described in the introduction, while specific results are presented in the next two sections. After discussing this analogy, let us point out three key issues (see also, e.g. 11,12 ).
First, instead of different energy and phase relaxation rates for a generic quantum two-level system, for a classical analogue they coincide: . Recall that we have considered identical oscillators, with equal m, k 0 , and γ. In general, all of these quantities should be different. Thus, the equivalence drawn between T 1 and T 2 for the classical system so far is not general. They are only equivalent here because it is assumed that the damping for both oscillators is the same, which in general will not be the case, particularly for different frequency mechanical oscillators. In general, T 1 and T 2 are not related for classical systems, as they would be for a quantum system. For example, 1/T 2 = 1/2T 1 + 1/T φ for a quantum system, but not in the case of two classical coupled oscillators. For two oscillators with only damping rates different, γ 1,2 , we would have Second, a qubit at zero temperature relaxes to the ground state, defined by the Bloch vector X = (0, 0, 1), while the classical analogue in equilibrium relaxes to the zero Bloch-type vector X = (0, 0, 0), as can be seen from Eq. (11). This difference originates from the absence of a classical analogue to the purely quantum process of quantum emission 12 .
Third, one should remember the approximations done: when neglecting the second time derivative we assumed that the classical system's characteristic evolution frequency ω is much smaller than Ω 0 , i. e. ω Ω  0 . We note that the specific parameters depend on the choice of the system, as it was outlined in the introduction. To show specific numbers, we can take the ones close to ref. 10 for a nanomechanical two-mode beam: These parameters give the following: Ω  6 0 MHz ⋅ 2π, which indeed satisfies γ Ω  0 , then, Δ  7 kHz ⋅ 2π, which makes driving with ω ∼ Δ feasible, and also Δ ∼ .
k 0 03 N/m, which allows to discuss the regime of strong driving, with the amplitude ω Δ  A , .

Solutions of the Schrödinger-like equation
Using the analogy between the dynamical equations for the two coupled mechanical resonators and a two-level system, one can rewrite results from the respective publications, e.g. from refs 39,48,49 . For convenience, some analytical results are written down below, while detailed solutions are illustrated in the next Section. For the stationary Hamiltonian (7c) with the time-independent bias ε = ε 0 , we can transform from the functions ψ i to ψ ± , which define the eigenstates, analogously to the diagonalization of the Hamiltonian for qubits, e.g. ref. 39 : The solution of these two uncoupled equations is Together with the ansatz (5), this gives eigenfrequencies for the oscillations: .) The eigen-functions ψ ± are the amplitudes of the respective eigen-modes. These are analogous to the energy-level occupation amplitudes for qubits. Accordingly, we will be interested in the "occupation" of one of the modes, namely, |ψ + | 2 , which can be related to a qubit upper-level occupation probability. In experiments such value can be probed as an amplitude of the oscillations of the respective mode 10 . For example, if initially one, say "−", mode is excited, the problem can be formulated in finding the amplitude of the other mode; in this sense, the value |ψ + | 2 can be interpreted as a transition probability, describing the transition from one mode to another.
Consider now several regimes for the dc + ac driving: 0 For the single passage of the avoided crossing, we have the Landau-Zener problem, for which the solution is given by the probability This coincides with the result for the LZ-problem in refs 9,10 with the linear driving: v = α = Aω. Analogously, for the double-passage problem, we have the Stückelberg oscillations: is a parameter varying from −π/4, for δ  1, to −π/2, for δ  1 15 . We note that the integral above can be estimated, first, by neglecting Δ, and, second, for ε 0 = 0, then the integral is given by a special function: a full elliptic integral of the second kind 50 .
For the multiple-passage problem, consider here only the fast-passage (small Δ) limit. Then multi-photon Rabi oscillations are envisaged; next to the k-th resonance, where ε ω ∼ k 0 , we have For large arguments, the Bessel function has the oscillating asymptote If necessary, the damping is described by adding the factor exp(−γt) in Eqs (16,17). The time-averaged probability distribution is given by the series of Lorentzians: This, when plotted as a function of ε 0 and A, could serve as a visualization of the LZSM interference, which is the main subject of the next Section.
An analysis of the above equations allows for interpretations of specific systems. In our example of two coupled classical oscillators, we have for the Rabi frequency: R R (1) c which means that the energy transfer between the two oscillators, or between the two modes of a mechanical resonator, appears with a rate proportional to the coupling strength 4 .

Dynamics and interference in the classical two-state system
The Schrödinger-type and Bloch-type equations were presented above for two coupled mechanical resonators. This was done with several assumptions, which allowed reducing the original Eq. (4) to Eq. (7a). From the latter, a qubit-like behaviour follows. In this section we confirm and demonstrate this, by solving numerically the original equations (4). We consider a driven system with the bias ε(t) = ε 0 + ε 1 (t), where ε 1 (t) is one of the following: (1) the sinusoidal function ε 1 (t) = A sin ωt with weak or strong amplitudes, so-called Rabi and LZSM regimes, respectively; (2) the rectangular driving with ε 1 (t) = A sgn(sin ωt), so-called latching modulation; and (3) a noisy signal, corresponding to random jumps between ε 1 = +A and −A, with the characteristic switching frequency χ.
Rabi oscillations. Let us first consider the Rabi regime with a weak-amplitude sinusoidal driving. For calculations, we here choose the resonant frequency, ω = ω 0 with ε 0 = 5Δ, weak amplitude A = 0.7ω, and significant relaxation γ = 0.006ω, in order to see the damping of the Rabi oscillations. In Fig. 2, the thick red curve shows the numerical solution of the exact Eq. (6), the thin black curve is for the numerical solution of the approximate Schrödinger-like Eq. (7a), and the dashed blue curve depicts the analytical solution, Eq. (17). Similarly to their quantum counterparts, classical oscillations appear at the same conditions, of weak resonant driving, and have a similar expression for the Rabi frequency, Eq. (17). An important distinction is that the oscillations relax to zero, in contrast to the quantum case, where resonant oscillations result in a steady state with nonzero population of the excited state. As we can conclude from Fig. 2, for non-trivial results we should average before the system relaxes. This means that the analogue simulation of the quantum system has to be realized as the dynamics of the classical system on time scales Δt < γ −1 . In particular, for the resonant excitation, like the one in Fig. 2, after averaging the oscillating dynamics, we obtain ψ ∼ . LZSM interferometry. This regime occurs when a two-level system is strongly driven by a sinusoidal signal and monitored when changing its parameters; see ref. 39 and references therein. In such processes one is interested in the excitation probability, which increases periodically due to LZSM transitions between the states. After time averaging, the excitation probability displays increased values in the vicinity of the multi-photon resonances, where the energy difference between the states is matched by multiples of the driving frequency. Changing the system parameters, one can plot such interference fringes. This approach is useful both for controlling the system state and for defining the system parameters and its coupling to the environment 49,51 . To demonstrate this regime, we now consider the classical-oscillators system in the strong-driving regime, where, for a driven qubit in the same regime, LZSM interference takes place. Figure 3 shows the time-averaged mode occupation ψ + 2 as a function of the bias ε 0 and the driving amplitude A. With the numerical solution of the exact Eq. (6) for relatively high and low frequencies  ω Δ ( ) , we obtain the interferograms in Fig. 3(a,b), respectively. We can observe an important feature in these graphs: they display that the excitations appear at the position |ε 0 | = kω, with an integer k, analogously to the multi-photon transitions for qubits, as illustrated by the arrows in Fig. 1(a). These correspond to the minima of the denominator in Eq. (18), while the narrowing of these resonance lines appears around the zeros of the Bessel functions, entering in the numerators in Eq. (18). Note that this results in the interruptions of the resonance lines, where there are no transitions, even though the resonance condition, |ε 0 | = kω, is fulfilled. This phenomena for a quantum system is known as the coherent destruction of tunneling [52][53][54] .
Note that the interferograms in Fig. 3(a,b) are analogous to the ones obtained for diverse qubit systems, see ref. 39 and references therein. In particular, the interferograms in Fig. 3(a,b) are analogous to Fig. 7(b) and Fig. 8(b) in ref. 39 , respectively. See also ref. 55 where an analogous interferogram was recently calculated for a quantum-dot-electromechanical device.
Furthermore, Fig. 3(c) demonstrates the case of relatively stronger damping, when the interference fringes appear mostly due to the neighboring two transitions of the avoided crossing in Fig. 1(a). In this case, the resonances in Fig. 3(c) form characteristic V-shaped lines. We can observe that the outer lines are inclined along A = ±ε 0 , so that there is no excitation beyond these lines. This is because at small driving amplitudes, A < |ε 0 |, the avoided crossing in Fig. 1(a) is not reached, which explains the absence of the excitation. This regime of relatively strong damping with the V-shaped fringes was called quasiclassical in ref. 56 , which was studied for superconducting and semiconducting qubits in refs 40,56 .
In Fig. 3 we have chosen the following parameters for calculations: relatively weak damping in (a) and (b), γ = 0.02⋅ω/2π and γ = 0.1 ⋅ ω/2π, respectively, and stronger damping in (c), γ = ω/2π. The initial occupation was zero, ψ + (t = 0) = 0, like in Fig. 2, and then we averaged for the time interval γ Δ ∼ − t 1 . Namely, we took Δt = 50, 8, 1 ⋅ 2π/ω for the three panels in Fig. 3, respectively. If we choose γΔ  t 1 or γΔ  t 1, we would obtain similar data as in Fig. 3, but with a maximum amplitude closer to 0.5 or 0, respectively. After discussing this here, for the rest of the interferograms below we will assume γΔ  t 1. Alternatively, in addition to the above interferograms, one may be interested in the dependence on the driving frequency ω, as in refs 41,42 . We present such diagram in Fig. 4. Here the driving amplitude A is considered constant Latching modulation. If the system is driven by rectangular pulses instead of a sinusoidal signal, it experiences periodic fast changes between the two states. Namely, if the bias is ε(t) = ε 0 + Asgn(sin ωt), this means that the system is abruptly latched between the two limiting states with the bias given by either ε 0 + A or ε 0 − A, where it stays in each for about half of the period 41 . A conceptual distinction from the sinusoidal driving is in that the avoided-level crossing, where LZSM transitions appear, is crossed rapidly. That is why the theory developed for smooth sinusoidal driving had to be revisited for the rectangular-pulse driving, which was done in ref. 41 both theoretically and experimentally for a superconducting qubit. Here, we study such a formulation for our classical system: two coupled classical oscillators.
To describe this regime, we have solved Eq. (6), as in the previous subsection, but now with a different bias. Figure 5 shows the time-averaged mode occupation ψ + 2 as a function of the bias ε 0 and the driving frequency ω. Such latching modulation displays interference fringes, different from the ones shown above, for a harmonic driving. Interestingly, at low frequencies the system is indeed mostly latched to the two states with the resonances around ε 0 = ±A, while for higher frequencies there is no trace of the latching, and the position of the resonances is described by the inclined resonance lines. This means that, due to the interference, our system latching between ε 0 = ±A displays resonances for any other values of ε 0 . Note that here, for the classical system, we obtained a remarkable agreement with the diagram obtained recently for the experimental qubit in ref. 41 .
Motional averaging. Differently from above, in the regime which we refer to as the motional averaging, the system is driven by noise rather than by a periodic signal. The rapid jumps between the two states appear stochastically. For this we follow ref. 42 , where it was demonstrated for a qubit that when increasing the frequency of the jumps, two separate spectral lines merge into one line. This observation allowed the authors of ref. 42 to speculate  2 is plotted when the system is driven by rectangular pulses. Note that the analogous interferogram for a qubit was experimentally observed in ref. 42  about the analogue simulation of this so-called motional averaging, originally observed in nuclear magnetic resonance spectroscopy.
Accordingly, we consider here the classical two-state system driven by a non-periodic signal. Namely, we take the bias ε(t) = ε 0 + ε 1 (t), admitting random jumps between ε 1 = +A and −A. The jumps are assumed to appear with the average jumping rate χ. In Fig. 6 we plot the time-averaged mode occupation ψ + 2 as a function of the detuning ε 0 and the characteristic switching frequency χ.
Our data shown in Fig. 6 are consistent with the result of ref. 42 in that, at low χ, there are two peaks, which merge into one, for increasing χ. To make this point explicit, in Fig. 6(b) we plot the two cross-sections of Fig. 6(a) along the horizontal white dashed lines. The red curve in Fig. 6(b) is plotted for a low switching frequency, and this displays two distinct peaks at ε 0 = ±A. The blue curve, shifted vertically for clarity, is plotted for a relatively high switching frequency. This displays a main peak aroung ε 0 = 0, which can be interpreted as the averaging due to the motion between the two states.

Conclusions
Remarkably, the classical system of two weakly coupled classical oscillators form the doppelgänger of a quantum two-level system. Namely, its equation of motion formally coincides with either the Schrödinger equation or with the Bloch equation in the cases when the relaxation is either ignored or taken into account, respectively. This means that the dynamical phenomena of the two-state classical system can be directly described by the ones already studied for the quantum two-level systems, and vice versa. This was known and studied some time ago, e.g. in refs 4,5,7,47 . However, diverse experimental mechanical resonators which are good enough for analogue simulations appeared only recently 11,13,43,44 . Such mechanical resonators have high quality factors and have reliable control of their inter-mode coupling. Moreover, there is recent interest in strongly-driven quantum systems, e.g. refs 39,41,42 . These two developments stimulated us to further consider the analogy between weakly and controllably coupled mechanical oscillators and the driven quantum few-level system. In particular, we demonstrated classical analogues of the effects recently studied for qubits: Landau-Zener-Stückelberg-Majorana interferometry, latching modulation, and motional averaging. Besides the pure interest of such dynamical phenomena linking classical and quantum physics, one may consider simulating some quantum phenomena with classical systems. The time-averaged mode occupation ψ + 2 as a function of the bias ε 0 and the jumping rate χ; this displays two peaks at around ε 0 = ±A for a slowly jumping signal with χ  A, while the faster jumping, with  χ A, results in the merging of the two peaks into one. (b) The time-averaged mode occupation as a function of the bias ε 0 for the two values of the jumping rate χ. Note that analogous dependencies were experimentally demonstrated for a superconducting qubit in ref. 42 .