Small-scale self-focusing (SSSF) is, actually, filamentation instability of a plane wave propagating in a medium with cubic nonlinearity. The SSSF manifests itself as a fast amplification of wave components with high spatial frequencies, which ultimately leads to a beam splitting into filaments, thus deteriorating the beam quality and significantly increasing the risk of optical elements breakdown. SSSF was first theoretically predicted by Bespalov and Talanov^[1^], experimentally confirmed in Refs. [2–6], and the theory was further developed in Refs. [7–12]. At the linear stage of development, SSSF is reduced to the energy transfer from a plane wave to spatial noise, that is, to the amplification of seed noise waves. The gain K depends on the angle θ = k_⊥/k between the wave vectors of the noise component k_⊥ and plane wave k = n₀k₀ (where n₀ is the linear index of refraction). In the linear theory, K has a maximum value at the angle as follows:

$$\begin{align}{\theta}_{\mathrm{max}}=\sqrt{\frac{2{n}_2I}{n_0}},\end{align}$$ (1)

$$\begin{align}K\left({\theta}_{\mathrm{max}}\right)=\cosh\ (2B),\end{align}$$ (2)

$$\begin{align}B={k}_0{n}_2 Il,\end{align}$$ (3)

In classical SSSF works, nanosecond pulses were considered, where the beam was usually split into filaments at B > 3. SSSF studies for femtosecond pulses^[13^–16^] showed that there is a significant difference between nano- and femtosecond pulses, which allows using beam self-filtering as a method for SSSF suppression. There are two mechanisms of beam self-filtering during propagation in free space. In the first (which we will call spatial self-filtering) the noise components propagate at an angle to the main beam and, if the beam passed a rather large length of free space, these components go outside the beam aperture. According to the second mechanism, called temporal self-filtering, the noise components lag behind the main beam, which, in the case of short pulses, allows reducing noise at the input to the nonlinear medium and, hence, suppressing the SSSF. The efficiency of both self-filtering methods depends on the angle θ_max. For nanosecond pulses, a typical θ_max value is about 1 mrad, which impedes self-filtering as it demands a too large propagation length. At the same time, for femtosecond pulses, the intensity is as a rule three orders of magnitude higher and, according to Equation (1), the typical values of θ_max are tens of mrad, which makes self-filtering an efficient mechanism of noise reduction for θ of the order of θ_max.

The linear theory of SSSF gives a correct result only at an early stage of instability development, when the power of the amplified noise is low. From a practical point of view, this stage is not so interesting, since because of small noise the beam retains an acceptable quality. Of major interest is the determination of the B-integral values at which hot spots leading to optical breakdown appear in the beam, or at which a significant part (e.g., 10%) of the plane wave energy is converted into noise. Thus, to quantify the efficiency of self-filtering for suppressing SSSF, the nonlinear regime should be considered. It is also worth mentioning that the experimental results demonstrated a significant suppression of SSSF, which cannot be explained by the existing theoretical concepts^[17^]. One of the goals of our paper is to propose such a concept.

This work is devoted to a detailed numerical simulation of the nonlinear stage of SSSF, taking into account various self-filtering parameters of the input pulse. The conditions under which self-filtering leads to significant suppression of the filamentation instability are found, and a new type of instability is discovered.

The propagation of the field in a medium with Kerr nonlinearity will be studied using the nonlinear Schrödinger equation (NSE)^[18^] of the following form:

$$\begin{align}\frac{\partial A}{\partial z}+\frac{i}{2k}{\Delta }_{\perp }A+i\frac{3\pi {\omega}_0{\chi}^{(3)}}{2n\left({\omega}_0\right)c}{\left|A\right|}^2A=0,\end{align}$$ (4)

At the entrance to the nonlinear medium (z = 0), the field is represented as a plane wave with superimposed noise. The noise is determined by the manufacturing quality of the optical elements through which the beam has passed, and is usually characterized by a 1D power spectral density PSD(k_x). The PSD(k_x) function typical of the optical elements is well approximated by the power dependence^[19^–21^]:

$$\begin{align}\mathrm{PSD}\left({k}_x\right)=\frac{\Phi_1}{k_x^{\beta}},\end{align}$$ (5)

$$\begin{align}\mathrm{PSD}2\left({k}_{\perp}\right)=\frac{\Phi_2}{k_{\perp}^{\beta +1}},\end{align}$$ (6)

The beam self-filtering is modeled using the expression for temporal self-filtering^[23^] at which the free space for a Gaussian pulse in the time domain is the filter of spatial frequencies (angles θ) with transmittance:

$$\begin{align}{T}_{\mathrm{time}}\left(\theta \right)=\exp \left(-{\left(\ln (2)\frac{L}{\lambda}\frac{\theta^2}{N}\right)}^2\right)=\exp \left(-\frac{\theta^4}{\theta_{\mathrm{thr}}^4}\right),\end{align}$$ (7)

Equation (7) was obtained for an ideal Gaussian pulse. In practice, the contrast of femtosecond pulses is limited, as a pedestal outside the main pulse restricts the filter contrast (Equation (7)), which was taken into consideration in the numerical simulations. It is worth noting that at spatial filtering the free space is a filter with transmittance T_space(θ), the shape of which^[16^] is close to Equation (7), but the contrast may be much higher, as there is almost no spatial pedestal. An important feature of both T_time(θ) and T_space(θ) is a very sharp decrease when θ exceeds the threshold value θ_thr.

When a laser pulse propagates through optical elements, inaccuracy in manufacturing the shape of optical surfaces gives rise to phase noise (the imaginary part of the complex field amplitude) and various scratches and dust particles to amplitude noise (the real part of the complex field amplitude). During the propagation in free space, the phase noise and amplitude noise transform into each other; therefore, after passing through a large number of optical elements, it can be assumed with a high degree of certainty that the noise is half phase and half amplitude. It is this noise that was used in the simulation of the SSSF. In this work, we will neglect the noise arising from the refraction at the input surface of a nonlinear medium.

An example of intensity and phase distribution in the beam at the input of the nonlinear element is given in Figure 1. The noise in Figure 1 has been filtered according to Equation (7) with the self-filtering threshold θ_thr = 4 mrad. The input noise power P_noise is about 0.02% of the power P₀ of the principal wave. This noise gives the ratio of the maximum intensity in the beam I_max to the mean intensity I₀ of about 1.07.

The main goal of the simulation was to see the difference in the development of instability in two cases: θ_max< θ_thr = 4 mrad and θ_max > θ_thr. The first case is typical for nanosecond and picosecond pulses, for which θ_max is about 1 mrad; therefore, self-filtering does not affect the noise components with maximal increment. The case with θ_max > θ_thr is typical for femtosecond pulses, for which θ_max is usually tens of mrad and the noise components with maximal increment are filtered at the input to the nonlinear medium, and only components with a much lower gain remain. Hence, we performed numerical simulations for four θ_max values: 1, 3, 10 and 30 mrad. If fused silica is used as the nonlinear medium, then the above values of θ_max correspond to intensity I₀ = 3, 27, 300 and 2700 GW/cm², respectively. The calculations were performed up to the values of the B-integrals, at which a significant power of the principal wave was converted into noise. With a further increase in the B-integral, the sizes of bright spots become smaller than the cell size of the computational grid, which varied from 2.5 to 20 μm in different calculations.

The calculated noise spectra S(θ) = PSD2(k_⊥ = kθ) are shown in Figures 2(a) and 3(a) for the original beam (B = 0) and the beam that has passed in a nonlinear medium distance corresponding to B = 3, 4, 5, 6; the self-filtering contrast is 10⁸. Figure 2(a) is for θ_max = 1 mrad and Figure 3(a) is for θ_max = 3 mrad. The dashed curves correspond to the gain K(θ) in the linear (without plane wave depletion) regime^[18^]:

$$\begin{align}K=1+\frac{2}{\frac{2{\theta}^2}{\theta_{\mathrm{max}}^2}-\frac{\theta^4}{\theta_{\mathrm{max}}^4}}{\cdot \sinh}^2\left(B\sqrt{\frac{2{\theta}^2}{\theta_{\mathrm{max}}^2}-\frac{\theta^4}{\theta_{\mathrm{max}}^4}}\right).\end{align}$$ (8)

For B = 3, 4, maximum gain is seen at θ = θ_max, which is quite expected. For B = 5, 6, the radiation collapses and the noise spectrum approaches white noise, with the noise power increasing by many orders of magnitude, even outside the instability region, that is, at θ > $\sqrt{2}$θ_max. The intensity distribution in the beam is illustrated in Figure 2(b) (for θ_max = 1 mrad, В = 5) and in Figure 3(b) (for θ_max = 3 mrad, В = 6). The collapse is seen in Figures 2(b) and 3(b). The radiation is focused into a large number of filaments of huge intensity: I_max/I₀ = 274 for θ_max = 1 mrad and I_max/I₀ = 187 for θ_max = 3 mrad. In this case p = 0.28 and 0.45, where p = P_noise /P₀, that is, 28% and 45% of the beam power has been converted into noise, respectively. The power in each filament is P_fil ≈ 2.25P_cr, where the following applies:

$$\begin{align}{P}_{\mathrm{cr}}=\frac{0.174{\lambda}^2}{n_2{n}_0}\end{align}$$ (9)

The effect of beam self-filtering is expected to manifest itself at θ_max > θ_thr. The noise spectra S(θ) for θ_max = 10 mrad are shown in Figures 4(a) and 5(a) for various values of the B-integral. Figure 4(a) corresponds to the case with self-filtering of the beam before it enters the nonlinear medium with a contrast of 10⁸, and Figure 5(a) corresponds to the case without self-filtering. In the latter case, the noise gain maximum is at θ = θ_max. There is no such gain maximum in the noise spectra after self-filtering, and the nonlinear stage looks qualitatively different, as is clearly demonstrated in Figure 4(b) for B = 11 (p = 0.32). It can be seen that, instead of a large number of bright dots, a random honeycomb structure has appeared. A typical cell size of this structure corresponds to λ/(n₀θ_thr). To the best of our knowledge, such an instability has never been observed before either in numerical simulation or in experiment. We will refer to it as honeycomb instability. The honeycomb structure is preserved with an increase in the B-integral, the intensity in the walls of the honeycombs increases and their thickness decreases until it becomes less than the cell size of the computational grid. For comparison, the intensity distribution without self-filtering is shown in Figure 5(b) for B = 6, that is, for a value significantly smaller than that in Figure 4(b). Here, 41% of the beam power was converted into noise (p = 0.41), and a classical collapse of the radiation into a large number of bright points with a typical transverse scale λ/(n₀θ_max) is visible. There are two videos of SSSF for θ_max = 10 mrad (the time scale corresponds to the beam propagation in a nonlinear medium and, accordingly, to the growth of the B-integral) in the supplementary materials for this paper. In one of the videos (Silin_supplementmovie1.avi), the calculation was performed taking into account beam self-filtering, which leads to the development of honeycomb instability. In the second video (Silin_supplementmovie2.avi), there is no beam self-filtering, and filamentation instability develops.

Now consider the case of θ_max = 30 mrad that is typical of femtosecond pulses. The noise spectra S(θ) for different values of the B-integral are shown in Figure 6(a). Despite the beam self-filtering one can see gain maxima at θ = θ_max. We attribute this to the fact that θ_max = 30 mrad is far from the filtering threshold θ_thr = 4 mrad; therefore, the amplification of the noise components at θ < 4 mrad is very low at the linear stage of the SSSF, and the noise components in the θ∼θ_max region are in an advantageous situation despite the self-filtering. This is what differentiates the case with θ_max = 30 mrad from the case with θ_max = 10 mrad, where the components with θ < θ_thr are advantageous, so the honeycomb instability develops faster than the filamentation one. The intensity distribution in the beam for B = 16 is shown in Figure 6(b). It can be seen that the beam collapsed into a large number of random points and there is no honeycomb structure (here p = 0.51). For the honeycomb structure to appear in the beam cross-section, the development of filamentation instability should be slowed down even more. The numerical simulations have shown that this can be achieved by increasing the noise filter contrast by up to 20 orders of magnitude or more. Naturally, such a level of noise suppression cannot be achieved in practice, even with spatial self-filtering. However, to demonstrate the physical aspect of the problem, we have performed the corresponding numerical simulations.

The noise spectra for θ_max = 30 mrad at a noise filter contrast of the 24th order are plotted in Figure 7(a), and the corresponding intensity distribution in the beam for B = 27 is presented in Figure 7(b) (here p = 0.59). In this case, the filamentation instability is inferior to the honeycomb instability. A honeycomb structure with characteristic size λ/(n₀θ_thr) is well seen in the beam.

At the end of this section, we will consider a benchmark to provide some verification of our numerical simulation. Let us calculate the noise gain K(θ) in the linear mode (there is no depletion of the principal wave). In this case, the gain must coincide with the analytical formula (Equation (8)). Letting θ_max = 10 mrad and the filter contrast be 10⁸, we observed the appearance of a honeycomb structure with these parameters (see Figure 4). However, the input noise level we used produces a discrepancy between the gain K(θ) obtained from the numerical simulation and Equation (8) even at B = 1. Therefore, we have reduced the input noise level by two orders of amplitude. In this case, numerical simulation also leads to the appearance of a honeycomb structure; however, at small values of B, the gain coincides with Equation (8) with good accuracy. This can be seen in Figure 8, which shows the noise gain for two values of the B-integral, B = 5 and B = 8, obtained from numerical simulation and Equation (8).

Beam self-filtering allows suppressing the noise components possessing maximum gain at the SSSF, thus increasing the value of the B-integral, at which the beam is not fit for most applications, that is, an appreciable part of the energy goes into noise, and/or I_max becomes significantly larger than I₀. To understand how beam self-filtering affects the permissible values of the B-integral, we will construct various beam parameters as a function of B.

The share of radiation power converted into noise p = P_noise /P₀ versus the B-integral is shown in Figures 9(a) and 9(b). The maximum intensity in the beam normalized to mean intensity I_max/I₀ versus the B-integral is plotted in Figures 9(c) and 9(d). The root mean square (RMS) intensity in the beam where the intensity is normalized to the mean intensity RMS_I = RMS(I/I₀–1) as a function of the B-integral is shown in Figure 9(e). Finally, the RMS phase in the beam relative to the mean phase in the beam RMS_φ = RMS(φ−<φ>) as a function of the B-integral is plotted in Figure 9(f). Figures 9(a), 9(c), 9(e) and 9(f) show the curves for the input noise level of 0.02% (this noise level was used for obtaining the results described in Section 3), and Figures 9(b) and 9(d) correspond to the input noise level of 0.002%. The curves in Figure 9 correspond to the four values θ_max = 1, 3, 10 and 30 mrad. The curves for θ_max = 1, 3 mrad are plotted for the case of self-filtering. However, since θ_max< θ_thr = 4 mrad, the self-filtering does not affect the result, so the curves in the absence of self-filtering are the same. For θ_max = 10 mrad, the curves are plotted both for self-filtering and without it at a contrast of 10⁸. For θ_max = 30 mrad, the curves are plotted for two values of noise filter contrast, 10⁸ and 10²⁴, as well as without self-filtering. The effect of collapse is clearly demonstrated in Figures 9(c) and 9(d): at a definite value of the B-integral, I_max increases ‘almost vertically’, despite the logarithmic scale.

From the plots presented in Figure 9, we can draw the following conclusions. Without noise filtering, an increase in θ_max leads to an insignificant decrease in the parameters p, RMS_I and RMS_φ, while I_max remains almost unchanged. Thus, from the point of view of the I_max growth, the admissible B-integral does not depend on θ_max. Beam self-filtering at small values of the B-integral (B < 4) practically does not affect p or RMS_φ, while from the point of view of I_max and RMS_I, self-filtering suppresses the SSSF. This means that the amplified noise is phase noise.

With self-filtering, all curves in Figures 9(a)–9(d) are notably shifted to the right, which points to a significant increase in the permissible B-integral. The increase in the noise filter contrast increases the admissible B-integral still more. Another important fact is that the decrease in the input noise level affects the increase in the permissible B-integral in different ways. Without self-filtering, as well as for θ_max = 30 mrad and filter contrast 10⁸, when a honeycomb structure is not formed, a decrease in the input noise by a factor of 10 increases the admissible B-integral by approximately 1.5. At the same time, during self-filtering, when a honeycomb structure appears, the admissible B-integral with the same decrease in the input noise increases by about 2–3 for θ_max = 10 mrad and filter contrast 10⁸ and by about 6–14 for θ_max = 30 mrad and filter contrast 10²⁴. The permissible values of the B-integral at θ_max = 30 mrad are given in Table 1 for different values of filter contrast and input noise level under the specified criteria I_max/I₀ < 1.5, p < 10%. It can be seen from the table that without filtering, the criterion I_max/I₀ < 1.5 is more stringent than p < 10%. However, with an increase in contrast, the allowable values of the B-integral in terms of the criterion I_max/I₀ < 1.5 grow faster and, with a high filter contrast, the energy transfer from the principal beam to noise comes to the fore. Note that even with a contrast of 10⁶ only, the admissible values of the B-integral become two-digit.

The obtained results may be physically interpreted as follows. Without self-filtering, the instability is developing primarily at the spatial frequencies k_⊥max= k₀n₀θ_max and on the spatial scales w_max = 2π/k_⊥max = λ/(n₀θ_max) corresponding to the highest increment. From the known expression for the length of self-focusing (the length at which a beam having diameter w and power ${P}_1$ collapses)^[24^],

$$\begin{align}{L}_{\mathrm{sf}}=\frac{0.0925{k}_0{n}_0{w}^2}{\sqrt{{\left(\sqrt{P_1/{P}_{\mathrm{cr}}}-0.852\right)}^2-0.03}},\end{align}$$ (10)

Let us try to explain the absence of 2D self-focusing. Shortly before the development of filamentation or honeycomb instability (B is 2–3 radians less than in the case of collapse), the intensity distribution in the beam looks similar in both cases. An example of such a distribution is demonstrated in Figure 10. For the sake of generality, the color scale is not shown in Figure 10, as it may be different, but I_max here exceeds I₀ by several tens of percent. The characteristic scale w (characteristic transverse sizes of the spots) corresponds to the scale w_thr if θ_max > θ_thr (at a sufficient level of self-filtering contrast) and w_max if θ_max < θ_thr, or there is no self-filtering. Here, the conditions used to obtain Equation (10) are not fully realized, since there are no isolated round spots for which the self-focusing length could be calculated. One can see in Figure 10 compact bright (red) spots, but they are located on the elongated spots with a lower intensity. Therefore, the nonlinear phase gradient on the bright spots varies in different directions, resulting in different self-focusing rates in different directions. Further SSSF development depends only on the ratio of the scales w and w_max (in other words, on the $\theta_{\max} / \theta_{\mathrm{thr}}$ ratio). If there is no self-filtering (or θ_max < θ_thr), then on the w = w_max scale the 2D self-focusing increment significantly exceeds the 1D self-focusing increment, as a result of which the elongated spots break up into separate round ones and classical 2D self-focusing occurs with the formation of filaments. At a high enough level of contrast, self-filtering (θ_max > θ_thr) results in the characteristic scale w = w_thr > w_max on which 2D self-focusing does not have such an advantage over 1D self-focusing. So, due to a larger nonlinear phase gradient in the transverse direction of the elongated spots than along them, 1D self-focusing is faster than 2D self-focusing. As a result, the elongated spots turn into honeycomb borders.

The honeycomb patterns in nonlinear optics were reported in 2002^[25^]. There are fundamental differences in the honeycomb structure presented in Ref. [25] and in our paper. Firstly, somewhat different equations are solved. We consider the propagation in a medium with cubic nonlinearity (corresponding to the term in the NSE ~|A|²A), while in Ref. [25] the propagation in a medium with saturation of the linear part of susceptibility is considered (the corresponding term in the equation ~A/(1+|A/A_sat|²)). Secondly, and most importantly, is that in our work there occurs 1D SSSF into the honeycomb borders, so the intensity of the honeycomb borders becomes significantly higher than the intensity of the honeycomb interiors (see Figures 4(b) and 7(b)), while in Ref. [25], in contrast, the honeycomb borders are darker than the honeycomb interior.

It should be noted that honeycomb instability occurs over a wide range of initial parameters. The above figures present the results when the degree of PSD decay in Equation (5) is β = 1.55 and the spatial frequency filter has the form in Equation (7). However, numerical simulation has shown that the honeycomb structure in the beam appears for the β values at least in the 0 ≤ β ≤ 3 range, and the sharpness of the spatial frequency filter can vary in the range at least from $\exp \left(\hbox{--} {\theta}^2/{\theta}_{\mathrm{thr}}^2\right)$ to a step function. A necessary condition for honeycomb instability is sufficient filter contrast, which differs for different parameters. The parameters considered in the paper (β = 1.55, filter Equation (7) at θ_thr = 4 mrad, input noise power about 0.02% of the power of the main wave) require a filter contrast of at least 10⁶ for θ_max = 10 mrad, while for θ_max = 30 mrad a filter contrast of at least 10²⁰ is required. In any case, honeycomb instability exists (or does not exist) irrespective of the random realization of spatial noise used as a boundary condition for Equation (4).

We have considered a stationary instability model, taking into account pulse duration only in the filter transmittance (Equation (7)). Obviously, allowance for nonstationarity in the NSE (Equation (1)) will lead to pulse spreading in the course of propagation and, consequently, to a decrease in the effective B-integral; see, for example, Ref. [17]. This means that the values of admissible B-integrals listed in Table 1 will be even higher for fs pulses. A detailed investigation of this issue will be the subject of another publication.

The results of the performed numerical simulation lead to the following conclusions on the efficiency of beam self-filtering in terms of instability suppression.

(1) At θ_max > θ_thr, the classical filamentation instability leading to the appearance of bright points in the beam develops slower, that is, at larger values of the B-integral.

(2) With a high self-filtering contrast, the filamentation instability does not develop at all, since it is inferior to the earlier unknown honeycomb instability, resulting in a honeycomb structure in which all energy is accumulated within the boundaries between the cells.

(3) Beam self-filtering at θ_max > θ_thr significantly increases the admissible values of the B-integral, at which the beam retains acceptable quality: the maximum intensity exceeds the average one by no more than 50% and the fraction of power in the noise does not exceed 10%. The higher the filtering contrast, the larger the admissible values of the B-integral. In the absence of self-filtering, the mentioned intensity increase occurs earlier than 10% of the beam power converting into noise, whereas in the case of self-filtering with a sufficiently high contrast, the opposite is true.

(4) The reduction of the fraction of noise power in the input beam with self-filtering increases the admissible B-integral by a much larger value than in the absence of self-filtering.

(5) This study may be a route to explain experimental results^[17^] that cannot be explained by the existing theoretical concepts.