Dynamical non-linear optical response in time-periodic quantum systems

Introduction— Recent research has significantly advanced our understanding of both linear and nonlinear responses in real systems Boyd (2020); Shen et al. (2018), particularly within the contexts of topological phases Weinberg et al. (2017); Harper et al. (2020); Lu et al. (2022); Kundu et al. (2016); Seshadri and Pereg-Barnea (2023); Kitagawa et al. (2010), twisted bilayer graphene Topp et al. (2019); Jiang et al. (2024), and open quantum systems Mori (2023). These findings have important implications for quantum engineering, the experimental identification of nonlinear dynamic systems, and the study of topological edge states Mochizuki et al. (2020) in periodically driven systems New (2011). Floquet systems present a rich framework for exploring new quantum phenomena Berdanier et al. (2018), offering non-thermal behavior Seetharam et al. (2018), supporting various topologically non-trivial phases Fulga and Maksymenko (2016), and even enabling the creation of Floquet time crystals Berdanier et al. (2018). Research in Floquet systems opens new opportunities for controlling and manipulating quantum systems, with significant implications for condensed matter physics, quantum optics, and photonics Park et al. (2022); Seetharam et al. (2019).

Although the Floquet theorem traditionally applies to linear response Rudner and Lindner (2020); Oka and Kitamura (2019a); Castro et al. (2022), researchers have extended its concepts to nonlinear formalism Morimoto and Nagaosa (2016); Aversa and Sipe (1995); Jotzu et al. (2014); McIver et al. (2020), especially in DC regimes. Nonlinear transport explores the phenomenon where applying a sufficiently intense external field causes the resulting current to deviate from the linear response. Understanding this is crucial for studying particle dynamics in crystals and has significant applications across physics. The nonlinear response in time-independent real materials Chu and Telnov (2004); Di Francescantonio et al. (2024) remains an active research area with important implications for controlling quantum dynamics. The interplay between periodic driving and nonlinearity can give rise to novel phenomena, such as stability transitions of topological edge states, with potential applications in ultrafast spintronics and strongly correlated electron systems Oka and Kitamura (2019b); Rudner and Lindner (2020).

There are two primary formulations in the literature for calculating nonlinear conductivity from the Hamiltonian of static systems: the velocity gauge and the length gauge, each with its own advantages and disadvantages Taghizadeh et al. (2017); Di Stefano et al. (2019). Recently, the velocity gauge approach has been successfully applied using Feynman diagrams Parker et al. (2019). However, this method often suffers from spurious divergences at low frequencies, which can be problematic in numerical calculations with a limited number of bands Taghizadeh et al. (2017). In contrast, the length gauge, which reveals intriguing effects such as gyration current and Fermi surface phenomena Watanabe and Yanase (2021), faces challenges due to the complexity of the matrix elements of the position operator.

The linear response of Floquet systems has been extensively studied in previous works Oka and Aoki (2009); Zhou and Wu (2011); Dehghani and Mitra (2015); Kumar et al. (2020). Notably, the time dependence of linear conductivity has been addressed within the velocity gauge formulation Zhou and Wu (2011); Kumar et al. (2020). Recent experimental observations have demonstrated coherent control and significant modulation of optical nonlinearity in van der Waals layered magnetic insulators Shan et al. (2021), alongside theoretical studies on spin dynamics Gao et al. (2023). There are several important questions unanswered: What is the dominant nonlinear response to time-dependent quantum systems? How is it influenced by the bands and different optical processes, and what information can it provide about interband processes? To address these questions, we utilize the density matrix evolution formalism to examine both linear and nonlinear responses of time-dependent quantum systems using the length gauge. The dynamical tensor conductivity depends on three different frequencies. This versatile approach can be extended to arbitrary orders of perturbation and reveals significant phenomena, such as divergent AC responses to DC electric fields. Intriguingly, our findings on tensor conductivity reduce to nonlinear optical responses in certain limits, aligning well with results from previous studies Kumar et al. (2020); Zhou and Wu (2011). Additionally, our approach captures various DC photocurrents Watanabe and Yanase (2021), including shift current, injection current, gyration current, Berry dipole contributions, and intrinsic Fermi surface effects in multiband systems.

Perturbation theory— We construct a perturbation theory to predict the density matrix after applying a potential $\lambda V(t)$ to the time-dependent system. The time evolution of the density matrix is derived from the von Neumann equation:

where symbol $[A,B]=AB-BA$ is commutator. For time periodic Bloch Hamiltonian $H(t+T)=H(t)$ (where momentum index is suppressed) with time period $T=2\pi/\Omega$ the Floquet theorem suggests the form of Schrödinger equation eigenvectors as (assuming $e=\hbar=1$ hereafter): Floquet (1883) $i{{\partial}_{t}}|{{\psi}_{\alpha}}(t)\rangle=H(t)|{{\psi}_{\alpha}}(t)\rangle,$ where $|{{\psi}_{\alpha}}(t)\rangle={{e}^{-i{{\epsilon}_{\alpha}}t}}|{{\phi}_{\alpha}}(t)\rangle$ and Greek letter $\alpha$ is the band index, $|{{\phi}_{\alpha}}(t)\rangle=|{{\phi}_{\alpha}}(t+T)\rangle$ are time periodic Floquet quasi modes and ${{\epsilon}_{\alpha}}$ is quasienergy which is defined modulo $\Omega$, i.e. ${{\epsilon}_{\alpha}}\cong{{\epsilon}_{\alpha}}+\Omega$ and the momentum dependence is implicit. We restrict ${{\epsilon}_{\alpha}}$ to be in the first Floquet Brillouin zone i.e. $-\Omega<{{\epsilon}_{\alpha}}\leq\Omega$. Assuming $\lambda$ as a small constant, we make an expansion for $\rho$ in powers of $\lambda$ as $\rho={{\rho}^{[0]}}+\lambda{{\rho}^{[1]}}+{{\lambda}^{2}}{{\rho}^{[2]}}+...$. Equating the terms with the same powers of $\lambda$ leads to

where $n$ is a positive integer number with ${\rho}^{[-1]}=0$. We assume that the density matrix is diagonal in Floquet basis i.e. ${{\rho}^{[0]}}=\sum\limits_{\eta}{{{f}_{\eta}}|{{\phi}_{\eta}}(t)\rangle\langle{{\phi}_{\eta}}(t)|}$ with occupation of Floquet states ${{f}_{\eta}}$ remaining constant in time. As shown in Not , by defining $\rho_{\alpha\beta}\equiv\langle{{\phi}_{\alpha}}(t)|{{\rho}}|{{\phi}_{\beta}}(t)\rangle$, and $\epsilon_{\alpha\beta}=\epsilon_{\alpha}-\epsilon_{\beta}$, there is a recursion relation between density matrix elements of different orders as:

In the following, we find the optical tensor conductivity in the length gauge i.e. $V(t)=\mathbf{r}\cdot\mathbf{E}e^{-i\omega_{1}t}$ where $\mathbf{r}$ and $\mathbf{E}$ are the position operator and amplitude of the electric field. It is important to note that the matrix elements of $[{\bf r},\rho]$ when evaluated based on Bloch states, require careful consideration. The ${\bf r}$ can be divided into interband and intraband contributions $\mathbf{r}=\mathbf{r}^{i}+\mathbf{r}^{e}$ Aversa and Sipe (1995) and the matrix elements of the position operator between Bloch states generalized to Floquet system are:

where we have revived the momentum index $\mathbf{k}$ and ${{\xi}_{\alpha\beta}}=\langle{{\phi}_{\alpha}}(\mathbf{k},t)|i{{\partial}_{\mathbf{k}}}|{{\phi}_{\beta}}(\mathbf{k},t)\rangle$ is the generalized Berry connection. It is obvious from Eq. (4) that the intraband components are not well-behaved and demonstrate singularities. To address the singularity issue Aversa and Sipe (1995), we decompose the matrix elements of the position operator’s commutator into intraband and interband components. We extend this method to Floquet systems as follows:

Notably, ${\xi}_{\alpha\alpha}$ is not uniquely defined, but the covariant derivative ${{\left({{\rho}_{\alpha\alpha}}\right)}_{;\mathbf{k}}}$ is.

First order conductivity— Now, we develop the first-order perturbation theory by deriving the first-order density matrix. We define the perturbation as: $V(t)=Ez{{e}^{-i{{\omega}_{1}}t}}$ where $z=\mathbf{r}\cdot\mathbf{E}/E=z^{e}+z^{i}$ in Eq. (3). The ${\mathbf{r}}^{e}$ and ${\mathbf{r}}^{i}$ lead to $\rho^{[1]e}$ and $\rho^{[1]i}$, which can be evaluated using Eqs. (3) and (5)

where $z_{\alpha\beta}^{e(-j)}={1}/{T}\int_{0}^{T}{{{e}^{-ij\Omega t}}}\langle{{\phi}_{\alpha}}(t)|z^{e}|{{\phi}_{\beta}}(t)\rangle dt$ and $f_{\beta\alpha}=f_{\beta}-f_{\alpha}$. We find the expectation value of current $\langle\mathbf{J}\rangle=-\mathcal{N}\langle{{\mathbf{v}}}\rangle$ (we set the number of charges in the volume $\mathcal{N}=1$ hereafter and ${{v}^{x}}$ stands for velocity operator along the probe) and assume an expansion $\mathbf{J}={{\mathbf{J}}^{[0]}}+\lambda{{\mathbf{J}}^{[1]}}+{{\lambda}^{2}}{{\mathbf{J}}^{[2]}}+...$ where $\langle{{\mathbf{J}}^{[1]}}\rangle={{\sigma}^{[1]}}E{{e}^{-i{{\omega}_{1}}t}}$, so ${{\sigma}^{[1]}}=-\text{Tr}({{v}^{x}}{{\rho}^{[1]}})/({Ee^{-i\omega_{1}t}})$. The calculation shows that at the first order of perturbation when a probe electric field with a frequency of ${{\omega}_{1}}$ is applied to the system, the response would be at a frequency of ${{\omega}_{1}}+n\Omega,\,\,n\in\mathbb{Z}$ so ${{\sigma}^{[1]}}={{\sum\limits_{n}{e}}^{-in\Omega t}}{{\sigma}^{[1](n)}}$. As shown in Not we find

where the sum over momentum $\mathbf{k}$ indicates the integral over Brillouin zone $\sum_{\mathbf{k}}=\int d^{n}\mathbf{k}/(2\pi)^{n}$. It is beneficial to express the matrix elements of the position operator in terms of the matrix elements of the velocity operator Not as:

This is advantageous for numerical calculations as it eliminates the need to compute the derivative of Floquet quasi-modes, which would otherwise require defining a smooth gauge over the Brillouin zone.

Eq. (7) can be written in terms of only velocity matrix elements using Eq. (9) and is consistent with the results obtained previously in Kumar et al. (2020); Dehghani and Mitra (2015); Dabiri et al. (2022). Eq. (7) for $n=0$ is very similar to the result of static systems with only a difference in a substitution ${{\epsilon}_{\alpha\beta}}\to{{\epsilon}_{\alpha\beta}}+j\Omega$ and ${\mathbf{v}_{\beta\alpha}}\to\mathbf{v}_{\beta\alpha}^{(j)}$ and its interpretation is quite simple in terms of optical transitions. Fig. 1(a1) illustrates the different quasienergy bands defined in extended space. The undriven Hamiltonian has two bands, denoted by $\alpha=1,2$. We define the optical transition associated with $j$ in Eq. (7) as a ”j-photon assisted” optical transition. In Fig. 1(a1), we represent the optical transitions for different values of $j$ (0, 1, and 2) between the $\alpha=1$ and $\alpha=2$ bands. The total linear response is obtained by summing over all $j$ values and different initial sidebands.

Based on the interpretation outlined above, Eq. (7) reveals that $\sigma_{x\neq z}^{[1]e(0)}$ represents the sum of Hall conductivities for occupied bands across all sidebands. For ideal occupation, this sum is quantized according to the TKNN formula Dehghani and Mitra (2015); Kumar et al. (2020); Dabiri et al. (2022). It is proportional to the Chern number of occupied bands, $\mathcal{C}_{occ}$, as: $\sigma_{x\neq z}^{[1]e(0)}=\frac{1}{2\pi}\sum\mathcal{C}_{occ}$.

Second order conductivity— We calculate the second-order conductivity by obtaining the second-order density matrix. We assume $V(t)=E_{2}y{{e}^{-i{{\omega}_{2}}t}}$ and since $y$ and $\rho^{(1)}$ encompass the interband and intraband components as $y=y^{i}+y^{e}$ and $\rho^{[1]}=\rho^{[1]i}+\rho^{[1]e}$, subsequently, there will be four terms corresponding to $y^{e},\rho^{[1]e}$ (interband-interband), $y^{e},\rho^{[1]i}$ (interband-intraband), $y^{i},\rho^{[1]e}$ (intraband-interband) and $y^{i},\rho^{[1]i}$ (intraband-intraband) terms. Each term can be calculated separately using Eqs. (3), (5) and (S11) to find the $\rho^{[2]}$ as shown in Not .

When two probe electric fields with frequencies of ${{\omega}_{1}},{{\omega}_{2}}$ are applied to the system, the second order response would be at a frequency of ${{\omega}_{1}}+{{\omega}_{2}}+n\Omega,\,\,n\in\mathbb{Z}$, therefore, ${{\sigma}^{[2]}}={{\sum\limits_{n}{e}}^{-in\Omega t}}{{\sigma}^{[2](n)}}$. Defining $\langle{{\mathbf{J}}^{[2]}}\rangle={{\sigma}^{[2]}}EE_{2}{{e}^{-i({{\omega}_{1}}+\omega_{2})t}}$, we will find an expression for second-order conductivity as: Not

The second-order response of Floquet systems at a frequency of $\omega_{1}+\omega_{2}$ resembles that of static systems, with one key difference: the optical transitions are replaced by $j_{1}$- and $j_{2}$-photon-assisted transitions. Eq. (10) is recast to the formulas for the nonlinear response of static systems Watanabe and Yanase (2021) in the limit of time-independent Hamiltonian. Notice that these frequency-dependent features of the linear and nonlinear optical conductivity in Floquet systems are intricately linked to the driving frequency, influencing band structure, topological properties, and the balance between different optical processes.

Special limits— We will consider various scenarios. i) AC response to a DC electric field: Assume that $\omega_{1}=\omega_{2}=0$, hence the full-intraband linear and second-order conductivities diverge based on Eqs. (8) and (10). This is a real divergence (not spurious one usually arising in velocity gauge Taghizadeh et al. (2017)) and we call it Drude peak at finite frequency, since the responses in Eqs. (8) and (10) are at frequencies $n\Omega,\,\,n\in\mathbb{Z}$. This process, which is essentially the reverse of rectification, is experimentally intriguing and could have significant implications for the field of optoelectronics.

ii) Zero response of Floquet topological insulators (FTIs): We refer to FTIs as the Floquet systems having finite gaps at quasienergies $n\Omega~{}n\in\mathbb{Z}$. Then at frequencies $\omega_{1}=n\Omega$, there are no resonant optical transitions as shown in Fig. 1(a1). Let us consider for instance the first order conductivity $\sigma^{[1]e}_{xx}$ in Eq. (7). Using Eq. (9) and adding an infinitesimal constant $\eta$ to frequency $\omega_{1}$ to account for the relaxation, it can be written as

where we have used the approximation $\frac{1}{x-i\eta}\approx\frac{x}{x^{2}+\eta^{2}}+i\pi\delta(x)$ with $\delta(x)$ the Dirac delta function. Since for FTIs and $\omega_{1}=n\Omega$ the delta function in Eq. (11) is always zero, then $\operatorname{Re}\sigma_{xx}^{[1]e(0)}$ would be zero as well. This is in agreement with previous results Kumar et al. (2020); Zhou and Wu (2011), albeit this effect may be suppressed when gaps are small and the broadening factor $\eta$ is large.

iii) DC photocurrent: The DC response to an AC electric field is usually called photocurrent. It is obvious from Eqs. (7) and (8) that when an electric field at frequency of $m\Omega,~{}m\in\mathbb{Z}$ is applied to the system, then $\sigma^{[1]e(-m)}$ and $\sigma^{[1]i(-m)}$ give DC photocurrents. This is in contrast to static systems where there is no photocurrent at the first order of perturbation theory.

At the second order of perturbation, when $\omega_{1}+\omega_{2}=0$, there are DC photocurrents and all previous formulas of photocurrents Watanabe and Yanase (2021) including shift current, injection current, gyration current, Berry dipole contribution, and intrinsic Fermi surface effect are generalized Not to Floquet systems $\sigma^{[2]}\rightarrow\sigma^{[2](0)}$ which is similar to formulas for multiband static systems.

When $\omega_{1}+\omega_{2}=m\Omega,~{}m\in\mathbb{Z}$, there are also DC photocurrents which are obtained by $\sigma^{[2](-m)}$ in Eq. (10). This means that the probe electric field can make resonance with the Floquet system giving rise to DC photocurrent. It should be noted that the higher orders of nonlinear response can be derived straightforwardly using Eqs. (3) and (5).

One-dimensional example— Here we apply our finding results for a driven one-dimensional two-band system. We take as a prototype of the Su-Schrieffer-Heeger model proposed for electronic states in polyacetylene Su et al. (1979) subjected to an oscillating field. The Hamiltonian in $\mathbf{k}$ space reads as:

where $\sigma_{i},~{}~{}i\in{x,y,z}$ are Pauli matrices. The quasienergy bands are shown for $A=0.3$ in Fig. 1(a1) where the color scale indicates the physical weights of subbands defined as $W^{n}_{\alpha}=\langle\phi_{\alpha}^{(n)}|\phi_{\alpha}^{(n)}\rangle\leq 1$. For Floquet systems the occupation of Floquet states is generally not easily described by Fermi-Dirac distribution and we assume a quench occupation Zhou and Wu (2011); Kumar et al. (2020); Dabiri et al. (2022, 2024) of states which is obtained by projecting the non-driven ground states $|g_{0}\rangle$ on the Floquet states ${{f}_{\alpha}}=\sum\nolimits_{n}{|\langle{{g}_{0}}|\phi_{\alpha}^{(n)}\rangle{{|}^{2}}}$ and describes the situation where the drive is suddenly turned on. The occupations are shown in Fig. 1(a2) where there are band inversions at two resonances in the Brillouin zone.

Figures 1(b1) and (b2) depict the real and imaginary parts of linear interband conductivities according to Eq. (7). The peaks of conductivities in Fig. 1(b1) are associated with the optical transition between bands with higher physical weights and density of states which are marked by arrows emphasizing $j$-photon-assisted transition origin.

The interesting response which was not appreciated enough in previous works is the divergent full intraband responses shown in Figs. 1(c1) and (c2) as a function of time which have DC and AC components and are in the same phase.

Figures 1(d1) and (d2) show the $\sigma^{[2]ei}$ which shares some similarities to linear interband response. According to Eq. (10) it contains the derivative of distribution and shows a divergent behavior at zero frequency. In Fig. 1(d2) we mark different peaks with arrows showing the relevant $j$-photon assisted transitions which can be traced back to Fig. 1(a1). They show a sharp peak at zero frequency and a long tail at higher frequencies for coherent Drude carriers. We also observe another peak near the $\Omega/2$, since the derivative of occupations has finite value only near the resonances (see Fig. 1(a2)).

The full interband conductivity is shown in Figs. 1(e1),(e2). The peaks are attributed to van Hove singularities in the density of states, and also the group velocity, which makes a shift in the frequencies of the peaks. Seemingly, in the middle frequency, the interband and intraband contributions to nonlinear optical conductivity have opposite effects as the frequency increases.

To calculate $\sigma^{[2]ie}$ we should find a smooth gauge for Floquet quasi modes because otherwise the derivatives in the third equation of Eq. (10) will not be well-defined. The result is depicted in Figs. 1(f1) and (f2) showing more peaks and dips with respect to $\sigma^{[2]ee},\sigma^{[2]ei}$ because the $2\omega$ resonances are enabled. High heterodyne $\sigma^{(\pm 1)}$ responses to homodyne response $\sigma^{(0)}$ at low frequencies is notable in Figs. 1(d-f) which originates from the transitions near the dynamical gap at quasienergy $\epsilon=\pm\Omega/2$.

Although we assume quench occupation of Floquet states in Fig. 1, other occupations are possible like ideal occupation Dehghani and Mitra (2015); Kumar et al. (2020); Dabiri et al. (2022) of states where the lower (upper) Floquet band is full (empty) and for which the derivative of occupation is zero so $\sigma^{[1]i},\sigma^{[2]ei},\sigma^{[2]ii}$ are zero. This case is also studied and results are presented in Not . On the other hand for mean-energy occupation Zhou and Wu (2011); Dabiri et al. (2022) $f_{\alpha}=\theta(-\bar{\epsilon}_{\alpha})$ where $\bar{\epsilon}_{\alpha}=\epsilon_{\alpha}+\sum_{n}n\Omega W^{n}_{\alpha}$, the derivative of occupation to momentum is singular at resonances. The real occupation of states in experimental setups depends on the switch-on protocol and also dissipation mechanisms. There are also proposals for on-demand control of occupation of Floquet states using multifrequency drives and engineered switch-on protocol Castro et al. (2023) or coupling to Fermi and Bose baths Seetharam et al. (2015).

Conclusion— We have developed a theoretical framework for calculating the nonlinear response of Floquet systems to perturbations of any order within the length gauge. This theory yields a complex and non-trivial response when applied to time-dependent quantum systems. We have identified a strong interband-interband current in response to external light, which is intimately tied to topological properties through the Berry connection and band topology. Given recent experimental advances in ultrashort laser pulses, high-speed pump-probe experiments, and cold atoms in shaken optical lattices, our method is now poised for practical validation.

Acknowledgments— R. A. received partial funding from the Iran National Science Foundation (INSF) under project No. 4026871. S. S. D. thanks the late Seyed Mahdi Dabiri for his valuable help and support.

This supplementary material provides detailed calculations referenced in the main text. We present the conductivities in the two-band model limit, along with the conductivities expressed in terms of velocity matrix elements. Additionally, we include a discussion of the covariant derivative and derive the DC photo-conductivities, incorporating the Berry curvature dipole term, injection current, and shift current for Floquet systems. Finally, we present the results of the optical conductivities for a one-dimensional driven model with the ideal occupation of Floquet states.