*Article* **Full Explicit Numerical Modeling in Time-Domain for Nonlinear Electromagnetics Simulations in Ultrafast Laser Nanostructuring**

**Enrique Moreno \*, Huu Dat Nguyen, Razvan Stoian and Jean-Philippe Colombier \***

Laboratoire Hubert Curien UMR5516, Institute of Optics Graduate School, CNRS, UJM-Saint-Etienne, University of Lyon, F-42023 St-Etienne, France; huu.dat.nguyen@univ-st-etienne.fr (H.D.N.); razvan.stoian@univ-st-etienne.fr (R.S.)

**\*** Correspondence: enrique@moreno.ws (E.M.); jean.philippe.colombier@univ-st-etienne.fr (J.-P.C.)

**Abstract:** The purpose of this paper is to present a new and accurate, fully explicit finite-difference time-domain method for modeling nonlinear electromagnetics. The approach relies on a stable algorithm based on a general vector auxiliary differential equation in order to solve the curl Maxwell's equation in a frequency-dependent and nonlinear medium. The energy conservation and stability of the presented scheme are theoretically proved. The algorithms presented here can accurately describe laser pulse interaction with metals and nonlinear dielectric media interfaces where Kerr and Raman effects, as well as multiphoton ionization and metal dispersion, occur simultaneously. The approach is finally illustrated by simulating the nonlinear propagation of an ultrafast laser pulse through a dielectric medium transiently turning to inhomogeneous metal-like states by local freeelectron plasma formation. This free carrier generation can also be localized in the dielectric region surrounding nanovoids and embedded metallic nanoparticles, and may trigger collective effects depending on the distance between them. The proposed numerical approach can also be applied to deal with full-wave electromagnetic simulations of optical guided systems where nonlinear effects play an important role and cannot be neglected.

**Keywords:** finite-difference time-domain method (FDTD); nonlinear propagation; Raman effect simulation; Kerr effect simulation; light propagation in a photoionizable media; plasma; Maxwell equations solver; laser pulse interaction; general vector auxiliary differential equation (GVADE)

#### **1. Introduction**

The finite difference time-domain (FDTD) method enables the computation and modeling of light propagation and scattering processes in linear and nonlinear dispersive materials possessing wavelength-dependent and intensity-dependent properties. This ranks the method among the most universal and powerful numerical tools for optics, electrodynamics, antennas, and waveguides theory. In this context, the general vector auxiliary differential equation (GVADE) method, which has been reported in [1] and improved in [2], is one of the most popular methods to deal with nonlinear media. However, the method suffers from some drawbacks that limit its speed and stability due to the semi-implicit updating relation for the current polarization densities. In summary, the method is not appropriate or recommended for a broad range of problems, especially those that address the designs with complex 3D topology, where iterative solutions of coupled nonlinear equations can be unstable and resource-demanding [3]. Some direct or explicit methods take into consideration the Kerr effect in the Born–Oppenheimer approximation [4]. The Kerr effect arises from the nonlinear third-order electric susceptibility. In [3], the authors take into consideration the Kerr effect in an explicit form. Other direct methods [5] require excessive computational burdens, especially for the consideration of complex geometries. When metal–dielectric interfaces are treated, especially at nanometric scales, plasmonic effects

**Citation:** Moreno, E.; Nguyen, H.D.; Stoian, R.; Colombier, J.-P. Full Explicit Numerical Modeling in Time-Domain for Nonlinear Electromagnetics Simulations in Ultrafast Laser Nanostructuring. *Appl. Sci.* **2021**, *11*, 7429. https:// doi.org/10.3390/app11167429

Academic Editor: Gennady M. Mikheev

Received: 9 July 2021 Accepted: 9 August 2021 Published: 12 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

appear [6,7] and generate high concentrations or high local densities of the electromagnetic field [8]. Since dielectrics experience nonlinear optical processes, the treatment of these regions in the interfaces is particularly challenging.

In this manuscript, we develop an algorithm that permits solving a nonlinear dielectric medium in which Kerr and Raman's effects [9] occur simultaneously. The Raman effect arises from the scattering caused by the conversion of photon energies into vibrational energy of molecules. The scattered photons have less energy than the incident photons corresponding, therefore lowering frequency fields [2]. We also assume that this takes place in an inhomogeneous dielectric environment that can potentially include metallic heterogeneities with specific optical properties. Along with the laser propagation, provided that the intensity is sufficient to ensure nonlinear absorption, the medium can also turn into a metallic state as a free electron population is generated during multiphoton absorption (MPA) or tunneling and avalanche ionization processes. It is well known that the electron transition from the valence band (VB) to the conduction band (CB) is followed by various processes, including heating in the CB, impact ionization, and various kinds of electron collisions. The purpose of this work consists of computing all these events in a full explicit form. The carrier generation rate *G*(*r*, *t*) has been computed owing to the BVkP approach implying Bloch–Volkov states [10]. The algorithms presented in this manuscript are especially suitable for simulations dealing with confinements and guiding high-energy pulsed beams. In summary, in this work, we coupled Maxwell equations and the electric charge carriers conservation equation in which the quantum ionization and recombination are governed by BVkP self-consistently.

In the next sections, by following and enriching some pioneers' works [11–13], we detail a fully explicit, stable algorithm based on GVADE FDTD, able to solve the curl Maxwell's equation in a frequency-dependent way for a nonlinear medium. In order to illustrate the combination of all these nonlinear effects occurring simultaneously in a simulation, we have modeled a Bessel beam traveling through dielectric fused silica, where some metal inclusions are modeled in a similar form as in Ref. [14]. In Section 2, we introduce the theoretical model, which describes the nonlinear effects. Section 3 presents the algorithms that discretize in time and space the main magnitudes, electromagnetic fields and electric charge. Appendix A is devoted to the stability condition assessment for these algorithms, which imposed numerical limits to the simulation performance. These limits are studied in Section 4, where we specify the validity framework considering natural limits imposed by physics. In Section 5 we present the simulation outcomes provided by the full numerical modelling and in Section 6 the conclusions are drawn. Finally, the stability conditions are derived in Appendix A and the used computational resources are presented in Appendix B.

#### **2. Modeling Electromagnetic Fields in Nonlinear Media**

The Equations (1)–(3) model the laser beam propagation by means of the electromagnetic fields *E*(*r*, *t*) and *H* (*r*, *t*). Throughout the light propagation path, fields generate and couple with the density of free charges as free electrons *nf e*<sup>−</sup> (*r*, *t*), which in turn interact with the local field. In the given form, nonlinearities are determined by the electric displacement field *D* (*r*, *t*) in the Ampère-Maxwell Equation (1). Meanwhile, the evolution of the free-electron density or medium conductivity is described by Equation (3), where *G*(*r*, *t*) is the carrier generation rate and *R*(*r*, *t*) is the carrier recombination rate.

$$\nabla \wedge \vec{H}(\vec{r},t) \quad = \begin{array}{c} \partial \vec{D}(\vec{r},t) \\ \frac{\partial}{\partial t} \end{array} \tag{1}$$

$$\nabla \wedge \vec{E}(\vec{r},t) \quad = \quad -\mu \frac{\partial \vec{H}(\vec{r},t)}{\partial t} \tag{2}$$

$$\frac{\partial n\_{f\varepsilon^{-}}(\vec{r},t)}{\partial t} = -G(\vec{r},t) - \mathcal{R}(\vec{r},t) \tag{3}$$

The electric displacement field is composed by the electric field contribution plus the polarization vectors *D* (*r*, *t*) = *ε*0*ε*∞ *E*(*r*, *t*) + ∑ *P*(*r*, *t*) [2]. The magnitude *ε*<sup>∞</sup> is the non-unity high-frequency relative permittivity of the media. On the right hand side of the displacement electric field, the sum of polarization vectors in the model that we are considering ∑ *P*(*r*, *t*) is the contribution of four terms, the metal polarization *PMetal*(*r*, *t*), the Kerr effect polarization *PKerr*(*r*, *t*), the Raman polarization *PRaman*(*r*, *t*) and the free charge polarization vector *Pf e*<sup>−</sup> (*r*, *<sup>t</sup>*) such as <sup>∑</sup> *P*(*r*, *t*) = *PMetal*(*r*, *t*) + *PKerr*(*r*, *t*) + *PRaman*(*r*, *t*) + *Pf e*<sup>−</sup> (*r*, *t*) with the detail of each term below:

$$\vec{P}\_{\text{Metal}}(\vec{r},t) \quad = \begin{array}{c} \varepsilon\_0 a \\ \hline j\omega + b \end{array} \vec{E}(\vec{r},t) \quad a, b \in \mathbb{C} \tag{4}$$

$$\vec{P}\_{\text{Kerr}}(\vec{r},t) = \epsilon \epsilon\_0 \chi\_0^{(3)} \vec{E}(\vec{r},t) \int\_{-\infty}^t \delta(t-\tau) |\vec{E}(\vec{r},\tau)|^2 d\tau \tag{5}$$

$$\vec{P}\_{\text{Raman}}(\vec{r},t) \quad = \,^t\Psi \vec{E}(\vec{r},t) \int\_{-\infty}^{t} e^{\frac{-(t-\tau)}{\tau\_2}} \sin(\frac{t-\tau}{\tau\_1}) \mathcal{U}(t-\tau) |\vec{E}(\vec{r},\tau)|^2 d\tau \tag{6}$$

$$\vec{P}\_{f\varepsilon^{-}}(\vec{r},t) \quad = \int\_{-\infty}^{t} \sigma\_{f\varepsilon^{-}}(\vec{r},t)\vec{E}(\vec{r},t)d\tau \tag{7}$$

where <sup>Ψ</sup> <sup>=</sup> *<sup>ε</sup>*0*χ*(3)(<sup>1</sup> <sup>−</sup> *<sup>α</sup>*) *<sup>τ</sup>*<sup>2</sup> <sup>1</sup> <sup>+</sup>*τ*<sup>2</sup> 2 *τ*1*τ<sup>b</sup>* <sup>2</sup> *yte*,*intent*(*in*)::*types*<sup>2</sup> is a coefficient gathering parameters depending on the considered medium. Expressions Equations (5) and (6) providing the Kerr effect polarization and Raman polarization, respectively, are taken from [9]. Here, *α* is a real-valued constant in the range *α* ∈ [0, 1] that parameterizes the relative strengths of Kerr and Raman contributions, *δ*(*t*) refers to the Dirac delta function that models the Kerr nonresonant virtual electronic transitions, which is medium dependent and in the order

of few femtoseconds for fused silica, *U*(*t*) is the Heaviside unit step function, *e* −*t <sup>τ</sup>*<sup>2</sup> sin( *<sup>t</sup> τ*1 ) models the impulse response of a single Lorentzian relaxation centered on the optical phonon frequency 1/*τ*<sup>1</sup> and having a bandwidth of 1/*τ*2. Finally *<sup>χ</sup>*(3) <sup>0</sup> is the strength of the third-order nonlinear electric susceptibility.

The upper limit of the integrations in Equations (5)–(7) extends only up to *t* because the response function must be zero for *τ* > *t* to ensure causality [3]. Moreover, freecharge polarization is introduced by using a free-charge density of current, which will be defined later. Finally, to model the metal, we consider a complex electrical permittivity (*ω*) = *<sup>a</sup> <sup>j</sup>ω*+*<sup>b</sup>* where *a* and *b* are complex numbers that fit experimental data [15]. To model the dispersive media, we assume an electric permittivity, which is frequency-dependent as the calculations are devoted to ultrafast laser propagation, with bandwidth-limited pulses that exhibit wavelength dispersion.

In ref. [2], the general vector auxiliary differential equation technique allows considering nonlinearities or dispersive characteristics of a medium. They are introduced into the Ampère–Maxwell equation by means of current densities. These currents correspond intrinsically to the temporal variation of the polarization vectors. Thereby, the following Equations (8)–(11) can be derived to correlate the currents *JMetal*(*r*, *<sup>t</sup>*) <sup>=</sup> *<sup>∂</sup> PMetal*(*r*,*t*) *<sup>∂</sup><sup>t</sup>* , *JKerr*(*r*, *<sup>t</sup>*) <sup>=</sup> *<sup>∂</sup> PKerr*(*r*,*t*) *<sup>∂</sup><sup>t</sup>* , *JRaman*(*r*, *<sup>t</sup>*) <sup>=</sup> *<sup>∂</sup> PRaman*(*r*,*t*) *<sup>∂</sup><sup>t</sup>* , *Jf e*<sup>−</sup> (*r*, *<sup>t</sup>*) <sup>=</sup> *<sup>∂</sup> Pf e*<sup>−</sup> (*r*,*t*) *<sup>∂</sup><sup>t</sup>* with the electric field:

$$\vec{f}\_{Metal}(\vec{r}\_\prime t) \quad = \quad \frac{a\varepsilon\_0}{b} \frac{\vec{E}(\vec{r}\_\prime t)}{\partial t} - \frac{1}{b} \frac{\vec{J}\_{Metal}(\vec{r}\_\prime t)}{\partial t} \tag{8}$$

$$\vec{f}\_{\rm Kerr}(\vec{r},t) \quad = \; \imath \epsilon \epsilon \wp\_0^{(3)} \frac{\partial}{\partial t} \left[ \vec{E}(\vec{r},t) \int\_0^t \delta(t-\tau) |\vec{E}(\vec{r},\tau)|^2 d\tau \right] \tag{9}$$

$$\vec{J}\_{\text{Raman}}(\vec{r},t) = \varepsilon \Psi \frac{\partial}{\partial t} \left[ \vec{E}(\vec{r},t) \int\_0^t e^{\frac{-(t-\tau)}{\tau\_2}} \sin(\frac{t-\tau}{\tau\_1}) \mathcal{U}(t-\tau) |\vec{E}(\vec{r},\tau)|^2 d\tau \right] \tag{10}$$

$$\vec{J}\_{f\varepsilon^{-}}(\vec{r},t) \quad = \quad \sigma\_{f\varepsilon^{-}}(\vec{r},t)\vec{E}(\vec{r},t) \tag{11}$$

We define a free-charge density of current as an Ohmic current, where the conductivity *σf e*<sup>−</sup> (*r*, *t*) = *qμf e*<sup>−</sup> *nf e*<sup>−</sup> (*r*, *t*) is given by the following frequency dependent expression:

$$\sigma(\vec{r},\omega) = \frac{q\mu\_{f\varepsilon^{-}}f\_d(|\vec{E}|)n\_{\text{ph}-\varepsilon^{-}}(|\vec{E}|)}{j\omega + f\_r(|\vec{E}|) - f\_a(|\vec{E}|) - f\_\eta(|\vec{E}|)}\tag{12}$$

where *q* is the electron charge, *μf e*<sup>−</sup> the free carriers of charge mobility, *nf e*<sup>−</sup> (*r*, *t*) the freecharge density and *nph*−*e*<sup>−</sup> (*r*, *<sup>t</sup>*) the photo-ionized electron density. This definition for the density of current Equation (11) is similar to Equation (9) in [10]. In Equation (13), we employ free-charge which is governed by Equation (3), or generation Equation (13) and recombination Equation (14) rates. In [10], the authors explain that electron density *nph*−*e*<sup>−</sup> (*r*, *<sup>t</sup>*) produced in the conduction band only describes the interaction of an electron with the laser electric field and accounts for the lattice periodicity, but does not account for possible collisions with phonons, ions, or other electrons. Due to these collisions, the coherence between the excited electron to the CB and its parent ion is destroyed and may introduce deviations from the above-predicted evolution of the electron density in the CB with respect to time. Indeed, without collisions, the produced electron density was shown to oscillate with time as a conduction electron may go back and forth to the VB through the action of the electromagnetic field [10,16]. Obviously, in an oscillation, the energy stays constant [17]. Hence, the authors of [10] defined a free state as a state where the electrons can be heated by inverse bremsstrahlung (IB) and are decorrelated from the parent ion. The density of charge associated with this free state is called here density of freecharge *nf e*<sup>−</sup> (*r*, *t*), and its evolution along the time is modeled by Equation (3). Nevertheless, the *μf e*<sup>−</sup> , the free carriers of charge mobility, then we employ experimental data, which allow the determination of this parameter. The method is explained in detail in Section 5.

$$G(\vec{r},t) = \begin{bmatrix} f\_d(\vec{r},t)n\_{ph-\varepsilon^-}(\vec{r},t) + n\_{f^{\varepsilon^-}}(\vec{r},t) \left[ f\_d(\vec{r},t) + f\_\eta(\vec{r},t) \right] \end{bmatrix} \tag{13}$$

$$R(\vec{r},t) \quad = \; f\_r(\vec{r},t)n\_{fc^-}(\vec{r},t) \tag{14}$$

In the carrier generation rate, the term associated with the photo-ionized electron density *nph*−*e*<sup>−</sup> (*r*, *<sup>t</sup>*) produced in the conduction band is a function of <sup>|</sup> *E*(*r*, *t*)|. From Equation (5) in [10] we know that *<sup>∂</sup> ∂t nph*−*e*<sup>−</sup> (*r*, *<sup>t</sup>*) − *<sup>N</sup>*<sup>0</sup> <sup>∑</sup> *c* |*Tcv*(*t*)| 2 = 0 and hence *nph*−*e*<sup>−</sup> (*r*, *<sup>t</sup>*) − *<sup>N</sup>*<sup>0</sup> <sup>∑</sup> *c* |*Tcv*(*t*)| <sup>2</sup> <sup>=</sup> *cte*. If this constant *cte* <sup>=</sup> 0, then *nph*−*e*<sup>−</sup> (*r*, *<sup>t</sup>*) <sup>=</sup> 0 when *E*(*r*, *t*) = 0¯. Combining this result with Equations (3) and (4) in [10], we calculate *nph*−*e*<sup>−</sup> (*r*, *<sup>t</sup>*) by Equation (15).

$$m\_{ph-c^{-}}\left(\vec{r},t\right) = \sum\_{\vec{\zeta}=:=0}^{W\_f} \left| \frac{-\sqrt{E\_{\vec{\chi}}} \int\_0^t \frac{|\vec{E}(\vec{r},\theta)|e^{-j\theta(E\_{\vec{\chi}}+\zeta\_i)}}{\left(1-\frac{j}{\Omega}\int\_0^{\theta}|\vec{E}(\vec{r},\theta)|d\phi\right)^2} d\theta \right|^2 \tag{15}$$

Finally the functions related to the generation by ionization *fd*(*r*, *t*), impact *fη*(*r*, *t*) and avalanche *fa*(*r*, *t*), as well as the function related to the recombination of the free-charge *fr*(*r*, *t*), are given by the Equations (16)–(19), respectively.

$$f\_d(\vec{r}, t) \quad = \begin{cases} \frac{1}{\tau\_d} & \text{if } |\vec{E}(\vec{r}, t)| \ge \frac{E\_{\text{max}}}{20} \\ \frac{20|\vec{E}(\vec{r}, t)|}{\tau\_d E\_{\text{max}}} & \text{if } |\vec{E}(\vec{r}, t)| < \frac{E\_{\text{max}}}{20} \end{cases} \tag{16}$$

$$\begin{array}{rcl} f\_{\eta}(\vec{r},t) & = & \begin{cases} \eta & \text{if } |\vec{E}(\vec{r},t)| \ge \frac{E\_{\text{max}}}{20} \\ \frac{20|\vec{E}(\vec{r},t)|\eta}{E\_{\text{max}}} & \text{if } |\vec{E}(\vec{r},t)| < \frac{E\_{\text{max}}}{20} \end{cases} \end{array} \tag{17}$$

$$f\_{\vec{a}}(\vec{r},t) \;= \begin{cases} \frac{\Theta(\vec{E}(\vec{r},t))^2}{\sqrt{\epsilon}} & \text{Harmonic field} \\ \Theta\left(\vec{E}(\vec{r},t) \wedge \vec{H}(\vec{r},t)\right) \cdot 2 & \text{No harmonic field} \end{cases} \tag{18}$$

$$f\_r(\vec{r}, t) \quad = \begin{cases} \frac{1}{\tau\_r} & \text{if } |\vec{E}(\vec{r}, t)| \ge \frac{E\_{\text{max}}}{20} \\ \frac{20|\vec{E}(\vec{r}, t)|}{\tau\_r E\_{\text{max}}} & \text{if } |\vec{E}(\vec{r}, t)| < \frac{E\_{\text{max}}}{20} \end{cases} \tag{19}$$

In Equations (16)–(19), we find that the factor of maximum electric field reduction, which is suggested twenty times by [10], is taken into consideration. Nevertheless, when the photoionization process is considered in a truncated region where the magnitude of | *<sup>E</sup>*(*r*, *<sup>t</sup>*)<sup>|</sup> is from the beginning higher than the threshold *Emax* <sup>20</sup> reference value, these logical conditions can be suppressed in the nested loops, and this increases the algorithm performance. In other words, this simplification is possible when we assume the description with nonlinear effects in a particular region of the computational domain. So, the electromagnetic field that arrives at that region has a magnitude higher than the threshold *Emax* <sup>20</sup> .

In Equation (18), we assume that <sup>|</sup> *E*(*r*,*t*)| 2 *<sup>ζ</sup>* is the fluence with *ζ* being the medium impedance. However, this is a particular case only valid when the electromagnetic field is a harmonic function in time. In general, and assuming a propagation along the *z*-axis, the power per surface is given by the Poynting vector.

#### **3. Algorithms**

The polarization vector due to the Kerr effect at the time *t* = Δ*<sup>t</sup> <sup>s</sup>* <sup>−</sup> <sup>1</sup> 2 can be written as *<sup>P</sup>s*<sup>−</sup> <sup>1</sup> 2 *Kerr* <sup>=</sup> *<sup>α</sup>*0*χ*(3) <sup>0</sup> <sup>|</sup> *Es*<sup>−</sup> <sup>1</sup> 2 | 2 *Es*<sup>−</sup> <sup>1</sup> <sup>2</sup> ; however, afterwards, at the instant Δ*<sup>t</sup> s* + <sup>1</sup> 2 , in order to define a full explicit approach, we have to make the approximation (*s*<sup>+</sup> <sup>1</sup> <sup>2</sup> )Δ*<sup>t</sup>* <sup>0</sup> *δ*(*t* − *<sup>τ</sup>*)| *<sup>E</sup>*(*r*, *<sup>τ</sup>*)|*d<sup>τ</sup>* (*s*<sup>−</sup> <sup>1</sup> <sup>2</sup> )Δ*<sup>t</sup>* <sup>0</sup> *<sup>δ</sup>*(*<sup>t</sup>* <sup>−</sup> *<sup>τ</sup>*)| *E*(*r*, *τ*)|*dτ*. This approximation was demonstrated in [3] and leads to:

$$\bar{P}\_{Kerr}^{s+\frac{1}{2}} \simeq \varkappa \epsilon\_0 \chi\_0^{(3)} \bar{E}^{s+\frac{1}{2}} \int\_0^{\left(s-\frac{1}{2}\right)\Lambda\_t} \delta(t-\tau) |\bar{E}(\vec{r},\tau)| d\tau$$

This approximation *<sup>P</sup>s*<sup>+</sup> <sup>1</sup> 2 *Kerr <sup>α</sup>*0*χ*(3) <sup>0</sup> <sup>|</sup> *Es*<sup>−</sup> <sup>1</sup> 2 | 2 *Es*<sup>+</sup> <sup>1</sup> <sup>2</sup> is valid under particular physical conditions explained in Section 4. To compute the nonlinear Raman effect, we make an analogous approximation. Considering that in this case there is no instantaneous cause–effect relationship, it is evident that the approximation is justified [18]. Therefore, by defining the scalar magnitude <sup>|</sup>*Q <sup>a</sup>*(*r*, *<sup>t</sup>*)<sup>|</sup> that is defined from the vector *<sup>Q</sup> <sup>a</sup>*(*r*, *<sup>t</sup>*), which could be considered as an instantaneous anisotropic Raman conductivity, we can discretize the Raman polarization vectors at the instants Δ*<sup>t</sup> <sup>s</sup>* <sup>−</sup> <sup>1</sup> 2 and Δ*<sup>t</sup> s* + <sup>1</sup> 2 as follows: *<sup>P</sup>s*<sup>−</sup> <sup>1</sup> 2 *Raman* <sup>=</sup> <sup>|</sup>*Q <sup>s</sup>*<sup>−</sup> <sup>1</sup> <sup>2</sup> *<sup>a</sup>* <sup>|</sup> *Es*<sup>−</sup> <sup>1</sup> <sup>2</sup> and *<sup>P</sup>s*<sup>+</sup> <sup>1</sup> 2 *Raman* <sup>=</sup> <sup>|</sup>*Q <sup>s</sup>*<sup>−</sup> <sup>1</sup> <sup>2</sup> *<sup>a</sup>* <sup>|</sup> *Es*<sup>+</sup> <sup>1</sup> <sup>2</sup> . Where <sup>|</sup>*Q <sup>s</sup>*<sup>+</sup> <sup>1</sup> <sup>2</sup> *<sup>a</sup>* <sup>|</sup> is updated from <sup>|</sup>*Q <sup>s</sup>*<sup>−</sup> <sup>1</sup> <sup>2</sup> *<sup>a</sup>* | by means of: <sup>|</sup>*Q <sup>s</sup>*<sup>+</sup> <sup>1</sup> <sup>2</sup> *<sup>a</sup>* <sup>|</sup> <sup>=</sup> <sup>|</sup>*Q <sup>s</sup>*<sup>−</sup> <sup>1</sup> <sup>2</sup> *<sup>a</sup>* <sup>|</sup> <sup>+</sup> *<sup>g</sup>*((*<sup>s</sup>* <sup>−</sup> *nt*(*k*))Δ*t*)| *<sup>E</sup>s*<sup>−</sup> <sup>1</sup> 2 | 2

where the Raman response function for a particular medium is modeled by *g*(*t*) = Ψ*e* − *t <sup>τ</sup>*<sup>2</sup> sin *t τ*1 *<sup>u</sup>*(*t*) and the time delay *nt*(*k*) = *nt*<sup>0</sup> <sup>−</sup> *Hc* (*k*−*k*0)Δ*znmedium <sup>c</sup>*0Δ*<sup>t</sup>* , being *nt*<sup>0</sup> = *Hc c*0Δ*<sup>t</sup> N* ∑ *i*=1 *ini*, depends on the refractive index of the media *ni* as well as the location of the electromagnetic field in the media (Figure 1). Here *u*(*t*) is the Heaviside step function:

$$u(t) = \begin{cases} 0 & \text{if } t \le 0 \\ 1 & \text{if } t > 0 \end{cases}$$

In this way, we obtain the Raman current density and the Kerr current density:

$$
\vec{J}\_{Kerr}^{\rm n} = \frac{\vec{P}\_{Kerr}^{s + \frac{1}{2}} - \vec{P}\_{Kerr}^{s - \frac{1}{2}}}{\Delta\_{\rm t}} = \frac{a\epsilon\_0 \chi\_0^{(3)} |\vec{E}^{s - \frac{1}{2}}|^2}{\Delta\_{\rm t}} \left(\vec{E}^{s + \frac{1}{2}} - \vec{E}^{s - \frac{1}{2}}\right)
$$

$$
\vec{J}\_{\rm Raman}^{\rm n} = \frac{\vec{P}\_{\rm Raman}^{s + \frac{1}{2}} - \vec{P}\_{\rm Raman}^{s - \frac{1}{2}}}{\Delta\_{\rm t}} = \frac{|\vec{Q}\_{a}^{s - \frac{1}{2}}|}{\Delta\_{\rm t}} \left(\vec{E}^{s + \frac{1}{2}} - \vec{E}^{s - \frac{1}{2}}\right)
$$

The current density updating associated to the metal is given by the equation:

$$J\_{Matral}^{\kappa + \frac{1}{2}} = g\_a J\_{Matal}^{\kappa - \frac{1}{2}} + g\_b \left( E^{\kappa + \frac{1}{2}} - E^{\kappa - \frac{1}{2}} \right),$$

where *ga* = R 2−*b*Δ*<sup>t</sup>* 2+*b*Δ*<sup>t</sup>* and *gb* = R *<sup>a</sup>* <sup>0</sup> 2+*b*Δ*<sup>t</sup>* . The operator R[] calculates the real part of these complex numbers. Finally, the electric field updating is given by the algorithm:

$$
\vec{E}^{s+\frac{1}{2}} = \vec{E}^{s-\frac{1}{2}} + \mathcal{C}\_{\mathcal{B}}^{s} \left( \nabla \wedge \vec{H}^{s} - \mathcal{g}\_{c} \vec{J}\_{\text{Metal}}^{s-\frac{1}{2}} \right),
$$

where the parameter *C<sup>s</sup> <sup>B</sup>* has different time-dependent values at the different locations:

$$\mathbf{C}\_{B}^{s} \quad = \begin{cases} \frac{\Delta\_{l}}{\epsilon\_{0}\epsilon\_{\infty} + Y\_{\overline{\xi}\_{b}^{\*}}\Delta\_{l}} & \partial\bar{\mathcal{O}}\\ \frac{\epsilon\_{0}\epsilon\_{\infty} + \nu\_{f\_{F}}^{n} + \epsilon\_{0}\chi\_{0}^{(3)}a|E^{s-\frac{1}{2}}|^{2} + |\overline{Q}\_{a}^{\*-\frac{1}{2}}|}{\epsilon\_{0}\epsilon\_{\infty} + \nu\_{f\_{F}}^{n} + \epsilon\_{0}\chi\_{0}^{(3)}a|E^{s-\frac{1}{2}}|^{2} + |\overline{Q}\_{a}^{\*-\frac{1}{2}}|} & D-\partial\bar{\mathcal{O}} \end{cases}$$

where *∂*-¯ implicates any region in metals or CPMLs, meanwhile *<sup>D</sup>* <sup>−</sup> *<sup>∂</sup>*-¯ represents any computational region filled by a dielectric medium. Finally we have two parameters, *gc* and *Υ*(*r*), that are defined as follows:

$$\begin{array}{rcl} \mathfrak{g}\_{\mathcal{C}} & = & \chi \frac{1+\mathcal{g}\_{\mathcal{C}}}{2} \\ \mathcal{Y}(\vec{r}) & = & \begin{cases} 1 & \vec{r} \in \text{Metal} \\ 0 & \vec{r} \notin \text{Metal} \end{cases} \end{array}$$

*Υ* permits to switch the algorithm from dielectric to metal and vice versa. *gc* is a coefficient gathering parameters for convenience in the metal current. Table 1 summarizes the algorithms presented in this section, and Figure 2 sketches the general flow chart. In the next section, we discuss the stability conditions for these algorithms.

**Table 1.** The list of the new explicit algorithms added to the classical FDTD scheme in order to deal with the nonlinear effects. The current list of algorithms is depicted in the form of a flow chart and it is complemented by the global integration view shown in Figure 2.


In Appendix A, we determine the stability conditions for the presented algorithms. Table 2 summarizes the stability conditions of all these effects together or separately.


**Table 2.** List of stability conditions. Note that *Min*[] selects the minimum value in a list.

Begin magnitudes


**Figure 2.** Global flow chart of the proposed GVADE-FDTD algorithm.

#### **4. Performance and Framework**

Some approximations made here are not always applicable and, in this section, we will establish the validity framework of the proposed algorithms. In particular, we should consider the two approximations (*s*<sup>+</sup> <sup>1</sup> <sup>2</sup> )Δ*<sup>t</sup>* <sup>0</sup> *<sup>δ</sup>*(*<sup>t</sup>* <sup>−</sup> *<sup>τ</sup>*)| *<sup>E</sup>*(*r*, *<sup>τ</sup>*)|*d<sup>τ</sup>* (*s*<sup>−</sup> <sup>1</sup> <sup>2</sup> )Δ*<sup>t</sup>* <sup>0</sup> *<sup>δ</sup>*(*<sup>t</sup>* <sup>−</sup> *<sup>τ</sup>*)| *E*(*r*, *τ*)|*dτ* and <sup>|</sup>*Q <sup>s</sup>*<sup>+</sup> <sup>1</sup> <sup>2</sup> *<sup>a</sup>* <sup>|</sup> <sup>=</sup> <sup>|</sup>*Q <sup>s</sup>*<sup>−</sup> <sup>1</sup> <sup>2</sup> *<sup>a</sup>* <sup>|</sup> <sup>+</sup> *<sup>g</sup>*((*<sup>s</sup>* <sup>−</sup> *nt*0(*k*))Δ*t*)| *Es*<sup>−</sup> <sup>1</sup> 2 | 2. Both expression are derived in [9] under the Born–Oppenheimer approximation [4] where the Kerr non-resonant virtual electronic transitions are in the order of 1fs, so it can be considered an instantaneous event, which is modeled by the Dirac delta function. However, while the time stepping keeps beneath 1 fs, we could model the medium reaction by using a time narrow Gaussian function with a standard deviation of 1 fs. In our simulations, Δ*<sup>t</sup>* ∼ 0.1 fs, so we are adding in a fraction of <sup>1</sup> <sup>10</sup> the standard deviation. In this particular situation, the approximation is assumed to be valid. We employ the full explicit form <sup>|</sup>*Q <sup>s</sup>*<sup>+</sup> <sup>1</sup> <sup>2</sup> *<sup>a</sup>* <sup>|</sup> <sup>=</sup> <sup>|</sup>*Q <sup>s</sup>*<sup>−</sup> <sup>1</sup> <sup>2</sup> *<sup>a</sup>* <sup>|</sup> <sup>+</sup> *<sup>g</sup>*((*<sup>s</sup>* <sup>−</sup> *nt*(*k*))Δ*t*)| *Es*<sup>−</sup> <sup>1</sup> 2 | 2 in the place of right implicit form <sup>|</sup>*Q <sup>s</sup>*<sup>+</sup> <sup>1</sup> <sup>2</sup> *<sup>a</sup>* <sup>|</sup> <sup>=</sup> <sup>|</sup>*Q <sup>s</sup>*<sup>−</sup> <sup>1</sup> <sup>2</sup> *<sup>a</sup>* | + *g*((*s* − *nt*(*k*))Δ*t*) *<sup>E</sup>s*<sup>+</sup> <sup>1</sup> <sup>2</sup> + *<sup>E</sup>s*<sup>−</sup> <sup>1</sup> 2 2 2 based 

on the same idea than before, Δ*<sup>t</sup> τ*<sup>1</sup> and Δ*<sup>t</sup> τ*<sup>2</sup> so *e* −(*t*+Δ*t*) *<sup>τ</sup>*<sup>1</sup> sin *<sup>t</sup>*+Δ*<sup>t</sup> <sup>τ</sup>*<sup>2</sup> *e* −(*t*) *<sup>τ</sup>*<sup>1</sup> sin *t τ*2 . This expression is valid for all instant ∀*t*. Hence, we fix the limit for this full explicit approach in the relation Δ*<sup>t</sup>* < 1 fs. The advantage of a fully explicit approach concerns the simulation time. Due to the updating connection in an implicit scheme, most of the algorithms should be re-updating to advance a time step. In classical FDTD, it is usually fixing the limits of the index in all nested loops. In a significant computational domain, for instance, 1.2 G cells, the updating demands 30 s. If the algorithm has to be repeated at least three times, the simulation should need around three times the the implicit scheme. Figure 3 shows the number of simulation status updates on the abscissa axis and on the ordinate axis the time needed in order to update these simulation states. Each step or update corresponds to executing the 14 algorithms shown by Figure 2. In addition to these algorithms that are necessary to perform a system status update, two numerical processes are employed. One exports simulation data or results (every 40 updates) and the other exports the complete simulation status (every 500 updates). The latter is done for two reasons. The first is to recover a state of a simulation in case of computer collapse (for example a failure in the computer electric power source). The second reason concerns the alteration of the computational domain (electromagnetic properties or/and size of the computational domain) without the need to restart the simulation from the first time stepping. In Appendix B the Table A1 lists the computational resources employed in the simulations. Under these conditions we reach a performance which can be summarized by the computational burdens and some performance parameters. Those relies on a computational domain size of 121.5 GBs in RAM, a maximum number of Yee's cell charged and computed of 1.2 Giga cells, a minimum number of Yee's cell computed of 13 Mega cells, a speed at minimum and maximum load of 6.5 Mcells/s and 38 Mcells/s respectability, a probe exportation time of 1.47 s and a computational domain backup time of 8.76 s. Therefore, Figure 3 depicts the time performance for a simulation carried out with dynamics nested loops. In our scheme, the size of the nested loops is provided by the information created. We know that our implementation violates the second rule of NASA's Ten-Rules for Developing Safety Critical Code [19]. In general, it is a good practice to arrange the loop indexes in the simulation beginning and keep these ranges. By doing this, we fix, from the start, the nested loop size. Nevertheless, the employment of a dynamic nested loop saves simulation time. Figure 3 depicts that we save more than forty hours. In a simulation with the same characteristics and static nested loop, we should have a constant duration per time step updating of 32 seconds which is equal to the asymptotic value in Figure 3. The total simulation duration is the accumulated time or area under the plot schematized in Figure 3. Hence, the triangle in gray color is the time saved by means of dynamic nested loops.

**Figure 3.** The figure shows the time per step in the simulation updating in a computational domain with dynamics limits (index k is a dynamical magnitude). Furthermore, the number of threads in the palatalization process is dynamics.

Finally, we close this section with some consideration on Equation (3). This equation is a derivation from the charge/electrons continuity equation [20], which has the general form *∂ρ*(*r*,*t*) *<sup>∂</sup><sup>t</sup>* <sup>=</sup> *<sup>q</sup>*(*G*(*r*, *<sup>t</sup>*) <sup>−</sup> *<sup>R</sup>*(*r*, *<sup>t</sup>*)) −∇·*J*(*r*, *<sup>t</sup>*), where we neglect the term ∇ ·*J*(*r*, *<sup>t</sup>*) *qnf e*<sup>−</sup> (*r*, *<sup>t</sup>*). We know from theory that *<sup>ρ</sup>*(*r*, *<sup>t</sup>*) = *qnf e*<sup>−</sup> (*r*, *<sup>t</sup>*) and *J*(*r*, *<sup>t</sup>*) = *<sup>q</sup>μf e*<sup>−</sup> *nf e*<sup>−</sup> (*r*, *<sup>t</sup>*) *E*(*r*, *t*). The key is the extremely low mobility *<sup>μ</sup>f e*<sup>−</sup> ∼ <sup>10</sup><sup>−</sup>22, which is deduced from [21]. We can consider a highly intense laser beam with <sup>|</sup> *E*(*r*, *t*)| <sup>2</sup> ∼ <sup>10</sup>13, which induces a free electron concentration in fused silica of *nf e*<sup>−</sup> (*r*, *<sup>t</sup>*) <sup>∼</sup> <sup>10</sup><sup>26</sup> and compare with ∇·*J*(*r*,*t*) *q* ∼ <sup>10</sup>−22×✄ *<sup>q</sup>*1026×1013/2 10−7 ✄ *<sup>q</sup>* <sup>∼</sup> <sup>10</sup>18. Therefore, we can conclude that *nf e*<sup>−</sup> (*r*, *<sup>t</sup>*) <sup>∼</sup> <sup>10</sup><sup>8</sup> ∇·*J*(*r*,*t*) *<sup>q</sup>* and Equation (3) is a valid approximation in the present context.

#### **5. Results**

An efficient technique to compute scattered fields in the context of FDTD modeling is the total-field/scattered-field (TF/SF) incident wave source, which is employed by most all current commercial FDTD solvers [2]. Fundamentally, the TF/SF technique is an application of the well-known electromagnetic field equivalence principle [22]. By this principle, the original incident wave of infinite extent and arbitrary propagation direction, polarization, and time-waveform is replaced by electric and magnetic current sources appropriately defined on a finite closed surface, called the Huygens surface [23], containing the object of interest. The reformulated problem confines the incident illumination to a compact total-field region and provides a finite scattered-field region external to the total-field region that is terminated by an absorbing boundary condition (ABC) to simulate the FDTD grid extending to infinity. In particular, we use convolutional perfectly matched layers, a technique already used in other works [24]. In this work, we have used the finite-difference time-domain discrete plane wave technique (FDTD-DPW), which is a numerical approach based on TF/SF technique that allows one to propagate plane waves quasi-perfectly isolated [25]. This means a propagation in the total field domain without reflections to the scattered field domain on the order of machine precision (∼300 dB) [26]. This technique is valid for any angle of propagation, for any grid cell aspect ratio and even for nonuniform grids [27]. Figure 4 illustrates the total field region where there are two important elements, an axicon lens that generates the Bessel beam from the initial plane waves and a block of fused silica. In this domain, the incident plane wave is generated by the FDTD-DPW scheme.

**Figure 4.** Computational domain.

The calculation domain contains a block of fused silica glass where the nonlinear effects apply. The interaction between the incident plane wave and the axicon lens results in a Bessel beam, which impinges the SiO2 block. We assume a pulsed laser source of wavelength (800 ± 11) nm and a time full width at half maximum of *σ<sup>t</sup>* = 40 fs. The pulsed incident time-waveform is modeled by *h*(*t*) = *A*(*t*)sech<sup>2</sup> *t*−t0 *σt* sin(*tω*), where *<sup>A</sup>*(*t*) = tanh(*tω*)+1 2 *ς* with *ς* ≥ 1, is an appropriate way to avoid numerical instabilities in the simulation beginning. The incident field six-dimensional array (IFA) that feeds the FDTD-DPW is quite simple *E*(*rIFA*,*t*) = *E*0*h*(*t*)*x*ˆ. In addition, there are other parameters that characterized the media properties and are summarize in Table 3.


**Table 3.** Media main parameters.

By using the developed irradiation strategy, the computational domain elements and their properties as well as the algorithm we carried out, the obtained results can be divided into three groups:


The power intensity *I* = *Max*(*I*(*r*, *t*)) is calculated as the maximum power that flows through a section of *σ*<sup>2</sup> *<sup>w</sup>*, being *σ<sup>w</sup>* the pulsed laser transverse waist, *σ* the pulsed length in the propagation direction and *zp* the central pulse location at the calculation time. Hence,

$$\tilde{I}(\vec{r},t) = \frac{\sigma\_w}{-\sigma\_w} \int d\mathbf{x} dy \int\_{z\_p - \sigma\_\ell}^{z\_p + \sigma\_\ell} dz \mathcal{Z} \cdot \frac{\vec{S}(\vec{r},t)}{\sigma\_w^2}.$$

To validate the algorithm, which includes nonlinear effects during light propagation, we have investigated the features of the propagation regime depending on the laser intensity. Figure 5a,b redraw the distribution evolution of the pulsed square electric field module for simulation performed for different laser intensity sources. Figure 5a plots the iso-slices that depicts <sup>|</sup> *E*| <sup>2</sup> inside the fused silica bulk for high laser power in the system. From top to bottom, the time evolution sequence reveals the defocusing process of the Bessel beam, which ends by forming some filamentation structures in the light distribution. At the time *t* = 5720, Δ*<sup>t</sup>* = 1.52 ps, the pulsed square electric field module begins to penetrate the interface air-fused silica. All successive time plots correspond to a propagation inside the fused silica block. Figure 5b illustrates, for the same time moments, but a higher initial power source, the square electric field module during the propagation inside SiO2 material.

We complement the picture with a simulation carried out in the medium power regime. From the instant of time, *t*<sup>1</sup> = 1.67 ps to the instant *t*<sup>18</sup> = 2.56 ps, Figure 5c illustrates the antagonistic trade-off between the focusing and the defocusing process experienced by the Bessel beam along with its travel through the bulk fused silica target. Figure 6b shows the history of the power density at the critical threshold separating medium and high power regimes, along with a few steps of time. In this figure we can identify the interface air-fused silica due to the rings which record the history of a stationary wave formed after the interaction. We can observe a strong focusing process in the entrance of the fused silica block followed by an alternation of refocusing and the defocusing events induced and governed by the photoionization and Kerr effect competition, and we know that these are not diffraction artifacts because in the air, at the same laser intensity, we do not observe them. At this stage, we have to give a clear definition of the pulsed power density and the power history. The magnitude of pulsed power density is defined as ∇ · *S*(*r*, *t*), where *S*(*r*, *t*) = *<sup>E</sup>*(*r*, *<sup>t</sup>*) <sup>∧</sup> *<sup>H</sup>* (*r*, *<sup>t</sup>*) is the Poynting vector [31]. Figure 6a sketches the pulsed power density. From this magnitude ∇ · *S*(*r*, *t*), we can calculate the delivered pulsed electromagnetic density of energy *<sup>t</sup>* <sup>0</sup> ∇ · *S*(*r*, *t*)*dt*. Besides this, it is interesting to consider the range of values in Figure 6a. We can identify a negative power density that requires clarification. As the Poynting vector represents the flux of power through a given surface, the divergence of this magnitude represents the power density that comes into/out a particular point in a given region. In particular, in FDTD calculations, we are describing the power flux through the computational nodes that discretized the computational domain. In the density power accounting, the reference system is placed inside the cell that encloses the nodes. In this way, a negative power density, in a particular instance of time, in a local region, means a local sink of power ∇ · *S*(*r*, *t*) < 0 in that region, and vice versa, we have a local source of power as ∇ · *S*(*r*, *t*) > 0. These oscillations in the density power magnitude indicate a redistribution that is able to support the pulsed density power and the pulsed density energy traveling.

**Figure 5.** Propagation of the pulsed beam given by a time sequence of instants which has been plotted before reaching the theoretical critical power, with an energy per pulse of 2.8 μJ allowing to reach an electron density of 1027 m<sup>−</sup>3, the beam splits, enclosing the charge and defocusing the beam. (**a**) Several instants of the pulsed square electromagnetic field when the power is beyond the critical power (High power regime). (**b**) Pulse propagation sequences during the travel inside the fused silica bulk for low power regime. (**c**) Sequence of frames that sketches the pulse along during its travel inside the fused silica bulk from the instance *t*<sup>1</sup> = 6290Δ*<sup>t</sup>* to *t*<sup>18</sup> = 9620Δ*<sup>t</sup>* (where <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> 2.65765 <sup>×</sup> <sup>10</sup>−<sup>16</sup> s and Courant number 0.69). The difference between two time-steps is 200Δ*<sup>t</sup>* in the sequence. The power of this sequence is in the threshold between low and high power regimes (critical power).

Turning to the matter at hand, we define the history of the power density by applying the absolute value to the divergence of the Poynting vector |∇ · *S*(*r*, *t*)|.

The bi-dimensional longitudinal cut of the history of the power density depicted in Figure 6b permits one to identify the air–fused silica interface. There, in that interface, we can appreciate the reflected beam and the strong focusing in the entrance due to the Kerr effect. We assume the same refractive index for both the axicon lens and the fused silica block in which the wave Bessel beam impinges.

Although we assume nonlinear media, in the simulation results presented up to this point in the manuscript, we have considered that the fused silica is isotropic and homogeneous. However, in order to study the effect of micro- and nanopores in the fused silica, 650 regions or voids have been introduced that range in size from 200 to 650 nanometers in radius. These spherical embedded cavities, that are plotted in Figure 7a, emulate a nanoporous fused silica bulk. The traditional method of linear combined congruence random number generators has been used to assign the size and location of these pores of air/void particles. The linear congruential generator yields a sequence of pseudorandomized numbers calculated with a discontinuous piecewise linear equation. The most important thing about the method is that, under the same seeds, the same random number distributions will be obtained, which allows us to reproduce both size distributions and particle locations. This will be useful to evaluate the effects resulting from another kind of media, such as metal particles, which will allow us to make a point-to-point comparison and infer conclusions from the results. In this setup, Figure 7b shows the history of the power in the region where the particles are located. The simulation is stopped when the nonlinear Bessel beam exits the porous fused silica region. Looking at Figure 7c, it appears that a correlation can be made between the distribution of photoionized electrons and the laser power density, in full consistency with our model. Figure 7c shows a history of the density of the ionized electrons as a consequence of the Bessel beam, which travels through the fused silica block. This result leads us to conjecture that both quantities look similar in terms of spatial structure.

It is relevant to study the distribution or history of power in a non-homogeneous dielectric, as is the present case. If we observe Figure 7b in detail, we see that the history of power reveals the history of field intensity as it passes through the nanoporous medium of SiO2. Due to a simple matter of impedance, the wavefront of the pulse is redistributed, avoiding entering inside the voids. As we can see in the longitudinal section shown by Figure 7b, the power vanishes inside the voids and the maximum power is not found at the interface fused silica-void. The maximum is reached in the glass medium between the pores. There is a visual explanation for this behavior. In some lenses, an anti-coating multi-layer that matches the impedance of the lens with the surrounding media (in general, glasses could be done for water or air), is coupled in order to remove the reflections. This gradient in the refractive index is necessary to facilitate the flowing energy from the lens to the environment. In the porous fused silica, there is an abrupt change in the refractive index that explains our result.

There is a Supplementary Material consisting in a video animation (see the link https://youtu.be/nvGdtxQ9E8o) that shows dynamically the ultrafast pulsed Bessel beam, depicted by employing the pulsed electric field module, that crosses the interface air-fused silica and travels through a fused silica block with air nanovoids.

**Figure 6.** (**a**) Subsample taken from the computational domain in which we plot the pulsed density of power and the history of the power density. (**b**) Power density power plotted by means of iso-slices for a pulsed Bessel beam in power regime close to the critical value. the selected time corresponds to the Bessel beam entrance inside the fused silica block. We consider all nonlinear effects in this simulation. (**c**) History of the power density for a pulsed Bessel beam. At the interface between the fused silica bulk and the air a modulation is observed in the history of power density.

(92,92,264)m (2,2,264)m (47,47,400)m

**Figure 7.** (**a**) Computational domain including clusters of spherical nanovoids, randomly located inside the fused silica bulk. (**b**) Power history of the pulsed Bessel beam propagating inside the nanoporous block. (**c**) History of the photo-ionized electrons density in the complete domain.

As stated above, we have followed identical reasoning and repeat the simulation, this time replacing nanovoids with metal particles. The results are shown in Figure 8. Figure 8a illustrates the power history under the same intensity of light. Along the photoexcitation path, the maximum power occurs around the metal– fused silica interface, enhancing and localizing the light coupling to the nanocomposite material. Even if the non-diffractive behavior is still present, it is less observable, as the beam is strongly attenuated by metal absorption. Moreover, the arrangement of the higher region of light absorption follows the randomly distributed composition of the metallic sphere. This indicates that contrary to dielectric–void interfaces, the dielectric-metal region fosters local

effects, preventing any collective response of the medium. This statement may differ if the particles are judiciously defined in terms of concentration and size in order to enhance resonance effects as plasmonic modes. In Figure 9, we observe a pattern in the electron photoionized history. In particular, Figure 9a depicts a weak structure of parallel planes, which is more intense in some regions. Figure 9b shows a stronger pattern. The present results, sketched by Figure 9a, reveal that ultrashort laser pulse interacting with distributed performed scattering nanovoids induces local field enhancement that results in enhanced nonlinear absorption and localized optical breakdown. At the void-SiO2 interface, the multiphotonic ionization process transforms the material into a thin layer of absorbing plasma with metallic-like properties. The light coupling initiates the growth of nanoplates of enhanced energy deposition, perpendicular to laser propagation direction, but parallel to the polarization. This kind of arranged structure exhibiting a periodicity approximately the wavelength of light in fused silica was commonly experimentally reported in conditions of tightly focused and high repetition rate [32,33]. They directly emerge from the spatial coherence of the waves scattered from the inhomogeneous material in the plane perpendicular to the beam propagation direction. Note the relatively used large nanovoids associated to the grid resolution do not allow to discriminate nanogratings of higher frequency perpendicular to the electric field. For propagation inside laser-induced nanoporous media, the fused silica glass remains underdense, below the critical plasma concentration, and the standing wave effect remains weak. However, the other simulation conditions, portrayed in Figure 9b), clearly indicate that for metallic nanoparticles, local field enhancement occurs near the metal-dielectric interface [34]. Upon electron–matrix energy coupling and void deformation, the light pattern may result in elongated nanostructures [35]. This shows that nanoplasmonic behavior dominates and opens the route for advanced applications as embedded micro-reflectors for a stronger near-field interaction regime.

**Figure 8.** (**a**) Power history of the pulsed Bessel beam propagating inside the fused silica block in which we place a cluster of spherical metal particles. (**b**) History of the photo-ionized electrons density in the complete domain.

Propagation and scattering through sphere inclusions with different optical properties demonstrate that the current algorithms are stable and can model efficiently nonlinear media with sharp dielectric interfaces among them. Moreover, nonlinear metal–dielectric interface investigations are far reaching for technology-based photonics. For instance, in a wider frequency spectrum, resonators, metasurfaces with designed cavities, metal– dielectric waveguides or even coating optical fibers are composed of metal–dielectric interfaces that may involve simultaneously nonlinear mechanisms of different nature. If we imagine for a sake of simplicity an harmonic electromagnetic pulse, which is guided, for instance, inside of a optical fiber, we should solve two coupled inhomogeneous Helmholtz's equations in the frequency domain <sup>∇</sup>2 *<sup>E</sup>*(*r*) <sup>+</sup> *<sup>ω</sup>*2*μ* (*r*) *<sup>E</sup>*(*r*) <sup>=</sup> −∇ *<sup>E</sup>*(*r*) · ∇(ln((*r*))) and <sup>∇</sup>2 *<sup>H</sup>*(*r*) <sup>+</sup> *<sup>ω</sup>*2*μ* (*r*) *<sup>H</sup>*(*r*)(*r*) = (∇ ∧ *H*(*r*)) ∧ ∇ ln((*r*)). In these equations, the refractive index *n*(*r*) =  *<sup>r</sup>*(*r*)*μr*, which is a dynamic magnitude along the pulsed beam propagation or simply a nonlinear magnitude, has to be taken into consideration. In this way, the algorithms present here can be used to deal with this kind of problems.

**Figure 9.** Pattern observed in the history of the photo-ionized electron history for (**a**) voids and (**b**) metal particles.

#### **6. Conclusions**

In this work, we have proposed a numerical approach for nonlinear laser beam propagation based on the time-domain discretization of Maxwell equations in a fully explicit scheme, including four nonlinear optical effects. The developed algorithms have been benchmarked in three different ways. First, we have studied the algorithm's stability conditions by combining the von Neumann method with the Routh–Hurwitz criterion. Second, we have established the physical framework for the algorithms, with special attention given to the approximation carried out in the algorithm derivation. Third, we have employed the algorithms in order to determine the transient optical properties of a photoexcited media when a laser beam is focused inside a fused silica bulk. In particular, we have considered a Bessel beam as a free-diffraction beam due to its stable and rebuildablelike energy distribution along with electromagnetic pulse travel path. This non-diffractive behavior has been observed numerically by tracking the power/energy history. This has been determined by evaluating the absolute value of the power |∇ · *S*(*r*, *t*)| and the density energy history as *<sup>t</sup>* <sup>0</sup> |∇ · *S*(*r*, *t*)|*dt*. The main advantage of the present method consists of its remarkable numerical performance (see Figure 3), that together with the tolerant stability conditions found out and the accuracy results, under an also permissive physical framework, made it valid for a wide range of applications (such as ultrafast laser processing

which has been extensively applied to fabricate 3D photonic devices based on embedded optical waveguides inscribed in glass materials by inducing permanent refractive index changes in the focal volumes [36], optical couplers and splitters [37], volume Bragg gratings [38], or even diffractive lenses [39] which have been efficiently inscribed by refractive-index changes during light propagation) in which optical/electromagnetic nonlinear effects have to be taken into consideration. Finally, our approach is generalizable to any nonlinear electromagnetic problem with dielectric–metal interfaces. The proposed work elucidates the simultaneous treatment of nonlinear effects and provides new routes toward the simulation of light–matter interaction with sub-wavelength features for the optimization of guided structures or loss-free metasurfaces.

**Supplementary Materials:** The following are available at https://youtu.be/nvGdtxQ9E8o.

**Author Contributions:** Conceptualization, E.M. and J.-P.C.; methodology, E.M.; programming, E.M.; validation, E.M.; formal analysis, E.M. and J.-P.C.; writing—original draft preparation, E.M. and J.-P.C.; writing—review and editing, all authors; visualization, E.M.; supervision, J.-P.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the FUI project 3D HYBRIDE financed by the French Single Interministerial Fund and by the LABEX MANUTECH-SISE (ANR-10-LABX-0075) of the Universite de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) and the INTRALAS project (ANR-19-CE30-0036) operated by the French National Research Agency (ANR).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data is provided in the manuscript.

**Acknowledgments:** The authors would like to thanks the free use of some tools under General Public License, that have facilitated this work: Gcc, Gfortran, Guile-OpenGL, FreeCAD, Paraview, Inkscape, SciDaVis, and Science Eclipse Parallel Tools Platform.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Stability**

In this appendix, we study the stability of the algorithm by means of the combination of the von Neumann method with the Routh–Hurwitz criterion [40]. The von Neumann method employs a Fourier series expansion of the error at the mesh nodes at a given time instant *t* = *s*Δ*t*. Due to linearity, only a single term of this expansion needs to be considered, i.e.,

$$\overrightarrow{E}^{s}|\_{i,j,k} = \overrightarrow{\overline{E}}\_{0} \overline{y}^{s} e^{j\left(i\Lambda\_{3}\overline{k}\_{3} + j\Lambda\_{3}\overline{k}\_{3} + k\Lambda\_{3}\overline{k}\_{z}\right)}\tag{A1}$$

$$\vec{D}^s|\_{i,j,k} = \widetilde{\vec{D}}\_0 y^s e^{j\left(i\Lambda\_3\overline{k}\_x + j\Lambda\_3\overline{k}\_y + k\Lambda\_4\overline{k}\_z\right)}\tag{A2}$$

where *E*<sup>0</sup> and *D*<sup>0</sup> are complex amplitudes, indexes *i*, *j*, *k* denote the position of the nodes in the mesh, Δ*x*, Δ*y*, Δ*<sup>z</sup>* are the sizes of the discretization cell, and *kx*, *ky* and *kz* are the wavenumbers of the discrete modes. The factor *y* in (A1) and (A2) is a complex variable, often called the amplification factor, which gives the growth of the error in a time updating. To ensure that a finite-difference algorithm is stable, the error must not grow as the time increases, thus, the condition |*y*| ≤ 1 must be satisfied. The stability condition for a particular FDTD scheme leads to a characteristic polynomial in *y*:

$$\Phi(y) = \sum\_{\ell=0}^{N\_p} a\_{\ell} y^{N\_p - \ell} \tag{A3}$$

Assuming that *y* are the roots of Φ(*y*), the condition for stability implicates that the roots *y* must lie inside or on the unit circle in the Y-plane. In this stability purpose, the Routh–Hurwitz criterion provides a method that ensures the roots *y* on the unit circle in the Y-plane. It establishes that the polynomial <sup>Φ</sup>(*r*) = <sup>∑</sup>*Np* =<sup>0</sup> *<sup>b</sup>yNp*− with constant real coefficients *b* has no roots in the right-half of the r-plane if all the entries of the first column of the Routh table are non-negative quantities. In [40], there is a error that makes it impossible to build the Routh table, so we explain here how to build up the Routh table. The first and second rows of the coefficients correspond to even and odd powers of r, respectively. For simplicity in the notation, and without loss of generality, it is usual assuming *Np* to be an even number, therefore *<sup>c</sup>*1, <sup>=</sup> *<sup>b</sup>*2, with = 0, 1, 2, . . . , *Np* <sup>2</sup> and *<sup>c</sup>*2, <sup>=</sup> *<sup>b</sup>*2+1, with = 0, 1, 2, . . . , *Np* <sup>2</sup> − 1. The remaining entries of the table are obtained by using the following expression *cm*, <sup>=</sup> *cm*−2,+1*cm*−1,0−*cm*−2,0*cm*−1,+<sup>1</sup> *cm*−1,0 .

The Routh–Hurwitz criterion can be used to determine if any root of the stability polynomial Φ(*y*) lies outside the unit circle in the y-plane. To this end, the bilinear transformation *<sup>y</sup>* <sup>=</sup> *<sup>r</sup>*+<sup>1</sup> *<sup>r</sup>*−<sup>1</sup> is carried out in <sup>Φ</sup>(*y*). This transformation maps the exterior of the unit circle in the y-plane onto the right-half of the r-plane. To satisfy the von Neumann stability condition for a different scheme, as a function of the parameters of interest, all the entries of the first column of the Routh table are forced to be non-negative quantities. Nevertheless, before beginning the search for stability conditions, some definitions are necessary. We define the following operators applied to Ξ*<sup>s</sup>* : ¯ *<sup>∂</sup>t*[Ξ*s*] = <sup>Ξ</sup>*s*<sup>+</sup> <sup>1</sup> <sup>2</sup> <sup>−</sup>Ξ*s*<sup>−</sup> <sup>1</sup> 2 <sup>Δ</sup>*<sup>t</sup>* , <sup>Δ</sup>¯ [Ξ*s*] = <sup>Ξ</sup>*s*<sup>+</sup> <sup>1</sup> <sup>2</sup> <sup>+</sup>Ξ*s*<sup>−</sup> <sup>1</sup> 2 <sup>2</sup> , ¯ *<sup>∂</sup>x*[Ξ*s*] = <sup>2</sup>*<sup>j</sup>* sin *kx*Δ*<sup>x</sup>* 2 Ξ*s* , ¯ *<sup>∂</sup>y*[Ξ*s*] = <sup>2</sup>*<sup>j</sup>* sin*ky*Δ*<sup>y</sup>* 2 Ξ*s* , and ¯ *<sup>∂</sup>z*[Ξ*s*] = <sup>2</sup>*<sup>j</sup>* sin *kz*Δ*<sup>z</sup>* 2 Ξ*s* , where Ξ*<sup>s</sup>* could be for instance *E<sup>s</sup>* or *D <sup>s</sup>* .

#### *Appendix A.1. Stability Condition on Metal*

For dielectric–metal interface, we deal with the wave equation *μ ∂*2*D* (*r*,*t*) *<sup>∂</sup>t*<sup>2</sup> − ∇2 *E*(*r*, *t*) = 0¯ in dielectric and the differential equation for the metal writes (*b* + *<sup>∂</sup> <sup>∂</sup>t*)*D* (*r*, *<sup>t</sup>*) <sup>=</sup> 0 *b* <sup>∞</sup> + *a* + <sup>∞</sup> *<sup>∂</sup> ∂t E*(*r*, *t*). After discretization, the wave equation is rewritten as [41],

$$
\mu \frac{\vec{\partial}\_t}{\Delta\_t^2} \vec{D}^s - \left(\frac{\vec{\partial}\_x^2}{\Delta\_x^3} + \frac{\vec{\partial}\_y^2}{\Delta\_y^3} + \frac{\vec{\partial}\_z^2}{\Delta\_z^3}\right) \vec{E}^s = \vec{0}
$$

and the auxiliary differential equation for the metal,

$$(b\vec{\Delta} + \mathfrak{d}\_{\mathfrak{l}})\vec{D}^s = \epsilon\_0 \left( (b\epsilon\_{\infty} + a)\vec{\Delta} + \epsilon\_{\infty}\mathfrak{d}\_{\mathfrak{l}} \right) \vec{E}^s.$$

Now, we replace the electric field *E<sup>s</sup>* and the electric displacement field *D <sup>s</sup>* in the last two equations by the Equations (A1) and (A2). Hence, from the wave equation we obtain:

$$(y-1)^2 \ddot{\vec{D}}\_0 + 4y \epsilon\_{\infty} \epsilon\_0 \nu^2 \ddot{\vec{E}}\_0 = 0$$

and from the differential equation for the metal:

$$(\Delta\_t b(y+1) + 2(y-1))\tilde{\vec{D}}\_0 = \Delta\_t \epsilon\_0 ((b\epsilon\_\infty + a)(y+1) + 2\epsilon\_0 \epsilon\_\infty (y-1))\tilde{\vec{E}}\_0$$

where *ν*<sup>2</sup> = <sup>√</sup> <sup>Δ</sup>*<sup>t</sup> μ* <sup>∞</sup><sup>0</sup> <sup>2</sup> *<sup>x</sup>*,*y*,*<sup>z</sup>* ∑ ℘ sin2 *<sup>k</sup>*℘Δ<sup>2</sup> ℘ 2 Δ2 ℘ . By combining both polynomials we determine the characteristic polynomial <sup>Φ</sup>(*y*) = *<sup>a</sup>*Δ*t*(*<sup>y</sup>* + <sup>1</sup>)(*<sup>y</sup>* − <sup>1</sup>)<sup>2</sup> + <sup>∞</sup> *y* <sup>4</sup>*ν*<sup>2</sup> + *<sup>y</sup>* − <sup>2</sup> + 1 (*b*Δ*t*(*y* + 1) + 2(*y* − 1)). In order to consider the Routh–Hurwitz criterion we apply the bilinear transformation *<sup>y</sup>* <sup>=</sup> *<sup>r</sup>*+<sup>1</sup> *<sup>r</sup>*−<sup>1</sup> and determine the characteristic polynomial <sup>Φ</sup>(*r*) = *<sup>b</sup>*Δ*<sup>t</sup>* <sup>∞</sup>*ν*2*r*<sup>3</sup> <sup>+</sup> 2<sup>∞</sup>*ν*2*r*<sup>2</sup> + Δ*<sup>t</sup> a* + *b* <sup>∞</sup> <sup>1</sup> − *<sup>ν</sup>*<sup>2</sup> *r* + 2<sup>∞</sup> <sup>1</sup> − *<sup>ν</sup>*<sup>2</sup> . The first column in the Routh table, which should be positive, are *c*1,0 = *b*Δ*<sup>t</sup>* <sup>∞</sup>*ν*2, *c*2,0 = 2<sup>∞</sup>*ν*2, *c*3,0 = *a*Δ*<sup>t</sup>* and *c*4,0 = 2<sup>∞</sup> <sup>1</sup> − *<sup>ν</sup>*<sup>2</sup> .

From *c*1,0 and *c*3,0 we find out that R[*a*] ≥ 0 and R[*b*] ≥ 0. The condition over *c*4,0 ≥ 0

means that 1 − *<sup>ν</sup>*<sup>2</sup> ≥ 0. This is the classic stability condition based on single-cubic FDTD [42] where <sup>Δ</sup>*<sup>t</sup>* <sup>≤</sup> *ni c*0 <sup>1</sup> Δ2 *x* + <sup>1</sup> Δ2 *y* + <sup>1</sup> Δ2 *z* , where *c*<sup>0</sup> is the speed of light in the vacuum and *ni* the refractive index of the media.

#### *Appendix A.2. Stability Condition on Kerr Effect*

In addition to the wave equation, we have to deal with the differential equation associated to the Kerr effect *<sup>∂</sup> <sup>∂</sup><sup>t</sup> <sup>D</sup>* (*r*, *<sup>t</sup>*) <sup>=</sup> <sup>0</sup><sup>∞</sup> *<sup>∂</sup> ∂t E*(*r*, *t*) + *<sup>∂</sup> ∂t PKerr*(*r*, *t*). The discretization of this equation by employing the operators results to

$$
\bar{\partial}\_t \vec{D}^s = \left( \epsilon\_0 (\epsilon\_\infty + \alpha \chi\_0^{(3)} |\vec{E}^{s-\frac{1}{2}}|^2) \right) \bar{\partial}\_t \vec{E}^s.
$$

It is evident that <sup>|</sup> *Es*<sup>−</sup> <sup>1</sup> 2 | <sup>2</sup> = *E*0 *E* ∗ <sup>0</sup> *y*2*s*−<sup>1</sup> and this result leads to the characteristic polynomial

$$\Phi(y) = (y - 1)^2 (\epsilon\_0(\epsilon\_{\infty} + a\chi\_0^{(3)} |\dot{\vec{E}}\_0|^2 y^{2s - 1})) + 4y \epsilon\_0 \epsilon\_{\infty} \nu^2$$

Comparing this characteristic polynomial with this one obtained from the wave difference differential equation Φ (*y*)=(*<sup>y</sup>* − <sup>1</sup>)<sup>2</sup><sup>0</sup><sup>∞</sup> + <sup>4</sup>*<sup>y</sup>* <sup>0</sup><sup>∞</sup>*ν*2, we see that they are similar when <sup>∞</sup> *αχ*(3) <sup>0</sup> <sup>|</sup> *E*0| <sup>2</sup>*y*2*s*−1. Obviously, we assume that |*y*| < 1, but notice here that we are not assuming the same thing that we want to demonstrate. We know that Φ (*y*), under the classical stability condition, leads to <sup>|</sup>*y*<sup>|</sup> <sup>&</sup>lt; 1 and in general for <sup>∀</sup>*<sup>s</sup>* <sup>≥</sup> <sup>1</sup> <sup>2</sup> we find that Φ(*y*) ≡ Φ (*y*) when <sup>∞</sup> *αχ*(3) <sup>0</sup> <sup>|</sup> *E*0| 2. We can identify the problem at a glance. The stability depends on the laser intensity <sup>|</sup> *E*0| through its propagation. In that case, we find the classical condition for stability already presented. Therefore, we can impose a condition to the laser intensity in order to preserve the stability; for instance, we can define an arbitrary factor *Df* such that <sup>∞</sup> *Df αχ*(3) <sup>0</sup> <sup>|</sup> *E*0| 2. Therefore, to ensure the stability | *E*0| < <sup>∞</sup> *Df αχ*(3) 0 . From the last expression, we see that stability is related to the medium, optical properties as well as laser intensity. Hence a heuristic search helps to establish an asymptotic limit for this factor *Df* ≥ <sup>10</sup>4.

#### *Appendix A.3. Stability Condition on Raman Effect*

Again, we consider the wave equation plus the below difference equation for the Raman effect,

$$\begin{split} & (1-\alpha)\epsilon\_{0}\chi\_{0}^{(3)}\frac{\mathfrak{r}\_{1}^{2}+\mathfrak{r}\_{2}^{2}}{\mathfrak{r}\_{1}\mathfrak{r}\_{2}^{2}}\Big(\frac{\partial}{\partial t}\vec{E}(\vec{r},t)\int\_{0}^{t}\overline{\mathfrak{F}}(t-\tau)|\vec{E}(\vec{r},\tau)|^{2}d\tau \\ & + \overline{\mathfrak{F}}(t)|\vec{E}(\vec{r},t)|^{2}\vec{E}(\vec{r},t)\Big)+\epsilon\_{\infty}\epsilon\_{0}\frac{\partial}{\partial t}\vec{E}(\vec{r},t)=\frac{\partial}{\partial t}\vec{D}(\vec{r},t) \end{split}$$

By applying the operators, we determine the discretized expression:

$$\begin{split} & \epsilon\_0 \ddot{\partial}\_t \dot{E}^s \left( \epsilon\_\infty + \Psi \int\_0^{s\Delta\_t} \overline{\overline{\xi}} (s\Delta\_t - \pi) |\dot{E}(\vec{r}, \tau)|^2 d\tau \right) \\ & + \epsilon\_0 (1 - \alpha) \chi\_0^{(3)} \frac{\tau\_1^2 + \tau\_2^2}{\tau\_1 \tau\_2^2} \overline{\overline{\xi}} (s\Delta\_t) \bar{\Delta} \left[ \dot{E}^s |\dot{E}^s|^2 \right] = \bar{\partial}\_t \dot{D}^s \end{split}$$

We calculate the characteristic polynomial replacing equations Equations (A1) and (A2) in the discretized wave equation and in the above expression. Then, we combine both derivation in the same formula,

$$\begin{split} \Phi(y) &= (y-1)^2 \Big( \epsilon\_{\infty} + y^{2s-1} \overline{\Psi} \left( y \int\_0^{s\Delta\_t} \overline{\overline{\xi}} (s\Delta\_t - \tau) d\tau \right) \\ &+ \frac{\Delta\_t \overline{\overline{\xi}} (s\Delta\_t)}{2} (y^2 + y + 1) \Big) |\overline{E}\_0|^2 \Big) + 4y \epsilon\_{\infty} \nu^2 \end{split}$$

Under an approach analogous to the one utilizes in the previous case we have to link the power laser to the medium. That means that the scheme should verify the inequality,

$$\epsilon\_{\infty} \gg y^{2s-1} \Psi \left( y \int\_{0}^{s\Delta\_{\ell}} \overline{\mathfrak{F}}(s\Delta\_{\ell} - \tau) d\tau + \frac{\Delta\_{t} \overline{\mathfrak{F}}(s\Delta\_{t})}{2} (y^{2} + y + 1) \right) |\widetilde{E}\_{0}|^{2} .$$

We know that *<sup>s</sup>*Δ*<sup>t</sup>* <sup>0</sup> *<sup>g</sup>*(*s*Δ*<sup>t</sup>* <sup>−</sup> *<sup>τ</sup>*)*d<sup>τ</sup>* <sup>Δ</sup>*tg*(*s*Δ*t*) <sup>2</sup> , hence, a reasoning similar to the previous section applies and we establish the stability condition <sup>|</sup> *<sup>E</sup>*0<sup>|</sup> <sup>&</sup>lt; <sup>∞</sup> *Df*(1−*α*)*χ*(3) 0 *τ*2 <sup>1</sup> <sup>+</sup>*τ*<sup>2</sup> 2 *<sup>τ</sup>*1*τ*<sup>2</sup> 2 *<sup>s</sup>*Δ*<sup>t</sup>* <sup>0</sup> *<sup>g</sup>*(*s*Δ*t*−*τ*)*d<sup>τ</sup>* . Finally for a reason of normalization we find that <sup>|</sup> *E*0| < <sup>∞</sup> *Df*(1−*α*)*χ*(3) 0 , together with the classic condition <sup>Δ</sup>*<sup>t</sup>* <sup>≤</sup> *ni c*0 <sup>1</sup> Δ2 *x* + <sup>1</sup> Δ2 *y* + <sup>1</sup> Δ2 *z* are the stability conditions in this case. Once again, we observe that stability is related to the medium optical properties as well as the laser

#### *Appendix A.4. Stability Condition on Photoionization Effect*

intensity when the Raman nonlinear effect is present.

Together the wave equation we take into consideration the differential equation *∂ <sup>∂</sup><sup>t</sup> <sup>D</sup>* (*r*, *<sup>t</sup>*) <sup>=</sup> <sup>∞</sup><sup>0</sup> *<sup>∂</sup> ∂t <sup>E</sup>*(*r*, *<sup>t</sup>*) + *<sup>σ</sup>f e*<sup>−</sup> (*r*, *<sup>t</sup>*) *E*(*r*, *t*), which by means of the defined operators can be discretized as follows,

$$
\bar{\partial}\_t \vec{D}^s = \epsilon\_\infty \epsilon\_0 \bar{\partial}\_t \vec{E}^s + \sigma\_{f\epsilon^-} \bar{\Lambda} \vec{E}^s
$$

Again we combine the equation with the discretized wave differential equation with the aim to determine the characteristic polynomial Φ(*y*) = 2<sup>0</sup><sup>∞</sup>(1 + *y* <sup>4</sup>*ν*<sup>2</sup> − <sup>2</sup> + *y*2) + Δ*tσf e*<sup>−</sup> *<sup>y</sup>*<sup>2</sup> − <sup>1</sup> in which we apply the bilinear transformation *<sup>y</sup>* <sup>=</sup> *<sup>r</sup>*+<sup>1</sup> *<sup>r</sup>*−<sup>1</sup> . Therefore, we calculate the characteristic polynomial Φ(*r*) = 2<sup>0</sup><sup>∞</sup>*ν*2*r*<sup>2</sup> + *σf e*<sup>−</sup> Δ*tr* + 2<sup>0</sup><sup>∞</sup> <sup>1</sup> − *<sup>ν</sup>*<sup>2</sup> for which the Routh-Hurwitz criterion can be applied. The first column of the Routh table are formed by the coefficients *c*1,0 = 2<sup>0</sup><sup>∞</sup>*ν*2, *c*2,0 = *σf e*<sup>−</sup> Δ*<sup>t</sup>* and *c*3,0 = 2<sup>0</sup><sup>∞</sup> <sup>1</sup> − *<sup>ν</sup>*<sup>2</sup> . The later coefficient provides the classical condition for stability. In addition, we know that *σf e*<sup>−</sup> must be a positive magnitude so as to address the stability in the algorithm updating.

#### **Appendix B. Computational Resources**

Table A1 lists the computational resources employed in the simulations. The operative system in which the simulations were carried out and runs on this hardware is Ubuntu Bionic Beaver.

**Table A1.** Computational resources.


#### **References**


## *Article* **Incident Angle Dependence of the Waveform of the Polarization-Sensitive Photoresponse in CuSe/Se Thin Film**

**Arseniy E. Fateev, Tatyana N. Mogileva, Vladimir Ya. Kogai, Konstantin G. Mikheev and Gennady M. Mikheev \***

> Institute of Mechanics, Udmurt Federal Research Center, Ural Branch of the Russian Academy of Sciences, Tatyana Baramzinoy Street, 34, 426067 Izhevsk, Russia; a.e.fateev@mail.ru (A.E.F.); mogileva@udman.ru (T.N.M.); vkogai@udman.ru (V.Y.K.); k.g.mikheev@udman.ru (K.G.M.)

**\*** Correspondence: mikheev@udman.ru; Tel.: +7-3412-21-66-11

**Abstract:** The results of studying the waveforms of longitudinal and transverse photocurrent pulses generated in thin, semitransparent CuSe/Se films as a function of the angle of incidence (α) of a femtosecond laser beam at linear and circular polarizations are presented. It has been established that the durations of unipolar longitudinal photocurrent pulses at linear and circular polarizations of laser pumping do not depend on the angle α. It is shown that the evolution of the temporal profile of the helicity-sensitive transverse photocurrent with a change in α strongly depends on polarization. At linear polarization, the shape of the unipolar pulses remains virtually constant; however, at circular polarization, the generation of unipolar and bipolar pulses is possible, with the waveforms strongly depending on the angle α. The influence of the incidence angle on the waveforms of transverse photocurrent pulses is explained by the transformation of linear and circular polarization into an elliptical upon the refraction of light at the air/semitransparent film interface and by the interplay of photocurrents arising due to linear and circular surface photogalvanic effects in the film. The presented findings can be utilized to develop polarization and incidence angle-sensitive photovoltaic devices.

**Keywords:** circular photocurrent; polarization; helicity; bipolar photovoltage; waveforms; femtosecond excitation; surface photogalvanic effect; thin films

#### **1. Introduction**

One of the interesting features of the interaction of polarized light [1–3] with semiconductors and metals is the generation of a polarization-sensitive photocurrent (PSPC) that depends on the polarization of the incident radiation according to harmonic laws [4–9]. There are longitudinal and transverse PSPCs at the oblique incidence of light on the surface of a material in which the currents flow along and perpendicular to the plane of incidence, respectively [10]. The transverse PSPC can be a combination of circular photocurrents (CPC) and linear photocurrents (LPC). The CPC depends on the direction of rotation of the electric field vector (the sign of circular polarization) of the incident radiation, while the LPC does not [11–14]. It is of interest to study the laws governing the PSPC in various media in terms of excitation of spin-polarized electrons [15,16], for the creation and development of optospintronic devices [17], laser polarization analyzers [18–20], photodiodes with high spatial resolution [21], as well as for photosensors intended for the direct recording of the polarization state of circularly polarized light [22–25]. The mechanisms of the PSPC generation include the photogalvanic effect (PGE) [26–28], the circular photogalvanic effect (CPGE) [11,29–34], the photon drag effect (PDE) [5,10,35–38], and the surface photogalvanic effect (SPGE) [5,39,40]. All of the aforementioned effects are nonlinear optical phenomena. It should be noted that the PDE and SPGE, as well as some other nonlinear optical phenomena [41,42] can be observed in any media, regardless of the type of symmetry of the medium [35]. The PDE and SPGE photocurrent pulses excited by pulsed laser pumping are

**Citation:** Fateev, A.E.; Mogileva, T.N.; Kogai, V.Y.; Mikheev, K.G.; Mikheev, G.M. Incident Angle Dependence of the Waveform of the Polarization-Sensitive Photoresponse in CuSe/Se Thin Film. *Appl. Sci.* **2022**, *12*, 6869. https://doi.org/10.3390/ app12146869

Academic Editor: Alejandro Pérez-Rodróhiez

Received: 12 May 2022 Accepted: 5 July 2022 Published: 7 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

usually unipolar (see, for example, [14,43–46]). The dependence of CPC and LPC excited in film and 2D structures due to the PDE and SPGE on the beam incidence angle α is described by an odd function [36,37,39,46–48]. Therefore, under pulsed laser pumping, the LPC and CPC unipolar pulses change their polarity when the sign of the incidence angle changes. With the simultaneous occurrence of longitudinal PDE and SPGE, due to their interplay, it is possible to generate bipolar PSPC pulses consisting of front and tail parts of opposite polarity, which change their polarities when the sign of the angle of incidence changes [40]. It was shown in our most recent work [48] that, at a given angle of incidence, the interference of the CPC and LPC generated due to the SPGE in a CuSe/Se nanocomposite film also leads to the generation of bipolar photocurrent pulses. The temporal profiles of those pulses strongly depend on the polarization ellipse of the incident femtosecond laser beam. However, despite the large number of articles on the topic of PSPC generation in various film materials (see, for example, [33,38,49–52]), studies of the influence of the incidence angle on the waveforms of photocurrent pulses arising due to the PDE and SPGE (and also CPGE) under pulsed laser pumping have not been carried out yet.

In this work, using the study of the SPGE in a CuSe/Se semitransparent thin film under femtosecond laser excitation as an example, it is shown for the first time that the waveforms of the PSPC pulses can strongly depend on the angle of incidence. In particular, it has been shown that with circularly polarized pulsed pumping and a fixed sign of the incidence angle, unipolar transverse photocurrent pulses of opposite polarity are generated at small and large angles of incidence, and bipolar pulses are excited in the intermediate range of incident angles, smoothly transforming into unipolar pulses of opposite polarity at the boundaries of this interval. The results obtained are explained by the change in the state of polarization in the refracted beam and by the interference of the LPC and CPC that arises in the subsurface layer of the semiconducting medium under study.

#### **2. Materials and Methods**

Thin CuSe/Se films with a thickness of 130 nm and dimensions of 15 × 35 mm were synthesized on a glass substrate by successive vacuum thermal deposition of Se and Cu, according to the procedure described in one of our previous publications [14]. To stabilize the residual Se [53,54], the synthesized film structure was annealed at a temperature of 140 ◦C for 30 min. As a result, the selenium from an unstable amorphous phase was transformed into a stable polycrystalline trigonal selenium (*t*-Se). It should be noted that currently there are some other methods of film copper selenides structure synthesis that have also been developed [55–60].

The phase composition of the films was studied at room temperature on a Bruker D2 PHASER diffractometer. We used CuKα radiation (λ = 0.154 nm). Using the DIFFRAC.EVA universal program, the diffraction curves were smoothed and the background due to X-ray scattering on a glass substrate was subtracted. To study the phase composition of the films, we additionally used a HORIBA HR800 Raman spectrometer with laser excitation at 632.8 nm at the radiation intensity of 14 kW/cm2. Sample surface images were taken using a scanning electron microscope (SEM) (Thermo Fisher Scientific Quattro S, Brno-Cernovice, ˇ Czech Republic). The elemental composition of the synthesized films was determined using an energy-dispersive microanalyzer based on an EDAX "Octane Elect Plus EDS System" spectrometer built into the SEM. The optical transmittance spectra of the films were recorded using a two-beam spectrophotometer (PerkinElmer Lambda 650, Shelton, WA, USA). The film thickness was measured using a stylus profilometer. The refractive index *n* and the extinction κ of the film were determined using an ellipsometer.

The X-ray diffraction pattern (see Figure S1, Supporting Information) and the Raman spectrum (see Figure S2, Supporting Information) demonstrate that the film has two phases and consists of CuSe and *t*-Se. The mass percentages of CuSe and *t*-Se in the synthesized film are 33.3 and 66.7 wt.%, respectively. The film surface consists of a number of flat petallike structures, which are CuSe crystallites predominantly oriented in radial directions

from their centers (see Figure S3, Supporting Information). The centers of these structures are located at distances of about 15–30 μm from each other. The film is semitransparent in the visible wavelength range (see Figure S4, Supporting Information). The optical transmittance, refractive index *n*, and extinction coefficient κ at 795 nm are 27.6%, 1.64, and 2.14, respectively. The synthesized film is electrically conductive with a sheet resistance of <sup>39</sup> <sup>Ω</sup>/-, which is orders of magnitude lower than the corresponding value of many known thin-film semiconductor structures (see, for example, [61,62]). It should be added that the sheet resistance of CuSe thin films with thicknesses from 50 to 500 nm synthesized by the chemical bath method varies between 23–50 <sup>Ω</sup>/- [63]. Meanwhile, the sheet resistance of 280 nm thick bismuth selenide (Bi2Se3) film obtained by the same method immediately after synthesis is 1012 <sup>Ω</sup>/-, and after its annealing in air or in the nitrogen medium, the sheet resistance decreases to 3.6 <sup>×</sup> 103 <sup>Ω</sup>/-[64].

To measure the laser-induced photocurrent (photovoltage), two film gold electrodes A and B were deposited on opposite short sides of the synthesized film by vacuum thermal deposition. The film under study was placed on a special goniometric device, which made it possible to smoothly change its spatial orientation relative to the incident laser beam. The surface of the film remaining after depositing the electrodes was sufficient to study the generation of photocurrent under the action of a narrow laser beam 1.5 mm in diameter in a wide range of incidence angles.

The photocurrent pulses in a CuSe/*t*-Se nanocomposite film were excited by pulses of a Ti:S femtosecond laser at a wavelength of 795 nm (pulse duration 120 fs, pulse repetition rate 1 kHz) at an oblique incidence of radiation on the film. All the experiments were carried out at laser pulse energy *E*in = 100 μJ. The polarization of the radiation incident on the film was controlled using half-wave and quarter-wave plates. The photocurrent pulses were registered and recorded using a digital oscilloscope with a bandwidth of 400 MHz and a rise time of 875 ps with averaging over multiple exciting laser pulses. The input impedance of the oscilloscope was *R*in = 50 Ω. In the experiments, the extreme values of the pulses of the longitudinal and transverse photocurrents *j*x and *j*y, respectively, were recorded, determined by the formulas *j*<sup>x</sup> = *U*x/*R*in and *j*<sup>y</sup> = *U*y/*R*in, where *U*<sup>x</sup> and *U*<sup>y</sup> are the extreme values of the voltage pulses at the experimental geometries when the measuring electrodes were placed perpendicular and parallel to the plane of incidence, respectively. Since *U*<sup>x</sup> and *U*<sup>y</sup> depend linearly on *E*in [65], it was convenient to introduce the conversion coefficients η<sup>x</sup> = *j*x/*E*in and η<sup>y</sup> = *j*y/*E*in, which characterize the efficiency of converting light into longitudinal and transverse photocurrents, respectively. It should be noted that when studying the waveforms of the photocurrent pulses, we subtracted the time-varying background signal due to noticeable electrical interference coming from the power supply unit of the femtosecond laser.

#### **3. Results and Discussion**

Figure 1a shows the dependences of the coefficients of light conversion into longitudinal photocurrent η<sup>x</sup> on the angle of incidence α for linear and circular polarizations, measured according to the sketch of experimental arrangements presented in the upper frame. It can be seen that the photocurrent is absent at normal incidence of the radiation on the film (α = 0). The photocurrent changes polarity when the sign of the angle of incidence changes, increases in absolute value with an increase in the modulus of α, takes its extreme values at angles |α| ≈ 66–68◦, and disappears at a grazing incidence of radiation on the film. Thus, the longitudinal photocurrent for linear and circular polarizations is described by an odd function of the angle α:

$$
\mathfrak{η}\_{\mathbf{x}}(\mathfrak{a}) = -\mathfrak{η}\_{\mathbf{x}}(-\mathfrak{a}).\tag{1}
$$

**Figure 1.** The light to longitudinal photocurrent conversion coefficient η<sup>x</sup> as a function of (**a**) the light incidence angle α on the film at linear (*p*-polarization) and circular polarizations (points denote experimental results, curves denote fitting with η<sup>x</sup> = 35.37 sin 2α/(1.76 cos α + 1) <sup>2</sup> and η<sup>x</sup> = 15.68 sin 2α/(1.75 cos α + 1) <sup>2</sup> for linear and circular polarization, respectively) and (**b**) the azimuth of the polarization plane Φ at α = 43.5◦ (dots denote experimental results, curve denotes fitting with ηx(Φ) = 6.7 cos2 Φ; orientations of the polarization plane are shown at the top). The insets show the sketches of the experimental setups: **k** and **E** are the wave vector and the electric field vector of the incident radiation, respectively; σ is the plane of incidence; α is an angle of incidence; **n** is the normal to the film surface; **A** and **B** are the measuring electrodes; x, y are the axes of the rectangular coordinate system; the x and x axes lie in the σ plane (x ⊥ **k**, the x axis is perpendicular to the electrodes **A** and **B**); Φ is the angle between the x and **E**.

With linear polarization, the dependence of η<sup>x</sup> on the polarization azimuth Φ (the angle between the plane of polarization and the radiation incidence plane on the film σ) at a fixed α is described by an even function:

$$
\eta\_{\rm tx}(\mathfrak{a}, \Phi) = \eta\_{\rm tx}(\mathfrak{a}, \Phi = 0) \cos^2 \,, \tag{2}
$$

where ηx(α, Φ = 0) is the conversion coefficient at a given α and Φ = 0. It should be noted that Φ = 0 corresponds to *p*-polarized radiation. All this is clearly seen from the dependence ηx(Φ) = ηx(α, Φ = 0) cos2 Φ obtained at α = 43.5◦, where ηx(α, Φ = 0) = 6.7 mA/mJ (see Figure 1b). From Equation (2) it follows that ηx(α = 43.5◦, Φ = 90◦) = 0, i.e., no photocurrent is generated at *s*-polarization. Experiments have shown that this equation is valid for any α. Figure 2a shows the dependence of the conversion coefficient of light into transverse photocurrent η<sup>y</sup> on the angle of incidence α for linearly polarized radiation at polarization azimuth Φ = −45◦, obtained using the experimental setup shown in the upper inset to this figure. It can be seen that, similarly to the angular dependence of the longitudinal photocurrent, the experimental dependence ηy(α) is described by an odd function, i.e.:

$$
\eta\_\mathbf{y}(\mathfrak{a}) = -\eta\_\mathbf{y}(-\mathfrak{a}).\tag{3}
$$

**Figure 2.** The light to transverse photocurrent conversion coefficient η<sup>y</sup> as a function of (**a**) the incidence angle α at the polarization azimuth Φ = −45◦ (triangles correspond to the experimental data while the solid line represents the result of the fitting with equation η<sup>y</sup> = 8.37 sin 2α/(0.87 cos α + 1) 2 ) and (**b**) the polarization azimuth Φ at α = 43.5◦ (circles denote experimental data, the solid curve represents the fitting with the equation η<sup>y</sup> = −3.1 sin 2 Φ; the orientations of the plane of polarization of the incident radiation are shown at the top). The insets show the sketches of the experimental setups.

Meanwhile, the measured dependence of the transverse photocurrent on the polarization azimuth Φ at a fixed α = 43.5◦ (see Figure 2b) is approximated by an odd function:

$$
\eta\_\text{y} \left( \mathfrak{a} = 43.5^\circ, \Phi \right) = \eta\_\text{y} \left( \mathfrak{a} = 43.5^\circ, \Phi = 45^\circ \right) \sin 2\Phi,\tag{4}
$$

where ηy(α = 43.5◦, Φ = 45◦) is the conversion coefficient of light into transverse photocurrent at α = 43.5◦ and Φ = 45◦.

The set of angular and polarization dependences described by Equations (1)–(4) shown in Figures 1 and 2 indicates that the generation of photocurrent in the studied CuSe/*t*-Se thin films occurs according to the SPGE mechanism [5,39,40,65].

Figure 3a shows the unipolar pulse waveforms of the longitudinal photovoltage (measured in the geometry of the experiment when the measuring electrodes are oriented perpendicular to the plane of incidence), normalized to their extreme values at *p*-polarization of the incident radiation at α = ±45◦ and α = ±78◦. It follows that when the sign of the incidence angle is changed, the pulses are inverted, which is in agreement with Equation (1). It can also be seen that the waveforms of the pulses are virtually independent of the angle of incidence. This is evidenced by the dependence curve of the pulse duration τ (FWHM) on α, presented in the inset to Figure 3a. Similar patterns were obtained when the longitudinal photovoltage was excited by circularly polarized radiation pulses (see Figure 3b). Thus, the waveforms of the longitudinal photovoltage pulses in the films under study for linear and circular polarizations are virtually independent from the angle of incidence.

**Figure 3.** The waveforms of the longitudinal photovoltage pulses normalized to their extreme values at (**a**) linear polarization (Φ = 0, *p*-polarization) and (**b**) circular polarization recorded at positive (solid lines) and negative (dashed lines) incidence angles α. The upper insets show the corresponding dependences of the recorded pulse durations τ on the angle of incidence. The bottom insets show the sketches of the experimental setups.

The waveforms of the transverse photovoltage pulses recorded at linear polarization (Φ = −45◦) and normalized to their extreme values at different α are shown in Figure 4. It can be seen that the waveforms of these pulses weakly depend on α. This is evidenced by the dependence of the pulse duration on the incidence angle, shown in the inset in the same figure.

**Figure 4.** The waveforms of the transverse photovoltage pulses normalized to their extreme values at linear polarization (Φ = −45◦), recorded at positive (solid lines) and negative (dashed lines) angles of incidence α. The upper inset shows the dependence of the recorded pulses duration τ on the angle of incidence. The bottom inset shows the sketch of the experimental setup.

Figure 5 shows the waveforms of the transverse photovoltage pulses recorded for the right-hand circularly polarized laser beam at various angles of incidence. The pulses are normalized to their extreme values. It can be seen that at 0 < α < 58.5◦, positive pulses are generated, the duration of which gradually increases with increasing α. In the range of angles α approximately 58.5 ≤ α ≤ 76.5◦ bipolar pulses are excited. Figure 5 shows that at α = 58.5◦, a small negative pulse appears on the leading edge of the positive pulse. The amplitude of the front negative pulse increases with increasing the angle of incidence, while the amplitude of the positive tail of the pulse decreases. At large angles of incidence (α = 79.5◦ and 82.5◦), the positive part of the pulse completely disappears and the photovoltage is generated in the form of a negative unipolar pulse. It should be added that when the direction of rotation of the electric field vector of the incident radiation changes (i.e., at left-hand circularly polarized laser beam) the polarities of the pulses shown in Figure 5 change to the opposite ones. Thus, the waveforms of the transverse photovoltage pulses at circular polarization strongly depend on the angle of incidence, which is not typical of the longitudinal photovoltage pulses at different polarizations (see Figure 3) and the transverse photovoltage pulses at linear polarization (see Figure 4).

The unusual dependence of the waveform of the transverse photovoltage pulses on the angle of incidence at circular polarization can be explained on the basis of the interaction of linear and circular photocurrents arising in the surface layer of the film as a result of the transformation of the incident circular polarization into elliptical polarization upon the refraction of light at the air/film interface.

**Figure 5.** The waveforms of the transverse photovoltage pulses normalized to their extreme values for circularly polarized incident radiation, recorded for different angles of incidence α. The inset shows the sketch of the experimental setup.

To determine the polarization of the refracted beam, one can use the complex Fresnel refractive indices *t*p, *t*s for *p*- and *s*-polarizations, respectively, at the air/semitransparent film interface given in [66]:

$$t\_{\rm P} = \frac{2\hat{\mu}\cos\alpha}{\hbar^2\cos\alpha + \sqrt{\hbar^2 - \sin^2\alpha}}, t\_{\rm s} = \frac{2\cos\alpha}{\sqrt{\hbar^2 - \sin^2\alpha} + \cos\alpha},\tag{5}$$

where *n*ˆ = *n* + *i*κ is the complex refractive index, *n* is the real refractive index that determines the phase velocity (*n* = 1.64), and κ is the absorption coefficient (κ = 2.14) that determines the attenuation of light in the film itself. From expression (5), the phase shift δ<sup>t</sup> between the components of the refracted beam with *p*- and *s*-polarizations is determined by the formula:

$$
\delta\_{\rm t} = -\delta\_{\rm tp} + \delta\_{\rm ts} + \delta\_{\rm 0} \tag{6}
$$

where:

$$\delta\_{\rm tp} = \text{Arg}(t\_{\rm P}) = \tan^{-1}\left(\frac{-\kappa^2 \cos \alpha + \chi \cos \xi + n^2 \cos \alpha}{\chi \sin \xi + 2\kappa n \cos \alpha}\right) - \tan^{-1}\left(\frac{n}{\kappa}\right),\tag{7}$$

$$\delta\_{\rm fb} = \mathrm{Arg}(t\_{\rm s}) = \tan^{-1}\left(\frac{(\cos\alpha + \chi\cos\xi)\csc\xi}{\chi}\right) - \frac{\pi}{2},\tag{8}$$

$$\chi = \sqrt[4]{\left(-\sin^2\alpha - \kappa^2 + n^2\right)^2 + 4\kappa^2 n^2} \tag{9}$$

$$\mathcal{K} = \frac{1}{2\left(\frac{\pi}{2} - \tan^{-1}\left(\frac{-\sin^2\alpha - \kappa^2 + n^2}{2\kappa n}\right)\right)'} \tag{10}$$

δtp, δts are the phase shifts of the *p*- and *s*-beam components, respectively, resulting from refraction at the interface between two media with a complex refractive index, and δ<sup>0</sup> is the initial phase shift between the *p*- and *s*-components before refraction. For the circular polarization of the incident radiation δ<sup>0</sup> = π/2, where the signs "+" and "−" are the opposite of the sign of circular polarization.

It follows from Equation (5) that the amplitude values of the complex transmission coefficients *t*p and *t*s can be determined using:

$$\left| \mathbf{t}\_{\mathbb{P}} \right| = \frac{2 \cos \alpha \sqrt{\kappa^2 + n^2}}{\sqrt{\left(-\kappa^2 \cos \alpha + \chi \cos \xi + n^2 \cos \alpha\right)^2 + \left(2 \kappa n \cos \alpha + \chi \sin \chi\right)^2}},\tag{11}$$

$$|t\_s| = \frac{2\cos\alpha}{\sqrt{\chi^2\sin^2\xi + \left(\cos\alpha + \chi\cos\xi\right)^2}}\,\tag{12}$$

The amplitudes of the *p*- and *s*-components of the electric field vector of the refracted beam **E**(**t**) describing the ellipse in the x"y coordinate system can be found using the equations:

$$T\_{\mathbf{x''}} = \begin{vmatrix} t\_{\mathbf{P}} \end{vmatrix} E\_{\mathbf{P}}; \ T\_{\mathbf{Y}} = |t\_{\mathbf{s}}| E\_{\mathbf{s}}. \tag{13}$$

respectively, where x"y is the coordinate plane perpendicular to the wave vector **k**t of the refracted beam, axis x" is in the plane of incidence σ, and the axis y is perpendicular to σ.

In this case, the equation of the polarization ellipse in the refracted beam depending on the ratio of the *p*- and *s*-components of the vector **E**(**t**) and the phase shift δ*<sup>t</sup>* can be written as follows:

$$\frac{\left(\mathbf{x''}\right)^2}{\left(T\_{\mathbf{x''}}\right)^2} + \frac{\mathbf{y}^2}{T\_{\mathbf{y}}^2} - \frac{2\mathbf{x''}\mathbf{y}\cos\delta\_{\mathbf{t}}}{T\_{\mathbf{x''}}T\_{\mathbf{y}}} = \sin^2\delta\_{\mathbf{t}}.\tag{14}$$

The polarization ellipses are characterized by the angle ψ between the semi-major axis *a* of the ellipse and the x" axis lying in the incidence plane σ of radiation on the film (see Figure 6a,b right insets), as well as the degree of circular polarization *P*cir = γ2*ab*/ *a*<sup>2</sup> + *b*<sup>2</sup> and degree of linear polarization *P*lin = *<sup>a</sup>*<sup>2</sup> − *<sup>b</sup>*<sup>2</sup> / *a*<sup>2</sup> + *b*<sup>2</sup> , where *b* is the semi-minor axis of the ellipse, γ is the sign of circular polarization, and γ = 1 and γ = −1 denote the rotation of the electric field vector to the right and to the left, respectively.

Expression (14) makes it possible to calculate the angle ψ, as well as the length of the minor *b* and major *a* semiaxes of the refracted beam polarization ellipse using the following Equations:

$$\psi = \frac{1}{2} \tan^{-1} \left( \frac{2T\_{\mathbf{x}''}T\_{\mathbf{y}} \cos \delta\_{\mathbf{t}}}{T\_{\mathbf{x}''} ^2 - T\_{\mathbf{y}} ^2} \right), \tag{15}$$

$$a = \sqrt{\frac{1}{2} \left( T\_{\mathbf{x''}} ^2 + T\_{\mathbf{y}} ^2 + \sqrt{T\_{\mathbf{x''}} ^4 + T\_{\mathbf{y}} ^4 + 2T\_{\mathbf{x''}} ^2 T\_{\mathbf{y}} ^2 \cos 2\delta\_{\mathbf{t}}} \right)},\tag{16}$$

$$b = \sqrt{\frac{1}{2} \left( T\_{\mathbf{x}''} ^2 + T\_{\mathbf{y}} ^2 - \sqrt{T\_{\mathbf{x}''} ^4 + T\_{\mathbf{y}} ^4 + 2T\_{\mathbf{x}''} ^2 T\_{\mathbf{y}} ^2 \cos 2\delta\_{\mathbf{t}}} \right)}. \tag{17}$$

Equations (16) and (17) were used to calculate *P*cir and *P*lin. Figure 6a,b shows the calculated values of ψ, *P*cir, and *P*lin for different angles α, for which the waveforms of the photovoltage pulses were recorded for circular and linear polarizations (Φ = −45◦) of the incident radiation, respectively. It can be seen that when the incident beam is circularly polarized at large angles α, there is a significant change in the parameters ψ, *P*cir, and *P*lin characterizing the polarization ellipse in the refracted beam. However, at α < 20◦, the refracted beam remains virtually circularly polarized, i.e., at small angles of incidence, the polarizations of the incident and refracted beams

virtually coincide. It should be added that at all α, the signs of the circular polarization of the incident and refracted beams coincide, i.e., for both beams γ = 1.

**Figure 6.** The angle ψ of the semi-major axis *a* of the refracted beam polarization ellipse relative to the plane of incidence σ, as well as the degrees of circular *P*cir and linear *P*lin polarizations of the refracted beam (insets) as a function of the angle of incidence α for (**a**) circularly and (**b**) linearly polarized (Φ = −45◦) radiation incident on a thin CuSe/*t*-Se film (points are calculated values of ψ, *P*cir and *P*lin for angles α, for which photovoltage pulses were recorded, curves denote smoothing functions). The graphical insets show the corresponding calculated polarization ellipses of the refracted beam at α = 64.5◦ in the xy coordinate system, where the x axis lies in the refraction plane coinciding with the plane of incidence σ and is perpendicular to the wave vector of the refracted beam.

It follows from Figure 6b that with linear polarization of the incident light at Φ = −45◦ the refracted beam is also elliptically polarized, but in this case γ = −1 (see upper inset). In addition, when changing α in the range of 0–85◦, the modulus *P*cir of the refracted beam changes insignificantly from 0 to 0.32. At small angles α, the changes in the polarization state of the refracted beam are minimal.

It is known that the transverse photocurrent *jy* for elliptical excitation beam polarization consists of a CPC (*j*y,cir ) and an LPC (*j*y,lin ), i.e., *j*<sup>y</sup> = *j*y,cir + *j*y,lin . According to [48] and taking into account [29], for positive angles of incidence, the transverse photocurrent amplitude in the CuSe/Se film structure for elliptical beam polarization can be represented as follows:

$$j\_\mathbf{y} = A\_{\mathrm{circ}} P\_{\mathrm{circ}} - A\_{\mathrm{lin}} P\_{\mathrm{lin}} \sin 2\psi\_\prime \tag{18}$$

where *A*cir and *A*lin are the positive coefficients of the circular and linear photocurrents at a given α, respectively. Taking into account Equation (18) and in accordance with the results of [37], the waveforms of the transverse photovoltage pulses *U*y(*t*,α) can be written as follows:

$$dI\_\mathbf{y}(t,\alpha) = B\_{\mathrm{circ}}(\alpha)P\_{\mathrm{circ}}f\_{\mathrm{circ}}(t) - B\_{\mathrm{lin}}(\alpha)P\_{\mathrm{lin}}\sin 2\psi \, f\_{\mathrm{lin}}(t),\tag{19}$$

where *B*cir(α), *B*lin(α) are positive coefficients (for positive α) depending on α and characterizing the circular and linear contributions, respectively, and *f*cir(*t*) and *f*lin(*t*) are the transverse photovoltage waveforms normalized to the maximum values, recorded at angles α close to zero for circular and linear polarizations, respectively. Equation (19) allows one to approximate the waveforms of the photovoltage pulses, recorded at different α, with two unknown parameters *B*cir(α) and *B*lin(α). For example, Figure 7 shows the approximations of the waveforms for four of the recorded pulses at different α. It can be seen that the obtained approximating curves satisfactorily describe the experiment. It should be noted that if the waveforms *f*cir(*t*) and *f*lin(*t*) coincide with each other, the generation of bipolar pulses is impossible, and the duration of unipolar pulses of the resulting photocurrent, defined by Equation (19), does not depend on the angle of incidence.

Figure 8 shows the dependences of the calculated coefficients *B*cir and *B*lin on α. It follows from the figure that both coefficients *B*cir and *B*lin approach zero at small and also at grazing angles of incidence. However, the dependences of *B*cir and *B*lin on α differ significantly from each other. It can be seen that the dependences of *B*cir and *B*lin acquire their extreme values at α = 54◦ and α = 62◦, respectively. The coefficient *B*lin prevails over *B*cir for the entire range of the incident angle change and the ratio *B*lin/*B*cir increases monotonically with increasing α.

**Figure 7.** The oscillograms of the photovoltage pulses arising in a thin CuSe/*t*-Se film with circular polarization of exciting radiation and angles of incidence α = 10.5, 58.5, 64.5, and 82.5◦ (black circles), and their approximations shown by red lines according to Equation (19).

It follows from Figure 6b that with linear polarization (Φ = −45◦) of the incident beam, the refracted beam becomes elliptically polarized at a negative polarization sign (γ = −1). This means that both terms on the right side of Equation (19) remain negative for any positive angle of incidence. In addition, when changing α in the range of 0−85◦, the parameters *P*cir and *P*lin of the refracted beam do not change significantly. Taking into account Equation (19), this leads to a weak dependence of the waveform of the photovoltage

pulse on the angle α for linearly polarized laser pumping. It should be noted that the transverse photocurrent arising in the medium due to the PDE or CPGE also consists of the circular and linear contributions. Therefore, the temporal profile of the transverse photocurrent pulses of the PDE or CPGE generated in the nonlinear optical medium when pumped by short laser pulses of circular polarization can also depend in a complex way on the angle at which the light falls on the surface of the absorbing medium.

**Figure 8.** The coefficients *B*lin (green crosses), *B*cir (blue dots) as a function of the angle of incidence α, which characterize the linear and circular contributions to the transverse photovoltage pulses, calculated from the experimental data (the green and blue curves represent the corresponding approximations by the equations *B*lin = 55.98 sin 2α/(1.15cos α + 1) 2 , *B*cir = 11.75 sin 2α /(0.5cos α + 1) 2 , respectively), and calculated dependence of the ratio *B*lin/*B*cir on α (pink curve).

Good agreement between the experimental data and the calculated dependences confirms that the angular dependence of the waveforms of transverse photocurrent pulses in a thin CuSe/Se film under circularly polarized femtosecond laser pumping originates from the transformation of circular polarization into elliptical polarization upon refraction of light in a semitransparent CuSe/Se film and the interaction of LPC and CPC having different relaxation times in the film structure. As mentioned above, despite the large number of publications on the topic of generation of polarization-sensitive transverse photocurrent in various materials, arising by various mechanisms, such a phenomenon has not been observed before. It is possible that this is due to the fact that in many studies of the transverse photocurrent, cw laser radiation or nanosecond laser radiation has been used (see, for example, [33,67–69]).

The results obtained in this work can be used in various applications in optoelectronics. For example, the experimental setup presented in the inset to Figure 5 can be used for the fast direct detection of the circular polarization state of light. If at angles of incidence 0 < α < 58.5◦ (for example, at α = 45◦) the photocurrent pulses generated in a CuSe/*t*-Se thin film have a positive polarity, then this means that the incident radiation is right-hand polarized (looking towards the light source). If, under the same experimental conditions, the photocurrent pulses have a negative polarity, then the incident radiation is left-hand polarized. This method of determining the state of circular polarization does not require the use of optical elements. Further, it is obvious that this photovoltaic property of a CuSe/Se thin film can be used to determine the fast and slow axes of a quarter-wave plate without using a reference quarter-wave plate and an optical light polarization.

#### **4. Conclusions**

In thin semitransparent CuSe/*t*-Se films synthesized by vacuum thermal deposition, the generation of nanosecond photocurrent pulses is studied as a function of the angle of incidence and polarization of exciting femtosecond laser pulses at 795 nm. It has been established that the dependences of the longitudinal and transverse photocurrents on the angle of incidence are described by odd relationships, which are characteristic of the SPGE nonlinear optical phenomenon. The relationships found for the longitudinal and transverse photocurrents as a function of pump polarization are also in agreement with the mechanism of photocurrent generation due to the SPGE. It is shown that the pulse durations of the longitudinal photocurrent for linear and circular polarizations, as well as the transverse photocurrent for linear polarization, are virtually independent of the incidence angle. However, the waveforms of the transverse photocurrent pulses at circular polarization with a given direction of rotation of the electric field vector of the incident radiation at a fixed sign of the angle of incidence strongly depend on the angle of incidence. At small and large angles of incidence, unipolar pulses of opposite polarity are generated, and in the intermediate range of incidence angles (58.5 ≤ α ≤ 76.5◦), bipolar photocurrent pulses are excited, smoothly transforming into unipolar pulses of opposite polarity at the boundaries of this interval. The obtained features of the waveform of the transverse photocurrent pulses at circular polarization of the incident radiation are due to the following: (i) the transformation of the circular polarization of the incident radiation into an elliptical one (without changing the sign of the circular polarization) upon the refraction of light at the air/semitransparent film interface and the appearance of a linear component of the photocurrent in the film structure, depending on the angle of incidence; (ii) the interaction of multidirectional linear and circular components of the photocurrent, which have different relaxation times and strongly depend on the angle of incidence.

The relationships found for the influence of the incidence angle on the waveform of pulses of longitudinal and transverse photocurrents that arise in the CuSe/*t*-Se film structure under polarized pulsed pumping due to the SPGE can be found in various materials in which the PSPC is excited by the PDE or by the CPGE.

The results obtained in this work can be used in optoelectronics, in particular, to create a high-speed detector capable of distinguishing left-handed and right-handed polarized light, as well as in the development of a technique that allows one to quickly determine the fast and slow axes of quarter-wave plates.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/app12146869/s1, Figure S1: X-ray diffractogram of the CuSe/*t*-Se nanocomposite film and diffraction patterns of CuSe (PDF 00-034-0171) and *t*-Se (PDF 00-042-1425) powders. A copper-based X-ray tube generating radiation at a wavelength of 0.1541 nm was used; Figure S2: Raman spectrum of the CuSe/*t*-Se nanocomposite film (black line). The red line represents the fitting of the combined Gaussian profiles at (blue line) 232 and (purple line) 237 cm−<sup>1</sup> corresponding to the *t*-Se and (green line) 262 cm−<sup>1</sup> to the CuSe nanocrystallites Raman resonances. A He-Ne laser radiation at a wavelength of 632.8 nm as excitation pumping was used; Figure S3: Scanning electron microscope image of the CuSe/*t*-Se semitransparent thin film surface; Figure S4: The optical transmittance spectrum of the CuSe/*t*-Se nanocomposite.

**Author Contributions:** Conceptualization, A.E.F. and G.M.M.; methodology, A.E.F., V.Y.K. and G.M.M.; validation, A.E.F., T.N.M., K.G.M. and G.M.M.; investigation, A.E.F., T.N.M., V.Y.K., K.G.M. and G.M.M.; resources, V.Y.K. and G.M.M.; writing—original draft preparation, A.E.F. and G.M.M.; writing—review and editing, all authors; supervision, G.M.M.; project administration, G.M.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Ministry of Education and Science of the Russian Federation (state registration number 1021032422167-7-1.3.2) and the Academy of Finland (Grant Nos. 323053, 340115).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors are grateful to Y.P. Svirko for discussing the results obtained. This study was performed using equipment of the Shared Use Center "Center of Physical and Physicochemical Methods of Analysis and Study of the Properties and Surface Characteristics of Nanostructures, Materials, and Products" UdmFRC UB RAS.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**

#### **Symbol**

	- σ Radiation incidence plane on the film
	- **n** Normal to the film surface

#### **Abbreviations**



#### **References**


**Mikhail V. Durnev and Sergey A. Tarasenko \***

Ioffe Institute, Politechnicheskaya 26, 194021 St. Petersburg, Russia

**\*** Correspondence: tarasenko@coherent.ioffe.ru

**Abstract:** Edges in two-dimensional structures are the source of nonlinear transport and optical phenomena which are particularly important in small-size flakes. We present a microscopic theory of the edge photogalvanic effect, i.e., the formation of DC electric current flowing along the sample edges in response to AC electric field of the incident terahertz radiation, for two-dimensional Dirac materials including the systems with massive and massless charge carriers. The edge current direction is controlled by the AC field polarization. The spectral dependence of the current is determined by the carrier dispersion and the mechanism of carrier scattering, as shown for single-layer and bilayer graphene as examples.

**Keywords:** edge currents; high-frequency nonlinear transport; photogalvanic effect; two-dimensional Dirac structures; massive and massless fermions

#### **1. Introduction**

The discovery of graphene and other two-dimensional (2D) crystals opened a new page in the physics of low-dimensional systems [1,2] and triggered the research aimed at the development of efficient sources and detectors of terahertz radiation based on 2D Dirac materials [3,4]. In small-size samples, e.g., flakes obtained by mechanical exfoliation, the important and sometimes decisive role in the formation of photoelectric response is played by edges and nearby regions [5,6]. At the edges, the translational and space inversion symmetries are naturally broken, which gives rise to edge-related mechanisms of the photogalvanic effect [5–11] and the second harmonic generation [12–14].

The photocurrents flowing along the sample edges (the edge photogalvanic effect) were observed and studied in single-layer and bilayer graphene samples excited by terahertz radiation [5,6,10], also in an external magnetic field in the conditions of cyclotron resonance [9] and in the regime of the quantum Hall effect [15]. It is found that the edge photocurrent is induced by both linearly and circularly polarized radiation. Moreover, the photocurrent direction is controlled by the polarization of the incident radiation: the electric field direction with respect to the edge for the linearly polarized field and the photon helicity for the circularly polarized field. The edge photogalvanic effect in 2D structures can be considered as a low-dimensional analog of the surface photogalvanic effect studied in bulk semiconductor crystals and metal films and recently in nanocomposite films [16–22].

The microscopic theory of the edge photogalvanic effect in the spectral range of intraband transport has been developed so far for 2D systems with parabolic energy spectrum of charge carriers [6,8]. Here, we generalize the theory to the class of 2D Dirac materials. We present a comprehensive theoretical study of the edge currents for 2D systems with arbitrary dispersion *ε*(*p*) and arbitrary type of electron scattering in the 2D bulk. We show that the edge current is determined by the dispersion of carriers and the relaxation times of the first and second angular harmonics of the distribution function and compare the results for single-layer and bilayer graphene, which are examples of 2D systems with parabolic and linear energy spectra.

**Citation:** Durnev, M.V.; Tarasenko, S.A. Edge Currents Induced by AC Electric Field in Two-Dimensional Dirac Structures. *Appl. Sci.* **2023**, *13*, 4080. https://doi.org/10.3390/ app13074080

Academic Editor: Gennady M. Mikheev

Received: 25 February 2023 Revised: 19 March 2023 Accepted: 20 March 2023 Published: 23 March 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **2. Microscopic Theory**

Consider an electromagnetic wave incident on the structure hosting a 2D electron gas occupying the (*xy*) half-plane at *x* ≥ 0, see Figure 1. The AC electric field of the incident wave has the form *E*(*t*) = *E* exp(−*iωt*) + c.c., where *E* and *ω* are the field amplitude and frequency, respectively, and the abbreviation "c.c." denotes the complex conjugation. The AC electric field causes the back-and-forth in-plane motion of electrons. At the edge of the structure (here, at *x* = 0), the AC motion of electrons gets distorted due to electron reflection from the edge and dynamic charge accumulation, which leads to an asymmetry of the high-frequency electron transport. This asymmetry results in the rectification of the AC current and, hence, the emergence of a DC current *Jy* flowing along the edge. As Figure 1 illustrates, the DC edge current can be excited by both linearly polarized and circularly polarized electromagnetic waves and the current direction is controlled by the wave polarization.

**Figure 1.** Illustration of the edge current formation. The back-and-forth motion of 2D carriers occupying the *x* ≥ 0 half-plane by linearly polarized (**a**,**b**) or circularly polarized (**c**,**d**) AC electric field results, due to electron scattering from the edge and dynamical charge accumulation, in the DC current *Jy* flowing along the edge. The edge current direction is controlled by the field polarization.

Now we present the microscopic theory of high-frequency non-linear electron transport and calculate the edge currently. We consider the classical range of the electromagnetic wave frequencies, i.e., *h*¯ *ω EF*, where *EF* is the Fermi energy of the 2D electron gas, and describe the electron kinetics by the Boltzmann equation

$$\frac{\partial f}{\partial t} + v\_x \frac{\partial f}{\partial x} + \epsilon \mathcal{E}\left(\mathbf{x}, t\right) \cdot \frac{\partial f}{\partial p} = \mathbf{I}\{f\} \,. \tag{1}$$

Here, *f* = *f*(*p*, *x*, *t*) is the electron distribution function, *p* and *ε*(*p*) are the electron momentum and energy, respectively, *v* = *dε*/*dp* = *vp*/*p* is the velocity, *v* = *dε*/*dp*, *e* is the electron charge, *E*(*x*, *t*) = *E*(*x*) exp(−*iωt*) + c.c. is the local electric field acting upon the electrons, and I{ *f* } is the collision integral. At this stage, we assume that the electron spectrum is isotropic in the 2D plane but do not specify the exact form of the dispersion *ε*(*p*). The field *E*(*x*) near the edge differs from the incident field *E* by the correction *δE*(*x*) *x* due to the screening produced by dynamical charge accumulation [6,23,24]. Therefore, E*<sup>y</sup>* = *Ey* whereas E*x*(*x*) ∝ *Ex*.

At equilibrium, the electron distribution is isotropic and homogenous at *x* ≥ 0 and is described by the Fermi–Dirac function *f*0(*ε*). In the presence of an AC electric field, the distribution function acquires corrections. We expand the resulting distribution function *f*(*p*, *x*, *t*) in the series in the electric field amplitude as follows

$$f(\mathbf{p}, \mathbf{x}, t) = f\_0 + \left[ f\_1(\mathbf{p}, \mathbf{x}) \exp(-i\omega t) + \text{c.c.} \right] + f\_2(\mathbf{p}, \mathbf{x}) + \dots \,\tag{2}$$

where *f*<sup>1</sup> is the first-order correction, which determines the linear (Drude) conductivity, and *f*<sup>2</sup> is the time-independent second-order correction. The second-order correction oscillating at 2*ω* is not considered here since it does not contribute to DC electric current.

The density of DC electric current *jy*(*x*) is determined by the asymmetric part of the correction *f*<sup>2</sup> and is given by

$$\dot{y}\_y(\mathbf{x}) = \text{eg}\sum\_{\mathbf{p}} v\_{\mathbf{y}} f\_2(\mathbf{p}, \mathbf{x})\,. \tag{3}$$

where *g* is the factor that takes into account possible spin and valley degeneracy (e.g., *g* = 2 for GaAs quantum wells and *g* = 4 for single-layer and bilayer graphene) and <sup>∑</sup>*<sup>p</sup>* = (2*πh*¯)−<sup>2</sup> *<sup>d</sup>*2*p*.

Equation (1) with the perturbation term *eE*(*x*, *t*) · *∂ f* /*∂p* yields the following differential equations for *f*<sup>1</sup> and *f*<sup>2</sup>

$$-i\omega f\_1 + v\_x \frac{\partial f\_1}{\partial x} + \epsilon \mathcal{E}(x) \cdot \frac{\partial f\_0}{\partial p} = \mathcal{I}\{f\_1\} \, , \tag{4}$$

$$\mathrm{c}\,v\_{\mathrm{x}}\frac{\partial f\_{2}}{\partial \boldsymbol{x}} + \left[\mathrm{e}\mathcal{E}(\boldsymbol{x}) \cdot \frac{\partial f\_{1}^{\*}}{\partial \boldsymbol{p}} + \mathrm{c.c.}\right] = \mathrm{I}\{f\_{2}\}\,. \tag{5}$$

We solve Equations (4) and (5) in the approximation of elastic electron scattering in the bulk of the 2D system and for specular reflection of electrons from the edge. The latter implies that *f*(*px*, *py*, 0, *t*) = *f*(−*px*, *py*, 0, *t*) which also ensures the lack of electric current across the edge. Multiplying Equation (5) by the velocity *vy* and averaging the resulting equation over the directions of *p* one obtains

$$
\left\langle v\_{\mathbf{x}} v\_{\mathbf{y}} \frac{\partial f\_2}{\partial \mathbf{x}} \right\rangle + \left\langle v\_{\mathbf{y}} \left( \epsilon \mathcal{E} \cdot \frac{\partial f\_1^\*}{\partial \mathbf{p}} + \text{c.c.} \right) \right\rangle = -\frac{\left\langle v\_{\mathbf{y}} f\_2 \right\rangle}{\tau\_1} \,, \tag{6}
$$

where the angular brackets ... stand for the averaging and *τ*<sup>1</sup> is the momentum relaxation time (relaxation time of the first angular harmonic) defined as 1/*τ*<sup>1</sup> = −*v I*{ *f* }/*v f*. Such a definition of *τ*<sup>1</sup> enables the consideration of its dependence on the electron energy *ε*(*p*). Equations (3) and (6) yield the equation for the current density

$$j\_y(\mathbf{x}) = -eg\sum\_p \tau\_1 v\_x v\_y \frac{\partial f\_2}{\partial \mathbf{x}} - \epsilon g \sum\_p \tau\_1 v\_y \left( e\mathcal{E} \cdot \frac{\partial f\_1^\*}{\partial \mathbf{p}} + \text{c.c.} \right). \tag{7}$$

After the integration of the second term by parts, Equation (7) assumes the form

$$\begin{split} \dot{\boldsymbol{\varepsilon}}\_{\mathcal{Y}}(\boldsymbol{x}) &= \ -\mathrm{e}\boldsymbol{\mathcal{g}}\sum\_{\mathcal{P}}\boldsymbol{\tau}\_{1}\boldsymbol{v}\_{\mathcal{X}}\boldsymbol{v}\_{\mathcal{Y}}\frac{\partial f\_{2}}{\partial \boldsymbol{x}} + \mathrm{e}^{2}\boldsymbol{\mathcal{g}}\sum\_{\mathcal{P}}\left(\frac{\boldsymbol{\tau}\_{1}}{m}\right)' m\boldsymbol{v}\_{\mathcal{X}}\boldsymbol{v}\_{\mathcal{Y}}(\mathcal{E}\_{\mathcal{X}}f\_{1}^{\*} + \text{c.c.}) \\ &+ \quad \boldsymbol{e}^{2}\boldsymbol{\mathcal{g}}\sum\_{\mathcal{P}}\left[\frac{\boldsymbol{\tau}\_{1}}{m} + \left(\frac{\boldsymbol{\tau}\_{1}}{m}\right)' \frac{m\boldsymbol{v}^{2}}{2} - m\left(\frac{\boldsymbol{\tau}\_{1}}{m}\right)' \frac{\boldsymbol{v}\_{\mathcal{X}}^{2} - \boldsymbol{v}\_{\mathcal{Y}}^{2}}{2}\right] (\boldsymbol{E}\_{\mathcal{Y}}f\_{1}^{\*} + \text{c.c.}) \end{split} \tag{8}$$

where (...) = *d*(...)/*dε* and *m* = *p*/*v* = *p*/(*dε*/*dp*) is the (energy-dependent) effective mass. In the case of parabolic dispersion *ε*(*p*) = *p*2/(2*m*∗) the effective mass *m* = *m*<sup>∗</sup> is energy-independent whereas for linear dispersion *ε*(*p*) = *v*<sup>0</sup> *p* the effective mass *m* = *ε*/*v*<sup>2</sup> 0 linearly depends on energy. Note that the mass *m* also determines the quasi-classical cyclotron motion. Equation (8) is valid for arbitrary dispersion and arbitrary boundary conditions at the edge. Obviously, the DC current vanishes for the AC field polarized along or perpendicularly to the edge and emerges only if the incident field *E* has both *x* and *y* components. Therefore, the corrections *f*<sup>1</sup> in the second and third terms on the right-hand side of Equation (8) should be calculated for the *y* and *x* components of the field, respectively.

For specular reflection of electrons from the edge, the second sum in Equation (8) vanishes since *f*<sup>1</sup> in response to *Ey* is an even function of *vx*. The third sum in Equation (8) can be rewritten via *∂ f* ∗ <sup>1</sup> /*∂x* using the equalities

$$i\omega \langle f\_1 \rangle\_{\!\!\! \! \! \! /} = \left\langle v\_x \frac{\partial f\_1}{\partial x} \right\rangle\_{\!\!\! \! \! /},\tag{9}$$

$$
\left< i\omega - \frac{1}{r\_2} \right> \left< \frac{v\_x^2 - v\_y^2}{2} f\_1 \right> \quad = \left< v\_x \frac{v\_x^2 - v\_y^2}{2} \frac{\partial f\_1}{\partial x} \right> , \tag{10}
$$

which follow from Equation (4). Here, *τ*<sup>2</sup> is the relaxation time of the second angular harmonic of the distribution function defined as 1/*τ*<sup>2</sup> = −(*v*<sup>2</sup> *<sup>x</sup>* − *<sup>v</sup>*<sup>2</sup> *<sup>y</sup>*) *<sup>I</sup>*{ *<sup>f</sup>* }/(*v*<sup>2</sup> *<sup>x</sup>* − *<sup>v</sup>*<sup>2</sup> *<sup>y</sup>*)*f*. Therefore, Equation (8) assumes the form

$$\begin{split} j\_{\mathcal{Y}}(\mathbf{x}) &= -\varepsilon g \sum\_{p} \tau\_{1} v\_{\mathcal{X}} v\_{\mathcal{Y}} \frac{\partial f\_{2}}{\partial \mathbf{x}} + \frac{\varepsilon^{2} g}{\omega \nu} \sum\_{p} v\_{\mathcal{X}} \left[ \frac{\tau\_{1}}{m} + \left( \frac{\tau\_{1}}{m} \right)' \frac{mv^{2}}{2} \right] \left( i \mathbb{E}\_{\mathcal{Y}} \frac{\partial f\_{1}^{\*}}{\partial \mathbf{x}} + \text{c.c.} \right) \\ &+ \varepsilon^{2} g \sum\_{p} mv\_{\mathcal{X}} \frac{\nu\_{\mathbf{x}}^{2} - v\_{\mathcal{Y}}^{2}}{2} \left( \frac{\tau\_{1}}{m} \right)' \left( \frac{\tau\_{2} E\_{\mathcal{Y}}}{1 + i \omega \tau\_{2}} \frac{\partial f\_{1}^{\*}}{\partial \mathbf{x}} + \text{c.c.} \right) . \end{split} \tag{11}$$

The current density *jy*(*x*) is determined by spatial derivatives of the distribution function and, as expected, vanishes in the 2D bulk where the electron distribution is homogenous.

The total electric current flowing along the edge is given by

$$J\_{\mathcal{Y}} = \bigcap\_{0}^{\infty} j\_{\mathcal{Y}}(\mathfrak{x})d\mathfrak{x}\,. \tag{12}$$

Integrating Equation (11) over *x* we obtain

$$\begin{split} f\_{\mathcal{Y}} &= -\varepsilon g \sum\_{p} \tau\_{1} v\_{x} v\_{y} [f\_{2}(p,\infty) - f\_{2}(p,0)] \\ &+ \left\{ i \frac{e^{2}g}{\omega} \sum\_{p} \left[ \frac{\tau\_{1}}{m} + \left(\frac{\tau\_{1}}{m}\right)' \frac{mv^{2}}{2} \right] v\_{x} [f\_{1}^{\*}(p,\infty) - f\_{1}^{\*}(p,0)] E\_{\mathcal{Y}} + \text{c.c.} \right\} \\ &+ \left\{ e^{2}g \sum\_{p} \left(\frac{\tau\_{1}}{m}\right)' \frac{mv\_{2}}{1 + i\omega \tau\_{2}} \frac{v\_{x}^{2} - v\_{y}^{2}}{2} v\_{x} [f\_{1}^{\*}(p,\infty) - f\_{1}^{\*}(p,0)] E\_{\mathcal{Y}} + \text{c.c.} \right\}, \end{split}$$

where *fn*(*p*, 0) and *fn*(*p*, ∞) are the corrections to the distribution function at the edge and far from the edge, respectively. For the particular case of specular reflection of electrons from the edge, the functions *f*1(*p*, 0) and *f*2(*p*, 0) are even in *px* and, therefore, the sums with *f*1(*p*, 0) and *f*2(*p*, 0) vanish. As a result, the edge current *Jy* is determined by the functions *f*1(*p*, ∞) and *f*2(*p*, ∞) in the bulk where the actual field *E* coincides with the incident field *E*. The sums with *f*1(*p*, ∞) and *f*2(*p*, ∞) can be readily calculated from Equations (4) and (5) neglecting spatial inhomogeneous terms and the electric field screening. Below, we do such calculations for the degenerate electron gas with the Fermi energy *EF*.

The first-order correction in the 2D bulk has the form

$$f\_1(p, \infty) = -\frac{e\tau\_1 f\_0'}{1 - i\omega\tau\_1}(\boldsymbol{v} \cdot \boldsymbol{E})\,. \tag{14}$$

Therefore, the second sum in Equation (13) is calculated as follows

$$\begin{aligned} \left[ \mathrm{ieg} \sum\_{p} \left[ \frac{\tau\_{1}}{m} + \left( \frac{\tau\_{1}}{m} \right)' \frac{m \upsilon^{2}}{2} \right] v\_{\mathcal{X}} f\_{1}^{\*} (p, \infty) E\_{\mathcal{Y}} + \mathrm{c.c.} &= \left[ \frac{\tau\_{1}}{m} + \left( \frac{\tau\_{1}}{m} \right)' \frac{m \upsilon^{2}}{2} \right]\_{E\_{\mathcal{Y}}} \left[ \mathrm{i} \upsilon^{\*} (\omega) E\_{\mathcal{X}}^{\*} E\_{\mathcal{Y}} + \mathrm{c.c.} \right] \\ &= \left[ \frac{\tau\_{1}}{m} + \left( \frac{\tau\_{1}}{m} \right)' \frac{m \upsilon^{2}}{2} \right]\_{E\_{\mathcal{Y}}} \mathrm{Re} \, \sigma(\omega) \left( \omega \tau\_{1} S\_{2} - S\_{3} \right), \end{aligned} \tag{15}$$

where *σ*(*ω*) is the conductivity,

$$
\sigma(\omega) = -\frac{e^2 g}{2} \sum\_{p} \frac{\tau\_1 v^2 f\_0'}{1 - i\omega \tau\_1} = \frac{ne^2}{m} \frac{\tau\_1}{1 - i\omega \tau\_1} \,\,\,\,\tag{16}
$$

all the values are taken at the Fermi level, *n* = *g* ∑*<sup>p</sup> f*<sup>0</sup> is the carrier density, and *S*<sup>2</sup> = *ExE*∗ *<sup>y</sup>* + *E*<sup>∗</sup> *xEy* and *S*<sup>3</sup> = *i*(*ExE*<sup>∗</sup> *<sup>y</sup>* − *E*<sup>∗</sup> *xEy*) are the Stokes parameters of the incident radiation. The third sum in Equation (13) is calculated as follows

$$\log \sum\_{p} \left(\frac{\tau\_1}{m}\right)' \frac{m\tau\_2}{1 + i\omega\tau\_2} \frac{\upsilon\_x^2 - \upsilon\_y^2}{2} \upsilon\_x f\_1^\*(p, \infty) E\_y + \text{c.c.} = \frac{1}{4} \left(\frac{\tau\_1}{m}\right)' \frac{m\upsilon^2}{1 + i\omega\tau\_2} E\_x^\* E\_y + \text{c.c.}$$

$$= \frac{1}{4} \left(\frac{\tau\_1}{m}\right)' \frac{m\upsilon^2}{1 + (\omega\tau\_2)^2} \frac{\tau\_2 \text{Re}\,\sigma(\omega)}{1 + (\omega\tau\_2)^2} \left[ (1 - \omega^2 \tau\_1 \tau\_2) S\_2 + \omega(\tau\_1 + \tau\_2) S\_3 \right]. \tag{17}$$

Lastly, the sum with *f*2(*p*, ∞) in Equation (13) can be expressed with the help of Equation (5) via the sum with *f*1(*p*, ∞) as follows

$$\sum\_{p} \mathbf{r}\_{1} \upsilon\_{x} \upsilon\_{y} f\_{2}(p, \infty) = -e \sum\_{p} \mathbf{r}\_{1} \mathbf{r}\_{2} \upsilon\_{x} \upsilon\_{y} \left( E \cdot \frac{\partial f\_{1}^{\*} (\infty)}{\partial p} + \text{c.c.} \right). \tag{18}$$

Integration of the right-hand side of Equation (18) by parts gives

$$\begin{split} \sum\_{p} \tau\_{1} v\_{x} v\_{y} f\_{2}(p, \infty) &= \left\{ c \sum\_{p} \left[ \frac{\tau\_{1} \tau\_{2}}{m} + \frac{m^{2} v^{2}}{2} \left( \frac{\tau\_{1} \tau\_{2}}{m^{2}} \right)' \right] (v\_{y} E\_{x}^{\*} + v\_{x} E\_{y}^{\*}) f\_{1}(p, \infty) + \text{c.c.} \right\} \\ &+ \left\{ c \sum\_{p} m^{2} \left( \frac{\tau\_{1} \tau\_{2}}{m^{2}} \right)' \frac{v\_{x}^{2} - v\_{y}^{2}}{2} (v\_{y} E\_{x}^{\*} - v\_{x} E\_{y}^{\*}) f\_{1}(p, \infty) \right\}. \end{split} \tag{19}$$

The above sums can be calculated similarly to the sums in Equations (15) and (17), which yields

$$\log \sum\_{p} \tau\_1 v\_x v\_y f\_2(p, \infty) = 2 \left[ \frac{\tau\_1 \tau\_2}{m} + \frac{m^2 v^2}{4} \left( \frac{\tau\_1 \tau\_2}{m^2} \right)' \right]\_{E\_F} \text{Re}\,\sigma(\omega) \,\text{S}\_2. \tag{20}$$

Finally, summing up all contributions to the edge current we obtain

$$J\_{\mathcal{Y}} = \frac{e\text{Re}\,\sigma(\omega)}{m} \left\{ \tau\_1 (\tau\_1 - 2\tau\_2) + \frac{m^2 \upsilon^2}{2} \left[ \frac{\tau\_1}{2} \left( \frac{\tau\_1}{m} \right)' - m \left( \frac{\tau\_1 \tau\_2}{m^2} \right)' \right] + \frac{m^2 \upsilon^2 (\tau\_1 + \tau\_2)}{4 [1 + (\omega \tau\_2)^2]} \left( \frac{\tau\_1}{m} \right)' \right\} S\_2$$
 
$$- \frac{e\text{Re}\,\sigma(\omega)}{m\omega} \left[ \tau\_1 + \frac{m^2 \upsilon^2 [2 + \omega^2 \tau\_2 (\tau\_2 - \tau\_1)]}{4 [1 + (\omega \tau\_2)^2]} \left( \frac{\tau\_1}{m} \right)' \right] S\_3. \tag{21}$$

Equation (21) represents the main result of this paper. It describes the DC edge current in 2D electron gas with an arbitrary electron dispersion *ε*(*p*) and arbitrary energy-dependent relaxation times.

For parabolic energy spectrum with *ε*(*p*) = *p*2/(2*m*∗), Equation (21) gives

$$\begin{split} J\_{\mathcal{Y}}^{(\text{par})} &= \frac{e \text{Re}\,\sigma(\omega)}{m^{\*}} \Big\{ \pi\_{1} (\tau\_{1} - 2\tau\_{2}) + \left[ \frac{\tau\_{1}\tau\_{1}^{\prime}}{2} - (\tau\_{1}\tau\_{2})^{\prime} \right] \varepsilon\_{F} + \frac{(\tau\_{1} + \tau\_{2})\tau\_{1}^{\prime}\varepsilon\_{F}}{2[1 + (\omega\tau\_{2})^{2}]} \Big\} \mathcal{S}\_{2} \\ &- \frac{e \text{Re}\,\sigma(\omega)}{m^{\*}\omega} \Big\{ \pi\_{1} + \frac{2 + \omega^{2}\tau\_{2}(\tau\_{2} - \tau\_{1})}{2(1 + \omega^{2}\tau\_{2}^{2})} \tau\_{1}^{\prime}\varepsilon\_{F} \Big\} \mathcal{S}\_{3} .\end{split} \tag{22}$$

This result was obtained previously in the approximation of a single energy-independent relaxation time (*τ*<sup>1</sup> = *τ*<sup>2</sup> = const) in Ref. [6] and in the form of Equation (22) in Ref. [8].

For linear energy spectrum *ε*(*p*) = *v*<sup>0</sup> *p*, Equation (21) gives

$$J\_{y}^{(\text{lin})} = \frac{ev\_{0}^{2}\text{Re}\,\sigma(\omega)}{2\varepsilon\_{F}} \left\{ \tau\_{1} \left(\frac{3\tau\_{1}}{2} - 2\tau\_{2}\right) + \left[\frac{\tau\_{1}\tau\_{1}^{\prime}}{2} - (\tau\_{1}\tau\_{2})^{\prime}\right] \varepsilon\_{F} + \frac{(\tau\_{1} + \tau\_{2})(\tau\_{1}^{\prime}\varepsilon\_{F} - \tau\_{1})}{2[1 + (\omega\tau\_{2})^{2}]} \right\} S\_{2}$$
 
$$= \frac{ev\_{0}^{2}\text{Re}\,\sigma(\omega)}{2\varepsilon\_{F}\omega} \left\{ 2\tau\_{1} + \frac{[2 + \omega^{2}\tau\_{2}(\tau\_{2} - \tau\_{1})](\tau\_{1}^{\prime}\varepsilon\_{F} - \tau\_{1})}{2[1 + (\omega\tau\_{2})^{2}]} \right\} S\_{3} \,\. \tag{23}$$

#### **3. Results and Discussion**

Now, we discuss the spectral and polarization dependence of the edge current in 2D systems with parabolic and linear dispersions, also for different scattering mechanisms. First, we note that the current *Jy* is proportional to the square of the electric field amplitude (*S*2, *S*<sup>3</sup> ∝ *E*2), i.e., to the intensity of the incident field. Therefore, it belongs to the class of photocurrents. Here, the photocurrent emerges due to the intraband (Drude-like) absorption of the electromagnetic field and is proportional to the real part of the high-frequency conductivity *σ*(*ω*).

The direction of the edge current (the polarity along the *y* axis) depends on the polarization state of the field via the Stokes parameters *S*<sup>2</sup> = (*ExE*<sup>∗</sup> *<sup>y</sup>* + *EyE*<sup>∗</sup> *<sup>x</sup>* ) and *S*<sup>3</sup> = *i*(*ExE*<sup>∗</sup> *<sup>y</sup>* − *EyE*<sup>∗</sup> *<sup>x</sup>* ). The contribution *Jy* ∝ *S*<sup>2</sup> is induced by linearly polarized radiation. This current is maximal if the radiation is polarized at the angle ±*π*/4 with respect to the edge and vanishes if the radiation is polarized along or perpendicular to the edge, Figure 1a,b. The contribution *Jy* ∝ *S*<sup>3</sup> describes the edge current induced by circularly polarized radiation and its direction is controlled by the radiation helicity, Figure 1c,d.

Equation (21) is quite general and can be applied to a wide class of 2D systems, including conventional III-V and II-VI quantum wells, bilayer graphene, and transition metal dichalcogenide monolayers with the parabolic spectrum, graphene and HgTe/CdHgTe quantum wells of the critical thickness with the linear spectrum, and narrow-gap 2D systems with the Dirac-like spectrum *ε*(*p*) = *v*<sup>0</sup> (*m*∗*v*0)<sup>2</sup> + *<sup>p</sup>*2. Below, we calculate *Jy* for single-layer and bilayer graphene.

Figure 2 shows the frequency dependence of the edge photocurrent *Jy* in a 2D system with the linear dispersion law *ε*(*p*) = *v*<sup>0</sup> *p*. The parameters used for calculations are given in the caption of Figure 2 and correspond to high-quality graphene [25] with the electron density *<sup>n</sup>* = <sup>5</sup> × <sup>10</sup><sup>11</sup> cm−2. The curves are calculated after Equation (23) for linearly polarized radiation with the electric field directed at the angle *π*/4 with respect to the edge (*S*2/*E*<sup>2</sup> = 1, *S*<sup>3</sup> = 0) and for circularly polarized radiation (*S*3/*E*<sup>2</sup> = 1, *S*<sup>2</sup> = 0). We consider two model types of scattering potential: (i) short-range scatterers, resulting in the energy dependence of the relaxation times *τ*<sup>1</sup> = 2*τ*<sup>2</sup> ∝ *ε*−1, and (ii) charged scatterers with the Coulomb potential, resulting in *τ*<sup>1</sup> = 3*τ*<sup>2</sup> ∝ *ε* [7]. For charged scatterers, the ratio *τ*1,2/*m* does not depend on energy and the edge current is given by *Jy* = *e*Re*σ*(*ω*)*v*<sup>2</sup> 0*τ*2 <sup>1</sup> (*EF*)/3*EF* for linearly polarized radiation and *Jy* = −*e*Re*σ*(*ω*)*v*<sup>2</sup> <sup>0</sup>*τ*1(*EF*)/*ωEF* for circularly polarized radiation. Hence, the frequency dependence of the edge currents excited by linearly and circularly polarized radiation are determined by Re*σ*(*ω*) ∝ 1/(1 + *ω*2*τ*<sup>2</sup> <sup>1</sup> ) and Re*σ*(*ω*)/*ω*, respectively. For short-range scatterers, all terms in Equation (21) contribute to the current, and the frequency dependence gets more complicated. In particular, the current induced by linearly polarized radiation is constant at low frequencies, changes its sign at an intermediate frequency, and decays as ∝ *ω*−<sup>2</sup> at high frequencies. The circular contribution

for short-range scatterers behaves as ∝ *ω* at low frequencies, in contrast to the diverging behavior <sup>∝</sup> *<sup>ω</sup>*−<sup>1</sup> for Coulomb scatterers. The current magnitude is *Jy* ∼ 10 nA for 1 W/cm2 of the radiation intensity.

**Figure 2.** Frequency dependence of the edge current in 2D systems with a linear dispersion of carriers. Solid curves correspond to short-range scatterers and linearly polarized field with *S*<sup>2</sup> = *E*<sup>2</sup> (green curve) and circularly polarized field with *S*<sup>3</sup> = *E*<sup>2</sup> (red curve). Dashed curves present the results for scattering by charged impurities. The curves are calculated after Equation (23) for the parameters corresponding to single-layer graphene: *<sup>v</sup>*<sup>0</sup> <sup>=</sup> 108 cm/s, *<sup>n</sup>* <sup>=</sup> <sup>5</sup> <sup>×</sup> 1011 cm−<sup>2</sup> (the Fermi energy *EF* ≈ 80 meV and the effective mass at the Fermi level *m* ≈ 0.013 *m*0), *τ*1(*EF*) = 1 ps, and *E* = 8 V/cm corresponding to the radiation intensity *I* = 1 W/cm2.

Figure 3 shows the frequency dependence of the edge photocurrent *Jy* in the 2D system with the parabolic dispersion law *ε*(*p*) = *p*2/2*m*∗. The parameters used for calculations are given in the caption of Figure 3 and correspond to high-quality bilayer graphene [6]. We use the same electron density *<sup>n</sup>* = <sup>5</sup> × 1011 cm−<sup>2</sup> as in the case of single-layer graphene. The curves are calculated after Equation (22) for linearly polarized radiation with the electric field directed at the angle *π*/4 with respect to the edge (*S*2/*E*<sup>2</sup> = 1, *S*<sup>3</sup> = 0) and for circularly polarized radiation (*S*3/*E*<sup>2</sup> = 1, *S*<sup>2</sup> = 0). For 2D systems with the parabolic spectrum and short-range scatterers, both the mass *m* and relaxation times *τ*<sup>1</sup> = *τ*<sup>2</sup> are independent of energy and Equation (21) yields *Jy* = −*e*Re*σ*(*ω*)*τ*<sup>2</sup> <sup>1</sup> (*EF*)/*m*<sup>∗</sup> for linearly polarized radiation and *Jy* = −*e*Re*σ*(*ω*)*τ*1(*EF*)/*ωm*<sup>∗</sup> for circularly polarized radiation. For Coulomb scatterers, *τ*<sup>1</sup> = 2*τ*<sup>2</sup> ∝ *ε*, and the frequency dependence of *Jy* is more complicated. Note that both the direction and magnitude of the current are determined to a great extent by the scattering mechanism. The calculated current magnitude for bilayer graphene is of the order of several nA for the radiation intensity 1 W/cm2 and is slightly smaller than that for monolayer graphene at the same electron density due to the larger effective mass in bilayer graphene.

Experimentally, edge photocurrents are detected as electric currents in short circuits or as voltage drops between contacts in open circuits [5,6,9,10]. If the sample is small and fully illuminated, the photocurrents are generated along all the edges and the resulting distribution of the photoinduced electric potential in the sample is determined by the edge photocurrents and the compensating drift currents in the sample. For linearly polarized radiation, the photocurrents generated along different edges are generally not equal because of different orientations of the edges with respect to the electric field polarization. For circularly polarized radiation, the edge photocurrents form a vortex whose winding direction depends on the radiation helicity. The photocurrents circulating around the sample produce, in turn, a magnetic field, which can be seen as a manifestation of the inverse Faraday effect.

**Figure 3.** Frequency dependence of the edge photocurrent in 2D systems with parabolic energy dispersion. Solid curves correspond to short-range scatterers and linearly polarized field with *S*<sup>2</sup> = *E*<sup>2</sup> (green curve) and circularly polarized field with *S*<sup>3</sup> = *E*<sup>2</sup> (red curve). Dashed curves present the results for scattering by charged impurities. The curves are calculated after Equation (22) for parameters corresponding to bilayer graphene: *<sup>m</sup>*<sup>∗</sup> <sup>=</sup> 0.03 *<sup>m</sup>*0, *<sup>n</sup>* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>11</sup> cm<sup>−</sup>2, *<sup>τ</sup>*1(*EF*) = 1 ps, and *E* = 8 V/cm corresponding to the radiation intensity *I* = 1 W/cm2.

#### **4. Conclusions**

To conclude, we have developed a kinetic theory of the edge photogalvanic effect for the intraband electron transport in two-dimensional materials. It is shown that the backand-forth motion of charge carriers by AC electric field of incident radiation is distorted at the edges of the sample resulting in direct electric currents flowing along the edges. The edge current direction is controlled by the radiation polarization while its spectral dependence is determined by the carrier dispersion and the mechanism of carrier scattering. We have obtained an analytical expression for the edge current valid for arbitrary dispersion law and scattering mechanism and analyzed the result for single-layer and bilayer graphene for electron scattering by short-range defects and Coulomb impurities. Considering the important role of edge regions in small-size samples such as flakes of two-dimensional crystals, one can expect that the edge photogalvanic effect will find applications in fast detectors of terahertz radiation and radiation polarization.

**Author Contributions:** M.V.D. and S.A.T. have contributed equality to the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the Russian Science Foundation (project 21-72-00047).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript: 2D Two-dimensional

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

MDPI St. Alban-Anlage 66 4052 Basel Switzerland www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

Disclaimer/Publisher's Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Academic Open Access Publishing

mdpi.com ISBN 978-3-0365-9371-5