2. Equations

Main considerations

In the optical propagation calculations we describe evolution of electromagnetic field from a state defined at some time and space region to another state, that this field takes after some time and in a different space region.

This problem can be described by the Maxwell equations:

$$\begin{aligned} & \nabla \times \mathbf{E} = - \partial_t \mathbf{B}\,,\\\ & \nabla \times \mathbf{H} = \partial_t \mathbf{D} + \mathbf{J}\,,\\\ & \nabla \mathbf{D} = \rho\,,\\\ & \nabla \mathbf{B} =0\,, \end{aligned}$$

accompanied by the constitutive equations

$$\begin{aligned} & \mathbf{D} = \epsilon_0 \mathbf{E} + \mathbf{P} \,,\\\ & \mathbf{H} = \cfrac{1}{\mu_0} \mathbf{B} - \mathbf{M} \,. \end{aligned}$$

In the following we will assume that all charges of the system are free, i.e. $\epsilon_0$ and $\mu_0$ are the electric permittivity and magnetic permeability of vacuum. Moreover, at the current stage of development we will only consider plasma processes and discard the polarization and magnetization of the media, thus setting $\mathbf{P}=0$ and $\mathbf{M}=0$.

Full equation

In the optical field the electric and magnetic components are consistent with each other, and in most practical situations it is enough to follow the electric component, for which the equation can be derived as:

$$\begin{aligned} & \nabla^2 \mathbf{E} - \frac{1}{c^2} \; \partial_t^2 \mathbf{E} = \mu_0\; \partial_t\mathbf{J} + \frac{1}{\epsilon_0}\nabla \rho\,. \end{aligned}$$

Here the left hand side is the wave equation which describes propagation of the initial field in vacuum (the optical beam whose dynamics we study). The terms in the right hand side are the currents and density modulations, which act as the sources for the field generated by the media as a response to the optical field.

Source terms

The two source terms are essentially different in nature. The first one is excited instantaneously as electrons generate electric currents being moved by the field. The second term constitutes the charge density modulations, and it presents a slower response and is typically considered to be much weaker. Let us demonstrate this by assuming the charge continuity in plasma and re-writing the source terms in the integral form:

$$\begin{aligned} & \mu_0 \; \partial_t\mathbf{J} + \frac{1}{\epsilon_0}\nabla \rho = \frac{1}{\epsilon_0 } \int_{-\infty}^{t} dt' \left(\frac{1}{c^2} \; \partial_t^2 \mathbf{J} - \nabla (\nabla \mathbf{J}) \right)\,. \end{aligned}$$

In the Fourier space, the terms under the integral appear as $k^2 \hat{\mathbf{J}}$ and $\mathbf{k}\cdot (\mathbf{k} \cdot \hat{\mathbf{J}})$, where $k=\omega/c$, and for the polarized field, it is enough to consider a single transverse component (e.g. $\mathbf{J} \parallel \mathbf{E} \parallel \mathbf{e}_x$) and discard others. The ratio between source terms therefore reads, $k^2/k_x^2$, and for the field wavelength $\lambda_0$, and transverse size $R$, it scales as $\propto (R/\lambda_0)^2$. It is easy to see, that in most practical situations this indeed allows to assume $\mu_0 \; \partial_t\mathbf{J} \gg \frac{1}{\epsilon_0}\nabla \rho$, and discard contribution of the charge density modulations compared to the induced currents.

Working equation

To conclude this part, in the calculations we will consider the equation for the electric field in the following form:

$$\begin{aligned} & \nabla^2 \mathbf{E} - \frac{1}{c^2} \partial_t^2 \mathbf{E} = \mu_0 \partial_t\mathbf{J} \,, & \qquad (1) \end{aligned}$$

and moreover we will use it as a scalar equation, for the field and current components along the laser field polarization.

Vacuum case

In a case when no sources are present in the system, i.e. field propagates in vacuum, the equation simplifies down to the uniform wave equation

$$\begin{aligned} & \left(\nabla^2 - \frac{1}{c^2} \; \partial_t^2\right) E = 0 \end{aligned}$$

Plasma dispersion

A simplest linear plasma response can be described with the help of non-relativistic motion equation:

$$\begin{aligned} & \mu_0 \partial_t J = -e \; \mu_0 \;n_{pe} \; \partial_t v = \omega_{pe}^2 / c^2 E \end{aligned}$$

where $n_{pe}$ is a number density of plasma electrons, and $\omega_{pe} = \sqrt{\frac{e^2 n_{pe} }{\epsilon_0 m_e}}$ is a plasma frequency.

Introducing this term in the optical equation leads to the following equation:

$$\left(c^2 \nabla^2 - \partial_t^2 \; - \; \omega_{pe}^2 \right) \; E = 0\,,$$

which if considered in the spatio-temporal Fourier domain leads to the well-known dispersion relation

$$\omega^2 = k^2 c^2 + \omega_{pe}^2$$

Field representations

When approaching the computational problem of optical propagation we typically need to solve Eq. (1) with an initial condition that represents the initial beam state. There are two standard ways to define and consider the field.

Spatial representation

In so-called spatial form, the field is defined at a time $t_0$ in a full space as $E_0(x,y,z)=E(x,y,z, t_0)$, and we look for the state after some time pass, so we calculate, $E_1(x,y,z) = E(x,y,z, t_0+\Delta t)$. The spatial domain where $E_1$ is localized is usually adjusted to the propagated distance, typically $z\to z+c\Delta t$ (we will always consider propagation along $z$-axis). This approach is typical for the particle-in-cell (PIC) codes, where electromagnetic solver is couples with the charges particles motion, that can lead to the complex charge and current density structures in the spatial domain.

Temporal representation

Here the field is defined in the spatial plane located at $z_0$ as a function of time $E_0(x, y, t)=E(x, y, z_0, t)$, and we look at its state at a different location $E_1(x,y,t) = E(x,y,z_0+\Delta z, t)$ at a later time when the field arrives there, e.g. $t \to t+\Delta z/c$. In this approach, media is presented by an oscillating plane, and it is more convenient be used with the temporal spectral decomposition of the field, while resolving the longitudinal dynamics of the plasma on a scale of the pulse length $c\tau$, may be more challenging.

While both approaches are mathematically equivalent, they would lead to different implementations, and usually serve different purposes. The spatial representation is used in PIC simulation as it is giving better description of the excited plasma structures and can be easily extended to the more rigorous account of the source terms. Temporal representation is more common for the optical propagation problems, and in the following we will always be considering it as a default one. In vacuum without source terms, the conversion between representations can be easily done using propagation techniques.

Propagation equation with temporal representation

Let us calculate the field $E_1(x,y,t) = E(x,y,z_0+\Delta z, t)$ for a known $E_0(x, y, t)=E(x, y, z_0, t)$, i.e. considering $z$ the propagation variable, and the field being defined in the temporal domain. For this we divide the Laplace operator in Eq(1) into longitudinal and transverse parts:

$$\begin{aligned} & \partial_z^2 \; \mathbf{E} = \frac{1}{c^2} \;\partial_t^2 \;\mathbf{E} \;-\; \nabla_\perp^2 \;\mathbf{E} + \mu_0 \; \partial_t \; \mathbf{J} \,, \qquad (2) \end{aligned}$$

where $\nabla_\perp^2 = \partial_x^2 + \partial_y^2$ for the problem defined in the Cartesian geometry $(x,y,t)$, and $\nabla_{\perp}^2 = \cfrac{1}{r}\; \partial_r (r \; \partial_r) + \cfrac{1}{r^2} \; \partial_\theta^2$ in the cylindrical coordinates $(r, \theta, t)$.

Spectral methods

All our computations will adopt the so-called spectral methods. In the optical problems this commonly involves presentation of the field and the source terms as the sums of the waves both temporally and in the transverse plane.

For the cartesian geometry such decomposition is well-known as the spatio-temporal Fourier series:

$$\begin{aligned} & A (x,y,t) = Re \left [ \sum_{ k_x, k_y, \omega} \; \hat{A}(k_x, k_y, \omega) \; \exp(ik_x x + ik_y y-i \omega t) \right ]\,, \end{aligned}$$

and for the cylindrical coordinates we define is as the Fourier-Bessel (or Fourier-Hankel) series:

$$\begin{aligned} & A (r,\theta,t) = Re \left [ \sum_{k_r, m, \omega} \; \hat{A}(k_r, m, \omega) \; J_m(k_r r) \; \exp(-i m \theta-i \omega t) \right ] \end{aligned}$$

The properties of such representations are defined by the sampling that we choose in the real and spectral spaces. For the case of spatio-temporal Fourier series, it is well known that uniform sampling in space

Both these representations transform the first two terms in the right hand side of Eq(2) into:

$$\frac{1}{c^2} \;\partial_t^2 \;\mathbf{E} \;-\; \nabla_\perp^2 \;\mathbf{E} \to -(\omega^2 - k_r^2) \hat{E}$$