Propagation Characteristics of High-Power Vortex Laguerre-Gaussian Laser Beams in Plasma

The propagation characteristics of high-power laser beams in plasma is an important research topic and has many potential applications in fields such as laser machining, laser-driven accelerators and laser-driven inertial confined fusion. The dynamic evolution of high-power Laguerre-Gaussian (LG) beams in plasma is numerically investigated by using the finite-difference time-domain (FDTD) method based on the nonlinear Drude model, with both plasma frequency and collision frequency modulated by the light intensity of laser beam. The numerical algorithms and implementation techniques of FDTD method are presented for numerically simulating the nonlinear permittivity model of plasma and generating the LG beams with predefined parameters. The simulation results show that the plasma has different field modulation effects on the two exemplified LG beams with different cross-sectional patterns. The self-focusing and stochastic absorption phenomena of high-power laser beam in plasma are also demonstrated. This research also provides a new means for the field modulation of laser beams by plasma.


Introduction
The complex interaction of high-power laser with plasma has become a popular research topic due to its various applications in science and technology, including laser-driven inertial confined fusion (ICF) [1], laser-driven accelerators [2] and laser machining [3].The main research areas have focused on the propagation characteristics of high-power laser beams in plasma [4][5][6].The Drude model of electrical conduction is extensively considered in this research field [7,8].Sharma et al. [9] have found that the propagation characteristics of an intense laser beam in plasma depend on the power and width of the beam and the ratio of plasma frequency and light wave frequency.Wang et al. [10] have investigated the propagation characteristics of a Gaussian laser beam in unmagnetized cold plasmas, based on the theory of ponderomotive nonlinearity.They also applied the Wentzel-Kramers-Brillouin (WKB) method to study the propagation characteristics of a Gaussian laser beam in cold plasma [11].However, since the above studies are conducted by only solving the scalar wave equation with complex eikonal function assumption and paraxial approximation, the obtained information about the dynamic laser-plasma interaction is limited.Additionally, the solutions are not accurate enough for tightly focused laser beams whose sizes are in the order of several wavelengths.
In recent decades, the finite-difference time-domain (FDTD) method has become a popular numerical technique in directly solving Maxwell's curl equations [12].The FDTD method also has well-developed techniques for accurately simulating the Drude model [13][14][15] and numerically generating any types of laser beams based on the total-field/scattered-field (TF/SF) technique [16].
Therefore, the current study evaluates the FDTD method in investigating the propagation characteristics and relevant phenomena of high-power Laguerre-Gaussian laser beams in plasma.

Theoretical Background
In physical optics, the complex-form expression of the electric field of a Laguerre-Gaussian beam propagating in the free space under cylindrical coordinates (r, φ, z) is given by: where E 0 is the nominal magnitude of light field of laser beam, k = 2π/λ L is the wave number of laser beam with wavelength λ L , w 0 is the radius of beam waist, w(z) is the beam size at distance z from the beam waist with Rayleigh range z R = πw 2 0 /λ, q(z) = z + iz R is the complex beam parameter, ϕ pl (z) = (2p + |l| + 1) tan −1 (z/z R ) is the Gouy phase shift, and L l p (•) are the associated Laguerre polynomials.Herewith, each Laguerre-Gaussian beam of LG l p mode is specified by two mode indices, the angular mode index l (l = 0, ±1, ±2, . ..) and the transverse radial mode index p (p = 0, 1, 2, . ..).It is noted that a Laguerre-Gaussian beam degenerates to the fundamental Gaussian beam when p = l = 0.In this work, we focus on the category of vortex Laguerre-Gaussian laser beams with l = 0. To maintain simplicity without losing generalizability, we investigate the propagation of two typical low-order Laguerre-Gaussian beams in plasma.The first one is the doughnut-shaped LG 1 0 laser beam with electric field written as: Figure 1 shows the field magnitude and phase distributions of the LG 1 0 laser beam propagating in free space, where E m (r, z) = E l p (r, z) and Φ(r, z) = arg E l p (r, z) are the magnitude and principal argument of the complex electric field E l p (r, z) of laser beam given by Equation (1), respectively.The size of beam waist is w 0 = 3λ L as depicted in Figure 1, where λ L is the wavelength of laser beam.The distinct doughnut-shaped feature makes it ideal for applications in many areas of optics, such as optical trapping [17] and optical tweezers [18].
The second laser beam is the double-ring-shaped Laguerre-Gaussian vortex laser beam of LG 2 1 mode with its electric field given by according to the associated Laguerre polynomial L 2 1 (x) = 3 − x. Figure 2 shows the magnitude and phase distributions of an electric field of the Laguerre-Gaussian laser beam of LG 2  1 mode propagating in free space.Here, a smaller size of beam waist w 0 = 2λ L is applied, which implies a larger diffraction effect as we can see from Figure 2a.The cross-sectional view of field magnitude distribution E m (r, z = 0)/E 0 in Figure 2b shows that the laser beam is double-ring-shaped.LG mode in free space: The dynamic interaction of laser beams with plasma is complex and dependent on the level of laser power [19].When the intensity of the incident laser beam is low, the charged particles such as ions and electrons in plasma driven by the light fields are oscillating around their nearly-fixed oscillation centers.Therefore, the spatial distribution of the electron density keeps nearly unchanged.However, when the light intensity is extremely high, the spatially-inhomogeneous light beam exerts a significant nonlinear Lorentz force, called the ponderomotive force, on the charged particles in the plasma.Since the ponderomotive force scales with the inverse of particle mass, the ponderomotive effect on ions is generally negligible with respect to that of electrons.Thus, the   LG mode in free space: The dynamic interaction of laser beams with plasma is complex and dependent on the level of laser power [19].When the intensity of the incident laser beam is low, the charged particles such as ions and electrons in plasma driven by the light fields are oscillating around their nearly-fixed oscillation centers.Therefore, the spatial distribution of the electron density keeps nearly unchanged.However, when the light intensity is extremely high, the spatially-inhomogeneous light beam exerts a significant nonlinear Lorentz force, called the ponderomotive force, on the charged particles in the plasma.Since the ponderomotive force scales with the inverse of particle mass, the ponderomotive effect on ions is generally negligible with respect to that of electrons.Thus, the The dynamic interaction of laser beams with plasma is complex and dependent on the level of laser power [19].When the intensity of the incident laser beam is low, the charged particles such as ions and electrons in plasma driven by the light fields are oscillating around their nearly-fixed oscillation centers.Therefore, the spatial distribution of the electron density keeps nearly unchanged.However, when the light intensity is extremely high, the spatially-inhomogeneous light beam exerts a significant nonlinear Lorentz force, called the ponderomotive force, on the charged particles in the plasma.Since the ponderomotive force scales with the inverse of particle mass, the ponderomotive effect on ions is generally negligible with respect to that of electrons.Thus, the electron density in plasma is manifestly changed by such a force until the balance is restored with the plasma pressure gradient force.According to the theory of ponderomotive nonlinearity, the redistributed electron density is exponentially dependent on the square of light-field magnitude E 2 m (r, z) and is given by [20]: where n e0 (r, z) is the initial spatial distribution of electron density before the presence of laser beam, E 2 m (r, z) is the square of light-field magnitude, T e is the electron temperature of plasma in unit of kelvin, k B is Boltzmann's constant, ω L is the angular frequency of laser beam, and the coefficient ) is defined for the sake of conciseness of following equations.On the other hand, the electromagnetic properties of plasma are usually described by the famous Drude model of permittivity: where ε 0 is permittivity of free space, ω p is the plasma frequency and ν c is the collision frequency.
Theoretical derivations [7,11] demonstrated that ω 2 p and ν c both are proportional to the electron density n e (r, z).Therefore, for high-power laser beam incidence, both the plasma frequency and collision frequency are modulated by the square of field magnitude as follows: and where ω p0 = e √ n e0 /ε 0 m e and ν c0 are the initial plasma frequency and collision frequency before the incidence of laser beam.Thus, the modified Drude model of permittivity for describing the interaction of high-power laser beams with plasma is given by As apparent in Equation ( 8), such a physical model of permittivity for plasma is related to the light intensity of laser beam.Therefore, it is a nonlinear and space-dependent model.According to the electromagnetic theory, the refractive index of a medium is n = √ ε r µ r , where ε r and µ r are the relative permittivity and relative permeability of the medium, respectively.Since the relative permeability of plasma is usually defaulted to unit one µ r = 1 in the regime of optical frequencies, the complex refractive index of plasma is accordingly expressed as which also reflects the dispersion and dissipation characteristics of the plasma.

Numerical Methodology
The Laguerre-Gaussian laser beams are numerically generated by utilizing the FDTD method based on the total field/scattering field (TF/SF) source condition [12,16].It is possible to generate any types of laser beams with predefined pattern and parameters using this technique.Within the framework of FDTD method, assuming that the TF/SF interface is perpendicular to the z axis and locates at z = (k s − 1/4)∆z, the x component of scattering magnetic field H (n+1/2) s,x on the plane z = (k s − 1/2)∆z at time instant t = (n + 1/2)∆t can be updated by: LG,y (i, j + where the first item on the right side of ( 10) is the FDTD updating equation for H and the last item is the y component of the electric field of the to-be-generated laser beam, LG,y (i, j + Similarly, the y component of scattering magnetic field: LG,x (i + where the last item is the x component of the electric field of the to-be-generated laser beam: LG,x (i + Typically, the TF/SF interface is set at or near the plane of beam waist for convenience purposes, where the electric field has the simplest expression and there is no z component of electric field. The bilinear transform (BT) approach [13] can be used to implement the complex dispersive model given by Equation ( 8) via the numerical implementation of the constitutive equation, D = ε(ω)E, within the frame work of FDTD method.The BT approach is not only accurate, but also stable with a stability limit that is equal to the Courant stability limit, ∆t C = ∆s/c √ m, where ∆s = min{∆x, ∆y, ∆z} is the shortest side of Yee cell, c is the speed of light, and m is the dimensionality of the simulated problem.In fact, we can first introduce the auxiliary quantity: so that: Subsequently, by applying the bilinear transformation jω∆t ⇔ 2(1 − Z −1 )/(1 + Z −1 ) to transform Equation ( 15) from the frequency domain to the Z domain and noting that E n Z −1 = E n−1 and S n Z −1 = S n−1 , we have the updating equation for the discrete time-domain relationship between S and E: where /(2ν c0 ∆tβ + 4), and C 3 = (−2ν c0 ∆tβ + 4)/(2ν c0 ∆tβ + 4) with β = exp −αE 2 m (r, z) .By substituting Equation ( 16) into Equation ( 15), we have the final updating equation for E: after some mathematical manipulation.However, the coefficients C 1 , C 2 , and C 3 are dependent on β that is related to E 2 m (r, z).Thus, we also need to extract the parameter E 2 m (r, z).According to the time-harmonic property of a coherent laser beam, E 2 m (r, z) can be numerically extracted by applying the composite trapezoidal rule within the discrete framework of FDTD method [20], where N T = T L /∆t is an even number of time steps per period of laser oscillation.

Simulation Results
The propagation characteristics of the x-polarized high-power LG 1 0 laser beam in plasma is first modeled and simulated.The laser beam used is as the one depicted in Figure 1, where the beam waist w 0 = 3λ L with the wavelength λ L = 351nm in the ultraviolet region.The nominal electric field of the laser beam is E 0 = 2.55 × 10 10 V/m, which is feasible and practicable in modern giant laser facilities [1].The total three-dimensional computational region for the following FDTD simulations is 12λ L × 12λ L × 28λ L in the x, y, and z directions and is surrounded by perfectly matched layers (PMLs).A cubic Yee cell with size ∆x = ∆y = ∆z = λ L /20 is applied for the space discretization.The initial electron density of simulated plasma is spatially homogeneous with n e0 = 3.25 × 10 21 cm −3 and the collision frequency is ν c0 = 0.05ω p0 .As the interaction process is simulated in a relatively short time t = 3000∆t = 3000∆x/2c = 8.78 × 10 −14 s for a total of 3000 time steps, a constant electron temperature T e = 1 × 10 4 K is also considered. Figure 3a,b show the simulation results for the longitudinal views of the redistributed light field magnitude of high-power doughnut-shaped LG 1 0 laser beam on the xOz and yOz planes, respectively, when it propagates in plasma.Figure 4a-f show the cross-sectional views of light field magnitude of LG 1 0 laser beam at several propagation distances z = 4λ L , 8λ L , 12λ L , 16λ L , 20λ L and 24λ L from the beam waist.The dynamic interaction process between the LG 1 0 laser beam and plasma are vividly recorded in Videos S1 and S2.See Supplementary Videos S1 and S2 for details.
to the time-harmonic property of a coherent laser beam, 2 m ( , ) E r z can be numerically extracted by applying the composite trapezoidal rule within the discrete framework of FDTD method [20], where is an even number of time steps per period of laser oscillation.

Simulation Results
propagation characteristics of the x -polarized high-power LG laser beam in plasma is first modeled and simulated.The laser beam used is as the one depicted in Figure 1, where the beam waist in the ultraviolet region.The nominal electric field of the laser beam is  LG laser beam on the xOz and yOz planes, respectively, when it propagates in plasma.Figure 4a-f show the cross-sectional views of light field magnitude of LG laser beam at several propagation distances LG laser beam and plasma are vividly recorded in Videos S1 and S2.See Supplementary Videos S1 and S2 for details.Secondly, the propagation characteristics of the x-polarized high-power LG 2  1 laser beam in plasma is modeled and simulated.The investigated laser beam is as the one depicted in Figure 2 and the parameters for wavelength and light-field intensity are the same as those of the LG 1 0 laser beam.The parameters for FDTD simulations and plasma parameters are also the same as the first numerical example.Figure 5a,b show the simulation results for the longitudinal views of the redistributed light field magnitude of high-power double-ring-shaped LG 2  1 laser beam on the xOz and yOz planes, respectively, when it is propagating in plasma.Figure 6a-f show the cross-sectional views of light field magnitude of LG 2  1 laser beam at several specified propagation distances z = 4λ L , 8λ L , 12λ L , 16λ L , 20λ L and 24λ L from the beam waist.The dynamic interaction process between the LG 2  1 laser beam and plasma are vividly recorded in Videos S3 and S4.See Supplementary Videos S3 and S4 for details.LG laser beam in plasma at the several specified propagation distances (a) Secondly, the propagation characteristics of the x -polarized high-power LG laser beam in plasma is modeled and simulated.The investigated laser beam is as the one depicted in Figure 2 and the parameters for wavelength and light-field intensity are the same as those of the LG laser beam.The parameters for FDTD simulations and plasma parameters are also the same as the first numerical example.Figure 5a,b show the simulation results for the longitudinal views of the redistributed light field magnitude of high-power double-ring-shaped LG laser beam on the xOz and yOz planes, respectively, when it is propagating in plasma.Figure 6a-f  LG laser beam at several specified propagation distances LG laser beam and plasma are vividly recorded in Videos S3 and S4.See Supplementary Videos S3 and S4 for details.LG laser beam in plasma at the several specified propagation distances (a) Secondly, the propagation characteristics of the x -polarized high-power LG laser beam in plasma is modeled and simulated.The investigated laser beam is as the one depicted in Figure 2 and the parameters for wavelength and light-field intensity are the same as those of the LG laser beam.The parameters for FDTD simulations and plasma parameters are also the same as the first numerical example.Figure 5a,b show the simulation results for the longitudinal views of the redistributed light field magnitude of high-power double-ring-shaped LG laser beam on the xOz and yOz planes, respectively, when it is propagating in plasma.Figure 6a-f  LG laser beam at several specified propagation distances LG laser beam and plasma are vividly recorded in Videos S3 and S4.See Supplementary Videos S3 and S4 for details.

Discussion
From Figures 3 and 5, it is noted that the self-focusing and absorption phenomena of high-power laser beams in plasma are evident as compared with those propagating in free space as illustrated in Figures 1 and 2. The physical model of permittivity with ponderomotive nonlinearity given by Equation ( 8) determines that the refractive index of plasma is smaller at the place where the light intensity is lower and larger where the light intensity is stronger.The ensuing inhomogeneous distribution of refractive-index attribute to the self-focusing effect and the laser beam is therefore focused when this effect overly counteracts the natural diffractive divergence of the laser beam.This is apparent from Figures 4 and 6 in which the cross-sectional profiles of the light fields of LG 1 0 and LG 2 1 vortex laser beams within the plasma first experienced self-focusing and then stochastic absorption.Video S2 and Video S4 vividly displayed the rotations and temporal evolutions of the cross-sectional light fields of LG 1 0 and LG 2 1 laser beams, which demonstrate the vortex features of the LG laser beams.We also note that the field distributions of the xOz and yOz planes are different from each other since the two laser beams are both linearly polarized in the x direction.
We also highlight the ribbons with half-wavelength intervals in Figures 3 and 5. Since the incident laser beams are assumed spatially coherent, the transient light intensity proportional to E 2 (r, z, t) is also harmonic with a half-wavelength period.In each half period of space, the electrons at the place with high transient light intensity are pushed to the place with low transient light intensity.According to Equations ( 4) and ( 9), this intensifies the difference of refractive indices at the two places and the half-wavelength electron-free elongated cavities are formed.The ensuing high reflection of light at the walls of cavities causes the ribbon phenomenon.Moreover, the speckles found in Figures 4 and 6 can be attributed to the Rayleigh-Taylor fluid mechanics instability [19] under the local oscillation electronic acceleration mechanism.Once the initial laser beam has a small-scale perturbation or is disturbed by the unevenly medium when it propagates in the plasma, a rapid light intensity growth will occur in the near region that forms a speckle in the cross-sectional patter of the laser beam.It can be estimated that different modulation effects of light fields can be achieved by using LG laser beams of other polarization states, such as the radial and azimuthal polarizations.Thus, the plasma can serve as a new medium for modulating the LG beams into the desired field patterns.

Conclusions
In this article, the propagation characteristics of Laguerre-Gaussian vortex laser beams in plasma is modeled and simulated under the given parameters.By using the FDTD method, the numerical methodology for modeling the nonlinear permittivity and generating the LG beams is proposed.The simulation results about the 3D dynamic interaction of high-power LG beams with plasma are illustrated and discussed.The obtained simulation results verify the self-focusing and stochastic absorption phenomena of high-power laser beam propagation in plasma.Our work presents the methodology for the numerical simulations of the complex interaction of high-power Laguerre-Gaussian vortex laser beams with plasma by using the FDTD method.It provides detailed information about the dynamic 3D evolution of both the light field of the laser beam and the spatial distribution of electron density and many special laser beams can be simulated based on the TF/SF technique.This is much superior to the state-of-the-art studies [9][10][11] that were conducted only by solving the scalar wave equation based on the WKB method with complex eikonal function assumption and paraxial approximation, where only a simple laser beam of large beam size can be calculated and only the limited information such as the variation of beam size can be obtained.This research embraced the potential merits of theory models and simulation techniques for many research applications such as the laser machining and laser-driven inertial confined fusion.

Figure 1 .
Figure 1.The field magnitude and phase distributions of the doughnut-shaped Laguerre-Gaussian vortex laser beam of

Figure 2 .
Figure 2. The field magnitude and phase distributions of the double-ring-shaped Laguerre-Gaussian vortex laser beam of (a) The longitudinal view of the normalized field magnitude distribution m 0 ( , ) / E r z E ; (b) the cross-sectional view of the normalized field magnitude distribution field of laser beam on the plane of beam waist.

Figure 1 .
Figure 1.The field magnitude and phase distributions of the doughnut-shaped Laguerre-Gaussian vortex laser beam of LG 1 0 mode in free space: (a) The longitudinal view of the normalized field magnitude distribution E m (r, z)/E 0 ; (b) the cross-sectional view of the normalized field magnitude distribution E m (r, z = 0)/E 0 on the plane of beam waist; (c) the cross-sectional view of principal argument distribution Φ(r, z = 0) of the complex field of laser beam on the plane of beam waist.

Figure 1 .
Figure 1.The field magnitude and phase distributions of the doughnut-shaped Laguerre-Gaussian vortex laser beam of

Figure 2 .
Figure 2. The field magnitude and phase distributions of the double-ring-shaped Laguerre-Gaussian vortex laser beam of (a) The longitudinal view of the normalized field magnitude distribution m 0 ( , ) / E r z E ; (b) the cross-sectional view of the normalized field magnitude distribution m field of laser beam on the plane of beam waist.

Figure 2 .
Figure 2. The field magnitude and phase distributions of the double-ring-shaped Laguerre-Gaussian vortex laser beam of LG 2 1 mode in free space: (a) The longitudinal view of the normalized field magnitude distribution E m (r, z)/E 0 ; (b) the cross-sectional view of the normalized field magnitude distribution E m (r, z = 0)/E 0 on the plane of beam waist; (c) the cross-sectional view of principal argument distribution Φ(r, z = 0) of the complex field of laser beam on the plane of beam waist.
is also considered.Figure3a,bshow the simulation results for the longitudinal views of the redistributed light field magnitude of high-power doughnut-shaped 1 0 beam waist.The dynamic interaction process between the 1 0

Figure 3 .
Figure 3. Simulation results for the propagation of the high-power doughnut-shaped

Figure 3 .
Figure 3. Simulation results for the propagation of the high-power doughnut-shaped LG 1 0 laser beam in plasma (see Video S1).The longitudinal views of the light-field magnitude of laser beam normalized to E 0 and shown in logarithmic scale (a) on the xOz plane and (b) on the yOz plane.

Figure 4 .
Figure 4.The cross-sectional views of the light-field magnitude profiles of the high-power doughnut-shaped waist (see Video S2).The colormap bars of all the sub-figures are the same as that of Figure3.
show the cross-sectional views of light field magnitude of2 1

Figure 5 .
Figure 5. Simulation results for the propagation of the high-power double-ring-shaped

Figure 4 .
Figure 4.The cross-sectional views of the light-field magnitude profiles of the high-power doughnut-shaped LG 1 0 laser beam in plasma at the several specified propagation distances (a) z = 4λ L , (b) z = 8λ L , (c) z = 12λ L , (d) z = 16λ L , (e) z = 20λ L and (f) z = 24λ L from the beam waist (see Video S2).The colormap bars of all the sub-figures are the same as that of Figure 3.

Figure 4 .
Figure 4.The cross-sectional views of the light-field magnitude profiles of the high-power doughnut-shaped Video S2).The colormap bars of all the sub-figures are the same as that of Figure3.
show the cross-sectional views of light field magnitude of2 1 the beam waist.The dynamic interaction process between the 2 1

Figure 5 .
Figure 5. Simulation results for the propagation of the high-power double-ring-shaped

Figure 5 .
Figure 5. Simulation results for the propagation of the high-power double-ring-shaped LG 2 1 laser beam in plasma (see Video S3).The longitudinal views of the light-field magnitude of laser beam normalized to E 0 and shown in logarithmic scale (a) on the xOz plane and (b) on the yOz plane.

Figure 6 .
Figure 6.The cross-sectional views of the light-field magnitude profiles the high-power double-ring-shaped LG 2 1 laser beam in plasma at the several specified propagation distances (a) z = 4λ L , (b) z = 8λ L , (c) z = 12λ L , (d) z = 16λ L , (e) z = 20λ L and (f) z = 24λ L from the beam waist (see Video S4).The colormap bars of all the sub-figures are the same as that of Figure 3.