Equations of macroscopic electrodynamics for two-dimensional crystals

Equations of the macroscopic electrodynamics of bulk (three-dimensional, 3D) materials are derived from the microscopic Maxwell equations

by their averaging over “physically infinitesimal” volume elements,¹ 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 field, magnetic field, charge density, and current density. As a result of the averaging, one defines the macroscopic electric field $E=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² the higher order susceptibilities, e.g., the third order tensor $χijkl(3)$ 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,⁶ both metallic (graphene), dielectric (e.g., BN) and semiconducting (e.g., MoS₂). 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 χ⁽³⁾ or the nonlinear refractive index n₂. 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, $ρ→ρ¯=−divP$⁠, which has the meaning of the dipole moment per unit volume¹ and has the dimension [e/cm²]. 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³/statvolt², respectively, in contrast to the corresponding 3D quantities which are dimensionless (χ⁽¹⁾) and measured in cm²/statvolt² (χ⁽³⁾). In general, the (linear and nonlinear) physical quantities χ^2D and σ^2D are complex functions of the frequency ω and the wave-vector 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 Eqs. (4)–(7), and their electrodynamic properties should be described by the 2D quantities χ^2D or σ^2D. Describing, for example, the 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 and 9. The use of unphysical (for 2D crystals) quantities ϵ, n, and n₂ is inappropriate since they cannot be rigorously defined for 2D materials. All above said refers to any material consisting of a single or a few atomic layers, including, e.g., tilted Dirac cone 2D systems^10–12 and thin films of topological insulators.¹³ 

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–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 and 32, under the name “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_sg_$v$ = 4) and

The frequency dependence of real and imaginary parts of the function $FΩ$ [Eq. (12)] is shown in Fig. 1. The susceptibility has a logarithmic divergence at the frequency ℏ|ω| = 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 $χω→0(1),2D=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.²⁸ 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|⁻². Its frequency dependence is shown in Fig. 2. At small frequencies, ω ≪ cΔ/e², the function T(ω) decreases quadratically with the growing frequency,

with the coefficient D ∝ 1/Δ² determined by the bandgap. At frequencies higher than 2Δ, the inter-band absorption is switched on and the transmission coefficient falls down. Notice that at the absorption edge ℏω ≃ 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%).

Now consider the static screening of a point charge Qδ(r) = Qδ(r_∥)δ(z) placed in the plane of a 2D dielectric. The Poisson equation for the electric potential has the form

Its solution,

is shown in Fig. 3; here f(t, θ) = e^−tcosθJ₀(t sin θ), J₀ 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 of Fig. 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 bandgap 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 ϵ, and 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 Union’s Horizon 2020 research and innovation programme Graphene Core 2 under Grant Agreement No. 785219.