Equations of macroscopic electrodynamics for two-dimensional crystals

The electrodynamics of two-dimensional (2D) dielectric and conducting layers cannot be described by such three-dimensional macroscopic quantities as the dielectric constant $\epsilon$ or the refractive index $n$. By means of the proper averaging of the microscopic Maxwell equations we derive general macroscopic electrodynamic equations for 2D crystals and discuss some of their consequences.

by their averaging over "physically infinitesimal" volume elements 1 , i.e. over dimensions small as compared to the macroscopic scales, such as the wavelength of radiation and the sample dimensions, but large as compared to the inter-atomic distances. Here e, h, ρ and j are the microscopic electric and magnetic fields, charge and current densities. As a result of the averaging one defines the macroscopic electric field E =ē, the polarization vector P = χE, which has the meaning of the dipole moment of a volume element, the electric induction D = E + 4πP, the dielectric susceptibility χ, the dielectric constant ǫ = 1 + 4πχ and other macroscopic quantities. The dimensionless quantities χ and ǫ (we use the more physical Gaussian system of units), which can in general be functions of the frequency ω and the wave-vector q of the electromagnetic field, fully determine the linear response of the medium to the electromagnetic field. In the nonlinear optics one defines 2 the higher order susceptibilities, e.g., the third order tensor χ (3) ijkl which determines the cubic term in the Taylor expansion of the function P(E).
The discovery of graphene 3-5 -a one-atom-thick layer of carbon atoms -triggered great interest to this, as well as other two-dimensional (2D) materials 6 , both metallic (graphene), dielectric (e.g. BN) and semiconducting (e.g. MoS 2 ). The notions of the dielectric susceptibility χ, the dielectric constant ǫ and the refractive index n = √ ǫ are no longer applicable to these materials since the averaging of microscopic fields over the physically infinitesimal volume elements is impossible in the direction perpendicular to the 2D layer. However, by analogy with 3D materials, many authors continue to characterize the electrodynamic response of graphene and other 2D crystals, especially their nonlinear properties, by 3D quantities such as the third susceptibility χ (3) or the nonlinear refractive index n 2 . A basic inadequacy of such an approach is evident: the refractive index n characterizes the change of the phase velocity of the wave propagating inside the material, but it makes no sense to talk about the propagation of waves inside a one-atom-thick layer.
Thus a fundamental question arises, how to write down the macroscopic electrodynamic equations for atomically-thin (in fact, two-dimensional) materials and which physical quantity should be used to characterize their linear and nonlinear electrodynamic and optical properties.
In a 3D dielectric the averaging of the microscopic charge density over the physically infinitesimal volume elements leads to the definition of the polarization vector P, ρ →ρ = −div P, which has the meaning of the dipole moment per unit volume 1 and has the dimension [e/cm 2 ]. In two dimensions such an averaging can be performed only over a physically infinitesimal surface element, where the subscripts and 2 indicate 2D vectors or operators, the vector P = (P x , P y ) is the dipole moment of a surface element (the dimension [e/cm]), and we consider the polarization in the direction parallel to the layer. The corresponding contribution to the polarization current (the dimension statampere/cm) is In the presence of external charges ρ ex (r, t) and currents j ex (r, t) (which can be three-dimensional) the macroscopic Maxwell equations for a 2D nonmagnetic medium assume the form Equations (4) -(7) should be completed by a relation between the 2D polarization vector P and the electric field E(z = 0). If we ignore a possible spontaneous ferroelectric polarization of a 2D crystal, predicted in Ref. 7 , such a relation for centrosymmetric dielectric crystals has the form where the Greek indexes take only the values {x, y} and the fields are taken at z = 0. For conducting crystals it is more convenient to use the current-field relation The 2D quantities χ 2D and σ 2D are similar to the corresponding 3D ones but are measured in different units: for example, the first-and third-order 2D susceptibilities χ (1),2D and χ (3),2D are measured in cm and cm 3 /statvolt 2 , respectively, in contrast to the corresponding 3D quantities which are dimensionless (χ (1) ) and measured in cm 2 /statvolt 2 (χ (3) ). In general, the (linear and nonlinear) physical quantities χ 2D and σ 2D are complex functions of the frequency ω and the wavevector q of the electromagnetic field. They are related to each other by formulas which can be obtained from (3).
The electrodynamics of 2D materials should thus be studied using the system of equations (4) - (7), and their electrodynamic properties should be described by the 2D quantities χ 2D or σ 2D . Describing, for example, Kerr effect in graphene one should directly relate the experimentally measured quantities to the linear and nonlinear components of χ 2D or σ 2D , see, e.g., Refs. 8,9 . The use of unphysical (for 2D crystals) quantities ǫ, n and n 2 is inappropriate since they cannot be rigorously defined for 2D materials. All said above refers to any material consisting of a single or a few atomic layers, including, e.g., tilted Dirac cone 2D systems [10][11][12] and thin films of topological insulators 13 .
Let us consider some consequences of the above equations. The linear and nonlinear dynamic conductivities of conducting graphene and carbon nanotubes have been theoretically studied in many papers, [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30] . Electrodynamic properties of dielectric 2D materials have been discussed to a lesser extent. Below we discuss some properties of a model dielectric 2D crystal characterized by the linear susceptibility χ (1),2D . If a 2D crystal has a hexagonal lattice like graphene, the susceptibility χ (1),2D can be calculated within the tight-binding approximation assuming different on-site energies of electrons sitting on atoms of different sublattices. This model well describes the hexagonal boron nitride (BN) and was theoretically studied, for example, in Refs. 31,32 under the name of "gapped graphene". The spectrum of electrons and holes near the Dirac points in such a model reads where 2∆ = E gap is the energy gap and v F is the effective Fermi velocity of electrons. Assuming that the Fermi energy lies in the gap and that the temperature is small, T ≪ ∆, one can get the following analytical expression for the function χ (1),2D (ω) where g s and g v are the spin and valley degeneracies (g s g v = 4) and The frequency dependence of real and imaginary parts of the function F (Ω), Eq. (12), is shown in Figure  1. The susceptibility has a logarithmic divergence at the frequencyh|ω| = 2∆ = E gap corresponding to the inter-band transition between the valence and conduction bands. At low frequencies Ω ≪ 1 the function F (Ω) tends to unity, F (Ω → 0) = 1, so that the 2D static susceptibility is It dramatically grows when the band gap decreases. If ∆ = 1 eV (E gap = 2∆ = 2 eV), the static susceptibility (13) equals χ (1),2D ω→0 = 0.153 nm; if the gap lies in the terahertz range it is several orders of magnitude larger.
The third-order nonlinear susceptibility χ (3),2D can be calculated within the tight binding approximation, in a similar way as for graphene. 28 Within the Dirac Hamiltonian approach it was done in Ref. 32 .
If an electromagnetic wave is normally incident on a 2D dielectric layer characterized by the susceptibility (11) the transmission coefficient is determined by the formula T (ω) = |1 − 2πiωχ (1),2D (ω)/c| −2 . Its frequency dependence is shown in Figure 2. At small frequencies, ω ≪ c∆/e 2 , the function T (ω) decreases quadratically with the growing frequency, with the coefficient D ∝ 1/∆ 2 determined by the band gap. At frequencies higher that 2∆ the inter-band absorption is switched on and the transmission coefficient falls down. Notice that at the absorption edgehω ≃ 2∆ a single dielectric 2D layer (e.g., a monolayer of BN) absorbs more than 4% of the incident radiation energy which is about twice as large as in graphene (≃ 2.3%).
Its solution, is shown in Figure 3; here f (t, θ) = e −t cos θ J 0 (t sin θ), J 0 is the Bessel function, R = r/2πχ (1),2D and θ is the angle between the vector r and the z-axis. At a large distance from the charge, R ≫ 1, Eq. (16) gives the unscreened Coulomb potential φ(r) ≈ Q/r. At smaller distances R < ∼ 1 the r-dependence of the screened Coulomb potential is logarithmic, and depends on the angle θ, In the directions perpendicular (θ = 0) and parallel (θ = π/2) to the 2D plane the integrals in (18) are calculated analytically, C(0) = −γ and C(π/2) = ln 2 − γ (γ = 0.577 . . . is the Euler constant). At arbitrary angles θ the function C(θ) is shown in the Inset to Figure 3. The 2D dielectric substantially screens the field of the point charge at the distance smaller than or of order of the susceptibility length χ (1),2D which lies between ∼ 0.1 µm and ∼ 0.1 nm if the band gap varies from several meV to several eV. To summarize, we have shown that the standard procedure of averaging microscopic electromagnetic fields which is admitted in the macroscopic electrodynamics of 3D materials is inapplicable to 2D crystals like graphene and graphene related materials. We have performed a proper averaging of microscopic Maxwell equations and derived the corresponding macroscopic electrodynamic equations suitable for the description of 2D crystals. We have shown that electrodynamic properties of such crystals are adequately described by two-dimensional quantities σ 2D or χ 2D . The three-dimensional quantities such as χ (3);3D , the dielectric function ǫ, the refractive index n = √ ǫ (linear and nonlinear) cannot be properly defined and should not be used in the electrodynamics of 2D crystals.
The work has received funding from the European Unions Horizon 2020 research and innovation programme Graphene Core 2 under Grant Agreement No. 785219.