Dynamical band gap tuning in Weyl semi-metals by intense elliptically polarized normal illumination and its application to $8-Pmmn$ borophene

TL;DR

This paper investigates how intense elliptically polarized light can dynamically modulate the band gap in Weyl semimetals, specifically applying the findings to $8-Pmmn$ borophene, using a non-perturbative Floquet approach.

Contribution

It introduces a non-perturbative method to predict dynamical gap formation in Weyl semimetals under elliptically polarized light, applicable to various materials.

Findings

Derived the quasienergy spectrum using Mathieu equations.

Identified conditions for dynamical gap formation and stability.

Estimated electromagnetic parameters for gap tuning in borophene.

Abstract

The Dynamical-gap formation in Weyl semimetals modulated by intense elliptically polarized light is addressed through the solution of the time-dependent Schr\"odinger equation for the Weyl Hamiltonian via the Floquet theorem. The time-dependent wave functions and the quasi-energy spectrum of the two-dimensional Weyl Hamiltonian under normal incidence of elliptically polarized electromagnetic waves are obtained using a non-perturbative approach. In it, the Weyl equation is reduced to an ordinary second-order differential Mathieu equation. It is shown that the stability conditions of the Mathieu functions are directly inherited by the wave function resulting in a quasiparticle spectrum consisting of bands and gaps determined by dynamical diffraction and resonance conditions between the electron and the electromagnetic wave. Estimations of the electromagnetic field intensity and frequency, as well as the magnitude of the generated gap are obtained for the $8-Pmmn$ phase of borophene. We provide with a simple method that enables to predict the formation of dynamical-gaps of unstable wave functions and their magnitudes. This method can readily be adapted to other Weyl semimetals.

Equations

$$\hat{H}=\mu \left(v_t{\hat{P}}_y{\hat{\sigma }}_0+v_x{\hat{P}}_x{\hat{\sigma }}_x+v_y{\hat{P}}_y{\hat{\sigma }}_y\right),$$

$$E_{\eta ,k}^{\mu }=\mu \hslash v_t{k}_y+\eta ϵ,$$

$$ϵ=\hslash \sqrt{v_x^2{k}_x^2+v_y^2{k}_y^2}.$$

$$\psi_{\eta,\,\bf{k}}^{\mu}(\mathbf{r})=\mu\frac{\exp(i\bf{k\cdot r})}{\sqrt{2}}\left[\begin{array}[]{lcc}1\\ \eta\exp(i\Theta)\end{array}\right],$$

$$\hat{H}=\left(\begin{array}[]{lcc}v_{t}\hat{\Pi}_{y}&v_{x}\hat{\Pi}_{x}-iv_{y}\hat{\Pi}_{y}\\ v_{x}\hat{\Pi}_{x}+iv_{y}\hat{\Pi}_{y}&v_{t}\hat{\Pi}_{y}\end{array}\right),$$

$$\hat{H}(\mathbf{r},t)\mathbf{\Psi }(\mathbf{r},t)=i\hslash \frac{\partial }{\partial t}\mathbf{\Psi }(\mathbf{r},t),$$

$$\mathbf{\Psi }(\mathbf{r},t)=\exp \left(i\mathbf{k}\cdot \mathbf{r}\right)\mathbf{\Phi }(t),$$

$$\text{mathbbm}H(t)\mathbf{\Phi }(t)=i\hslash \frac{d}{dt}\mathbf{\Phi }(t),$$

$$\mathbf{\Phi }(t)=\exp \left[-\frac{i}{\hslash }{\int }^{t}ds{\alpha }_{\mathbf{k}}(s)\right]\mathbf{\chi }(t),$$

$$\frac{d^2}{d{t}^2}\mathbf{\chi }(t)+\text{mathbbm}F(t)\mathbf{\chi }(t)=0,$$

$$\text{mathbbm}F(t)=-\frac{i}{\hslash }\hat{\mathbf{\sigma }}\cdot \frac{d\mathbf{S}}{dt}+\frac{1}{{\hslash }^2}{\left[\widetilde{\text{mathbbm}H}(t)\right]}^2,$$

$$\mathbf{A}=\frac{1}{\Omega }\left(E_x\cos (\Omega t),E_y\sin (\Omega t)\right),$$

$$\varphi =\Omega t,$$

$${\mathbf{\chi }}^{\prime \prime }(\varphi )+\text{mathbbm}F(\varphi )\mathbf{\chi }(\varphi )=0,$$

$$\begin{gathered} \text{mathbbm}F(\varphi )=i\left(\frac{{\zeta }_x}{\hslash \Omega }{\hat{\sigma }}_x\sin \varphi -\frac{{\zeta }_y}{\hslash \Omega }{\hat{\sigma }}_y\cos \varphi \right) \\ +{\left(\frac{1}{\hslash \Omega }\right)}^2[{ϵ}^2-2{\mathbf{\kappa }}_{\mathbf{v}}\cdot \mathbf{S}+\frac{1}{2}\left({\zeta }_x^2+{\zeta }_y^2\right) \\ +\frac{1}{2}\left({\zeta }_x^2-{\zeta }_y^2\right)\cos (2\varphi )], \end{gathered}$$

$$\mathbf{S}$$

$${\zeta }_x$$

$${\zeta }_y$$

$${\mathbf{\chi }}^{\prime \prime }(\varphi )+[a-2q\cos (2\varphi )]\mathbf{\chi }(\varphi )=0.$$

$$q$$

$$a$$

$$\sqrt{a}=\frac{{\Omega }_0}{\Omega },$$

$$\mathbf{\chi }(\varphi )=\mathcal{C}(a,q,\varphi )=\exp [-ir(a,q)\varphi ]{\mathbf{f}}_C(\varphi ),$$

$$\mathbf{\chi }(\varphi )=\mathcal{S}(a,q,\varphi )=\exp [-ir(b,q)\varphi ]{\mathbf{f}}_S(\varphi ),$$

$$\frac{{\Omega }_0}{\Omega }=1,2,3,...,$$

$$\boldsymbol{\Psi}(\bm{r},t)=\mathcal{N}\exp\left[i\left(\boldsymbol{k}\cdot\boldsymbol{r}-v_{t}k_{y}t-\frac{v_{t}\zeta_{y}}{v_{y}\hbar\Omega}\cos(\Omega t)\right)\right]\\ \times\left[\,\mathcal{C}(a,q,\Omega t)\pm i\eta\mathcal{S}(a,q,\Omega t)\,\right]\left(\begin{array}[]{lcc}1\\ \eta\exp(i\Theta)\end{array}\right),$$

$$a\ge {\left(\frac{{\zeta }_y}{\hslash \Omega }\right)}^2-2q.$$

$$a=({\zeta }_y/\hslash \Omega {)}^2-2q,$$

$${\left(\frac{v_x}{v_F}\right)}^2{\tilde{k}}_x^2+{\left(\frac{v_y}{v_F}\right)}^2{\tilde{k}}_y^2=a+2q-{\left(\frac{{\zeta }_y}{\hslash \Omega }\right)}^2,$$

$$\Delta E(q)/2=\hslash \Omega \sqrt{c_r(q)-\left[{\left(\frac{{\zeta }_y}{\hslash \Omega }\right)}^2-2q\right]},$$

$$\mathbf{\Psi }(\mathbf{r},t)=\mathcal{N}\exp \left(-i\frac{\mathcal{E}t}{\hslash }\right)\mathbf{f}(\Omega t).$$

$$\mathcal{C}(a,q,\Omega t)-i\eta \mathcal{S}(a,q,\Omega t)=\exp [-i\eta r(a,q)\Omega t]f(\Omega t),$$

Full text

Keywords

Weyl semi-metals; dynamical gap; Mathieu equation

Suggested keywords

I Introduction

Recently, the so-called Dirac and Weyl materials have received considerable attention due to their possible implementation into next-generation electronic devices Novoselov et al. (2012); Peng et al. (2016); Ferrari et al. (2015); Naumis et al. (2017). The main physical properties of these materials were observed for the first time in graphene, an allotrope of carbon consisting of a monolayer of atoms in a honeycomb lattice with an electron linear dispersion near the Dirac points. As a result of this characteristics, the charge carriers in graphene behave like massless Dirac fermions Geim (2009); Neto et al. (2009); Oliva-Leyva and Naumis (2014); Naumis and Roman-Taboada (2014); Oliva-Leyva and Naumis (2015a, b, c); Naumis et al. (2017); Carrillo-Bastos and Naumis (2018). Thereafter, a wide variety of two-dimensional materials with similar properties has been discovered Wehling et al. (2014). Examples of these are: silicene Jose and Datta (2013); Liu et al. (2013), germanene Behera and Mukhopadhyay (2011); Zhang et al. (2016a), stanene Zhu et al. (2015); Chen et al. (2016) and artificial graphene Soltan-Panahi et al. (2011); Gomes et al. (2012). Borophen, a two dimensional allotrope of boron, also falls in this category. The chemical similarity between boron and carbon atoms has triggered the search for stable two-dimensional boron structures and synthesis techniques to produce themZhang et al. (2017). Since their theoretical prediction Boustani (1997), many different allotropes of borophene have been experimentally confirmed Mannix et al. (2015). Among its many different phases, the orthorhombic $8-Pmmn$ is one of the most energetically stable structures Peng et al. (2016), having a ground state energy lower than its analogs Verma et al. (2017). Borophene, in contrast with graphene, shows a highly anisotropic crystalline structure, which causes high optical anisotropy and transparency Peng et al. (2016); Verma et al. (2017); Champo and Naumis (2019). It is thus a strong candidate for flexible electronics, display technologies and in the design of smart windows where minimal photon absorption and reflection are required Zhang et al. (2016b); Peng et al. (2016); Zhang et al. (2017).

Despite the many useful and fascinating properties of graphene, borophene and Weyl materials in general, their lack of an electronic band gap has stimulated the search for either other two-dimensional materials with semiconducting properties or techniques to induce them artificially. Among other proposals to circumvent this problem, one of the most promising ideas is generating a light induced dynamical-gap. As the electromagnetic field is a periodical function of time, this technique has been termed Floquet gap engineering. High intensity electromagnetic waves interacting with graphene have been studied using perturbative approaches López-Rodríguez and Naumis (2010); Higuchi et al. (2017). However, it has been shown that light induces a renormalization of the electronic spectrum of Dirac materials not captured by simple perturbation techniques López-Rodríguez and Naumis (2008, 2010); Kibis (2010); Kristinsson et al. (2016a); Oliva-Leyva and Naumis (2016); Kristinsson et al. (2016b); Kibis et al. (2017). In this regard, borophene brigns interesting possibilities to study the light-matter interaction in Weyl semimetals due to its asymmetric spectrum. As graphene, borophene has a honeycomb lattice with two nonequivalent sublattices. However, its peculiar structure give rise to a tilted anisotropic cone in the vicinity of the Dirac points Islam and Jayannavar (2017); Verma et al. (2017); Champo and Naumis (2019), as opposed to graphene whose spectrum is completely isotropic in $K$ space.

Recently, the formation of energy gaps in borophene subject high-intensity linearly polarized light was studied beyond the perturbative approximation Champo and Naumis (2019). It was found that borophene, when interacting with light, acquires a complex band structure from the stability conditions of the solutions of the Mathieu differential equation. Among other effects, the interaction with light produces a gap in the vicinity of borophene’s Dirac point. The effects of an intense circularly polarized electromagnetic field have, nevertheless, not been discussed for borophene yet.

In this paper, we address the general problem of a particle that obeys the Weyl Hamiltonian subject to an intense elliptically polarized electromagnetic field. As it is schematized in Fig. 1, our solutions can be applied to the particular case of electrons in borophene under a strong elliptically polarized field. We report the wave functions, the quasi-energy spectrum, and the magnitude of the dynamic gap opening. The case of linearly polarized light, addressed by us previously Champo and Naumis (2019), is proven to be fundamentally different from the elliptically polarized one studied in this work. Our analysis mainly focuses on the stability and instability of the time-dependent wave functions. The results presented here display an interesting interplay between the tilted anisotropy and the relative orientation of the light-polarization ellipse. Moreover, we show that the gaps may be tuned by changing the orientation of the elliptical polarization profile of light.

The paper is organized as follows. In Sec. II we introduce the low-energy effective two-dimensional Weyl Hamiltonian under an arbitrary electromagnetic field. Subsequently, in Sec. III, we determine the time-dependent wavefunction of electrons in borophene subject to an elliptically polarized electromagnetic field. In this same section we analyze the stability of the solutions inherited from Mathieu functions in the strong electromagnetic field or long wavelength regimes. We workout the time-dependent wave functions and the solutions’ stability chart. To get an insight into the gap structure, in Sec. III.2, the stability and instability regions are projected onto the tilted Dirac cones of the free Weyl electrons. In Sec. III.2, we extract the quasi-energy spectrum from the time-dependent wave function and prove that it consistently shows a similar gap structure to that of the projected chart. Finally, we summarize and conclude in Sec. IV.

II Weyl electrons subject to electromagnetic

fields

II.1 The Weyl Hamiltonian

The single-particle low-energy effective Weyl Hamiltonian is given by Zabolotskiy and Lozovik (2016); Islam and Jayannavar (2017); Verma et al. (2017); Champo and Naumis (2019),

[TABLE]

where $\hat{P}_{x}$ and $\hat{P}_{y}$ are the momentum operators, $\hat{\sigma}_{i}$ are the Pauli matrices, and $\hat{\sigma}_{0}$ is the $2\times 2$ identity matrix. For $8-Pmmn$ borophene, the three velocities in the Weyl Hamiltonian (1) are $v_{x}=0.86v_{F}$, $v_{y}=0.69v_{F}$ and $v_{t}=0.32v_{F}$ in units of the Fermi velocity $v_{F}=10^{6}\,m/s$ Zabolotskiy and Lozovik (2016). The two Dirac points $\mathbf{k}=\pm\mathbf{k}_{D}$ are given by the valley index $\mu=\pm 1$. The first term in Eq. (1) gives rise to the tilting of the Dirac cones and the last ones correspond to the kinetic energy.

The previous Hamiltonian results in the energy dispersion relation Islam and Jayannavar (2017)

[TABLE]

where

[TABLE]

The corresponding free Weyl electron wave function is,

[TABLE]

where $\eta=\pm 1$ is the band index, $\Theta=\tan^{-1}(v_{y}k_{y}/v_{x}k_{x})$ and the two-dimensional momentum vector is given by $\mathbf{k}=(k_{x},k_{y})$.

II.2 The Weyl Hamiltonian in the presence of

an electromagnetic wave

Now we consider a charge carrier described by the Weyl Hamiltonian subject to an electromagnetic wave that propagates along a direction perpendicular to the surface of the crystal. From Eq. (1) and using the minimal coupling we obtain,

[TABLE]

where $\hat{\boldsymbol{\Pi}}=\hat{\boldsymbol{P}}-e{\boldsymbol{A}}$, with $\boldsymbol{A}=(A_{x},A_{y})$ being the vector potential of the incident electromagnetic wave. Calculations are considerably simplified by choosing a gauge in which the vector potential is only a function of time. The Schrödinger equation for charge carriers in a Weyl semimetal is thus given by

[TABLE]

where, in the two dimensional spinor $\boldsymbol{\Psi}(\mathbf{r},t)=\left(\Psi_{A}(\mathbf{r},t),\Psi_{B}(\mathbf{r},t)\right)^{\top}$, $A$ and $B$ label the two sublattices.

To deduce the explicit form of the wave function ${\bf{\Psi}}(\mathbf{r},t)$ from Eq. (6) we make the following ansatz

[TABLE]

where $\mathbf{\Phi}(t)=\left(\Phi_{A}(t),\Phi_{B}(t)\right)^{\top}$. Substituting (7) reduces Eq. (6) into

[TABLE]

where the matrix $\mathbbm{H}(t)$ is defined in Appendix A.

The diagonal terms of $\mathbbm{H}(t)$ can be lifted by explicitly adding a time-dependent phase to the wave function

[TABLE]

with $\alpha_{\mathbf{k}}(t)=\hbar v_{t}k_{y}-ev_{t}A_{y}(t)$ and $\bm{\chi}(t)=(\chi_{A}(t),\chi_{B}(t))^{\top}$. Following the procedure shown in the Appendix A, Eq. (8) can be recast in the form of a second-order ordinary differential equation as

[TABLE]

where the function $\mathbbm{F}(t)$ is defined by

[TABLE]

with $[\widetilde{\mathbbm{H}}(t)]^{2}=|\boldsymbol{\kappa}_{\boldsymbol{v}}-\boldsymbol{S}|^{2}$, $\boldsymbol{\kappa}_{\boldsymbol{v}}=\hbar(v_{x}k_{x},v_{y}k_{y})$ and $\boldsymbol{S}=e(v_{x}A_{x},v_{y}A_{y})$. In the last expression, the vector $\boldsymbol{\kappa}_{\boldsymbol{v}}$ is the directional energy flux of the electrons, and the components of $\boldsymbol{S}$ represent the work done by the electromagnetic wave along the $x$ and $y$ directions.

III Elliptically polarized waves

Let us now study the case of an elliptically polarized electromagnetic wave characterized by the vector potential

[TABLE]

where $E_{x}$ and $E_{y}$ are constants and $\Omega$ is the frequency of the electromagnetic wave. The vector potential (12) corresponds to the electric field $\bm{E}=-\partial\boldsymbol{A}/\partial t=(E_{x}\sin(\Omega t),-E_{y}\cos(\Omega t))$.

Rewriting Eq. (10) in terms of the phase

[TABLE]

yields the Hill equation Magnus and Winkler (2013)

[TABLE]

where $\mathbbm{F}(\phi)$ is

[TABLE]

and

[TABLE]

The unitless parameter $\epsilon/\hbar\Omega$ is the ratio of the electron energy to the photon energy. Similarly, the parameter $\zeta_{x}/\hbar\Omega$ ($\zeta_{y}/\hbar\Omega$) is the ratio of the work done by the electromagnetic wave along the $x$ ($y$) direction to the photon energy.

The determination of the stability regions of the differential Eq. (14) is quite challenging mainly due to the imaginary part in the first term of the right-hand side of Eq. (15). While the real part gives rise to the Whittaker-Hill equationUrwin and Arscott (1970), the imaginary term yields a Mathieu-like equation with complex characteristic values, rarely discussed in literature Ziener et al. (2012). Fortunately, in the intense electric field or long wavelength regimes the imaginary part is negligible. Other limits are treatable by perturbation theory López-Rodríguez and Naumis (2010); Higuchi et al. (2017).

Here we focus on the intense electric field regime. We thus assume that $\zeta_{i}/\hbar\Omega\gg 1$ with $i=x,\,y$, which is equivalent to $ecE_{x}/\hbar\Omega^{2}\gg 349$ and $ecE_{y}/\hbar\Omega^{2}\gg 435$. This corresponds to electric fields $E_{x}\gg 1.91\,V/m$ and $E_{y}\gg 2.39\,V/m$. Thereby, we can neglect the linear terms of $\zeta_{i}/\hbar\Omega$ in Eq. (14) that yield the imaginary terms. The obtained expression, best-known for describing the dynamics of the parametric pendulumAldrovandi and Ferreira (1980); Baker and Blackburn (2005), is the Mathieu differential equation

[TABLE]

The purely real parameters $q$ and $a$ are given by

[TABLE]

The characteristic value of the Mehtieu equation

[TABLE]

with $\Omega_{0}=\sqrt{2\epsilon^{2}+\zeta^{2}_{x}+\zeta^{2}_{y}}/\sqrt{2}\hbar$ is the ratio between the fundamental frequency $\Omega_{0}$ and the frequency of the electromagnetic wave $\Omega$. The characteristic value $a$ parametrizes the family of ellipses in the $k_{x}-k_{y}$ plane that are characterized by the eccentricity $[1-(v_{y}^{2}/v_{x}^{2})]^{1/2}$. Stated differently, each value of the parameter $a$ corresponds to a particular elliptical section of the Dirac cone. However, not all the ordered pairs in the $q-a$ plane produce stable solutions of the Mathieu equation. Consequently, in the presence of an intense electromagnetic radiation not all the elliptical sections of the Dirac cones correspond to stable solutions. In fact, the interaction with light induces elliptical sections of the Dirac cone that alternate between forbidden (unstable) and allowed (stable) solutions. As can be seen in Fig. 2, the stability chart in the $q-a$ plane consists of tongue-like stable regions (light blue) that neighbor with unstable regions (white). The Mathieu equation might have either even or odd stable solutions. Even stable solutions of (19) have the form

[TABLE]

where $\mathcal{C}(a,q,\phi)$ is the even Mathieu function, $\boldsymbol{f}_{C}(\phi)$ is an even function with period $\pi$ and $a=a_{r}(q)$ is the Mathieu even characteristic value. Conversely, odd stable solutions have the form

[TABLE]

where $\mathcal{S}(a,q,\phi)$ is the odd Mathieu function, $\boldsymbol{f}_{S}(\phi)$ is an odd function with period $\pi$ and $b=b_{r}(q)$ is the Mathieu odd characteristic value. When $r$ is a non-integer rational number, inside the stable regions, the even and odd characteristic values are identical, namely $a_{r}(q)=b_{r}(q)$. The rational function $r(a,q)$ depends on the Mathieu characteristic value $a$ and the parameter $q$. On the boundaries between the stable and unstable regions (solid and dashed blue lines Fig. 2) $r$ takes an integer value and inside the stable regions $r$ is a non-integer rational number. Thus, inside the stability regions the even and odd Mathieu functions have the same characteristic value. For the particular situation in which $q=0$, (19) reduces to the differential equation of an harmonic oscillator whose solutions are $\cos(\sqrt{a}\phi)$ and $\sin(\sqrt{a}\phi)$ McLachlan (1951). Evidently, in this case $r=\sqrt{a}$. Moreover, in the special case where $r=\sqrt{a}\in\mathbbm{Z}$, resonant states are generated for which

[TABLE]

and therefore when $q\to 0$ two contiguous stability zones are connected.

III.1 Wave function and stability spectrum

The general solution of Eq. (19) is the superposition of the even and odd Mathieu functions $\mathcal{C}(a,q,\phi)$ and $\mathcal{S}(a,q,\phi)$. The wave function is then given by

[TABLE]

where $\mathcal{N}$ is a normalization constant, $\Theta=\tan^{-1}\left(v_{y}k_{y}/v_{x}k_{x}\right)$, and $\eta=\pm 1$ denotes the conduction and valence bands, respectively. The wave function (26) reduces to the free-particle wave function (4) when the electric field vanishes.

Since the time-dependent wave function is expressed in terms of the Mathieu functions, its stability is governed by the stability chart in 2 that we discussed previously. Indeed, the structure of the dynamical gaps of Weyl electrons, generated in the presence of an intense electromagnetic radiation, is inherited from the properties of the characteristic values of the Mathieu functions.

The chart can be divided into two key regions according to the shape of the electromagnetic wave: $q<0$ ($\zeta_{x}>\zeta_{y}$) and $q>0$ ($\zeta_{y}>\zeta_{x}$) and $q=0$ ($\zeta_{x}=\zeta_{x}$). If $\zeta_{x}>\zeta_{y}$ ($E_{x}v_{x}>E_{y}v_{y}$) the work done by the electromagnetic wave on the electrons is higher along the $x$ axis. Conversely, if $\zeta_{x}<\zeta_{y}$ ($E_{x}v_{x}<E_{y}v_{y}$) the work is higher along $y$. Finally, if $\zeta_{x}=\zeta_{y}$ ($E_{x}v_{x}=E_{y}v_{y}$) the electromagnetic wave contributes with equal amounts of work in each direction. Nevertheless, the electron state can not access any point $(q,a)$ in the stability chart shown in Fig. 2; Eq. (21) imposes an extra constraint. Defining $\Gamma=\epsilon/\hbar\Omega=\sqrt{v_{x}^{2}k_{x}^{2}+v_{y}^{2}k_{y}^{2}}/\Omega$, Eq. (21) takes the form of a straight line $a=\Gamma^{2}+(\zeta_{y}/\hbar\Omega)^{2}-2q$ in the $q-a$ plane. Hence, for any state to be accessible to the electron, the ordered pair $(q,a)$ must satisfy the inequality

[TABLE]

In Fig. 2 the solid purple and solid red lines illustrate the limiting case

[TABLE]

for $\zeta_{x}>\zeta_{y}$ and $\zeta_{y}>\zeta_{x}$ respectively. Naturally, any of these points should also fall on the stable regions allowed by the Mathieu equation in order to produce a stable solution of the wave function.

For fixed $E_{x}$, $E_{y}$ and $\Omega$, $q$ is constant, and therefore the allowed states should be located on the vertical line $q=\mathrm{const}$ (see for example the green dotted line or the gray dotted line in Figs. 2 and 3). Along this lines, the ranges of stable and unstable states alternate producing the appearance of bands separated by dynamical energy gaps. The opening of these gaps is due to the space-time diffraction of electrons in phase with the electromagnetic field, and effect akin to the magnetoacoustic diffraction of electrons in phase with acoustic waves Landau (1946); Davydov (1980); Champo and Naumis (2019).

To further comprehend the connection between the Mathieu stability chart and the consequent wave function gap structure, it is illustrating to project the stable and unstable regions of Fig. 2 on the surface of the tilted Dirac cones that arise from the free particle Weyl equation. To this end, we explicitly express the normalized energy dispersion $\widetilde{E}=(1/\hbar\Omega)E^{+}_{\pm 1,k}$ from Eq. (2) in terms of the normalized wave vector components $(\widetilde{k}_{x},\widetilde{k}_{y})$ and the parameter $a$ in the Eq. (21) obtaining

[TABLE]

where $\widetilde{k}_{x}=(v_{F}/\hbar\Omega)k_{x}$ and $\widetilde{k}_{y}=(v_{F}/\hbar\Omega)k_{y}$. Elliptical rings of allowed and forbidden states form in the $(\widetilde{k}_{x},\widetilde{k}_{y})$ plane or on the surface of the Dirac cone for fixed values of $E_{x}$ and $E_{y}$ (or fixed values of $q$ and $\zeta_{y}/\hbar\Omega$). The dressed Dirac cones, the Dirac cones over whose surfaces the allowed and forbidden states have been projected, are shown in Fig. 4. The light blue portion of the surface represents the allowed states and the white rings are the forbidden ones. The first correspond to the stable regions and the latter to the unstable regions of Fig. 2.

In Figs. 3 (a) and (b) we plot the vertical line $q=-1.77$ (green dotted line) and the line (28) (solid purple line) superimposed to a zoom of the stability chart for typical electromagnetic field values $E_{x}=5.5$ V/m, $E_{y}=2.6$ V/m and $\Omega=50$ GHz. The crossing between these two lines, seen in Fig. 3 (a), is the starting point for the search of stable solutions. However, in the immediate region above the crossing we observe a gap of unstable solutions, that projected on to the Dirac cone produces the appearance of forbidden states at the tip, forming a gap. At higher energies, we observe the rings corresponding to the third order gap as can be appreciated in Figs. 2 and 3 (b). This gap yields and energy range $\Delta^{L}_{3}=13.22\,\mu$eV (see Appendix B) of forbidden states. It should be noted that the origins of the first gap at the Dirac point and the following ring-like forbidden regions are essentially the same. Both of them are generated in points that comply with the inequality (27), and as a result of inherent instabilities of the Mathieu solutions.

When the parameters are chosen to fall on the opposite side of the stability chart ($\zeta_{x}<\zeta_{y}$) the arrangement of the gaps is quite different. In Fig. 3 (c) we observe the crossing of the vertical line $q=0.87$ and the limiting line (28) for electromagnetic field values $E_{x}=2.6$ V/m, $E_{y}=5.5$ V/m and $\Omega=50$ GHz. In contrast to the previous case, above the crossing we find ourselves well inside a stability region. Hence, the tip of the Dirac cone is dressed entirely with allowed states and the forbidden rings appear well above it as can be seen in Fig. 4. At high energies the line $q=0.87$ crosses the second order gap as can be seen in Fig. 3 (d). A ring of unstable states with an energy gap of $\Delta^{R}_{2}=19.45\,\mu$eV is projected on to the Dirac cone until the line reaches the next stability zone [see Fig. 4 (b)].

To systematize the search of gaps in the $K$ point of the Dirac cone we define the indicator

[TABLE]

where $r\in\mathbbm{Z}$ and $c_{r}(q)$ is either $a_{r}(q)$ or $b_{r}(q)$, depending on which one is at the bottom of the allowed band. This indicator corresponds to the energy difference between the lower allowed band edge and the limiting case of the inequality (27) given by (28). The integer $r$ is chosen so that the purple line in Fig. 2 is situated directly below the top band edge associated with the stable region. Therefore, if for a given value of $q$ the purple line falls on a forbidden region of states then $\Delta E(q)>0$. If, on the other hand, the purple line falls on an allowed band $\Delta E(q)$ is a pure imaginary number. The domain where $\Delta E(q)$ is a pure real number corresponds, thus, to a gap of forbidden states. Hence, the function $E(q)$ provides with a clear-cut criterion to detect the formation of gaps in the surroundings of the Dirac point: $\Delta E(q)\in\mathbbm{R}$.

In Fig. 5 we analyze the gap formation through the behaviour of $\Delta E(q)$ in the case $\zeta_{x}>\zeta_{y}$ for fixed $E_{y}=2.6$V/m. The domains of $q$ where $\Delta E(q)>0$ are shown as solid brown lines in Fig. 5 (a). The point given by $E_{x}=5.5$ V/m, $q=-1.77$, examined previously in Fig. 4 (a), is shown in Fig. 5 (a). We notice that this value of $q$ falls inside one of the domains where $\Delta E(q)>0$, therefore indicating the presence of a gap opening at the tip of the cone. Moreover, another useful property of $\Delta E(q)$ is that it gives the energy of the gap. In this example $\Delta E(-1.77)/2=17\mu$eV which corresponds to the gap shown in Fig. 4(a). In contrast, for $E_{x}=2.6$ V/m, $E_{y}=5.5$ ($\zeta_{x}<\zeta_{y}$) and $q=0.89$, Fig. 5 (b) shows that there is no gap formation as expected from Fig. 4.

These results clearly show that if $E_{x}>E_{y}$, the function $\Delta E(q)/2$ presents more domains that yield forbidden gaps close to the Dirac points than $E_{y}>E_{x}$. This strongly suggests that the shape and position of the gaps is strongly influenced by the relative orientation of the minor and major axes of the elliptical profiles of the radiation and the free-Weyl particle Dirac cone. If the major and minor axes of the ellipses arising from the Dirac cone are perpendicular to the major and minor axes of the ellipses of the electromagnetic profile, a gap at the Dirac point is more likely to form. Otherwise, if the ellipses are oriented in the same direction, gaps are more improbable. It is important to emphasize that this does not imply that the two ellipses necessarily have to have the the same proportions. Nevertheless, the only way to correctly predict the formation of a gap in the Dirac point is to determine if the value of $q$ falls inside one of the domains were $\Delta E>0$, as we discussed above.

III.2 Quasi-energy spectrum

The Hamiltonian in Schrödinger equation (6) is a periodic function of time, therefore the Floquet theorem must hold and, consequently, the wave function must be of the form

[TABLE]

The phase $\mathcal{E}$ is usually termed the quasi-energy and $\boldsymbol{f}(\Omega t)=\boldsymbol{f}(\Omega t+2\pi)$ is a periodic function of time with the same period $T=2\pi/\Omega$ as the Hamiltonian. Using Eqs. (23) and (24) we can rearrange the Mathieu functions as

[TABLE]

where $r(a,q)$ is purely rational and therefore $a_{r}=b_{r}$ as we discussed before. Substituting Eq. (32) into (26) and comparing with (31) we find the explicit expression for the quasi-energy

[TABLE]

To better visualize the shape of the spectrum we may use $r\approx\sqrt{a}$ for $r\gg 1$ to approximate the quasi-energy. Using the definition of $a$ (21) and taking the limit for large quasi-momenta we get

[TABLE]

Hence, the quasi-energy spectrum asymptotically approaches the Dirac cone for large quasi-momentum values. This feature is clearly seen in Fig. 6. The light blue straight lines are the section of the Dirac cone cut by the $k_{x}=0$ plane and the solid gray lines are the quasi-energy spectrum in the approximation $r=\sqrt{a}$. The blue, orange, green and red lines correspond to the different bands arising from the exact expression of the quasi-energy (33). We readily confirm that both, the exact and the approximated quasi-energy spectrum lines, asymptotically come close to the Dirac cone. This is also seen in Fig. 7, where the full quasi-energy surface is shown and the Dirac cone is depicted as a light blue semi-transparent surface. The most striking characteristic of these plots is the formation of gaps. In Fig. 6 (a) a gap appears close to the Dirac point for $E_{x}=5.5$ V/m, $E_{y}=2.6$ V/m and $q=-1.77$. It is located inside the same region as the forbidden states in the dressed Dirac cone shown in Fig. 4 (a). In the full spectrum of Fig. 7 (a) this gap translates into a disc-shaped vacuum of states in the tip of the quasi-energy spectrum. Instead, for $E_{x}=2.6$ V/m, $E_{y}=5.5$ V/m and $q=0.87$ there is no gap formation close to the Dirac point, though, a ring-shaped gap appears around the tip of the quasi-energy spectrum as it is shown in Figs. 6 (b) and 7 (b). These results are consistent with the ones found for the dressed Dirac cones. Newly, the relative orientation of the minor and major axes of the radiation and the Dirac cone determine the location and structure of the gaps.

IV Conclusions

We have systematically investigated the wave function stability of charge carries in a Weyl Hamiltonian under an intense elliptically polarized electromagnetic radiation. To this end, we have worked out the time-dependent wave function from the Schrödinger equation in the limit of strong electric field (or long wavelength). We have proven that the stability properties of the wave functions are inherited from the Mathieu functions, in terms of which they are expressed. The analysis of the stability chart of the Mathieu functions projected onto the Dirac cones shows the formation of gaps of unstable states for certain domain regions of the quasi-momentum space. We have shown that the structure of the gaps strongly depends on the alignment between the minor and major axes of the elliptical profiles of the radiation and the Dirac cone of the free-Weyl particle spectrum. In summary, the formation of a gap at the Dirac point is more likely if the radiation and Dirac cone axes do not match, or are perpendicular. Otherwise, if the axes are aligned, ring-shaped gaps form for higher energies. The quasi-energy spectrum extracted from the phase of the wave function consistently reproduces the position and shapes of these gaps. Magnitude estimations of the electromagnetic fields and gap values were presented for the $8-Pmmn$ borophene phase.

V Acknowledgements

This work was supported by DCB UAM-A grant numbers 2232214 and 2232215, and UNAM DGAPA PAPIIT IN102717. V.G.I.S and J.C.S.S. acknowledge the total support from DGAPA-UNAM fellowship.

Appendix A

In this Appendix, we derive Eqs. (6) and (10) from the ansatz (7) and (9). We start from the Dirac equation (6) and applying the solution (7), we obtain Eq. (8), where

[TABLE]

whose matrix elements are given by

[TABLE]

In the previous equatins, the vector components $A_{x}$ and $A_{y}$ only depend on the time variable.

Now, we substitute Eq. (9) in Eq. (8), to get

[TABLE]

where

[TABLE]

From Eq. (38) we obtain

[TABLE]

where

[TABLE]

In the equations above we used the vectors $\hat{\bm{\sigma}}=(\hat{\sigma}_{x},\hat{\sigma}_{y})$ and $\bm{S}=e(v_{x}A_{x},v_{y}A_{y})$. Therefore, Eq. (40) can be reduced into

[TABLE]

with

[TABLE]

and

[TABLE]

where $\boldsymbol{\kappa}_{\boldsymbol{v}}=\hbar(v_{x}k_{x},v_{y}k_{y})$.

Appendix B

To estimate the gap size, first, we calculate the value of the parameter $q$ from the Eq. (20), and subsequently use the expression López-Rodríguez and Naumis (2010); Champo and Naumis (2019)

[TABLE]

where $a_{r}(q)$ and $b_{r}(q)$ for $r\in\mathbbm{Z}$ are the boundaries of a forbidden region as it is shown in Fig. 2.

For a microwave with frequency $\Omega=50$ GHz, we estimate the gap in the following cases: a) For $E_{x}=5.5$ V/m and $E_{y}=2.6$ V/m ( $\zeta_{x}>\zeta_{y}$), we find $q=-1.77$ and $\Delta^{L}_{3}=13.22\,\mu$ eV. b) For $E_{x}=2.6$ V/m and $E_{y}=5.5$ V/m ($\zeta_{x}<\zeta_{y}$), we find $q=0.87$ and $\Delta^{R}_{2}=19.45\,\mu$ eV.