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₀ from the source to the first branching point follows the scaling $d_0\propto l_c{v}_0^{-2/3}$,1,15,21–24 where l_c is the correlation length of the medium’s unevenness, and v₀ is its strength parameter (to be defined later).

With the development of chirped pulse amplification technology,25 lasers with intensity I ranging from 10¹⁴ to 10²⁰ W/cm² 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¹⁴ W/cm², 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¹⁸ W/cm², 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¹⁴–10¹⁷ W/cm²) 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.

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₂ plasma (with Si²⁺ 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²⁺, O⁺, and electrons are 0.02n_c, 0.04n_c, and 0.08n_c, respectively, where n_c ∼ 1.1 × 10²¹${\text{cm}}^{-3}/{\lambda }_0^2$ is the critical density. The coefficient of density variation is ∼30%. An intense circularly polarized Gaussian laser pulse of central wavelength λ₀ = 1.06 µm, peak intensity I₀ = 10¹⁴ W/cm², and FWHM spot size r₀ = 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 simulations29 (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 clusters32 and in backward Raman amplification schemes for intense-laser pulse compression.33

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²⟩/⟨I⟩²) − 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.2n_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¹⁴ W/cm², 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 $|k_l|=\sqrt{k_x^2+k_y^2}\sim \sqrt{1-{({\omega }_{pe}/{\omega }_0)}^2}{\omega }_0/c\sim 0.96k_0$, and the field strength in Fourier space, ${\tilde{E}}_z$, decreases with increasing laser spread angle defined as ϕ = arctan(|k_y/k_x|), as can be seen in Fig. 2(a). Here, ω₀ and k₀ 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.96k₀ 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 ∼ ω₀ 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 λ₀, so that the geometric optics limit is satisfied.

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 $v_0=\sqrt{⟨{\left(\overline{{\eta }^2}-{\eta }_{\text{eff}}^2\right)}^2⟩}/2\overline{{\eta }^2}$, or unevenness, of the plasma by the action of a laser with intensity ranging from 10¹⁴ to 10²⁰ W/cm², as obtained from PIC simulations. Here, $\overline{{\eta }^2}=⟨{\eta }_{\text{eff}}^2⟩$ is the mean-square value of the plasma refractivity η_eff. For laser intensities below the relativistic threshold (I ≲ 10¹⁷ W/cm²), v₀ reaches a plateau after a sharp increase by photoionization. However, in the relativistic regime (I ≳ 10¹⁸ W/cm²), v₀ 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/2c_s, where $c_s=\sqrt{⟨Z^2/m_i⟩T_e/⟨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_ec² is the bulk electron temperature,44$\gamma =\sqrt{1+a_0^2}$ for circular polarization, a₀ = e|E|/m_eω_lc is the normalized laser electric field E, and 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₀ 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₀.

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 equation17$$-{\nabla }^2\mathit{E}+k_0^2(\overline{{\eta }^2}-{\eta }_{\text{eff}}^2)\mathit{E}=k_0^2\overline{{\eta }^2}\mathit{E},$$(1)where ${\eta }_{\text{eff}}=\sqrt{1-n_e(I)/\gamma n_c}$ is the effective refractivity, and $\overline{{\eta }^2}=⟨{\eta }_{\text{eff}}^2⟩$, 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 $v_0=\sqrt{\overline{\delta {(I)}^2}}/2[n_c-n_e(I)/\gamma ]$, where $\delta (I)=⟨n_e(I)/\gamma ⟩-n_e(I)/\gamma$ is the local fluctuating strength of the electron density, and γ ∼ 1 for nonrelativistic lasers. We see that since v₀ 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₀, the branching is also enhanced by photoionization.

Comparison with simulations for laser intensity at 10¹⁶ W/cm² 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¹⁶ W/cm², 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.24n_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₀(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₀ 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 ${\tilde{E}}_z(\varphi )=\left|\int {\mathrm{e}}^{-\mathrm{i}{k}_{l}x\cos \varphi }{\mathrm{e}}^{-\mathrm{i}{k}_{l}y\sin \varphi }{E}_z(\varphi )\mathrm{d}x\mathrm{d}y\right|$. Here, E_z is used instead of E_y, to exclude self-generated fields. As shown in Fig. 4(d), the dependence of ${\tilde{E}}_z(\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 ${\tilde{E}}_z$ 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.

For a relativistic laser with I > 1.37 × 10¹⁸ W/cm², 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 front35–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₀(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.

Figure 6(a) for the spread angle Θ of the laser branches for different initial laser intensities shows that Θ increases with I₀ until I₀ ≲ 10¹⁷ W/cm², and then it decreases as I₀ increases further. This is in good agreement with the dependence of the potential strength v₀ 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₀ 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₀ 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₀ is $d_0\propto l_c{v}_0^{-2/3}$. For nonrelativistic picosecond lasers, where plasma homogenization can be ignored and $v_0=\sqrt{\overline{\delta {(I)}^2}}/2[n_c-n_e(I)/\gamma ]$, one can obtain$$d_0\propto l_c{⟨\delta {(I)}^2⟩}^{-1/3}{⟨n_c-\frac{n_e(I)}{\gamma }⟩}^{2/3}.$$(2)Note that the laser intensity appears in the scaling, which is due to the nonlinear effect of photoionization. We see that d₀ decreases with increasing strength of density fluctuations given by $⟨\delta {(I)}^2⟩$, and increases with $⟨n_c-n_e(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¹⁸ W/cm², 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₀ 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¹⁶ W/cm², 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.)

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.

Acknowledgment. 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.