Universal patterns in multifrequency-driven dissipative systems

The dynamics of strongly driven systems may be rather complex, in particular when the driving consists of various components with different frequencies. A particular and well-studied case is a higher harmonic added to the driving with a fundamental frequency. This type of bichromatic driving can be used to control spatio-temporal symmetries via the relative phase between the two monochromatic fields [1–5]. It allows one to induce directed motion by the action of an oscillating field with zero mean, giving rise to the celebrated ratchet effect [6,7]. Moreover, one may use bichromatic driving for quantum state preparation [8,9].

By contrast, there exists considerably less work on driving forces with two incommensurable frequencies, i.e., frequencies whose proportion is an irrational number. An intriguing feature of such drivings is that while each of its components is time periodic, the system as a whole lacks discrete time-translation symmetry. In spite of this, the response may have higher symmetry than in the commensurable case [5,10]. The reason for this is that, in the incommensurable case, the long-time average is equivalent to the average over the relative phases among the driving components [11,12]. In ref. [10] the difference between commensurable and incommensurable drivings has been demonstrated both theoretically and experimentally for the electron transport through a double quantum dot.

At first sight, the distinction between drivings with commensurable and incommensurable frequencies seems surprising, since irrational numbers can be approximated to any degree of accuracy by rational numbers. This apparent paradox has motivated several studies of the response of multi-chromatically driven systems —both dissipative [11–14] and Hamiltonian ones [15,16]—as a function of one driving frequency while keeping the others constant.

In this perspective, we shed light on how the discrete time translation symmetries of each periodic driving component provides resonance peaks with a generic shape. These results are illustrated with some paradigmatic examples of classical and quantum dissipative systems. Moreover, we demonstrate that the phase average of the long-time response can be computed with a matrix continued-fraction method originally developed for incommensurable frequencies [10]. Possible implementations in quantum optical systems are also suggested together with an outlook for further studies.

We are interested in systems whose dynamical equations depend on time through N time-periodic functions of the form

$$\begin{equation} f_j(t)=\epsilon_j \cos\left(\Omega_j t+\varphi_j\right), \end{equation} \tag{ 1 }$$

where $j=1,\dots,N$, and $\epsilon_j$, $\Omega_j$, and $\varphi_j$ denote, respectively, the amplitude, the angular frequency, and the initial phase of f_j (t). The state of the system at time t will be denoted as $\mathscr{S}(t)$. Depending on the case, $\mathscr{S}$ may represent the values of a finite number of state variables characterizing a classical deterministic system, the density operator of a quantum system, the one-time probability density of a classical stochastic system, etc.

We focus on generic properties that do not depend on the precise nature of the functions f_j (t) and the specific details of the underlying dynamics. Our only assumption is that there exists a unique steady state, $\smash{\mathscr{S}^{\mathrm{st}}(t)}$, to which the system converges in the long-time limit —an assumption that holds for a wide class of dissipative systems. In general, the steady state will depend on the specific values taken by the parameters appearing in the functions f_j (t). When necessary, this dependence will be made explicit by the notation $\smash{\mathscr{S}^{\mathrm{st}}(t, \bm{\epsilon}, \bm{\Omega}, \bm{\varphi})}$, where $\bm{\epsilon}$, $\bm{\Omega}$, and $\bm{\varphi}$ are N-dimensional vectors with components $\epsilon_j$, $\Omega_j$, and $\varphi_j$, respectively.

Since we are assuming a unique steady state, its time evolution must be uniquely determined by the dynamical equations. Hence, the steady state shares the symmetry properties of the dynamical equations. To be specific, the set of functions in eq. (1) is invariant under the N + 1 transformations

$$\begin{eqnarray} \mathcal{T}^{(j)}&{\!\!:} \{t,\bm{\epsilon},\bm{\Omega},\bm{\varphi}\}\mapsto\{t,\bm{\epsilon}^{(j)},\bm{\Omega},\bm{\varphi}+\pi \bm{u}^{(j)}\}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \mathcal{T}&{\!\!:} \{t,\bm{\epsilon},\bm{\Omega},\bm{\varphi}\}\mapsto \{t+\tau,\bm{\epsilon},\bm{\Omega},\bm{\varphi}-\tau \bm{\Omega}\}, \end{eqnarray} \tag{ 3 }$$

where $\smash{\bm{\epsilon}^{(j)}}$ is the vector of amplitudes with the sign of component j inverted, while $\smash{\bm{u}^{(j)}}$ is the j-th canonical basis vector. More formally, $\smash{\epsilon^{(j)}_k=(1-2 \delta_{j,k})\epsilon_k}$ and $\smash{u_k^{(j)}=\delta_{j,k}}$, with $\delta_{j,k}$ being the Kronecker delta. Since the only explicit dependence of the dynamical equations on t, $\bm{\epsilon}$, $\bm{\Omega}$, and $\bm{\varphi}$ comes from the functions f_j (t), the steady state will also be invariant under these N + 1 transformations, i.e.,

$$\begin{align} \mathscr{S}^{\mathrm{st}}(t,\bm{\epsilon},\bm{\Omega},\bm{\varphi})&=\mathscr{S}^{\mathrm{st}}(t,\bm{\epsilon}^{(j)},\bm{\Omega},\bm{\varphi}+\pi \bm{u}^{(j)}) \end{align} \tag{ 4 }$$

$$\begin{align} &=\mathscr{S}^{\mathrm{st}}(t+\tau,\bm{\epsilon},\bm{\Omega},\bm{\varphi}-\tau \bm{\Omega}). \end{align} \tag{ 5 }$$

Let $Q=Q(\mathscr{S})$ represent a certain (physical) quantity that depends on the state of the system. In particular, in the steady state, the dependence of Q on t, $\bm{\epsilon}$, $\bm{\Omega}$, and $\bm{\varphi}$ is $\smash{Q^{\mathrm{st}}(t,\bm{\epsilon},\bm{\Omega},\bm{\varphi})\equiv Q[\mathscr{S}^{\mathrm{st}}(t,\bm{\epsilon},\bm{\Omega},\bm{\varphi})]}$. By applying eq. (4) twice, it follows that $\smash{ Q^{\mathrm{st}}(t,\bm{\epsilon},\bm{\Omega},\bm{\varphi})}$ is $2\pi$-periodic in all the components of the vector $\bm{\varphi}$, i.e., $\smash{Q^{\mathrm{st}}(t,\bm{\epsilon},\bm{\Omega},\bm{\varphi}+2\pi \bm{u}^{(j)})=Q^{\mathrm{st}}(t,\bm{\epsilon},\bm{\Omega},\bm{\varphi})}$ for $j=1,\dots,N$. In addition, taking in eq. (5) $\tau=-t$, one obtains that $\smash{Q^{\mathrm{st}}(t,\bm{\epsilon},\bm{\Omega},\bm{\varphi})= Q^{\mathrm{st}}(0,\bm{\epsilon},\bm{\Omega},\bm{\varphi}+ \bm{\Omega}t)}$, i.e., the time evolution of $\smash{Q^{\mathrm{st}}(t,\bm{\epsilon},\bm{\Omega},\bm{\varphi})}$ admits a description in terms of a time-dependent phase vector of the form $\bm{\varphi}+ \bm{\Omega} t$. Taking into account these two properties and performing a Fourier expansion in $\bm{\varphi}$, it is easy to see that the time average of $Q^\text{st}$ from 0 to T reads

$$\begin{equation} \overline{Q}_{T}(\bm{\epsilon},\bm{\Omega},\bm{\varphi})=\sum_{\bm{k}\in \mathbb{Z}^N}q_{\bm{k}}(\bm{\epsilon},\bm{\Omega})e^{i\bm{k}\cdot( \bm{\varphi}+\bm{\Omega}T/2)}\mathrm{sinc}\left(\frac{\bm{k}\cdot\bm{\Omega}T}{2}\right), \end{equation} \tag{ 6 }$$

where $\mathrm{sinc}(x)\equiv\sin(x)/x$ denotes the unnormalized sinus cardinalis, while the centered dot denotes the usual scalar product of N-dimensional vectors, and

$$\begin{equation} q_{\bm{k}}(\bm{\epsilon},\bm{\Omega})=\int_{-\pi}^{\pi}\dots\int_{-\pi}^{\pi}e^{-i\bm{k}\cdot \bm{\varphi}}Q^{\mathrm{st}}(0,\bm{\epsilon},\bm{\Omega},\bm{\varphi})\prod_{j=1}^N\frac{{\rm d}\varphi_j}{2\pi}. \end{equation} \tag{ 7 }$$

Note that, according to eqs. (4) and (5), the Fourier coefficients $q_{\bm{k}}$ satisfy the symmetry property $q_{\bm{k}}(\bm{\epsilon}^{(j)},\bm{\Omega})=(-1)^{k_j}q_{\bm{k}}(\bm{\epsilon},\bm{\Omega})$ and, hence, can be written in the form

$$\begin{equation} \smash{q_{\bm{k}}(\bm{\epsilon},\bm{\Omega})\equiv \gamma_{\bm{k}}(\bm{\epsilon},\bm{\Omega})\prod_{j=1}^N\epsilon_j^{|k_j|}}, \end{equation} \tag{ 8 }$$

where the functions $\smash{\gamma_{\bm{k}}(\bm{\epsilon},\bm{\Omega})}$ are even in each of the arguments $\epsilon_j$. Under quite general conditions, it can be shown that the functions $\smash{\gamma_{\bm{k}}(\bm{\epsilon},\bm{\Omega})}$ admit a Taylor expansion in the amplitudes $\epsilon_j$ [11]. In practice, for sufficiently small values of the amplitudes $\epsilon_j$ expressed in suitable dimensionless units, this expansion can be truncated after a few terms [11]. The dependence of the Fourier coefficients on the amplitudes in eq. (8) can also be obtained by a functional expansion on the driving [17,18].

In the limit $T\to \infty$, the sinc functions appearing in eq. (6) vanish unless the resonance condition $\bm k\cdot\bm\Omega = 0$ is fulfilled. Therefore, eq. (6) leads to

$$\begin{equation} \overline{Q}_{\infty}(\bm{\epsilon},\bm{\Omega},\bm{\varphi})=\sum_{\bm{k}\in S_{\bm\Omega}^{\perp}}q_{\bm{k}}(\bm{\epsilon},\bm{\Omega})e^{i\bm{k}\cdot \bm{\varphi}}, \end{equation} \tag{ 9 }$$

where $S_{\bm\Omega}^{\perp}$ is the set of vectors $\bm k$ that have integer components and are orthogonal to $\bm \Omega$. In practice, the limit T → ∞ can be calculated only approximately by taking a sufficiently large value of T. If we consider the vicinity of a resonance at a fixed frequency vector $\bm{\Omega_0}$ and focus on driving frequencies $\bm \Omega = \bm\Omega_0 +\bm{\delta\omega}$, with $\bm{\delta\bm\omega}$ of the same order of magnitude as $T^{-1}$, then the asymptotic behavior of eq. (6) for $T\to \infty$ is given by [12]

$$\begin{eqnarray}\overline{Q}_{T}(\bm{\epsilon},\bm{\Omega_0}+\bm{\delta \omega},\bm{\varphi}) &\sim\sum_{\bm{k}\in S_{\bm\Omega_0}^{\perp}}q_{\bm{k}}(\bm{\epsilon},\bm{\Omega_0}) e^{i\bm{k}\cdot( \bm{\varphi}+\bm{\delta \omega}T/2)}\nonumber\\ &\quad \times\mathrm{sinc}\left(\frac{\bm{k}\cdot\bm{\delta \omega}T}{2}\right). \end{eqnarray} \tag{ 10 }$$

Rather importantly, owing to the orthogonality condition $\bm k\cdot \bm\Omega_0 = 0$, the only dependence on $\bm\Omega_0$ is contained in the Fourier coefficients $q_{\bm{k}}$.

A non-trivial solution of the equation $\bm k\cdot\bm\Omega = 0$ requires that the N components of $\bm \Omega$ are commensurable, i.e., that one of the frequencies can be expressed as a linear combination of the others with rational coefficients. Otherwise, the set $S_{\bm\Omega}^{\perp}$ reduces to the trivial solution $\bm{k} =\bm{0}$, and $q_{\bm{0}}(\bm{\epsilon},\bm{\Omega})$ is the only non-vanishing term in $\overline{Q}_{\infty}(\bm{\epsilon},\bm{\Omega},\bm{\varphi})$. Since $q_{\bm{0}}(\bm{\epsilon},\bm{\Omega})$ is independent of the phases in $\bm{\varphi}$, it provides a smooth background for sinc-shaped peaks in $\overline{Q}_{T}(\bm{\epsilon},\bm\Omega,\bm\varphi)$. Interestingly, this background vanishes if the dynamical equations are invariant under a mapping that involves the phases $\bm{\varphi}$ and inverts the sign of the observable Q. Then, $\overline{Q}_{T}(\bm{\epsilon},\bm{\Omega},\bm{\varphi}) = -\overline{Q}_{T}(\bm{\epsilon},\bm{\Omega},\tilde{\bm{\varphi}})$, which links the response for the phase $\bm{\varphi}$ and the transformed phase $\tilde{\bm\varphi}$. If, in addition, the Jacobian of the phase transformation is unity (which is fulfilled for any phase inversion and phase shift), one can conclude from eq. (7) that $q_{\bm 0}(\bm\epsilon,\bm\Omega) = -q_{\bm 0}(\bm\epsilon,\bm\Omega) = 0$.

Equation (9) contains a phase factor $e^{i\bm{k}\cdot\bm{\varphi}}$ which is invariant under a phase shift $\bm{\varphi} \to \bm{\varphi}+\delta\bm{\varphi}$ provided that

$$\begin{equation} \bm{k}\cdot \delta\bm{\varphi} = 2\pi \end{equation} \tag{ 11 }$$

(or any other multiple of $2\pi$), for all the vectors $\bm k$ orthogonal to $\bm \Omega$ and with integer components. In the case that the N components of $\bm \Omega$ are pairwise commensurable (i.e., if there exist a frequency $\Omega_0$ and an N-dimensional vector $\bm{n}$ with positive integer components such that $\bm{\Omega} =\Omega_0 \bm{n}$), the orthogonality condition $\bm k\cdot \bm\Omega = 0$ becomes equivalent to the Diophantine equation $\bm k\cdot \bm n = 0$. This Diophantine equation together with eq. (11) implies invariance of the response under non-trivial phase shifts $\delta\varphi_j<2\pi$.

While for multi-chromatic driving, the general solution of the Diophantine equation $\bm k\cdot \bm n = 0$ may be complicated (see, e.g., refs. [11,12,19,20]), for bichromatic driving, it can be derived explicitly. Setting $\bm{\Omega} = (q,p)\Omega_0$ with p and q coprime, the general solution is $(p,-q) \ell$, with ℓ being any integer. Therefore, only terms with $\bm{k}$ being integer multiples of $(p,-q)$ contribute to $\overline{Q}_\infty$ as discussed above. Then, condition (11) becomes $p\delta\varphi_1 - q\delta\varphi_2 = 2\pi$. Since only relative phases of the drivings f_j (t) matter, we can set one phase to zero, such that we can conclude invariance of the response $\overline{Q}_\infty$ for the phase shifts [5,12,13]

$$\begin{align} \varphi_1 \to{}& \varphi_1 + 2\pi/p, \end{align} \tag{ 12 }$$

$$\begin{align} \varphi_2 \to{}& \varphi_2 + 2\pi/q. \end{align} \tag{ 13 }$$

Notice that this is a kind of cross-relation, since the integer q or p that defines the frequency of one of the drivings appears in the phase invariance of the other driving. In ref. [12], it has been demonstrated numerically that the equivalence of these phases holds (approximately) in a whole vicinity of a (q, p)-resonance.

To illustrate the features discussed so far, we provide explicit numerical results for a classical random walk and a dissipative quantum mechanical two-level system.

Classical system

Several classical models have been considered in the literature to analyze generic properties of dissipative dynamical systems under multifrequency drivings. For example, a one-dimensional model consisting of a Brownian particle, moving in a periodic potential, under the influence of a biharmonic force has been used to study some general asymptotic properties of driven nonlinear dissipative systems in the long-time limit [13]. This same model has also been considered to elucidate the connection between irrationality and quasiperiodicity in these kinds of systems [14]. The generality of these results has been revealed by replacing the periodic potential with a double-well potential (see the supplemental material in ref. [14]). In ref. [12], the generic shape of the resonance peaks in the vicinity of commensurable frequencies has been illustrated using a classical random walk model. Given the simplicity of this latter model, it will be the one considered here.

As an example for a classical stochastic process, we thus employ an infinite one-dimensional chain with thermal nearest-neighbor hopping with the forward and backward rates

$$\begin{equation} r_\pm(t) = r_0 e^{-\beta [E_0 \pm\Delta E(t)]}, \end{equation} \tag{ 14 }$$

where β denotes the inverse thermal energy $1/k_BT$ [12]. The energy difference between two adjacent sites with distance a contains a static contribution E₀ and a time-dependent one, $\Delta E(t) = f(t)/\beta$ with

$$\begin{equation} f(t) = A_1\cos(\Omega_1 t+\varphi_1) + A_1\cos(\Omega_2 t+\varphi_2). \end{equation} \tag{ 15 }$$

It can be shown that the stationary state of the corresponding master equation reads

$$\begin{equation} v(t) = a[r_+(t) - r_-(t)] = v_0 \sinh[f(t)], \end{equation} \tag{ 16 }$$

with $v_0 = 2ar_0 \exp(-\beta E_0)$. For details of the calculation, see ref. [12].

To evaluate the long-time average of the velocity, it is convenient to decompose v(t) into a Taylor series in the amplitudes A_i . Then the time integration of each term can be evaluated analytically, while for the summation of the resulting terms, we resort to numerics [12]. The result for $\Omega_2$ in the vicinity of $\Omega_1/4$, i.e., close to the (4, 1) resonance, is depicted in fig. 1. It nicely shows that while the response for equivalent phases defined in the previous section may be different for small averaging times, all the curves become indistinguishable for sufficiently large T. Moreover, once convergence is practically reached, the curves exhibit the sinc shape proposed for their enveloping function.

Zoom In Zoom Out Reset image size

Dissipative two-level system

In ref. [12], the shape of the resonance peaks has been investigated also for the two-level Hamiltonian as an example for dissipative quantum systems. Here, we consider a quantum mechanical system defined by the Hamiltonian

$$\begin{eqnarray} H(t) =\,& \frac{\Delta}{2}(\sigma_z\cos\theta + \sigma_x\sin\theta)\nonumber\\ & {}+ A_1\sigma_x\cos(\Omega_1 t) + A_2\sigma_z\cos(\Omega_2 t+\varphi), \end{eqnarray} \tag{ 17 }$$

with the Pauli matrices $\sigma_{x,z}$ and the angle θ which allows the control of the symmetry, as we will see below. Dissipation is provided by a Lindblad form such that the quantum master equation for the density operator reads

$$\begin{equation} \dot\rho = -\frac{i}{\hbar}[H(t),\rho] +\Gamma(2 \sigma_\downarrow\rho\sigma_\uparrow - \sigma_\uparrow\sigma_\downarrow\rho -\rho \sigma_\uparrow\sigma_\downarrow), \end{equation} \tag{ 18 }$$

with the dissipation rate Γ and the Lindblad operator $\sigma_\downarrow \equiv |\phi_0\rangle\langle\phi_1| = \sigma_\uparrow^\dagger$, which is the projector to the ground state $|\phi_0\rangle$ for given angle θ.

Figure 2(a) depicts the behavior of the long-time solution for which we predicted the result in eq. (6). The grey lines show the long-time average $\overline{\langle\sigma_z\rangle}_T$ for various initial phases. They are computed via straightforward numerical propagation of the Lindblad master equation (18). The curves possess resonance peaks at rational values of $\Omega_1/\Omega_2$. Their enveloping function clearly exhibits the shape of a sinc function, which implies that the series in eq. (6) is governed by a single coefficient with index ${\bm k}\neq{\bm 0}$. The smooth background corresponds to the only coefficient for which the sinc becomes equal to unity, namely $q_{\bm 0}({\bm\Omega})$. Figure 2(b) visualizes how the shape of the enveloping function of the response for different phases emerges. For sufficiently large T, firstly the sinc shape is assumed. Then with T increasing further, the sinc becomes ever narrower and eventually shrinks to a single discontinuity located at $\Omega_2 = (p/q)\Omega_1$.

Zoom In Zoom Out Reset image size

As a consequence of spatio-temporal symmetries, the background may vanish. For example, when $\theta=\pi/2$ the Hamiltonian H(t) is invariant under unitary transformation with $\sigma_x$ accompanied by a phase shift $\varphi\to\varphi+\pi$, while our observable $\sigma_z$ acquires a minus sign. Moreover, since the Lindblad dissipator is defined via the eigenstates of the time-independent part of the Hamiltonian, it is invariant under this transformation as well. Since for incommensurable frequencies, the phase is irrelevant in the limit T → ∞, the time-averaged response is equal to its negative value and, hence, must be zero.

An established technique for treating periodically driven systems is the Floquet theory. It is based on the discrete time translation by the period of the driving T. Under this symmetry, linear differential equations $\dot\psi = L(t)\psi$ possess a complete set of solutions of the form $\psi(t) = e^{-i\mu t}\phi(t)$, where $\phi(t) = \phi(t+T)$ shares the time periodicity of the driving [21]. Then the Floquet function $\phi(t)$ is an eigen solution of $L(t)-\partial_t$ in a Hilbert space extended by a periodic time coordinate [22,23].

For bichromatic driving with incommensurable frequencies, the discrete time translation symmetry gets lost. Nevertheless, one can employ a Floquet ansatz extended by a further Fourier index that reflects the periodicity of the second driving. Hence the solutions are still of the form $\psi(t) = e^{-i\mu t}\phi(t)$, but now with the modified Floquet function [24,25]

$$\begin{equation} \phi(t) = \sum_{\bm k} e^{-i{\bm k}\cdot{\bm\Omega}t} \phi_{\bm k}({\bm\Omega}), \end{equation} \tag{ 19 }$$

where ${\bm k} = (k_1,k_2)$. It corresponds to a two-dimensional Fourier ansatz which for very strong driving may be numerically expensive. Notice that in contrast to the ansatz for the long-time solution (6), the general solution of the dynamical system, $\psi(t)$, contains an exponential prefactor.

For dissipative equations of motion, an efficient numerical method has been developed in ref. [10]. It starts from the observation that the Floquet index μ of the long-time solution must vanish. Then one readily obtains a set of coupled homogeneous equations for the Fourier coefficients $\phi_{\bm k}({\bm\Omega})$. The idea is now to derive for one Fourier index, say k₂, a recurrence equation which can be solved with the matrix continued fraction (MCF) method with a numerical effort that grows only linearly with the size of the cutoff index. Finally, this provides the Fourier coefficient $\phi_{\bm 0}(\bm\Omega)$ which contains all information about the time-averaged long-time solution. Hence, it is equivalent to the component $q_{\bm 0}(\bm\Omega)$ of the long-time solution (6).

While the ansatz (19) looks rather natural, it has to be handled with care, because for commensurable frequencies, it is overcomplete. Technically this may lead to divergences in the MCF iteration. In practice, dissipation generally cures this problem, but its emergence cannot be ruled out.

Hence, for commensurable frequencies, the MCF algorithm provides the phase average of the long-time average (6). In fig. 2, we verify this numerically for the case of the two-level system defined in eq. (17). The red curve is computed with the MCF iteration and indeed provides the smooth background of the resonance peaks. A further test may be performed with the mixing angle $\theta=\pi/2$ for which the symmetry considerations above predict $\overline{\langle\sigma_z\rangle_T} = 0$.

The dissipative dynamics of eq. (18) can be carried out straightforwardly with a trapped ion quantum platform, as has been demonstrated both theoretically [26] and experimentally [27]. We consider a two-ion system, where the first ion will encode the two-level system under study, and the second ion will be an ancillary qubit that will provide the dissipative part. The unitary part given by H(t) amounts to a single-qubit time-dependent operation, which can always be decomposed onto single-qubit drivings with appropriate laser intensities and frequencies [28]. With respect to the dissipative, Lindblad-form term of eq. (18), one can couple the previous two-level system of the first ion with a second two-level system of the ancillary ion, and perform a digital decomposition of the Lindblad dynamics, as described in ref. [26]. In each digital step, one would implement the Kraus operators, of the form (considering only the dissipative part for simplicity, while the unitary part would be carried out with a subsequent digital step)

$$\begin{equation} \rho(t)=E_0\rho(0)E_0^\dag + E_1\rho(0)E_1^\dag, \end{equation} \tag{ 20 }$$

where

$$\begin{equation} E_0= \begin{pmatrix} \sqrt{1-\gamma'} & 0 \\ 0 & 1 \end{pmatrix},\quad\,\,\, E_1= \begin{pmatrix}0 & 0 \\ \sqrt{\gamma'} & 0 \end{pmatrix} \end{equation} \tag{ 21 }$$

are 2× 2 matrices acting on the two-level subspace of the considered system with exponentially decaying $1-\gamma'=\exp(-2 \Gamma t)$. These Kraus operators correspond to the dissipative channel in the basis $|\phi_0\rangle$, $|\phi_1\rangle$, providing, in each small time step, the Lindblad part of the dynamics of eq. (18) [29]. To carry out these operations in a digital quantum simulator with the two-ion system, for each digital step, one would initialize the ancilla qubit in state $|0\rangle$ and apply a two-qubit gate U₂ such that $\langle 0|U_2|0\rangle=E_0$ and $\langle 1|U_2|0\rangle=E_1$ are the required matrix elements in the ancillary qubit basis (corresponding to single-qubit gates in the system qubit). A U₂ fulfilling these requirements can always be obtained via at most three CNOT gates combined with single-qubit gates [29]. Subsequently, as described in refs. [26,27], one would apply optical pumping to the ancillary qubit, to map it to state $|0\rangle$, providing the entropy increase that realizes the dissipation. Finally, one would carry out the unitary part of H(t). The complete master equation dynamics would be provided by the subsequent iteration of this digital step for n total steps. The long-term solution would be obtained for sufficiently large n, and the measurement $\overline{\langle\sigma_z\rangle}_T$ can be straightforwardly carried out with the trapped ion system via resonance fluorescence [28].

In this article we have reviewed the generic behavior of multi-chromatically driven dissipative systems in the classical as well as in the quantum mechanical dynamics. Most prominently, in the time-averaged signal as a function of one driving frequency, one observes phase-dependent peaks with a sinc-shaped envelope on top of a smooth phase-independent background. The width of the peaks diminishes with the averaging time. By contrast, the background does not depend on the phases and converges rather rapidly to its asymptotic value. Symmetries may suppress the background, while features of the peaks remain. It is worth mentioning that the width of these peaks is generally smaller than the one predicted by the Fourier inequality [30], where the difference can be expressed as a factor determined by the driving frequencies [13]. Therefore, our results may have a direct and practical application for the identification of dissipative systems displaying sub-Fourier resonances [13,14,30].

To observe these peaks, one has to leave the linear response limit and enter the regime of harmonic mixing. Then with an increasing amplitude an increasing number of resonances emerges, but "simple resonances" such as 1/1, 1/2, or 2/3 dominate the overall behavior. A challenging open problem is the question whether there is any rule for the relative magnitude of the peaks as a function of the "simplicity" of the frequency ratio.

For the computation of the numerical examples, we have used simple propagation schemes which, however, may be rather time consuming. For the phase-averaged response, it turned out that a two-frequency Floquet theory provides reliable results for practically all frequency ratios, despite that its convergence is guaranteed only for incommensurable frequencies. This allows one to employ a computational method based on matrix continued fractions, which is numerically rather efficient.

Finally, future explorations of the magnitudes of the resonances for further systems may deepen our understanding of multi-chromatically driven systems and may open perspectives in science, technology, and industry.

We acknowledge funding by the Junta de Andalucía (P20-00617 and US-1380840) and by the Spanish Ministry of Science, Innovation, and Universities under grant Nos. FIS2017-86478-P, PGC2018-095113-B-I00, PID2019-104002GB-C21, PID2019-104002GB-C22, PID2020-117787GB-I00 (MCIU/AEI/FEDER, UE), and via the CSIC Research Platform on Quantum Technologies PTI-001.

Data availability statement: The data that support the findings of this study are available upon reasonable request from this author.