Nonlinear Optical Conductivity in Graphene and other 2-Band 2-D Materials

TL;DR

This paper develops a quantum master equation framework to analyze the nonlinear optical conductivity of two-dimensional materials like graphene, silicene, and MoS2, highlighting their unique electronic responses to time-dependent electric fields.

Contribution

It introduces a first-principles quantum master equation approach to study nonlinear optical responses in 2D Dirac materials, providing detailed analysis for graphene, silicene, and MoS2.

Findings

Derived nonlinear optical conductivity expressions for graphene, silicene, and MoS2.

Compared theoretical results with existing data, confirming the model's validity.

Highlighted the impact of microscopic interactions on nonlinear responses.

Abstract

Graphene, Silicene, $\mathrm{MoS}_2$ and other similar two-dimensional structures have unusual electronic properties that lend themselves to exotic device applications. These properties emanate from the fact that the electrons are endowed with Dirac fermion-like attributes. Thus these materials are not only characterized by certain fundamental principles, they also have amazing practical uses. Our emphasis here is on one such basic property concerning nonlinear response to time-dependent electric fields. We set up a first principle quantum master equation for the underlying density operator which is based on microscopic interactions between the Dirac electron with phonons and other electrons. While such an equation has general applicability to a variety of non-equilibrium phenomena in two-band systems, we focus onto the case of nonlinear optical conductivity. The derived results are separately analyzed for graphene, silicene and $\mathrm{MoS}_2$, and comparison made with other known results.

Figures

Equations

$${\mathcal{H}}_t={\mathcal{H}}_0(\sigma ,s)+{\mathcal{H}}_I(\sigma ,s)\sin (\omega t)$$

$${\mathcal{H}}_0=-J\sum _{(iA,jB)}({|\uparrow ⟩}_i{⟨\downarrow |}_j+{|\downarrow ⟩}_i{⟨\uparrow |}_j),$$

$${\mathcal{H}}_0=-\sum _k[J(k)|\uparrow ,k⟩⟨\downarrow ,k|+J^{\ast }(k)|\downarrow ,k⟩⟨\uparrow ,k|].$$

$${\mathcal{H}}_0=v_F(\vec{\sigma }.\vec{k})=v_{F}k(\tau {\sigma }_x\cos {\chi }_k+{\sigma }_y\sin {\chi }_k),$$

$$|c_{k\tau }⟩=\frac{1}{\sqrt{2}}\left[\begin{array}{c}\tau e^{-i\tau {\chi }_k/2}\\ e^{i\tau {\chi }_k/2}\end{array}\right]$$

$$|v_{k\tau }⟩=\frac{1}{\sqrt{2}}\left[\begin{array}{c}\tau e^{-i\tau {\chi }_k/2}\\ -e^{i\tau {\chi }_k/2}.\end{array}\right]$$

$${\mathcal{H}}_t=v_F(\vec{\sigma }.\vec{k})+\tau {\sigma }_x(e{v}_{F}E/\omega )\sin (\omega t),$$

$${\mathcal{H}}_I=\tau \frac{e{v}_{F}E}{\omega }{\sigma }_x.$$

$${\mathcal{H}}_0=v_F(\tau k_x{\sigma }_x-k_y{\sigma }_y)-\tau {\lambda }_{SO}{s}_z{\sigma }_z+l{E}_z{\sigma }_z,$$

$${\mathcal{H}}_0^s=v_F(\tau k_x{\sigma }_x-k_y{\sigma }_y)-s\tau {\lambda }_{SO}{\sigma }_z+l{E}_z{\sigma }_z.$$

$$|c_{s\tau }⟩=\frac{1}{\sqrt{2|{ϵ}_{s\tau }|{\eta }_{s\tau }}}\left[\begin{array}{c}\tau k{v}_F{e}^{i\tau {\chi }_k}\\ {\eta }_{s\tau }\end{array}\right],$$

$$|v_{s\tau }⟩=\frac{1}{\sqrt{2|{ϵ}_{s\tau }|{\bar{\eta }}_{s\tau }}}\left[\begin{array}{c}\tau k{v}_F{e}^{i\tau {\chi }_k}\\ -{\bar{\eta }}_{s\tau }\end{array}\right],$$

$${\mathcal{H}}_0=v_F(\tau k_x{\sigma }_x+k_y{\sigma }_y)+\frac{\Delta }{2}{\sigma }_z-\lambda \tau \frac{{\sigma }_z-1}{2}{s}_z,$$

$$|c_{s\tau }⟩=\frac{1}{\sqrt{2|{ϵ}_{s\tau }|{\eta }_{s\tau }}}\left[\begin{array}{c}\tau k{v}_F{e}^{-i\tau {\chi }_k}\\ {\eta }_{s\tau }\end{array}\right],$$

$$|v_{s\tau }⟩=\frac{1}{\sqrt{2|{ϵ}_{s\tau }|{\bar{\eta }}_{s\tau }}}\left[\begin{array}{c}\tau k{v}_F{e}^{-i\tau {\chi }_k}\\ -{\bar{\eta }}_{s\tau }\end{array}\right],$$

$$\frac{\partial \rho }{\partial t}=-i[{\mathcal{H}}_t,{\rho }_t].$$

$$\frac{\partial \rho (t)}{\partial t}=\frac{\partial {\rho }_S(t)}{\partial t}+\frac{\partial {\rho }_D(t)}{\partial t},$$

$${\mathcal{H}}_0=\left(\frac{{ϵ}_c+{ϵ}_v}{2}\right)I+\frac{{ϵ}_c-{ϵ}_v}{2}(|c⟩⟨c|-|v⟩⟨v|),$$

$${\mathcal{H}}_0=k{v}_F(|c⟩⟨c|-|v⟩⟨v|).$$

$${\tilde{\rho }}_S(t)=e^{i{\mathcal{H}}_{0}t}{\rho }_S(t)e^{-i{\mathcal{H}}_{0}t}$$

$$\frac{\partial {\tilde{\rho }}_S}{\partial t}=-i[{\tilde{\mathcal{H}}}_I,{\tilde{\rho }}_S(t)]\sin (\omega t),$$

$${\tilde{\mathcal{H}}}_I=\tau \frac{e{v}_{F}E}{\omega }{e}^{i{\mathcal{H}}_{0}t}{\sigma }_x{e}^{-i{\mathcal{H}}_{0}t}.$$

$${\mathcal{H}}_{eff.}(t)=\tau \frac{e{v}_{F}E}{4\omega }[i⟨c|{\sigma }_x|v⟩|c⟩⟨v|e^{-i(\omega -{\omega }_{cv})t}+h.c.].$$

$$\frac{\partial {\tilde{\rho }}_S(t)}{\partial t}=-i[{\mathcal{H}}_{eff.}(t),{\tilde{\rho }}_S(t)],$$

$${\mathcal{H}}_{eff}(t)=\frac{{\Omega }_k}{2}|c⟩⟨v|e^{-i{\Delta }_{k}t}+h.c.,$$

$$\frac{\partial {\tilde{n}}_k(t)}{\partial t}$$

$$\frac{\partial {\tilde{\Pi }}_k(t)}{\partial t}=\frac{1}{4}⟨v|{\mathcal{H}}_{\mathcal{I}}|c⟩{\tilde{n}}_k(t)e^{i{\Delta }_{k}t}.$$

$${\tilde{n}}_k(t)={\tilde{\rho }}_{c_k{c}_k}(t)-{\tilde{\rho }}_{v_k{v}_k}(t),$$

$${\tilde{\Pi }}_k^{\ast }(t)={\tilde{\rho }}_{c_k{v}_k}(t),{\tilde{\Pi }}_k(t)={\tilde{\rho }}_{v_k{c}_k}(t).$$

$$⟨v|{\mathcal{H}}_I|c⟩$$

Full text

I 1. Introduction

Two dimensional electronic materials such as graphene, silicene, $\mathrm{MoS_{2}}$ etc., notwithstanding their myriad technological and device applications, are a gold mine for testing theoretical concepts of great contemporary interest Geim1 ; Geim2 ; katsnelson ; Houssa . These involve novel ideas of the chemistry of hybridized carbon orbitals saito , ultra high mobility bolotin ; du , spin-orbit interaction ezawa , Andreev reflection and Klein tunneling beenakker1 ; beenakker2 ; stander , magneto-resistance and weak localization falko ; tikhonenko ; chen , quantum Hall effect kim ; novoselov ; du2 ; bolotin2 , spintronics han ; inoue , valleytronics mak and so on. Much of the extraordinary properties of these materials manifest themselves on account of their electrons behaving as Dirac Fermions in the low energy physics sector. Indeed, all the interesting physics emanates from the underlying Dirac equation, though the electrons travel with not the speed of light but much smaller Fermi velocity. Thus these materials have turned out to be a laboratory for the realization of relativistic quantum mechanics, albeit in a non-relativistic setting of solid state physics katsnelson .

Given this background to all the extraordinary phenomena in graphene, and graphene-like materials, one aspect that is going to particularly occupy our attention in the present study is the breakdown of linear response theory in connection to a transport property such as the frequency-dependent dynamic conductivity mishchenko . The necessity to go beyond the usual ‘Kubo regime’ hinges on the fact that our interest is in the systems for which the Fermi energy $\epsilon_{F}$ is vanishingly small. In our calculation of the nonlinear electrical conductivity we follow a density matrix approach for a general two band model mishchenko ; ishikawa ; agarwal . But, in a departure from the existing treatments in which dissipative terms are incorporated phenomenologically, we employ a first-principles approach wherein these terms occur in a natural manner from the considerations of probability conservation and detailed-balance. Further, our approach is couched in the framework of quantum dissipation of two-state systems in contact with phonon and electronic baths weiss11 ; SDSP . Though the quantum baths are not treated in detail here, the indicated method opens up possibilities for an elaborate analysis of dissipation in graphene and graphene-like systems.

An outline of the paper is as follows. In section 2, we introduce the relevant Hamiltonian for graphene, silicene and $\mathrm{MoS}_{2}$ within the product Hilbert space of two distinct spin ‘variables’ – one operating within the two sub-lattice states of the material, and the other within the actual spin-states of the electron. The form of the Hamiltonian is illustrated by citing three different materials of grapheme, silicene and $\mathrm{MoS}_{2}$, in sub-sections 2a, 2b and 2c, respectively. Section 3 deals with the von-Neuman – Liouville equation for the density operator. That equation has two parts, one governed by systematic Hamiltonian dynamics in the presence of an external electric field, and the other concerned with dissipative dynamics triggered by the surrounding heat bath. These issues are separately treated in sections 3a and 3b respectively. While discussing the systematic dynamics in Section 3a, we employ a ‘Rotating Wave Approximation’ (RWA) in which rapidly oscillating terms that have no influence in the non-transient regime, are dropped. Section 3b is devoted to standard rate equation method in which the need for probability conservation and detailed-balance automatically emerge from microscopic considerations. While section 3 contains a general framework for dealing with dissipative effects in the 2-d materials at hand, we specialize in section 4 to the calculation of nonlinear optical conductivity. Once again, separate results are presented for the three systems of graphene, silicene and $\mathrm{MoS}_{2}$. These results are analyzed in section 5 as a prelude to offering certain concluding remarks in section 6 .

II 2. The Dirac Hamiltonian

The generic Hamiltonian of our focus is an explicitly time dependent one given by  mishchenko ,

[TABLE]

where $\mathcal{H}_{0}({\sigma},{s)}$ is an unperturbed Hamiltonian that lives in the product space of the two distinct dynamical variables $\sigma$- a pseudo spin and $s$- the real spin of the electron. The time $t$ is inserted as a suffix to underscore the fact that the underlying time dependence is only a parametric (and not a dynamic) one. The meaning of $\sigma$ and $s$ is clarified in the context of three different electronic materials, below. The time dependent interaction term arises from an externally applied electronic field of amplitude $E$ and frequency $\omega$. Evidently, $\mathcal{H_{I}}$ being off-diagonal in the representation in which $H_{0}$ is diagonal, would cause quantum transitions amongst the eigenstates of $H_{0}$. These transitions would imply gain or loss in energy of the unperturbed system which would have to be compensated by corresponding transitions in the surrounding heat bath (treated in section 3 below), thus leading to dissipation. One measurable quantity that characterizes dissipation is a non-linear electrical conductivity which has to be calculated to all orders in $E$ and $\omega$. The concrete case study of Eq. (1) can be effected in the context of three different materials, as described in the following sub-sections.

II.1 2a. Graphene

As our first example, we consider graphene which has a 2-dimensional honeycomb structure as shown in Fig. (1), It is evident that there are two inter-penetrating triangular lattices marked A and B. If only nearest neighbour tunneling is allowed, an electron can only jump from an A-site to a B-site and vice-versa. The corresponding tight-binding Hamiltonian may be written as

[TABLE]

where $J$ is the (constant) tunneling (or hopping) energy, $\ket{\uparrow}(\ket{\downarrow})$ denotes presence (absence) at the site A (B). The above can be rewritten in the Bloch-Fourier space as

[TABLE]

Here the wave-vector dependent tunneling $J(k)$ is complex that assumes a linear form near the $K$ and $K^{\prime}$-points of the Brillouin zone in the reciprocal space (Fig.  1). Thus the Hamiltonian for a given $k(=p,\,\hbar=1)$ has the following $2\times 2$ representation in the space of the pseudo-spin $\sigma$

[TABLE]

where $\ket{\uparrow}$ and $\ket{\downarrow}$ are the eigenstates of $\sigma_{z}$ with eigenvalues + 1 and -1 respectively.

In the above equation, $\tau=\pm 1$ for $K/K^{\prime}$ valleys, $\chi_{k}$ is the angle made by the 2-dimensional vector $\vec{k}$ with the $X$-axis and $v_{F}(=3Jd/{2\hbar})$ is the Fermi velocity that is estimated to be of the order of $10^{6}\text{m/sec}$, $d$ being the lattice parameter (Fig. 1) . Comparing Eq. (4) with Eq. (1) we note that $\mathcal{H}_{0}$ is independent of the spin degree of freedom s, for graphene. Thus the eigenstates of $\mathcal{H}_{0}$ are the electronic band states associated with the conduction band ($\ket{c_{k}}$) and the valence band ($\ket{v_{k}}$), thus

[TABLE]

and

[TABLE]

The equations (5) and (6) constitute our generic two state system(for fixed momentum, valley index and spin index). It is the linearity in the momentum (or equivalently the wave vector) of the Hamiltonian that justifies the epithet: ‘massless Dirac Fermions’ for the electrons in graphene, as stated earlier. Following mishchenko we now add to Eq. (4) a term associated with a sinusoidal electric field of amplitude $E$ and frequency $\omega$ in the x-direction,

[TABLE]

$e$ being the electron charge. Comparing the above equation with Eq. (1), we have

[TABLE]

II.2 2b. Silicene

Our next example is that of silicene which is made up of a single layer of silicon atoms forming a buckled honeycomb structure (see Fig. 2).

The low energy Hamiltonian about the $K/K^{\prime}$ points has the following form ezawa ,

[TABLE]

where $\tau=\pm 1$ again represent $K/K^{\prime}$ valleys and $\sigma$- and $s$- denote Pauli matrices corresponding to the sub-lattice and spin degrees of freedom, respectively. The first term of Eq. (9) is the linear Dirac term while the second term is due to the spin-orbit interaction (with coupling constant $\lambda_{SO}$) which is responsible for opening up a gap in the spectrum. The third term of the Hamiltonian serves as a gap tuning parameter which arises when an external electric field $E_{z}$ is applied perpendicular to the surface, thus generating different potentials on the A and B lattice sites. Note that the only relevant component of the electronic spin which is to be reckoned with lies along the $z$-axis. However, as $s_{z}$ is a good quantum number, the static Hamiltonian simply splits into two disconnected sectors associated with spin-up and spin-down states (with $s_{z}\rightarrow\pm 1$), thus

[TABLE]

The eigenvalues are $\epsilon_{s\tau}=\pm\sqrt{v_{F}^{2}k^{2}+\Delta_{s\tau}^{2}}$, where $\Delta_{s\tau}=lE_{z}-s\tau\lambda_{SO}$. The corresponding eigenstates are,

[TABLE]

and

[TABLE]

where $\eta_{s\tau}=\sqrt{v_{F}^{2}k^{2}+\Delta_{s\tau}^{2}}-\Delta_{s\tau}$ and $\bar{\eta}_{s\tau}=\sqrt{v_{F}^{2}k^{2}+\Delta_{s\tau}^{2}}+\Delta_{s\tau}$. The perturbation $\mathcal{H}_{I}$ (i.e. the second term in Eq. 7) is taken to be of same form, as the direct coupling of the electric field to the spin $s$ is small.

II.3 2c. $\text{MoS}_{2}$

The third example considered is a monolayer of transition-metal dichalcogenide $\text{MoS}_{2}$ which is composed of an Mo atom caged-in by six neighboring S atoms that are arranged in a prismatic geometry. The top view of monolayer $\text{MoS}_{2}$ has a hexagonal structure with the alternate sites occupied by Mo and S atoms (see Fig. 3). Unlike graphene, the inversion symmetry is explicitly broken in the monolayer of $\text{MoS}_{2}$ making it amenable to valleytronics and valley-dependent optical excitations yao1 . It also supports a strong spin-orbit interaction making it an ideal 2D-system to study spintronics min ; yao2 .

The low-energy unperturbed Hamiltonian about the $K/K^{\prime}$ points is given by yao3 ,

[TABLE]

where, once again, $\tau=\pm 1$ represent $K/K^{\prime}$ valley and $\sigma$ and $s$ are Pauli matrices corresponding to the sub-lattice and spin degrees of freedom, respectively. Note that the last term in Eq. (13) represents the spin-orbit interaction. As before, spin is a good quantum number and the Hamiltonian can be decoupled in to the $s=\pm 1$ sectors. The eigenvalues are $\epsilon_{s\tau}=\pm\sqrt{v_{F}^{2}k^{2}+\Delta_{s\tau}^{2}}+\frac{\lambda\tau s}{2}$, where $\Delta_{s\tau}=\frac{\Delta-\lambda\tau s}{2}$. The corresponding eigenstates are,

[TABLE]

and

[TABLE]

where $|{\epsilon}_{s\tau}|=\sqrt{v_{F}^{2}k^{2}+\Delta_{s\tau}^{2}}$, $\eta_{s\tau}=\sqrt{v_{F}^{2}k^{2}+\Delta_{s\tau}^{2}}-\Delta_{s\tau}$ and $\bar{\eta}_{s\tau}=\sqrt{v_{F}^{2}k^{2}+\Delta_{s\tau}^{2}}+\Delta_{s\tau}$. Once again, our goal will be to study the non-linear current response of $\text{MoS}_{2}$ subject to an oscillating electric field term in the Hamiltonian given by $\tau\sigma_{x}(ev_{F}E/w)\sin(\omega t)$.

III 3. Quantum Master Equation for the Density Operator

With the Hamiltonian (Eq. (1)) at hand, applicable to nonlinear electrical conductivity measurements in graphene, silicene and $\mathrm{MoS_{2}}$ etc., the Von-Neumann-Liouville density operator equation can be written as

[TABLE]

In order to arrive at a master equation for the density operator one has to add to the Hamiltonian $\mathcal{H}_{t}$ in Eq. (1), the coupling to a heat bath. The later, in the present instance, comprises phonon and electron interactions. Upon averaging over the bath in the spirit of Lindblad weiss11 , the averaged density operator can be split as:

[TABLE]

where the subscript S denotes systematic component whereas the subscript D signifies the dissipative part. Needless to say, the heat bath has no influence on the systematic part $\rho_{S}$. We deal with these two terms separately in the following two sub-sections.

III.1 3a. Systematic Dynamics in RWA

Given that our natural basis states in the two-band mode are the conduction and valence states we prefer to diagonalize the unperturbed Hamiltonian in that representation. Further, in order to avoid clumsy notation we simply re-express the conduction and valence states by $\ket{c}$ and $\ket{v}$ respectively. Thus

[TABLE]

wherein all the indices $k$, $s$ and $\tau$ have been subsumed for the sake of brevity. For instance, in graphene,

[TABLE]

We note from the discussion above that though the structure of $\mathcal{H}_{0}$ is different in the three distinct cases of graphene, silicone and $\mathrm{MoS}_{2}$ (because of the additional presence of the spin s) the form of $\mathcal{H}_{I}$ has been assumed to be the same in all the three cases.

As our strategy is to work in the representation in which $\mathcal{H}_{0}$ is diagonal (in the space of conduction and valence states of $\ket{c}$ and $\ket{v}$), it is natural for us to go to the ‘interaction picture’ in which,

[TABLE]

that yields,

[TABLE]

where,

[TABLE]

Rewriting the above in the basis representation of $\ket{c}$ and $\ket{v}$, and ignoring all the rapidly oscillating terms such as $e^{\pm i\omega t}$, $e^{\pm i(\omega+\omega_{cv})t}$ (here, $\omega_{cv}=\epsilon^{0}_{c}-\epsilon^{0}_{v}$), we arrive at an ‘effective’ interaction Hamiltonian:

[TABLE]

Therefore, as in  walls we shall work with an effective Von-Neumann Liouville equation

[TABLE]

It is instructive to spell out what $\mathcal{H}_{eff.}$ is, say for graphene, using the form of $\mathcal{H}_{I}$ given by the second term of Eq.(7) and the fact that $\epsilon^{0}_{c}=kv_{F}$ and $\epsilon^{0}_{v}=-kv_{F}$. We find

[TABLE]

where $\Delta_{k}=\omega-\omega_{cv}$ and $\Omega_{k}=(\frac{ev_{F}E}{2\omega})\sin\chi_{k}$ is the ‘Rabi frequency’.

Using Eqs.(23) and (24) we derive:

[TABLE]

and

[TABLE]

In the above, the ‘population inversion’ $n_{k}$ is defined as

[TABLE]

whereas ‘polarization or inter-band coherence’ is given by

[TABLE]

As we have taken $\mathcal{H}_{I}=(ev_{F}E/\omega)\sigma_{x}$, we obtain the required off-diagonal matrix elements of $\mathcal{H}_{I}$ ($\bra{c}\mathcal{H_{I}}\ket{v}$ and $\bra{v}\mathcal{H_{I}}\ket{c}={\bra{c}\mathcal{H_{I}}\ket{v}}^{*}$), using the form of conduction and valance band states, ($\ket{c}$ and $\ket{v}$ ) corresponding to the case of graphene, silicene and $\mathrm{MoS_{2}}$, given by Eqs. (5), (6), (11), (12),  (14) and  (15). We find, for graphene:

[TABLE]

where $\Omega_{vc}=\Omega_{k}=\frac{ev_{F}E}{2\omega}\sin{\chi_{k}}$, is the so-called ‘Rabi frequency’. Note that these matrix elements are independent of $s$ and $\tau$ i.e. are same for $s_{z}=\pm\frac{1}{2}$ and $\tau=\pm 1$ (K and K’ valleys).

For silicene :

[TABLE]

where $\Omega_{vc}=\alpha_{s\tau}(-\Omega_{k}\,|\epsilon_{s\tau}|+i\,\Omega^{\prime}_{k}\,\Delta_{s\tau})$ and $\Omega_{cv}=\Omega^{*}_{vc}$.

For $\mathrm{MoS_{2}}$:

[TABLE]

Here $\Omega_{vc}=\alpha_{s\tau}(\Omega_{k}\,|\epsilon_{s\tau}|+i\,\Omega^{\prime}_{k}\,\Delta_{s\tau})$ for $\mathrm{MoS_{2}}$ case.

In all three cases

[TABLE]

$\alpha_{s\tau}=\frac{kv_{F}}{|\epsilon_{s\tau}|\sqrt{|{\epsilon_{s\tau}}|^{2}-{\Delta^{2}_{s\tau}}}}$, $\Omega_{k}=\frac{ev_{F}E}{2\omega}\sin{\chi_{k}}$ and $\Omega^{\prime}_{k}=\frac{ev_{F}E}{2\omega}\cos{\chi_{k}}$. Furthermore, $|\epsilon_{s\tau}|$ and $\Delta_{s\tau}$ have been previously defined in the respective sub-sections on silicene and $\mathrm{MoS_{2}}$.

III.2 3b. Dissipative Dynamics

In this sub-section we introduce the relaxation effects which are present in the system because of its interaction with the surrounding heat bath. These effects account for dissipation that characterizes transition of the system from one equilibrium configuration to another. Suffice it to say, the underlying effects are intrinsic to the system, ever present even in the absence of the external perturbation (the electric field, in the present case) – because of spontaneous statistical fluctuations. Thus, our strategy is to examine the dissipative phenomenon in the absence of the external field, by expanding the unperturbed Hamiltonian to a form that incorporates the system-plus-bath coupling. Once we have a master equation for the density operator at hand (with dissipative or relaxational terms present) we will add that to the contribution arising from the systematic component, as dealt with in section 3a.

The system-plus-bath Hamiltonian will have to account for two distinct physical interactions between the Dirac electron and the surrounding phonons and other electrons. The phonons that describe lattice relaxations can cause transitions between the conduction and valence bands thus maintaining the system in thermal equilibrium governed by the Fermi-Dirac distribution. Thus these transitions involve ‘depopulation’ - a terminology well known in vibrational relaxation in molecular spectra dattagupta – of the valence and conduction levels. Evidently, the coupling to the phonons will have to occur via a term that is purely off-diagonal in the basis states of conduction and valence bands, thus yielding:

[TABLE]

Here, $g_{q}$ is a coupling constant and the last term denotes a ‘free’ phonon Hamiltonian, the $b$’s and $b^{\dagger}$‘s being the annihilation and creation operators for phonons.

On the other hand, there exists a much stronger interaction with the other electrons which however do not trigger any transition between the two bands but simply cause ‘de-phasing’ thereby leading to what is known as ‘dissipation-less de-coherence’  SDSP . The corresponding interaction can be represented by

[TABLE]

Once again, $G_{p}$ parametrizes the coupling with the electrons. Evidently, the coupling term now commutes with the unperturbed Hamiltonian and hence, does not influence any energy exchange between the system and the surrounding electronic bath. The last term contains the free electron interaction in terms of electron creation and annihilation operators. The two level system in coupling with a fermionic bath has been treated in detail in chang . As shown there, the bath(fermionic) can be ‘bosonized’ as long as we are only interested in the electron-hole excitations near the fermi surface, for what is called ‘Ohmic dissipation’. We will adopt that interpretation here, since we are only considering low energy Hamiltonian and excitations near the K/K’ valley points.

The full system-plus-bath Hamiltonian can then be written as

[TABLE]

The Hamiltonian given above in Eq. (36) can form the basis of a detailed study of quantum dissipative phenomena for the two-band materials at hand, following the methods reviewed in weiss11 . While relegating that investigation to the future we study here a simpler model in which the phonon and electronic bath operators are replaced by effective fields and follow standard treatments of master equations (chapter 1) of  SDSP . We should note however that because of the harmonic nature of the phonon and electron baths and the assumed linear coupling, as is prevalent also in Feynman-Vernon, Caldeira-Leggett treatments (See weiss11 ), only Gaussian fluctuations are relevant. Thus, Cumulant Expansion methods, as detailed in Appendix I.A.I of  SDSP , are the most suitable because all cumulants beyond the second order vanish.

The simplified system-bath Hamiltonian can then be rewritten as

[TABLE]

In the above, g and G are the phonon and electron coupling constants (in which their momentum-dependence has been ignored), B and C are the classical version of phonon and electron ‘displacement operators’, and $\mathcal{H}_{ph}$ and $\mathcal{H}_{el}$ are the free phonon and electron Hamiltonians, respectively.

Following the results given in I.7 of SDSP , the different components of the dissipative part of the density operator, in the interaction picture of the unperturbed Hamiltonian $\mathcal{H}_{0}$, can be shown to satisfy

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The Eq. (39) is consistent with the conservation of probability: $\bra{c}\tilde{\rho}_{D}\ket{c}+\bra{v}\tilde{\rho}_{D}\ket{v}=1$. In the following we treat $\lambda$ and $\gamma$ as lumped parameters which have been listed by couplings as $\lambda\approx 10^{12}$ Hz while $\lambda+\gamma=\Upsilon\approx 10^{14}$ Hz for graphene. Here, we identify that, $\Upsilon=\lambda+\gamma$ is the total rate due to the combined phonon and electron bath.

It turns out that SDSP , for phonon bath,

[TABLE]

where plus/minus sign in the exponential is for valence/conduction band index and $<\textbf{B}(\tau)\textbf{B}(0)>$ is an auto-correlation of the phonon bath operators defined in the interaction picture as

[TABLE]

Note that the relation between the Fermi distribution $p_{c}$ and $p_{v}$(detailed balance) follows from the spectral properties of the correlation function. The equilibrium population is then, $n^{eq}_{k}=p_{c}-p_{v}$. For electron bath, the relaxation rate $\gamma$, on the other hand, is given by

[TABLE]

where $<\textbf{C}(\tau)\textbf{C}(0)>$ is an auto-correlation of the electron bath operators defined in the interaction picture as,

[TABLE]

The coupling constants $g$ and $G$ (on which the rates $\lambda$ and $\gamma$ depend) can be calculated from microscopic considerations of electron-phonon scattering and electron-electron scattering for different materials e.g. graphene as discussed in grapheneEP1 ; grapheneEP2 ; grapheneEP3 ; grapheneEP4 and grapheneEE0 ; grapheneEE1 ; grapheneEE2 ; grapheneEE3 , silicene SiliceneEP1 ; SiliceneEP2 and $\mathrm{MoS}_{2}$ MoS2EP1 ; MoS2EP2 .

In the above analysis we have obtained the dissipative evolution part and the systematic evolution part of the Eq. (17) as given by the Eq. (26) and Eq. (27), where the matrix elements $\bra{v}\mathcal{H_{I}}\ket{c}=2i\Omega_{vc}$ and $\bra{c}\mathcal{H_{I}}\ket{v}=-2i\Omega_{cv}$ and expressions for $\Omega_{vc}$ (and $\Omega_{cv}=\Omega^{*}_{vc}$) have been stated in the previous section corresponding to our three cases. We are now ready to write down the equation (17) in component form (in the interaction picture), thereby yielding

[TABLE]

[TABLE]

Let us define a following transformation for the elements of the density matrix,

[TABLE]

Then the above equations becomes

[TABLE]

[TABLE]

The steady state solutions for these equations are given by

[TABLE]

[TABLE]

IV 4. Nonlinear Optical Conductivity

For calculating the conductivity we need the cycle-averaged power absorbed, defined by

[TABLE]

The general expression for the d-dimensional system is given by

[TABLE]

where $g_{s}\,(g_{v})$ is the spin (valley) degeneracy factor(in our case both are 2) and momentum dependent component of particle current density

[TABLE]

and $x$-direction component of this,

[TABLE]

A general expression for the $\vec{J}_{\vec{k}x}(t)$ can be calculated, given the Hamiltonian and the valence and conduction eigenstates for any 2-band system, which in our case

[TABLE]

In order to convert this expression in terms of transformed density matrix’s elements ($\Pi^{\prime}_{k},\,\Pi^{\prime*}_{k},\,n^{\prime}_{k}$) we use Eqs. (20) and  (48) to get $\Pi_{k}(t)=\Pi^{\prime}_{k}\,e^{i\omega t}$, $\Pi^{*}_{k}(t)=\Pi^{\prime*}_{k}\,e^{-i\omega t}$ and $n_{k}(t)=n^{\prime}_{k}$. When we substitute Eq. (57) in to the Eq. (53) and integrate over a complete cycle of the applied field, the first term of Eq. (57) yields zero while the second term gives non-zero contribution. Therefore we will consider from Eq. (57) only the second term, i.e.,

[TABLE]

Using the solution for $\Pi^{\prime}_{k}$ from Eq. (52) in Eq. (58) yields,

[TABLE]

where $\Omega_{cv}$ and $\Omega^{*}_{vc}=\Omega_{cv}$ have been given in the previous sections corresponding to the three cases, while $n^{\prime}_{k}(t)$ is given by Eq. (51) and $\Delta_{k}=\omega-\omega_{cv}$. It should be noticed that only the term proportional to $\cos{\omega t}$ (out of phase with the applied term of Eq. (7) ) is responsible for dissipation that contributes to the conductivity and the other term proportional to $\sin{\omega t}$ will be dropped from further consideration. Therefore,

[TABLE]

The general expression for the longitudinal non-linear optical conductivity is given by,

[TABLE]

We next write down the expressions for $\textbf{J}_{\vec{k}x}(t)$ for the three cases of interest. For graphene, we have $\Omega_{vc}=\Omega_{k}=\frac{ev_{F}E_{0}\sin{\chi_{k}}}{2\omega}=\Omega_{cv}$, $\Delta_{k}=\omega-\omega_{cv}$ and $\omega_{cv}=2kv_{F}$, plugging these into the Eq. (60) one gets,

[TABLE]

For silicene, we have $\Omega_{cv}=\alpha_{s\tau}(-\Omega_{k}|\epsilon_{s\tau}|+i\Omega^{\prime}_{k}\Delta_{s\tau})$, $\Omega^{*}_{vc}=\Omega_{cv}$, $\Omega_{cv}=2|\epsilon_{s\tau}|=2\sqrt{k^{2}v^{2}_{F}+\Delta^{2}_{s\tau}}$ , $\Omega_{k}=\frac{ev_{F}E_{0}\sin{\chi_{k}}}{2\omega}$, $\Omega^{\prime}_{k}=\frac{ev_{F}E_{0}\cos{\chi_{k}}}{2\omega}$ and $\alpha_{s\tau}=\frac{kv_{F}}{|\epsilon_{s\tau}|\sqrt{|\epsilon_{s\tau}|^{2}-\Delta^{2}_{s\tau}}}$. putting these in the Eq. (60) gives us,

[TABLE]

Finally, for $\mathrm{MoS_{2}}$ we have $\Omega_{cv}=\alpha_{s\tau}(\Omega_{k}|\epsilon_{s\tau}|+i\Omega^{\prime}_{k}\Delta_{s\tau})$ and $\Omega^{*}_{vc}=\Omega_{cv}$. Plugging it in Eq. (60) yields a term identical to Eq. (63). Thus the longitudinal optical conductivity in all the three cases can be calculated by plugging the corresponding $\textbf{J}_{\vec{k}x}(t)$ in Eq. (61) and integrating over the allowed range of momentums (near the vicinity of K/K’ points) for which the Dirac Hamiltonian form is valid.

V 5. Results and Analysis

Employing $\Omega_{vc}=\frac{\tau ev_{F}E_{0}}{\omega}\sigma^{x}_{vc}$, then Eq. (60) can be cast in the following form

[TABLE]

where,

[TABLE]

And $n^{\prime}_{k}$ is given by the Eq. (51), which can also be written in the following form after defining dimensionless parameter $\xi=\frac{ev_{F}E_{0}}{\omega\sqrt{\lambda\Upsilon}}$ characterizing linear regime ($\xi\ll 1$) and non-linear regime($\xi\geq 1$),

[TABLE]

In our cases, spin and valley degeneracy factors $g_{s},g_{v}=2$ and putting the expression for $\textbf{J}_{\vec{k}x}(t)$ in Eq. (61) from Eq. (64), we get,

[TABLE]

where $|\sigma^{x}_{vc}|^{2}$ and $n^{\prime}_{k}$ are respectively given by Eq. (65) and Eq. (66). Note that the Eq. (67) is a simplified version of Eq. (22) of  agarwal . The latter authors in agarwal have made further discussion of the longitudinal optical conductivity in detail for different cases-linear and non-linear in different regimes-clean and dirty. We can derive results using Eq. (67) for silicene and $\mathrm{MoS}_{2}$ which however appear to match those of Agarwal et al. for ‘gapped’ graphene. We do not repeat that analysis here.

VI 6. Concluding Remarks

We have calculated the non-linear optical conductivity for gap-less(graphene) and gapped (silicene and $\mathrm{MoS_{2}}$) two band systems. We follow a density matrix approach to study the time evolution of these systems. The treatment presented in section 3 is the core of the paper that offers a relaxation model for describing dissipative response of the 2-D materials at hand, to an externally applied frequency dependent field. After obtaining the time evolution equations for different components of the density matrix including damping terms, we solve for the case of steady state and get analytical expressions for population inversion and inter-band coherence given by Eq. (51) and Eq. (52). A comment is in order on clarifying the difference in our treatment from that of mishchenko . The phenomenological equations(10) of  mishchenko incorporate a damping parameter $\Gamma_{p}$ in however the wave function equations. The problem with this approach is that the total probability, in the presence of dissipation, is not conserved. Our method circumvents this issue by looking at a fully quantum master equation for density operator, besides taking care of two time-scales (inverse relaxation rates) associated with electron-phonon and electron-electron interactions.

As discussed in agarwal , we may categorize the four regimes depending on two parameters $\xi$ (defines non-linearity) and ${\Upsilon}/{\omega}$ (characterizes clean or dirty system) as follow: (a) linear clean limit (in which $\xi\ll 1$ and $\Upsilon/\omega\ll 1$ ), (b) linear dirty limit (in which $\xi\ll 1$ and $\Upsilon/\omega\geq 1$ ), (c) non-linear clean limit (in which $\xi\geq 1$) and $\Upsilon/\omega\ll 1$ and (d) non-linear dirty limit (in which $\xi\geq 1$ and $\Upsilon/\omega\geq 1$). The parameter $\xi$ first appears naturally in the steady state solutions for population inversion and inter-band coherence and then in the conductivity expressions in various cases. If we assume $\xi\ll 1$, the results match the linear response results, that’s why $\xi$ specifies the strength of optical non-linearity in the system. Furthermore, for $\Upsilon/\omega\ll 1$, the term $\frac{\Upsilon}{(\omega_{k}-\omega)^{2}+\Upsilon^{2}}\rightarrow\delta(\omega_{k}-\omega)$ in the momentum dependent conductivity expression (Eq. (64)), which means this non-zero contribution comes only from the frequency $\omega=\omega_{k}$, that’s why $\Upsilon/\omega\ll 1$ is called ‘clean limit’. When $\Upsilon/\omega\geq 1$, then nearby frequencies to $\omega=\omega_{k}$ can also contribute to the optical conductivity and hence is called ‘dirty limit’.

Finally, with the steady state solutions for density matrix components in hand, we can calculate the optical conductivity in these four regimes and obtain analytical solutions for most of the cases in $T\rightarrow 0$ limit. Some specific cases need numerical evaluation for momentum dependent conductivity integrals. One may note that the linear response theory (based on the Kubo formula) for conductivity is only valid in the high frequency limit ($\Upsilon/\omega\ll 1$).

In conclusion, our ability to derive the known results for non-linear optical conductivity in gapless and gapped graphene, based on a simplified treatment of the Hamiltonian Eq. (36), may spur further and more extensive study of what is generally known as Spin-Boson Hamiltonian (Eq. (36)) in the context of quantum dissipative systems weiss11 .

VII 7. Acknowledgement

S.D. is grateful to the Indian National Science Academy for support through their Senior Sciencist scheme and to IISER, Bhopal for kind hospitality. We thank Suhas Gangadharaiah for useful discussions in the early stages of the work.