1. Introduction

Filamentation damage limits the output of high-power laser systems, and thus understanding the condition of its formation is key to the design and operation of these systems [1,2]. Filamentary damage in an optic’s bulk is a result of self-focusing, a process in which the Kerr nonlinearity exceeds diffraction and light self-traps into a rapidly tightening spot until the intensity is high enough to generate damage. The damage formation process, once the required intensity has been exceeded, involves multi-photon ionization and generation of a carrier density in the bulk which in turn limits the intensity growth. The physics of damage formation is discussed elsewhere [3–10] and is beyond the scope of this work. Here, we focus on the process only until self-focusing is reached, and thus the underlying physics can be described by the nonlinear Schrödinger wave equation (NLSE). Limiting the physics to this more simplified form will also allow examination of the validity of the rules previously derived under this framework.

The condition leading to filamentation formation depends on the input beam shape in addition to the governing equation, and therefore this study focuses on perturbation growth on a large top-hat like beam shape typical of high-power lasers [11]. The sensitivity of filamentation to beam shape is known and can be traced back to the two-dimensional (2D) nature of the NLSE [[10] chapter 17, [12]]. The volume of the literature on filamentation is too large to be covered here, and we refer the reader to comprehensive reviews of the subject [10]. At large, great efforts were focused on whole beam collapse (especially of Gaussian and super-Gaussian beams), whereas the initial condition of the beam treated here as perturbations to the ideal top-hat shape, is less studied. There are several differences in the dynamics between these two excitation cases; to name few: the scale and shape of the self-collapsing features are different, resulting in a different self-focusing and diffraction balance; for the large beam case, the power in the beam usually exceeds the critical power required for collapse, and therefore the research question becomes what would be the distance to collapse, rather than whether the beam will collapse; the peak intensity in the whole beam case is replaced by two parameters in the large beam case: the beam base intensity and the perturbation amplitude. Even though we refer in this study to material system and dimensions relevant to high power lasers, the conclusions, at large, may also apply to other cases of perturbation collapse in other Kerr nonlinear media, such as atmospheric propagation for free space communications and directed energy applications.

An empirical observation that serves as a powerful design tool is that the product of the beam intensity, I, and the distance to collapse, L, (known also as the self-focusing length) has a constant value, commonly referred to as the IL rule. The basis for this rule (to the best of our knowledge) is given in [13,14] and is based on the asymptotic behavior above the critical power for whole beam collapse studied through numerical simulations. Adaptation of this notion of the IL rule to perturbation growth on a large beam was reported by Milam based on experiments [15,16]. In this context a large beam is incident on a thick fused silica glass rod after passing through a metal wire that generates a spatial perturbation in the beam intensity. The product of the distance to the onset of the observed damage and the input intensity is measured and reported to be about 25 GW/cm at 351nm wavelength. This sets a very powerful design rule, since it allows a prediction of the highest permissible intensities per a given optic’s thickness to avoid filamentation damage. This releases the laser designer / operator from the requirement to control the finer details of the amplitude and phase of the beam entering the optics, which is usually challenging in complex laser systems. However, it raises multiple questions to the limits of the validity of this empirical observation, the interpretation of the intensity (is it the peak local power or the average quantity, and how does it relate to the perturbation amplitude?), and the role that other parameters not included in the expression play (e.g., perturbation characteristics).

Another model of relevance is the Bespalov-Talanov (BT) gain theory for perturbation growth in NLSE type system [17,18]. The classic BT gain expression is a linearized perturbative theory which describes the growth of a small harmonic perturbation due to nonlinear phase accumulation. BT gain has been suggested in the past to describe the dependence of perturbation growth and collapse as a function of its spatial period [20,21]. One unresolved difficulty though, is that during the collapse process, the initial ‘perturbation’ substantially grows beyond the base beam intensity and substantially reshapes, diverting from its original harmonic profile, which puts into question the base assumptions of the BT gain model. In addition, we will also extend the model from its classic 1D harmonic perturbation to a more general 2D arbitrary shape.

Here we will study the validity of the IL rule, BT gain and the relations between them in the excitation settings relevant to high power laser systems using numerical simulations of the NLSE. We derive expressions based on modified and generalized BT gain model and compare them to the simulation result. The goal of this study is to examine the validity of the practiced IL rule, and uncover what additional physics beyond it required to properly design, operate and monitor high power laser systems.

This paper is organized as follows: in section 2, we rederive the BT gain model for a generalized 2D perturbation and modify it to predict the self-focusing length. We then implement this generalized concept for cases of harmonic and Gaussian amplitude perturbations as well as for phase perturbations. In section 3, we compare the IL rule and the modified BT gain model expressions to the numerical simulations results to find which of the properties are captured well by these models. In section 4, we examine a practical implication of the BT gain concept in the form of spatial filtering for suppression of the filamentation problem. In section 5, we examine the effect of the base beam edge on perturbation growth and the local IL rule. In section 6, we examine the effect of two-color cross-phase modulation (XPM) on the BT gain model and the IL rule which is relevant in frequency converted systems when unconverted light is present.

2. Rederivation of a generalized BT gain-based model

In this section, we extend the classic small signal BT gain model beyond 1D harmonic perturbations to a more general perturbation shape and calculate the self-focusing length, L_fil. In principle, since BT gain derivation is a linearized approach, a Fourier decomposition to 2D harmonics should apply. However, since the perturbation growth is only the first step leading to a collapse, which is highly nonlinear, and since the initial conditions for reaching the collapse point depends on the weight of the different frequencies components (intensities, phase, and the resulting interference), it is somewhat harder to think of this problem in terms of a decomposition problem, and it is more illustrative to extend the perturbation spatial shape beyond infinite harmonic wave. This generalization will turn insightful with respect to the local source term of filamentation formation.

The derivation starts with the governing nonlinear wave equation for the phasor representation of the E-field (harmonic in time, $E\propto \exp (i\omega \cdot t-ik\cdot z)$), assuming paraxial propagation along z direction [19]:

The small-signal gain expression suggests that the perturbation growth relates to the local curvature, providing the specific functional dependence (represented by $\varphi ={\nabla }^{2}f/f$). We note that using a 1D harmonic perturbation in Eq. (6) will result in the classic BT gain theory as will be shown later. Assuming the perturbation amplitude grows exponentially with the propagation distance and the small-signal gain, L_fil is reached when the intensity reaches some arbitrary value $I_{fil}=(1+a^2\exp (2g{L}_{fil}))I_0$(to reach filamentation conditions the perturbation amplitude has to substantially exceed the initial perturbation size and thus the cosh takes the simpler form of an exponent). Since the filamentation intensity is substantially larger than the input one, the ‘1’ term can be neglected, and the self-focusing length expression can be simplified:

The gain dependence on the intensity local curvature suggests that the gain maximizes for certain perturbation shape parameters. Heuristically, by adjusting the local curvature of the perturbation (represented by the functional ϕ(f)), the balance between diffraction and self-focusing is altered: too strong of a curvature will diffract strongly, while too weak of a curvature will result in weak lensing. Assuming the perturbation functional ϕ depends on some spatial measure of the perturbation, K, the perturbation will receive a maximal BT gain when ${\partial }_{K}g=0$and ${\partial }_{KK}g<0$. Using calculus, detailed in Appendix A, the maximum gain is:

The IL rule could be interpreted from the modified BT gain model as the damage onset for an excitation of a uniform spatial frequency spectrum of perturbations. Without any a-priori knowledge of the spectral content of the perturbations on the base beam, we choose to assume a uniform spatial frequency spectrum. In this case the damage onset will result only from the perturbation with spatial measure K that sees the maximal BT gain. Under these conditions the expression for the maximal gain is given in Eq. (8a), and Eq. (7) takes the shape of a modified IL rule after some algebra:

Before comparing the model predictions with the numerical simulations results, the equations above will be examined for some particular cases of interest.

2.1 Phase perturbation

The extension from amplitude perturbation to a phase one is almost trivial using a Taylor series for the source term. A harmonic phase perturbation on a base beam expression is, then, $A=A_0\exp \left(ia\cdot f\left(x,y\right)\right)\approx A_0\left(1+ia\cdot f\left(x,y\right)\right)$, and therefore, the BT gain derivation is the same as described before, where the equation of interest is for the imaginary part, $a_I$, rather than the real part, $a_R$. Nevertheless, both have the same small signal gain expression. Therefore, the BT gain expression is the same for phase and amplitude perturbation terms.

Even though the BT gain coefficient expression is the same for both amplitude and phase perturbations, they represent a growth in different physical properties, and thus carry a different weight. While the amplitude growth results in secondary index change and self-lensing through intensity (i.e., through n₂I), the phase growth is with respect to the absolute value of 2π and affects directly phase bending (i.e., through n₀). Therefore, phase effects are expected to be more substantial, even though harder to measure and monitor on complex laser systems.

2.2 2D harmonic amplitude perturbation

The analysis of the tendency of 2D harmonic perturbations to filament is insightful since spatial frequency decomposition is a commonly used concept in optical systems. At least under the linearized case assumption, a spectrum of harmonic perturbations with different spatial frequencies could be treated separately, and we could infer which of the spatial frequencies grows faster than the other. Furthermore, since the BT gain is determined locally (ϕ(f) is local operator), it is the same for infinite and finite number of periods, whereas the latter is more representative of the practical amplitude perturbation observed on a base beam.

In the case of 2D harmonic function, the perturbation function takes the form:

2.3 2D Gaussian amplitude perturbation

To exemplify how the BT-gain expression changes with the postulated spatial shape of the perturbation, we examine a 2D Gaussian shape, which is different from 2D harmonic shape, being non-periodic with varying spatial curvature.

The perturbation in this case takes the form:

 

Fig. 1 numerical calculation of self-focusing length (L_fil) for 2D harmonic perturbation on a flat base beam: (a) illustration of the four parameters varied: I₀, a, W_x, and W_y. (b) four snap shots of the intensification (local intensity I, normalized by the initial base beam intensity I₀) at different propagation length (detailed on their title: 0, 3.043, 4.7247, 5.3623 cm) (see Visualization 1). Run-parameters for this case are: I₀ = 10 GW/cm², a = 0.05, W_x = W_y = 0.15 mm, and the self-focusing length was found to be L_fil = 5.47 cm. (colors throughout this work are added as a visualization aid for the z-axis values of surface plots).

Download Full Size | PDF

3. Comparison of the analytic model to NLSE numerical calculations

In this section we will compare the analytic model derived in the previous section to numerical calculations with respect to the predicted self-focusing length and its dependence on the excitation beam parameters.

3.1 Method of the numerical calculation

The numerical calculation is based on a beam propagation method (BPM) implementation using a split-step algorithm [22–25] of the NLSE (Eq. (1)). The propagation step size is adjusted along the propagation and is kept much smaller than the linear and nonlinear lengths [23]. The numerical lateral resolution is chosen small enough (~few microns) to not introduce noticeable numerical error, and the spatial spectrum is monitored to stay well below the spatial Nyquist frequency.

The excitation is of a 2D harmonic amplitude perturbation on a spatially invariant intensity, representing the top-hat region of a base beam, in accordance to the mathematical representation given in Eq. (2) and (10). For each simulation the parameters that are being changed are: the base beam intensity I₀, the perturbation amplitude a, the perturbation period in x and y, W_x and W_y accordingly (W = 2π/K), as illustrated in Fig. 1(a). The two axis periods are not exclusively the same; however, we use the convention that the longer period is in the y axis. The parameters that are kept constant throughout this paper are: the free space wavelength λ₀ = 351 nm, and the optical material constants at this wavelength (from tabulated data of fused silica glass): n = 1.4767, n₂ = 3.6 × 10²⁰ m²/W [26,27]. For each simulation the domain is fitted to the size of one cycle of the perturbation, and since the boundary conditions are periodic, the scenario calculated is for an infinitely periodic harmonic perturbation on a flat base beam, in accordance to the excitation used for the analytic model. The numerical calculation was validated with the model (Eq. (12)) for the linearized case of a small perturbation growth, and was found to be in excellent agreement (see Fig. 7, Appendix C).

For each simulation, a set of run-parameters are chosen {I₀, a, W_x, W_y}, and the propagation calculation stops when the peak intensity reaches an arbitrary high value, determining L_fil for this case. We choose here an arbitrary value of I = 500 GW/cm² which is more than an order of magnitude higher than the initial base beam intensities examined here. Since the self-focusing in the NLSE governing equations (Eq. (2)) is not limited (unlike the practical case, where carrier density forms at high intensity which limit the intensity growth through deflection and absorption), the growth is exponential, and when it reaches intensities substantially higher than the initial one, the growth is very rapid, and the choice of L_fil has small sensitivity to the particular arbitrary threshold selection. This methodology is a common practice in this field [10,13]. The calculation for one prototypical case is illustrated in the Fig. 1(b) with propagation snapshots (movie online).

3.2 Self-focusing length dependence on the run-parameters

The self-focusing length calculated for different run parameters, shows a strong dependence on the perturbation amplitude and period, which are beyond the physics detailed by the IL rule. The calculation of L_fil, as described in section 3.1, was conducted for different run-parameters: base beam intensities, I₀ = 6 and 10 GW/cm²; perturbation amplitudes, a = 0.05, 0.15, 0.25; x-periods, W_x = 100 μm – 600 μm (in 50μm steps); y-to-x period ratio, W_y / W_x = 1 - 1.5 (in 0.1 steps). The results are presented in Figs. 2(a) – 2(c), as a function of W_x, W_y, and W_BT accordingly. The results clearly, show that for each given base beam input intensity the self-focusing length depends also on the perturbation amplitude and period. Also, if considering peak input intensity (i.e., given combinations of I₀ and a) to be used in the IL rule instead of the classically used base beam intensity, the results still strongly depends on additional parameters, being the x and y periods.

 

Fig. 2 Dependence of self-focusing length on the 2D harmonic perturbation run-parameters: L_fil vs: (a) x-period, (b) y-period, and (c) BT-period; (d) normalized BT model prediction to L_fil vs period curve based on Eq. (12); subset of the Lfil data presented in (a-c) for I₀ = 10 GW/cm², a = 0.05, and W_x≥300μm, vs: (e) y-period; and (f) BT-period.

Download Full Size | PDF

The dependence observed for L_fil on the perturbation period has the prototypical shape of the BT gain model. Before getting into a more quantitative comparison between the model and the numerical simulations, we first quantitatively compare L_fil vs W curve shape. The BT model prediction for the curve shape, based on Eq. (12) and (7) and simplifying to normalized units, is:

The 2D combined period based on the BT model agrees well with the numerical simulation data. The numerical simulations where conducted for different W_x and W_y values. When presenting L_fil as a function of each of them separately, as depicted in Figs. 2(a) and 2(b), it is clear from the y-values spread of the data (i.e., uncertainty) that each one of them separately does not contain all the required information to determine L_fil. The BT model, derived in section 2, suggests that the perturbation growth depends on a combination of the two periods, in a functional shape detailed in Eq. (11). When using this BT combined period representation, the uncertainty is substantially reduced, suggesting that this combined period describes well the filamentation dependencies on the spatial periods.

To further illustrate the reduction of the uncertainty when using the BT representation, we replot the data for one of the curves separately for different x-to-y periods ratio. In Fig. 2(e) and 2(f), the L_fil values for the curve of I₀ = 10 GW/cm² and a = 0.05 is plotted in the range of W_x values between 300 – 600 μm, vs W_y and W_BT accordingly. The different colors represent the different x-to-y period ratios, and clearly shows that the uncertainty vanishes when using the BT representation. Additionally, we have observed some slight differences between the curves for x-to-y period ratios of 1.4 and above, and especially for BT periods closer to the minima period and at {I₀, a} scenarios of leading to the longer L_fil. This suggests that a breakdown of x-y symmetry in the initial excitation conditions might reduce the accuracy of the BT model in describing the filamentation formation.

We emphasis that the agreement we have presented so far is a non-trivial conclusion, since the BT model is based on a linearized perturbative assumption that is clearly disrupted and violated in the highly nonlinear collapse process beyond the initial perturbation growth step.

A quantitative comparison shows that a substantial linear transformation is required for the modified BT model to meet the numerical data. For the more quantitative comparison, additional {I₀, a} run-parameter cases where added, but focused on period range around the gain maxima with finer resolution and for a smaller x-to-y difference. The examination set included: I₀ = 6, 8 and 10 GW/cm²; perturbation amplitudes, a = 0.05 – 0.25 (in 0.05 steps); x-periods, W_x = 100 μm – 300 μm (in 12.5μm steps); y-to-x period ratio, W_y / W_x = 1 - 1.2 (in 0.1 steps). For each run-parameters set a linear function of the y = Ax + B form was applied to the BT model prediction for the L_fil vs Λ_BT curve (based on Eqs. (12) and (7)), and fitting algorithm was used to find the {A, B} constants that minimize the sum of absolute difference between the two. In Figs. 3(a) and 3(b), two prototypical run-parameter cases show the comparison of the numerical data and the fitted modified BT model along with the fitting parameters, showing good agreement with the data.

 

Fig. 3 comparison of the modified BT model and simulations data: (a,b) comparison of self-focusing length data and fitted BT model for two cases ({I₀, a} and fitting parameters are given in plots); (c) fitting parameter A vs a; (d) fitting parameter B vs max input intensity; (e) error in max gain period between model and data; (f) product of input intensity and minimal L_fil for each{I₀, a} case vs a functional dependence on a given in Eq. (9).

Download Full Size | PDF

The fact that only a linear function is required to fit the model is non-trivial. It suggests that the growth of the different perturbation periods under BT model, and the collapse resulting from the different perturbation periods are related, and slower BT growth of perturbation will result in a longer collapse distance. However, it also suggests that there is additional physics that is required to describe the filamentation distance beyond the BT model. The relation between the different fitting parameters was plotted against the different run-parameters {I₀, a}, and significant correlation could be found mainly between A and the perturbation amplitude a, and between B and the input peak intensity I₀^max (Figs. 3(c) and 3(d)). The fact that the slope coefficient values are large (i.e., about 20-40) makes it unlikely for the discrepancy to result from the arbitrary value choice of the I_fil. For example, in a case where A = 30, being an increase of 30 times in L_fil, the discrepancy has to be explained by an increase of exp(30) times in the I_fil (see Eq. (7)). The linear fitting function and its coefficients’ functional dependence on {I₀, a} may point to the following two-step picture of the process, where the perturbation grows at first step into Townes profile [28] and collapse into a filament in the second step. The dependencies suggest that the multiplier A relates to the perturbation amplitude a, perhaps relating to the growth of the perturbation into Townes profile at the first step, where larger initial amplitude requires less filamentation distance (therefore less A). Similarly, the multiplier B is related to the base beam intensity I₀, perhaps, since the Townes profile starting point for the second step of the collapse is common to all the cases, and since less distance (therefore less B) is required for more intense base beam.

The BT model is in reasonable agreement with the numerical data for the peak growth spatial period. The error between the period of minimal L_fil (from the numerical data) and the peak gain period (from non-fitted BT model, i.e., Eq. (13)) is depicted in Fig. 3(e). The error values are in the range 8-15%, and larger for larger a. These error values suggest the BT gain model gives a good practical prediction of the period most likely to filament. Nevertheless, the non-negligible error and its dependence on the perturbation amplitude, both suggest that a physics beyond the BT model is required to describe it more accurately.

The modified IL rule corrected to the dependence on a (Eq. (9)) describes well the damage onset for uniform spatial frequencies spectrum of same amplitude perturbations. For the numerical data for each {I₀, a} case, the product of the minimal L_fil and I₀ is depicted vs the predicted function dependence on a from Eq. (9). For each perturbation amplitude, the IL product is invariant in the input intensity. Furthermore, the dependence on the perturbation amplitude is consistent with the predicted in Eq. (9). This conclusion gives a more rigorous context to the conditions for which the IL rule could be used.

Similar BT gain like behavior was also observed for the 2D harmonic phase and 2D amplitude Gaussian perturbation spatial functions. The simulation results for these cases are given in Fig. 8, Appendix D.

4. Effect of spectral filtering on self-focusing length

The results above suggest that spatial filtering might suppress filamentation. The self-focusing length depends on the perturbation period, such that periods longer than the BT gain max period the L_fil increases monotonically (~linearly) with the period. For example, for {I₀, a} = {6 GW/cm², 0.05}, L_fil changes from about 10 cm to about 15 cm as the spatial period changes from 200 μm to 400 μm (see Fig. 2(c)). Therefore, spatial low pass filtering (LPF) could increase significantly L_fil. This hypothesis is based on the analysis and numerical calculations above that studied the filamentation dynamics of each harmonic separately. Therefore, studying the effect of spectral filtering on L_fil, has the added value of concurring that this hypothesis agrees with simultaneous excitation of continuous spatial spectrum of harmonic perturbations. That is, off course, beyond the practical utility of the understanding developed here to improving laser system using tailored spatial filtering.

The numerical calculations agree with the hypothesis, showing that tighter spatial filter pinhole will result in larger L_fil, as illustrated in Fig. 4(a). To study the spatial filtering, a 2D white noise was generated with uniform random distribution of spatial spectrum (on a base beam with given intensity, and perturbation amplitude defined with respect to the base beam amplitude, similarly to the harmonic perturbation), it was low pass filtered to remove the periods shorter than the designated BT period (W_x = W_y was assumed), and then used as the excitation input to the simulation described in section 3.1. The size of the calculation grid was kept the same for all the cutoff periods, 2 mm on 2 mm. Because of the random nature of the excitation, each run-parameter set was repeated 9 times.

 

Fig. 4 filamentation suppression with tighter spatial filtering: L_fil vs LPF cutoff W_BT (a) without and (b) with renormalization of the input signal rms to a = 0.15. W_x = W_y = 0.1mm ÷ 0.5mm, I = 8GW/cm², repeat of each parameter set 9 times for statistics. (c) illustration of the filamentation formation dynamics for one renormalization case LPF cutoff W_BT = 0.25mm (see Visualization 2, and for LPF cutoff W_BT = 0.5mm see Visualization 3).

Download Full Size | PDF

The increase in L_fil with spatial filtering is a result of both, the removal of the more prone to filament periods (as discussed above) as well as the reduction of the spatial signal rms (root mean square amplitude). According to the Parseval’s theorem the total power in the Fourier spectrum equal the spatial one, and the rms amplitude of the random signal in space will be reduced due to the power reduction resulting of the filtering. To try determining which of the two discussed effects impacts the increase in L_fil, the signal was renormalized to the same amplitude (a = 0.15) after filtering, which eliminates the reduction of spatial rms effect (even though this setting represents less accurately the function of a practical low pass filter). The rest of the methodology was the same. The results are presented in Fig. 4(b), elucidating that both effects have a substantial contribution. In Fig. 4(c), the initial excitation and the filamentation dynamics is illustrated with three propagation snap shots (and movie online). The minima in Fig. 4(b) is likely a result of the method of normalization and the BT max gain shape: for cases where the cutoff period is lower than the BT max gain results in less power at the particular period of max gain with respect to the case of cut off equals the BT max gain, which dominantly determines L_fil. The broadening of the L_fil uncertainty for larger cutoff periods is a result of statistics: the larger the periods, the fewer the peaks that could potentially filament for a given calculation grid (which also represent the situation for real optics with a given size).

5. Effect of base beam shape on self-focusing length

The base beam shape could also change the conditions for the formation of the filamentation. The study in the previous section focused on a perturbation on top of a spatially invariant base beam. That situation represents an approximation of the flat main central region of a top-hat beam. Another spatial region in the top-hat beam shape is the roll-off region at the beam edge. These regions in the beam are observed to have a tendency for filamenting thus motivating a further examination [2] (we note that the beam roll-off segments used here are on a mm scale, much smaller than the ~cm ones used on the NIF beam).

To study the beam role off regions, a 2D Gaussian perturbation on a 1D top-hat like base beam profile was used as the initial field excitation. The Gaussian perturbation was used to make the perturbation finite. The base beam profile is invariant in 1D (x-axis) and has a hyperbolic tangent profile (y-axis) to represent the prototypical shape of the beam roll-off segment, where periodic and absorbing boundary conditions are used in x and y axis, accordingly. The prototypical normalized intensity profile is illustrated in Fig. 5(a).

 

Fig. 5 filamentation dependence on its location on the base beam edge roll-off: (a) prototypical normalized input intensity; (b) base beam amplitude profile vs initial perturbation central location; (c) propagation snapshots for I₀ = 10 GW/cm², W_x = 25 μm, a = 0.15, initial location −0.25mm (see Visualization 4, and see Visualization 5 as a reference for initial location 0); (d) localized IL factor, being the product of intensity at the initial filamentation location (normalized by the peak intensity for the base beam) and L_fil (normalized by L_fil value at the base beam center). Base beam roll-off width is 1mm.

Download Full Size | PDF

The effect of the perturbation location on its growth and filamentation formation is studied by repeating the simulation for spatially shifted perturbation on the roll-off profile. The location of the Gaussian perturbation of a given peak amplitude (a refers to the base beam peak intensity, not the local intensity at its location) with reference to top roll-off position, as illustrated in Fig. 5(b). The parameter set includes: Gaussian waist W_x = W_y = 25μm, 50μm, 75μm (see reference values for Gaussian perturbation in Fig. 9, Appendix E); peak base beam intensity I₀ = 6 and 10 GW/cm²; a = 0.05, 0.15, 0.25; and, the beam roll-off is 1 mm. The filamentation formation dynamics for one case is illustrated in Fig. 5(c), in snapshots along the propagation direction (movie online).

The IL product reduces in the roll-off regions, even though the self-focusing length increases. The simulation was repeated with the perturbation located at different initial locations at shifts of 50μm between the roll-off starting point (location = 0) and half the base beam intensity (location = −0.5 mm). The self-focusing length increases monotonically with reduction in the intensity profile for all the parameter set conditions used (see Fig. 10, Appendix F). When normalizing the L_fil at each location by the L_fil value at the beam center for this parameter set {W_x, I₀, a}, and multiplying it by the intensity reduction compared with the beam center, we get a location modified IL product, depicted in Fig. 5(d). The reduction in the location modified IL product suggests that the tendency to filament increases moving into the roll-off – that is, when comparing cases of same base beam local intensity at the different locations. This is possibly a result of the base beam curvature increasing the local curvature of the perturbation, and thus the filamentation tendency, in agreement with the model detailed in section 2. Nevertheless, the intensity profile drop is more dominant than the reduction in the local IL product, resulting in an overall increase in L_fil moving into the roll-off. Perhaps, a further study of this effect could explain what conditions of the roll-off profile might result in enhanced tendency to filament at the beam edge.

6. Effect of cross phase modulation from additional frequency harmonic on L_fil

Some of the high-power laser systems frequency-convert the beam at the final optics, and due to finite conversion efficiency, two frequency harmonics with similar spatial shape co-propagate. For example, in the national ignition facility (NIF), the beam propagates through most of the beam path in 1053nm wavelength (i.e., 1ω), and is mostly converted (~75%) only at the final optics assembly to its third harmonic, (i.e., 3ω, 351nm wavelength) [11]. The two harmonics are co-propagating coupled by cross-phase modulation (XPM) [9,23], in addition to exhibiting the Kerr self-focusing (i.e., self-phase modulation, SPM) and diffraction as described in Eq. (1). This motivates the study of the effect of the XPM on the filamentation formation. It is interesting for practical reasons to try evaluating if the IL rule using the intensity (I) as the total intensity of the two harmonics could be used as an approximation or even as an asymptotic limit.

To study this effect, we use a set of two NLSEs coupled by XPM, and a split step algorithm implementation as described in previous sections where γ₁ is the Kerr coefficient for 3ω, γ₂ is the Kerr coefficient for 1ω, and the XPM term linking the two are γ₁₂ and γ₂₁.

We have found that cross-phase modulation could enhance or inhibit the filamentation formation, depending on the perturbation period. We study L_fil dependence on W_x with and without XPM, as presented in Table 1. Two cases are compared: when no additional harmonics is presented (η_1ω = 0) and 25% power is in the additional harmonics at excitation (η_1ω = 25%), both for I₀ = Ι_1ω + Ι_3ω = 6 GW/cm², (Ι_1ω = η_1ω⋅I_0, Ι_3ω = η_3ω⋅I₀), a = 0.25 (whereas, the effect for other harmonics ratios is illustrated in Table 2, Appendix G). The table data elucidates that as the period is increased above about 200 μm the XPM further enhance the filamentation (rows shown in red), while as the period is reduced below about 200 μm the XPM further suppress the filamentation formation. Note that the minimal L_fil period increases when unconverted light is present, from about 200 μm (when η_1ω = 0) to about 250 μm (when η_1ω = 25%).

Table 1. L_fil vs perturbation period, with (η_1ω = 25%) and without (η_1ω = 0) XPM; I₀ = 6GW/cm², a = 0.25. Rows in red show cases where cross-phase modulation of converted 1ω light makes filamentation worse.

View Table | View all tables in this article

Table 2. L_fil vs XPM; W_x = 750μm, I₀ = 6GW/cm², a = 0.25

View Table | View all tables in this article

Examination of the simulation dynamics reveals that enhanced diffraction at the infrared harmonics may explain the L_fil dependence on XPM. In Fig. 6, propagation snapshots (and movies online) are shown for very short period (i.e., 150 μm) and for very long period (i.e., 750 μm), with and without XPM, after 5 cm of propagation. At the shorter perturbation period, the diffraction of the additional harmonic dominates its dynamics (having a longer wavelength) and breaks apart the central peak, which creates a self-defocusing effect on the main harmonic, thus suppressing filamentation formation. At larger period sizes, both frequencies tend to self-focus, and thus, mutually enhance the self-focusing of the coupled system. This results in a shorter self-focusing length. As a result, the IL product of one harmonic system could be increased or reduced substantially in the presence of two frequencies, even if the total power / intensity is kept the same. The effect will depend on the feature size of the perturbation, where smaller features will inhibit the self-focusing effect, while larger features will enhance it.

 

Fig. 6 filamentation dependence on XPM for (a) short period of 150 μm and (b) long period of 750 μm. Propagation snapshot for the two harmonics are presented after 5 cm, without XPM η_1ω = 0 (center column), and with XPM η_1ω = 25% (right column). Left column shows the intensity distribution at z = 0. I₀ = 6 GW/cm², a = 0.25.

Download Full Size | PDF

7. Conclusion

This study examines and generalizes the well-known BT-gain and IL-rule for filamentation under realistic conditions associated with high power lasers: filamentation seeded by both amplitude and phase perturbations on a large, flat-top beam, and the impact of cross-phase modulation from unconverted light in UV frequency-converted lasers. Studying this excitation sensitive problem using a high-power lasers prototypical input beam (unlike the more commonly studied whole beam collapse), provides a more directly linked interpretations to designing and operating these systems. The validity of the useful practical IL rule is examined with NLSE numerical calculations. We also derive expressions for a generalized BT gain model and compare them to the simulations results, which suggest that the perturbation growth relates to its local intensity curvature. This generalization addresses spatial shaped perturbation beyond the originally BT gain described harmonic functions, and therefore the role of practical perturbation shape deviation from ideal harmonic shape. The IL rule is further derived from the BT gain model under the context of damage onset for a uniform spatial frequency spectrum of perturbations. This newly derived form of IL rule establishes a framework for its validity scope and provides the dependence on the perturbation amplitude. We exemplify the general BT gain expressions for specific 2D perturbations amplitude shapes (harmonic and Gaussian), and phase perturbations. We quantitatively examine the ability of this BT-gain based IL-rule to predict the numerical simulations results, which is key to high-power lasers design. The simulations and model show that there are parameters beyond the commonly-used IL rule, such as the perturbation amplitude and period content, however it does apply to the damage onset of a system with a given perturbation spectrum. The BT model presents a fair description of the tendency of spatial periods to filament, but not of the quantitative self-focusing length, unless an empirical linear transformation is applied. The fact that only a linear transformation is required, and the correlation we report to the fitting parameters on the excitation parameters might suggest that some future modification to the model could perhaps improve the model accuracy. In concert with this description, spatial filtering of short periods is shown to suppress filamentation, due to both, the removal of the more prone to filament periods, as well as the reduction of the spatial intensity amplitude rms. At the edge of a top hat beam we find that the IL product reduces in the roll-off regions, even though the self-focusing length increases. When adding a co-propagating harmonic, we find that the cross-phase modulation (XPM) could enhance or inhibit the filamentation formation, depending on the perturbation period. This study highlights the underlying physics involved in the creation of filamentation in high power and high energy laser systems and provides guidance to designing and operating such systems and the physical parameters to monitor in order to control this significant problem.

Appendix A conditions for maximal BT gain

We rewrite Eq. (6) expression for the gain as:

(18)$$g=\frac{\sqrt{\phi }}{2k};\phi =-\varphi \cdot \left(C+\varphi \right);C=2\gamma I_0;$$

The maximum gain is obtained when ∂_Λg = 0. Since:

(19)$${\partial }_{K}g={\partial }_K\phi \cdot \left(\frac{1}{4k\sqrt{\phi }}\right)$$

The maximum gain condition is equivalent to ∂_Λφ=0. This becomes ill-defined when φ=0 (ϕ=0 or -C), however, at this case the gain is identical zero, and thus not interesting case. By substituting φ value, the condition for extremum point in gain translates to:

(20)$${\partial }_{K}g=0\iff {\partial }_K\varphi \left(C+2\varphi \right)=0(\varphi \ne 0,-C)$$

For the extremum point to be maxima ∂_KKg should be negative. Applying an additional derivative to Eq. (19):

(21)$${\partial }_{KK}g={{\partial }_{KK}\phi \cdot \left(\frac{1}{4k\sqrt{\phi }}\right)-{\left({\partial }_K\phi \right)}^2\cdot \left(\frac{1}{8k\phi \sqrt{\phi }}\right)|}_{{\partial }_K\phi =0}={\partial }_{KK}\phi \cdot \left(\frac{1}{4k\sqrt{\phi }}\right)$$

Therefore, the sign of ∂_KKg equals the sign of ∂_KKφ, which reads:

(22)$${\partial }_{KK}\phi =-\left(C+2\varphi \right){\partial }_{KK}\varphi -2{\left({\partial }_K\varphi \right)}^2$$

Examining the sign of the second derivative for the two conditions at the extrema point (based on Eq. (20)):

(23)$$\begin{array}{l}{{\partial }_{KK}\phi |}_{\left(C+2\varphi \right)=0}=-2{\left({\partial }_K\varphi \right)}^2<0\to \text{max}\text{point}\\ {{\partial }_{KK}\phi |}_{{\partial }_K\varphi =0}=-\left(C+2\varphi \right){\partial }_{KK}\varphi \end{array}$$

For prototypical perturbation function, the curvature function ϕ dependence on the scale variable Λ is such that ∂_Kϕ=0 only when K=0. For example, for harmonic function, ϕ=-K², and thus ∂_Kϕ=-2K which vanishes only with -K, which could be considered a non-practical case. Limiting the conclusion to the maxima detailed by the first condition, the conditions for maximal BT-gain are:

(24)$$\varphi =-\gamma I_0$$

And the max gain expression is:

(25)$$g_{\max }=\frac{\gamma I_0}{2k}$$

However, we should note that if the perturbation function results in ∂_Kϕ=0 where ∂_KKϕ≠0 and where K≠0, another gain maxima point may exist.

Appendix B agreement of Eq. (12) with classic BT-gain expression from literature

Starting with 12 and using K instead of K_BT to refer to the 1D case, writing explicitly the expression for γ, and after some algebra the BT-gain expression is:

(26)$$g=K\sqrt{\frac{\left(4k^2{n}_2{I}_0/n-K^2\right)}{4k^2}}=K\sqrt{\frac{n_2{I}_0}{n}-\frac{K^2}{4k^2}}$$

Which is agrees with the expression for classic BT-gain, given as Eq. (1) in [20].

Appendix C validation of numerical calculations with Eq. (12) for perturbation growth (linear regime)

 

Fig. 7 The increase in an harmonic perturbation amplitude after a distance L = 1 cm, as a function of the perturbation period: numerical simulation results (stars) compared to Eq. (12) (dashed line). The examination here is in the regime of small signal growth: initial perturbation amplitude is a = 0.01 (2 orders of magnitude smaller than the base beam amplitude), and the growth in the perturbation amplitude is smaller than 1.25 times its original size.

Download Full Size | PDF

Appendix D dependence of the linear fit coefficients on the parameter set

 

Fig. 8 Dependence of the fit parameters (A, B) on the parameter set: perturbation amplitude, a, and peak input intensity I_in^max. The 25 points on the plots represents the different initial settings: a = 0.05 – 0.25 (in 0.05 steps), and I₀ = 6 – 10 (in 1 steps) GW/cm². The fit coefficients of the y-values linear function were optimized to reduce the error between the BT gain based model and the data based on the simulations.

Download Full Size | PDF

Appendix E parameter set plots for 2D harmonic phase and amplitude 2D Gaussian perturbation

 

Fig. 9 Dependence of the L_fil vs W_BT for different initial parameter sets {a,I₀}for: (a) 2D harmonic phase, and (b) 2D amplitude Gaussian spatial perturbation functions (note that a Gaussian beam waist resembles about ¼ a period of an harmonic function, which agrees with the scaling of the x-axis).

Download Full Size | PDF

Appendix F dependence of L_fil on perturbation location on the beam edge

 

Fig. 10 (a) dependence of the L_fil on the location of the Gaussian perturbation on the edge. (b) same as a, but each curve of L_fil vs perturbation location is normalized by the L_fil for this curve at the beam center.

Download Full Size | PDF

Appendix G effect of changing XPM on the filamentation

Funding

US Department of Energy.

Acknowledgments

The authors would like to acknowledge the useful conversations with C. Widmayer, C. W. Carr and W.H. Williams. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. LLNL-JRNL-764047.

References

1. K. R. Manes, M. L. Spaeth, J. J. Adams, M. W. Bowers, J. D. Bude, C. W. Carr, A. D. Conder, D. A. Cross, S. G. Demos, J. M. G. Di Nicola, S. N. Dixit, E. Feigenbaum, R. G. Finucane, G. M. Guss, M. A. Henesian, J. Honig, D. H. Kalantar, L. M. Kegelmeyer, Z. M. Liao, B. J. MacGowan, M. J. Matthews, K. P. McCandless, N. C. Mehta, P. E. Miller, R. A. Negres, M. A. Norton, M. C. Nostrand, C. D. Orth, R. A. Sacks, M. J. Shaw, L. R. Siegel, C. J. Stolz, T. I. Suratwala, J. B. Trenholme, P. J. Wegner, P. K. Whitman, C. C. Widmayer, and S. T. Yang, “Damage mechanisms avoided or managed for NIF large optics,” Fus. Sci. Technol. 69(1), 146–249 (2016).

2. J. M. Di Nicola, T. Bond, M. Bowers, L. Chang, M. Hermann, R. House, T. Lewis, K. Manes, G. Mennerat, B. MacGowan, R. Negres, B. Olejniczak, C. Orth, T. Parham, S. Rana, B. Raymond, M. Rever, S. Schrauth, M. Shaw, M. Spaeth, B. Van Wonterghem, W. Williams, C. Widmayer, S. Yang, P. Whitman, and P. Wegner, “The National Ignition Facility: Laser performance status and performance quad results at elevated energy,” Nucl. Fusion 59(3), 032004 (2019). [CrossRef]  

3. G. N. Steinberg, “Filamentary tracks formed in transparent optical glass by laser beam self-focusing. I. Experimental investigation,” Phys. Rev. A 4(3), 1182–1194 (1971). [CrossRef]  

4. A. V. Smith and B. T. Do, “Bulk and surface laser damage of silica by picosecond and nanosecond pulses at 1064 nm,” Appl. Opt. 47(26), 4812–4832 (2008). [CrossRef]   [PubMed]  

5. L. Sudrie, M. Franco, B. Prade, and A. Mysyrowicz, “Study of damage in fused silica induced by ultra-short IR laser pulses,” Opt. Commun. 191(3-6), 333–339 (2001). [CrossRef]  

6. L. Sudrie, A. Couairon, M. Franco, B. Lamouroux, B. Prade, S. Tzortzakis, and A. Mysyrowicz, “Femtosecond laser-induced damage and filamentary propagation in fused silica,” Phys. Rev. Lett. 89(18), 186601 (2002). [CrossRef]   [PubMed]  

7. D. G. Papazoglou, I. Zergioti, S. Tzortzakis, G. Sgouros, G. Maravelias, S. Christopoulos, and C. Fotakis, “Sub-picosecond ultraviolet laser filamentation-induced bulk modifications in fused silica,” Appl. Phys., A Mater. Sci. Process. 81(2), 241–244 (2005). [CrossRef]  

8. S. Tzortzakis, L. Sudrie, M. Franco, B. Prade, A. Mysyrowicz, A. Couairon, and L. Bergé, “Self-guided propagation of ultrashort IR laser pulses in fused silica,” Phys. Rev. Lett. 87(21), 213902 (2001). [CrossRef]   [PubMed]  

9. Y. R. Shen, The Principles of Nonlinear Optics, (Wiley, 2002)

10. R. W. Boyd, S. G. Lukishova, and Y. R. Shen, Self-focusing: past and present, (Springer, 2009)

11. M. L. Spaeth, K. R. Manes, D. H. Kalantar, P. E. Miller, J. E. Heebner, E. S. Bliss, D. R. Spec, T. G. Parham, P. K. Whitman, P. J. Wegner, P. A. Baisden, J. A. Menapace, M. W. Bowers, S. J. Cohen, T. I. Suratwala, J. M. Di Nicola, M. A. Newton, J. J. Adams, J. B. Trenholme, R. G. Finucane, R. E. Bonanno, D. C. Rardin, P. A. Arnold, S. N. Dixit, G. V. Erbert, A. C. Erlandson, J. E. Fair, E. Feigenbaum, W. H. Gourdin, R. A. Hawley, J. Honig, R. K. House, K. S. Jancaitis, K. N. LaFortune, D. W. Larson, B. J. Le Galloudec, J. D. Lindl, B. J. MacGowan, C. D. Marshall, K. P. McCandless, R. W. McCracken, R. C. Montesanti, E. I. Moses, M. C. Nostrand, J. A. Pryatel, V. S. Roberts, S. B. Rodriguez, A. W. Rowe, R. A. Sacks, J. T. Salmon, M. J. Shaw, S. Sommer, C. J. Stolz, G. L. Tietbohl, C. C. Widmayer, and R. Zacharias, “Description of the NIF Laser,” Fus. Sci. Technol. 69(1), 25–145 (2016). [CrossRef]  

12. T. D. Grow, A. A. Ishaaya, L. T. Vuong, A. L. Gaeta, N. Gavish, and G. Fibich, “Collapse dynamics of super-Gaussian Beams,” Opt. Express 14(12), 5468–5475 (2006). [CrossRef]   [PubMed]  

13. E. L. Dawes and J. H. Marburger, “Computer Studies in Self-Focusing,” Phys. Rev. 179(3), 862–868 (1969). [CrossRef]  

14. H. Bercegol, L. Lamaignere, B. Le Garrec, M. Loiseau, and P. Volto, “Self-focusing and rear surface damage in a fused silica window at 1064 nm and 355 nm,” SPIE 4932, 0277 (2003).

15. D. Milam, J. T. Hunt, K. R. Manes, and W. H. Williams, “Modeling of filamentation damage induced in silica by 351-nm laser pulses,” SPIE 2966, 425–428 (1996).

16. W.H. Williams, K. R. Manes, J. T. Hunt, P.A. Renard, S. Milam, and D. Eimerl, “Modeling of self-focusing experiments by beam propagation codes,” UCRL-LR-105821–96–1 (1995).

17. V. I. Baspalov and V. I. Talanov, “filamentary structure of light beams in nonlinear liquids,” JETP Lett. 3, 307–310 (1966).

18. J. Trenholme, “Review of small signal theory,” Laser Program Annual Report 1974 (UCRL-50021–74), 179–191 (1974).

19. J. T. Hunt, K. R. Manes, and P. A. Renard, “Hot images from obscurations,” Appl. Opt. 32(30), 5973–5982 (1993). [CrossRef]   [PubMed]  

20. E. S. Bliss, D. R. Speck, J. F. Holzrichter, J. H. Erkkila, and A. J. Glass, “Propagation of a high-intensity laser pulse with small-scale intensity modulation,” Appl. Phys. Lett. 25(8), 448–450 (1974). [CrossRef]  

21. P. J. Wegner, M. D. Feit, J. A. Fleck, and D. Eimerl, “Measurement of the Bespalov-Talanov gain spectrum in a dispersive medium with large n₂,” SPIE 2633, 661–667 (1997).

22. E. Feigenbaum and T. A. Laurence, “Filament damage formation in fused silica glass as a result of 1-50 ps near-infrared laser pulses,” Appl. Opt. 56(13), 3666–3672 (2017). [CrossRef]   [PubMed]  

23. G. P. Agrawal, Nonlinear fiber optics, (Academic, 1989)

24. R. A. Fisher and W. Bischel, “The role of linear dispersion in plane-wave self-phase modulation,” Appl. Phys. Lett. 23(12), 661–663 (1973). [CrossRef]  

25. M. D. Feit and J. A. Fleck Jr., “Light propagation in graded-index optical fibers,” Appl. Opt. 17(24), 3990–3998 (1978). [CrossRef]   [PubMed]  

26. D. Milam, “Review and assessment of measured values of the nonlinear refractive-index coefficient of fused silica,” Appl. Opt. 37(3), 546–550 (1998). [CrossRef]   [PubMed]  

27. E. D. Palik, Handbook of Optical Constants of Solids, (Academic, 1985)

28. K. D. Moll, A. L. Gaeta, and G. Fibich, “Self-similar optical wave collapse: observation of the Townes profile,” Phys. Rev. Lett. 90(20), 203902 (2003). [CrossRef]   [PubMed]