Two-Step Asymmetric Perfectly Matched Layer Model for High-Order Spatial FDTD Solver of 2D Maxwell’s Equations

Abdelrahman Mahdy

Department of Physical Sciences, College of Science, Jazan University, Jazan, Kingdom of Saudi Arabia.

DOI: 10.4236/jamp.2025.132030                

Abstract

We implemented a two-step Asymmetric Perfectly Matched Layer (APML) model in High-Order Finite Difference Time Domain (FDTD) algorithm for solving two-dimensional Maxwell’s equations. Initially, we applied the APML method to the standard second-order FDTD algorithm to derive a two-step time-staggered APML (APML-2SS) and a two-step time-centered APML (APML-2SC) formulation for these equations, afterwards, we extended these formulations in high-order FDTD algorithm in order to derive a APML high-order FDTD (APML-HOFDTD) formulation for our Maxwell’s equations. To examine the performance and check out the accuracy of APML model, we conducted a numerical study using a 2D fluid where the three derived formulations were to analyze selected phenomena in terahertz radiation production by the filamentation of two femtosecond laser beams in air plasma. Numerical results illustrated that the two-step APML model is sufficiently accurate for solving our 2D Maxwell’s equations in high-order FDTD discretization and it demonstrated a great performance in studying the THz radiation production.

Keywords

Perfectly Matched Layer, The Finite-Difference-Time-Domain, Terahertz Radiation Production, Filamentation of Femtosecond Laser, Maxwell’s Equations Solver

1. Introduction

The Perfectly Matched Layer (PML) [1] is an Absorption Boundary Condition (ABC) algorithm that is presented to shorten the unbounded spatial domain in the computational electrodynamics [2]. The PML idea [3] is to block the reflection from the computational domain boundaries and at the same time, absorb the incident waves using an absorbing layer that is installed to absorb the incident waves without reflection. The PML is derived [4] by splitting the problem fields before the lower terms that simulate the absorption are added. The standard PML was initially introduced to solve Maxwell’s equations and Helmholtz equation [5] to simulate the electromagnetic wave propagation in unbounded media using the Finite-Difference-Time-Domain (FDTD) method [6]. In actual fact, the PML formulation is simple and its implementation is straightforward [3] in multi-dimensional Cartesian and cylindrical space, and it has been proven to be extremely efficient [7] in comparison with the classical Mur absorption algorithms. The standard PML has been developed [3] to solve 3D Maxwell’s equations to study the time-dependent wave propagation in free and arbitrary space, later, it has been upgraded for the real and complex spatial coordinates (CPML). The upgraded CPML at stretching conditions [7] offers an absorption for the seismic waves [8] and a solution for its second-order wave equation [9]. At that moment, a CPML frequency-shift (CFS-PML) [10] was implemented by shifting the frequency dependence off the real axis to improve the absorption of the EMW propagation at low grazing angles.

As an essential procedure in the standard PML algorithms, the filed equations must be splitted and then separately solved. To avoid the field splitting, Wang et al. [11] proposed a non-splitted PML method (NPML), in this NPML method, the conventional integration is replaced with the explicit time-matching scheme based on the trapezoidal integration. Following that, Ramadan [7] proposed another un-splitted field implementation for the standard PML by using auxiliary differential equation (ADE-PML) to update the auxiliary memory variables [12]. Thereupon, alternative unsplitted field PML formulations were reported for both loss and lossless media by the PML in Helmholtz scheme [13] [14], after which new unsplitted PML formulation was presented for absorbing seismic wave energy and a non-splitting split CPML for second order equation in elastic media [15]. It turns out that the unsplitted PML algorithms result in 25% memory saving [16] over the essential Berenger’s split field, as well, the Arbitrarily Wide-angle Wave Equation (AWWE) unsplitted CMPL requires less memory storage than the revised second-order AWWE with the split PML [12]. Over and above, the unsplitted CPML with a general CFS factor has proven to be more effective in absorbing the evanescent wave in media [9], beside it is highly absorbing for surface wave in elastic media.

In this article, we implemented a two-step Asymmetric Perfectly Matched Layer (APML) model [17] in the High-order FDTD algorithm in order to solve 2D Maxwell’s equations for terahertz radiation production by the filamentation of two femtosecond laser beams in air plasma. First of all, we imposed the APML boundary conditions in our basic equations, and then we applied the two-step APML method to derive two-step time-staggered APML (APML-2SS) and two-step time-centered APML (APML-2SC) formulation in the standard second-order FDTD algorithm. Subsequently, we extended these formulations for the high-order FDTD to derive the high-order finite-difference time-domain APML (APML-HOFDTD) formulation of our Maxwell’s equations. The article is organized as follows: In Section 2, we present our Maxwell’s equations with their assumptions and derivation conditions. In Section 3, we conduct numerical procedures to derive the APML-2SS and the APML-2SC formulation in the standard second-order FDTD algorithm, and the APML-HOFDTD formulation in the high-order FDTD discretization. In Section 4, we carried out a numerical study for fundamental phenomena in THz radiation production where the three formulations are applied in order to examine the performance and check out the accuracy of the two-step APML model.

2. Basic Maxwell’s Equations

The production of terahertz radiation by the filamentation of two femtosecond laser beams in air plasma is our main research interest [18]. Two theories govern the THz production in this medium, namely the Four-Wave Mixing rectification (FWM) [19] and the transient photocurrent model (PC) [18]. To carry out this research, a 2D fluid code that is referred to the PC theory is employed. The moment equation and the ionization equation coupled with a 2D Maxwell’s equations are the basic equations of this code. The moment equation is the standard equation of motion [20] for non-relativistic electrons dynamics in collisional cold plasma, while the ionization equation is the classical ADK tunneling ionization rate formula [21] that maintains the time-dependent ionization process after the threshold ionization. With respect to our central equations; the Maxwell’s equations where the two-step APML model is applied, it 2D equations that maintain the filamentation of a fundamental First Harmonic (FH) laser pulse (${\omega }_0$ ) and its corresponding Second Harmonic (SH) pulse ($2{\omega }_0$ ) in dielectric, dispersive, lossless, and isotropic media that is unchanged during this filamentation.

It is worthy noted that in our code Maxwell’s equations are derived [22] by taking the divergence of Gauss’s laws in combination with Faraday and Ampere-Maxwell’s laws to obtain these field equations:

$\nabla \times E\mathrm{<mo>=</mo>}-\frac{\partial B}{\partial t}$ , Faraday’s Law

$\nabla \times B\mathrm{<mo>=</mo>}\mu ϵ\frac{\partial E}{\partial t}$ , Ampere’s Law

where $E$ and $B$ are the electric and magnetic filed, and $ϵ\mathrm{,}\mu$ are the dielectric permittivity and the magnetic susceptibility, receptively. Concerning the air plasma, in our code it is an anisotropic plasma that is filed in a bounded 2D slab; along x and y-direction, where the electric and magnetic in this slab are lossy filed with standard boundary that fulfilled the Perfectly Matched Layer (PML) conditions.

As a main physical consideration in our study, the incident beam is a Transverse Electric (TE) polarized plane wave with field components ($E_x\mathrm{,}{E}_y\mathrm{,}{B}_z$ ). The magnetic field $B_z$ of this beam is along the direction of the propagation (z-direction) and it is normal to the slab interfaces (x and y). Inside this slab the electric fields $\left(E_x\mathrm{,}{E}_y\right)$ are calculated, while the magnetic field $B_z$ is derived. However, this magnetic filed is assumed to have anisotropic dependence on the intersections along these interfaces, thus this magnetic field can be re-presented as $B_z\to \left(B_{zx}\mathrm{,}{B}_{zy}\right)$ . Therefore, under this consideration and assumption, our Maxwell’s equations can be re-written as:

$\frac{\partial E_x}{\partial t}=\mu ϵ\frac{\partial B_z}{\partial y}\mathrm{,}$ (2.1)

$\frac{\partial E_y}{\partial t}=-\mu ϵ\frac{\partial B_z}{\partial x}\mathrm{.}$ (2.2)

$\frac{\partial B_{zx}}{\partial t}=-\frac{\partial E_y}{\partial x}$ (2.3)

$\frac{\partial B_{zy}}{\partial t}=\frac{\partial E_x}{\partial y}$ (2.4)

Equations (2.1)-(2.4) are our basic Maxwell’s equations where the APML model is going to be applied.

3. The Asymmetric Perfectly Matched Layer (APML) Model

The Asymmetric Perfect Matched Layer model (APML) is the standard PML boundary conditions [2] method with a particular constrains that are imposed on Maxwell’s equations that are then modified by extra added coefficients and more extracted terms. To apply the APML in our Maxwell’s equations, the field solutions in our problem are assumed to be explicit and linear [23]. Herein, these solutions; the electric field as an example, can be written as:

$\begin{array}{c}{E}_j^{i+1}={\alpha }_1{E}_j^i+{\beta }_{11}{B}_{j+1/2}^{i+1/2}-{\beta }_{12}{B}_{j-1/2}^{i+1/2}+\cdots +{\beta }_{1n}{B}_{j+n/2}^{i+1/2}\\ +{\alpha }_2{E}_j^{i-1}+{\beta }_{21}{B}_{j+1/2}^{i-1/2}-{\beta }_{22}{B}_{j-1/2}^{i+1/2}+\cdots +{\beta }_{2n}{B}_{j+n/2}^{i-1/2}\\ +\cdots \\ +{\alpha }_m{E}_j^{i-m+1}+{\beta }_{m1}{B}_{j+1/2}^{i-m+3/2}-{\beta }_{m2}{B}_{j-1/2}^{i-m+3/2}+\cdots +{\beta }_{mn}{B}_{j+n/2}^{i-m+3/2}\end{array}$

where ${\alpha }_m$ and ${\beta }_{mn}$ are constants. Then, with a particular constrains that are added to the electric filed in order to make its solution minimal and symmetric, the above electric filed series can be expressed in the from

$E_j^{i+1}=\alpha E_j^i+{\beta }_{py}{B}_{j+1/2}^{i+1/2}-{\beta }_{my}{B}_{j-1/2}^{i+1/2}\mathrm{.}$ (3.1)

$\alpha$ , ${\beta }_{py}$ , and ${\beta }_{my}$ are three constants, where ${\beta }_{py}={\beta }_{my}=c\Delta t/\Delta x$ , and

$\begin{array}{ll}\alpha =1+{\beta }_{py}-{\beta }_{my}\hfill & \text{for}\text{forward}\text{traveling}\text{waves}\mathrm{,}\hfill \\ \alpha =1-{\beta }_{py}+{\beta }_{my}\hfill & \text{for}\text{backward}\text{traveling}\text{waves}\mathrm{,}\hfill \end{array}$

and ${\beta }_{py}$ , and ${\beta }_{my}$ are determined by particular ways described in [23] [24].

In another alternative, the electric field in Equation (3.1) can be also re-expressed as:

$\frac{E_j^{i+1}-E_j^i}{\delta t}+{\sigma }_{E_y}\frac{E_j^{i+1}+E_j^i}{2}=C_{E_x}\frac{B_{j+1/2}^{i+1/2}-B_{j-1/2}^{i+1/2}}{\delta x}+{\sigma }_{B_y}\frac{B_{j+1/2}^{i+1/2}+B_{j-1/2}^{i+1/2}}{2}\mathrm{,}$ (3.2)

where,

$\begin{array}{ccc}{\sigma }_{E_y}=\frac{2}{\delta t}\frac{1-\alpha }{1+\alpha },& C_{E_x}=\frac{\delta y}{\delta }\frac{{\beta }_{py}+{\beta }_{my}}{1+\alpha },& {\sigma }_{B_y}=\frac{2}{\delta t}\frac{{\beta }_{py}-{\beta }_{my}}{1+\alpha }\end{array}$

At the infinitesimal limit [23], the alternative electric filed solution in Equation (3.2) becomes

$\frac{\partial E_x}{\partial t}+{\sigma }_{E_y}{E}_x=C_{E_x}\frac{\partial B_z}{\partial y}+{\overline{\sigma }}_{B_y}{B}_z\mathrm{.}$ (3.3)

Equation (3.3) is Equation (2.1) of our basic Maxwell’s equations where the APML boundary conditions are imposed. Similar procedures are carried out to the remaining equations; Equations (2.2)-(2.4), which become

$\frac{\partial E_y}{\partial t}+{\sigma }_{E_x}{E}_y=-C_{E_y}\frac{\partial B_z}{\partial x}+{\overline{\sigma }}_{B_x}{B}_z,$ (3.4)

$\frac{\partial B_{zx}}{\partial t}+{\sigma }_{B_x}{B}_{zx}=-C_{B_x}\frac{\partial Ey}{\partial x}+{\overline{\sigma }}_{E_x}{E}_y,$ (3.5)

$\frac{\partial B_{zy}}{\partial t}+{\sigma }_{B_y}{B}_{zy}=C_{B_y}\frac{\partial E_x}{\partial y}+{\overline{\sigma }}_{E_y}{E}_x\mathrm{,}$ (3.6)

where,

$\begin{array}{lll}{\sigma }_{E_x}=\frac{2}{\delta t}\frac{1-\alpha }{1+\alpha },\hfill & C_{E_x}=\frac{\delta x}{\delta t}\frac{{\beta }_{px}+{\beta }_{mx}}{1+\alpha },\hfill & {\sigma }_{B_x}=\frac{2}{\delta t}\frac{{\beta }_{px}-{\beta }_{mx}}{1+\alpha }\hfill \\ {\sigma }_{B_x}=\frac{2}{\delta t}\frac{1-\alpha }{1+\alpha },\hfill & C_{B_x}=\frac{\delta y}{\delta t}\frac{{\beta }_{px}+{\beta }_{mx}}{1+\alpha },\hfill & {\sigma }_{E_x}=\frac{2}{\delta t}\frac{{\beta }_{px}-{\beta }_{mx}}{1+\alpha }\hfill \\ {\sigma }_{B_y}=\frac{2}{\delta t}\frac{1-\alpha }{1+\alpha },\hfill & C_{B_y}=\frac{\delta y}{\delta }\frac{{\beta }_{py}+{\beta }_{my}}{1+\alpha },\hfill & {\sigma }_{E_y}=\frac{2}{\delta t}\frac{{\beta }_{py}-{\beta }_{my}}{1+\alpha }\hfill \end{array}$

After all, Equations (3.3)-(3.6) are our fundamental equations where the APML boundary conditions are applied.

3.1. Standard Second-Order FDTD Discretization

After we have applied the APML boundary conditions on our basic Maxwell’s equations, in this section, we implement the two-step APML model in our fundamental equations. For this purpose, the discretization of Equations (3.3)-(3.6) is essential, in this regards, the “Yee” staggered-grid Finite-Difference Time Domain (FDTD) discretization in space-time domain is commonly used. Among various FDTD discretization orders, we started with the second-order FDTD discretization. In this context, using the exponential time-stepping (${\text{e}}^{-{\sigma }_{E_x}},{\text{e}}^{-{\sigma }_{E_y}},{\text{e}}^{-{\sigma }_{B_x}},{\text{e}}^{-{\sigma }_{B_y}}$ ) for the electric and magnetic filed, with the assumption that the coefficients ${\overline{\sigma }}_{E_x}={\overline{\sigma }}_{E_y}={\overline{\sigma }}_{B_x}={\overline{\sigma }}_{B_y}=0$ , Equations (3.3)-(3.6) in the standard second-order FDTD discretization are written as:

${B_{zx}|}_{i+1/2,j+1/2}^{n+1/2}={\text{e}}^{-{\sigma }_{B_x}\Delta t}{B_{zx}|}_{i+1/2,j+1/2}^{n-1/2}-\frac{1-{\text{e}}^{-{\sigma }_{B_x}\Delta t}}{{\sigma }_{B_x}\Delta x}\frac{C_{B_x}}{c}\left({E_y|}_{i+1/2,k+1/2}^{n+1}-{E_y|}_{j,k+1/2}^{n+1}\right),$ (3.7)

${B_{zy}|}_{i+1/2,j+1/2}^{n+1/2}={\text{e}}^{-{\sigma }_{B_y}\Delta t}{B_{zy}|}_{i+1/2,j+1/2}^{n-1/2}-\frac{1-{\text{e}}^{-{\sigma }_{B_y}\Delta t}}{{\sigma }_{B_y}\Delta y}\frac{C_{B_y}}{c}\left({E_x|}_{i+1/2,j+1}^{n+1}-{E_x|}_{i+1/2,j}^{n+1}\right),$ (3.8)

${E_x|}_{i+1/2,j}^{n+1}={\text{e}}^{-{\sigma }_{E_y}\Delta t}{E_x|}_{i+1/2,j}^n+\frac{1-{\text{e}}^{-{\sigma }_{E_y}\Delta t}}{{\sigma }_{E_y}\Delta y}\frac{C_{E_y}}{c}\left({B_z|}_{i+1/2,j+1/2}^{n+1/2}-{B_z|}_{i+1/2,j-1/2}^{n+1/2}\right),$ (3.9)

${E_y|}_{i,j+1/2}^{n+1}={\text{e}}^{-{\sigma }_{E_x}\Delta t}{E_y|}_{i,j+1/2}^n-\frac{1-{\text{e}}^{-{\sigma }_{E_x}\Delta t}}{{\sigma }_{E_x}\Delta y}\frac{C_{E_x}}{c}\left({B_z|}_{i+1/2,j+1/2}^{n+1/2}-{B_z|}_{i-1/2,j+1/2}^{n+1/2}\right).$ (3.10)

where, $\Delta x\mathrm{,}\Delta y$ are the spatial grid size along $x$ and $y$ direction, respectively, $\Delta t$ is the time step, $i\mathrm{,}j$ are the spatial notation at the center position of the relevant cell, and $n$ is the time step.

Equations (3.7)-(3.10) are the two-step APML formulation in the second-order FDTD algorithm for Equations (3.3)-(3.6).

With the following special choices for the $C_{B_x}$ and $C_{B_y}$ coefficients:

$C_{B_x}=c{\text{e}}^{-{\sigma }_{B_x}\Delta t}\frac{{\sigma }_{B_x}\Delta t}{1-{\text{e}}^{-{\sigma }_{B_x}\Delta t}},$

$C_{B_y}=c{\text{e}}^{-{\sigma }_{B_y}\Delta t}\frac{{\sigma }_{B_y}\Delta t}{1-{\text{e}}^{-{\sigma }_{B_y}\Delta t}}.$

Equations (3.3)-(3.7) are re-written as:

${B_{zx}|}_{i+1/2,j+1/2}^{n+1/2}={\text{e}}^{-{\sigma }_{B_x}\Delta t}\left[{B_{zx}|}_{i+1/2,j+1/2}^{n-1/2}-\frac{\Delta t}{\Delta x}\left({E_y|}_{i+1,j+1/2}^n-{E_y|}_{j,j+1/2}^n\right)\right],$ (3.11)

${B_{zy}|}_{i+1/2,j+1/2}^{n+1/2}={\text{e}}^{-{\sigma }_{B_y}\Delta t}\left[{B_{zy}|}_{i+1/2,j+1/2}^{n-1/2}+\frac{\Delta t}{\Delta y}\left({E_x|}_{i+1/2,j+1}^n-{E_x|}_{i+1/2,j}^n\right)\right],$ (3.12)

${E_x|}_{i+1/2,j}^{n+1}={\text{e}}^{-{\sigma }_{E_y}\Delta t}\left[{E_x|}_{i+1/2,j}^n+\frac{\Delta t}{\Delta y}\left({B_z|}_{i+1/2,j+1/2}^{n+1/2}-{B_z|}_{i+1/2,j-1/2}^{n+1/2}\right)\right],$ (3.13)

${E_y|}_{i,j+1/2}^{n+1}={\text{e}}^{-{\sigma }_{E_x}\Delta t}\left[{E_y|}_{i,j+1/2}^n-\frac{\Delta t}{\Delta x}\left({B_z|}_{i+1/2,j+1/2}^{n+1/2}-{B_z|}_{i-1/2,j+1/2}^{n+1/2}\right)\right].$ (3.14)

Equations (3.11)-(3.14) are a novel two-step APML formulation in the second-order FDTD algorithm for Equations (3.3)-(3.6).

3.2. Two-Step Time-Staggered APML (APML-2SS)

When we push the electric and magnetic field into the vacuum and then reshaping the electric and magnetic field by the exponential coefficients ${\text{e}}^{-{\sigma }_{E_u}\Delta t}$ and ${\text{e}}^{-{\sigma }_{b_u}\Delta t}$ along $u$ -direction, respectively, the following modified formulation of Equations (3.11)-(3.14) are explored:

${{\tilde{B}}_{zx}|}_{i+1/2,j+1/2}^{n+1/2}={B_{zx}|}_{i+1/2,j+1/2}^{n-1/2}-\frac{\Delta t}{\Delta x}\left({E_y|}_{i+1,j+1/2}^n-{E_y|}_{i,j+1/2}^n\right),$ (3.15)

${{\tilde{B}}_{zy}|}_{i+1/2,j+1/2}^{n+1/2}={B_{zy}|}_{i+1/2,j+1/2}^{n-1/2}+\frac{\Delta t}{\Delta y}\left({E_x|}_{i+1/2,j+1}^n-{E_x|}_{i+1/2,j}^n\right),$ (3.16)

${B_{zx}|}_{i+1/2,j+1/2}^{n+1/2}={\text{e}}^{-{\sigma }_{B_x}\Delta t}{{\tilde{B}}_{zx}|}_{i+1/2,j+1/2}^{n+1/2},$ (3.17)

${B_{zy}|}_{i+1/2,j+1/2}^{n+1/2}={\text{e}}^{-{\sigma }_{B_y}\Delta t}{{\tilde{B}}_{zy}|}_{i+1/2,j+1/2}^{n+1/2},$ (3.18)

${{\tilde{E}}_x|}_{i+1/2,j}^{n+1}={E_x|}_{i+1/2,j}^n+\frac{\Delta t}{\Delta y}\left({{\tilde{B}}_z|}_{i+1/2,j+1/2}^{n+1/2}-{{\tilde{B}}_z|}_{i+1/2,j-1/2}^{n+1/2}\right),$ (3.19)

${{\tilde{E}}_y|}_{i,j+1/2}^{n+1}={E_y|}_{i,j+1/2}^n-\frac{\Delta t}{\Delta x}\left({{\tilde{B}}_z|}_{i+1/2,j+1/2}^{n+1/2}-{{\tilde{B}}_z|}_{i-1/2,j+1/2}^{n+1/2}\right),$ (3.20)

${E_x|}_{i+1/2,j}^{n+1}={\text{e}}^{-{\sigma }_{E_y}\Delta t}{{\tilde{E}}_x|}_{i+1/2,j}^{n+1},$ (3.21)

${E_y|}_{i,j+1/2}^{n+1}={\text{e}}^{-{\sigma }_{E_x}\Delta t}{{\tilde{E}}_y|}_{i,j+1/2}^{n+1}.$ (3.22)

Equations (3.15)-(3.22) are the two-step time-staggered APML (APML-2SS) formulation in the second-order FDTD algorithm for Equations (3.3)-(3.6).

3.3. Two-Step Time-Centered APML (APML-2SC)

The two-step time-centered APML (APML-2SC) is an alternative formulation of the (APML-2SS), where the field components are time-centered pushing into the vacuum. In this formulation, Maxwell’s equations are solved in one-time step and the field components are updated by multiplying each of this component by its corresponding damping coefficient ${\text{e}}^{-{\sigma }_n\Delta t}$ . Under these considerations, this alternative formulation becomes:

${{\tilde{B}}_{zx}|}_{i+1/2,j+1/2}^{n+1/2}={B_{zx}|}_{i+1/2,j+1/2}^n-0.5\frac{\Delta t}{\Delta x}\left({E_y|}_{i+1,j+1/2}^n-{E_y|}_{i,j+1/2}^n\right),$ (3.23)

${{\tilde{B}}_{zy}|}_{i+1/2,j+1/2}^{n+1/2}={B_{zy}|}_{i+1/2,j+1/2}^n-0.5\frac{\Delta t}{\Delta y}\left({E_x|}_{i+1/2,j+1}^n-{E_x|}_{i+1/2,j}^n\right),$ (3.24)

${{\tilde{E}}_x|}_{i+1/2,j}^{n+1}={E_x|}_{i+1/2,j}^n+\frac{\Delta t}{\Delta y}\left({{\tilde{B}}_z|}_{i+1/2,j+1/2}^{n+1/2}-{{\tilde{B}}_z|}_{i+1/2,j-1/2}^{n+1/2}\right),$ (3.25)

${{\tilde{E}}_y|}_{i,j+1/2}^{n+1}={E_y|}_{i,j+1/2}^n-\frac{\Delta t}{\Delta x}\left({{\tilde{B}}_z|}_{i+1/2,j+1/2}^{n+1/2}-{{\tilde{B}}_z|}_{i-1/2,j+1/2}^{n+1/2}\right),$ (3.26)

${{\tilde{B}}_{zx}|}_{i+1/2,j+1/2}^{n+1}={{\tilde{B}}_{zx}|}_{i+1/2,j+1/2}^{n+1/2}-0.5\frac{\Delta t}{\Delta x}\left({{\tilde{E}}_y|}_{i+1/2,j+1/2}^{n+1}-{{\tilde{E}}_y|}_{i,j+1/2}^n\right),$ (3.27)

${{\tilde{B}}_{zy}|}_{i+1/2,j+1/2}^{n+1/2}={{\tilde{B}}_{zy}|}_{i+1/2,j+1/2}^{n+1/2}+0.5\frac{\Delta t}{\Delta y}\left({{\tilde{E}}_x|}_{i+1/2,j+1}^n-{{\tilde{E}}_x|}_{i+1/2,j}^n\right),$ (3.28)

${B_{zx}|}_{i+1/2,j+1/2}^{n+1}={\text{e}}^{-{\sigma }_{B_x}\Delta t}{{\tilde{B}}_{zx}|}_{i+1/2,j+1/2}^{n+1},$ (3.29)

${B_{zy}|}_{i+1/2,j+1/2}^{n+1}={\text{e}}^{-{\sigma }_{B_y}\Delta t}{{\tilde{B}}_{zy}|}_{i+1/2,j+1/2}^{n+1},$ (3.30)

${E_x|}_{i+1/2,j}^{n+1}={\text{e}}^{-{\sigma }_{E_y}\Delta t}{{\tilde{E}}_x|}_{i+1/2,j}^{n+1},$ (3.31)

${E_y|}_{i,j+1/2}^{n+1}={\text{e}}^{-{\sigma }_{E_x}\Delta t}{{\tilde{E}}_y|}_{i,j+1/2}^{n+1}.$ (3.32)

Equations (3.23)-(3.32) are the two-step time-centered APML (APML-2SC) formulation of the second-order FDTD algorithm for Equations (3.3)-(3.6).

3.4. Two-Step High-Order Spatial FDTD (APML-HOFDTD)

The two-step APML High-order Spatial FDTD (APML-HOFDTD) is an extension formulation for the APML-2SS and APML-2SC where the spatial derivatives of our fundamental Maxwell’s equations are highly-order Finite-Difference Time Domain (FDTD) discretized in the “Yee” staggered spacial domain. Thus, the two-step APML formulation of the high-order FDTD algorithm (APML-HOFDTD) for Equations (3.3)-(3.6) are given by

${B_{zx}|}_{i+1/2,j+1/2}^{n+1/2}={B_{zx}|}_{i+1/2,j+1/2}^{n-1/2}-\frac{\Delta t}{\Delta x}{\displaystyle \sum _{l=1}^{p/2}}{c}_l^p\left({E_y|}_{i+l,j+1/2}^n-{E_y|}_{i-\left(l-1\right),j+1/2}^n\right),$ (3.33)

${B_{zy}|}_{i+1/2,j+1/2}^{n+1/2}={B_{zy}|}_{i+1/2,j+1/2}^{n-1/2}+\frac{\Delta t}{\Delta y}{\displaystyle \sum _{l=1}^{p/2}}{c}_l^p\left({E_x|}_{i+1/2,j+1+l}^n-{E_x|}_{i+1/2,j-\left(l-1\right)}^n\right),$ (3.34)

${E_x|}_{i+1/2,j}^{n+1}={E_x|}_{i+1/2,j}^n+\frac{\Delta t}{\Delta y}{\displaystyle \sum _{l=1}^{p/2}}{c}_l^p\left({B_z|}_{i+1/2+\left(l-1\right),j+1/2}^{n+1/2}-{B_z|}_{i+1/2-l,j-1/2}^{n+1/2}\right),$ (3.35)

${E_y|}_{i,j+1/2}^{n+1}={E_y|}_{i,j+1/2}^n-\frac{\Delta t}{\Delta x}{\displaystyle \sum _{l=1}^{p/2}}{c}_l^p\left({B_z|}_{i+1/2,j+1/2+\left(l-1\right)}^{n+1/2}-{B_z|}_{i-1/2,j+1/2-l}^{n+1/2}\right).$ (3.36)

where $c_{\ell }^p\mathrm{,}\ell =\mathrm{1,2,}\cdots ,p/2$ are the Fornberg coefficients [25] for any arbitrary discretization order $p$ of the splited Maxwell’s equations in vacuum. These coefficients are calculated by the Taylor expansion-based closed-form formulas [26], and given on the staggered grid by [27]

$c_{\ell }^p=\frac{{\left(-1\right)}^{\ell +1}{16}^{1-p/2}\left(p-1\right){\mathrm{<mo>!</mo>}}^2}{{\left(2\ell -1\right)}^2\left(p/2+\ell -1\right)\mathrm{<mo>!</mo>}\left(p/2-\ell \right)\mathrm{<mo>!</mo>}\left(p/2-\ell \right){\mathrm{<mo>!</mo>}}^2}\mathrm{,}\ell =\mathrm{1,2,}\cdots ,p/2\mathrm{.}$

4. Numerical Study

After we derived three different two-step APML formulations for our Maxwell’s equations, it is necessary to examine the performance and check out the accuracy and the reliability of these formulations. Therefore, in this section, using a 2D fluid code; where these formulations are applied, we carried out a numerical study for selected fundamental phenomena in THz radiation production. In this study, the initial input TE fs beam profile is given by [28]

$E\left(t\right)=E_L\left(\sqrt{1-\zeta }\sin \left({\omega }_{0}t\right){\text{e}}^{-\frac{t^2}{2{\tau }_0^2}}+\sqrt{\zeta }\sin \left(2{\omega }_{0}t+\varphi \right){\text{e}}^{-\frac{t^2}{{\tau }_0^2}}\right)\mathrm{.}$ (4.1)

where $E_m=\sqrt{2I_0/{ϵ}_{0}c}$ is the amplitude, $I_0$ is the initial input intensity, $\zeta$ is the energy fraction factor of SH pulse, $\varphi$ is the relative phase between the FH and SH pulse, ${\omega }_0\left(=2\pi c/{\lambda }_0\right)$ is the fundamental frequency, and ${\tau }_0$ is the initial pulse duration. In addition, in this study the 2D air plasma slab has equal dimension $x=y=50\text{-}100\mu \text{m}$ , and it is filled with plasma at initial electron density $n_e=3.7\times {10}^{19}{\text{cm}}^{-3}$ .

4.1. Optimizing the Incident Beam Parameters

As clearly listed in the initial beam profile given in Equation (4.1), in whole, the fundamental frequency ${\omega }_0$ , the initial input intensity $I_0$ , the initial pulse duration ${\tau }_0$ , the fractional energy factor $\zeta$ , and the relative phase $\varphi$ are the key parameters of this profile. For an efficient THz production, suitable initial values have to be given for these parameters. In this regards, ${\omega }_0=800\text{nm}$ value is necessary given for best coupling between the FH and the SH pulse [29], and $I_0\approx {10}^{14}~{10}^{17}\text{W}/{\text{cm}}^{\text{2}}$ period is normally selected, since this period is above the threshold intensity of the air gases ionization and below the relativistic intensity to avoid the relativistic dynamics, and ${\tau }_0=40\text{fs}$ is the pulse duration value where the yield and the efficiency of the produced THz radiation is enhancing [20] in the presence of the pulse chirping effects. With respect to $\zeta$ and $\varphi$ , in actual fact, the suitable values of these two parameters can’t be given, but it should be optimized in accordance with the conditions of the problem of our interest. In our previous studies [30] and other similar researches [31] [32], $\zeta =0.5$ value is acknowledged to avoid the energy exchange between the two fs beams and arrest the Modulation Instability (MI) effects [33], and $\varphi =\pi /2$ is identified for an efficient THz radiation within the selected input intensity period.

To examine the performance of the derived formulations, in this section, we re-optimized the suitable fractional energy factor $\zeta$ and the relative phase $\varphi$ value for an efficient THz production. In this respect, we simulated the excitation energy ${\mu }_{THz}$ of produced THz radiation as function of $\zeta$ and $\varphi$ , the possible energy fraction factor $\zeta =0\text{-}1$ and relative phase $\varphi =0\text{-}\pi$ values are considered in this simulation, and an initial input intensity value withing the selected period $I_0={10}^{15}\text{W}/{\text{cm}}^{\text{2}}$ is given, the result of this simulation is displayed in Figure 1.

Figure 1. The excitation energy ${\mu }_{THz}$ of the generated THz radiation using the APML-2SS (left), the APML-2SC (middle), and the APML-HOFDTD formulation (right).

As clearly noted in the ${\mu }_{THz}$ spectrum demonstrated in Figure 1, the maximum THz radiation is approximately obtained at $\zeta \approx \pi /2$ and $\varphi \approx \pi /2$ for the three formulations. It should be emphasized that, although the same acknowledged $\zeta$ and identified $\varphi$ value are also obtained for the three formulations in this figure, the APML-HOFDTD spectrum around these values is more refined than APML-2SS and APML-2SC formulation.

4.2. The Conversion Efficiency of the Produce THz Radiation

The optical-to-THz radiation conversion efficiency ${\eta }_{THz}$ is the ratio between the intensity of the produced THz radiation ${\left|E_{THz}\right|}^2$ and the input beam intensity ${\left|E_0\left(r\mathrm{,}t\right)\right|}^2$ [34],

${\eta }_{THz}={\displaystyle {\int }_0^{THz}}{\left|E_{THz}\left(\omega \right)\right|}^2\text{d}\omega /{\displaystyle {\int }_0^t}{\left|E_0\left(r\mathrm{,}t\right)\right|}^2\text{d}t\mathrm{.}$

For an efficient THz radiation production, a considerable conversion rate is crucially required. For this purpose, numerous researches have been carried out [29] [35] in order to calculate the conversion efficiency and analyze its behavior at different input beam parameters and air plasma structures. ${\eta }_{THz}={10}^{-3}\text{-}{10}^{-4}$ conversion rate is achieved in these researches, which is recognized as standard conversion efficiency value in THz radiation production topic.

In order the check the accuracy of the three derived formulations, herein we carried out a numerical comparison between the conversion efficiency as function of the selected input intensity period obtained by theses formulations and the standard ${\eta }_{THz}$ conversion [36] [37]. The optimum energy fraction factor $\eta =0.5$ and relative phase $\varphi =\pi /2$ , in addition to the suitable initial pulse duration ${\tau }_0=40\text{fs}$ are considered in this comparison, the result of this comparison is shown in Figure 2.

Figure 2. Comparison between the conversion efficiency of the produced THz radiation ${\eta }_{THz}$ as function of the input intensity obtained using the APML-2SS, APML-2SC, and the APML-HOFDTD formulations, and the standrad conversion efficiency.

As clearly seen in Figure 2, at low input intensity range, the three formulations have almost exact behavior to the standard conversion efficiency, moreover the maximum efficiency value obtained by APML-HOFDTD is approximately identical with the standard value, while the maximum values obtained by the APML-2SS and APML-2SC are little lower. At the high input intensity range, a discrepancy between the three formulations behavior and the standard one is emerged, although of that, the APML-HOFDTD has an insignificant discrepancy, while the APML-2SS and APML-2SC have a growing discrepancy with increasing the input intensity in this range.

4.3. The Spectral Broadening of the Produced THz Radiation

The Self-Phase Modulation (SPM) is the variation of the pulse phase as it is propagating [38]. Due to this variation, an addition time-dependent phase is added ($\Delta {\varphi }_{SPM}$ ). In nonlinear optical medium, the SPM is originated due to the Kerr-nonlinearity effect [38] where intensity is temporarily varying by the own pulse and due to the plasma induced SPM effect [20] where time-varying electron density $n\left(t\right)$ is evolved due the ionization. In both effects, and time dependent refractive index $n\left(r\mathrm{,}t\right)=n_0+\Delta n\left(t\right)$ is established. The presence of the SPM leads to spectral broadening of the pulse due to the generation of extra frequencies $\Delta \omega$ that are added to the central frequency (${\omega }_0$ ) of the pulse $\omega \left(t\right)={\omega }_0+\Delta \omega \left(t\right)$ . In this part, we checked the reliability of the three applied formulations in monitoring this broadening. In this context, we presented in Figure 3 a time-evolution for the produced THz radiation envelope where the three formulations are applied. Physical assumptions and considerations are assumed where a laser-envelope evolution which has a linear phase modulation with a symmetric frequency variation from the central frequency of the pulse is proceeding.

Figure 3. The time-evolution for the induced THz radiation envelope using the APML-2SS (left), the APML-2SC (middle), and the APML-HOFDTD formulation (right).

As illustrated in Figure 3, in a comparison of the time-evolution among the three envelopes, it clearly noted a small spectral broadening is ongoing with the APML-2SS and APML-2SC, while a clear spectral broadening is proceeding with the APML-HOFDTD formulation.

5. Conclusion

We successively implemented a two-step Asymmetric Perfect Matched Layer (APML) model with the High-order FDTD algorithm for solving 2D Maxwell’s equations. In the constrains of the APML boundary conditions, we derived three different formulations for these equations: namely, the two-step time-staggered APML (APML-2SS) and two-step time-centered APML (APML-2SC) formulation with the standard second-order FDTD algorithm, and the high-order FDTD APML (APML-HOFDTD) formulation with the high-order FDTD discretization. In the numerical study for selected fundamental phenomena in THz radiation production by the filamentation of two fs laser beams in air plasma where the three formulations are applied, the numerical results affirmed the validity of the two-step APML model for solving our Maxwell’s equations and the reliability in studying the THz radiation production. Although the three derived formulations were equivalently presented a proper performance in this study, the APML-HOFDTD formulation demonstrated a better performance and a more refining behavior in these results in comparison with the APML-2SS and APML-2SC formulations.

Acknowledgements

The author acknowledges the use of the High-Performance Computing Facility at Bibliotheca Alexandrina, and the associated support services in the completion of this work.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Berenger, J.-P. (1996) Three-Dimensional Perfectly Matched Layer for the Absorption of Electromagnetic Waves. Journal of Computational Physics, 127, 363-379. https://doi.org/10.1006/jcph.1996.0181
[2]	Brenger, J.-P. (2007) Perfectly Matched Layer (PML) for Computational Electromagnetics. Springer International Publishing. https://doi.org/10.1007/978-3-031-01696-7
[3]	Cummer, S.A. (2003) A Simple, Nearly Perfectly Matched Layer for General Electromagnetic Media. IEEE Microwave and Wireless Components Letters, 13, 128-130. https://doi.org/10.1109/LMWC.2003.810124
[4]	Holland, R. and Williams, J.W. (1983) Total-Field Versus Scattered-Field Finite-Difference Codes: A Comparative Assessment. IEEE Transactions on Nuclear Science, 30, 4583-4588. https://doi.org/10.1109/TNS.1983.4333175
[5]	Ji, D.Q., Wang, Q., Chen, Z.X., Harding, T.G., Dong, M.Z., Ge, C.Q., Li, Z.P. and Lai, F.P. (2020) A Property-Dependent Perfectly Matched Layer with a Single Additional Layer for Maxwells Equations in Finite Difference Frequency Domains. Computer Methods in Applied Mechanics and Engineering, 372, Article ID: 113355. https://doi.org/10.1016/j.cma.2020.113355
[6]	Taflove, A., Susan, C., Hagness, and Piket-May, M. (2005) 9-Computational Electromagnetics: The Finite-Difference Time-Domain Method. In: Chen, W.-K., Ed., The Electrical Engineering Handbook, Academic Press, Burlington, 629-670. https://doi.org/10.1016/B978-012170960-0/50046-3
[7]	Ramadan, O. (2003) Auxiliary Differential Equation Formulation: An Efficient Implementation of the Perfectly Matched Layer. IEEE Microwave and Wireless Components Letters, 13, 69-71. https://doi.org/10.1109/LMWC.2003.808706
[8]	Festa, G. and Vilotte, J.-P. (2005) The Newmark Scheme as Velocity Stress Time-Staggering: An Efficient PML Implementation for Spectral Element Simulations of Elastodynamics. Geophysical Journal International, 16, 789-812. https://doi.org/10.1111/j.1365-246X.2005.02601.x
[9]	Alan Roden, J. and Gedney, S.D. (2000) Convolution PML (CPML): An Efficient FDTD Implementation of the CFS-PML for Arbitrary Media. Microwave and Optical Technology Letters, 27, 334-339. https://doi.org/10.1002/1098-2760(20001205)27:5<334::AID-MOP14>3.0.CO;2-A
[10]	Chew, W.C. and Weedon, W.H. (1994) A 3d Perfectly Matched Medium from Modified Maxwell’S Equations with Stretched Coordinates. Microwave and Optical Technology Letters, 7, 599-604. https://doi.org/10.1002/mop.4650071304
[11]	Wang, T.S.l. and Tang, X.M. (2003) Finite? Difference Modeling of Elastic Wave Propagation: A Nonsplitting Perfectly Matched Layer Approach. Geophysics, 68, Article ID: 1620648. https://doi.org/10.1190/1.1620648
[12]	Chen, H.M., Zhou, H. and Li. Y.Q. (2014) Application of Unsplit Convolutional Perfectly Matched Layer for Scalar Arbitrarily Wide-Angle Wave Equation. Geophysics, 79, T313-T321. https://doi.org/10.1190/geo2014-0103.1
[13]	Fang, J.Y. and Wu, Z.h. (1995) Generalized Perfectly Matched Layer—An Extension of Berenger’S Perfectly Matched Layer Boundary Condition. IEEE Microwave and Guided Wave Letters, 5, 451-453. https://doi.org/10.1109/75.481858
[14]	Xu, F., Zhang, Y.l., Hong, W., Wu, K. and Cui, T.-J. (2003) Finite-Difference Frequency-Domain Algorithm for Modeling Guided-Wave Properties of Substrate Integrated Waveguide. IEEE Transactions on Microwave Theory and Techniques, 51, 2221-2227. https://doi.org/10.1109/TMTT.2003.818935
[15]	Connolly, D.P., Giannopoulos, A. and Forde, M.C. (2015) A Higher Order Perfectly Matched Layer Formulation for Finite-Difference Time-Domain Seismic Wave Modeling. Geophysics, 80, T1-T16. https://doi.org/10.1190/geo2014-0157.1
[16]	Petropoulos, P.G., Zhao, L. and Cangellaris, A.C. (1998) A Reflectionless Sponge Layer Absorbing Boundary Condition for the Solution of Maxwell’s Equations with High-Order Staggered Finite Difference Schemes. Journal of Computational Physics, 139, 184-208. https://doi.org/10.1006/jcph.1997.5855
[17]	Shapoval, O., Vay, J.-L. and Vincenti, H. (2019) Two-Step Perfectly Matched Layer for Arbitrary-Order Pseudo-Spectral Analytical Time-Domain Methods. Computer Physics Communications, 235, 102-110. https://doi.org/10.1016/j.cpc.2018.09.015
[18]	Kim, K.Y. Taylor, A.J., Glownia, J.H. and Rodriguez, G. (2008) Coherent Control of Terahertz Supercontinuum Generation in Ultrafast Lasergas Interactions. Nature Photonics, 2, 605-609. https://doi.org/10.1038/nphoton.2008.153
[19]	Cook, D.J. and Hochstrasser, R.M. (2000) Intense Terahertz Pulses by Four-Wave Rectification in Air. Optics Letters, 25, 1210-1212. https://doi.org/10.1364/OL.25.001210
[20]	Nguyen, A., González de Alaiza Martínez, P., Thiele, I., Skupin, S. and Bergé, L. (2018) THz Field Engineering in Two-Color Femtosecond Filaments Using Chirped and Delayed Laser Pulses. New Journal of Physics, 20, Article ID: 033026. https://doi.org/10.1088/1367-2630/aaa470
[21]	Ammosov, M.V., Delone, N.B. and Krainov, V.P. (1986) Tunnel Ionization of Complex Atoms and Atomic Ions in Electromagnetic Field. Soviet Physics-JETP, 64, 138-141. https://doi.org/10.1117/12.938695
[22]	Jackson, J.D. (1998) Classical Electrodynamics. 3rd Edition, Wiley.
[23]	Vay, J.-L. (2002) Asymmetric Perfectly Matched Layer for the Absorption of Waves. Journal of Computational Physics, 183, 367-399. https://doi.org/10.1006/jcph.2002.7175
[24]	Vay, J.-L. (2001) An Extended FDTD Scheme for the Wave Equation: Application to Multiscale Electromagnetic Simulation. Journal of Computational Physics, 167, 72-98. https://doi.org/10.1006/jcph.2000.6659
[25]	Fornberg, B. (1990) High-Order Finite Differences and the Pseudospectral Method on Staggered Grids. SIAM Journal on Numerical Analysis, 27, 904-918. https://doi.org/10.1137/0727052
[26]	Khan, I.R. and Ohba, R. (1999) Closed-Form Expressions for the Finite Difference Approximations of First and Higher Derivatives Based on Taylor Series. Journal of Computational and Applied Mathematics, 107, 179-193. https://doi.org/10.1016/S0377-0427(99)00088-6
[27]	Vincenti, H. and Vay, J.-L. (2016) Detailed Analysis of the Effects of Stencil Spatial Variations with Arbitrary High-Order Finite-Difference Maxwell Solver. Computer Physics Communications, 200, 147-167. https://doi.org/10.1016/j.cpc.2015.11.009
[28]	Mahdy, A.I. and Eltayeb, H.A. (2024) The Efficiency of Terahertz Radiation Generated by Two Chirped Femtosecond Laser Pulses at Different Pulse Durations. AIP Advances, 14, Article ID: 015345. https://doi.org/10.1063/5.0185372
[29]	Zhang, Z., Panov, N., Andreeva, V., Zhang, Z.l., Slepkov, A., Shipilo, D., Thomson, M.D., Wang, T.-J., Babushkin, I., Demircan, A., Morgner, U., Chen, Y.P., Kosareva, O. and Savel’ev, A. (2018) Optimum Chirp for Efficient Terahertz Generation from Two-Color Femtosecond Pulses in Air. Applied Physics Letters, 113, Article ID: 241103. https://doi.org/10.1063/1.5053893
[30]	Mahdy, A.I. and Eltayeb, H.A. (2024) The Conversion Efficiency of Terahertz Radiation Produced by the Filamentation of Two Femtosecond Laser Beams in Highly Charged Noble Gases. Physics of Plasmas, 31, Article ID: 122108. https://doi.org/10.1063/5.0233676
[31]	Babushkin, I., Kuehn, W., Köhler, C., Skupin, S., Berge, L., Reimann, K., Woerner, M., Herrmann, J. and Elsaesser, T. (2010) Ultrafast Spatiotemporal Dynamics of Terahertz Generation by Ionizing Two-Color Femtosecond Pulses in Gases. Physical Review Letters, 105, Article ID: 053903. https://doi.org/10.1103/PhysRevLett.105.053903
[32]	Kim, K.Y., Glownia, J.H., Taylor, A.J. and Rodriguez, G. (2007) Terahertz Emission from Ultrafast Ionizing Air in Symmetry-Broken Laser Fields. Optics Express, 15, 4577-4584. https://doi.org/10.1364/OE.15.004577
[33]	Ding, P.J., Guo, Z.Q., Wang, X.S., Cao, Y., Sun, M.Z., Zhao, P.X., Shi, Y.C., Sun, S.H., Liu, X.L. and Hu, B.T. (2013) Energy Exchange between Two Noncollinear Filament-Forming Laser Pulses in Air. Optics Express, 21, 27631-27640. https://doi.org/10.1364/OE.21.027631
[34]	Houard, A., Liu, Y., Mysyrowicz, A. and Leriche, B. (2007) Calorimetric Detection of the Conical Terahertz Radiation from Femtosecond Laser Filaments in Air. Applied Physics Letters, 91, Article ID: 241105. https://doi.org/10.1063/1.2821371
[35]	González de Alaiza Martínez, P., Davoine, X., Debayle, A., Gremillet, L. and Berge, L. (2016) Terahertz Radiation Driven by Two-Color Laser Pulses at Near-Relativistic Intensities: Competition between Photoionization and Wakefield Effects. Scientific Reports, 6, Article No. 26743. https://doi.org/10.1038/srep26743
[36]	Oh, T.I., You, Y.S. Jhajj, N., Rosenthal, E.W., Milchberg, H.M. and Kim, K.Y. (2013) Intense Terahertz Generation in Two-Color Laser Filamentation: Energy Scaling with Terawatt Laser Systems. New Journal of Physics, 15, Article ID: 075002. https://doi.org/10.1088/1367-2630/15/7/075002
[37]	Wang, T.-J., Yuan, S., Chen, Y.P., Daigle, J.-F., Marceau, C., Thberge, F., Chteauneuf, M., Dubois, J. and Chin, S.L. (2010) Toward Remote High Energy Terahertz Generation. Applied Physics Letters, 97, Article ID: 111108. https://doi.org/10.1063/1.3490702
[38]	Chen, G., Kawayama, I., Murakami, H., Teramoto, T. and Tonouchi M. (2021) Intensity-Dependent Self-Induced Dual-Color Laser Phase Modulation and its Effect on Terahertz Generation. Scientific Reports, 11, Article No. 498. https://doi.org/10.1038/s41598-020-80105-7