Nonlinear dynamics of topological Dirac fermions in 2D spin-orbit coupled materials

The graphene family materials are two-dimensional staggered monolayers with a gapped energy band structure due to intrinsic spin-orbit coupling. The mass gaps in these materials can be manipulated on-demand via biasing with a static electric field, an off-resonance circularly polarized laser, or an exchange interaction field, allowing the monolayer to be driven through a multitude of topological phase transitions. We investigate the dynamics of spin-orbit coupled graphene family materials to unveil topological phase transition fingerprints embedded in the nonlinear regime and show how these signatures manifest in the nonlinear Kerr effect and in third-harmonic generation processes. We show that the resonant nonlinear spectral response of topological fermions can be traced to specific Dirac cones in these materials, enabling characterization of topological invariants in any phase by detecting the cross-polarized component of the electromagnetic field. By shedding light on the unique processes involved in harmonic generation via topological phenomena our findings open an encouraging path towards the development of novel nonlinear systems based on two-dimensional semiconductors of the graphene family.


Results
Topological properties of the graphene family. Let us consider the following Dirac-like Hamiltonian describing staggered GFM in the low-energy regime 19 , where d η s = {η v F k x , v F k y , � η s } describes a meron structure in momentum space (Fig. 1c,d), p = (k x , k y ) is the particle momentum, τ i are the Pauli matrices, η, s = ±1 are valley and spin indexes, v F is the Fermi velocity, and the mass term � η s = ηs SO − E − η L + s M corresponds to half energy band gap. Here, SO represents the intrinsic spin-orbit coupling energy and the last three terms in � η s account for interactions with external fields, which allow for tailoring the Dirac mass gap. Indeed, the second term E corresponds to the potential difference between sub-lattices in the buckled structure in the presence of a static electric field E z applied perpendicularly to the plane of the monolayer 41 . The third component L describes the anomalous quantum Hall effect and arises due to the coupling between an off-resonant circularly polarized laser and the GFM 45,52 . The final term M depicts the staggered antiferromagnetic exchange interaction 53 . For simplicity, we have neglected Rashba couplings due to their small magnitude compared to the other contributions. It is worth mentioning that the Hamiltonian in Eq.
(1) H η s = τ · d η s = v F (ηk x τ x + k y τ y ) + � η s τ z , www.nature.com/scientificreports/ (1) applies more generally to other 2D systems, including antiferromagnetic manganese chalcogenophosphates (MnPX 3 , X = S, Se) 54 and perovskites 55 . The formalism developed here is also valid for these systems. Within the Dirac picture the electronic properties of the GFM can be fully characterized through a set of four topological invariants 19 , namely the Chern C = η,s C η s , spin Chern C s = η,s s C η s /2 , valley Chern C η = η,s η C η s , and spin-valley Chern C sη = η,s ηs C η s /2 numbers. Here C η s = η sign[� η s ]/2 is the Pontryagin number 10,11 , a topological quantity that counts how many times the vector d η s /|d η s | wraps a sphere in momentum space (Fig. 1c,d). By varying the external parameters E , L , and M the monolayer can be driven through various phase transitions as depicted in Fig. 1b, where we show the topological phases in the planes E = 0 , L = 0 , and M = 0 . The electronic states quantum spin Hall insulator (QSHI), anomalous quantum Hall insulator (AQHI), band insulator (BI), antiferromagnetic insulator (AFI), and polarized-spin quantum Hall insulator (PS-QHI) have non-zero mass gaps for all four Dirac cones. The system undergoes a topological phase transition when the band gap of at least one Dirac cone changes sign. Hence, the boundaries between insulating states are determined by the condition � η s = 0 , which defines the so called single Dirac cone phases (diagonal solid lines in Fig. 1b). At the intersection between boundary lines there are two Dirac gaps closed, and the system is either in a spin, valley, or spin-valley polarized semimetal state. Finally, note that there are four points in the 3D phase diagram, where three Dirac cones simultaneously close, corresponding to unique electronic states that remain unexplored to date.
Nonlinear light-matter interactions in the GFM. The nonlinear dynamics of charge carriers in a given Dirac cone interacting with a linearly polarized electromagnetic plane wave impinging normally on the monolayer can be described via the density matrix ρ(t) , which satisfies the equation of motion Here, Ĥ i (t) = ep · A(t)/mc is the field-monolayer interaction Hamiltonian in the velocity gauge 56 , where A(t) = A(t)x is the electromagnetic vector potential. We assume that at large times the system relaxes to the equilibrium density matrix ρ (0) = l,k f lk |lk��lk| with a phenomenological decay rate Ŵ , where f lk is the Fermi-Dirac distribution for fermions in the valence or conduction bands ( l = ±1 ) with momentum k , and |lk� are the corresponding eigenstates of the unperturbed Hamiltonian H η s (see Methods). In general, the equation of motion cannot be solved analytically due to the time dependence of the interaction Hamiltonian, and numerical methods are often employed. Here, instead, we use perturbation theory by expanding the density matrix ρ(t) = nρ (n) (t) in powers of the vector potential amplitude and iteratively solve for ρ (n) (t) ∝ A(t) n . We use the standard expression for total current j(t) = −e T r[ρ(t)(v + e mcÂ (t))] , which contains contributions arising from all orders in the perturbative expansion 57,58 . Here, v = −1 ∇ kĤ η s is the velocity operator. The n-th order ac electric current is given by Although the last term in the previous equation does not vanish in general, one can demonstrate that it is exactly canceled by another contribution arising from the first term and, therefore, its contributions do not appear in the final expressions for the conductivities. The cross-canceling of contributions involving the first and second terms in Eq. (2) is actually critical to reobtain the well know expression for the conductivity in the linear regime, given by Kubo's formalism (Methods). The n-th order optical conductivity is a tensor of rank n + 1 , accounts for both intra and interband transition contributions as well as topological effects emerging from the Berry connection 19 , and can be directly computed after Fourier transforming the above equation and expressing the current as a product between the electric field and conductivity.
The longitudinal σ (1)η,s xx (ω) and transverse σ (1)η,s yx (ω) linear optical conductivities per Dirac cone obtained from Eq. (2) are in full accordance with previous studies, and have been investigated in details in the M = 0 plane in Refs. 33,41,47,49 . For completeness, we explore in Fig. 2 the charge, spin, valley, and spin-valley transverse linear conductivities (summed over valley and spin indexes as shown in the Methods) in the four phases with only one Dirac cone open, since these have been previously overlooked in the literature. By measuring the dynamic linear conductivities at these unique points in phase space one can investigate topological properties emerging from each Dirac cone individually. Indeed, the nonzero Dirac gap in each of the cases considered leads to Chern numbers C , C η = ±1/2 , and C s , C sη = ±1/4 , which clearly manifest in the zero frequency limit of the transverse conductivities for neutral monolayers, as shown in Fig. 2. Note that in panels 2a,b the charge (spin) and the valley (spin-valley) currents are identical, which implies that the nonzero Dirac gap belongs to the K valley, while in panels 2c,d the charge (spin) and valley (spin-valley) conductivities are the negative of each other, hence the nonzero Dirac gap belongs to K ′ valley. Finally, the sign of the spin or spin-valley linear conductivities in the quasi-static regime allows to identify the spin of the charge carriers.
The second order response of the system described by the Hamiltonian in Eq. (1) is zero. The first nonlinear contribution comes from the third-order response and can be described by a rank four conductivity tensor σ (3)η,s αβγ δ (ω 1 , ω 2 , ω 3 ) , which is invariant under simultaneous permutations of the indexes β, γ , δ and frequencies ω 1 , ω 2 , ω 3 , where α, β, γ , δ = x, y (see Methods). For simplicity we will concentrate in cases where β = γ = δ = x and we will assume that the incident electromagnetic radiation is monochromatic, i.e., A(t) = A 0 cos(ωt)x , allowing us to separated the third-order conductivity into terms oscillating with frequencies ω and 3ω . The first term adds to the linear conductivity response and gives a correction to the fundamental harmonic that is quadratic in the vector potential amplitude (Kerr effect), while the second one describes the nonlinear process of third harmonic generation. Therefore, up to third-order in perturbation theory, the optoelectronic conductivity of the system at the fundamental and third-harmonics due to each Dirac mass gap can be cast as (3)η,s αxxx (ω, ω, ω) is a shorthand notation. In the next section we discuss topological fingerprints embedded in these nonlinear conductivities.
Topological signatures in the GFM nonlinear response. Let us first look into topological signatures buried in the phase of the optical field by investigating the polarization state (helicity) of light resulting from third-harmonic generation processes. We consider the low temperature regime for a neutral monolayer, i.e., the Fermi level lies in the middle of the band gap, and assume ω = 0.2 SO and Ŵ = 0.05 SO . In Fig. 3a we plot the phase diagram associated with the difference h 3ω = (I (3) yx (3ω) between intensities I   Finally, we mention that the variations of h 3ω across a phase transition depend on whether the initial and final phases have trivial charge transport properties (Chern number C equal to zero). For instance, in the M = 0 we observe that h 3ω is a symmetric function around E / SO = 1 when we consider the topological phase transition between the QSHI and BI phases, both with C = 0 . On the other hand, a remarkably asymmetric behavior takes place around L / SO = 1 (which is energetically equivalent to the previous case) when one considers a topological phase transition involving the QSHI and the AQHI ( C = −2 ) phases.
The phase diagram of h 3ω includes contributions from both third-harmonic longitudinal and transverse conductivities. However, only the latter encodes topological features of the GFM and deserves to be independently studied. To this end one could, for instance, use a linear polarizer filter to block the co-polarized scattered field and detect only its cross-polarized component, which is proportional to σ (3) yx . In Fig. 4 we investigate the charge conductivity σ In panel 4a we plot the conductivities along the E axes as we move from the QSHI to the BI phase by closing two Dirac gaps when E / SO = 1 . In this case the Chern numbers ( C, C s , C η ) change from (0, 1, 0) to ( 0, 0, −2 ) and we observe that the charge conductivity is always zero. The resonant peaks and dips near E / SO = 1 in the spin and valley conductivities divulge information about the change in C s and C η across the phase transition. For instance, as the system is driven from the QSHI to the BI phase the resonance at 3 ω = 2� yx (3ω) flips from negative to positive values, which is a signature indicating a decrease in the corresponding Chern numbers across the phase transition. Also, the magnitude of the resonance peaks and dips in the valley conductivity is twice as those of the spin conductivity, which coincides with the fact that the change in C η is two times that in C s . These features are general and can be applied to any paths in the phase diagram regardless of the nature of the transition involved.
In panel 4b we plot the conductivities as the system moves from the AQHI to the BI phase via the PS-QHI phase. Along this path the monolayer crosses two phase boundaries, with Chern numbers ( C, C s , C η ) changing as (−2, 0, 0) → ( −1, 1/2, −1)→ ( 0, 0, −2 ). Note that the signs of the resonances at 3 ω = 2� η s for the valley conductivity are opposite to those for the charge conductivity, while their strengths are the same, meaning that the variations in C and C η across the phase transition are equal in magnitude. The sign of the third harmonic resonance for the charge (valley) conductivity changes from negative (positive) in the AQHI phase to positive (negative) in the PS-QHI phase, which agrees with an increase (decrease) in C ( C η ). The same feature repeats near the phase transition between the PS-QHI and BI phases, signaling another increase (decrease) in C ( C η ). www.nature.com/scientificreports/ Moreover, the behavior of the spin conductivity in Fig. 4b reflects the fact that C s increases along the first phase boundary and then decreases along the second transition, further confirming that the nonlinear resonant behavior of the conductivities encode topology fingerprints.
In the third panel of Fig. 4 the system is driven along the phase boundaries highlighted in path (c), where at least one of the Dirac gaps is closed all the time and the Chern numbers change according to (−3/2, 1/4, −1/2) → ( −1/2, 3/4, −1/2 ) → ( −1/2, 1/4, −3/2 ). It is important to notice that C does not change near E / SO = 1 and C η remains the same across the transition at L / SO = 1 ( E = 0 in the parametrized curve). This fact is clearly reflected in Fig. 4c, where the third harmonic resonances associated with charge (valley) conductivity have the same signs as the system is driven through E / SO = 1 ( L / SO = 1 ). Finally, in Fig. 4d we plot the conductivities for a path (not shown in Fig. 4e) going from the origin of the phase diagram to the point { E / SO , L / SO , M / SO } = {2, 2, 2} across the diagonal E = L = M . Once again, the behavior of the conductivities and the change in Chern numbers agrees well with signatures we explained above.
Next we show how one can identify signatures of the Chern numbers in each topological phase by investigating the nonlinear spectral response of the GFM. In Fig. 5 we present the frequency dispersion of the real part of the third-order optoelectronic longitudinal and transverse conductivities for both fundamental and third harmonics at four different points in phase space. These correspond to cases where all gaps are open (Fig. 5a) or one (Fig. 5b), two (Fig. 5c), or three (Fig. 5d) Dirac gaps are closed. The corresponding curves for the imaginary part (not shown) of the conductivities can be obtained directly from our formalism or via generalized Kramers-Kronig relations for nonlinear systems 59 . We expect that the nonlinear conductivities should have non-negligible contributions when at least one Dirac gap matches one or three times the energy of the incident photons. Indeed, Fig. 5 shows that Kerr σ  Note that near resonances the longitudinal Kerr and third harmonic conductivities always shows normal and anomalous dispersion, respectively, regardless of the resonant Dirac cone. This is due to the fact that both of these conductivities depend only on the absolute value of the resonant mass gap, hence not allowing to directly distinguish between energetically equivalent but topologically different points in the phase space. On the other hand, the transverse Hall conductivity is sensitive to the sign of the resonant � η s and, therefore, encodes information of the corresponding Pontryagin number C η s . Note in Fig. 5, for instance, that σ yx (3ω) precisely at the resonant frequencies, accounting appropriately for degeneracy among mass gaps, and multiplying the result by +1/2 or −1/2 , respectively. Consider, for example, the case in Fig. 5a. Note that in this case the third-order conductivity at the fundamental (third) harmonic evaluates to negative (positive) values in three of the four resonances. Thus, the sum of their signs equals −2 ( +2 ) and results in a Chern number C = −1 , as expected for the PS-QHI phase in the M = 0 plane. Note that this argument can be extended to the nonlinear spin, valley, and spin-valley Hall currents, allowing for obtaining the full set of Chern

Conclusions
In summary, we developed a comprehensive investigation of nonlinear light-matter interactions in spin-orbit coupled graphene materials, with focus on third harmonic generation and nonlinear corrections to the fundamental harmonic (Kerr effect). We have shown that by controlling various external parameters the monolayer can be driven through different phase transitions to explore the interplay between nonlinear effects and topological properties of two-dimensional materials. Specifically, we demonstrated that the helicity of the nonlinear scattered field encodes information about the location of the phase transitions boundaries, the dependence of the cross-polarized component of the field (proportional to the transverse conductivity) on the external parameters near resonances divulges details about the topological phases involved in the transition, and the dispersive nonlinear response of the system enables characterization of the Chern number in each phase. The recent progress in synthesis of topological semiconductors of the graphene family materials [12][13][14][15] and well established nonlinear characterization photonic techniques [23][24][25][26][27][28] suggest that our findings can be accessed experimentally with current technologies. We envision that the the effects predicted here will greatly impact research at the crossroads between nonlinear optics, topological materials, spintronics, and valleytronics.

Methods
In order to compute the nonlinear conductivity we need the unperturbed eigenvectors and the velocity matrix elements associated to the Hamiltonian in Eq. (1). In the following we will omit the valley and spin indexes η and s whenever possible, and the reader should keep in mind that, unless otherwise stated, all results presented are valid for a single Dirac cone. The total conductivities are given by summing the contributions of all mass gaps in the band structure. The eigenvector associated to the Hamiltonian in Eq. (1) has the form www.nature.com/scientificreports/ where |k| = k 2 x + k 2 y , ǫ l = lǫ = l 2 v 2 F |k| 2 + � 2 , l is the conduction ( l = +1 ) and valence ( l = −1 ) band index, and θ = tan −1 (k y /k x ) . The matrix elements �lk|v α |l ′ k� of the velocity operator v α = −1 ∂Ĥ η s /∂k α are given by We give a brief derivation for both the linear and third order optical conductivity terms for a particular Dirac cone. To this end, the equation of motion i ∂ρ(t)/∂t = Ĥ η s +Ĥ i (t),ρ(t) for the density matrix operator is solved perturbatively by expanding the density matrix in powers of the potential vector amplitude ρ(t) = nρ (n) (t) with ρ (n) (t) ∝ A n . To leading order in the potential vector, the first order density matrix upon Fourier transformation to the frequency domain can be cast as where E α (ω) = iωA α (ω)/c . Substituting ρ (1) (ω) in equation (2) and making use of the identity , we eliminate the first term in the current ∝ Tr(ρ (0) ) and re-obtain the well know Kubo formula after expressing the current as the product of conductivity and electric field where σ 0 = αc/4 and the phenomenological relaxation rate is accounted for by replacing ω → ω + iŴ . Equation (7) accounts for both interband and intraband transitions. The interband contribution is straightforward and follows directly from (7) by setting l ′ � = l . The intraband contribution can be derived by enforcing l ′ = l , in which case (f l ′ − f l )/(ǫ l ′ − ǫ l ) → ∂f lk /∂ǫ l . The total charge, spin, and valley linear conductivities are then computed as σ For the third order conductivities we first solve the equation of motion for third order correction ρ (3) (t) , which is expressed in the Fourier space as Substituting the above expression in Eq. (2) and using the same identity as before with the replacement ω → ω 1 + ω 2 + ω 3 , we obtain the third order optical conductivity per Dirac cone as The phenomenological relaxation rate can be accounted in a similar manner as in the linear case. Note that, Eq. (9) includes contributions from Berry connection, as well as intra and interband transitions 56,57 . The nonlinear contributions to the fundamental harmonic associated to the Kerr effect follow by setting ω 1 = ω 2 = −ω 3 , while the third harmonic conductivity is can be derived for ω 1 = ω 2 = ω 3 . Similarly to the linear case the total charge, spin, valley, and spin-valley conductivities due to all Dirac cones are given by σ �lk|ρ (3) (ω)|l ′ k� = e 3 E α (ω 1 )E β (ω 2 )E γ (ω 3 )/ω 1 ω 2 ω 3 i(ǫ l ′ − ǫ l + ω 1 + ω 2 + ω 3 ) l ′′ ,l ′′′ �lk|v δ |l ′′′ k��l ′′′ k|v γ |l ′′ k��l ′′ k|v β |l ′ k� ǫ l ′ − ǫ l ′′′ + ω 1 + ω 2 �lk|v δ |l ′′′ k��l ′′′ k|v γ |l ′′ k��l ′′ k|v β |l ′ k� ǫ l ′′ − ǫ l + ω 2 + ω 3 σ (3) αβγ δ (ω 1 , ω 2 , ω 3 ) = −e (ω 1 + ω 2 + ω 3 ) E α (ω 1 )E β (ω 2 )E γ (ω 3 ) × × k,l,l ′ �lk|v α |l ′ k� ǫ l ′ − ǫ l �l ′ k|ρ (3) (ω)|lk�.

Data availability
The data that support the plots within this paper and other findings of this study are available from the corresponding author on reasonable request.