Many experiments conducted in multi-kilojoule laser facilities, whether they concern astrophysical phenomena, high-energy-density physics,1 or inertial confinement fusion (ICF),2–4 bring critical insights into matter under extreme conditions. In these facilities, energetic laser beams are used to heat and compress the matter to millions of degrees and mega- to gigabar pressures. The predictability of these experiments requires a precise understanding of laser–plasma interactions, such as the mechanisms responsible for the energy deposition and those responsible for the growth of various instabilities such as stimulated Raman or Brillouin scattering,5–8 cross-beam energy transfer,9 two-plasmon decay,10,11 collective scattering,12–14 and self-focusing.15

Furthermore, optical smoothing techniques combining random-phase plates (RPPs) and smoothing by spectral dispersion (SSD), available at the Laser Méga Joule (LMJ) facility, the National Ignition Facility (NIF), the Laboratory for Laser Energetics (LLE), and SG-III-class lasers, improve the laser intensity profile.16,17 These techniques reduce the spatial coherence (through RPPs) and temporal coherence (through SSD) of the light, resulting in intensity fluctuations on the scale of a few wavelengths and lasting a few picoseconds: the so-called speckles. Polarization smoothing (PS), which consists in splitting the pulse into two uncorrelated superimposed beams with perpendicular polarizations, may also be used. Together, these smoothing techniques will in turn affect the laser–plasma interactions on macroscopic scales, such as the energy deposition region,18,19 the scattering direction of the light wave,20–25 and the amount of expected reflectivity.26–29 Accounting for the hot-spot dynamics thus requires to reconcile the sub-micrometer and sub-picosecond physics of the RPP/SSD beam with the millimeter size and nanosecond duration of the experiments. In this context, a hydrodynamic description of the plasma can be coupled with an approximated Maxwell solver.30–34 However, this formalism fails to capture self-consistently the Landau damping of acoustic waves or other kinetic effects and may turn out to be numerically unsuitable whenever multidimensional effects, solid-density physics, or radiative phenomenons arise.

In such systems, hydrodynamic codes rest on a rough description of light, such as the classical ray tracing scheme,35–38 which, in its native form, only captures the wave refraction and basic energy deposition. Accurate description, by the ray tracing scheme, of simple light properties such as the intensity of the wave and its local spectrum, is still an active area of research.39–42 Although great efforts have been made to improve these schemes by including back or side scattering of the light caused by wave mixing processes,41,43 modeling of the speckle-scale physics remains largely unexplored.44–47

Although the flow-induced deviation may be entangled with other wave mixing processes,25 this study addresses the beam bending of a laser beam48–51 under the assumption that other instabilities, such as filamentation and forward/backward Brillouin/Raman scattering, are negligible. The beam bending occurs when the laser-driven density fluctuations are advected by a flow, resulting, owing to a wave-guide effect, to deflection of the electromagnetic wave from its original axis of propagation, toward the flow direction. As this effect may take place in a perfectly homogeneous plasma, the beam deflection may add up to the well-known refraction of light caused by density gradients.

In Sec. II, we first briefly recall the kinetic and fluid modeling of Ref. 52, which predicts the deflection angle of a Gaussian laser pulse, in both the transient and asymptotic regimes of the plasma density response. The model is then extended in Sec. III to the centroid deviation of a spatially (using RPP) and temporally (with SSD) incoherent laser pulse, with the effect of polarization smoothing by fits coming from three-dimensional (3D) Parax53 simulations. Our model is incorporated into a ray tracing description of the beam propagation and compared successfully with our numerical results in Sec. IV. Then, the laser propagation under realistic indirect drive ICF plasma conditions (NIF Shot N181209 from the Hybrid B campaign54–57) is addressed and evidences a sensible impact of the flow-induced deviation on the light energy path and energy deposition. Our concluding remarks and perspectives are gathered in Sec. V.

International system of units (SI) units are used throughout this paper, the Boltzmann constant is dropped, temperatures are given in eV, and vectors are denoted by bold symbols.

In Ref. 50, it is proved that in the small-angle limit, the transverse averaged beam wavevector ⟨k⊥⟩k⊥ can be related to the plasma electron density fluctuations δn_e through1k0d〈k⊥〉k⊥dx⋅vd|vd|=−12nenc∇⊥δnene⊥,(1)X⊥=∫dr⊥X(r)I(r)∫dr⊥I(r),(2)Xk⊥=∫dk⊥X(k)I(k)∫dk⊥I(k),(3)where we have introduced k₀ = 2π/λ₀, n_c, and v_d, the main laser wavevector, critical density, and flow velocity, respectively. Moreover, r_⊥ and k_⊥ are the position and wavevector transverse to the main laser x direction. The transverse averages in real and Fourier space, ⟨·⟩_⊥ and ⟨⋅⟩k⊥, respectively, are defined through Eqs. (2) and (3). We assume hereinafter a perfectly homogeneous plasma. The dependence of the normalized density fluctuations on time and position, δn_e(t, r)/n_e, will be assumed to follow the linearized hydrodynamic equations. In the plasma rest frame and after a transverse Fourier transform [defined hereinafter as f(ω, k) = ∫f(t, r)e^−·+drdt], the relation between the density fluctuations and the laser intensity I(k) = I₀g(k) satisfies∂2∂t2+2γ0|k|cs∂∂t+cs2k2δnene=A|k|cs2k2I0ncvgTeαkαfg(k)eivd⋅kt.(4)Here, we have introduced the sound speed cs=[(ZiTe+3Ti)/mi]1/2, the electron/ion mass and temperature m_e/i and T_e/i, the normalized acoustic Landau damping rate γ₀ = ν/|k|c_s, the laser group velocity vg=c1−ne/nc, and the light speed in vacuum c. The ratio α_k/α_f, whereαk[M0cos(θ)]=−Z′(ζe)2∑iZ′(ζi)ZiTeTiZineniZ′(ξe)+∑iZ′(ζi)ZiTeTiZineni,(5)ζe/i=−k⊥⋅vd|k⊥|me/i2Te/i,(6)αf[M0cos(θ)]=κ1−M02cos2θ+2iγ0|M0|cosθ,(7)cos(θ)=k⊥⋅vd|k⊥||vd|,(8)κ=ZiTemics2,(9)introduced in Ref. 52, is an ad hoc correction that ensures convergence to the kinetic asymptotic limit for highly Landau-damped plasmas. We have also introduced the Mach number M₀ = v_d/c_s and Z′, the first-order derivative of the plasma dispersion function.58 Finally, the factor A_k on the right-hand side of Eq. (4) accounts for nonlocal thermal effects on the density fluctuation amplitude59,60 and can be written asAk(u)=12+Zi0.074u2+0.88u4/7+2.541+5.5u2Ω,u=|k⊥|λmfpZi,Ω=1if |k⊥⋅vd|/νei<1,0elsewhere,(10)where λ_mfp is the electron mean free path and where A_|k| = 1/2 in the collisionless limit. As mentioned in Ref. 61, this correction can be used provided the acoustic wave frequency |ω_s| = |kv_d| remains smaller than the electron–ion collisional frequency ν_ei (momentum exchange). This condition is contained in the factor Ω, which ensures that A_k recovers its collisionless value when the condition is not fulfilled. Although easily introduced in our theoretical estimates, this factor remains far more challenging to implement in a self-consistent numerical modeling where other laser–plasma effects are accounted for, in which case Ω is set to unity (see Sec. II B). Likewise, the Landau damping factor γ₀ in Eqs. (7) and (12) may be expressed in the collisionless limit (Landau) [Eq. (10) in Ref. 52] as in Sec. III. When collisional effects are not negligible, the fit introduced in Ref. 62 will be used, leading to a dependence of γ₀ on |k|, as in Sec. IV B.

In the laboratory frame, the solution of Eq. (4) for γ₀ < 1, with the initial conditions δn_e(t = 0) = ∂_tδn_e(t = 0) = 0, readsδne(k,t)ne=αkA|k|I0g(k)ncvgTef(k,t),(11)f(k,t)=1+a+eg+cs|k|t−a−eg−cs|k|t,(12)a±=−iM0cos(θ)−γ0∓i1−γ022i1−γ02,(13)g±=−γ0±i1−γ02−iM0cos(θ).(14)Here, M₀ = |M₀|, where M₀ and v_d are assumed to be normal to the main laser x axis. The combination of the above density fluctuations, written in Fourier space, with Eq. (1) implies an inverse Fourier transform and leads todθdx=1k0d〈k⊥〉⊥dx⋅vd|vd|=−12nencI0ncvgTe∫d2k(2π)2ik⋅vd|vd|αkA|k|g(k)f(k,t)〈eik⋅r⊥〉⊥.(15)

Hence, the estimation of ⟨eik⋅r⊥⟩⊥ with Eq. (2) leads to the final beam bending rate, and this depends on the Fourier transform of the transverse intensity profile, g(k). References 48 and 50 assume Gaussian beams, while Ref. 63 addresses RPP beams, in the asymptotic regime only. To avoid any assumption about the SSD beam temporal spectrum, as in Refs. 64 and 65, we will simply relate the averaged SSD beam centroid deflection to the speckle contributions using Gaussian transient regime estimates.52

Hence, using g(k) = exp(−k²σ²/8) with σ = λ₀f_♯ and assuming a y-aligned drift velocity, the deflection rate may be recast asdθdx=−nencI02vgncTe1σG(t),(16)withG=σ38πI∫−∞∞dkk2Ake−k2σ2/4×∫0πdθαk/f[M0cos(θ)]cos(θ)f(k,t).(17)

Figure 1(a) illustrates the temporal evolution of a single 3D Gaussian speckle of f-number f_♯ = 8 and wavelength λ₀ = 0.35 μm in a fully ionized carbon plasma as predicted by Eqs. (16) and (17). As expected,52 the deflection is larger for a flow with M₀ = 1 (red lines) in the asymptotic limit t > 100 ps (where dθ/dx ≃ 0.2 rad/mm, red solid line) than for M₀ = 0.8 or 1.2 (where dθ/dx ≃ 0.17 rad/mm and ∼0.02, blue and black solid lines, respectively). In the transient, however, the bending at resonance remains slightly weaker than for the case M₀ = 0.8 (black lines) or 1.2 (blue lines, t < 10 ps). As a corollary, the resonance of the hot spot bending rate around a Mach number of the order of unity is not noticeable during the transient.

The propagation of a spatially and temporally smoothed laser beam in a flowing plasma, as shown in Ref. 44, exhibits a deflection of its centroid that can be related to the contribution of hot spots. Assuming that the speckles are independent and of Gaussian form, we may now relate the time-averaged beam centroid deflection rate to the contribution of the speckles (subscript s) throughdθSSDdx=∑sdθsdxt,(18)where ⟨X⟩_t is the average of X over a 2π/ω_m period (where ω_m is the main modulation frequency of the SSD device). The speckle bending rate dθ_s/dx stems from Eq. (16). Owing to diffraction, σ and I₀ become σ[1+(x−xs)2/zc2]1/2 and Is/[1+(x−xs)2/zc2]3/2, respectively. The speckle reaches its maximum intensity I_s at the position x_s and has a Rayleigh length of zc=πf♯2λ0, leading todθSSDdx=−Sn0ncI02vgncTe⟨G⟩τSSDσ,(19)where I₀ now denotes the beam-averaged intensity and whereS=∑sIs/I0[1+(x−xs)2/zc2]3/2.(20)The factor S introduced in Eqs. (19) and (20) incorporates the speckle-dependent factors such as the speckle intensities, spatial profiles, and positions. As S takes account of the complex RPP-induced speckle dynamics, we consider it as a parameter to be fitted from the paraxial simulations. This procedure further smooths out the hot spot contributions over the whole pulse. This approximation neglects the speckle disparities resulting from the variability of their intensity. Furthermore, we neglect any significant modification of the speckle distribution during the beam propagation, since this could affect the value of S. To mimic the effect of the SSD on the beam bending rate, the speckle pattern is assumed to last a time τ_SSD, proportional to the coherence time T_SSD. The latter is related to the modulation depth Δ (defined hereinafter at 3ω) and to the modulation frequency ω_m of the SSD device through T_SSD = 2π/ω_m(2Δ + 1). On introduction of the factor β, the parameter τ_SSD, referred to as the speckle coherence time, can be written asτSSD=β2πωm(2Δ+1)=βTSSD.(21)This second factor, of the order of unity, represents the efficiency of SSD in beam bending and will be deduced from the paraxial simulations.

We introduce the temporal mean of a function X(t) over an SSD period as⟨X⟩τSSD=1τSSD∫0τSSDX(t)dt,(22)Equation (17) leads to⟨G⟩τSSD≃σ38πI∫−∞∞dkk2Ake−k2σ2/4×∫0πdθαk/f[M0cos(θ)]cos(θ)×1+a+eg+cs|k|τSSD−1g+cs|k|τSSD−a−eg−cs|k|τSSD−1g−cs|k|τSSD.(23)For the sake of simplicity, we propose in Appendix D a fully analytical fit of the above function, obtained in the fluid framework and valid in the regime 0.007 ≤ γ₀ ≤ 0.05 and |M₀| < 3 for a monospecies plasma, illustrated in Figs. 7(a)–7(d).

The mean deflection angle is related to the set of Eqs. (18), (19), and (23). The dependences of the above predictions on the Mach number are illustrated in Fig. 1(b) for He²⁺ (blue lines), C⁶⁺ (red lines), and Au⁵⁰⁺ (black lines) 10% critical plasmas (with T_e = 2.5 keV and T_i = 1 keV) and for a speckle coherence time τ_SSD = 7 ps. The deflection rate is an increasing function of the ionization number, owing to the thermal correction of Eq. (10). By virtue of Sec. II B, the deflection rate is not peaked around M₀ = 1, owing to the short coherence time compared with the transient. A higher density causes the mean free path to decrease, thus increasing the local energy deposition and thereby enhancing the hole boring. The regime of validity of Eq. (10) (ω = kv_d < ν_ei), incorporated in the factor Ω, may be assessed by estimating the critical Mach number M_c = ν_ei/[c_sk₀/(2f_♯)] above which the thermal correction is no longer valid. When Z_i = 2 (we obtain M_c ∼ 0.7), the factor A_k is ∼20% above its collisionless value, indicating a poor impact of thermal effects whatever the value of Ω. However, the carbon deflection rates exhibit slight inflections located around M₀ = M_c ≃ 1.7, due directly to the condition ω = kv_d < ν_ei. Additionally, the large value of M_c ≃ 25 obtained for the Au⁵⁰⁺ case demonstrates the validity of the thermal corrections of Eq. (10) in this material.

Because of the very large values of Z_iT_e/T_i in the gold case (Z_iT_e/T_i = 125), the kinetic and fluid 3D deflection rates (black and red lines respectively) are in agreement.52 For the same reason, accounting for an accurate distribution of gold ionization levels under non-local-thermodynamic-equilibrium (NLTE) conditions as predicted by DEDALE66 [in the multiple-ion-species formulation of Eq. (5)] leads to only a small reduction (by ∼10%) of the final bending rates compared with the single-fluid average ionization number model as used in this study.

Finally, we performed five 3D hydrodynamic simulations with a paraxial solver (Hera32,67) of a Gaussian beam propagating in a C⁶⁺ 10%-critical plasma with the laser and plasma parameters of Fig. 1(b) and for various Mach numbers. The simulation parameters and the method of measuring the deflection rates are detailed in Appendix A. The thermal corrections of Eq. (10) (with Ω set to 1) have been accounted for. The average over the first 7 ps of the beam centroid deflection angles [⟨⟨dθ/dx⟩⊥⟩t, Eqs. (2) and (22)] at the exit plane of the simulations are illustrated as red circles in Fig. 1(b) and agree very well with the corresponding theoretical prediction (obtained with S = β = 1) represented by the dashed red line. This confirms the validity of Eqs. (19) and (23) for a Gaussian beam. In the following, the time averaged centroid deflection rate of a SSD beam will be related to the above Gaussian speckle predictions by fitting the two scalar variables S and β using dedicated 3D paraxial numerical results.

We aim at finalizing and validating our speckle-scale SSD beam bending model in light of paraxial 3D Parax53 hydrodynamic simulations. We will thus use in this section the fluid plasma response in our model [Eq. (7)]. More realistic predictions, such as in Sec. IV B, will then be made using the kinetic deflection rates. Ideally, the kinetic flow-induced deflection should be fitted with kinetic codes (such as particle-in-cell codes)—a kind of simulation that is still out of reach with present supercomputers. In the Parax code, the linear Landau damping operator is computed in Fourier space,50,68,69 and the light is propagated through a paraxial solver. Parax uses a linearized hydrodynamics module that is applied transversely to the main laser direction and that will moderate the subsequent required numerical load. This approach is consistent with the fact that the density perturbation amplitudes considered in our model are always weak. For the sake of simplicity and for the purposes of comparison, the bremsstrahlung energy deposition is also neglected. Note that as the pulse deflection is reduced by its loss of temporal coherence, large simulation domains have to be considered in order to obtain clean quantitative comparisons. The thermal corrections of Eq. (10) are accounted for in our 3D simulations.

Parax has specific field injections and diagnostics to correctly simulate and characterize the different smoothing techniques. To reach more realistic ICF conditions than those in Ref. 52, the choice has been made to use a laser wavelength of λ₀ = 0.35 μm. The focal spot is located at the center of the simulation domain, x_foc = 1 mm, with a f-number f_♯ = 8.88. The spectral dispersion has a modulation frequency ω_m = 2π × 14.25 GHz. The laser beam with I₀ = 2 × 10¹⁴ W/cm² propagates through a 2 mm flowing plasma with v_d/c_s = 0.9 and composed of helium with n_e = 0.1n_c, T_e = 2 keV, and T_i = 500 eV while accounting for the thermal corrections of Eq. (10). The simulation details are given in Appendix B. Different smoothing configurations have been simulated.

Owing to the phase modulation, the center of the beam oscillates as soon as it has been injected, with a period corresponding to the modulation frequency, 2π/ω_m ≃ 70 ps. This effect leads to periodic motion of the centroid at the simulation exit plane, as shown in Fig. 2. Note that the SSD beam bending model that we propose only predicts the temporal average of the beam centroid deflection, thus neglecting the motion of the spatial envelope.

The dependence of the beam’s deflection angle on the f-number f_♯ is illustrated in Fig. 3(a). Note that we varied the f-number, keeping the mean intensity I₀ constant. Within the ansatz of Gaussian beam deflection, the latter should be proportional to f_♯ [see the integration of a Gaussian beam bending deflection rate over a Rayleigh length, Eq. (31) of Ref. 52]. However, under the assumption that S can be replaced by a constant, Eq. (19) predicts θSSD∝G(τSSDcs/f♯λ0)/f♯=o(1/f♯), as illustrated by the solid line in Fig. 3(a). Furthermore, although the speckle-scale contribution to the deflection increases with f_♯, the inter-speckle distance and the transient regime amplitude both depend on f_♯. Hence, smoothed over the whole SSD beam, the resulting centroid deviation decreases faster than 1/f_♯ [see the dashed line in Fig. 3(a)]. Figure 3(a) also indicates that if the beam’s effective f-number f_♯ were to decrease, as is to be expected when plasma smoothing occurs,70 the resulting flow-induced deviation should be exacerbated.

We may now fit our model with the most relevant smoothing techniques in the light of the paraxial simulations performed with the Parax code.53 In the case of 1D transverse smoothing by spectral dispersion (TSSD) oriented normally to the flow direction, our numerical results, represented by the red triangles in Fig. 3(b), confirm that the increase in SSD modulation depth (i.e., the decrease in T_SSD) reduces the amount of beam bending. Moreover, using S = 1.7 and β = 2 (solid lines) correctly reproduces the centroid deviations. Equation (45) of Ref. 65, represented by the red dashed line, predicts the deflection of a beam assuming a Gaussian frequency spectrum. Its spectral width is thus narrower than for a SSD laser beam. Consequently, the dependence of beam bending on the speckle coherence time is too steep and the predictions overestimate our numerical results for τ_SSD ≳ 5 ps.

For the case in which the TSSD direction is parallel to the flow, Ref. 65 has found a significantly smaller deflection rate. The corresponding Parax simulations [crosses in Fig. 3(b)] exhibit a deflection rate roughly five times smaller than in the perpendicular configuration. In this case, the theoretical predictions with S = 1.7 and β = 1 give good agreement (dot-dashed line). Those two cases confirm that for a beam smoothed by 1D transverse SSD, such as on the NIF, different deviations are to be expected, depending on the direction of the flow relative to the smoothing direction. Indeed, when the direction of the spectral dispersion is perpendicular to the flow, the beam bending rate depends exclusively on the value of the flow. By contrast, when the direction of spectral dispersion is parallel to the flow, the speckles are moving along the flow with supersonic velocities. The relevant parameter to assess beam bending is the difference between the flow velocity and the transverse velocity of the speckles, Δv. In TSSD, the time evolution of the speckle transverse velocity is roughly sinusoidal.71 For example, for Δ = 18, the maximum transverse velocity reaches v_y/c ∼ 0.004, which corresponds to M ∼ 3.3. Δv then evolves rather quickly as a function of time, and is most of the time far from c_s, where the maximum of the deflection rate is located. This explains both the small time-averaged value of the beam bending and the important centroid oscillations with time that can be seen in Fig. 2, in which we represent the temporal evolution of the beam centroid for Δ = 5.1 in the four simulated configurations for beam smoothing.

Additionally, the impact of PS has also been addressed for flows oriented normally to the transverse SSD direction. In our simulations, the polarizations in the near field are in the diagonal configuration, as for the inner cones on NIF.72 Illustrated by black triangles in Fig. 3(b), our 3D Parax simulations agree fairly well with our model for (S, β) = (0.8, 3). They also suggest a weak influence of the polarization smoothing on the bending rate when the coherence time is small (T_SSD < 4 ps). For longer coherence time (Δ = 5.1), beam bending appears to be reduced when PS is used in the simulation. This reduction can be attributed either to the uncorrelated spatial speckle distribution between the two polarizations or to another nonlinear effect. Indeed, when PS is not used, beam spreading is noticeable in the simulation (the angular aperture is increased by ∼20%, and the effective f_♯ decreases from 9 to 7.5). The increased beam bending in the absence of PS could be explained by the decrease in f_♯, as shown in Fig. 3(a). Given the weak influence of PS, the latter will be neglected when the flow is parallel to the transverse SSD direction.

The dependence of the bending rate on the angle between the flow component normal to the main laser axis and the SSD direction, ϕ, is outside the scope of this study. However, the value of β evolves from 1 to 2 when ϕ = 0 (TSSD ‖v_d) and ϕ = π/2 (TSSD ⊥v_d), respectively. The formula β(ϕ) = 1 + |sin(ϕ)| could be well suited to generalized the model to 3D.

Finally, another set of Parax simulations were performed with longitudinal temporal smoothing (LSSD), as used on the LMJ facility. In this configuration, the phase shift imposed between the different frequencies is radial, causing the speckles to move mainly longitudinally. The results are represented by the red squares in Fig. 3(b). The deflection rates are bounded by the two cases SSD ⊥v_d and SSD ‖v_d, with the advantage of presenting a value independent of the flow direction. The use in our model of (S, β) = (6, 0.8) leads to fair agreement.

LSSD is more efficient than TSSD for small coherence times only. The LMJ optimum modulation depth is close to 15 (defined for λ₀ = 0.35 μm) for ω_m = 2π × 14.25 GHz, leading to τ_SSD ≃ 2.26 ps, and so for the parameters of Fig. 3(b), the overall LMJ deflection over 2 mm is ∼1 µm. The NIF currently uses τ_SSD ≃ 6.5 ps (ω_m = 2π × 17 GHz and Δ ∼ 4) with PS and thus corresponds to a larger deflection than that for the LMJ configuration (∼10 and ∼3 µm for flows perpendicular and parallel to the TSSD direction, respectively). As a result, the beam bending mitigation in high-energy laser experiments is efficiently done by increasing the laser’s temporal spectral width. However, a large spectral width, apart from the associated FM–amplitude modulation (AM) conversion effect on laser performance,73,74 may lead to anisotropic deflection rates when TSSD is used.

Before including the present theoretical calculations in a ray tracing scheme, we gather in Table I the values of the fitting parameters as obtained in the different geometries and SSD configurations.

The time-averaged centroid deviation rate may naturally be implemented in a ray tracing scheme as a correction to the well-known refraction effect arising from the eikonal equation, by deflecting each ray accordingly to Eq. (19). The ray direction and position k and r, respectively therefore satisfydkds=−k0∇ne2ncη−k0Sn0ncI02vgncTe⟨G⟩τSSDσv⊥|v⊥|,drds=kvgω0,v⊥=vd−vd⋅k|k|2k,(24)where s is the ray abscissa and η the optical index.

The above model has been implemented in a new ray tracing module of the hydrodynamic code Hera.32,67 For solving the ray trajectory, we split the rectangular meshes of the 2D Hera hydrodynamic module into four triangles using their diagonals. In each triangle, the density profile is assumed to be linear, allowing us to have a continuous density profile throughout the simulation domain. Hence, solving the eikonal equation in each triangle allows us to split the ray trajectory into a sequence of parabolas, ideal for fast computing. Summing the contribution of each ray that passes through a given mesh following Ref. 43, we obtain the local averaged intensity I₀ used in Eq. (24). To circumvent the nonlinearity, we propose to use the local intensity calculated at the previous hydrodynamic time-step.

To mimic the averaged intensity profile at and off best focus, we define a lens position, hereinafter at x_lens = −10 m from which rays are shot. The rays’ departure y positions on the lens are chosen randomly, and their initial direction points toward a randomly distributed position at the focal plane (at x = x_foc), following the measure given by the chosen intensity profile at best focus. Finally, the rays are propagated in vacuum from the lens to the left boundary condition of the simulation (at x = x_BC). To avoid statistical bias, the positions and directions of the rays are shot at the beginning of each time steps. Note that this procedure, as used in Ref. 38, prevents nonphysical thermal filamentation arising from systematic statistical errors induced by the fluctuations in laser ray number per cell.

Up to now, we have assumed v_d to be perpendicular to the main laser axis. For any flows, the axis along which the beam bending contribution in Eq. (24) lies results from the projection v_⊥ of the fluid velocity v_d on the plane perpendicular to the beam main direction. Importantly, this modification to the classical ray tracing scheme [Eq. (24)] remains minimal, since it does not modify the main algorithm for ray propagation (described above), apart from adding a dependence on the beam intensity of the rays’ trajectories. Equation (24) may be applied in a 2D ray tracing simulation performed in the (k₀, v_d) plane (where k₀ is the main laser axis) with the 3D beam bending deflection rates [and thus using ⟨G⟩τSSD, Eq. (23), instead of its 2D equivalent]. Referred to subsequently as a 2.5D geometry, this will allow us to estimate the bending rate in a realistic geometry in a cheap 2D simulation.

To validate the implementation of the beam bending model in the ray tracing module of Hera, we performed four 2.5D ray tracing Hera simulations. A mesh size of dx = 7.8 µm and dy = 4 µm was used, along with 10³ rays with various values of f_♯. All other simulation characteristics (equations of state, boundary conditions, and simulation domain) were identical to those of the paraxial simulations, as described in Sec. III. Illustrated as black triangles in Fig. 3(a), the centroid deviations resulting from the modified ray tracing module of Hera [using the fluid beam bending rate, Eqs. (7), (10), (19), and (23)] reproduce correctly the temporally averaged theoretical (solid lines) and paraxial (circles) predictions.

It is worth noting that restricting the model to the beam centroid deviation implies to neglect the contributions to beam spreading coming from the different deflections caused by the various speckle intensities and lifetimes. Likewise, the forward Brillouin instability45,46,75 and associated plasma smoothing effects may contribute to reducing the effective speckle waist, thereby affecting the final beam centroid deviation.70Figure 3(a) demonstrates that the bending rate decreases faster than 1/f_♯, and so plasma smoothing effects may greatly worsen the impact of beam bending on the propagation of the laser. Moreover, our crude modeling of the speckle dynamics does not account for the oscillations of the coherence time and speckle velocity as characterized in Ref. 76, which also contribute to the fluctuations in the beam direction.

To quantify the beam bending level under realistic conditions, we performed a 2D-axisymmetric hydrodynamic simulation with the Troll code (see Appendix C for details)38 of Hybrid B NIF Shot N181209,55–57 i.e., with a diamond ablator and a low-gas-fill hohlraum. No significant amounts of backscattering were measured during these shots, and so the beam bending may be addressed here independently from the stimulated backward Brillouin/Raman instability. One representative time was retained, 6 ns, which corresponds to the main power drive, as illustrated in Fig. 4(a). We will focus here on the beam bending of the inner cone, whose main axis lies on y = 0 and propagates from left to right in Figs. 4(b)–4(e). The beam bending deflection rate is calculated using the local plasma and laser parameters (f_♯ = 8), with a speckle coherence time τ_SSD = 7 ps. Use is made of the 3D theoretical kinetic predictions of Eqs. (5), (19), and (23), together with the nonlocal correction of Eq. (10). Note that we do not use the analytical fit of ⟨G⟩τSSD presented in Appendix D, but we solve numerically the quadrature of Eq. (23). The TSSD direction is vertical in the NIF chamber,72 and therefore its direction lies in the (r, z) plane of the hydrodynamic simulation. Owing to the weak effect of the PS, we will thus apply our deflection model with (S, β) = (1.7, 1), according to Table I.

Figures 4(b)–4(e) show the local ⟨Z_i⟩, T_e, n_e/n_c, and y component of the flow velocity, respectively. Note that the color map in (b) is saturated to ⟨Z_i⟩ = 10, and thus the high-Z gold medium is represented by the black areas. The map of the y component of velocity, normalized to the local acoustic speed, suggests that at 6 ns, helium, carbon, and gold, with electron densities ranging from n_e ∼ 0.05 to n_e ∼ 0.2, may be found on the inner beam path with M₀ ∼ 1. Indeed, as the helium gas is being compressed by the expanding ablator [the green region in Fig. 4(b)] and gold wall [the black region in Fig. 4(b)], the electron density remains in the 10% critical density range. Likewise, the electron temperature, higher than 3 keV in these regions, leads to qualitatively small Landau damping rates and therefore to a bending not so far from its large asymptotic limit. Note that at 6 ns, the density of the expanded window is of the order of n_e ∼ 10⁻³n_c [see Fig. 4(b)], which is too small for speckle-scale beam bending to occur.

We performed two 2.5D ray tracing Hera simulations using the plasma profiles from the Troll simulation at 6 ns. These simulations did not solve the hydrodynamic equations. We used a domain size of L_x × L_y = 8 × 3 mm², a mesh size of dx = 9.7 µm and dy = 11.7 µm, an f-number of f_♯ = 8, and a beam waist of 1 mm. The beam was composed of 10³ rays, and all the other simulation characteristics were identical to those of the paraxial simulations. Figures 5(a) and 5(c) present the intensity map and the inverse bremsstrahlung power of the inner cone with the modified scheme of Eq. (24). These are to be compared with the corresponding classical ray propagation results in Figs. 5(b) and 5(d). The intensity profiles evidence a weak impact of the beam bending dynamics on the pointing direction. When account is taken of the beam deviation [Fig. 5(a)], the expanding gold bubble around x = 2 mm tends to deflect the laser toward the ablator. Indeed, Fig. 5(e) shows the local 2.5D bending rate and reveals a slight deflection of ∼−70 mrad/mm located in the gold material. Likewise, the compressed helium and a small part of the ablator tend to deviate the other half of the beam (by ∼+50 mrad/mm) in the direction of the gold wall (upward), thus enhancing the already present focusing effect due to refraction. The temporal incoherence of the laser is thus not strong enough to totally suppress the beam bending and compensate for the thermal effects, especially in the gold medium, leading to significant deflection rates, as illustrated in Fig. 5(e). Quantitatively, the energy deposition is 0.5% weaker in the gold medium as extracted from our ray tracing simulation with our beam bending model than without it, and 10% less power is deposited onto the carbon ablator when the deviation is taken into account. However, the profile of the laser absorption is modified, being much more peaked around (x, y) ≃ (4.5, 0) mm when the beam bending is included. Quantitatively, ∼0.05% of the incoming laser power is able to leave the gold bubble in the presence of beam bending, which greatly contrasts with the ∼2.8 × 10⁻⁴% obtained without flow-induced refraction. These numbers, highly sensitive to the beam bending model, seem to be at least qualitatively consistent with the measurements presented in another experiment.77

Importantly, the part of the beam that is specularly refracted by the gold wall is more energetic according to our model than when the classical ray tracing scheme is used [see Figs. 5(a) and 5(b) around (x, y) ≃ (6, −1) mm]. The small flow-induced deflections in the gold bubble seem to be sufficient, at grating incidence, to significantly modify the penetration depth of the laser in the dense plasma and its absorption. The resulting deflected beam, of intensity in excess of 10¹³ W/cm², could reach the opposite laser entrance hall and therefore could lead to possible wave mixing with opposite-side beams. The striking qualitative similarity with experimental measurements ascribed to stimulated Brillouin side scatter or the so-called glint77–79 indicates an urgent need to address the impact of beam bending on the implosion dynamics in ICF experiments.

On the basis of previous work, we have addressed the beam bending physics of a realistic (RPP, SSD, and PS) high-energy laser beam propagating under ICF plasma conditions. By mean of large homogeneous paraxial hydrodynamic simulations, we have constructed a 3D ray-tracing model that quantitatively captures the temporally averaged centroid deviation of realistic beams. The use of our ray tracing model under ICF Hybrid B conditions reveals a weak modification of the beam propagation properties as a consequence of speckle-scale beam bending. Our results indicate that flow-induced deviation in the gold bubble is a good candidate for explaining the anomalous refraction measurements: the so-called glint.78,79 Our model may thus be used to predict and assess the impact of this anomalous refraction on implosion dynamics and wave mixing processes. The impact of beam bending in other ICF scenarios is left for future work.

Speckle-scale physics may thus be included in a coarse description of laser light, such as ray tracing-based schemes, opening the way for a more realistic laser energy deposition model in radiative hydrodynamic codes. Obtaining more accurate predictions of the implosion symmetry and bang time without resorting to nonphysical artifacts such as laser power multipliers80 may allow smart ICF designs81,82 and possible control of deleterious laser–plasma effects.

Our ray tracing model with spectral dispersion has been compared with homogeneous ideal simulations, and so its validity can be ensured only in smooth density gradients. Furthermore, most wave-mixing processes are absent from our analysis, although these might possibly affect the beam deflection and thus bring further complexity to the system.25 According to Fig. 3(a), the speckle-scale beam bending may be greatly amplified by plasma smoothing effects.70,83 The ICF estimates in Figs. 5(a), 5(c), and 5(e) could thus represents a lower bound on the realistic flow-induced deflection. Since our study is limited to the beam centroid deviation, further analyses, including forward scattering and beam spreading,46,48,84 are required to obtain a deeper understanding of all the physical mechanisms that can affect the beam pointing and subsequent laser energy deposition region inside the hohlraum.

Furthermore, our theory and simulations in realistic geometries allow a comparison of the efficiency of various temporal smoothing techniques that can be applied to achieve control of the beam pointing direction. Compared with LSSD, we have confirmed that 1D TSSD results in an anisotropic deflection rate. The latter can be substantially reduced by increasing the SSD modulation depth to reach field coherence times below 2 ps. Although outside the scope of this study, alternative smoothing strategies such as multicolored beams85–87 could also be of interest with regard to the beam bending dynamics.

Finally, a ray tracing-based simulation incorporating our beam bending model would require evaluation of the quadrature of Eq. (23) each time a ray exits a cell. The exact calculations seem numerically too demanding unless one resorts to a tabulation of the deflection rates or to an accurate fit that would significantly alleviate the computational load. For that purpose, a fit is presented in Appendix D that should facilitate implementation of a ray tracing beam bending model in hydrodynamic codes. A radiative hydrodynamics simulation of a whole-capsule implosion using our beam bending model is left for future work.

Acknowledgment. We acknowledge fruitful discussions with S. Laffite, D. Hinkel, and W. Farmer and important inputs from D. Penninckx and J. M. Di Nicola. We also acknowledge F. Gilleron and R. Piron for the use of the DEDALE model, D. Dureau for his advice on the ray tracing algorithm O. Morice does not need to be in the Acknowledgments if he is in the author list We also acknowledge the HERA team for the development and maintenance of the simulation code. This work has been done under the auspices of CEA-DAM, and the simulations were performed using HPC resources at TGCC/CCRT and CEA-DAM/TERA.