Generic shape of multichromatic resonance peaks

Depósito de investigación de la Universidad de Sevilla https://idus.us.es/ “This version of the article has been accepted for publication, after peer review (when applicable) and is subject to Springer Nature’s AM terms of use, but is not the Version of Record and does not reflect post-acceptance improvements, or any corrections. The Version of Record is available online at: https://doi.org/10.1140/epjb/e2020-100595-0 ” Eur. Phys. J. B manuscript No. (will be inserted by the editor) Generic shape of multichromatic resonance peaks Mar´ıa Laura Olivera-Atencio1, Jes´us Casado-Pascual1, and Sigmund Kohler2 1F´ısica Te´orica, Universidad de Sevilla, Apartado de Correos 1065, 41080 Sevilla, Spain 2Instituto de Ciencia de Materiales de Madrid, CSIC, 28049 Madrid, Spain January 15, 2020 Abstract. We investigate dissipative dynamical systems under the influence of an external driving with two or more frequencies. Our main quantities of interest are long-time averages of expectation values which turn out to exhibit universal features. In particular, resonance peaks in the vicinity of commensurable frequencies possess a generic enveloping function whose width is inversely proportional to the averaging time. While the universal features can be derived analytically, the transition from the specific short-time behavior to the long-time limit is illustrated for the examples of a classical random walk and a dissipative two-level system both with biharmonic driving. In these models, the dependence of the time-averaged response on the relative phase between the two driving frequencies changes with increasing integration time. For short times, it exhibits the 2πperiodicity of the dynamic equations, while in the long-time limit, the period becomes a fraction of this value. PACS. 1 Introduction Dynamical systems driven by time-dependent forces represent paradigmatic models for non-equilibrium effects in the realm of both classical and quantum mechanics. There one may find counter-intuitive phenomena such as stochastic resonance [1,2], synchronization [3–7], and the ratchet effect by which directed currents emerge despite the absence of any net force [8–10]. The common physics behind these phenomena is an interplay of non-linearities and non-equilibrium. Driven systems have been studied in the Hamiltonian limit [11–13], as well as in the steady state of dissipative classical [8,9,14,15] and quantum mechanical models [16–18]. Even far from equilibrium, spatio-temporal symmetries may inhibit the emergence of a dc response such as a ratchet current. Many of these symmetries are based on the fact that a sinusoidal driving force changes its sign when time is shifted by half a period. For bichromatic or multichromatic forces, the periods of the various forces are different and, thus, the symmetry analysis has to be revised. Typically, one has to distinguish two cases, namely those of commensurable and incommensurable frequencies. In the former case, the symmetry depends on the phase between the components of the driving, which has been verified with transport experiments in quantum dots [19, 20]. In reference [19], it has also been demonstrated theoretically and experimentally that for incommensurable frequencies, the symmetry may be higher, despite that the driving force is only quasi periodic and may not possess any symmetry or anti-symmetry. Here, we work out generic properties of dissipative dynamical systems under multi-frequency driving, in particular the shape of the resulting resonance peaks in the longtime limit. Moreover, we study how this limit emerges. The article is organized as follows. In Section 2, we formulate the problem and in Section 3 derive how the independence on the initial time provides generic properties. In Section 4, we consider finite-time effects for bichromatic driving and in Sections 5 and 6 show how the limits are approached for a classical random walk model and for a quantum mechanical two-level system, respectively. Finally, we summarize and conclude in Section 7. 2 Formulation of the problem Suppose that the dynamical equations governing the time evolution of the system under consideration depend on time through Ntime-periodic functions of the form fj(t) = Φj(Ωjt+ϕj),(1) with j= 1, . . . , N. In the above expression, Φjare 2πperiodic functions [i.e., Φj(θ+ 2π) = Φj(θ)∀θ∈R] and Ωjand ϕjdenote the angular frequency and the phase of fj(t), respectively. The precise nature of the functions fj(t) is irrelevant for our purposes, e.g., they may represent external oscillatory forces, modulating amplitudes, etc. Assume that the system has been prepared in a state s0 at an initial time t0. Depending on the case, s0may represent the system’s density operator (in the case of quantum 2 Mar´ıa Laura Olivera-Atencio et al.: Generic shape of multichromatic resonance peaks systems), the one-time probability density (in the case of classical stochastic systems), the values of a finite number of state variables or classical fields (in the case of classical deterministic systems), etc. For our purpose, the only relevant fact is that all the physical properties of the system at any subsequent time t≥t0are uniquely determined by this initial preparation and the dynamical equations. In particular, the (expectation) value at time t≥t0of a certain physical quantity Qof this system, denoted by Qt, will depend on the specific values taken by the parameters appearing in the functions fj(t), as well as, on the initial preparation. When necessary, this dependence will be made explicit by the notation Qt(Ω,ϕ|s0, t0) with the vector notation Ω= (Ω1, . . . , ΩN) and ϕ= (ϕ1, . . . , ϕN). The transformation property of Qt(Ω,ϕ|s0, t0) is determined by the fj, as for any time shift τ, the set of periodic functions {fj(t)}j=1,...,N is invariant under the transformation {t, ϕ} 7→ {t+τ, ϕ−Ωτ}. Consequently, since the only explicit time-dependence of the dynamical equations comes from the functions fj(t), it readily follows that Qt(Ω,ϕ|s0, t0) = Qt+τ(Ω,ϕ−Ωτ|s0, t0+τ).(2) In particular, taking τ=−t, one obtains Qt(Ω,ϕ|s0, t0) = Q0(Ω,ϕ+Ωt|s0, t0−t).(3) Henceforth, we assume that, after a certain transient time (or, more formally, in the limit t0→ −∞), the observable under investigation reaches a stationary value Qst t(Ω,ϕ) = lim t0→−∞Qt(Ω,ϕ|s0, t0).(4) This assumption has been found in many dissipative systems. From equation (3), it then follows that Qst t(Ω,ϕ) = Qst 0(Ω,ϕ+Ωt).(5) Therefore, in the stationary regime, the time evolution of Qadmits a description in terms of a time-dependent phase vector of the form ϕ+Ωt. Here, we are interested in the generic properties of the time average QT(Ω,ϕ) = 1 TZT 0 dt Qst t(Ω,ϕ),(6) and more specifically on the infinite-time average Q∞(Ω,ϕ) = lim T→∞QT(Ω,ϕ).(7) Since the functions Φj(Ωjt+ϕj) are 2π-periodic in the phases ϕj, so will be Qst t(Ω,ϕ). Therefore, taking into account equation (5), it can be Fourier expanded as Qst t(Ω,ϕ) = X k∈ZN qk(Ω)eik·(ϕ+Ωt),(8) where the centered dot denotes the usual scalar product in RNand qk(Ω) = Zπ −π···Zπ −π e−ik·ϕ (2π)NQst 0(Ω,ϕ)dNϕ.(9) By inserting equation (8) into equation (6), we obtain QT(Ω,ϕ) = X k∈ZN qk(Ω)eik·(ϕ+ΩT/2)sinc k·ΩT 2, (10) where sinc(x) = sin(x)/x denotes the unnormalized sinus cardinalis. Since limT→∞ sinc (αT/2) equals 1 if α= 0 and 0 otherwise, from equations (7) and (10) we get Q∞(Ω,ϕ) = X k∈SΩ qk(Ω)eik·ϕ,(11) where SΩis the set of all ordered N-tuples of integers orthogonal to Ω, i.e., SΩ=k∈ZN:k·Ω= 0. If the Ncomponents of the frequency vector Ωare incommensurable (i.e., it is not possible to express one of them as a linear combination of the others with rational coefficients), then the set SΩreduces to the single element k=0. In this case, from equation (11), it readily follows that Q∞(Ω,ϕ) = q0(Ω) (12) and, therefore, the infinite-time average is independent of ϕ. An experimental verification of this general statement has been reported in reference [19]. In contrast, if the Ncomponents of Ωare commensurable (i.e., it is possible to express one of them as a linear combination of the others with rational coefficients), then the set SΩcontains additional elements other than k=0 and, according to equation (11), Q∞(Ω,ϕ) depends on the phases in ϕ. In particular, if the frequencies are pairwise commensurable, then there exists a frequency Ωsuch that Ω=Ωn, where n= (n1, . . . , nN), with njbeing positive integers. Consequently, the condition k·Ω= 0 becomes equivalent to the Diophantine equation k·n= 0. The general solution of the latter equation can be expressed as an integer linear combination of a set of N−1 generating vectors, g(1),...,g(N−1), each of which satisfies the equation g(j)·n= 0 [15, 21, 22]. Thus, equation (11) can be written as Q∞(Ω,ϕ) = X `∈ZN−1 qk(`)(Ω)eik(`)·ϕ,(13) where k(`) = PN−1 j=1 `jg(j). Let WCand WIbe the sets of all Ωwhose components are commensurable and incommensurable, respectively. From the above results, it can easily be shown that the infinite-time average Q∞(Ω,ϕ) is discontinuous on the set WC. Indeed, since the set WIis dense in RN, for any Ωc∈ WCone can find a sequence {Ω(n)}n∈N⊂ WI such that limn→∞ Ω(n)=Ωc. If Q∞(Ω,ϕ) were continuous at Ωc, then Q∞(Ωc,ϕ) = lim n→∞Q∞(Ω(n),ϕ) = lim n→∞q0(Ω(n)) (14) and, as a result, Q∞(Ωc,ϕ) would be independent of ϕ. Obviously, this contradicts the fact that, since Ωc∈ WC, Q∞(Ωc,ϕ) depends on ϕ. Hence, we can conclude that Mar´ıa Laura Olivera-Atencio et al.: Generic shape of multichromatic resonance peaks 3 Q∞(Ω,ϕ) is discontinuous at Ωcand, thus, on the set WC. Taking into account that WCis also dense in RN, this last result implies that the infinite-time average Q∞(Ω,ϕ) is a highly discontinuous function of the frequency vector Ω. In this context, the following questions arise: (i) Can this discontinuity actually be observed? (ii) How does it manifest itself in practice? 3 Long-time asymptotic behavior of QT In real situations, the physical quantity Qis known only in a finite time-interval. Therefore, the infinite time-average Q∞can be calculated only approximately by taking a sufficiently large value of T. Suppose we are interested in analyzing the dependence of QTon Ωnear a fixed frequency vector Ω0. More precisely, we focus on values of Ωsuch that |Ω−Ω0|is of the same order of magnitude as T−1, where |Ω−Ω0|= [PN j=1(Ωj−Ω0,j)2]1/2. Then, as Tincreases, the size of the region of interest becomes smaller in inverse proportion to T. By defining the dimensionless frequency vector δ˜ω= (Ω−Ω0)T, the finite time-average can be brought to the form QT(Ω,ϕ) = QT(Ω0+δ˜ω/T, ϕ) =Qas(Ω0, δ˜ω,ϕ) + RT(Ω0, δ˜ω,ϕ),(15) where Qas(Ω0, δ˜ω,ϕ) = lim T→∞QT(Ω0+δ˜ω/T, ϕ) (16) is the leading-order of the asymptotic behavior of the function QT(Ω0+δ˜ω/T, ϕ) for T→ ∞, while δ˜ωis held constant. For the rest RT(Ω0, δ˜ω,ϕ) = QT(Ω0+δ˜ω/T, ϕ)− Qas(Ω0, δ˜ω,ϕ), it readily follows that lim T→∞RT(Ω0, δ˜ω,ϕ)=0,(17) and together with equation (7), that Qas(Ω0,0,ϕ) = Q∞(Ω0,ϕ).(18) To obtain an expression for Qas(Ω0, δ˜ω,ϕ), we set Ω=Ω0+δ˜ω/T in equation (10), and insert the resulting expression into equation (16). Then, using again limT→∞ sinc (αT/2) = δα,0, one obtains Qas(Ω0, δ˜ω,ϕ) = X k∈SΩ qk(Ω0)eik·(ϕ+δ˜ω/2)sinc k·δ˜ω 2. (19) Using equation (11), one sees that equation (19) can be written in the more compact form Qas(Ω0, δ˜ω,ϕ) = Z1 0 dλ Q∞(Ω0,ϕ+λδ˜ω).(20) Therefore, in the limit T→ ∞, the asymptotic behavior of QTin a neighborhood of Ω0of size proportional to T−1is completely determined by the infinite-time limit Q∞at Ω0. Notice that, according to the above results, no matter how large T, it is always possible to find a neighborhood of Ω0within which QTis continuous. Thus, the discontinuity mentioned in the previous section is a mathematical idealization that, owing to the necessarily finite measurement time, cannot be observed in reality. From equation (20) there immediately follows an interesting conclusion. If Ω0∈ WI(i.e., all Ncomponents are incommensurable), then Q∞(Ω0,ϕ) is independent of ϕand Qas(Ω0, δ˜ω,ϕ) = Q∞(Ω0,ϕ) for all δ˜ω. Consequently, from equations (7) and (16), one obtains that lim T→∞QT(Ω0+δ˜ω/T, ϕ)−QT(Ω0,ϕ)= 0 (21) for all δ˜ω. In contrast, if Ω0∈ WC(i.e., its Ncomponents are commensurable), then Q∞(Ω0,ϕ) depends on ϕand, in general, lim T→∞QT(Ω0+δ˜ω/T, ϕ)−QT(Ω0,ϕ)6= 0.(22) This difference in behavior will be apparent in the numerical calculations presented below. 4 Bichromatic driving The analytic considerations made so far provide the longtime limit of the non-linear response to a multi-frequency forcing. In order to investigate numerically how this limit is approached, we focus on bichromatic driving, i.e., on dynamic equations which contain terms of the form f1(t) = A1cos(Ω1t+ϕ1) and f2(t) = A2cos(Ω2t+ϕ2), where our line of reasoning still holds if the cosines are replaced by any other 2π-periodic functions. 4.1 Commensurable frequencies and periodicity in ϕ While the 2π-periodicity of {f1(t), f2(t)}in ϕ1and ϕ2 is evident, discrete time translation symmetry is present only when Ω1and Ω2are commensurable, i.e., for rational values of Ω2/Ω1. Indeed, let us assume that f1(t) and f2(t) have a common fundamental period T. From the periodicity of the cosine function, it then follows that there must exist two integers qand psuch that Ω1T= 2πq and Ω2T= 2πp. This implies that Ω1=qΩ and Ω2=pΩ, with Ω= 2π/T, and consequently, that Ω2/Ω1=p/q. In addition, the integers qand pmust be coprime because otherwise T/gcd(q, p)<T, would be the common fundamental period of f1(t) and f2(t), with gcd(q, p) denoting the greatest common divisor of qand p. Let us now focus on a frequency vector Ω0≡(q, p)Ω, with some coprime integers pand q, henceforth referred to as (q, p)-resonance. Using equation (8), it can easily be seen that the stationary value Qst t(Ω0,ϕ) is a T-periodic function of time. A further interesting fact is that the infinite-time average Q∞(Ω0,ϕ) displays a higher symmetry in ϕthan the obvious ϕ1→ϕ1+ 2πand ϕ2→ 4 Mar´ıa Laura Olivera-Atencio et al.: Generic shape of multichromatic resonance peaks ϕ2+ 2π[14]. To see that this is so, note that in the present case the condition k·Ω0= 0 becomes equivalent to k1q+k2p= 0. The general solution of this Diophantine equation is `g, with the generating vector g≡(−p, q) and any integer `. Thus, in this case, equation (13) becomes Q∞(Ω0,ϕ) = X `∈Z q`g(Ω0)ei`g·ϕ.(23) Since g·ϕ=−pϕ1+qϕ2, from equation (23) it readily follows that Q∞(Ω0,ϕ) is 2π/p-periodic in ϕ1and 2π/qperiodic in ϕ2. 4.2 The neighborhood of commensurable frequencies Let us now turn our attention to the vicinity of a (q, p)- resonance. More specifically, and following the general approach outlined in Section 3, we focus on values of Ω such that |Ω−Ω0|= [(Ω1−qΩ)2+ (Ω2−pΩ)2]1/2is of the same order as T−1. Then, introducing the notation δω=δ˜ω/T =Ω−Ω0, from equation (15), (17), and (19), it follows that QT(Ω0+δω,ϕ)∼X `∈Z q`g(Ω0)ei`g·(ϕ+T δω/2) ×sinc `Tg·δω 2,(24) provided that Tis large enough so that the term RTcan be neglected. Two conclusions can be drawn from equation (24). First, the time average QT(Ω0+δω,ϕ) is 2π/p-periodic in ϕ1and 2π/q-periodic in ϕ2. There is, however, an important difference to the exact periodicity of the infinite-time average Q∞(Ω0,ϕ) presented in Section 4.1. Here the periodicity holds only for sufficiently large values of Tand for Ωsufficiently close to Ω0—to be precise, for |δω|= O(T−1). Second, under these same restrictions, the relative height ∆QT(Ω0+δω,ϕ)≡QT(Ω0+δω,ϕ)−q0(Ω0) vanishes when all sinc functions in equation (24) are zero for all `6= 0, i.e., when g·δωis a nonzero integer multiple of 2π/T. In particular, if we set δω1= 0 and vary δω2, ∆QT(Ω0+δω,ϕ) vanishes when δω2is a nonzero integer multiple of 2π/(qT). Analogously, if we set δω2= 0 and vary δω1,∆QT(Ω0+δω,ϕ) vanishes when δω1is a nonzero integer multiple of 2π/(pT). Hence, in the vicinity of the (q, p)-resonance, 2π/(qT) and 2π/(pT) represent the frequency scales on which ∆QT(Ω0+δω,ϕ) varies. The main aim of the next two sections is to provide quantitative information on how these asymptotic properties can actually be observed in practice. Note, for example, that a further limitation may stem from the fact that the set of commensurable frequencies vectors is dense in R2and, consequently, any (q, p)-resonance may be disturbed by another resonance that lies arbitrarily close to it. It will turn out, however, for all cases investigated, practically only resonances with rather small qand pmatter. Moreover, for sufficiently large values of T, the shape of the resonance peak is mainly governed by the term with `= 1 in equation (24) and, therefore, their enveloping function appears as a rather clean sinc function. 5 Classical random walk model In the model presented in this section, the motion of a particle in a periodic substrate is given by a random walk on a one-dimensional lattice. The lattice sites are located at xn=na, where nis any integer and athe distance between two neighboring sites. The evolution of the probabilities pn(t) that the particle is at site nis governed by the master equation [23] ˙pn(t) = −[r+(t) + r−(t)] pn(t) +r+(t)pn−1(t) + r−(t)pn+1(t),(25) where r+(t) and r−(t) are the transition rates from site n to site n+1 and n−1, respectively. They are assumed to be independent of nand to follow the Van’t Hoff-Arrhenius law [24] r±(t) = r0e−β[E0±∆E(t)],(26) where r0is a prefactor with the dimension of an inverse time, β= (kBΘ)−1is the inverse temperature, and E0+ ∆E(t) and E0−∆E(t) are, respectively, the activation energies for the forward and backward steps. These activation energies oscillate around a constant value E0with time-dependent amplitudes ∆E(t) and −∆E(t), respectively. In this section, we will consider that the role of Qtis played by the mean particle velocity Vt, which is defined as the time derivative of the mean particle position Xt= aPn∈Znpn(t). Using equations (25) and (26), it can be seen easily that Vt=a[r+(t)−r−(t)] = vsinh [f(t)] ,(27) where v= 2ar0e−βE0and f(t) = −β∆E(t). Henceforth, we will assume that f(t) = A1cos(Ω1t+ϕ1) + A2cos(Ω2t+ϕ2),(28) with A1and A2being two dimensionless constants, which can be taken as positive by suitable choice of the phases ϕ1 and ϕ2. Note that, in the present model, there is no difference between the stationary Vst tand Vtbecause the mean particle velocity is independent of the initial preparation. In addition, from equations (27) and (28), it immediately follows that Vtsatisfies the symmetry properties V−t(Ω,−ϕ) = Vt(Ω,ϕ),(29) Vt(Ω,ϕ+π) = −Vt(Ω,ϕ),(30) where π≡(π, π). Now the Fourier expansion in equation (8) takes the form Vt(Ω,ϕ) = X k∈Z2 vkeik·(ϕ+Ωt),(31) Mar´ıa Laura Olivera-Atencio et al.: Generic shape of multichromatic resonance peaks 5 with vk=vZπ −πZπ −π e−ik·ϕ (2π)2sinh (A1cos ϕ1+A2cos ϕ2)d2ϕ. (32) Since the Fourier expansion is unique, using equations (29) and (31), it is easy to see that v−k=vk. For the same reason, from equations (30) and (31), it readily follows that vk=−ei(k1+k2)πvk. In addition, from equation (32), it is clear that v∗ k=v−k, where the asterisk denotes complex conjugation. In conclusion, all the coefficients vkare real and vanish when k1+k2is an even integer. With the above results, let us now examine the dependence of the infinite-time average velocity on Ω. If Ω1 and Ω2are incommensurable, from equation (12) one concludes that V∞(Ω,ϕ) = v0= 0. If, by contrast, Ω1and Ω2are commensurable, taking into account that vk=v−k and that v2k= 0, it is easy to see from equation (23) that V∞(Ω0,ϕ)=2 ∞ X `=0 v(2`+1)gcos [(2`+ 1)g·ϕ].(33) Using the definition of gand the fact that vkvanishes when k1+k2is even, it can be verified with equation (33) that V∞(Ω0,ϕ) = 0 when q−pis even. Therefore, as a consequence of the symmetry property (30), only the (q, p)-resonances with q−podd are present in this model. The integral in equation (32) can be evaluated explicitly by expanding the hyperbolic sine into a power series. Then, after some calculations, we obtain vk=v∞ X j=0 2j+1 X `=0 ` X m=0 2j+1−` X n=0 δk1,2m−`δk2,2n+`−2j−1 ×A` 1A2j+1−` 2 22j+1m!n!(`−m)!(2j+ 1 −`−n)!.(34) It can be verified that the coefficients vkgiven by equation (34), as could not be otherwise, satisfy the conditions discussed above. In addition, since A1and A2are positive, all the coefficients vkare clearly non-negative. Thus, from equation (33), it follows that the maximum value of V∞(Ω0,ϕ) is V∞,M(Ω0)=2 ∞ X `=0 v(2`+1)g,(35) and it occurs when g·ϕis an integer multiple of 2π, i.e., when qϕ2−pϕ1= 2πn (36) with nbeing any integer. Assuming that ϕ1and ϕ2satisfy this condition of maximum resonance, from equation (24) together with v−k=vk, it follows that VT(Ω0+δω,ϕ)∼2∞ X `=0 v(2`+1)gsinc [(2`+ 1)Tg·δω], (37) Fig. 1. Dependence of the dimensionless average velocity VT/v on Ω2/Ω1for A1=A2= 1, ϕ1=ϕ2= 0, and T= 104/Ω1. The resonance peaks corresponding to the fractions 0 (only partially shown for better visibility of the remaining peaks) and 1/2 are clearly visible. The ones corresponding to 1/4 and 2/3 are also visible but considerably smaller. Other resonances, such as the ones corresponding to 2/5 and 6/7 shown in Figure 3 cannot be appreciated on this scale. The insets show a zoomed-in view around Ω2/Ω1= 1/4 (left inset) and Ω2/Ω1= 2/3 (right inset) in terms of δω/Ω1=Ω2/Ω1−p/q. provided that |δω|=O(T−1) and that Tis large enough such that the term RTcan be neglected. Note that, in this case, VT(Ω0+δω,ϕ) vanishes whenever g·δωis a nonzero integer multiple of π/T. Thus, if we set δω1= 0 and vary δω2,VT(Ω0+δω,ϕ) vanishes when δω2is a nonzero integer multiple of π/(qT), whereas if we set δω2= 0 and vary δω1,VT(Ω0+δω,ϕ) vanishes if δω1is a nonzero integer multiple of π/(pT). To illustrate our theoretical results, we have calculated VT(Ω,ϕ) using equation (10) with QT(Ω,ϕ) = VT(Ω,ϕ), and qk(Ω) = vk. The coefficients vkappearing in that equation have been evaluated using equation (34). To ensure that both drivings have roughly the same influence, we have focused on the most symmetric case of equal driving amplitudes. Specifically, in all the figures of this section we have taken A1=A2= 1. In addition, to maximize the height of the resonance peaks, we have restricted our analysis to values of ϕ1and ϕ2that satisfy the condition of maximum resonance in equation (36). In Figure 1, we plot the dimensionless time-average velocity VT/v as a function of Ω2/Ω1for ϕ1=ϕ2= 0 and T= 104/Ω1. According to our theoretical results, there should emerge resonance peaks when Ω2/Ω1is equal to a rational number. Furthemore, only the resonances corresponding to the fractions 0, 1/4, 1/2, and 2/3 are visible. Other resonances, such as the ones corresponding to 2/5 and 6/7 studied below, cannot be appreciated in the figure. This absence of resonances is not due to the use of a finite averaging time. In fact, these peaks would be imperceptible even in the limit T→ ∞, because they are very small compared to the smallest peak visible in Fig- 6 Mar´ıa Laura Olivera-Atencio et al.: Generic shape of multichromatic resonance peaks ure 1. Indeed, using equation (35), it is easy to verify that at Ω2/Ω1= 1/4 (the smallest peak visible in Figure 1), V∞,Mis approximately equal to 3 ×10−3v, whereas at Ω2/Ω1= 2/5 and Ω2/Ω1= 6/7, V∞,Mis approximately equal to 7 ×10−5vand 7 ×10−11v, respectively. To analyze in more detail the behavior of VTin the vicinity of a (q, p)-resonance, we now consider that one frequency, say Ω1, is kept fixed, while the other frequency, Ω2, is varied around the value pΩ1/q. In terms of our previous notation, this corresponds to setting δω1= 0 and δω2=Ω2−pΩ1/q. Henceforth, for notational simplicity, we will write δω instead of δω2. In addition, to facilitate comparison with the asymptotic expression (37), we will use the dimensionless variable qδωT, which here is nothing but Tg·δω. In Figures 2 and 3, we depict the dependence of VT/v on qδωT for the resonances (q, p) = (4,1) (left column in Figure 2), (q, p) = (3,2) (right column in Figure 2), (q, p) = (5,2) (left column in Figure 3), and (q, p) = (7,6) (right column in Figure 3). To analyze how these peaks emerge as the averaging time increases, different values of Thave been considered, which are indicated in the panels. In addition, in each panel, qcurves have been plotted, corresponding to the values ϕ2= 2πn/q, with n= 0, . . . , q −1. Since ϕ1= 0, all these values satisfy the condition of maximum resonance in equation (36). The results in Figures 2 and 3 corroborate the theoretical prediction that, for sufficiently large values of T, the q curves converge to the asymptotic result in equation (37). Surprisingly, the averaging times necessary to reach the asymptotic regime are huge in comparison to the period of the driving T= 2πq/Ω1. 6 Quantum mechanical two-level system Let us now consider a dissipative quantum mechanical model with the Hamiltonian H(t) = H0+HD(t), where H0= 2(σxcos θ+σzsin θ) (38) and the bichromatic driving HD(t) = A1σxcos(Ω1t+ϕ1) + A2σzcos(Ω2t+ϕ2). (39) To ensure that both drivings have roughly the same impact, we focus on the most symmetric case θ=π/4 and equal driving amplitudes, A1=A2≡A. For the consideration of dissipation, one may start from a system-bath model to obtain an equation of motion for the reduced density operator of the dissipative system. Then one can show that generally dissipation is quantitatively affected by the driving [25, 26]. Here however, we are interested in the generic response to bichromatic drivings and, thus, we follow a less involved path which allows an efficient numerical solution for rather long propagation times. In doing so, we employ a Lindblad master equation for the density operator [27], ˙ρ=−i[H0+HD(t), ρ]+ Fig. 2. Dependence of the dimensionless average velocity VT/v on the dimensionless variable qδωT for the resonances (q, p) = (4,1) (left column) and (q, p) = (3,2) (right column), and the values of the averaging times displayed in the panels. For all the curves ϕ1= 0 and A1=A2= 1. In each panel, there are qcurves corresponding to the values ϕ2= 2πn/q, for n= 0,...,q−1 (from blue to red), which are obtained from the condition of maximum resonance in equation (36). As expected from the theoretical analysis, with increasing the value of T, the qcurves converge to the asymptotic curve given by equation (37). Fig. 3. The same as in Figure 2 but for the resonances (q, p) = (5,2) (left column) and (q, p) = (7,6) (right column). γD(ρ) (in units with ~= 1) with a dissipator [27] Dρ= ˜σ−ρ˜σ+−1 2˜σ+˜σ−ρ−1 2ρ˜σ+˜σ−,(40) where ˜σ−=|ϕ0ihϕ1|=σ† +induces dissipative decay from the excited state |ϕ1iof the undriven Hamiltonian H0to the corresponding ground state |ϕ0i. As an observable Q, we may choose any combination of Pauli matrices. Generic features, however, will not depend on the particular choice such that without loss of generality, we consider Qst t≡ hσzit. To ensure independence Mar´ıa Laura Olivera-Atencio et al.: Generic shape of multichromatic resonance peaks 7 0.5 0.6 0.7 0.8 0.911.1 −0.5 −0.25 0 0.25 Ω2/Ω1 hσziT T= 500/Ω1 T= 1000/Ω1 1 1 1 2 4 73 5 2 3 5 7 3 4 5 6 Fig. 4. Time-averaged response of the bichromatically driven two-level system after a transient stage averaged over various times T. When Ω2/Ω1is close to a rational p/q, peaks with generic shape emerge (labeled by the corresponding fraction). Parameter values are: =√2Ω1,θ=π/4, A1=A2=Ω1, and γ= 0.2Ω1. For graphical reasons, the curve for T= 500/Ω1is vertically shifted. of details, we have verified all numerical results by using also slightly different setups and observables. The calculations are performed by numerically integrating the Lindblad master equation starting at time t= 0 in the ground state of H0to obtain the density operator ρtand, thus, Qt≡tr(ρtσz). Since we are interested in stationary expectation values, we exclude the transient stage and compute the time averages in an interval of duration Tstarting at time 10/γ. We start by sketching the global picture of the response (Figure 4) which shows hσziTfor two different averaging times Tas a function of Ω2and for various phases ϕ2 (in this section, we always take ϕ1= 0). As expected, hσziTexhibits resonance peaks at simple rational values of Ω2/Ω1which sharpen with increasing T. The size of the peaks as well as the shape of the background depend on details. Generally it is such that the (1,1) resonance is rather prominent, which can be understood by frequency mixing: For two equal driving frequencies, a zero-frequency response can be obtained already in second-order perturbation theory. For all other values of (q, p), one has to go to higher order. Henceforth we focus on the universal features in narrow regions around pronounced peaks. As smaller peaks tend to be less compromised by components with `6= 1, we study the emergence of a peak with increasing propagation time Tfor (q, p) = (5,3). Figure 5 shows hσziTfor the 20 equally spaced phases ϕ2= 2πn/20 with n= 0,1,...,19. For the relatively short propagation time T= 150/Ω1, all curves are significantly different from each other, while a clear peak structure is still missing. With increasing Tand staring at the center δω = 0, curves for phases that differ by 2π/q coincide, such that eventually 20/q = 4 groups of curves emerge. This reflects the 2π/q periodicity in ϕ2derived above for δω = 0. Moreover, it confirms the generalization conjectured from equation (24), namely that the periodicity to −4π−2π0 2π4π −0.5 −0.475 q δω T hσziT T= 5000/Ω1 −0.5 −0.475 hσziT T= 1000/Ω1 −0.5 −0.475 hσziT T= 750/Ω1 −0.5 −0.475 hσziT T= 250/Ω1 −0.5 −0.475 hσziT T= 150/Ω1 Fig. 5. Resonance peak (q, p) = (5,3) for the phases ϕ1= 0 and ϕ2= 2πn/(4q), with n= 0,1, .., 4q−1 (from blue to red), showing the transition from 2π-periodicity to 2π/q-periodicity. Notice that the abscissa is scaled with the averaging time T. All other parameters are as in Figure 4. a good approximation holds in a whole neighborhood of the (q, p)-peak. To underline this result, we also plot the curves within one 2π/q period and those for ϕ2equal to multiples of 2π/q separately, but now for the (q, p) = (7,5) resonance, see Figure 6. The left column contains the results for values of ϕ2in the range [0,2π/q]. They show that only for a sufficiently large T, the enveloping function is dominated by the `= 1 component and resembles the sinc conjectured in equation (24). Moreover, the first and the last curve smoothly connect to each other, which depicts how the 2π/q-periodicity emerges. Accordingly, the curves for ϕ2at multiples of 2π/q eventually coincide, as can be appreciated in the right column of Fig. 6 As a remnant of finite propagation time T, we witness in all panels of Figure 6 an inclination of the resonance peak, which to a smaller extent is noticeable also in Figure 5. It stems from the global background of hσziTvisible in Figure 4. Owing to the scaling of the abscissa, it diminishes with increasing Tand, in accordance with equation (24), it eventually vanishes. 7 Conclusions We have studied the asymptotic limit of multichromatically driven, dissipative dynamical systems. It turned out 8 Mar´ıa Laura Olivera-Atencio et al.: Generic shape of multichromatic resonance peaks −5π0 5π −0.43 −0.42 q δω T hσziT T= 2000/Ω1 −5π0 5π q δω T T= 2000/Ω1 −0.43 −0.42 hσziT T= 500/Ω1T= 500/Ω1 −0.43 −0.42 hσziT ϕ2= 0 ϕ2= 2π9 10q T= 250/Ω1 ϕ2= 0 ϕ2= 2πq−1 q T= 250/Ω1 Fig. 6. Resonance peak (q, p) = (7,5) for the averaging times displayed in the graphics, while all other parameters are as in Figure 4. Left column: Result for 10 equally spaced phases ϕ2= 2πn/(10q) for n= 0,1,...,9 showing that for sufficiently large T, the curve for ϕ2= 0 (n= 0, blue) smoothly connects to the one for ϕ= 2π/q. Right column: The same but for ϕ= 2πn/q,n= 0,1,...,q−1. With increasing T, the curves start to coincide. that a strict distinction between commensurable and incommensurable frequencies requires an infinite propagation time. Nevertheless, there exits a noticeable difference between the two cases, namely that only for commensurable frequencies the phases of the driving fields may matter. Moreover, the phase dependence generally has a lower periodicity than the naively expected 2π. While this implies non-generic features of the resonances, the resulting peaks upon variation of the phase exhibit a generic form given by sinc functions. While the limiting behavior can be derived analytically, we have performed numerical studies to see how the limits are approached. For a classical random walk on a lattice with bichromatically time-dependent transition rates, the velocity has been obtained analytically up to the numerical computation of a sum. The response as a function of the two driving frequencies shows how resonances emerge around rational values of Ω2/Ω1. The case of a dissipative two-level system has been treated fully numerically. It revealed how with increasing propagation time, the generic features of resonance peaks emerge, namely the sinc shape and the sub 2πperiodicity in the phase shift. The width of the resonance peaks at rational frequency quotients shrinks with increasing averaging time, such that the background eventually appears flat and the peaks become pronounced. The value of the response depends on the relative phase of the two drivings, while in its vicinity, the response becomes phase independent. An important point in practical calculations is that commensurable frequency ratios with rather large numerator or denominator imply large periods. Then owing to the necessarily finite propagation time, this periodicity may still not be manifest in the result. In other words, up to such finite time, the system behaves as if it were quasiperiodic, i.e., as if the frequencies were incommensurable. However, in particular for the random-walk model, it may take even considerably longer until the peaks assume their generic shape. Quantitative statements about this issue still represent a challenge for future investigations. This work was supported by the Spanish Ministry of Science, Innovation, and Universities through the CSIC Research Platform on Quantum Technologies PTI-001 and via grants No. MAT2017-86717-P and FIS2017-86478-P. Author contribution statement JCP has derived the analytical results, MLOA and SK have performed the numerical calculations for the random walk model and the two-level system, respectively. All authors have contributed to writing the manuscript. References 1. L. Gammaitoni, P. H¨anggi, P. Jung, F. Marchesoni, Rev. Mod. Phys. 70, 223 (1998) 2. J. Casado-Pascual, J. G´omez-Ord´o˜nez, M. Morillo, P. H¨anggi, Phys. Rev. Lett. 91, 210601 (2003) 3. V. Anishchenko, A. Neiman, A. Astakhov, T. Vadiavasova, L. Schimansky-Geier, Chaotic and Stochastic Processes in Dynamic Systems (Springer, Berlin, 2002) 4. J.A. Freund, L. Schimansky-Geier, P. H¨anggi, Chaos 13, 225 (2003) 5. B. Lindner, J. Garcia-Ojalvo, A. Neiman, L. SchimanskyGeier, Phys. Rep. 392, 321 (2004) 6. J. Casado-Pascual, J. G´omez-Ord´o˜nez, M. Morillo, J. Lehmann, I. Goychuk, P. H¨anggi, Phys. Rev. E 71, 011101 (2005) 7. I. Goychuk, J. Casado-Pascual, M. Morillo, J. Lehmann, P. H¨anggi, Phys. Rev. Lett. 97, 210601 (2006) 8. P. Reimann, Phys. Rep. 361, 57 (2002) 9. P. H¨anggi, F. Marchesoni, Rev. Mod. Phys. 81, 387 (2009) 10. D. Cubero, F. Renzoni, Brownian Ratchets: From Statistical Physics to Bio and Nano-motors (Cambridge University Press, Cambridge, 2016) 11. S. Flach, O. Yevtushenko, Y. Zolotaryuk, Phys. Rev. Lett. 84, 2358 (2000) 12. D. Cubero, F. Renzoni, Phys. Rev. E 97, 062139 (2018) 13. D. Cubero, G.R. Robb, F. Renzoni, Phys. Rev. Lett. 121, 213904 (2018) 14. J. Casado-Pascual, D. Cubero, F. Renzoni, Phys. Rev. E 88, 062919 (2013)