Exploring the Optical States for Black Phosphorus: Anisotropy and Bandgap Tuning
Andrii Iurov Liubov Zhemchuzhna Godfrey Gumbs Danhong Huang





aiurov@unm.edu

Center for High Technology Materials, University of New Mexico,
1313 Goddard SE, Albuquerque, NM, 87106, USA


Department of Physics and Astronomy, Hunter College of the City University of New York,
695 Park Avenue, New York, NY 10065, USA

Department of Physics and Astronomy, Hunter College of the City
University of New York, 695 Park Avenue, New York, NY 10065, USA
Donostia International Physics Center (DIPC),
P de Manuel Lardizabal, 4, 20018 San Sebastian, Basque Country, Spain

Air Force Research Laboratory, Space Vehicles Directorate,
Kirtland Air Force Base, NM 87117, USA
Center for High Technology Materials, University of New Mexico,
1313 Goddard SE, Albuquerque, NM, 87106, USA





The dressed states arising from the interaction between electrons and holes, and off-resonant electromagnetic
radiation have been investigated for recently fabricated gapped and anisotropic black phosphorus.  Our calculations
were carried out for the low-energy electronic subbands near the Γ point. States for both linear 
and circular polarizations of the incoming radiation have been computed. However, our principal emphasis 
is on linearly polarized light with arbitrary polarization since this case has not been given much attention
for dressing fields imposed on anisotropic structures. We have considered various cases for 
one- and few-layer phosphorus, including massless Dirac fermions with tunable in-plane anisotropy. Initial 
Hamiltonian parameters are renormalized in a largely different way compared to those for previously 
reported for gapped Dirac structures and, most importantly, existing anisotropy which could be modified 
in every direction. 


78.67.-n, 78.67.Wj, 81.05.Xj, 73.22-f.anisotropic

Exploring the Optical States for Black Phosphorus: Anisotropy and Bandgap Tuning 
    Danhong Huang
    December 30, 2023
=================================================================================





§ INTRODUCTION




Black phosphorus (BP) is a layered structure of buckled atomic phosphorus for which the layers are 
connected by weak van der Waals forces. <cit.> This is the most stable phosphorus-based crystal 
at room  temperature and for a wide range of pressure values. Crystalline BP is a semiconductor, for 
which P atoms are covalently bonded with three adjacent atoms.  Such a structure exhibits strong 
anisotropy for an arbitrary number of layers, as well as hybrid electron and hole 
states in the vicinity of the band edge in phosphorene featuring both Dirac cone and a conventional 
Schrödinger electron behavior. While bulk BP is a semiconductor with a small 
bandgap of 0.3 eV , its monlayer counterpart is a semiconductor with a large direct bandgap 
(∽ 2 eV) and is referred to as phosphorene in analogy to graphene and could be 
exfoliated in a mechanical way.



It is not surprsing that reliable information for the band structure of BP with 
a specified number of layers has been appreciated as being extremely important for device applications,
including a field-effect transistor. <cit.> This has stimulated a large number of first-principles
calculations based on the density functional theory as well as the tight-binding model  <cit.>,
group-theoretical calculations  <cit.> and the continuum model. <cit.>
 of BP structures, which recently received a significant amount of attention for their success 
in analyzing experimental data. <cit.> The quality of all such 
models  <cit.> depends on the accuracy of the exchange-correlation approximation. 
Over time, we  are going to witness numerous attempts to engineer new types of anisotropic electronic 
bandstruture of BP-based devices using various mechanisms such as electron-photon interaction with a 
dressing field. The corresponding effect from an imposed electrostatic field 
was addressed in Ref. [elfield].




 The possibility of an electronic topological transition and electrically-tunable Dirac cone was 
theoretically predicted for multi-layer BP. <cit.> The fact that the bandgap of BP 
is mainly determined by the number of layers was confirmed in experiment,<cit.> i.e., the 
energy gap is strongly decreased for a larger number of layers and could be effectively neglected 
for N_L > 5. Yet demonstrating a Dirac cone with no effective mass or energy 
bandgap, such electrons still posses strong anisotropic properties for pristine black 
phosphorus.  <cit.>     Consequently, non-symmetric Klein tunneling could be 
observed  <cit.> with a shifted transmission peaks similar to electron tunneling 
for graphene in the presence of magnetic field. <cit.>

 

Phosphorene is one of the most recently discovered members of a sizable group of low-dimensional structures 
with great potential for nanoelectronics applications. The most famous part of this family 
is graphene, fabricated in 2004. Because of its unique properties, <cit.> graphene
has initiated a new direction in electronic devices.  A subsequent crucial advance was the discovery 
of the buckled honeycomb lattices such as silicene and germanene. Their most distinguished feature 
is the existence and tunability of two non-equivalent spin-orbit and sublattice-asymmetry bandgaps. 
The latter one, Δ_z is directly modified by an external 
electrostatic field. <cit.> 
This is a result of the out-of plane buckling, which is due to a larger radius 
of Si and Ge atoms compared to carbon and sp^3 hybridization. The most 
recently fabricated germanene possesses similar buckling features but with different bandgaps and
Fermi velocity. <cit.>



Another important class of such materials are the transition metal dichalcogenides 
(TMDC's) structures such as MC_2, where M  denotes a metal such as Mo, W,  
and C is a chalcogen atom (S, Se, or Te). Molybdenum disulfide 
MoS_2, a typical representative  <cit.> of TMDC's, has exhibited a large energy 
bandgap ⋍ 1.9  eV, and broken symmetry between the electron and hole 
subbands so that for all experimentally accessible electron density values only 
one hole subband is doped. As a result, all the electronic, collective and transport 
properties vary significantly for the electron and hole types of doping. However, all 
these low-dimensional materials exhibit almost complete (with a slight 
deviation for MoS_2) isotropy in the x-y plane. Consequently, phosphorene 
is a totally new material with its very unusual properties, so that complete 
and thorough studies of these characteristics open a new important chapter in 
low-dimensional science and technology. 







With the newest achievements in laser and microwave science, it has become possible to 
achieve substantial control and tunability of the electronic properties of low-dimensional
condensed-matter materials by subjecting them to strong off-resonant high-frequency 
periodic fields (so-called "Floquet engineering“, schematically shown in 
Fig. ]refFIG:1 (a)).<cit.> If the electron-phonon coupling 
strength is high, such a bound system could be considered as a single, holistic 
object and has been investigated using quantum optics and mechanics. These 
electrons with substantially modified energy dispersions, referred to as  
“dressed states"', became a commonly used model in present-day
low-dimensional physics. <cit.>
One of the first significant achievements has been the demonstration of a
metal-insulator transition in graphene  <cit.>, which drastically 
affected the electron tunneling and the Klein 
paradox. <cit.> Important collective properties such as exchange and 
correlation energies are also affected by the presence of an energy gap, <cit.> 
and spin dynamics on the surface of a three-dimensional topological insulator  <cit.> is also modified.

The rest of our paper is organized in the following way. In Sec. <ref>, 
we present our model, the low-energy Hamiltonian and the energy dispersions 
for phosphorene, i.e., single-layer black phosphorus with strong in-plane 
anisotropy and a large energy bandgap. The electron-photon dressed states 
for phosphorene are presented and discussed in Sec. <ref>.  Section
   <ref> is devoted to calculating the dressed states for few-layer
phosphorus, in which the electrons are anisotropic massless Dirac 
fermions without a gap. Mathematical generalizations of such a model with 
both on- and off-diagonal bandgaps is considered, and the corresponding 
dressed states are also obtained. Concluding remarks are provided in Sec. <ref>.





§ LOW-ENERGY MODEL FOR PHOSPHORENE



Our calculations utilize the model for BP presented in Refs. [GGBpPRL, GGBpPRB].
Being somewhat similar to the truly two-dimensional hexagon structure of graphene, the atomic 
arrangement for single layer BP results in a packered surface due to the  sp^3 
hybridization composed of the 3s and 3p orbitals. For silicene, such hybridization 
is responsible for out-of-plane “buckling" displacement of the Si atoms. 



The continuouum k-dependent Hamiltonian is usually based on the tight-binding model.
Close to the Γ point, approximated up to second order in the wave vector components, 
it is given as


    H_ph^Δ  ( k) = ( 
    E_i + ∑_i = x,yη_i k_i^2
    ) Î_2× 2+
    ( 
    Δ_O + ∑_i = x,yγ_i k_i^2
    ) Σ̂_x-
    χ k_y Σ̂_y   ,

or in the matrix form, 


    H_ph^Δ  ( k) = [
    [            E_i + ∑_i = x,yη_i k_i^2 Δ_O 
    + ∑_i = x,yγ_i k_i^2 - i χ k_y;  Δ_O + ∑_i = x,yγ_i k_i^2 + i χ k_y            E_i + ∑_i = x,yη_i k_i^2 ]]   .

This Hamiltonian clearly displays significantly different structure and properties, 
compared to that for graphene. First, there are no linear k- terms, except 
± i χ k_y. Furthermore, there are no linear k_x elements. As one of the 
most evident consequences of this structure, we note that circularly polarized 
irradiation with x- and y-components being equally important, couples to 
such electrons only in the ⋍ A^2 level.




Secondly, the energy bandgap is presented in a Σ̂_x, off-diagonal 
form, contributing to the asymmetry between the electron and hole states in 
contrast to the Σ̂_z-type gap. These properties, coming directly 
from the Hamiltonian structure, are new and have not been encountered previously. 
The energy dispersions are 


    ε_Δ^± ( k) = E_i + ∑_i = x,yη_i k_i^2 ±[( Δ_O  + ∑_i = x,yγ_i k_i^2 )^2 
    + χ^2 k_y^2  ]^1/2  ,

where γ_e,h = ± 1 corresponds to the electron or hole solution. For    
small values of the wave vector, these dispersions are approximated as


    ε_Δ^± ( k) ⋍E_i ±Δ_O + ( η_x ±γ_x ) k_x^2 +
    [   η_y ±(
    γ_y + χ^2/2 Δ_O)
    ] k_y^2  .

The effective masses given by  <cit.>


    m^(e,h)_ x =  ħ^2/2 ( η_x ±γ_x )  , 
           m^(e,h)_ y =  ħ^2 / 2/η_y ±(
    γ_y + χ^2 / (2 Δ_O) )

are anisotropic, and this anisotropy is different for the electron and hole 
states as ∽χ^2/Δ_O. 








§ ELECTRON DRESSED STATES FOR IN SINGLE LAYER



In this Section, we calculate electron-light dressed states for phosphorene.  
As far as circularly polarized irradiation is concerned, one must consider second
order coupling in order to see how both components of the wave vector are modified.
Such consideration is critical for the vector potential given by Eq.(<ref>),
but is clearly beyond the scope of conventional analytical methods. The mere presence of an 
off-diagonal energy gap Δ_O means that there is no electron/hole symmetric solution,
obtained in Refs. [kibissrep,kibisall]. This situation also leads us to 
conclude that the Hamiltonian parameters, such as the energy gap, are affected at 
lower order than such parameters for Dirac fermions. Consequently, for   monolayer
BP, we focus on the case for linearly polarized irradiation.



 §.§ Linear polarization of the dressing field and induced anisotropy


Since electron energy dispersion relations and their effective masses are intrinsically
anisotropic for phosphorenes, the direction of the dressing field polarization 
now plays a substantial role. We define this direction by an arbitrary angle 
θ_0 from the x-axis and generalize the vector potential used in 
Ref. [kibisall] so that it now has both x- and y- non-zero components 


    A^L(t) =  {E_0/ωcos ω tcosθ_0;  E_0/ωcos ω t sinθ_0 }  .




The renormalized Hamiltonian for the dressed states is obtained by the canonical substituion 
k_x,y⇒ k_x,y - e /ħ A_x,y,  where e stands for the electron charge, yielding


    Ĥ(k) = H_ph^Δ + ĥ_0 + ĥ_int  ,

where H_ph^Δ is the “bare", non-interacting electron Hamilotonian, 
given by Eq.(<ref>). The zeroth-order, k-independent 
interaction Hamiltonian may be expressed as


    ĥ_o = χ E_0 e/ħω sinθ_0  cosω t  Σ̂_y = 
    c_0  (
    [         0 -i cosω t;  i cosω t         0 ])  ,

where c_0 = χ sinθ_0    E_0 e/(ħω), e = | e |. 
We conclude that the vector potentail of linearly polarized light (<ref>)  
must have a non-zero y- component in order to enable ∽ A^L coupling. 
We now turn our attention to the interaction term which is linear in k_x,y and given by


    ĥ_int =
     2 e/ħ( 
     ∑_i = x,yη_i k_i   A_i^L(t)  I_2 × 2 + 
     ∑_i = x,yγ_i k_i   A_i^L(t)  Σ̂_x
     )   ,

 or, introducing the following simplifying notations 
ϵ_α = √(α_x^2 k_x^2 + α_y^2 k_y^2), 
 ϕ^(α) = tan^-1[ α_y k_y / (α_x k_x) ] 
for α = η, γ and 
c^(2) = (2 e E_0) / (ħω), we express it as


    ĥ_int = c^(2) cos ω t  [[ ϵ_ηcos (ϕ^(η) - θ_0) ϵ_γcos (ϕ^(γ) - θ_0); ϵ_γcos (ϕ^(γ) - θ_0) ϵ_ηcos (ϕ^(η) - θ_0) ]]   .

As the fisrt step, we solve a time-dependent Schrödinger equation 
for k = 0:


    i ħ∂Ψ_0(t)/∂ t = ĥ_oΨ_0(t)   .

We obtain the solution in a  straightforward way to be


    Ψ_0^β = ± 1(t) = 1/√(2) [
     [    1; β  i ]] exp{
     - i β c_0/ħω sinω t
     }  .

It is noteworthy that even if both energy bandgaps Δ_O and the initial 
energy shift E_i are included so that the Hamiltonian takes the form


    h_o = E_i  I_2 × 2 + Δ_O  Σ̂_x + c_0  cosω t  Σ̂_y = 
     c_0  (
     [            E_i  Δ_O -i cosω t; Δ_O + i cosω t            E_i ])   ,

the solution could still could be determined analytically as 


    Ψ_0^β = ± 1(t) = 1/√(2) [
     [    1; β  i ]] exp{
     - βi/ħ[
     c_0/ω sinω t - (Δ_O - βE_i)   t
     ]
     }  .

We note that such a trivial solution could only be obtained for equal diagonal energies 
E_0, i.e., for an introduced non-Σ_z type of energy gap. If the diagonal energies
are not equivalent, given by E_1,2 = E_i ±Δ_D, complete symmetry between 
the two components of the wave function no longer exists and such a solution could 
not be determined analytically. Consequently, we will use basic set in 
Eq.(<ref>) for the rest of the present calculation. Such a situation (non-zero diagonal 
gap Δ_D) appears if there is a relatively small non-zero vertical electrostatic 
field component and the symmetry between the vertically displaced phosphorus atoms
is broken similar to silicene. We will often use this mathematical generalization in the rest
of our work.



For finite wave vector, we present the eigenfunction as an expansion in terms of
a basis set  <cit.>


    Ψ_ k (t) = F^⇑ Ψ_0^+(t) + F^⇓ Ψ_0^-(t)   ,

where F^⇑,⇓ = F^⇑,⇓(k_x,k_y  |  t) are scalar, time-dependent coefficients with anisotropic k-dependence. 
This equation immediately results in the two following indentities


    i ħd  F^ ⇑,⇓/dt = 
     ⟨Ψ_0^±(t) |δH_ph^γ( k) |Ψ_0^+(t) ⟩ F^ ⇑ + ⟨Ψ_0^±(t) 
    |δH_ph^γ( k) |Ψ_0^-(t) ⟩ F^ ⇓  ,

where 


    δH_ph^γ( k) = H_ph^γ + ĥ_int

is the bandgap and wave vector dependent portion of the total Hamiltonian (<ref>). 
This system becomes


    i d  F^ ⇑,⇓/dt       = 
    [
    E_i + η_x k_x^2 + ( η_y k_y ±χ ) k_y + c^(2)ϵ_ηcos (ϕ^(η) - θ_0 )  cosω t] 
     F^ ⇑, ⇓±
          ±   
    i [
    ∓ i Δ_D + Δ_O + γ_x k_x^2 + γ_y k_y^2 
    + c^(2)ϵ_γcos (ϕ^(γ) - θ_0 )  cosω t]
     F^ ⇓, ⇑ exp[
    ± 2  i   c_0/ħω  .
     sinω t ]


The quisiparticle energy dispersion relations ε_d ( k) 
are calculated by using the Floquet theorem from the following substitution


    F^⇑, ⇓(t) = exp[ - i ε_d ( k) /ħ  t ] 
    ∑_λ = - ∞^∞ f_λ^⇑, ⇓ e^i λω t  ,

where f^⇑, ⇓(t) = f^⇑, ⇓(t + 2 π/ ω) are 
time-dependent periodic functions. The nested exponential dependence is usually 
simplified using the Jacobi-Anger identity


    exp[
    ± 2  i   c_0/ħω  sinω t 
    ] = 
    ∑_ν = - ∞^∞e^i νω t  J_ν(
    ± 2 c_0/ħω) .

The orthonormality of Fourier expansions results in the following system of 
2 μ, μ⇒∞ equations


    ε_d ( k)   f^⇑, ⇓_μ = 
    [
    μ ħω + E_i ±η_x k_x^2 + ( η_y k_y ±χ ) k_y
    ]  f^⇑, ⇓_μ +     [ c^(2)/2ϵ_ηcos (ϕ^(η) -θ_0) ]
     ( f^⇑, ⇓_μ+1 + f^⇑, ⇓_μ-1) 
    ∑_λ = -∞^∞
     f^⇓, ⇑_λ×
          × {[ Δ_D 
     ±  i  ( Δ_O + γ_x k_x^2 + γ_y k_y^2 )
     ]    J_μ-λ( ± 2 c_0/ħω) + [ c^(2)/2ϵ_γcos (ϕ^(γ) - θ_0) ]  ∑_α = ± 1
     J_μ+ α -λ( ± 2 c_0/ħω) }  .
  
In our consideration, the frequency of the off-resonant dressing field is high 
enough so that only diagonal elements in the eigenvalue equation (<ref>) are 
retained. However, if we need to include the first-order electron-field coupling terms
∽ c^(2), we must keep the summations with λ = μ± 1. 



In the simplest case, where only diagonal elements are kept, the quasiparticle 
energy dispersion relations are    
 

    ε_d ( k) = E_i + ∑_i=x,yη_i k_i^2 ±{χ^2 k_y^2 + [Δ_D^2 + (Δ_O + ∑_j=x,yγ_j k_j^2 )^2 ]   J_0^2  (
    2 c_0/ħω) }^1/2  .

This result is valid only if c_0 = χ sinθ_0    E_0 e/(ħω) 
≠ 0, or there is a finite y-component of the polarization direction 
of the dressing field.


    ε_d ( k) = E_i ±{(1 - α_c^2) [
     Δ_D^2 + Δ_O^2
     ]
    }^1/2 + [
    η_x ±γ_x  √(1 - α_c^2) Δ_O/√(Δ_D^2 + Δ_O^2)]   k_x^2 + [
    η_y ±γ_y  Δ_O  (1 - α_c^2) + χ^2/2/√((1 - α_c^2 ) [ Δ_D^2 + Δ_O^2 ])]  k_y^2   ,

where α_c  = 2 c_0 / (ħω) is a dimensionless coupling coefficient. 
The electron effective masses are now readily obtained. If there is no Σ_z-type 
energy bandgap Δ_D, the expressions  are simplified as


    m^(e,h)_ x =  ħ^2/2 ( η_x ±α̃_c  γ_x )  , 
           m^(e,h)_ y =  ħ^2 / 2/η_y ±[ α̃_c  γ_y + χ^2 / (2 Δ_O α̃_c  ) ]  ,

where α̃_c = √(1 - α_c^2)∽ 1 - α_c^2 / 2. 
The obtained energy dispersion relations are presented in Fig. <ref>. 
It is interesting to note that both energy bandgaps are renormalized by the 
electron-photon interaction, showing a substantial decrease, while the diagoanl 
terms for the initial effective masses of a “bare" electron η_x,y
are unchanged. 








§ ANISOTROPIC MASSLESS FERMIONS IN FEW-LAYER PHOSPHORUS




The central property, as well as the research focus of phosphorene, is the electron 
dispersion relation and effective mass anisotropy. At the same time, BP-based materials
have a band gap which is determined by the number of layers, which varies from 
0.6 eV for five layers to 1.5 eV for a single layer. Specifically, we consider 
anisotropic massless Dirac particles, which could be observed in special few-layer 
(N_L > 5) black phosphorus superlattices for a narrow range of energies. 




This anisotropic Dirac Hamiltonian 


    H_ml^γ_0 = ħ v_F' ( k_x Σ̂_x  + γ_0 k_y  Σ̂_y  )

leads to the following energy dispersion


    ε_γ_0^± ( k) = ±ħ v_F' √(k_x^2 + (γ_0 k_y)^2)  .

For such fermions interacting with light having linear polarization in arbitrary 
direction θ_0 described by Eq.(<ref>), we obtain the following Hamiltonian


    Ĥ( k) =  H_ml^γ_0( k) + ĥ_0   ,

where


    ĥ_0 = e v_F   Σ̂· A^L(t) = e v_F E_0/ω (
     [                0 e^ - i θ_γcosω t;   e^ i θ_γcosω t                0 ])

is the k = 0 portion of the total Hamiltonian. Here, we also introduced 
θ_γ = tan^-1[γ_0 tan(θ_0)] so that 
e^± i θ_γ = cosθ_0 ± i γ_0 sinθ_0. 
We define c_0 = e v_F E_0 / ω to be the electron-photon
interaction coefficient with the dimension of energy and for a given range of frequency 
c_0 ≪ħω  - dressing field which cannot be absorbed by electrons.
Traditionally, we first need to solve the time-dependent Schrödinger equation 
for k = 0 and Hamiltonian ĥ_0. The eigenfunction is obtained in a 
straightforward way as



    Ψ_0^β = ± 1(t) = 1/√(2) [
     [        1; βe^i θ_γ ]] exp{
     - i β c_0/ħω sinω t
     }  .


In order to determine the solution for a finite wave vector we once again employ
the expansion (<ref>) and solve Eq. (<ref>) for the time-dependent coefficients
F^⇑,⇓ = F^⇑,⇓(k_x,k_y  |  t). 
In our case, this leads to  
    

    i d  F^ ⇑,⇓/dt = ± v_F' k_γcos(ϕ^(γ) - θ_γ)  F^ ⇑,⇓± i   v_F'k_γsin(ϕ^(γ) - θ_γ)  F^ ⇓,⇑ exp[± 2  i   c_0/ħω  sinω t 
     ]   ,

where ϕ^(γ) = tan^-1[γ_0 k_y/k_x] or ϕ^(γ) 
= tan^-1[γ_0 tan(ϕ_0)] and tanϕ_0 = k_y/k_x. We also denote
k_γ = √(k_x^2 + (γ_0 k_y)^2).

Now we also use the Floquet theorem to extract the quasiparticle energy  ε_d ( k) 
and expand the remaining time-periodic functions as a Fourier series. The result is 
again a system of 2 μ, μ⇒∞ equations


    ε_d ( k) f^⇑, ⇓_μ = ∑_λ = -∞^∞{δ_μ, λ [
      μ ω±ħ v_F' k_γcos (ϕ^(γ) - θ_γ)
      ]  f^⇑, ⇓_λ±  i  ħ v_F' k_γsin (ϕ^(γ) - θ_γ)    J_μ-λ(
    ± 2 c_0/ħω)   f^⇓, ⇑_λ}  .







In the region of interest, i.e., for large frequency ω≫ v_F' k and 
ω≫ϵ ( k) we approximate  
f^⇑, ⇓_μ≠ 0⋍ 0. Finally, the eigenvalue 
equation becomes ε_d ( k) f^⇑, ⇓_0 = K(⇑, ⇓ | γ_0,  k) 
× f^⇑, ⇓_0, where


    K(⇑, ⇓ | γ_0,  k) = [
     [                                  ħ v_F' k_γcos (ϕ^(γ) - θ_γ) i γ_0  ħ v_F' k_γsin (ϕ^(γ) - θ_γ)  
        J_0 [
    2 c_0/(ħω)
     ];     - i ħ v_F' k_γsin (ϕ^(γ) - θ_γ)    J_0 [
    - 2 c_0/(ħω)
     ]                                - ħ v_F' k_γcos (ϕ^(γ) - θ_γ) ]]

The energy eigenvalues are given by


    ε_d ( k) = ±ħ v_f' √(k_x^2 + (γ_0 k_y)^2) {cos^2 (ϕ^(γ) - θ_γ) +[ sin (ϕ^(γ) - θ_γ)   J_0 (
    2c_0/ħω) ]^2   
     }^1/2 .

For small light intensity c_0 ≪ħω the zeroth-order Bessel function of
the first kind behaves as


    J_0^c_0 ≪ħω[
    2c_0/ħω] ⋍ 1 - c_0^2/(ħω)^2 + c_0^4/4 (ħω)^4 - ...

and we have approximately for the energy dispersion 


    ε_d ( k) = ±ħ v_f' {[1 - 2c_0^2/(ħω)^2sin^2 θ_γ]^2 k_x^2 + γ_0^2  [1 - 2c_0^2/(ħω)^2cos^2 θ_γ]^2 k_y^2 + 2γ_0 c_0^2/(ħω)^2 k_x k_y
    sin (2 θ_γ) 
    }^1/2  .

If the light polarization is directed along the x-axis, then θ_0 = θ_γ = 0 
and 


    ε_d ( k) = ħ v_f' √(
     k_x^2 + γ_0^2  [1 - 2 c_0^2/(ħω)^2]^2 k_y^2 )  .

Angular dependence of the dressed states dispersions is shown in Fig. <ref>.
We notice that initially existing in-plane anisotropy is affected for all in-plance 
angles, depending on the dressing field polarization direction. For small intensity 
of the incoming radiation, polarized along the x- axis, the anisotropy 
coefficient is simply renormalized. 




 §.§ Circular Polarization




For circular polarization of the dressing radiation, the vector potential is 


    A^C(t) = { A_0,x,   A_0,y} = E_0/ω{cos ω t,   sin ω t }  .

Being completely isotropic, this type of field is known to induce the metal-insulatron 
transition in graphene,<cit.> resulting in the creation of a non-zero energy bandgap. 
If such a gap already exists, it could be increased or decreased depending
on its initial value. <cit.> At the same time, the slope of the Dirac dispersions, 
known as Fermi velocity and the in-plane isotropy are not changed. The situation cannot be 
the same for an initially anisotropic Dirac cone for AMF's.

The total Hamiltonian for the interacting quasiparticle now becomes


    Ĥ( k) =  H_ml^γ_0( k) + ĥ_0^(c)  ,

where the k = 0 interaction term is


    ĥ_0^(c) = e v_F E_0/ω Σ̂· A^C (t) = 
     c_0/2 [
     [                            0 ∑_α = ± 1 (1 - αγ)  e^iαγω t; ∑_α = ± 1 (1 + αγ)  e^iαγω t                            0 ]]  .

It seems rather surprising, although physically justified, that this problem is 
mathematically identical to the isotrpic Dirac cone interacting with elliptically 
polarized light addressed in Refs. [kibisall, goldprx]. The interaction 
term also could be presented as 


    ĥ_0^(c) = Ŝ_γ e^i ω t 
    + Ŝ_γ^† e^ -i ω t
          Ŝ_γ = c_0/2∑_α = ± 1 (1 + αγ) Σ̂_α  ,

where Σ̂_± = 1/2 (Σ̂_x ± i Σ̂_y ).
This Hamiltonian represents an example of a wide class a periodically driven quantum 
systems. <cit.> Such problems are generally solved perturbatively, in powers of 
c_0/(ħω), if the electron-field coupling is weak. The effective Hamiltonian 
for such problem has been shown to be
 

    Ĥ_eff( k) ⋍H_ml^γ_0( k) + 1 / (ħω)  [  Ŝ_γ Ŝ^†_γ ] +  1/2 (ħω)^2 { [ [  Ŝ_γ, H_ml^γ_0( k)  ]  Ŝ^†_γ ] +  h.c. } + ⋯

Evaluating this expression, we obtain 

    Ĥ_eff( k) =   ħ v_F' k_x (
      1 -γ_0   c_0^2/2 (ħω)^2)  Σ̂_x + ħ v_F'  γ_0 k_y (
      1 -c_0^2/2 (γ_0  ħω)^2)  Σ̂_y  -c_0^2/ħω γ_0  Σ̂_z   .

Finally, the energy dispersion is given by
   

    ε_d ( k) = ±{(γ_0   c_0^2/ħω)^2 + ħ^2 v_F'^2 [
     (
      1 -γ_0   c_0^2/2 (ħω)^2)^2 k_x^2 + γ_0^2 (
      1 -c_0^2/2 (γ_0  ħω)^2)^2 k_y^2 
     ]   }^1/2  .

This result is an approximation. As we mentioned, in the case of an isotropic 
Dirac cone, electrons in graphene interacting with a circularly-polarized dressing 
field, the energy gap was found to be  <cit.>


    Δ_g/2 = √(ħ^2 ω^2 + 2   c_0^2 ) - ħω⋍c_0/ħω - 1/4( c_0/ħω)^2 
     + ...     ,

while the Fermi velocity v_F is unaffected.



 §.§ Gapped anisotropic fermions


We now present a generalization of previously considered massless Dirac particles 
with a finite energy bandgap. Two different gaps added to both on- (Δ_D) and 
off-diagonal terms (Δ_O) of the Hamiltonian. Here, an anisotropic Dirac cone 
is combined with the energy gaps attributed to phosphorene, a single-layer structure.  
Even though this model does not exactly describe any of the fabricated black phosphorus 
structures, we conider it as an interesting mathematical generalization of the anisotropic 
Dirac fermions case, which may become relevant from a physical point of view. 
Apart from that, this is an intermediate case between phosphorene and few-layer 
gapless materials, which is expected to approximate the electronic properties of 
a system with a small number of phospohrus layers.  We have 


    H_g =        E_i  Î_2× 2 +
    ( ħ v_F' k_x + Δ_O ) Σ̂_̂x̂  + ħ v_F' k_y  Σ̂_y + Δ_D 
    Σ̂_z   = 
    
     =       [
    [ E_i + Δ_D         0;         0 E_i - Δ_D ]]  + Δ_O
    [  [ 0 1; 1 0 ]]  +
     ħ v_F' [
     [               0 k_x - i γ_0 k_y; k_x + i γ_0 k_y               0 ]]   .

The corresponding energy dispersion is given by 


    ϵ_γ_0^±(k) = E_0 ±{Δ_D^2 + [Δ_O + ħ v_F' k_x ]^2 + (γ_0  ħ v_F k_y)^2
    }^1/2 .

We now address the interaction of these gapped Dirac electrons with the dressing 
linearly polarized field. The vector potential here
is again specified by Eq.(<ref>) and the new Hamiltonian is


    Ĥ( k) =  H_g( k) + ĥ_0   ,

where ĥ_0 is identical to Eq.(<ref>) since both equations share 
similar k-dependent terms. Following this approach, we expand the wave function 
for a finite wave vector Ψ_ k (t) over the basis (<ref>) and obtain 
the following equations for the expansion coefficients 
F^ ⇑,⇓ (k_x, k_y  |  t)


    i Ḟ^ ⇑,⇓       = ±{±E_i + (ħ v_F' k_x + Δ_O)cosθ_γ + γ_0  ħ v_F' k_y sinθ_γ} F^ ⇑,⇓  + 
           + i {
    - i Δ_D ±γ_0  ħ v_F' k_y cosθ_γ∓ (ħ v_F' k_x + Δ_O) sinθ_γ} F^ ⇓,⇑ exp[ ± 2 i  c_0/ħω  sinω t 
     ]   .

Similar to our previous case, we introduce the following simplifying notation 


    ϕ^(O) = tan^-1{γ_0  ħ v_F' k_y
     /ħ v_F' k_x + Δ_O
     }
          ϵ_O = √((ħ v_F' k_x + Δ_O)^2 + (γ_0  ħ v_F' k_y)^2) .

After performing the Floquet theorem substitution and the expansions similar
to the previously adopted procedure and we obtain


    ε_d ( k) f^⇑, ⇓_μ =         [
      μ ω + E_i ±ϵ_Ocos (ϕ^(O) - θ_γ)
      ]  f^⇑, ⇓_λ + 
      
      +       ∑_λ = -∞^∞[ Δ_D 
     ±  i  ϵ_Osin (ϕ^(O) - θ_γ) 
     ]
       J_μ-λ(
    ± 2 c_0/ħω)   f^⇓, ⇑_λ  .

It is interesting to note that the diagonal bandgap Δ_D is now on the 
main diagonal and is affected by the Bessel function. The energy dispersions are now
given by


    ε_d ( k)        = E_i ±{[ϵ_O cos (ϕ^(γ) - θ_γ)]^2 +[ Δ_D^2 +
      ϵ_O^2  sin^2 (ϕ^(γ) - θ_γ)]   J_0^2 (
    2c_0/ħω)    
     }^1/2  , 
     
    ϵ_O       = √((ħ v_F' k_x + Δ_O)^2 + (γ_0  ħ v_F' k_y)^2)   .

If the electron-photon interaction is small, then α_c = 2 c_0/(ħω) ≪ 1,
the energy dispersion relation is approximated by


    ( ε_d ( k) - E_i )^2        ⋍ (1 - α_c^2)  Δ_D^2 +
    [1 - α_c^2 sin^2 θ_γ]^2 (ħ v_F' k_x + Δ_O)^2 + (ħ v_F')^2  γ_0^2  [1 - α_c^2 cos^2 θ_γ]^2 k_y^2 + 
           +    2 γ_0  α ħ v_F' (ħ v_F' k_x + Δ_O)
    k_y  sin (2 θ_γ)   .

From Eqs. (<ref>) and (<ref>)  we note that the diagonal Δ_D bandgap 
is decreased as ∽α_c^2, similar to that for gapped graphene or the 
transition metal dichalcogenides, <cit.> while the off-diagonal gap 
modification drastically depends on the direction of the radiation polarization. 
This behavior has no analogy in the previously considered structures. At the same time, 
the Fermi velocity components are modified similarly to that for anisotropic massless 
fermions. 



§ CONCLUDING REMARKS





We have derived closed-form analytic expressions for electron-photon dressed states in one- and few layer 
black phosphorus. The energy gap is determined by the number of layers comprising the system, reaching its largest
value for a single layer (phosphrene) and effectively vanishes for a few-layer structure. The latter case gives rise to   
anisotropic massless fermions, which exhibits an anisotropic Dirac cone. Since the anisotropy is the most significant 
and common property. Of all cases of BP, we focused on linearly polarized dressing field in an arbitrary direction.
As a result, we demonstrated the Hamiltonian parameters are modified in an essentially different way compared to
all isotropic Dirac systems. Anisotropy of the energy dispersion is modified in all directions, as well as all the 
electron effective masses.   If both diagonal and off-diagonal gaps are present, the latter one remains unaffected, 
but only for a specific light polarization direction. For AMF's interacting with circularly polarized light, 
the problem is mathematically identical with Dirac electrons irradiated by the field with elliptical polarization. 
In that case we found that initially absent energy bandgap is created and the non-equivalent Fermi velocities 
in various directions are renormalized. These results are expected to be of high importance for electronic 
device applications based on recently discovered and fabricated black phosphorus. 

. 
D.H. would like to thank the support from the Air Force Office of Scientific Research (AFOSR).



