Floquet oscillations in periodically driven Dirac systems

TL;DR

This paper demonstrates that periodically driven Dirac systems can exhibit Floquet oscillations, analogous to Bloch oscillations, with tunable periods and richer patterns, observable in materials like graphene under circularly polarized light.

Contribution

It introduces the concept of Floquet oscillations in Dirac systems and shows how they can be generated and analyzed using different driving protocols.

Findings

Floquet oscillations occur in driven Dirac systems with linear dispersion.

Oscillation patterns are richer than traditional Bloch oscillations.

Oscillation period can be tuned via driving frequency.

Abstract

Electrons in a lattice exhibit time-periodic motion, known as Bloch oscillation, when subject to an additional static electric field. Here we show that a corresponding dynamics can occur upon replacing the spatially periodic potential by a time-periodic driving: Floquet oscillations of charge carriers in a spatially homogeneous system. The time lattice of the driving gives rise to Floquet bands that take on the role of the usual Bloch bands. For two different drivings (harmonic driving and periodic kicking through pulses) of systems with linear dispersion we demonstrate the existence of such oscillations, both by directly propagating wave packets and based on a complementary Floquet analysis. The Floquet oscillations feature richer oscillation patterns than their Bloch counterpart and enable the imaging of Floquet bands. Moreover, their period can be directly tuned through the driving frequency. Such oscillations should be experimentally observable in effective Dirac systems, such as graphene, when illuminated with circularly polarized light.

Equations

$$H(\mathbf{k},t)=H_0(\mathbf{k})+V_T(t)=H(\mathbf{k},t+T).$$

$$H_0(\mathbf{k})=\hslash v_F\mathbf{k}\cdot \mathbf{\sigma },$$

$$V(x)=-eEx,e>0,$$

$$k(t)=k_i-(eE/\hslash )t.$$

$$T_B=2\pi /(a|eE|\hslash ).$$

$$x_{\alpha \alpha }(k(t))=\frac{\hslash }{eE}\left[{ϵ}_{\alpha }(k(t))-{ϵ}_{\alpha }(k_i)\right].$$

$$T_F=\frac{\hslash \omega }{v_F|eE|}.$$

$$V_T(t)=\sum _{l=1}^{\infty }\Theta (t-(lT-\Delta t))\Theta (lT-t)M{\sigma }_z$$

$$\tilde{\Psi }(k,0)=\sqrt{\frac{1}{\sqrt{\pi }\Delta k}}\exp \left(-\frac{1}{2\Delta k^2}(k-{\bar{k}}_i{)}^2\right)\cdot \left(\begin{array}{c}1\\ 1\end{array}\right),$$

$$v(t)=⟨\hat{v}(t)⟩=\frac{i}{\hslash }⟨\Psi (t)|\left[H(t),\hat{x}\right]|\Psi (t)⟩,$$

$$U_k(t,0)=\sum _{\alpha =\pm }\exp \left(-\frac{i}{\hslash }{ϵ}_{\alpha }(k)t\right)|{\varphi }_{\alpha ,k}(t)⟩⟨{\varphi }_{\alpha ,k}(0)|.$$

$$|\Psi (0)⟩=\sum _{k_i}{c}_{\alpha ,k_i}|{\varphi }_{\alpha ,k_i}(0)⟩,$$

$$v(t)=$$

$$\times \exp \left(-\frac{i}{\hslash }\left[{ϵ}_{\beta }\left(k(t)\right)-{ϵ}_{\alpha }\left(k(t)\right]\right)t\right).$$

$$x(t)=⟨\hat{x}(t)⟩=\int |\Psi (x,t){|}^{2}x\text{d}x$$

$$v(t)=⟨\hat{v}(t)⟩=\frac{\text{d}}{\text{d}t}⟨\hat{x}(t)⟩.$$

$$V_T(t)=\mathbf{A}(t)\cdot \mathbf{\sigma }=A\left(\begin{array}{c}\cos (\omega t)\\ \sin (\omega t)\end{array}\right)\cdot \mathbf{\sigma }.$$

$$H(\mathbf{k})=\hslash v_F\mathbf{k}\cdot \mathbf{\sigma }-\mu \left[\left(k_y^2-k_x^2\right){\sigma }_x+2k_x{k}_y{\sigma }_y\right],$$

$${\mathcal{H}}_F=H_0(\mathbf{k})+V_T(t)-i\hslash {\partial }_t,$$

$${\Phi }_{\alpha }(\mathbf{k},t)=\sum _{n=-\infty }^{\infty }{\mathbf{u}}_{\alpha }(\mathbf{k},n){\text{e}}^{in\omega t}.$$

$$\begin{array}{rl}\sum _{n=-\infty }^{\infty }& \underset{{\mathcal{H}}_{F,mn}(\mathbf{k})}{\underbrace{\left({\mathcal{H}}_{0F,mn}(\mathbf{k})+V_{F,mn}\right)}}{\mathbf{u}}_{\alpha }(\mathbf{k},n)\\ & ={ϵ}_{\alpha }(\mathbf{k}){\mathbf{u}}_{\alpha }(\mathbf{k},m).\end{array}$$

$${\mathcal{H}}_{0F,mn}(\mathbf{k})=(\hslash v_F\mathbf{k}\cdot \mathbf{\sigma }+m\hslash \omega ){\delta }_{mn}$$

$$V_{F,mn}=\frac{1}{T}{\int }_0^T\text{d}t{V}_T(t){\text{e}}^{i(n-m)\omega t}$$

$$V_{F,mn}=\frac{M}{2\pi i(n-m)}\left(1-{\text{e}}^{-2\pi i(n-m)\frac{\Delta t}{T}}\right){\sigma }_z.$$

$$A(k)=-\frac{i}{\sqrt{1+{\eta }^2}}\sin \left(\mu \sqrt{1+{\eta }^2}\right),$$

$${\omega }_{ZB}=\left({ϵ}_{\beta }(k)-{ϵ}_{\alpha }(k)\right)/\hslash .$$

$$\text{i}\hslash \frac{\partial }{\partial t}|\psi ⟩=\hat{H}|\psi ⟩,$$

$$U(t,t_0)=\mathcal{T}\exp \left(-\frac{\text{i}}{\hslash }{\int }_{t_0}^t\hat{H}(t^{\prime })\text{d}{t}^{\prime }\right),$$

$$U(t,t_0)=U(t,t^{\prime })U(t^{\prime },t_0),$$

$$|\psi (t)⟩=U(t,t_0)|\psi (t_0)⟩.$$

$$U(t_0,t)=\exp \left(-\text{i}\frac{\hat{H}}{\hslash }\cdot (t-t_0)\right).$$

Full text

Floquet, Bloch, Floquet-Bloch oscillations, time periodic systems

I Introduction

In the early days of quantum mechanics, F. Bloch and C. Zener Bloch (1929); Zener (1934) predicted that electrons in a periodic potential, when accelerated by a constant external electric field, perform a time-periodic motion, by now well known as Bloch oscillation Wacker (2002). It took about 60 years until this phenomenon was observed in biased semiconductor superlattices Keldysh (1963); Feldmann et al. (1992); Waschke et al. (1993); Ignatov et al. (1994). Since then Bloch oscillations or analogs of them have been found in various systems ranging from cold atom gases Ben Dahan et al. (1996); Wilkinson et al. (1996) to classical optical Pertsch et al. (1999); Morandotti et al. (1999) and acoustic waves Sanchis-Alepuz et al. (2007), to name a few. In 2014, Bloch oscillations due to the crystal lattice of a biased bulk semiconductor were eventually observed O. Schubert et al. (2014).

In the meantime, scientific interest in tuning Bloch bands by means of external time-periodic driving has rapidly grown, especially since the proposal of so-called Floquet topological insulators N. H. Lindner et al. (2011) demonstrating the powerful influence external driving can exert on the properties of a crystal. Recent experiments also showed that Floquet band engineering allows for switching Bloch oscillations on and off Fujiwara et al. (2019). Moreover, additional driving can immensely increase the amplitude of conventional Bloch oscillations, giving rise to “super” Bloch oscillations Alberti A. et al. (2009); Haller et al. (2010); Arlinghaus and Holthaus (2011); Kudo and Monteiro (2011).

Here we propose to consider the opposite limit of a time-periodically driven system without any spatial lattice, but still subject to a constant external electric field. We demonstrate that, most notably, still spatially periodic motion of the charge carriers can appear. We call this type of dynamics Floquet oscillations since they arise from the periodic repetitions of Floquet quasi-energy bands. So far very few works have addressed Bloch-type oscillations in the absence of an external lattice. One interesting prediction refers to Bloch oscillations of light, *i.e. *frequency oscillations of photons Yuan and Fan (2016). Further Bloch-type oscillations were predicted theoretically for interacting 1d spinor gases Gangardt and Kamenev (2009) and recently observed in an atomic Bose liquid Meinert et al. (2017). In these settings interactions lead to the dynamical formation of periodic structures, which can yield oscillations à la Bloch. They do not, however, involve external drivings, and hence are of different nature than the Floquet oscillations predicted here.

Specifically, we show that such periodic modes can emerge in spatially uniform systems governed by effective Dirac Hamiltonians, where the linear dispersion converts the energy periodicity of the Floquet spectrum into approximately $k$-periodic bands. Most notably, Floquet oscillations are a quantum phenomenon distinctly different from the classical oscillatory motion of a charge in an ac field. Instead of following the driving frequency, they exhibit a frequency inversely proportional to it. Moreover, Floquet oscillations offer a possibility to directly image the Floquet quasi-band structure. Interestingly, they additionally show zitterbewegung features. We support our predictions by numerical calculations for two experimentally relevant prototypes of external driving, a periodic pulse sequence and circularly polarized radiation.

II General concept of Floquet oscillations

Consider a system $H_{0}(\mathbf{k})$ (with momentum operator $\mathbf{k}$) subject to a time-periodic driving $V_{T}(t)$ with period $T$ and frequency $\omega=2\pi/T$ described by the Hamiltonian

[TABLE]

Via Floquet theory Floquet (1883); Ince (1956); Magnus and Winkler (2004) the problem is transformed to finding the quasi-energy eigenvalues $\epsilon(\mathbf{k})$ of the Floquet Hamiltonian $\mathcal{H}_{F}(\mathbf{k})=H(\mathbf{k},t)-\textrm{i}\hbar\frac{\partial}{\partial t}$. The quasi-energies $\epsilon(\mathbf{k})$ extend to infinity in $\mathbf{k}$-space in absence of a spatial lattice. However, they are periodic in quasi-energy $\epsilon\in[-\hbar\omega/2,\hbar\omega/2]$ forming a sequence of Floquet replica which are the analogue of usual Bloch bands – the latter being periodic in quasi-momentum $k\in[-\pi/a,\pi/a]$ due to spatial periodicity with lattice constant $a$.

For the undriven system, we take an (effective) two-dimensional Dirac Hamiltonian

[TABLE]

with constant Fermi velocity $v_{F}$ and $\boldsymbol{\sigma}=(\sigma_{x},\sigma_{y})$ the vector of Pauli matrices. The spectrum of $H_{0}(\mathbf{k})$ is composed of two energy branches $E_{\pm}=\pm\hbar v_{F}k$ with $k=|\mathbf{k}|$. The Floquet bands $\epsilon_{\pm}$ that emerge when adding the time-periodic driving $V_{T}(t)$ to the system are sketched in Fig. 1. Due to the radial symmetry of the bands, we only show a cut along $k_{y}=0$. In the limit of infinitely small driving $V_{T}(t)$, the bare Dirac cone $E_{\pm}$ (black) is accompanied by a mesh of intersecting replica (gray) that are shifted in energy by multiples of $\hbar\omega$. The blue dashed lines mark the corresponding first Brillouin zone (BZ) in energy. For finite $V_{T}(t)$ band gaps open at replica crossings and separated Floquet bands emerge (red dashed curves, see Appendix A for an introduction to the relevant aspects of Floquet theory). The band gaps sketched in Fig. 1 emerge since different replica are coupled. Under the influence of the drift field $E$, a particle follows the band adiabatically along $k$ and successively absorbs/emits an increasing number $n$ of photons $\hbar\omega$ at each avoided crossing between two replica of the (undriven) Dirac cone with relative energy shift $n\hbar\omega$. Hence, the farther away from $k=0$, the more photons are needed to bridge the energy difference required for adiabatic dynamics associated with Floquet oscillations. As we will show below, time-periodic potentials $V_{T}$ with pronounced higher-order (in $\omega$) Fourier components will correspondingly open band gaps farther away from $k=0$ compared to, e.g., single frequency harmonic driving. Due to the underlying linear (Dirac) dispersion, the Floquet bands are approximately $k$-periodic, implying particularly pronounced Floquet oscillations with a well-defined frequency, in close analogy to Bloch oscillations 111General nonlinear dispersions $\epsilon_{\pm}(k)$, *e.g. *two hyperbolic branches, give rise to sequences of intersections that are not equidistant in $k$ and thereby deny a constant Floquet oscillation period..

To gain insight into the Floquet oscillation dynamics let us be definite and consider the time evolution of a wave packet (WP) of Dirac electrons under the influence of an additional, constant electric field $\mathbf{E}=E\mathbf{e}_{x}$. Here we consider the 1d motion along the field direction, generalizations to higher dimensions are straight forward. Due to the drift potential

[TABLE]

the WP is accelerated and its initial wave vector $k_{i}$ changes linearly in time Arlinghaus and Holthaus (2011):

[TABLE]

For ordinary Bloch bands, the BZ with period $\Delta k_{osc}\!=\!2\pi/a$ is traversed in the time

[TABLE]

While crossing the BZ, a charge carrier changes its velocity according to the change in slope of the $k$-periodic band structure, leading to a Bloch oscillation in $k(t)$ with frequency $\omega_{B}=2\pi/T_{B}$.

For Floquet systems the velocity operator is given by $\hat{v}=\textrm{d}\mathcal{H}_{F}/\textrm{d}k$ Arlinghaus and Holthaus (2011). The diagonal terms of a WP’s velocity read $\langle\hat{v}_{\alpha\alpha}\rangle\!=\!v_{\alpha\alpha}\!=\!\textrm{d}\epsilon_{\alpha}/\textrm{d}k$, with $\alpha\!=\!\pm$, and the corresponding position expectation value $\langle\hat{x}_{\alpha\alpha}\rangle=x_{\alpha\alpha}$ is obtained by time integration of $v_{\alpha\alpha}$. Using Eq. (4) this can be substituted by an integration over $k$ leading to

[TABLE]

These diagonal contributions to $\langle\hat{x}\rangle(t)$ encode features of the Floquet band structure into the WP position. In particular, analogously to conventional Bloch oscillations, the WP is expected to perform oscillations for Floquet bands similar to the ones sketched in Fig. 1 (red dashed lines).

During one (Floquet) oscillation $k$ changes by the period $\Delta k_{osc}$ of the Floquet bands (Fig. 1). Hence, Eq. (4) implies the corresponding period $T_{F}=(\hbar/|eE|)\Delta k_{osc}$. Due to the linear dispersion $E_{\pm}(k)=\pm\hbar v_{F}k$, we have $\hbar v_{F}\Delta k_{osc}=\hbar\omega$ so that the Floquet oscillation period reads

[TABLE]

$T_{F}$ is proportional to the inverse period $1/T=\omega/(2\pi)$ of the driving in Eq. (1). While its analogue, the Bloch period $T_{B}$, Eq. (5), is determined by the inverse (super-)lattice constant $1/a$, usually fixed in experiment, $T_{F}$ can be tuned through $\omega$ in a range such that $T_{F}>T$.

III Results for representative driving protocols

In the following we consider two representative types of driving well suited to generate Floquet oscillations: a periodic sequence of short pulses and a circularly polarized light field.

III.1 Short pulse sequence

For the first driving protocol we use

[TABLE]

in Eq. (1) (with Heaviside function $\Theta$ and Pauli matrix $\sigma_{z}$). This pulse train periodically couples the two branches of the Dirac Hamiltonian $H_{0}({\mathbf{k}})$, Eq. (2), by opening a mass gap of strength $M$ and duration $\Delta t$, see lower inset in Fig. 2(a). The resulting Floquet spectrum can be tuned either by $M$ or by the ratio $\Delta t/T$. To be definite we choose the normalized pulse strength $M/\hbar\omega=4.4/\pi$ and $\Delta t/T=0.09$. The resulting Floquet band structure for this representative set of parameters is shown in Fig. 2(a). The driving opens a sequence of gaps around $\epsilon\!=\!0$ at the intersections of the original Dirac spectrum and its replicas. A detailed analysis of these Floquet bands and the $k$-dependence of the gaps is given in Appendix B.

A WP with initial quasi-energy and wave number as marked by the blue dot in the zoomed in area in Fig. 2(a) (red box, $\bar{k}_{i}v_{F}/\omega=0.25$) will undergo Floquet oscillations (blue arrows) when driven through the bands by a static electric field $E<0$, Eq. (3), such that $k$ increases in time, Eq. (4). Notably, since the Floquet band maxima and hence the band width vary over $k$-space, the resulting Floquet oscillation is expected to change its amplitude but not its frequency. In our simulations we choose the field strength $E$ such that, according to Eq. (7), $T_{F}/T=2\pi\hbar/T^{2}v_{F}|eE|\simeq 20.8\pi$. We took initially Gaussian WPs of the form

[TABLE]

with width $\Delta k$ and initial central mode $\bar{k}_{i}$. We employ two complementary approaches to compute and analyze Floquet oscillations: Floquet theory and direct time-integration of the full time-dependent effective Dirac equation including the $\bf E$-field.

To compute the WP velocity within Floquet theory we start from Ehrenfest’s theorem

[TABLE]

where $\left[H(t),\hat{x}\right]=-i\hbar v_{F}\sigma_{x}$ for the Dirac case. $|\Psi(t)\rangle$ is obtained via the time-evolution operator of a Floquet system Hänggi (1998) that for a single $k$-mode reads

[TABLE]

Here $\epsilon_{\alpha}$ are the Floquet quasi-energies and $|\phi_{\alpha,k}(t)\rangle$ the corresponding eigenstates of ${\cal H}_{F}$, including replicas $n\hbar\omega$ (see Eq. (20)), each consisting of two branches $\alpha\!=\!\pm$ from the linear dispersion. The additional electric field induces a linear change of $k$, which we account for by adjusting $k(t)$ according to Eq. (4). Applying $U_{k}(t,0)$ to an initial (WP) state

[TABLE]

where $|c_{\alpha,k_{i}}|^{2}=|\langle\phi_{\alpha,k_{i}}(0)|\Psi(0)\rangle|^{2}$ describes the initial occupation of branch $\alpha$, and plugging Eq. (11) into Ehrenfest’s theorem gives

[TABLE]

Here $k(t)$ is given by Eq. (4). The occupation $|c_{\alpha,k_{i}}|^{2}$ is time-independent as long as different Floquet bands are far enough apart for Landau-Zener interband transitions Landau (1932a); Zener (1932); Landau (1932b); Stueckelberg (1933); Breid et al. (2006); Krueckl and Richter (2012) to be neglected. This is the case for the time scales $t\leq 200T$ considered below.

For the periodic pulse sequence, Eq. (8), we numerically compute $v(t)$ by means of Eq. (13) and $x(t)=\int_{0}^{t}\textrm{d}t^{\prime}v(t^{\prime})$. Due to the rectangular pulse shape we must include up to $n=500$ Floquet modes to achieve sufficient convergence. The results are shown in the left panel of Fig. 2(b) as red and green lines for $x(t)$ and $v(t)$, respectively. They indeed show distinct Floquet oscillations, as predicted, with period $T_{F}\simeq 20.5\pi T$, matching the expected value $T_{F}\simeq 20.8\pi T$ from Eq. (7). The off-diagonal velocity term ($\alpha\neq\beta$) in Eq. (13) encodes the interference of states living in different Floquet bands, giving rise to an additional feature, zitterbewegung Schrödinger (1930); Zawadzki and Rusin (2011a) caused by the Floquet Dirac band structure (see Appendix C).

Our second, complementary method to compute Floquet oscillations is based on the WP propagation algorithm “Time-dependent Quantum Transport” (TQT) Krückl (2013), see also Appendix D. The WP is discretized on a rectangular grid and the time-evolution operator for the full Hamiltonian including the drift field, $H(t)\!=\!H_{0}(\mathbf{k})\!+\!V_{T}(t)\!-\!Ex$, is computed. Then the Lanczos method Lanczos (1950) is used to evaluate the action of the time-evolution operator on the WP and to compute $\Psi(x,t)$. The time-dependent position and velocity expectation values are then calculated through

[TABLE]

and

[TABLE]

The resulting TQT data is shown in Fig. 2(b), left panel, as black and blue curves for $x(t)$ and $v(t)$, respectively. They approximately coincide with those computed within Floquet theory, also showing zitterbewegung on top of the Floquet oscillations. Moreover, the WP position $x(t)$ computed with TQT precisely reflects characteristic features of the Floquet quasi bands, shown in the red box in Fig. 2(a), namely the increasing amplitude and the sharpening of the maxima and minima although TQT directly integrates the time-dependent Schrödinger equation without using Floquet formalism.

While the static electric drift field enters into the full Hamiltonian governing the numerically exact TQT time evolution, within the Floquet approach its effect is included via the acceleration theorem (4) into the time evolution, Eqs. (11, 13). The latter approximation, together with residual numerical errors from the cutoff in the Floquet quantum number, could explain the slight deviations between Floquet and TQT data in the left panel of Fig. 2(b).

Finally, in the right panel of Fig. 2(b) we present snapshots of the absolute square $|\Psi(x,t)|^{2}$ taken at the turning points (marked as black crosses) of the red curve in the left panel. They show clear-cut Floquet oscillations of the full WP in configuration space generated for the setting of a periodic pulse sequence.

III.2 Circularly polarized light

The experimental realization of Floquet oscillations could be easier to achieve in an alternative setup, employing circularly polarized light as periodic driving. The associated vector potential $\mathbf{A}$ enters (linearly) the Dirac Hamiltonian (1) via the minimal coupling

[TABLE]

The Floquet quasi-bands of graphene illuminated by circularly polarized light have already been studied extensively Oka and Aoki (2009); Kibis (2010); Zhou and Wu (2011); Savel’ev and Alexandrov (2011). Recently, also transport Atteia et al. (2017) and topological Kelardeh et al. (2016); Delplace et al. (2013); Usaj et al. (2014); Perez-Piskunow et al. (2015) properties were investigated. Instead, here we focus on generating Floquet oscillations for realistic experimental parameters. To be closer to measurements, we explicitly treat the 2d case with an initial, radially symmetric Gaussian WP analogous to Eq. (9), with $k$ and $\bar{k}_{i}$ replaced by $\mathbf{k}$ and $\bar{\mathbf{k}}_{i}=(\bar{k}_{i},0)$. Using again TQT we simulate the WP dynamics (see Fig. 4) in presence of the electric field $\mathbf{E}=E\mathbf{e}_{x}$, where $T_{F}/T=2\pi\hbar/T^{2}v_{F}|eE|\simeq 140$.

We calculate the Floquet band structure and the WP dynamics for two different dimensionless driving strengths $\tilde{A}=v_{F}eA/(\hbar\omega)$ to demonstrate the field amplitude influence on the Floquet oscillation frequency. Figure 3(a) shows the radially symmetric Floquet band structure along the direction of the electric field. For small enough driving strengths, the distance of local band minima is independent of $k$ (black). For larger driving amplitudes, multiple-photon processes significantly alter the formation of band gaps.

As displayed in Fig. 3(b), our TQT simulations of the position expectation values, Eq. (14), of two WPs with scaled initial momenta $\bar{k}_{i}v_{F}/\omega=-2.07$ and $-1.22$ for $\tilde{A}=$ 1.1 and 0.5, respectively (marked by blue dots in panel (a)) clearly show Floquet oscillations, nicely reflecting the shape of the underlying Floquet band structure as expected. Since the gaps between unperturbed Dirac cones open in a smaller $k$-range than for the periodic pulse train (8), there are less cycles of regular Floquet oscillations. At longer times Landau-Zener transitions to neighboring Floquet bands substantially alter the WP motion. Nevertheless, Fig. 3 shows Floquet oscillations with 4 full periods. In the famous experiments of Bloch oscillations in superlattices Feldmann et al. (1992); Waschke et al. (1993); Ignatov et al. (1994) and bulk semi-conductors O. Schubert et al. (2014), their detection was possible even though only $1-3$ periods could be achieved. Moreover, oscillations involving many periods exist for the case of periodic kicking which can be experimentally realized through laser pulse trains.

In analogy to Fig. 2(b), Fig. 4(a) shows the position expectation value, Eq. (14), along the electric field, Eq. (3), and the 2d shape of the WP. Here we choose $\bar{k}_{i}v_{F}/\omega=-0.64$ and $T_{F}/T=2\pi\hbar/T^{2}v_{F}|eE|\simeq 14$ without shifting the initial time $t_{0}$ and including orange crosses to mark when the snapshots of the WP are taken. The snapshots displayed on the right show the absolute square of the WP in 2d space and the oscillation of its center of mass around $x=0$. Note that according to Eq. (7) $T_{F}$ is inversely proportional to $E$, therefore in Fig. 4 the timescale of the Floquet oscillations is a factor of 10 smaller than in Fig. 3. Here, the oscillations induced by the circulating electric field, Eq. (16), are resolved on top of the slower Floquet oscillations. In Fig. 4(b), we plot the trajectory of the WP’s center of mass in the $x$-$y$-plane. The red curve is meant to serve as a guide for the eye. In the trajectory, one can recognize the Floquet oscillations but also some additional features. The smaller oscillations have the same origin as the tiny oscillations on top of the comparably slow Floquet oscillations in Fig. 4(a): They arise from the momentum change due to the circulating electric field of the pulse, Eq. (16), and are modulated by zitterbewegung. The drift of the WP in $y$-direction however has a less intuitive explanation. We assume that it is caused by the anomalous velocity, *i.e. *the Berry curvature of the Floquet bands, since it is orthogonal to the electric driving field. This assumption will be tested in future studies.

III.3 Including trigonal warping and introducing experimental parameters

Finally, we want to give a few more details on the experimental realizability of Floquet oscillations in graphene. First of all we numerically verify that when choosing the correct energy scales trigonal warping effects do not play a role. Therefore we extend our Hamiltonian $H_{0}(\mathbf{k})$, Eq. (2), to

[TABLE]

with $v_{F}=10^{6}\,\rm m/s$ and $\mu=3a^{2}t/8$, where $a=1.42\,\rm\AA$ is the nearest-neighbor distance and $t=2.7\,\rm eV$ the nearest-neighbor hopping strength of graphene Katsnelson (2012). We set the light amplitude $A=45\,\rm nVs/m$ and the frequency $\omega/2\pi=10\,\rm THz$. We simulate the wave packet motion for $\bar{k}_{i}=0.013\,\rm 1/\AA$. The described parameters are equivalent to the unitless values chosen for Fig. 3. In Fig. 5 we compare the Floquet oscillations obtained with and without the trigonal warping term and find no qualitative differences. However, for experimental realization in graphene the transport relaxation time $\tau$ has to be higher than the period $T_{F}$ of the Floquet oscillations, Eq. (7), typically $\tau=1-20\,\rm ps$ Wang et al. (2013); Berdyugin et al. (2019). As $T_{F}\simeq 13.5\,\rm ps$ for $E=30\,\rm V/cm$ (Fig. 5 (a)), good quality samples would already allow for the observation of one period. If we increase the electric field by a factor of 10 (Fig. 5 (b)), $T_{F}\simeq 1.3\,\rm ps$. Then, all 4 Floquet oscillation cycles could be detected. Note that for the higher electric field the Floquet oscillations become slightly altered. The additional oscillations are caused by Zitterbewegung and the momentum change induced by the circularly polarized light field and are resolved here since the Floquet oscillations occur on the time-scale of pico-seconds.

IV Conclusions

The above analysis and simulations constitute a proof of principle for generating Floquet oscillations in systems with an effective Dirac dispersion. Concerning possible experimental realizations, graphene Castro Neto et al. (2009), topological insulators Hasan and Kane (2010); Qi and Zhang (2011) and cold atoms in artificial honeycomb lattices L. Tarruell et al. (2012) appear suitable Kolovsky and Bulgakov (2013); Kolovsky (2018). The latter have the advantages that one can additionally tune $v_{F}$ entering the Floquet frequency and that relaxation through disorder or interaction effects can be avoided. Other effective Dirac systems, *e.g. *for plasmons Weick et al. (2013) or polaritons Jacqmin et al. (2014), could also be considered. In the following, we will quantitatively focus on realizations for charge carriers in real monolayer graphene.

The radiation frequency $\omega$ must be chosen such that several Floquet BZs lie in the energy range governed by the linear dispersion. In graphene, a frequency $\nu=2\pi\omega$ of $1-10\,$THz is small enough to accommodate roughly 50 to 5 Floquet replicas over the energy range of $400\,$ meV, for which the linear Dirac cone is a very good approximation. Additionally, $\omega>\omega_{F}\gg 2\pi/\tau$ is required, where $\tau$ denotes a typical relaxation time of charge carriers. To realize the oscillations shown in Fig. 5(b) with $\tilde{A}=1.1$ at a radiation frequency of $\omega/2\pi=10\,\rm THz$, the moderate intensity (avoiding sample heating) of $I=\frac{c\epsilon_{0}}{2}A^{2}\omega^{2}\simeq 1\frac{\rm MW}{\rm cm^{2}}$ is needed. Here $c$ is the speed of light and $\epsilon_{0}$ the vacuum permittivity. Then an electric field $E\simeq 0.3\frac{\rm kV}{\rm cm}$ is sufficient to generate Floquet oscillations of frequency $\omega_{F}/2\pi\gtrsim 1\,\rm THz$. Hence $\omega_{F}\gtrsim 2\pi/\tau$ for typical inverse transport relaxation times $1/\tau=0.05-1\,\rm THz$ of clean hexagonal boron nitride-encapsulated graphene Wang et al. (2013); Berdyugin et al. (2019). Thus, Floquet oscillations could in principle be observed, opening an alternative way to generate THz radiation Ganichev and Prettl (2005); Lewis (2014).

To conclude we showed that free particles in a static electric drift field and obeying a linear Dirac-type dispersion can perform spatially periodic motion, Floquet oscillations, when subject to time-periodic driving. The Floquet time lattice takes on the role of the spatial lattice required for conventional Bloch oscillations. Such Floquet oscillations feature zitterbewegung and characteristic amplitude modulations that could provide a tool to experimentally map the Floquet quasi-bands. A closer consideration of Landau-Zener transitions between different Floquet bands and the question of how the topology of Floquet bands N. H. Lindner et al. (2011) is reflected in corresponding Floquet oscillations opens interesting perspectives for future research.

Acknowledgements.

We thank Simon Maier for support with calculating Floquet band structures at an early stage of this work and Sergey Ganichev for helpful discussions and careful reading of the manuscript. We further thank an anonymous referee for suggesting to interprete Floquet oscillations in terms of multiple photon exchange. We acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 314695032 – CRC 1277 (subproject A07).

Appendix A Basic relations in Floquet theory

In the following we describe how to obtain the Floquet band structure for a finite driving strength $V_{T}(t)$ without the electric field, which is later included by a shift of $k$ when computing Floquet oscillations. Generally, the Floquet operator $\mathcal{H}_{F}$ is given by

[TABLE]

where $H_{0}(\mathbf{k})$ is the Hamiltonian of the time-independent system. Since the eigenstates of $\mathcal{H}_{F}$, $\Phi_{\alpha}(\mathbf{k},t)$, are periodic in time, they can be expanded in a Fourier series:

[TABLE]

The dimension of the coefficients $\mathbf{u}_{\alpha}(\mathbf{k},n)$ is equal to the number of branches $\alpha$ of the time-independent Hamiltonian $H_{0}(\mathbf{k})$, *i.e. *two in our case. In order to find these coefficients and the corresponding eigenenergies $\epsilon_{\alpha}$, the Floquet equation (18) is multiplied by $\textrm{e}^{-im\omega t}$, $m\in\mathbb{Z}$, and averaged over one period $T$ to end up with Shirley (1965); Sambe (1973)

[TABLE]

Here the contributions from $V_{T}(t)$ are denoted by $V_{F,mn}$ and the contributions from $H_{0}(\mathbf{k})$ and $-i\hbar\partial_{t}$ by $\mathcal{H}_{0F,mn}(\mathbf{k})$. Note that the resulting eigenvalue problem is time-independent and all the dynamics of a WP in the system are incorporated in the Floquet basis states Arlinghaus and Holthaus (2011). To account for this when projecting a WP from the time-dependent basis to the Floquet basis, the relation $\langle k_{o}(t)\rangle_{t}=k_{Floquet}$ is used, where $\langle k_{o}(t)\rangle_{t}=\frac{1}{T}\int_{0}^{T}k_{o}(t)\ \textrm{d}t$. The time dependence of $k_{o}(t)$ is introduced by the time-periodic driving $V_{T}(t)$. Hereafter, we refer to $k_{Floquet}$ when talking of wave numbers and suppress its subscript to simplify the notation. The Floquet Hamiltonian without coupling

[TABLE]

is diagonal and describes the Dirac bandstructure shown in Fig. 1. The driving Hamiltonian

[TABLE]

on the other hand couples different Floquet modes $\mathbf{u}_{\alpha}(\mathbf{k},n)$ and thus can lead to the opening of band gaps in the originally linear spectrum (see Fig. 1). For the numerical evaluation, the resulting infinite matrix $\mathcal{H}_{F}$ has to be truncated at a finite value $\pm n$ that corresponds to the number of Floquet replicas taken into account. When performing numerical calculations, one has to make sure that the results are converged for the chosen value of $n$.

Appendix B Floquet quasi-band structure of a Dirac system with periodically opened mass gap

Here we give a more detailed analysis of the Floquet band structure of a Dirac system with a periodically opened mass gap. The potential describing this time-dependent pulsing is given in Eq. (8). For the mass gap, we find with Eq. (22)

[TABLE]

Qualitatively, the two time scales $\Delta t$ and $T$ involved are reflected in $k$-space. While the longer scale $T$ makes for the high-frequency oscillations of the Floquet bands due to the replicas, the smaller time scale $\Delta t$ is responsible for the slow modulation, *i.e. *the different gap sizes as function of $k$ shown in Fig. 6.

On a more quantitative level, the dependence of the Floquet bands on the pulse amplitude $M$ and pulse duration $\Delta t$ can be best understood when studying the influence of one pulse on a WP in the static Dirac cone. This has been done extensively in Ref. Reck et al. (2017) and will only be summarized here. Let us consider a WP initially occupying states in the upper cone. The opening and closing of the gap causes a redistribution of the WP to both cones, leading to a new superposition state. The amplitude of the part now occupying the other cone - in our example the lower one - can easily be calculated analytically Reck et al. (2017):

[TABLE]

where $\eta=E_{k,\pm}/M$ and $\mu=M\Delta t/\hbar$. This transition amplitude only depends on the initial energy $E_{k,\pm}$ of the state and the pulse parameters.

A numerical comparison of the corresponding transition probability $P(k)=|A(k)|^{2}$ and the Floquet band structure reveals that the gaps that open at the intersections of the repetitions of the original cone are directly related to the transition probability at that $k$-value: The larger $P(k)$, the larger the band gap. The reason for this dependence can be motivated in the following way. $P(k)$ describes for a single pulse the proportion of the WP that is transferred to the other band, *i.e. *how much a single pulse couples upper and lower band of the Dirac cone. On the other hand, the band gap of the Floquet bands is due to this coupling of (initially) linear band replicas. Therefore, it is not surprising that $P(k)$ and the band gap width are directly related.

We show this for two exemplary band structures in Fig. 6. In panel (a) we set $\mu=6.3$ and $\Delta t/T=0.5$. Since $P(k)=0$ around $kv_{F}/\omega=0$ the original Dirac cone is preserved. For larger $k$ band gaps open. The resulting complex band structure is a perfect example of how nicely the Floquet band structure can be tuned based on the transition probability $P(k)$. The band structure shown in panel (b) is the same as the one shown in Fig. 2(a) of the main paper bur for a larger range of $k$-values. There, a wide area around $kv_{F}/\omega=0$ is gapped. For an appropriate choice of parameters (as in panel (b)) a large $k$-window, in which Landau-Zener transitions are suppressed, can be chosen to support Floquet oscillations. As a rule of thumb, the smaller $\Delta t/T$, the more band gaps open and thus allow for more periods of Floquet oscillations before Landau-Zener transitions diminish them.

Appendix C Zitterbewegung in a Dirac system with periodically opened mass gap

Zitterbewegung (ZB) was originally predicted by Schrödinger for the Dirac equation Schrödinger (1930) but an analogue is also visible in multiband systems Schliemann et al. (2005); Gerritsma et al. (2010); Zawadzki and Rusin (2011b); Qu et al. (2013); LeBlanc et al. (2013); Reck et al. (2020). The reason is the interference of particle and antiparticle contributions in a WP, or respectively the contributions of different bands, in cases when the velocity operator does not commute with the Hamiltonian. The corresponding term is the offdiagonal term of the expectation value of the velocity operator and the frequency is given by the difference of the energies. In our case, we also have an effective two band model due to the Floquet bands with an offdiagonal term of the velocity as seen in Eq. (13). The frequency of the corresponding “Floquet zitterbewegung” is given by the energy gap between the two Floquet quasi-energies,

[TABLE]

In Fig. 7 we show the Fourier transform $\tilde{v}(\omega_{v})$ of the velocity of a WP starting at different $k$-values in the Floquet band structure to analyze the frequency spectrum. We denote the variable of the Fourier transform of the velocity by $\omega_{v}$. For these calculations no static electric field was applied, such that $k$ and therefore $\omega_{ZB}$ stay constant in time. Before performing the Fourier transform to investigate the frequency of the oscillations, the mean velocity value was subtracted of the corresponding data to avoid a peak at $\omega_{v}=0$.

The dashed lines mark $\omega_{ZB}$ as calculated by Eq. (25). Their good agreement with the spectrum confirms that the off-diagonal velocity in the Floquet picture describes ZB caused by the interference of states occupying different Floquet bands. Since the velocity has a rectangular shape, the peaks are repeated at higher harmonics. This rectangular shape can be explained in our example by the fact that the velocity can only change during the mass gap, which means that the harmonic oscillation of the ZB is effectively sampled with the driving frequency $\omega$.

Appendix D The c++ library ”Time-dependent Quantum Transport” (TQT)

To propagate a quantum state $|\psi\rangle$, one has to solve the time-dependent Schrödinger equation,

[TABLE]

with the Hamilton operator $\hat{H}$, which depends in general on time. Formally, it can be solved using the time-evolution operator

[TABLE]

which is unitary and fulfills

[TABLE]

where $t_{0}$ is the initial time, $t$ is some arbitrary later time and $t^{\prime}$ is a time in between. The time-evolution of a state then yields

[TABLE]

Moreover, for time-independent Hamiltonians, the time-evolution operator simplifies to

[TABLE]

On the other hand, any function can be approximated by step-wise constant functions – the smaller the steps, the better the approximation. Thus, the time-ordered exponential of Eq. (27) can be estimated by

[TABLE]

where the Hamiltonian is made step-wise constant for the time duration $\delta t$. The advantage is that instead of the time-ordered product of Eq. (27), a rather easy multiplication can be performed. Of course, one has to make sure that the time step $\delta t$ is small enough, such that the numerical result is converged.

Using the time-evolution operator, we shifted the problem of solving the differential equation in Eq. (26), to having the Hamiltonian operator in an exponential, which is defined by its (infinite) series expansion. The publicly available c++ library “Time-dependent Quantum Transport” (TQT) by Viktor Krueckl Krückl (2013) takes care of this expansion as efficient as possible for 1d or 2d systems. The expanded time-evolution operator acts on a numerically defined initial state in real space. Since an sufficiently smooth function can be approximated by its values at discrete points, the space is discretized by a grid, in 1d with $N_{x}$ and in 2d with $N_{x}\times N_{y}$ points. Thus, the wave function becomes a complex valued $N_{x}$-component vector or $(N_{x}\times N_{y})$-matrix, respectively. For this work we usually take $N_{x}=8192$ for 1d and $N_{x}\times N_{y}=8192\times 256$ for 2d systems.

The Hamiltonian can be either given as tight-binding Hamiltonian or as mixed position and momentum space representation, *i.e. *a function of both, the position and momentum operator, the latter being the Hamiltonian used in most cases for analytical calculations. In the mixed representation, instead of using the spatial derivative, the momentum operator acts in momentum space, *i.e. *the wave function is transformed by a fast Fourier transform, then the momentum operator acts as factor, and finally the inverse fast Fourier transformation is applied to get back to position space. The reason for using the Fourier transformation instead of the derivative is the numerical instability of the latter. Since the momentum operator acts several times (in higher orders $k_{i}^{n}$) in each small time step, the errors add up quickly. In this paper, only the mixed representation of position and momentum operator is used.

Due to the explicitly time-dependent Hamiltonian in our problem, a Lanczos method is used to expand the time-evolution operator Nauts and Wyatt (1983); Park and Light (1986) instead of a Chebyshev expansionTal-Ezer and Kosloff (1984); Kosloff (1994). The difference here is that instead of expanding in a fixed set of polynomials, the time-evolution operator is expanded in terms of the wave function $\psi$ itself and powers of the Hamiltonian acting on the wave function $\hat{H}^{n}\psi$. The thereby spanned subspace is a $N$-dimensional Krylov subspace $\mathcal{K}=\operatorname{span}\{\psi,\hat{H}\psi,\dots,\hat{H}^{N-1}\psi\}$, which is orthonormalized to get the basis vectors $u_{n}$ by a Gram-Schmidt procedure during the recursive creation for better numerical stability:

[TABLE]

with the overlaps $\alpha_{n}=\langle u_{n}\mid\hat{H}\mid u_{n}\rangle$ and $\beta_{n-1}=\langle u_{n-1}\mid\hat{H}\mid u_{n}\rangle$. Note that $u_{n}$ is a linear combination of powers of $\hat{H}$ acting on $\psi$, with highest order $n$.

The truncated Hamiltonian in this subspace becomes tridiagonal

[TABLE]

which can be diagonalized by conventional algorithms and enables the calculation of approximate eigenvalues of the operator $\hat{H}$ Lanczos (1950). With the matrix of eigenvectors $\mathbf{T}$ and eigenvalues $\mathbf{E}$ of the Hamiltonian in the reduced Krylov space $H_{\mathcal{K}}$, the time-evolution of one small time step is given by

[TABLE]

The expansion in the Krylov subspace is faster than a Taylor expansion Hochbruck and Lubich (1997) and for the Krylov space, a dimension $N$ in the range $10$–$40$ is usually enough. It turned out that for the calculations in this paper, the dimension of the Krylov space of $N=15$ is sufficient.

With the thus obtained time-dependent state on our discrete timeline, an arbitrary (observable) quantity like the position expectation value can be obtained as a function of the time for the propagation, which yields in our case the Floquet oscillations.