High-Order Spatial FDTD Solver of Maxwell’s Equations for Terahertz Radiation Production

Abdelrahman Mahdy

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

DOI: 10.4236/jamp.2024.124063                

Abstract

We applied a spatial high-order finite-difference-time-domain (HO-FDTD) scheme to solve 2D Maxwell’s equations in order to develop a fluid model employed to study the production of terahertz radiation by the filamentation of two femtosecond lasers in air plasma. We examined the performance of the applied scheme, in this context, we implemented the developed model to study selected phenomena in terahertz radiation production, such as the excitation energy and conversion efficiency of the produced THz radiation, in addition to the influence of the pulse chirping on properties of the produced radiation. The obtained numerical results have clarified that the applied HO-FDTD scheme is precisely accurate to solve Maxwell’s equations and sufficiently valid to study the production of terahertz radiation by the filamentation of two femtosecond lasers in air plasma.

Keywords

The Finite-Difference-Time-Domain, Terahertz Radiation Production, Filamentation of Femtosecond Laser, Maxwell’s Equations Solution

1. Introduction

The finite-difference-time-domain (FDTD) is the principal numerical method in the electromagnetic waves (EMWs) propagation in materials research [1] . Although the fourth stage Runge-Kutta integrator, the ADI-FDTD [2] [3] , the septicemic FDTD [4] , the discontinuous Crank-Nicholson [5] , and the quasi-static FDTD (QS-FDTD) [6] are widely considered, the conventional second order in both of space and time algorithm FDTD (2, 2) was the most predominant FDTD algorithm for solving Maxwell’s equations [7] in this research. As a result of laser invention and due to development of new materials, there has been an increasing demand to study the EMW or laser propagation in ultrashort time and spatial domain in a complex geometry for a long propagation time. In this new computational arrangement, the conventional FDTD (2, 2) has shown a negative feedback in simulating such changeling computational problems that are represented with more refinement propagation conditions and complex materials structure. It has been determined [8] that the FDTD (2, 2) generates errors in form of computational dispersion anisotropy, and these errors are accumulating and enlarging as computational time is proceeding. For that reason, the FDTD (2, 2) is considered insufficiently accurate, accordingly, new high-order FDTD algorithms have to be developed.

The high-order FDTD (HO-FDTD) means higher than two orders in both space and time. From this perspective, the explicit staggered FDTD (2, 4) was first developed for frequency-dependent medium (dispersive medium), then the FDTD (3, 4) was added, and later a FDTD (4, 4) for the frequency-independent medium (non-dispersive) was proposed [9] . After it has been verified [10] and demonstrated notable advancements over the conventional FDTD (2, 2), many HO-FDTD schemes based on these three intrinsic HO-FDTD [11] algorithms have been developed. For instance, some HO-FDTD schemes were at first introduced for dispersive mediums that are associated with nonlinear laser-mater interaction effects (Kerr-effects) [12] , others were efficiently applied to solve Maxwell’s equations for perfect conduction periodic exploration and interpolation periodic boundaries, then more accurate HO schemes were more developed for the perfectly matched layer [13] , which is independent on the frequency, polarization and the incidence/reflection conditions. In principle, the developed HO-FDTD schemes are engaged to study the contemporary physics problems, as they can take any spatial grid sizes of the modern physical problems (1/10, 1/20, ....) for the wavelength of the applied EMW. Although these developed HO-FDTD schemes require more computational time than the conventional FDTD (2, 2) at the same grid size, it present a compromise solution between the computational cost and accuracy with an adequate dispersion for the real problems.

In this article, we applied a HO spatial FDTD scheme (HO-FDTD) [14] to develop a 2D fluid model employed to study the terahertz radiation production by the filamentation of two femtosecond laser beams in air plasma. In this model, we solved the fluid equations that include the equation of motion, the ionization equation, and the Maxwell’s equations, note, that the HO-FDTD scheme is mainly applied to solve Maxwell’s equations of this model. The article is organized as follows: Section 2 is devoted to describing the basic equations of our fluid model with a main focus on Maxwell’s equations and its physical assumptions and considerations. In Section 3, we apply the HO-FDTD scheme to solve Maxwell’s equations, the numerical procedures of this solution are described in this section, and in addition, the initial and boundary conditions used in this solution are also listed in this section. In Section 4, we examine the performance of the applied scheme, as we implement the developed 2D model to study selected phenomena and processes in the terahertz radiation production, numerical results for this phenomenon and processes are presented in this section.

2. Basic Equations

In the production of terahertz radiation by two femtoseconds (fs) laser beams in air plasma, the filamentation of a fundamental first harmonic (FH) laser pulse ( ${\omega }_0$ ) and its corresponding second harmonic (SH) pulse ( $2{\omega }_0$ ) is the fundamental mechanism of the production. Two theories govern the THz production in connection with this mechanism, namely the four-wave mixing rectification (FWM) [15] and the transient photocurrent theory (PC) [16] . Herein, we study the production of THz radiation by the filamentation of two femtosecond lasers in air plasma based on the PC theory [17] . In that respect, a 2D fluid model is implemented to conduct this study, the moment equation and the ionization equation coupled with a 2D Maxwell’s equation are the basic equations of this model. In these equations, the moment equation is the standard equation of motion [18] for non-relativistic electrons dynamics in collisional cold plasma, and the ionization equation is the classical ADK tunneling ionization rate formula [19] that maintains the time-dependent ionization process after the threshold ionization, with respect to our central equations; the Maxwell’s equation, it 2D equations that are occupied to describe the fs beam propagation and filamentation in dielectric, dispersive, lossless, and isotropic media that is unchanged during the filamentation.

As a normal procedure in the theoretical calculations for both of the electric and magnetic fields in the ultra short laser plasma interaction using Maxwell’s equations [20] , we take the divergence of Gauss’s laws in combination with Faraday and Ampere-Maxwell’s laws to obtain these fields. According to the charge conversation nature of filamentation process of the femtosecond beam in air plasma, Gauss’s laws always satisfy the initial conditions in this process, and at this nature, Faraday and Ampere-Maxwell’s laws are fully described the solution of these fields that are remained always satisfied by Gauss’s laws. At this circumstance, solution of Faraday and Ampere’s laws are sufficient to describe the electric and magnetic fields evolution in this interaction, upon that our basic Maxwell’s equations are limited to the following equations:

$\nabla \times E=-\frac{\partial B}{\partial t}$ , Faraday’s Law (1)

$\nabla \times B=\mu \epsilon \frac{\partial E}{\partial t}$ , Ampere’s Law (2)

where E and B are the electric and magnetic field, and $\epsilon \mathrm{,}\mu$ are the dielectric permittivity and the susceptibility of the air plasma, receptively.

TE Maxwell’s Equations

To simplify applying the HO-FDTD scheme to solve our modified Maxwell’s equations, we considered an incident transverse electric (TE) polarization plane wave with the field components ( $E_x\mathrm{,}{E}_y\mathrm{,}{B}_z$ ) in our study. Under this consideration, our basic Maxwell’s equations can be re-written as:

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

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

$\frac{\partial B_z}{\partial t}=\left(\frac{\partial E_x}{\partial y}-\frac{\partial E_y}{\partial x}\right)$ , Faraday’s Law (5)

3. The High-Order Spatial FDTD Scheme

In our fluid model employed to study the production of THz radiation by the filamentation of two fs beams in air plasma, the moment equation is numerically solved using the standard leapfrog integrator [21] and the ionization equation is numerically solved using standard explicit particle-in-cell method [22] , for the Maxwell’s equations given in Equations ((3), (4)) and Equation (5), the HO-FDTD scheme is applied to solve these equations, the numerical procedures of the solution are fully described in this section.

3.1. The Initial Conditions

In the scope of our investigating for the production of THz radiation in air plasma, the appropriate initial conditions are stated by the initial values of the applied fields $\left(E^0\mathrm{,}{B}^0\right)$ at $t=0$ for any $\left(x\mathrm{,}y\right)$ value in the 2D spatial domain $\Omega =\left[\mathrm{0,}a\right]\times \left[\mathrm{0,}b\right]$ , where $a=I{\Lambda }_x$ , $b=J{\Lambda }_y$ , and ${\Lambda }_x\mathrm{,}{\Lambda }_y$ are the spatial step, $I\mathrm{,}J$ are the maximum spatial step along $x\mathrm{,}y$ , respectively, thus the initial conditions in our study are given by:

$E\left(x\mathrm{,}y\mathrm{,0}\right)=E^0\left(x\mathrm{,}y\right)=\left(E_x^0\left(x\mathrm{,}y\right)\mathrm{,}{E}_y^0\left(x\mathrm{,}y\right)\right)\mathrm{,}$ (6)

and

$B_z\left(x\mathrm{,}y\mathrm{,0}\right)=B_z^0\left(x\mathrm{,}y\right)\text{for}\left(x\mathrm{,}y\right)\in \Omega \mathrm{.}$ (7)

where $E_x^0\mathrm{,}{E}_y^0$ are the initial electric field components and $B_z^0$ is the initial magnetic field component of the applied pulse amplitude. Regardless of the initial conditions of the applied fields in the time or the spatial domain, these fields should be constrained with the conservation laws as:

$\nabla \cdot E=0\text{and}\nabla \cdot {\partial }_{t}E=0,$

$\nabla \cdot B=0\text{and}\nabla \cdot {\partial }_{t}B=0.$

3.2. The Numerical Procedures

Starting at the time step $n=0$ , and based on the initial values given for $E_x^0$ and $B_z^0$ in Equations ((6), (7)).

3.2.1. Compute $E_x^{\ast }$ and $B_z^{\ast }$

We obtain $E_x^{\mathrm{<mo>\; *\; </mo>}}$ and $B_z^{\mathrm{<mo>\; *\; </mo>}}$ as following:

· For the interior nodes, $i=0,1,2,\cdots ,I-1$ ,

$\frac{E_{x_{i+\frac{1}{2},j}}^{*}-E_{x_{i+\frac{1}{2},j}}^n}{\Delta t}=\frac{\mu \epsilon }{4}{\Lambda }_y\left\{B_{z_{i+\frac{1}{2},j}}^{*}+B_{z_{i+\frac{1}{2},j}}^n\right\},j=2,3,\cdots ,J-2$ (8)

$\frac{B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{*}-B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^n}{\Delta t}=\frac{1}{4}{\Lambda }_y\left\{E_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^{*}+E_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^n\right\},j=1,2,\cdots ,J-2$ (9)

· For the near boundaries nodes, $i=0,1,2,\cdots ,I-1$ ,

$\frac{E_{x_{i+\frac{1}{2},j}}^{*}-E_{x_{i+\frac{1}{2},j}}^n}{\Delta t}=\frac{\mu \epsilon }{4}{\bar{\Lambda }}_y\left\{B_{z_{i+\frac{1}{2},j}}^{*}+B_{z_{i+\frac{1}{2},j}}^n\right\},j=1\text{and}j-1$ (10)

$\frac{B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{*}-B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^n}{\Delta t}=\frac{1}{4}{\bar{\Lambda }}_y\left\{E_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^{*}+E_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^n\right\},j=0\text{and}j-1$ (11)

· Where, for $i=0,1,2,\cdots ,I-1$ and $j=\mathrm{2,}\cdots ,J-2$ ,

${\Lambda }_y{B}_{z_{i+\frac{1}{2},j}}^{*}=\frac{1}{8}\left(9{\delta }_y-{\delta }_{2,y}\right)B_{z_{i,j}}^{*},{\Lambda }_y{B}_{z_{i+\frac{1}{2},j}}^n=\frac{1}{8}\left(9{\delta }_y-{\delta }_{2,y}\right)B_{z_{i,j}}^n,$

and for $i=0,1,2,\cdots ,I-1$ and $j=1,\mathrm{2,}\cdots ,J-2$ ,

${\Lambda }_y{E}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^{*}=\frac{1}{8}\left(9{\delta }_y-{\delta }_{2,y}\right)E_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^{*},{\Lambda }_y{E}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^n=\frac{1}{9}\left(9{\delta }_y-{\delta }_{2,y}\right)E_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^n,$

· While for $i=0,1,\cdots ,I-1$ and $j=J-1$ ,

${\bar{\Lambda }}_y{B}_{z_{i+\frac{1}{2},j}}^{*}=\frac{1}{8}\left(9{\delta }_y-{\bar{\delta }}_{2,y}\right)B_{z_{i,j}}^{*},{\Lambda }_y{B}_{z_{i+\frac{1}{2},j}}^n=\frac{1}{8}\left(9{\delta }_y-{\delta }_{2,y}\right)B_{z_{i,j}}^n,$

and

${\bar{\Lambda }}_y{E}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^{*}=\frac{1}{8}\left(9{\delta }_y-{\bar{\delta }}_{2,y}\right)E_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{*},{\Lambda }_y{E}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^n=\frac{1}{9}\left(9{\delta }_y-{\delta }_{2,y}\right)E_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^n,$

· In general, we defined the operators ${\delta }_t\mathrm{,}{\delta }_x\mathrm{,}{\delta }_y,{\delta }_{\mathrm{2,}x}$ and ${\delta }_{\mathrm{2,}y}$ as:

${\delta }_t{U}_{i,j}^n=\frac{U_{i,j}^{n+\frac{1}{2}}-U_{i,j}^{n-\frac{1}{2}}}{\Delta t},{\delta }_x{U}_{i,j}^n=\frac{U_{i+\frac{1}{2},j}^n-U_{i-\frac{1}{2},j}^n}{\Delta x},{\delta }_y{U}_{i,j}^n=\frac{U_{i,j+\frac{1}{2}}^n-U_{i,j-\frac{1}{2}}^n}{\Delta y},$

and

${\delta }_{2,x}{U}_{i,j}^n=\frac{U_{i+\frac{3}{2},j}^n-U_{i-\frac{3}{2},j}^n}{3\Delta x},{\delta }_{2,y}{U}_{i,j}^n=\frac{U_{i,j+\frac{3}{2}}^n-U_{i,j\text{}--\frac{3}{2}}^n}{3\Delta y},{\delta }_u{\delta }_u{U}_{i,j}^n={\delta }_u\left({\delta }_u{U}_{i,j}^n\right)$

From Equations ((8), (9)), for the interior nodes, $i=0,1,2,\cdots ,I-1$ ,

$E_{x_{i+\frac{1}{2},j}}^{*}=\left(\frac{\mu \epsilon }{4}{\Lambda }_y\left\{B_{z_{i+\frac{1}{2},j}}^{*}+B_{z_{i+\frac{1}{2},j}}^n\right\}\right)\Delta t+E_{x_{i+\frac{1}{2},j}}^{n}j=2,3,\cdots ,J-2$ (12)

and

$B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{*}=\left(\frac{1}{4}{\Lambda }_y\left\{E_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^{*}\text{}\text{}+E_{x_{i+\frac{1}{2},j\text{}\text{}+\frac{1}{2}}}^n\right\}\right)\Delta t+B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^n,j=1,2,3,\cdots ,J-2$

While Equations ((10), (11)), for the near boundaries nodes, $i=0,1,2,\cdots ,I-1$ ,

$E_{x_{i+\frac{1}{2},j}}^{*}=\left(\frac{\mu \epsilon }{4}{\Lambda }_y\left\{B_{z_{i+\frac{1}{2},j}}^{*}+B_{z_{i+\frac{1}{2},j}}^n\right\}\right)\Delta t+E_{x_{i+\frac{1}{2},j}}^n,j=1\text{and}J-1$ (13)

and

$B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{*}=\left(\frac{1}{4}{\bar{\Lambda }}_y\left\{E_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^{*}+E_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^n\right\}\right)\Delta t+B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^n,j=0\text{and}J-1$ (14)

3.2.2. Compute $E_y^{n+1}$ and $B_z^{\ast \ast }$

After we obtained $E_x^{\mathrm{<mo>\; *\; </mo>}}$ and $B_z^{\mathrm{<mo>\; *\; </mo>}}$ in the previous step, we calculate $E_y^{n+1}$ and $B_z^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}$ by discretize Equations ((4), (5)) as the following:

· For the interior nodes, $j=\mathrm{0,1,2,}\cdots ,J-1$ ,

$\frac{E_{y_{i,j+\frac{1}{2}}}^{n+1}-E_{y_{i,j+\frac{1}{2}}}^n}{\Delta t}=-\frac{\mu \epsilon }{2}{\Lambda }_x\left\{B_{z_{i,j+\frac{1}{2}}}^{**}+B_{z_{i,j+\frac{1}{2}}}^n\right\},i=2,3,\cdots ,I-2$ (15)

$\frac{B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{**}-B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{*}}{\Delta t}=\frac{1}{2}{\Lambda }_x\left\{E_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+E_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^n\right\},i=1,2,\cdots ,I-2$ (16)

· For the near boundary, $j=\mathrm{0,1,2,}\cdots ,J-1$ ,

$\frac{E_{y_{i,j+\frac{1}{2}}}^{n+1}-E_{y_{i,j+\frac{1}{2}}}^n}{\Delta t}=-\frac{\mu \epsilon }{2}{\bar{\Lambda }}_x\left\{B_{z_{i,j+\frac{1}{2}}}^{**}+B_{z_{i,j+\frac{1}{2}}}^n\right\},i=1\text{and}i=I-1$ (17)

$\frac{B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{**}-B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{*}}{\Delta t}=\frac{1}{2}{\bar{\Lambda }}_x\left\{E_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+E_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^n\right\},i=0\text{and}i=I-1$ (18)

· Where for $i=0,1,2,\cdots ,I-1$ ; $j=2,\cdots ,J-2$ ,

${\Lambda }_x{B}_{z_{i,j+\frac{1}{2}}}^{**}=\frac{1}{9}\left(9{\delta }_x-{\delta }_{2,x}\right)B_{z_{i,j}}^{**},{\Lambda }_x{B}_{y_{i,j+\frac{1}{2}}}^{*}=\frac{1}{9}\left(9{\delta }_x-{\delta }_{2,x}\right)B_{z_{i,j}}^{*}$

and

${\Lambda }_x{E}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1}=\frac{1}{9}\left(9{\delta }_x-{\delta }_{2,x}\right)E_{y_{i+\frac{1}{2},j}}^{n+1},{\Lambda }_x{E}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^n=\frac{1}{9}\left(9{\delta }_x-{\delta }_{2,x}\right)E_{y_{i\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^n$

· While for $i=1$ and $j=0,1,2,\cdots ,J-1$ ,

${\bar{\Lambda }}_x{B}_{z_{i,j+\frac{1}{2}}}^{**}=\frac{1}{9}\left(9{\delta }_x-{\bar{\delta }}_{2,x}\right)B_{z_{i,j+\frac{1}{2}}}^{**},{\bar{\Lambda }}_x{B}_{y_{i,j+\frac{1}{2}+\frac{1}{2}}}^{*}=\frac{1}{9}\left(9{\delta }_x-{\bar{\delta }}_{2,x}\right)B_{z_{i,j+\frac{1}{2}}}^{*}$

· And for $i=0$ and $i=I-1$ ,

${\bar{\Lambda }}_x{E}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1}=\frac{1}{9}\left(9{\delta }_x-{\bar{\delta }}_{2,x}\right)E_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1},{\bar{\Lambda }}_x{E}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^n=\frac{1}{9}\left(9{\delta }_x-{\bar{\delta }}_{2,x}\right)E_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^n$

· From Equations ((15), (16)), for the interior nodes, $j=0,1,2,\cdots ,J-1$ ,

$E_{y_{i,j+\frac{1}{2}}}^{n+1}=\left(-\frac{\mu \epsilon }{2}{\Lambda }_x\left\{B_{z_{i,j+\frac{1}{2}}}^{**}+B_{z_{i,j+\frac{1}{2}}}^n\right\}\right)\Delta t+E_{y_{i,j+\frac{1}{2}}}^n,i=2,3,\cdots ,I-1$ (19)

$B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{**}=\left(\frac{1}{2}{\Lambda }_x\left\{E_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+E_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^n\right\}\right)\Delta t+B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{*},i=1,3,\cdots ,I-2$ (20)

• From Equations ((17), (18)) for the near boundary, $j=\mathrm{0,1,2,}\cdots ,j-1$ ,

$B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{**}=\left(\frac{1}{2}{\bar{\Lambda }}_x\left\{E_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+E_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^n\right\}\right)\Delta t-B_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{*},i=0\text{and}i=I-1$ (21)

$E_{y_{i,j+\frac{1}{2}}}^{n+1}=\left(-\frac{\mu \epsilon }{2}{\bar{\Lambda }}_x\left\{B_{z_{i,j+\frac{1}{2}}}^{**}+B_{z_{i,j+\frac{1}{2}}}^n\right\}\right)\Delta t+E_{y_{i,j+\frac{1}{2}}}^n,i=1\text{and}i=I-1$ (22)

3.2.3. Compute $E_x^{n+1}$ and $B_z^{n+1}$

Finally, we obtain the $E_x^{n+1}$ and $B_z^{n+1}$ based on the calculated variables $E_y^{n+1}$ and $B_z^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}$ in the previous step, as:

· For the interior nodes, $i=0,1,2,\cdots ,I-1$ ,

$\frac{E_{x_{i+\frac{1}{2}},j}^{n+1}-E_{x_{i+\frac{1}{2},j}}^n}{\Delta t}=\frac{\mu \epsilon }{4}{\Lambda }_x\left\{B_{z_{i+\frac{1}{2},j}}^{n+1}+B_{z_{i+\frac{1}{2},j}}^{**}\right\},j=2,3,\cdots ,J-2$ (23)

$\frac{B_{z_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{n+1}-B_{z_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{**}}{\Delta t}=\frac{1}{4}{\Lambda }_x\left\{E_{x_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+E_{x_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{*}\right\},j=1,2,\cdots ,J-2$ (24)

· For the near boundary $i=0,1,2,\cdots ,I-1$ ,

$\frac{E_{x_{i+\frac{1}{2}},j}^{n+1}-E_{x_{i+\frac{1}{2},j}}^n}{\Delta t}=\frac{\mu \epsilon }{4}{\bar{\Lambda }}_x\left\{B_{z_{i+\frac{1}{2},j}}^{n+1}+B_{z_{i+\frac{1}{2},j}}^{**}\right\},j=1\text{and}j-1$ (25)

$\frac{B_{z_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{n+1}-B_{z_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{**}}{\Delta t}=\frac{1}{4\mu }{\bar{\Lambda }}_x\left\{E_{x_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+E_{x_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{*}\right\},j=0\text{and}j-1$ (26)

· From Equations ((23), (24)), for the interior nodes, $i=0,1,2,\cdots ,I-1$ ,

$E_{x_{i+\frac{1}{2}},j}^{n+1}=\left(\frac{\mu \epsilon }{4}{\Lambda }_x\left\{B_{z_{i+\frac{1}{2},j}}^{n+1}+B_{z_{i+\frac{1}{2},j}}^{**}\right\}\right)\Delta t+E_{x_{i+\frac{1}{2},j}}^n,j=2,3,\cdots ,J-2$ (27)

$B_{z_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{n+1}=\left(\frac{1}{4}{\Lambda }_x\left\{E_{x_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+E_{x_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{*}\right\}\right)\Delta t+B_{z_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{**},j=1,2,3,\cdots ,J-2$ (28)

· From Equations ((25), (26)), for the near boundary nodes, $i=0,1,2,\cdots ,I-1$ ,

$E_{x_{i+\frac{1}{2},j}}^{n+1}=\left(\frac{\mu \epsilon }{4}{\bar{\Lambda }}_x\left\{B_{z_{i+\frac{1}{2},j}}^{n+1}+B_{z_{i+\frac{1}{2},j}}^{**}\right\}\right)\Delta t-E_{x_{i+\frac{1}{2}},j}j=1\text{and}j-1$ (29)

$B_{z_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{n+1}=\left(\frac{1}{4}{\bar{\Lambda }}_x\left\{E_{x_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+E_{x_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{*}\right\}\right)\Delta t+B_{z_{i+\frac{1}{2}+\frac{1}{2},j+\frac{1}{2}}}^{**},j=0\text{and}j-1$ (30)

We increment the time step n by one, and repeat the procedures of Sections 3.2.1-3.2.3 to the end of the simulation time T to obtain the electric and magnetic fields over the total time domain of the problem.

3.3. Boundary Conditions

Due the physical nature of our problem, some considerations have to be taken into account for the boundary conditions in the applied numerical procedures. Initially, the air plasma target in an air-dielectric-air 2D slab, this slab is filled with a uniform air plasma gas, in addition, the air-dielectric boundary is perfectly conducting interface where the applied electric field is permitted to propagate without absorption, while the dielectric-air boundary is reflected interface where the electric field is sum of the reflected and transmitted field, moreover, the energy is conserved in both interfaces. Based on these considerations, the time-derivatives $\left(\partial E_x/\partial t\mathrm{,}\partial E_y/\partial t\mathrm{,}\partial B_z/\partial t\right)$ and the tangential derivatives $\left(\partial E_x/\partial y\mathrm{,}\partial E_y/\partial x\mathrm{,}\partial B_z/\partial y\mathrm{,}\partial B_z/\partial x\right)$ of the fields are continuous at these interfaces.

4. Numerical Examples

After we have point by point explained the numerical procedures of solving our Maxwell’s equations using the HO-FDTD scheme, sequentially, it will be necessary to examine the performance of this scheme. In this section, we employed the developed 2D fluid model to study selected fundamental phenomena and processes in THz radiation production by the filamentation of two fs beams in air plasma, in particular optimizing the laser parameters for an efficient THz production, the conversion efficiency of the produced THz radiation in some noble gases, a comparison for the efficiency evolution behavior among these gases, and finally the chirping of the input laser pulse and its influence on the properties of the produced THz radiation. In our examination, a 2D plasma slab of equal dimension $x=y=50\text{-}100\mu \text{m}$ is considered, this slab is filled with plasma at initial electron density $n_e=3.7\times {10}^{19}{\text{cm}}^{-3}$ . The applied laser beam is a TE plane wave that is normally incident on the entrance interface at $x=0$ , the electric fields $\left(E_x\mathrm{,}{E}_y\right)$ of this wave are calculated in this slab, while the magnetic field $B_z$ is derived. The diagnostic and the measurements for the selected fundamental phenomena and processes are carried out after the the applied beam pass the exit interface at $x=a$ by $d_0=3\text{-}10\mu \text{m}$ . The initial combined two-beams profile in this investigation is given by [23] :

$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{.}$ (31)

where $E_m=\sqrt{2I_0/{\epsilon }_{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.

4.1. Optimizing the Suitable Beam Parameters

It necessary to bear in mind that the initial values of the applied beam parameters defined in Equation (31) are key players in our research. Concerning these values, it should be aware that in the most of our research parts, the initial wavelength ( ${\lambda }_0=800\text{nm}$ ) and the initial pulse duration ${\tau }_0=35\text{fs}$ are fundamentally given, while the initial input intensity is the dependent parameter where most of the variables are investigated and analyzed as function on, its given a period of $I={10}^{14}\text{-}{10}^{17}$ W/cm². With respect to the relative phase $\varphi$ and energy fraction factor $\zeta$ , truly the initial values of these two parameters are the investigated problem dependent, in other words, $\varphi$ and $\zeta$ can’t be given, but its suitable values should be estimated and optimized for the problem of our interest. In this section, we optimized the suitable relative phase and energy fraction factor, in this context, we simulated the spectrum of the excitation energy ${\mu }_{THz}$ of the produced THz radiation, this simulation swapped over the whole relative phase $\varphi =0\text{-}\pi$ and energy fraction factor $\zeta =0\text{-}1$ possible values, and in order to cover the given initial intensity period, we presented this spectrum for three initial intensity values, $I_0=8.0\times {10}^{15}$ W/cm², $I_0=3.0\times {10}^{16}$ W/cm² and $I_0=1.0\times {10}^{17}$ W/cm², the results of the ${\mu }_{THz}$ spectrum at these intensity values are displayed in Figure 1.

Figure 1. The excitation energy ${\mu }_{THz}$ of the produced THz radiation at different input intensities.

As clearly illustrated in this spectrum, the maximum ${\mu }_{THz}$ is induced at $\varphi \approx \pi /2$ and $\zeta \approx 0.5$ for $I_0=8.0\times {10}^{15}$ W/cm² and $I_0=3.0\times {10}^{16}$ W/cm², although the maximum ${\mu }_{THz}$ is also induced at $\varphi =\pi /5$ and $\zeta =0.7$ for $I_0=1.0\times {10}^{17}$ W/cm², $\zeta =0.7$ value is excluded in our research to avoid the possible propagation instabilities that are emerged at the unequal input intensity. The suitable $\varphi$ and $\zeta$ obtained in this section are typically compatible with the standard suitable ( $\varphi =\pi /2$ , $\zeta =0.5$ ) values determined in the preliminary THz research [16] , therefore, this typical compatibility highlights on the good performance of the applied HO-FDTD scheme to study the production of the THz radiation.

4.2. The Conversion Efficiency of the Produce THz Radiation

The optical-to-terahertz conversion is the substantial process in THz radiation production, and the efficiency of this conversion is the main scale that characterizes the properties of produced THz radiation. The efficiency

${\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$ , is defined [24] as the ratio of the generated terahertz intensity ${\left|E_{THz}\left(\omega \right)\right|}^2$ that accumulated over the THz period to the applied beam intensity ${\left|E_0\left(r\mathrm{,}t\right)\right|}^2$ that is summed over the total simulation time. To study the conversion efficiency, we calculated this efficiency ${\eta }_{THz}$ as function of the input intensity of the applied beam over the given intensity period. Three noble gases; which are Neon (Ne), Argon (Ar), and Xenon (Xe) at different binding energies, are selected to conduct this calculation, the binding energies of these gases [25] are listed in Table 1 and the results of the calculation is presented in Figure 2.

As clearly demonstrated in Figure 2, within the intensity period of few 10¹⁴ W/cm², the efficiency measurements is initiated, then it is commonly enlarging with increasing the intensity for all of the gases. At $I_0=2.0\times {10}^{15}\text{W}/{\text{cm}}^{\text{2}}$ , the efficiency reaches its maximum value for Xe, at $I_0=4\times {10}^{15}\text{W}/{\text{cm}}^{\text{2}}$ for Ar, and at $I_0=9.0\times {10}^{15}\text{W}/{\text{cm}}^{\text{2}}$ for Ne. Shortly after, the efficiency is commonly decaying for all the gases with increasing the intensity to reach its lowest value within the input intensity period $I_0=2.9\text{-}4\times {10}^{16}\text{W}/{\text{cm}}^{\text{2}}$ , afterwards the efficiency is saturating for all of these gases to the end of the given intensity period.

The efficiency-evolution behavior with its associated sequential dynamics, i.e. the measurements initiation, enlarging, maximize, decaying, and saturation, demonstrated in Figure 2 is a standard evolution behavior for the efficiency as function on the input intensity of the applied beams in the THz production researches [26] . As this behavior and its associated dynamics is similarly behaved in our research, it can be concluded that the applied HO-FDTD scheme is valid

Table 1. The binding energy (eV) of some noble gases.

Figure 2. The efficiency of the produced THz as function of the input intensity for different noble gases.

to study the conversion efficiency.

The efficiency-evolution behavior is not the only valid phenomena that is demonstrated in Figure 2. As depicted in this figure, in a comparison analysis for the efficiency values obtained for the selected gases, it is clearly noted that the efficiency value obtained for Xe; of the lowest binding, energy is greater than the one obtained for Ar which is higher than the value obtained for Ne. As a matter of fact, increasing the efficiency with decreasing the binding energy is a predominant event in the filamentation of fs in plasma [27] . This occurrence is initially observed in an experimental study [28] and then it is approved in a theoretical calculation [29] . Now, it is additionally affirmed in our study to ensure the validity of the applied HO-FDTD scheme for our THz study.

4.3. Pulse Chirping Influence on the Produced THz Radiation

Laser pulse chirping [30] is changing the frequency/wavelength during the duration of this pulse. The pulse chirping is classified [31] into: 1) positive pulse chirping which implies increasing the frequency of the propagating pulse during its duration, and 2) negative chirping which points to decreasing the frequency of the pulse during this duration. In the filamentation of fs in air plasma, the pulse chirping phenomena offers a complete control on the spatially-temporal properties of the optical pulse, thus in accordance with this control, more efficient optical-to-terahertz conversion is achieved. To examine the performance of the HO-FDTD scheme in studying the pulse chirping, essentially we replaced the initial profile given in Equation (31) with a particular initial combined pulse profile [18] where the pulse chirping with its positive and negative properties can be established, following that we calculated the induced THz yield for positively and negatively pulse chirping at different pulse durations, the results of this calculation are displayed in Figure 3.

As seen in Figure 3, the induced THz yield reaches its maximum value at chirped pulse duration of around 125 fs for both the positive and negative chirping, meanwhile, in this figure, the maximum value obtained for the positive chirping is higher than for the negative chirping pulse. The disadvantage of the negative chirping in reaching a higher THz yield than the positive chirping is widely-recognized characteristic in the filamentation process [26] that is now verified in Figure 3. The interpretation for this disadvantage is that the filamentation for the negative chirping stars earlier, which gives rise to some destructive filament competition.

5. Conclusion

We developed a 2D fluid model to study the terahertz radiation production by the filamentation of two femtosecond laser beams in air plasma. A spatially high-order finite-difference-time-domain (HO-FDTD) scheme was applied to solve the Maxwell’s equations of this model. Extensive numerical results have approved the performance of the applied scheme to study the terahertz radiation production. As an illustration, in the first place, the scheme has exactly estimated the suitable values of some beam parameters of an efficient THz production, in addition, it has precisely monitored the conventional efficiency-evolution behavior

Figure 3. The THz radiation yield as function of the pulse duration of +ve and −ve pulse chirping.

and its associated dynamics in some noble gases, moreover, the applied scheme has justified the experimental and theoretical event of increasing the efficiency of the produced THz radiation with decreasing the binding energies of these gases. Lastly, the HO-FDTD scheme has closely demonstrated the well-known influence of the positive and negative chirping on the yield of the produced THz radiation.

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]	Taflove, A. and Hagness, S.C. (2005) Computational Electrodynamics: The Finite-Difference Time-Domain Method. Artech House, Norwood. https://doi.org/10.1016/B978-012170960-0/50046-3
[2]	Namiki, T. (1999) A New FDTD Algorithm Based on Alternating-Direction Implicit Method. IEEE Transactions on Microwave Theory and Techniques, 47, 2003-2007. https://doi.org/10.1109/22.795075
[3]	Zhen, F.H., Chen, Z.Z. and Zhang, J.Z. (2000) Toward the Development of a Three-Dimensional Unconditionally Stable Finite-Difference Time-Domain Method. IEEE Transactions on Microwave Theory and Techniques, 48, 1550-1558. https://doi.org/10.1109/22.869007
[4]	Hirono, T., Lui, W., Seki, S. and Yoshikuni, Y. (2001) A Three-Dimensional Fourth-Order Finite-Difference Time-Domain Scheme Using a Symplectic Integrator Propagator. IEEE Transactions on Microwave Theory and Techniques, 49, 1640-1648. https://doi.org/10.1109/22.942578
[5]	Li, J.C. and Hesthaven, J.S. (2014) Analysis and Application of the Nodal Discontinuous Galerkin Method for Wave Propagation in Metamaterials. Journal of Computational Physics, 258, 915-930. https://doi.org/10.1016/j.jcp.2013.11.018
[6]	Kim, M., Park, S. and Jung, H.-K. (2018) An Advanced Numerical Technique for a Quasi-Static Electromagnetic Field Simulation Based on the Finite-Difference Time-Domain Method. Journal of Computational Physics, 373, 917-923. https://doi.org/10.1016/j.jcp.2018.07.032
[7]	Yee, K. (1966) Numerical Solution of Initial Boundary Value Problems Involving Maxwell’s Equations in Isotropic Media. IEEE Transactions on Antennas and Propagation, 14, 302-307. https://doi.org/10.1109/TAP.1966.1138693
[8]	Sha, W., Huang, Z.X., Chen, M.S. and Wu, X.L. (2008) Survey on Symplectic Finite-Difference Time-Domain Schemes for Maxwell’s Equations. IEEE Transactions on Antennas and Propagation, 56, 493-500. https://doi.org/10.1109/TAP.2007.915444
[9]	Fang, J. (1989) Time Domain Finite Difference Computation for Maxwell’s Equations. Ph.D. Dissertation, University of California, Berkeley.
[10]	Hwang, K.-P. (2003) Computational Efficiency of Fangs Fourth-Order FDTD Schemes. Electromagnetics, 23, 89-102. https://doi.org/10.1080/02726340390159450
[11]	Hesthaven, J.S. (2003) High-Order Accurate Methods in Time-Domain Computational Electromagnetics: A Review. In: Advances in Imaging and Electron Physics, Vol. 127, Elsevier, Amsterdam, 59-123. https://doi.org/10.1016/S1076-5670(03)80097-6
[12]	Bokil, V.A. and Gibson, N.L. (2011) Analysis of Spatial High-Order Finite Difference Methods for Maxwell’s Equations in Dispersive Media. IMA Journal of Numerical Analysis, 32, 926-956. https://doi.org/10.1093/imanum/drr001
[13]	Wei, X.-K., Shao, W., Ou, H.Y. and Wang, B.-Z. (2016) An Efficient Higher-Order PML in WLP-FDTD Method for Time Reversed Wave Simulation. Journal of Computational Physics, 321, 1206-1216. https://doi.org/10.1016/j.jcp.2016.06.032
[14]	Liang, D. and Yuan, Q. (2013) The Spatial Fourth-Order Energy-Conserved S-FDTD Scheme for Maxwell’s Equations. Journal of Computational Physics, 243, 344-364. https://doi.org/10.1016/j.jcp.2013.02.040
[15]	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
[16]	Kim, K.Y., Taylor, A.J., Glownia, J.H. and Rodriguez, G. (2008) Coherent Control of Terahertz Supercontinuum Generation in Ultrafast Laser-Gas Interactions. Nature Photonics, 2, 605-609. https://doi.org/10.1038/nphoton.2008.153
[17]	Kim, K.-Y. (2009) Generation of Coherent Terahertz Radiation in Ultrafast Laser-Gas Interactions. Physics of Plasmas, 16, Article ID: 056706. https://doi.org/10.1063/1.3134422
[18]	Nguyen, A., González de Alaiza Martnez, 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
[19]	Ammosov, M.V., Delone, N.B. and Krainov, V.P. (1986) Tunnel Ionization of Complex Atoms and of Atomic Ions in an Alternating Electromagnetic Field. Soviet Physics JETP, 64, 1191-1194.
[20]	Jackson, J.D. (1925-2016) Classical Electrodynamics. 3rd Edition, Wiley, New York.
[21]	Hockney, R.W., Goel, S.P. and Eastwood, J.W. (1974) Quiet High-Resolution Computer Models of a Plasma. Journal of Computational Physics, 14, 148-158. https://doi.org/10.1016/0021-9991(74)90010-2
[22]	Chen, M., Cormier-Michel, E., Geddes, C.G.R., Bruhwiler, D.L., Yu, L.L., Esarey, E., Schroeder, C.B. and Leemans, W.P. (2013) Numerical Modeling of Laser Tunneling Ionization in Explicit Particle-in-Cell Codes. Journal of Computational Physics, 236, 220-228. https://doi.org/10.1016/j.jcp.2012.11.029
[23]	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
[24]	Xi, T.T., Zhao, Z.J. and Hao, Z.Q. (2014) Filamentation of Femtosecond Laser Pulses with Spatial Chirp in Air. Journal of the Optical Society of America B, 31, 321-324. https://doi.org/10.1364/JOSAB.31.000321
[25]	González De Alaiza Martnez, P., Babushkin, I., Bergé, L., Skupin, S., Cabrera-Granado, E., Köhler, C., Morgner, U., Husakou, A. and Herrmann, J. (2015) Boosting Terahertz Generation in Laser-Field Ionized Gases Using a Sawtooth Wave Shape. Physical Review Letters, 114, 183-901. https://doi.org/10.1103/PhysRevLett.114.183901
[26]	Zhao, H., Zhang, L.L., Huang, S.X., Zhang, S.J. and Zhang, C.L. (2018) Terahertz Wave Generation from Noble Gas Plasmas Induced by a Wavelength-Tunable Femtosecond Laser. IEEE Transactions on Terahertz Science and Technology, 8, 299-304. https://doi.org/10.1109/TTHZ.2018.2820425
[27]	Debayle, A., González De Alaiza Martnez, P., Gremillet, L. and Bergé, L. (2015) Nonmonotonic Increase in Laser-Driven Thz Emissions through Multiple Ionization Events. Physical Review A, 91, Article ID: 041801. https://doi.org/10.1103/PhysRevA.91.041801
[28]	Chen, Y.Q., Yamaguchi, M., Wang, M.F. and Zhang, X.-C. (2007) Terahertz Pulse Generation from Noble Gases. Applied Physics Letters, 91, Article ID: 251116. https://doi.org/10.1063/1.2826544
[29]	Moradi, S., Ganjovi, A., Shojaei, F. and Saeed, M. (2015) Influences of Different Gases on the Terahertz Radiation Based on the Application of Two-Color Laser Pulses. Physics of Plasmas, 22, Article ID: 103115. https://doi.org/10.1063/1.4931168
[30]	Wang, T.-J., Chen, Y.P., Marceau, C., Thberge, F., Chateauneuf, M., Dubois, J. and Chin, S.L. (2009) High Energy Terahertz Emission from Two-Color Laser-Induced Filamentation in Air with Pump Pulse Duration Control. Applied Physics Letters, 95, Article ID: 131108. https://doi.org/10.1063/1.3242024
[31]	Wang, W.-M., Sheng, Z.-M., Wu, H.-C., Chen, M., Li, C., Zhang, J. and Mima, K. (2008) Strong Terahertz Pulse Generation by Chirped Laser Pulses in Tenuous Gases. Optics Express, 16, 16999-17006. https://doi.org/10.1364/OE.16.016999