1 Introduction
Present day numerical methods (in particular, the finite-difference time-domain (FDTD) method (Yee Reference Yee1966; Taflove & Hagness Reference Taflove and Hagness2005)) allow one to simulate the propagation of an electromagnetic signal through plasmas in real time. While the propagation of the short-term pulses is always transient, for the long-term pulses there is a transition time to self-consistent oscillations of plasma particles in the external electromagnetic field. The electromagnetic wave propagation through plasmas with such self-consistent oscillations is well described by the plasma theory and is the basis of classical plasma electrodynamics. But the transient plasma electromagnetic field is still not fully described analytically, making it difficult to predict the initial stage of the external electromagnetic wave interaction with the plasma and, especially, to analyse the data obtained in the numerical simulations. In particular, the analytics can predict in advance if the simulated electromagnetic field is transient or not.
The transient plasma electromagnetic filed was described already both analytically and numerically (Pavlenko et al. Reference Pavlenko, Girka, Trush, Melnyk and Velizhanina2019). But the task was somewhat artificial. A metal antenna with alternating surface current was assumed to be placed inside the plasma. As a result, if the plasma permittivity goes to zero, the amplitude of the plasma electric field grows to infinity. That is why the duration of the transient processes was very large for plasmas with small permittivities. The correlation between the analytical and numerical data was demonstrated to be very good. The developed analytics are applied in the present paper in order to study the normal incidence of a single frequency plane wave on a semi-bounded uniform, cold, unmagnetized and collisionless plasma which is initially unperturbed.
The problem of the normal wave incidence on the plasma–vacuum interface is widely used to validate FDTD codes (Luebbers, Hunsberger & Kunz Reference Luebbers, Hunsberger and Kunz1991; Young Reference Young1996; Gamliel Reference Gamliel2017). Therein, the wave magnetic field is an undifferentiable function of the coordinate at the plasma–vacuum interface. This fact always introduces a numerical error in comparison with the analytical description because the FDTD codes represent the differential operators in a numerical way. This error was additionally emphasized in the Conclusion of Young (Reference Young1996). It is negligible if the FDTD codes simulate short-term pulses or are validated by short-term pulses. However, in general, special attention should be paid in the numerical simulations to the sharp interface between the vacuum and the plasma. This problem is discussed in the present paper in detail.
The rather extensive review of the relevant publications can be found in the Introduction of Pavlenko et al. (Reference Pavlenko, Girka, Trush, Melnyk and Velizhanina2019). However, the most important published developments are outlined below.
The proposed problem was studied many years ago as a part of the topic related to the electromagnetic wave incidence on plasmas. It was solved by the Laplace transform method (Case Reference Case1964) where the transmitted field was written as the complicated double integral of equation (61). A similar problem was considered by Schmitt (Reference Schmitt1964, Reference Schmitt1965), but the author did not study the transmitted signal in detail. The transient signal propagation in isotropic plasmas (Haskell & Case Reference Haskell and Case1967) was described by the saddle-point method and by the high-frequency expansion technique. And although a signal of the step sinusoidal wave was launched inside the plasma (Haskell & Case Reference Haskell and Case1967), the transient propagation had much in common mathematically with the problem of the incident wave.
Pulse propagation of an arbitrary time duration in dispersive media was studied by Oughstun (Reference Oughstun1991). The asymptotic description with the integral representation of the propagated plane-wave pulse was used in the study. The same technique was used to study the propagation of the step sinusoidal wave (Cartwright & Oughstun Reference Cartwright and Oughstun2009). Although the obtained result was valid for all input carrier frequencies, special attention was paid to the carrier frequency below the plasma frequency. A more detailed description of the method and its results can be found in Oughstun (Reference Oughstun2009).
The propagation of a transverse electromagnetic plane wave of the step sinusoidal signal variety was studied by Zablocky & Engheta (Reference Zablocky and Engheta1993). The transient behaviour of the electric field components at different locations and at various times of observation was analysed using the standard Fourier transform technique. In the paper the time-varying electric field was applied at some fixed coordinate in the medium.
The transient electric signal propagation in a dispersive medium (Lam Reference Lam1974) was simulated numerically as a direct solution of the integro-differential equation. In the case of a lossless medium the equation was reduced, of course, to a differential equation, but both coordinate and time discretization was still required.
Progress in the numerical description of transient signal propagation in plasmas was demonstrated as the frequency-dependent finite-difference time-domain formulation (Luebbers et al. Reference Luebbers, Hunsberger and Kunz1991). The method is based on the recursive convolution technique and was applied effectively to calculate the reflection and transmission coefficients for plane-wave incidence on a plasma slab. Later, the method was used successfully by many authors to study electromagnetic signal propagation in dispersive media.
The electromagnetic field reflected from the plasmas was studied in many other papers (see, for example, Wait (Reference Wait1969), Kalluri (Reference Kalluri1988) and references in the books Wait (Reference Wait1970), Kalluri (Reference Kalluri2018)), which is why it is not discussed in the present paper. Only the transmitted plasma electromagnetic field is described. A general present day view of transient fields in dispersive media can be found in the books Felsen (Reference Felsen1976), Orfanidis (Reference Orfanidis2016) and Kalluri (Reference Kalluri2018).
In the present paper, the transmitted plasma electromagnetic field is written for the cold, unmagnetized, collisionless and underdense plasmas (the incident wave frequency is larger than the plasma frequency) in a simple integral form obtained by the Fourier integral method (Budak, Samarskii & Tikhonov Reference Budak, Samarskii and Tikhonov1964). The case when the frequencies are equal is studied in detail. The analytics of the problem are presented in § 2. The comparison between the analytical and numerical results is provided in § 3. Practical recommendations for analytical and numerical studies of similar problems are given in the conclusions.
2 Analytical description of the transmitted plasma field
A uniform hydrogen plasma is assumed to occupy the semi-space 
$x>0$ (figure 1). The model of a sharp boundary of the homogeneous plasma is widely used to study the propagation characteristics of the electromagnetic waves both analytically and numerically (Thoma et al. Reference Thoma, Rose, Miller, Clark and Hughes2009; Yuan, Zhou & Sun Reference Yuan, Zhou and Sun2010; Zheng et al. Reference Zheng, Zhao, Liu, Xing and Chen2014). The electromagnetic wave from the vacuum semi-space 
$x<0$ is incident normally upon the plasma–vacuum interface. The 
$y$ axis is directed along the wave electric field and the 
$z$ axis is directed along the wave magnetic field. The wave is partially reflected from the interface and partially transmitted into the plasma. Both the magnetic and electric fields of the incident wave have the following time–coordinate dependence in the vacuum:
$$\begin{eqnarray}B_{z},E_{y}=\sin (\unicode[STIX]{x1D714}t-kx),\end{eqnarray}$$ where 
$\unicode[STIX]{x1D714}$ is the angular wave frequency, 
$k=\unicode[STIX]{x1D714}/c$ is the vacuum wavenumber and 
$c$ is the speed of light in vacuum. The wave amplitude is normalized to be unity and is small enough to neglect the contributions of the plasma density perturbations in the plasma current densities and the plasma frequencies (see § 3). The objective is to describe in detail the transmitted electromagnetic field.
Figure 1. Schematic of the problem.
Both the magnetic and electric fields are governed by the wave equation in the vacuum
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}v}{\unicode[STIX]{x2202}x^{2}}=\frac{1}{c^{2}}\frac{\unicode[STIX]{x2202}^{2}v}{\unicode[STIX]{x2202}t^{2}}\end{eqnarray}$$and by the telegraph equation in the plasma
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}u}{\unicode[STIX]{x2202}x^{2}}=\frac{1}{c^{2}}\frac{\unicode[STIX]{x2202}^{2}u}{\unicode[STIX]{x2202}t^{2}}+\frac{\unicode[STIX]{x1D714}_{pl}^{2}}{c^{2}}u,\end{eqnarray}$$ where 
$\unicode[STIX]{x1D714}_{pl}^{2}=\sum _{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D6FC}}^{2}$ and 
$\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D6FC}}$ is the plasma frequency of plasma species 
$\unicode[STIX]{x1D6FC}$.
Both the magnetic and electric fields should be continuous at the plasma–vacuum interface
$$\begin{eqnarray}\displaystyle & \displaystyle B_{z}^{vac}(0,t)=B_{z}^{pl}(0,t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle E_{y}^{vac}(0,t)=E_{y}^{pl}(0,t). & \displaystyle\end{eqnarray}$$The continuity of the coordinate derivative of the electric field at the plasma–vacuum interface follows from the continuity of the magnetic field (2.4)
$$\begin{eqnarray}\left.\frac{\unicode[STIX]{x2202}E_{y}^{vac}}{\unicode[STIX]{x2202}x}\right|_{x=0}=\left.\frac{\unicode[STIX]{x2202}E_{y}^{pl}}{\unicode[STIX]{x2202}x}\right|_{x=0}.\end{eqnarray}$$The solution of (2.2) for the electromagnetic field in the vacuum region should be in the form (Jackson Reference Jackson1998)
$$\begin{eqnarray}v(x,t)=\sin (\unicode[STIX]{x1D714}t-kx)+f(\unicode[STIX]{x1D714}t+kx),\end{eqnarray}$$ where 
$f(\unicode[STIX]{x1D714}t+kx)$ is any twice differentiable functional of time-coordinate dependence 
$\unicode[STIX]{x1D714}t+kx$ which describes the reflected wave. Then the boundary conditions for the electric field (
$u(x,t)\equiv E_{y}^{pl}(x,t)$) at the plasma–vacuum interface become
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\sin \unicode[STIX]{x1D714}t+f(\unicode[STIX]{x1D714}t)=u(0,t),\\ -k\cos \unicode[STIX]{x1D714}t+kf^{\prime }(\unicode[STIX]{x1D714}t)=u_{x}(0,t),\end{array}\right\}\end{eqnarray}$$ where 
$u_{x}$ means the partial derivative with respect to the coordinate and 
$f^{\prime }(g)$ means the derivative of 
$f(g)$ with respect to the argument.
Excluding the unknown functional 
$f(g)$ from (2.8) one can get the boundary condition for the electric field in the plasma
$$\begin{eqnarray}\left\{\frac{1}{c}u_{t}(x,t)-u_{x}(x,t)\right\}_{x=0}=2k\cos \unicode[STIX]{x1D714}t,\quad t\geqslant 0.\end{eqnarray}$$Thus, the problem which has to be solved for the electric field in the plasma consists of (2.3) with the boundary condition (2.9) and two initial conditions
$$\begin{eqnarray}\displaystyle & \displaystyle u(x,0)=0,\quad x\geqslant 0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u_{t}(x,0)=0,\quad x\geqslant 0. & \displaystyle\end{eqnarray}$$ Let the function 
$\unicode[STIX]{x1D719}(x,t)$ be any solution of (2.3) which can be written (see chapter II in Budak et al. (Reference Budak, Samarskii and Tikhonov1964) and Pavlenko et al. (Reference Pavlenko, Girka, Trush, Melnyk and Velizhanina2019)) in the form
$$\begin{eqnarray}\unicode[STIX]{x1D719}(x,t)=\int _{x}^{ct}\unicode[STIX]{x1D713}(g-x)J_{0}\left(\frac{\unicode[STIX]{x1D714}_{pl}}{c}\sqrt{c^{2}t^{2}-g^{2}}\right)\text{d}g,\end{eqnarray}$$ where 
$J_{0}(x)$ is the Bessel function of the first kind. The solution of (2.3) can be written as (2.12) because 
$J_{0}(\unicode[STIX]{x1D714}_{pl}\sqrt{(ct-\unicode[STIX]{x1D701})^{2}-(x-\unicode[STIX]{x1D709})^{2}}/c)$ is Riemann’s function of the telegraph equation (Tikhonov & Samarskii Reference Tikhonov and Samarskii2011, chapter II, § 5).
However, the solution of the problem (2.3), (2.9)–(2.11) is sought in the form
$$\begin{eqnarray}u(x,t)=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}+c\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}x},\end{eqnarray}$$ which is obviously also a solution of (2.3) if the required derivatives of the function 
$\unicode[STIX]{x1D719}(x)$ can be calculated.
Then the boundary condition for the function 
$\unicode[STIX]{x1D719}(x,t)$ (which is not the electric field) becomes
$$\begin{eqnarray}\unicode[STIX]{x1D719}(0,t)=-2\frac{\unicode[STIX]{x1D714}}{\unicode[STIX]{x1D714}_{pl}^{2}}\cos \unicode[STIX]{x1D714}t.\end{eqnarray}$$Luckily, the problem defined by (2.3), (2.10), (2.11) with the boundary condition (2.14) has been solved already (Pavlenko et al. Reference Pavlenko, Girka, Trush, Melnyk and Velizhanina2019), therefore
$$\begin{eqnarray}\unicode[STIX]{x1D719}(x,t)=-2\frac{k}{\unicode[STIX]{x1D714}_{pl}^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\int _{x}^{ct}\cos (k\sqrt{\unicode[STIX]{x1D716}}(g-x))J_{0}\left(\frac{\unicode[STIX]{x1D714}_{pl}}{c}\sqrt{c^{2}t^{2}-g^{2}}\right)\text{d}g,\end{eqnarray}$$ where 
$\unicode[STIX]{x1D716}=1-\unicode[STIX]{x1D714}_{pl}^{2}/\unicode[STIX]{x1D714}^{2}$ is the plasma permittivity.
Using the notation
$$\begin{eqnarray}I_{1}(x,t)\equiv \int _{x}^{ct}\cos (k\sqrt{\unicode[STIX]{x1D716}}(g-x))J_{0}\left(\frac{\unicode[STIX]{x1D714}_{pl}}{c}\sqrt{c^{2}t^{2}-g^{2}}\right)\text{d}g,\end{eqnarray}$$the solution of the problem for the electric field due to (2.13) becomes
$$\begin{eqnarray}u(x,t)=-2\frac{k}{\unicode[STIX]{x1D714}_{pl}^{2}}\left\{\frac{\unicode[STIX]{x2202}^{2}I_{1}}{\unicode[STIX]{x2202}t^{2}}+c\frac{\unicode[STIX]{x2202}^{2}I_{1}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}t}\right\}.\end{eqnarray}$$It is convenient to rewrite the second derivative with respect to time in (2.17) using (2.3)
$$\begin{eqnarray}u(x,t)=-2\frac{k}{\unicode[STIX]{x1D714}_{pl}^{2}}\left\{c^{2}\frac{\unicode[STIX]{x2202}^{2}I_{1}}{\unicode[STIX]{x2202}x^{2}}+c\frac{\unicode[STIX]{x2202}^{2}I_{1}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D714}_{pl}^{2}I_{1}\right\}.\end{eqnarray}$$ Calculating the partial derivatives with respect to the coordinate 
$x$ in (2.18) and introducing a new convenient notation
$$\begin{eqnarray}I_{2}(x,t)\equiv \int _{x}^{ct}\sin (k\sqrt{\unicode[STIX]{x1D716}}(g-x))J_{0}\left(\frac{\unicode[STIX]{x1D714}_{pl}}{c}\sqrt{c^{2}t^{2}-g^{2}}\right)\text{d}g,\end{eqnarray}$$the exact solution for the electric field becomes
$$\begin{eqnarray}\displaystyle E_{y}(x,t) & = & \displaystyle -\frac{2}{(1-\unicode[STIX]{x1D716})^{1/2}}\frac{ct-x}{\sqrt{c^{2}t^{2}-x^{2}}}J_{1}\left(\frac{\unicode[STIX]{x1D714}_{pl}}{c}\sqrt{c^{2}t^{2}-x^{2}}\right)\nonumber\\ \displaystyle & & \displaystyle +\,\frac{2}{1-\unicode[STIX]{x1D716}}\left\{kI_{1}-\frac{\sqrt{\unicode[STIX]{x1D716}}}{c}\frac{\unicode[STIX]{x2202}I_{2}}{\unicode[STIX]{x2202}t}\right\}.\end{eqnarray}$$The magnetic field is obtained easily from (2.17) and (2.20) using Maxwell equations
$$\begin{eqnarray}\displaystyle B_{z}(x,t) & = & \displaystyle \frac{2}{(1-\unicode[STIX]{x1D716})^{1/2}}\frac{ct-x}{\sqrt{c^{2}t^{2}-x^{2}}}J_{1}\left(\frac{\unicode[STIX]{x1D714}_{pl}}{c}\sqrt{c^{2}t^{2}-x^{2}}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\frac{2}{1-\unicode[STIX]{x1D716}}\left\{\unicode[STIX]{x1D716}kI_{1}-\frac{\sqrt{\unicode[STIX]{x1D716}}}{c}\frac{\unicode[STIX]{x2202}I_{2}}{\unicode[STIX]{x2202}t}\right\}.\end{eqnarray}$$ The expressions for the electric and magnetic fields differ in sign and multiplier 
$\unicode[STIX]{x1D716}$ in the first term of the second row.
For the obtained electromagnetic field in the plasmas (2.20), (2.21) one can repeat the asymptotic analysis of the integrals 
$I_{1}$ and 
$I_{2}$ in Pavlenko et al. (Reference Pavlenko, Girka, Trush, Melnyk and Velizhanina2019). In the approximations 
$\unicode[STIX]{x1D714}_{pl}t\gg 1$ and 
$x^{2}\ll (ct)^{2}$ (for plasmas with 
$\unicode[STIX]{x1D716}\geqslant 0.5$) or 
$(1-\unicode[STIX]{x1D716})x^{2}\ll \unicode[STIX]{x1D716}(ct)^{2}$ (for plasmas with 
$\unicode[STIX]{x1D716}\leqslant 0.5$) the plasma electric and magnetic fields become
$$\begin{eqnarray}\displaystyle E_{y}(x,t) & = & \displaystyle \frac{2}{1+\sqrt{\unicode[STIX]{x1D716}}}\sin (\unicode[STIX]{x1D714}t-k\sqrt{\unicode[STIX]{x1D716}}x)\nonumber\\ \displaystyle & & \displaystyle -\,2\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{(1-\unicode[STIX]{x1D716})^{1/4}}{\unicode[STIX]{x1D716}(\unicode[STIX]{x1D714}t)^{1/2}}\left(\frac{kx}{\unicode[STIX]{x1D714}t}-\frac{(kx)^{2}}{(\unicode[STIX]{x1D714}t)^{2}}\right)\sin (\unicode[STIX]{x1D6E5}-\unicode[STIX]{x03C0}/4),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle B_{z}(x,t) & = & \displaystyle \frac{2\sqrt{\unicode[STIX]{x1D716}}}{1+\sqrt{\unicode[STIX]{x1D716}}}\sin (\unicode[STIX]{x1D714}t-k\sqrt{\unicode[STIX]{x1D716}}x)\nonumber\\ \displaystyle & & \displaystyle -\,2\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{(1-\unicode[STIX]{x1D716})^{1/4}}{\unicode[STIX]{x1D716}(\unicode[STIX]{x1D714}t)^{1/2}}\left(\frac{(kx)^{2}}{(\unicode[STIX]{x1D714}t)^{2}}-\frac{(kx)^{3}}{(\unicode[STIX]{x1D714}t)^{3}}\right)\sin (\unicode[STIX]{x1D6E5}-\unicode[STIX]{x03C0}/4),\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6E5}\equiv \unicode[STIX]{x1D714}_{pl}\sqrt{c^{2}t^{2}-x^{2}}/c$.
The electric field becomes steady state when the amplitude of the second term in (2.22) is much lower than the amplitude of the first term. Similarly, the magnetic field becomes steady state when the amplitude of the second term in (2.23) is much lower than the amplitude of the first term. But, in general, one can say that the steady state is reached when the amplitude deviation of the Poynting flux from the steady state one becomes less than some small value 
$\unicode[STIX]{x1D702}$. Here, the Poynting flux means the component of the Poynting vector 
$S_{x}=cE_{y}B_{z}/4\unicode[STIX]{x03C0}$, since the problem is uniform in cross-section to the direction 
$x$.
From (2.22), (2.23), the dominant terms of the normalized Poynting flux are
$$\begin{eqnarray}\displaystyle S_{x}(x,t) & = & \displaystyle \frac{4\sqrt{\unicode[STIX]{x1D716}}}{(1+\sqrt{\unicode[STIX]{x1D716}})^{2}}\sin ^{2}(\unicode[STIX]{x1D714}t-k\sqrt{\unicode[STIX]{x1D716}}x)\nonumber\\ \displaystyle & & \displaystyle -\,4\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{(1-\unicode[STIX]{x1D716})^{1/4}}{\unicode[STIX]{x1D716}(1+\sqrt{\unicode[STIX]{x1D716}})}\frac{1}{(\unicode[STIX]{x1D714}t)^{1/2}}\nonumber\\ \displaystyle & & \displaystyle \times \left(\sqrt{\unicode[STIX]{x1D716}}\frac{kx}{\unicode[STIX]{x1D714}t}+(1-\sqrt{\unicode[STIX]{x1D716}})\frac{(kx)^{2}}{(\unicode[STIX]{x1D714}t)^{2}}\right)\sin (\unicode[STIX]{x1D714}t-k\sqrt{\unicode[STIX]{x1D716}}x)\sin (\unicode[STIX]{x1D6E5}-\unicode[STIX]{x03C0}/4).\end{eqnarray}$$Therefore, the Poynting flux becomes steady state when the amplitude of the second term in (2.24) is much lower than the amplitude of the first term
$$\begin{eqnarray}\frac{x}{\unicode[STIX]{x1D706}}<\frac{1}{2}\frac{t}{T}\frac{\sqrt{\unicode[STIX]{x1D716}}}{1-\sqrt{\unicode[STIX]{x1D716}}}\left\{\sqrt{1+4\unicode[STIX]{x1D702}\unicode[STIX]{x03C0}\frac{\sqrt{\unicode[STIX]{x1D716}}(1-\sqrt{\unicode[STIX]{x1D716}})^{3/4}}{(1+\sqrt{\unicode[STIX]{x1D716}})^{5/4}}\left(\frac{t}{T}\right)^{1/2}}-1\right\},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D706}$ is the vacuum wavelength and 
$T$ is the wave period. The duration of the transient processes of the magnetic and electric field formation can be obtained easily from (2.22), (2.23) in the manner presented in Pavlenko et al. (Reference Pavlenko, Girka, Trush, Melnyk and Velizhanina2019), therefore these relations are not written here.
The formulas (2.22), (2.23) are written already without the small terms which are proportional to the small parameter 
$1/\unicode[STIX]{x1D714}_{pl}t$ and do not have a coordinate dependence. These terms can be neglected usually inside the plasma where 
$\unicode[STIX]{x1D714}_{pl}t\gg (ct)^{2}/x^{2}$ but they dominate near the plasma–vacuum interface. Therefore exactly at the plasma–vacuum interface in the aforementioned approximations the electromagnetic field becomes
$$\begin{eqnarray}\displaystyle & \displaystyle E_{y}(0,t)=\frac{2}{1+\sqrt{\unicode[STIX]{x1D716}}}\sin \unicode[STIX]{x1D714}t+2\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{1}{\unicode[STIX]{x1D716}(1-\unicode[STIX]{x1D716})^{1/4}}\frac{1}{(\unicode[STIX]{x1D714}t)^{3/2}}\cos (\unicode[STIX]{x1D714}_{pl}t-\unicode[STIX]{x03C0}/4), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle B_{z}(0,t)=\frac{2\sqrt{\unicode[STIX]{x1D716}}}{1+\sqrt{\unicode[STIX]{x1D716}}}\sin \unicode[STIX]{x1D714}t-2\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{1}{\unicode[STIX]{x1D716}(1-\unicode[STIX]{x1D716})^{1/4}}\frac{1}{(\unicode[STIX]{x1D714}t)^{3/2}}\cos (\unicode[STIX]{x1D714}_{pl}t-\unicode[STIX]{x03C0}/4). & \displaystyle\end{eqnarray}$$The formulas (2.26), (2.27) can be helpful for an analysis of the reflected field.
The time of the transient processes for the magnetic field (dotted line), electric field (dashed line) and Poynting flux (thick solid line) are presented in figure 2 as a function of the distance from the plasma–vacuum interface. The curves have been built for a plasma with 
$\unicode[STIX]{x1D716}=0.16$ and show when the ratio of the second term amplitude in (2.23), (2.22), (2.24) to the steady state amplitude becomes 
$\unicode[STIX]{x1D702}=0.02$ for the magnetic field, electric field and Poynting flux respectively. The values 
$\unicode[STIX]{x1D716}$ and 
$\unicode[STIX]{x1D702}$ are chosen because they were used when studying the case of the metal antenna inside the plasma (Pavlenko et al. Reference Pavlenko, Girka, Trush, Melnyk and Velizhanina2019). This gives the opportunity to compare the time of the transient processes in the following two cases: the electromagnetic wave is incident from the vacuum (the present approach) and the electromagnetic wave is launched by the antenna inside the plasma (presented in Pavlenko et al. (Reference Pavlenko, Girka, Trush, Melnyk and Velizhanina2019)). The time of the transient processes is smaller if the electromagnetic wave is incident from the vacuum. This is discussed in the next section.
Figure 2. Time of transient processes for plasma with 
$\unicode[STIX]{x1D716}=0.16$ as a function of the distance from the plasma–vacuum interface. Above each curve, the amplitude deviation of the corresponding value from the steady state one is less than 2 %. The approximations used are not valid below the thin solid line.
2.1 Formation of electrostatic oscillations in plasmas
There is a practical interest in the possibility of both analytical and numerical descriptions of the formation of the electrostatic oscillation in plasmas when the frequency of the incident wave is equal to the plasma frequency (or, in a more general sense, when the incident wave packet includes the frequency which is equal to the plasma frequency). Of course, thermal effects and collisions are very important for the processes near 
$\unicode[STIX]{x1D714}_{pl}$. However, sometimes they can be left in the background (Kylychbekov et al. Reference Kylychbekov, Song, Kwon, Ra, Yoon, Chung, Yu, Yoffe, Ersfeld and Jaroszynski2020). Here, the transient propagation at 
$\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{pl}$ is studied as a limit case of the developed cold collisionless model. The thermal effects and the collisions can be included in a separate study to bring the model closer to the experimental behaviour.
The structure of the electromagnetic field in a plasma in this case (initially the wave field is still electromagnetic and will become electrostatic in the large time limit only) can be obtained from (2.20), (2.21) if one puts 
$\unicode[STIX]{x1D716}=0$
$$\begin{eqnarray}\displaystyle & \displaystyle B_{z}(x,t)=2\frac{ct-x}{\sqrt{c^{2}t^{2}-x^{2}}}J_{1}(k\sqrt{c^{2}t^{2}-x^{2}}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle E_{y}(x,t)=-2\frac{ct-x}{\sqrt{c^{2}t^{2}-x^{2}}}J_{1}(k\sqrt{c^{2}t^{2}-x^{2}})+2k\int _{x}^{ct}J_{0}(k\sqrt{c^{2}t^{2}-g^{2}})\,\text{d}g. & \displaystyle\end{eqnarray}$$The first term in both equations describes the electromagnetic field which decays with time. The term is referred to further as the ‘electromagnetic term’. It describes a vanishing magnetic field. The second term of the electric field dependence (2.29) can be rewritten with the help of the table integral 6.517 in Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2007)
$$\begin{eqnarray}\displaystyle E_{y}(x,t) & = & \displaystyle -2\frac{ct-x}{\sqrt{c^{2}t^{2}-x^{2}}}J_{1}(k\sqrt{c^{2}t^{2}-x^{2}})\nonumber\\ \displaystyle & & \displaystyle +\,2\sin \unicode[STIX]{x1D714}t-2k\int _{0}^{x}J_{0}(k\sqrt{c^{2}t^{2}-g^{2}})\,\text{d}g.\end{eqnarray}$$ Here the second term describes a Langmuir oscillation (or plasma oscillation) (Krall & Trivelpiece Reference Krall and Trivelpiece1973, § 4.2), (Chen Reference Chen2012, § 4.3) with zero wavenumber in the limit 
$t\rightarrow \infty$ and the third term corrects it at intermediate time. Usually a Langmuir oscillation with zero wavenumber is called a plasma dipole oscillation (Kylychbekov et al. Reference Kylychbekov, Song, Kwon, Ra, Yoon, Chung, Yu, Yoffe, Ersfeld and Jaroszynski2020) because there is no fluctuation of the plasma density in this case. But it is an electrostatic oscillation since there is no fluctuating magnetic field. The amplitude of the electrostatic oscillation in a steady state is twice the amplitude of the incident electric field. The third term is referred to below as the ‘electrostatic term’.
Initially the plasma is electrodynamically similar to a vacuum, therefore, at 
$t=0$ all incident power penetrates into the plasma. Build-up of the plasma particle oscillations near the plasma–vacuum interface reduces the plasma’s ability to transmit the incident electromagnetic power. Therefore, over time, an increasing fraction of the incident power is reflected back into the vacuum. Finally, in the limit 
$t\rightarrow \infty$, the incident power is fully reflected. One can calculate the time when the fraction of the power penetrating into the plasma becomes less than 
$\unicode[STIX]{x1D702}$.
The electromagnetic field at the plasma–vacuum interface is
$$\begin{eqnarray}\displaystyle & \displaystyle B_{z}(0,t)=2J_{1}(\unicode[STIX]{x1D714}t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle E_{y}(0,t)=-2J_{1}(\unicode[STIX]{x1D714}t)+2\sin \unicode[STIX]{x1D714}t. & \displaystyle\end{eqnarray}$$Therefore, the normalized Poynting flux from the vacuum to the plasma is
$$\begin{eqnarray}S_{x}(0,t)=4J_{1}(\unicode[STIX]{x1D714}t)\sin \unicode[STIX]{x1D714}t-4J_{1}^{2}(\unicode[STIX]{x1D714}t).\end{eqnarray}$$ The Bessel function 
$J_{1}(\unicode[STIX]{x1D714}t)$ can be expanded in the large time limit 
$\unicode[STIX]{x1D714}t\gg 1$. Therefore, the Poynting flux amplitude is less than 
$\unicode[STIX]{x1D702}$ if the following inequality is valid:
$$\begin{eqnarray}4\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{1}{(\unicode[STIX]{x1D714}t)^{1/2}}<\unicode[STIX]{x1D702}.\end{eqnarray}$$ However, the power transfer through the plasma–vacuum interface is described by the Poynting flux averaged over the wave period. Therefore, the ratio of the penetrated power to the incident power becomes less than 
$\unicode[STIX]{x1D702}$ when
$$\begin{eqnarray}\frac{t}{T}>\frac{8}{\unicode[STIX]{x03C0}^{2}}\frac{1}{\unicode[STIX]{x1D702}^{2}}.\end{eqnarray}$$For example, the power flux in plasma becomes less than 5 %, 2 % and 1 % of the incident power flux after approximately 325, 2025 and 8100 wave periods respectively.
Now one can easily separate the ranges of the electromagnetic and quasi-electrostatic fields on the coordinate–time plane. The electromagnetic term (and the magnetic field itself) can be expanded in the large time limit 
$\unicode[STIX]{x1D714}t\gg 1$ as
$$\begin{eqnarray}B_{z}(x,t)=2\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{(\unicode[STIX]{x1D714}t-kx)^{1/4}}{(\unicode[STIX]{x1D714}t+kx)^{3/4}}\sin (k\sqrt{c^{2}t^{2}-x^{2}}-\unicode[STIX]{x03C0}/4).\end{eqnarray}$$ It is convenient to introduce the useful notation 
$\unicode[STIX]{x1D707}\equiv kx/\unicode[STIX]{x1D714}t$ which is not a small parameter in this subsection. The wave field can be considered as a quasi-electrostatic one if the magnetic field amplitude does not exceed a small ratio 
$\unicode[STIX]{x1D702}$ of the incident magnetic field amplitude
$$\begin{eqnarray}2\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{(1-\unicode[STIX]{x1D707})^{1/4}}{(1+\unicode[STIX]{x1D707})^{3/4}}\frac{1}{(\unicode[STIX]{x1D714}t)^{1/2}}<\unicode[STIX]{x1D702}.\end{eqnarray}$$Thus, the wave field is quasi-electrostatic in the range
$$\begin{eqnarray}\frac{t}{T}>\frac{4}{\unicode[STIX]{x03C0}^{2}}\frac{(1-\unicode[STIX]{x1D707})^{1/2}}{(1+\unicode[STIX]{x1D707})^{3/2}}\frac{1}{\unicode[STIX]{x1D702}^{2}}.\end{eqnarray}$$ The last relation is not resolved for both coordinate and time. However, it is convenient for graphical analysis. The analysis is presented in figure 3. The amplitude of the magnetic field in plasma is shown by the contour lines on the coordinate–time plane. The lines show where the magnetic field amplitude is 5 %, 4 %, 3 % and 2 % of the amplitude of the incident magnetic field. One can see that the magnetic field amplitude never reaches 5 %, 4 %, 3 % and 2 % at the distances from the plasma–vacuum interface 
$d>32\,\unicode[STIX]{x1D706}$, 
$d>49\,\unicode[STIX]{x1D706}$, 
$d>87\,\unicode[STIX]{x1D706}$ and 
$d>195\,\unicode[STIX]{x1D706}$ respectively. Similar data can be obtained from (2.38) for any value of the magnetic field amplitude. Also, the magnetic field amplitude in plasma near the plasma–vacuum interface drops to 5 %, 4 %, 3 % and 2 % at times 
$162\,T$, 
$253\,T$, 
$450\,T$ and 
$1013\,T$ respectively.
Figure 3. Magnetic field amplitude in plasma with 
$\unicode[STIX]{x1D716}=0$ normalized by the magnetic field amplitude of the incident wave as the contour lines on the coordinate–time plane (see text for details). The solid thin line is 
$x=ct$.
And, finally, let us look at the electrostatic term which defines the electric field amplitude of the forming electrostatic oscillation as a function of coordinate and time. Lommel’s expansion (Watson Reference Watson1944) allows one to integrate the term
$$\begin{eqnarray}\int _{0}^{x}J_{0}(k\sqrt{c^{2}t^{2}-g^{2}})\,\text{d}g=\mathop{\sum }_{m=0}^{\infty }\frac{1}{m!(2m+1)}\left(\frac{\unicode[STIX]{x1D714}x^{2}}{2c^{2}t}\right)^{m}xJ_{m}(\unicode[STIX]{x1D714}t).\end{eqnarray}$$ In the large time limit 
$\unicode[STIX]{x1D714}t\gg 1$, the sum on the right-hand side, after multiplying by 
$k$, can be rewritten through the Fresnel integrals 
${\mathcal{C}}(g)$ and 
${\mathcal{S}}(g)$ defined as 8.250 in Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2007), therefore
$$\begin{eqnarray}\displaystyle E_{y}(x,t) & = & \displaystyle -2\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{(\unicode[STIX]{x1D714}t-kx)^{1/4}}{(\unicode[STIX]{x1D714}t+kx)^{3/4}}\sin (k\sqrt{c^{2}t^{2}-x^{2}}-\unicode[STIX]{x03C0}/4)\nonumber\\ \displaystyle & & \displaystyle +\,2\sin \unicode[STIX]{x1D714}t-2\sqrt{2}\cos (\unicode[STIX]{x1D714}t-\unicode[STIX]{x03C0}/4){\mathcal{C}}\left(\frac{kx}{\sqrt{2\unicode[STIX]{x1D714}t}}\right)\nonumber\\ \displaystyle & & \displaystyle -\,2\sqrt{2}\sin (\unicode[STIX]{x1D714}t-\unicode[STIX]{x03C0}/4){\mathcal{S}}\left(\frac{kx}{\sqrt{2\unicode[STIX]{x1D714}t}}\right).\end{eqnarray}$$At small distances from the plasma–vacuum interface
$$\begin{eqnarray}\frac{kx}{\sqrt{2\unicode[STIX]{x1D714}t}}\ll 1\end{eqnarray}$$the electrostatic term becomes
$$\begin{eqnarray}-\frac{2\sqrt{2}kx}{\sqrt{\unicode[STIX]{x03C0}\unicode[STIX]{x1D714}t}}\cos (\unicode[STIX]{x1D714}t-\unicode[STIX]{x03C0}/4).\end{eqnarray}$$ It’s amplitude can be compared with the amplitude of the steady state electrostatic oscillation (which is equal to 2) or with the amplitude of the incident electric field (which is equal to 1). An amplitude deviation of the electric field 
$\unicode[STIX]{x1D6FE}$ can be introduced as
$$\begin{eqnarray}\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{kx}{\sqrt{\unicode[STIX]{x1D714}t}}<\unicode[STIX]{x1D6FE}.\end{eqnarray}$$ Here, if 
$\unicode[STIX]{x1D6FE}$ is equal to 0.1, 0.3 or 0.5, then the relation (2.43) gives the coordinate–time range where the amplitude of the electric field oscillation exceeds 90 %, 70 % or 50 % of the steady state electric field amplitude respectively. It should be noted once again that 
$\unicode[STIX]{x1D6FE}=0.5$ gives the range where the amplitude of the electric field oscillation in the plasma exceeds the amplitude of the incident electric field. The results of the electric field analysis are demonstrated in figure 4. The amplitude of the electric field oscillation is shown as the contour lines on the coordinate–time plane. At the plasma–vacuum interface, the amplitude of the electric field oscillation is almost twice the amplitude of the incident electric field. At the time 
$t=325\,T$, for example, when the power flux into plasma becomes 5 % of the incident power flux, the amplitude of the electric field oscillation is 1.8, 1.4 and 1.0 of the amplitude of the incident electric field at the distances of 0.9, 2.7 and 4.5 wavelengths from the plasma–vacuum interface respectively. For all these distances the magnetic field amplitude is approximately 3.5 % of the incident wave amplitude. That is why the field oscillation is almost electrostatic (quasi-electrostatic). If the approximation (2.41) is violated (below the thin solid line of the graph), the amplitude has to be calculated from the formula (2.40).
Figure 4. Electric field amplitude in plasma with 
$\unicode[STIX]{x1D716}=0$ normalized by the electric field amplitude of the incident wave as the contour lines on the coordinate–time plane (see text for details). The approximation used is not valid below the thin solid line.
3 Numerical simulations of the transmitted plasma field
The obtained analytical results are compared with the results of numerical simulations by the FDTD code. The FDTD code solves Maxwell’s equations and the momentum equation for each plasma species 
$\unicode[STIX]{x1D6FC}$
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D735}\times \boldsymbol{E}=-\frac{1}{c}\frac{\unicode[STIX]{x2202}\boldsymbol{B}}{\unicode[STIX]{x2202}t},\\ \displaystyle \unicode[STIX]{x1D735}\times \boldsymbol{B}=\frac{1}{c}\frac{\unicode[STIX]{x2202}\boldsymbol{E}}{\unicode[STIX]{x2202}t}+\frac{4\unicode[STIX]{x03C0}}{c}\mathop{\sum }_{\unicode[STIX]{x1D6FC}}\boldsymbol{J}_{\unicode[STIX]{x1D6FC}},\\ \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\boldsymbol{J}_{\unicode[STIX]{x1D6FC}}=\frac{\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D6FC}}^{2}}{4\unicode[STIX]{x03C0}}\boldsymbol{E}+[\boldsymbol{J}_{\unicode[STIX]{x1D6FC}}\times \unicode[STIX]{x1D734}_{\unicode[STIX]{x1D6FC}}].\end{array}\right\}\end{eqnarray}$$ where 
$\boldsymbol{E}$ is the electric field in the medium, 
$\boldsymbol{B}$ is the magnetic field in the medium, 
$\boldsymbol{J}_{\unicode[STIX]{x1D6FC}}$ is the current density of plasma species 
$\unicode[STIX]{x1D6FC}$ based on the unperturbed plasma density, 
$\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D6FC}}^{2}\equiv 4\unicode[STIX]{x03C0}n_{\unicode[STIX]{x1D6FC}}q_{\unicode[STIX]{x1D6FC}}^{2}/m_{\unicode[STIX]{x1D6FC}}$ is the square of the plasma frequency of plasma species 
$\unicode[STIX]{x1D6FC}$ based on the unperturbed plasma density 
$n_{\unicode[STIX]{x1D6FC}}$, 
$\unicode[STIX]{x1D734}_{\unicode[STIX]{x1D6FC}}\equiv q_{\unicode[STIX]{x1D6FC}}\boldsymbol{B}_{0}/cm_{\unicode[STIX]{x1D6FC}}$ is the cyclotron frequency of plasma species 
$\unicode[STIX]{x1D6FC}$, 
$\boldsymbol{B}_{0}$ is the external magnetic field, 
$q_{\unicode[STIX]{x1D6FC}}$ is the charge of the plasma species 
$\unicode[STIX]{x1D6FC}$, 
$m_{\unicode[STIX]{x1D6FC}}$ is the mass of the plasma species 
$\unicode[STIX]{x1D6FC}$. Other features of the code have been described in Pavlenko et al. (Reference Pavlenko, Girka, Trush, Melnyk and Velizhanina2019). Therefore these details are not repeated here.
For the purpose of the present analysis the code uses the model of non-magnetized plasmas which consist of electrons and hydrogen ions. The pointwise hard source (Taflove & Hagness Reference Taflove and Hagness2005) is realized numerically for the electric field with the time dependence from (2.1). It is far away (in the vacuum) from the plasma–vacuum interface, therefore the wave reflected by the plasma does not reach it. The second boundary of the computational domain is also far away from the plasma–vacuum interface, so the transmitted plasma wave is not reflected from it. The incident electromagnetic wave in the vacuum at 
$t=0$ is prescribed in the code as an initial condition. Therefore, numerical modelling starts from the moment when the incident electromagnetic wave just reaches the plasma–vacuum interface.
3.1 Case of electromagnetic steady state field (
$\unicode[STIX]{x1D716}>0$)
The wave frequency (1 GHz) and the plasma density (
$1.0414\times 10^{10}~\text{cm}^{-3}$) are chosen as in similar calculations in Pavlenko et al. (Reference Pavlenko, Girka, Trush, Melnyk and Velizhanina2019) to enable a comparison of the transient propagations. The analytical formulas (2.22), (2.23), (2.25) are used to calculate the times when the amplitude deviations at the distance of 5 plasma wavelengths (12.5 vacuum wavelengths) from the plasma–vacuum interface are 5 %. For the plasma with 
$\unicode[STIX]{x1D716}=0.16$ these times are 48, 67 and 87 wave periods for the magnetic field, electric field and Poynting flux respectively. The calculated electromagnetic field for 
$t=48\,T$ is shown in figure 5. The shaded areas show the ranges of 5 % deviation from the steady state amplitudes. The vertical solid line at five plasma wavelengths divides the coordinates into the left side, where the amplitude of the magnetic field is inside the range (the magnetic field is close to the steady state), and the right side, where its amplitude goes out of the range. However, the electric field amplitude leaves the range on both the left and right sides (the electric field is far enough from the steady state on both sides).
In fact, both the magnetic field data and the electric field data in figure 5 consist of three different curves. Two of them were calculated by the FDTD code using the different spatial resolutions of the numerical grid (50 and 100 points per vacuum wavelength). The third curve was calculated directly from the exact formulas (2.20) (for the electric field) or (2.21) (for the magnetic field). These curves are indistinguishable from each other in figure 5. Therefore, the coordinates near the plasma–vacuum interface of figure 5 are zoomed in figure 6. It is seen in figure 6 that the data of the FDTD code with the spatial resolution of 50 points per vacuum wavelength differ essentially from other numerical data. But, at the same time, the curve from the FDTD code with the spatial resolution of 100 points per vacuum wavelength and the curve from the exact formulas are indistinguishable from each other. The spatial resolution of 50 points per vacuum wavelength is not enough to reproduce the exact formula data because the sharp plasma–vacuum interface introduces a numerical error which is discussed in the next subsection. The time resolution of the numerical grid (number of the points per wave period) is twice as large as the spatial resolution in all calculations in order to ensure numerical stability (Cummer Reference Cummer1997).
Figure 5. Electromagnetic field as a function of the distance from the plasma–vacuum interface for 
$t=48\,T$ and for a plasma with 
$\unicode[STIX]{x1D716}=0.16$. The magnetic field data are marked by the letter ‘B’. The electric field data are marked by the letter ‘E’. The shaded areas show the ranges of 5 % amplitude deviation from the steady state field.
Figure 6. Electromagnetic field as a function of the distance from the plasma–vacuum interface for 
$t=48\,T$ and for a plasma with 
$\unicode[STIX]{x1D716}=0.16$. The magnetic field data are marked by the letter ‘B’. The electric field data are marked by the letter ‘E’. The thin horizontal lines indicate the amplitudes of the steady state field. The graph consists of the curves obtained from the different numerical sources (see text for details).
One can see that the duration of the transient processes is shorter in the case of the wave incidence from the vacuum in comparison with the case of the metal antenna inside the plasma (Pavlenko et al. Reference Pavlenko, Girka, Trush, Melnyk and Velizhanina2019). This is easy to explain. According to the wave theory of steady state plasmas, the amplitude of the magnetic field in the plasma should be less (by a factor 
$\sqrt{\unicode[STIX]{x1D716}}$) than the amplitude of the electric field. For example, in the simulated case the magnetic field amplitude is 2.5 times smaller, since 
$\sqrt{\unicode[STIX]{x1D716}}=0.4$. Since the metal antenna inside the plasma sets the amplitude of the wave magnetic field, the electric field amplitude has to increase with time (initially they are equal to each other). If in the numerical simulations 
$\unicode[STIX]{x1D716}$ goes to 0 the electric field amplitude has to increase to infinity. The latter requires an unrestricted amount of power. Therefore, the steady state cannot be reached for 
$\unicode[STIX]{x1D716}=0$ at all, as was shown in Pavlenko et al. (Reference Pavlenko, Girka, Trush, Melnyk and Velizhanina2019).
However, in the case of the wave incidence from the vacuum the magnetic field amplitude in the plasma is not fixed by the antenna. For the underdense plasmas (
$\unicode[STIX]{x1D714}>\unicode[STIX]{x1D714}_{pl}$ or 
$0<\unicode[STIX]{x1D716}<1$) the magnetic field amplitude decreases and the electric field amplitude increases with time. If the incident wave amplitude is 1, then in the steady state the amplitude of the plasma magnetic field is 
$2\sqrt{\unicode[STIX]{x1D716}}/(1+\sqrt{\unicode[STIX]{x1D716}})$ (belonging to the range 
$[0;1]$) and the amplitude of the plasma electric field is 
$2/(1+\sqrt{\unicode[STIX]{x1D716}})$ (belonging to the range 
$[1;2]$). In other words, the electric field amplitude in the plasma tends to 2 when 
$\unicode[STIX]{x1D716}\rightarrow 0$. Therefore the build-up of the plasma particle oscillations requires much less power in comparison with the case of the metal antenna inside the plasma (for small 
$\unicode[STIX]{x1D716}$).
3.2 Case of electrostatic steady state field ($\unicode[STIX]{x1D716}=0$)
The formation of the electrostatic oscillation in the plasma is also studied numerically by the FDTD code. Since the frequency of the incident wave is chosen to be 1 GHz, the density of the hydrogen plasma is taken as 
$1.23976753\times 10^{10}~\text{cm}^{-3}$ to provide the case 
$\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{pl}$ exactly.
The analytical formulas (2.28), (2.30) describe exactly the magnetic and electric fields at any distance from the plasma–vacuum interface for any time. It is known from the discussions in the previous sections that the power flux into the plasma becomes less than 5 % of the incident power flux after 325 wave periods. That is why the FDTD code is used to simulate the electromagnetic field inside the plasma at this time but with the shift by a quarter of the wave period to get the maximum of the electric field. Thus, the electromagnetic field in the plasma is simulated for 
$t=325.25\,T$. Initially, the simulations were carried out with a spatial resolution of 50 points per vacuum wavelength and with the time resolution of 100 points per wave period but the results were not accurate enough. The reason for this is clear. The magnetic field as the exact analytical solution of the problem is an undifferentiable function of the coordinate at the plasma–vacuum interface, while the FDTD code operates with the differential operators in a numerical way. Therefore, the numerical differentiation at the plasma–vacuum interface introduces a computational error. If this explanation is correct then the computational error can be decreased by increasing the numerical grid resolution. That is why the resolution was increased several times. The results of the calculations are presented in figure 7 (for the magnetic field) and in figure 8 (for the electric field). The insets in the figures zoom in on the range of the small coordinates to see the differences in the results. The solid line presents the exact analytical solutions (2.28) and (2.30) respectively. The lines with numbers 100 (dash-dotted), 200 (dotted) and 500 (dashed) correspond to the simulations with spatial resolutions of 100, 200 and 500 grid points per vacuum wavelength respectively. The time resolution of the numerical grid was twice as large as the spatial one. One can see that the line 100 is sufficiently far from the analytical data but the line 500 is in good agreement with them. No doubt that a subgridding near the plasma–vacuum interface can essentially accelerate the numerical calculations; however, it is not realized in the current version of the FDTD code.
Figure 7. Magnetic field in plasma with 
$\unicode[STIX]{x1D716}=0$ normalized by the magnetic field amplitude of the incident wave as a function of the distance from the plasma–vacuum interface for the time 
$t=325.25\,T$. The solid line represents the exact analytical solution of the problem. Other lines are the results of numerical simulations by the FDTD code with different resolutions of the numerical grid (see text for details).
Figure 8. Electric field in plasma with 
$\unicode[STIX]{x1D716}=0$ normalized by the electric field amplitude of the incident wave as a function of the distance from the plasma–vacuum interface for the time 
$t=325.25\,T$. The solid line represents the exact analytical solution of the problem. Other lines are the results of numerical simulations by the FDTD code with different resolutions of the numerical grid (see text for details).
4 Conclusions
The exact formulas (2.20), (2.21) for the electromagnetic field transmitted into cold, unmagnetized and collisionless plasmas are written in the integral form for the case when the electromagnetic wave is incident normally upon the plasma half-space. The simplified formulas (2.22), (2.23) describe the electromagnetic field near the steady state propagation of the electromagnetic wave in the plasmas when the plasma particle oscillations are almost self-consistent with the electromagnetic field. The steady state propagation of the electromagnetic wave through the plasmas is known from plasma electrodynamics and does not require numerical simulation. The developed analytics make it possible to analyse the numerical data at transient times. The time of the transient processes (the time of the plasma transition from the initially quiet state to the state of the self-consistent plasma particle oscillations in the external electromagnetic field) can be easily calculated from (2.25). The analytics were developed for the sinusoidal incident wave defined by the equation (2.1). Of course, a change in time dependence of the incident signal can have a significant effect on the physics of the transient processes. But any long-term periodic signal can easily be expanded into the sinusoidal waves to be analysed by the developed analytics.
The formulas (2.20), (2.21) are also applicable when the frequency of the incident wave is equal to the plasma frequency. In this case they describe the build-up of the electrostatic oscillation in plasmas (2.28), (2.30). The last formulas provide one with the possibility to calculate the amplitude of the vanishing magnetic field and the amplitude of the forming electrostatic oscillation as functions of time and the distance from the plasma–vacuum interface.
If the frequency of the incident wave is equal to the plasma frequency, the power flux into the plasma drops below 1 % of the incident flux after approximately 8100 wave periods. The build-up of the plasma particle oscillations near the plasma–vacuum interface reduces the ability of the plasma to transmit the incident electromagnetic power, therefore, over time, an increasing part of the incident power is reflected back into the vacuum. As a result, the power flux into the plasma is not sufficient to build up the electrostatic oscillations at large distances from the plasma–vacuum interface. According to the theory (2.43), the distance where the amplitude of the electric field oscillation in the plasma is equal to the electric field amplitude of the incident wave, is proportional to the square root of the incidence time. Therefore, for example, these amplitudes are equal at the distance of 4 wavelengths after 250 wave periods, at the distance of 8 wavelengths after 1000 wave periods and so on.
Any plasma electromagnetic code operating in the time domain should be able to reproduce the analytical results. In the present paper the FDTD code is used to simulate two cases: the propagation of the electromagnetic wave through plasmas and the build up of the electrostatic oscillation in plasmas. In the first case the agreement between the analytical and numerical results is very good. However, the second case is much more difficult for the numerical simulations. Since the exact solution of the problem is an undifferentiable function of the coordinate at the plasma–vacuum interface, the numerical data of the FDTD code always contain the related numerical error. However, it is shown that the analytical and numerical results are the same if the resolution of the numerical grid is sufficiently large. One can conclude that an electromagnetic wave’s propagation through a sharp boundary between two media is a good challenge for the FDTD codes. It can be simulated, but additional attention should be paid to the boundary.
Acknowledgements
The authors would like to thank A. Kostenko from the University of Toronto for his assistance in preparing the paper for publication.
