Branched flow is an interesting phenomenon that can occur in diverse systems. It is usually linear in the sense that the flow does not alter the properties of the medium. Branched flow of light on thin films has recently been discovered. It is therefore of interest to know whether nonlinear light branching can also occur. Here, using particle-in-cell simulations, we find that in the case of an intense laser propagating through a randomly uneven medium, cascading local photoionization by the incident laser, together with the response of freed electrons in the strong laser fields, triggers space–time-dependent optical unevenness. The resulting branching pattern depends dramatically on the laser intensity. That is, the branching here is distinct from the existing linear ones. The observed branching properties agree well with theoretical analyses based on the Helmholtz equation. Nonlinear branched propagation of intense lasers potentially opens up a new area for laser–matter interaction and may be relevant to other branching phenomena of a nonlinear nature.

## I. INTRODUCTION

Waves propagating in a randomly uneven medium with correlation length greater than the wavelength *λ* can form filaments in a tree-branch-like manner.^{1} This phenomenon, known as branched flow, has been observed in diverse systems with different length scales.^{2–20} For example, instead of smoothly diffusing or spreading, an electron beam passing through a two-dimensional (2D) electron gas can form branching strands that become successively narrower,^{2–6} tsunami waves are attributed to random unevenness on the ocean floor,^{7–10} and branching of microwave radiation emitted by pulsars is attributed to interstellar dust clouds.^{12–14} Recently, branched flow of light (specifically continuous-wave laser light) has been found on thin soap films.^{17,18} This light branching is attributed to variations in the film’s refractivity *η*, which bend and bundle the light rays at favorable locations and form caustics.^{17–20} Despite its complexity, branched flow is usually a linear (passive) phenomenon in the sense that the flow does not alter the properties of the medium. That is, the properties of the medium remain space–time-independent, despite the presence of the flow. In this regime, the distance *d*_{0} from the source to the first branching point follows the scaling $d0\u221dlcv0\u22122/3$,^{1,15,21–24} where *l*_{c} is the correlation length of the medium’s unevenness, and *v*_{0} is its strength parameter (to be defined later).

With the development of chirped pulse amplification technology,^{25} lasers with intensity *I* ranging from 10^{14} to 10^{20} W/cm^{2} are readily available. The propagation of such intense lasers through matter initiates new phenomena and provides opportunities for a wide range of applications.^{26} When *I* ≳ 10^{14} W/cm^{2}, the atomic Coulomb barrier is suppressed by the strong laser electric field, electrons are set free, and the affected medium is ionized into a plasma,^{27} whose optical properties then become dominated by electron dynamics. Moreover, at higher laser intensities *I* ≳ 1.37 × 10^{18} W/cm^{2}, namely, when the laser intensity is above the relativistic threshold, in addition to photoionization, the laser ponderomotive force and relativistic plasma motion can significantly modify the original unevenness in the density as well as the local refractivity, thereby affecting the laser propagation.^{26} Whether branched flow of intense laser light can occur, and, if so, how it evolves, in such a space–time-dependent medium involving complex nonlinear effects remains unknown.

In this article, we present an investigation of nonlinear branched flow of intense laser light in uneven media. Since the laser propagation behavior determines the interaction region and can greatly influence the energy coupling efficiency, such an investigation may provide some fundamental insights into laser–matter interaction. Our particle-in-cell (PIC) simulations show that laser branching can occur at moderate laser intensities (*I* ∼ 10^{14}–10^{17} W/cm^{2}) in an inhomogeneous plasma with a randomly uneven density distribution. In contrast to linear light branching, in this regime, the branching depends crucially on the laser intensity. In particular, photoionization induced by the strong laser electric field raises the density unevenness along the laser paths and enhances branching. However, relativistic lasers can suppress branching by smoothing the unevenness and thus the local refractivity of the plasma. An analysis of the branching process and the resulting properties consistent with the simulation results is also given. Our work sheds light on the nonlinear dynamics of branched flow in space–time-dependent complex media. Branched flow of intense laser light, as a new mode for laser–matter coupling, may be of use in many fields, such as generation of intense quasi-monochromatic partially spatially incoherent light for nonlinear optics and soliton science.

## II. SIMULATION METHOD

Light branching is usually three-dimensional (3D). However, if the irregularity of the uneven background medium is isotropic, branching can be angular independent with respect to the laser axis. This enables us to conduct 2D simulations of a certain plane containing the laser axis without loss of generality. The 3D effects will be discussed later. In the PIC simulations, the initial background medium [see Fig. 1(a)] is a weakly pre-ionized SiO_{2} plasma (with Si^{2+} and O^{+} ions) with uneven density distribution located in 0 *µ*m < *x* < 215 *µ*m, −55 *µ*m < *y* < 55 *µ*m. The density unevenness has an isotropic correlation length *l*_{c} = 4.8 *µ*m, as obtained from the autocorrelation function (ACF)^{28} shown in Fig. 1(d) (see the Appendix for details). The average densities of Si^{2+}, O^{+}, and electrons are 0.02*n*_{c}, 0.04*n*_{c}, and 0.08*n*_{c}, respectively, where *n*_{c} ∼ 1.1 × 10^{21} $cm\u22123/\lambda 02$ is the critical density. The coefficient of density variation is ∼30%. An intense circularly polarized Gaussian laser pulse of central wavelength *λ*_{0} = 1.06 *µ*m, peak intensity *I*_{0} = 10^{14} W/cm^{2}, and FWHM spot size *r*_{0} = 16.65 *µ*m is incident normally from the left boundary. The laser pulse duration is 1.5 ps, with a flat-top time profile for simplicity. The effects of photoionization during laser–plasma interaction are self-consistently included in the simulations^{29} (see the Appendix and supplementary material for details). The uneven plasma can be created by ionizing appropriate chemically fabricated low-density foams with an extra preceding hundred-picosecond laser.^{30,31} We note that similar plasmas with uneven density distribution occur widely in intense laser–matter interaction, for example, in laser ablation on gas clusters^{32} and in backward Raman amplification schemes for intense-laser pulse compression.^{33}

## III. BRANCHED FLOW PATTERNS FOR INTENSE LASER LIGHT AT MODERATE INTENSITIES

Figure 1(b) for the distribution of the laser intensity at *t* = 1095 fs clearly shows light branching: the laser breaks up into several filaments after the first caustics at *x* ∼ 40 *µ*m. As they propagate, the filaments break up further into narrower and weaker ones in a cascade manner (see the movie in the supplementary material). To the best of our knowledge, such bifurcations of intense lasers have not been reported before. A useful quantity for characterizing the branching pattern is the scintillation index Σ = (⟨*I*^{2}⟩/⟨*I*⟩^{2}) − 1, which measures the relative strength of intensity fluctuations.^{11} For statistical accuracy, here Σ is obtained from simulations using a reference plane wave with otherwise identical interaction parameters. We see in Fig. 1(e) that Σ increases from 0 at the plasma front surface to a local maximum of 0.42 at *x* ∼ 40 *µ*m. The dependence of Σ on *x* agrees well with the branching pattern of the laser intensity shown in Fig. 1(b).

It is worth mentioning that the laser power here is below the self-focusing threshold, namely, *P*_{cr}[GW] ∼ 16.2*n*_{c}/*n*_{e}, with *n*_{e} the electron density.^{34–37} Thus, the conventional laser filamentation triggered by nonlinear modulation instabilities cannot occur.^{38–41} Indeed, Fig. 1(c) for the simulation results of laser propagation in a homogeneous plasma background at the same average density as the uneven one shows that the laser light maintains its initial Gaussian profile well, with only weak diffraction. In addition, the corresponding scintillation index Σ remains effectively zero, as can be seen in Fig. 1(e). The distinct laser dynamics in homogeneous and uneven plasma backgrounds indicate that the laser branching observed here does indeed originate from the plasma density fluctuations and is thus quite different from laser filamentation. That is, light branching can occur in purely linear systems, such as a continuous-wave laser in a soap film. The term “nonlinear light branching” used here refers to our finding that branched flow patterns become critically dependent on the laser intensity for *I* ≳ 10^{14} W/cm^{2}, as shown in the movie in the supplementary material.

Figure 2 shows the laser light spectrum. As the laser propagates through the uneven plasma, successive random weak scatterings by the density fluctuations cause the laser to branch into other directions. Accordingly, the (*k*_{x}, *k*_{y}) spectrum gradually forms a quarter-arc pattern centered at $|kl|=kx2+ky2\u223c1\u2212(\omega pe/\omega 0)2\omega 0/c\u223c0.96k0$, and the field strength in Fourier space, $E\u0303z$, decreases with increasing laser spread angle defined as *ϕ* = arctan(|*k*_{y}/*k*_{x}|), as can be seen in Fig. 2(a). Here, *ω*_{0} and *k*_{0} are the laser frequency and wavenumber in vacuum, *ω*_{pe} is the plasma frequency, and *c* is the speed of light in vacuum. By contrast, in a homogeneous plasma, the spectrum is simply at *k*_{x} ∼ 0.96*k*_{0} and *k*_{y} ∼ 0, as expected. Figure 2(c) for the time Fourier transform of *E*_{z} shows that the laser frequency maintains its vacuum value well at *ω*_{l} ∼ *ω*_{0} during branching. We can therefore consider that the laser pulse evolves from a coherent Gaussian beam to quasi-monochromatic partially spatially incoherent light by interaction with the uneven plasma. It is of interest to note that Fig. 2(c) also indicates that laser–plasma instabilities (e.g., stimulated Raman/Brillouin scatterings) and the associated plasma heating are negligible for the interaction parameters considered here. Thus, laser filamentation (and spraying) triggered by laser–plasma instabilities is also unlikely to happen here, even though the threshold for these instabilities is much lower than that for self-focusing.^{42,43}

In addition, the robustness of intense laser light branching has been verified by varying the laser polarization, spot size, pulse duration, and profile, as well as the material and unevenness profile of the background medium. Branch formation is clearly observed under the condition that the medium’s correlation length *l*_{c} is larger than the laser wavelength *λ*_{0}, so that the geometric optics limit is satisfied.

## IV. BACKGROUND PLASMA AS A SPACE–TIME-DEPENDENT UNEVEN MEDIUM

A feature of intense-laser branching is the nonlinear response of the background plasma medium, which changes the refractivity and the optical unevenness along the laser paths. In Fig. 3(a), we show the evolution of the potential strength $v0=\eta 2\u0304\u2212\eta eff22/2\eta 2\u0304$, or unevenness, of the plasma by the action of a laser with intensity ranging from 10^{14} to 10^{20} W/cm^{2}, as obtained from PIC simulations. Here, $\eta 2\u0304=\eta eff2$ is the mean-square value of the plasma refractivity *η*_{eff}. For laser intensities below the relativistic threshold (*I* ≲ 10^{17} W/cm^{2}), *v*_{0} reaches a plateau after a sharp increase by photoionization. However, in the relativistic regime (*I* ≳ 10^{18} W/cm^{2}), *v*_{0} after peaking decreases continuously to even smaller values than the initial one. This decrease is due to plasma homogenization by the laser interaction, on a time scale estimated to be *τ* = *l*_{c}/2*c*_{s}, where $cs=Z2/miTe/Z$ is the ion acoustic speed, *Z* and *m*_{i} are the charge number and rest mass of each ion species, *T*_{e} ∼ (*γ* − 1)*m*_{e}*c*^{2} is the bulk electron temperature,^{44} $\gamma =1+a02$ for circular polarization, *a*_{0} = *e*|** E**|/

*m*

_{e}

*ω*

_{l}

*c*is the normalized laser electric field

**, and**

*E**m*

_{e}and −

*e*are the electron rest mass and charge. One can see from Fig. 3(b) that for nonrelativistic lasers,

*τ*is much larger than the pulse duration. Therefore, the plasma unevenness

*v*

_{0}can be considered as quasistatic after the rapid increase by photoionization, consistent with the plateaus shown in Fig. 3(a). By contrast, for relativistic lasers,

*τ*becomes comparable to or smaller than the pulse duration. In this case, removal of the density unevenness by the laser interaction becomes dominant. The background plasma loses its original uneven quality along the laser path, resulting in the pronounced decrease in

*v*

_{0}.

## V. BRANCHING ENHANCEMENT BY PHOTOIONIZATION

For nonrelativistic lasers, since plasma homogenization, as well as laser frequency shift, during the laser interaction can be neglected and the refractivity of the uneven plasma can be considered as slowly varying (see the supplementary material for a detailed discussion), the propagation of the laser can be described by the Helmholtz equation^{17}

where $\eta eff=1\u2212ne(I)/\gamma nc$ is the effective refractivity, and $\eta 2\u0304=\eta eff2$, with the average taken over the laser spot area. Note that the electron number density *n*_{e} is an explicit function of the instantaneous laser intensity *I*, since electrons are produced by photoionization. In this case, one obtains $v0=\delta (I)2\u0304/2[nc\u2212ne(I)/\gamma ]$, where $\delta (I)=ne(I)/\gamma \u2212ne(I)/\gamma $ is the local fluctuating strength of the electron density, and *γ* ∼ 1 for nonrelativistic lasers. We see that since *v*_{0} includes the effect of photoionization, it increases with *I* owing to the increased ionization rate, in agreement with the simulation results in Fig. 3(a). Since laser branching is directly related to *v*_{0}, the branching is also enhanced by photoionization.

Comparison with simulations for laser intensity at 10^{16} W/cm^{2} where photoionization is switched on/off further confirms the above analysis. As shown in Fig. 4(a), a considerable increase in the electron density *n*_{e} along the laser paths is observed if photoionization is included. At a laser intensity *I* = 1 × 10^{16} W/cm^{2}, tunneling ionization of both the Si and O ions to the +4 state can occur (see the supplementary material for a detailed discussion of photoionization).^{27} For our simulation parameters, the average *n*_{e} around the laser axis is ∼0.24*n*_{c}, which is three times the initial value. In addition, since plasma homogenization can be ignored in the nonrelativistic regime, the electrons have the initial disordered distribution of the ions (see the supplementary material for electron and ion dynamics). Therefore, the ACF of the electron density at *t* = 1095 fs remains almost the same as the initial one. Likewise, the local potential strength *v*_{0}(*x*) retains its initial correlations (indicated by the vertical streaks) after the rapid increase due to photoionization, as shown in Fig. 4(c). Note that the correlation length *l*_{c} of 4.8 *µ*m is much smaller than the laser spot size, i.e., the effect of the Gaussian intensity profile across the beam on *l*_{c} can be ignored. The increase in *v*_{0} leads to stronger (but still relatively weak) scattering of the laser, resulting in enhancement of the branching pattern, as shown in Fig. 4(b).

To further characterize the branching, we consider the angular dependence of the laser electric field in Fourier space, defined by $E\u0303z(\varphi )=\u222be\u2212iklx\u2061cos\u2061\varphi e\u2212ikly\u2061sin\u2061\varphi Ez(\varphi )dxdy$. Here, *E*_{z} is used instead of *E*_{y}, to exclude self-generated fields. As shown in Fig. 4(d), the dependence of $E\u0303z(\varphi )$ on *ϕ* at *t* = 35.3 fs is quite similar for both cases, indicating that most of the laser energy still flows in the *x* direction. However, at *t* = 1095 fs, as a result of many successive weak scatterings in the uneven plasma, a large amount of laser light is branched into other directions. The spread angle Θ of the light branches after the first caustics can be defined as that when $E\u0303z$ drops to one-quarter of its maximum. We find Θ ∼ 2*π*/9 when photoionization is included, which is about two times larger than that without photoionization. This result further demonstrates that photoionization enhances the unevenness of the refractivity and thus the branching.

## VI. BRANCHING SUPPRESSION BY RELATIVISTIC EFFECTS

For a relativistic laser with *I* > 1.37 × 10^{18} W/cm^{2}, most of the electrons on the outer shells of the ions are freed, and they can be accelerated to light speed by the laser fields within a single cycle. In this case, further ionization becomes marginal and relativistic laser–plasma interaction effects become significant. The plasma homogenization time *τ* ∼ 217 fs becomes much smaller than the pulse duration. The local refractivity along the laser path now changes simultaneously as the laser propagates, and Eq. (1) becomes inapplicable. In fact, Figs. 5(a) and 5(b) show that the unevenness in the initial electron density distribution vanishes right behind the laser pulse front. Rapid plasma homogenization leads to a decrease in *l*_{c}, and electron resonance in the laser fields causes longitudinal modification of the density distribution and the ACF, as can be seen in Fig. 5(b). In addition, the strong laser ponderomotive force expels the affected electrons, resulting in the formation of a plasma channel behind the laser front^{35–38} and further reduction of the density unevenness, as shown in Fig. 5(c). Figure 5(e) shows that the corresponding local potential strength *v*_{0}(*x*) decreases to much less than the initial one after the rapid increase caused by photoionization. The initial correlation of the unevenness also vanishes. As shown in Fig. 5(d), branching of the laser is suppressed, and its spread angle Θ remains small at 2*π*/67.

## VII. SCALING WITH LASER INTENSITIES

Figure 6(a) for the spread angle Θ of the laser branches for different initial laser intensities shows that Θ increases with *I*_{0} until *I*_{0} ≲ 10^{17} W/cm^{2}, and then it decreases as *I*_{0} increases further. This is in good agreement with the dependence of the potential strength *v*_{0} on the laser intensity at *t* = 1095 fs shown in Fig. 3(a). At moderate intensities, as a result of strong laser branching, the spread angle Θ in an uneven plasma is considerably larger than that in a homogeneous plasma, even though the latter also increases with *I*_{0} owing to ionization-induced self-defocusing.^{45} For relativistic lasers, since the plasma density unevenness is greatly reduced by laser interaction, one can expect that the laser dynamics, as well as the spread angle Θ, become rather similar in uneven and homogeneous plasmas in this regime. Such large spread angles for moderately intense lasers in uneven plasmas can be considered as evidence of nonlinear laser branching in experiments.

Another parameter for characterizing flow branching is the distance *d*_{0} from the boundary (where the flow enters) of the uneven medium to the first branching point. In the linear case, the flow has no influence on the medium, and a universal scaling law for *d*_{0} is $d0\u221dlcv0\u22122/3$. For nonrelativistic picosecond lasers, where plasma homogenization can be ignored and $v0=\delta (I)2\u0304/2[nc\u2212ne(I)/\gamma ]$, one can obtain

Note that the laser intensity appears in the scaling, which is due to the nonlinear effect of photoionization. We see that *d*_{0} decreases with increasing strength of density fluctuations given by $\delta (I)2$, and increases with $nc\u2212ne(I)/\gamma $, indicating that the higher the effective plasma density *n*_{e}/*γ*, the earlier branching occurs. The quasilinear scaling agrees well with the branched flow of nonrelativistic lasers, as shown in Fig. 6(b). However, for relativistic lasers with *I* > 10^{18} W/cm^{2}, laser branching becomes suppressed owing to plasma homogenization. The first caustics, instead of being a branching point, now mark the location where self-focusing starts. Figure 6(b) shows that *d*_{0} now increases with the laser intensity. It is of interest to note that Eq. (2) still agrees fairly well with the simulation results, even though Eq. (1) no longer holds in this regime.

Figure 7 shows that laser branching is also robustly observed in 3D. Here, the incident laser intensity is 10^{16} W/cm^{2}, and the interaction parameters are kept the same as those in 2D, except for the dimension. One can see that a pronounced angle-independent branching pattern forms, which resembles that in Fig. 4(b). At the branching point, the intensity pattern shows intense caustic surfaces separated by bubble-like voids.

We have found that the scintillation index Σ and spread angle Θ that describe the branch pattern can be applied successfully in 3D. As shown in Fig. 6(a), the dependence of the laser spread angle Θ on intensity shows a similar tendency to that in 2D. This confirms that our analysis and the nonlinear nature of intense laser branching are robust in 3D systems, even though the absolute values of Θ are reduced compared with those in 2D, owing to considerably greater number of degrees of freedom involved in a 3D uneven medium. (More detailed discussions of dimensional effects can be found in the supplementary material.)

## VIII. CONCLUSION

In conclusion, we have shown that an intense laser propagating through an uneven plasma can form complex light branches in a nonlinear manner. The nonlinearity originates from the space–time dependence of the plasma refractivity resulting from the laser action. In particular, photoionization can increase the unevenness in the density and thus enhance branch formation. However, relativistic effects of too-intense lasers can suppress branch formation by smoothing the plasma unevenness. These regimes can potentially be verified by experiments based on laser interaction with pre-ionized low-density fibrous or foamy materials, or gas clusters. Our work extends existing studies of optical branching to the nonlinear regime and can thus induce investigations of nonlinear branching in other areas, such as branched wind and ocean waves (sources of freak waves) and fluid flow branching in frangible porous media. Furthermore, branched flow of intense laser light should also open a new area for laser–matter interaction and should be of interest in optical communications, nonlinear optics, and strong-field physics, as well as in laser interactions with foam or turbulent plasmas.

## SUPPLEMENTARY MATERIAL

See the supplementary material for discussions of PIC simulations incorporating photoionization models, electron and ion dynamics, the validity of using the Helmholtz equation, and dimensional effects. The supplementary movie shows the results of simulations of intense laser light (at intensities from 10^{14} to 10^{20} W/cm^{2}) propagation in space–time-dependent media.

## ACKNOWLEDGMENTS

This work is supported by the National Natural Science Foundation of China (Grant Nos. 12205201, 12175154, 11875092, and 12005149), the Natural Science Foundation of Top Talent of SZTU (Grant Nos. 2019010801001 and 2019020801001), and GCS Jülich (Project No. QED20) in Germany. The EPOCH code is used under a UK EPSRC contract (Grant Nos. EP/G055165/1 and EP/G056803/1). K.J. would like to thank Q. Wang, X. Luo, X. F. Shen, H. Peng, and T. Y. Long for useful discussions.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

K.J. and T.W.H. developed the theoretical work. K.J. conducted the simulations. K.J., T.W.H., and C.N.W. analyzed the data and produced the figures. H.Z., S.Z.W., H.B.Z., and A.P. helped review and interpret the data. K.J., T.W.H., and M.Y.Y. wrote the article. T.W.H., C.T.Z., and S.C.R. supervised the work. All authors have reviewed, discussed, and agreed to the complete article.

**K. Jiang**: Conceptualization (equal); Data curation (lead); Formal analysis (equal); Investigation (lead); Methodology (lead); Validation (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **T. W. Huang**: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Project administration (equal); Supervision (equal); Validation (equal); Writing – original draft (supporting); Writing – review & editing (equal). **C. N. Wu**: Formal analysis (supporting); Visualization (supporting). **M. Y. Yu**: Investigation (supporting); Writing – original draft (supporting); Writing – review & editing (equal). **H. Zhang**: Data curation (supporting); Formal analysis (supporting); Resources (equal); Software (equal); Validation (supporting). **S. Z. Wu**: Data curation (supporting); Resources (supporting); Software (supporting); Validation (supporting). **H. B. Zhuo**: Formal analysis (supporting); Investigation (supporting); Validation (supporting). **A. Pukhov**: Formal analysis (supporting); Investigation (supporting); Resources (supporting); Supervision (supporting); Validation (supporting). **C. T. Zhou**: Data curation (supporting); Formal analysis (supporting); Funding acquisition (equal); Investigation (supporting); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (supporting). **S. C. Ruan**: Resources (supporting); Supervision (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

### APPENDIX: SIMULATION DETAILS

Our PIC simulations are conducted using the epoch code.^{29} For the simulations of the Gaussian laser pulse, “*simple_laser*” and “*simple_outflow*” longitudinal boundaries and open lateral boundaries are used. Periodic lateral boundaries are used for the reference simulations of plane waves. Unless otherwise stated, “*field_ionization*” is turned on, which accounts for different modes of photoionization, including multiphoton ionization, tunneling ionization, and barrier-suppression ionization. In 2D (3D) simulations, the simulation box is −5 *µ*m < *x* < 215 *µ*m, −55 *µ*m < *y* (and *z* in 3D) < 55 *µ*m, with 2200 × 1100 (1100 × 550 × 550) grid cells and 30 (2) macroparticles per cell for each species. Note that more macroelectrons will be produced by photoionization. For the interacting medium, we choose a weakly pre-ionized SiO_{2} plasma (Si^{2+} and O^{+}) with uneven density distribution, which is located in 0 *µ*m < *x* < 215 *µ*m, −55 *µ*m < *y* (and *z* in 3D) < 55 *µ*m. The average densities of Si^{2+}, O^{+}, and electrons are 0.02*n*_{c}, 0.04*n*_{c}, and 0.08*n*_{c}, respectively, where *n*_{c} ∼ 1.1 × 10^{21} $cm\u22123/\lambda 02$ is the critical density. The uneven density distribution is generated by assigning random points at a minimum separation of 2 *µ*m within the plasma area with values taken from a set satisfying a Gaussian random distribution between zero and three times the average densities. The full plasma density map with *C*^{2} continuity is obtained by applying biharmonic spline interpolation based on these random points. The overall average densities and a coefficient of density variation of 30% are used as constraint conditions. The obtained density unevenness is of an isotropic correlation length *l*_{c} = 4.8 *µ*m, defined as the shortest distance in which the value of the corresponding autocorrelation function (ACF) drops to 10% of the zero-shift value.^{28} The ACF is calculated as $ACF(i,j)=\u2211m=0M\u22121\u2211n=0N\u22121ne(m,n)ne(m+i,n+j)$, where *n*_{e} is the initial electron density, *M* and *N* are the numbers of grid cells occupied by the electrons in the longitudinal and lateral directions, respectively, and 0 ≤ *i* < 2*M* − 1 and 0 ≤ *j* < 2 *N* − 1. The input laser pulse is simulated as a Gaussian beam with circular polarization incident normally from the left boundary. Its central wavelength *λ*_{0} = 1.06 *µ*m, peak intensity *I*_{0} = 10^{14}–10^{20} W/cm^{2}, and FWHM spot size *r*_{0} = 16.65 *µ*m. The laser pulse duration is 1.5 ps, with a flat-top time profile for simplicity. The reference simulations of plane waves use identical parameters.