## Abstract

In focusing Kerr media, small-scale filamentation is the major obstacle to imaging at high light intensities. In this article, we experimentally and numerically demonstrate a method based on statistical averaging to reduce the detrimental effects of filamentation on the reconstructed images. The experiments are performed with femtosecond optical pulses propagating through a nonlinear liquid (toluene). We use digital holography to capture the transmitted optical image. The reverse propagation of the captured field is numerically performed using a numerical solution of the nonlinear Schrödinger equation. Because of their intrinsic sensitivity to measurement noise, filaments fail to propagate back on their initial trajectories and parasitic filaments form. The principle of the method is the introduction of artificial perturbations on the measurement, which spatially displace the parasitic filaments. By averaging the reconstruction over many realizations of the artificial perturbation, we show that the reconstruction improves the quality of the images. Finally, in order to identify the different regimes of optical power for which the filaments are time reversible, we also derive an analytical estimate for the condition number of the nonlinear propagator.

© 2015 Optical Society of America

## 1. Introduction

In recent years, researchers have displayed interest in optical imaging through nonlinear media [1–6]. Such activity is motivated by the potential improvement of optical resolution through mode mixing [5] and the possibility of extracting more information from noisy signals than what linear systems would allow [3]. In focusing Kerr media, small-scale filamentation [7] plays a dominant role. This phenomenon has been extensively studied [8–10] and some attempts have been made to transmit an optical image through Kerr media in the presence of small-scale filaments [1, 6].

In a previous work [6], we showed that small-scale filamentation is the main limitation for imaging through focusing Kerr media for two main reasons. The first one is the limitation of the experimental detection. In particular, imaging experiments are carried out using pulsed laser light, which is subject to significant material dispersion. As a result, measurements of the pulses in the spatio-temporal domain are required, something that is experimentally quite challenging. The second reason is the intrinsic sensitivity of the system to perturbations, which is the focus of the present work. As it will be shown below, even if we could acquire perfectly accurate holograms of the femtosecond pulses, experimental noise would still significantly impair the formation of images in the presence of filaments. Such sensitivity of filaments to noise has been recently studied and experimentally demonstrated in self-focusing bulk media [11, 12].

In this work, we show that the reconstruction error from a single noisy measurement can be reduced by adding digitally a spatial perturbation on the measurement and by averaging over several realizations of this perturbation. We will show that each perturbation induces a relative motion of the filaments with respect to the image. The collective effect of the filaments over different perturbation realizations is averaged out and leads to reconstructed images of better quality than the image reconstructed from the original noisy measurement with no digital perturbation. The optical power at which the experiment is performed is of crucial importance, as there are regimes where the effects of the filaments cannot be reduced by averaging. These power regimes can be systematically identified by analytically estimating the condition number of the nonlinear propagation operator. The experiments were performed with a femtosecond laser propagating in a nonlinear liquid.

## 2. Numerical method

Under the paraxial approximation, the wave dynamics of scalar optical beams in Kerr media, can be described by the nonlinear Schrödinger equation (NLSE) in physical units:

*,*$y$the transverse plane coordinates,$z$the propagation distance, ${\lambda}_{0}$ is the free space wavelength, ${k}_{0}=2\pi /{\lambda}_{0}$, $k={k}_{0}{n}_{0}$, where ${n}_{0}$is the background linear refractive index, and ${\tilde{n}}_{2}$ is the nonlinear refractive index (Kerr coefficient). Higher order nonlinear terms and dispersion can be added to the NLSE to refine the physical model, however, the simple form of Eq. (1) is sufficient to demonstrate the validity of the averaging method.

The procedure of propagation of an image through a nonlinear medium and its digital reconstruction (reverse propagation [13]) is schematically depicted in Figs. 1(a) and 1(b), respectively. In Fig. 1(a), a complex optical field${u}_{0}(x,y)$carrying the pictorial information about an object is incident on the input window of a nonlinear medium. The object field is transformed into an output field${u}_{\text{L}}(x,y)$, at a distance $z=\text{L}$, at the output window of the nonlinear medium through a physical nonlinear propagation operator${\widehat{F}}_{\text{L}}$_{.} The digital reconstruction of the object, i.e. the image, is performed by applying operator${F}_{\text{L}}^{-1}$to the output field,${F}_{\text{L}}$representing the digital model of ${\widehat{F}}_{\text{L}}$. The output field is complex-valued and includes both phase and intensity information. We use a non-paraxial split-step beam propagation method (SSF-BPM) [14] to perform the reconstruction. In the actual experiments, a field${u}_{\text{L}}+\delta {u}_{\text{L}}$is detected, where$\delta {u}_{\text{L}}$represents the experimental noise and any systematic error. In Fig. 1(b), an image${u}_{\text{0}}+\delta {u}_{0}$of the object is reconstructed from${u}_{\text{L}}+\delta {u}_{\text{L}}$by applying operator${F}_{\text{L}}^{-1}$, $\delta {u}_{0}$ representing the distortion in the reconstruction. When applying the averaging method described in Section 2.2, a digital perturbation ${\eta}_{\text{L}}$ is added to the output${u}_{\text{L}}+\delta {u}_{\text{L}}$. The induced contribution to the reconstruction, ${\eta}_{\text{0}}$, is simply defined from the equation ${F}_{\text{L}}^{-1}({u}_{\text{L}}+\delta {u}_{\text{L}}+{\eta}_{\text{L}})={u}_{0}+\delta {u}_{0}+{\eta}_{\text{0}}$.

#### 2.1 Sensitivity to noise and the condition number

In order to quantify how images are affected by experimental perturbations, we define a condition number for the nonlinear propagation operator *F*_{L} [15]:

*x-y*plane for each

*z*plane. The condition number of the corresponding linear operator is unity, since the linear operator is unitary. Equation (3) confirms the intuitive notion that the condition number increases in the high intensity regions of the field. Filaments form in the most intense regions in the field in a focusing Kerr medium and are responsible for the increase of the condition number. Note that ${\widehat{K}}_{\text{L}}$ is always monotonically growing and is generally a large overestimate of $\gamma $because it is calculated using the maximum intensity in the

*x-y*plane. This amounts to assume that, at each

*z*plane, the perturbation affecting the field is the one that is going to lead to the largest error in the reconstruction. However, this is very unlikely in reality, which is why the observed value of$\gamma $is generally much smaller than${\widehat{K}}_{\text{L}}$. It is also noteworthy that $\gamma $can decrease locally. Assuming a filament has formed at a distance$z={z}_{f}$, the intensity may drop again for$z>{z}_{f}$ if the filament is not sustained, which is generally the case for of small-scale filamentation (see for example [16]). The error on the reconstruction is the largest (large $\gamma $) when an experimental measurement is made in the $z={z}_{f}$plane, because the intensity is large. Besides being sensitive to the phase of the measured field directly, the filaments are also sensitive to the noise in intensity, because intensity translates into a phase shift through the Kerr effect. The reconstruction quality critically depends on the phase in the filament, hence a large value of $\gamma $. If the measurement is taken at $z>{z}_{f}$, in a plane where no high intensity features are present, that is not high enough to induce a significant nonlinear phase shift, the value of is $\gamma $generally smaller. During reverse propagation, the noise or perturbation originating from a plane $z>{z}_{f}$far from the filament tends to smooth out through diffraction before it reaches the plane$z={z}_{f}$. The induced disturbance on the filament is thus less, on average, than in the case the noise is added in the $z={z}_{f}$plane.

The relationship between the light intensity and the error can be observed both in experiments (Fig. 2(a)) and in simulations (Fig. 2(b)). The nonlinear image shown in Fig. 2(a) has been reconstructed while filaments were present in the measured field. In Fig. 2(b), we illustrated the sensitivity of the system with a dynamical simulation by adding noise on the simulated measurement. An initial object field was digitally propagated forward over a distance of 20mm. At each propagation step, a virtual measurement was taken and random noise was added to it. The power of the noise was 10^{10} times smaller than the power of the measurement. Then, for each virtual measurement, a reverse propagation was carried out to reconstruct the object. Figure 2 shows both the noise growth factor$\gamma $(calculated according to Eq. (2)) and the largest field amplitude as a function of propagation distance. The simulation was performed 100 times in order to get a statistical maximum for$\gamma $. Each time the field locally collapses to a filament, the maximum of the field intensity increases and a corresponding surge in the value of$\gamma $is observed.

#### 2.2 The averaging method

In this paragraph, we describe a general method that reduces the effect of small-scale filaments on reconstructed images of an object seen through a self-focusing Kerr medium. The method relies on the observation that, for moderate power, most of the nonlinear distortion is spatially concentrated in and around the filaments. As we showed in the previous section, once the field has undergone filamentation, propagation becomes difficult to invert. Newly formed filaments, which were not present in the forward propagation, tend to form in the reverse propagation. However the spatial configuration of these parasitic filaments is sensitive to the noise on the measurement. Therefore, if an artificial numerical perturbation${\eta}_{\text{L}}$is added on top of the experimental noise (see Fig. 1(b)), the parasitic filaments can be spatially displaced. For many realizations of${\eta}_{\text{L}}$, the filament-induced distortion can be reduced by averaging over the reconstructions. The perturbation ${\eta}_{\text{L}}$ is a band-limited signal with random Fourier components within the bandwidth, with uniformly distributed amplitudes, and is thus a smooth fully differentiable function. Depending on the scale of the image and on the power, the bandwidth and the relative strength of ${\eta}_{\text{L}}$ with respect to the measurement can be empirically adjusted.

The fundamental assumption behind the averaging method is that the weak linear pictorial background can be separated from the filaments for moderate power. In the reverse propagation process, the bright regions in the picture lead to the formation of filaments, whereas the weaker pictorial background only undergoes a negligible nonlinear phase shift. In addition, large spatial frequencies escape from filamentation and the corresponding parts of the image spectrum get restored with minimal distortion.

The imaging process follows the next steps: 1- The object field propagates in the nonlinear medium and filaments are formed. 2- The output field is detected and is generally affected by detection noise. 3- In the reverse propagation, the large phase curvature in the filaments is perturbed by the noise and prevents the filament to retrace its path on the way back. As a result, the area corresponding to the filament in the reconstructed image receives less light and the energy is distributed in a ring around the affected region. Even more importantly, parasitic filaments that were not initially present in the forward propagation can form on the way back. The position of these parasitic filaments is affected by the noise. Averaging over multiple noise realizations smoothes out the effects of filamentation. We provide more details on the averaging method in Section 5.2.

## 3. Experimental method and results

#### 3.1 Method

The experimental apparatus is shown in Fig. 3. The light source is a Spectra-Physics Ti:Sapphire regenerative laser amplifier (TSA Spitfire) with a pulse duration of 150fs, a central wavelength of 800nm and a maximal pulse energy of 2mJ. In the experiments shown below, each image is recorded with a single laser pulse. The imaging detector is an Apogee Neo monochrome sCMOS camera with a resolution of 2560 by 2160 pixels, 6.5x6.5 microns each. The dynamic range is 14bits and the camera was cooled down to −30°C. As in a previous work, we measure the output field using off-axis digital holography [4, 6]. The experimental apparatus consists essentially of a Mach-Zehnder interferometer and is depicted in Fig. 3. A tunable neutral density filter was placed at the laser output to adjust the pulse energy. The beam is sampled with a beam splitter (BS1) and the actual pulse energy is monitored with an amplified InGaAs photodiode. The photodiode output voltage was first calibrated against the pulse energy measured using a Molectron EPM1000 energy meter, which was placed at the input window of the nonlinear medium for the calibration measurement. The output compressor of the regenerative amplifier was detuned to produce a positively chirped pulse of 1.5ps in order to attenuate the combined effect of self-phase modulation and dispersion. A beam splitter (BS2) reflects 10% of the light toward a pulse compressor (with negative group delay dispersion) that recreates a 150fs transform-limited pulse, which is used as a reference beam. The reference pulse duration was optimized using a frequency resolved optical gating (FROG) device. The reference beam is passed trough a spatial filter to clean the spatial modes and then sent to beam splitter BS3. A piston (delay line) was inserted in the reference path to match the optical path length of the other arm of the interferometer. The initial beam that is transmitted through BS2 is reduced by a 1/2-magnification telescope and sent onto an amplitude mask that served as the object. The object consists of an array of absorbing stripes and was placed against the input window of a 100mm silica cell. The cell was filled with toluene as a Kerr medium, with a nonlinear refractive index${n}_{2}=1.3\pm 0.6\times {10}^{-19}$m^{2}/W at 800nm and 110fs [17] and a linear refractive index${n}_{0}=1.486$ [18] at 800nm). The corresponding Kerr coefficient for the amplitude of the electric field is${\tilde{n}}_{2}=1.03\times {10}^{-21}$m^{2}/V^{2} and the third order electrical susceptibility is${\chi}^{(3)}=1.02\times {10}^{-21}$m^{2}/V^{2} for toluene. The critical power for self-focusing, as given in [7] is: ${P}_{\text{cr}}=\pi {(0.61)}^{2}{\lambda}_{0}^{2}/(8{n}_{0}{n}_{2})$MW. A measured pulse energy of 1μJ corresponds to an average intensity in the bright area of the object (e.g. bright stripes of object in Fig. 1(c) of $1.66\times {10}^{12}$ W/m^{2}. For toluene this amounts to a nonlinear change in the refractive index $\Delta {n}_{\text{NL}}={n}_{2}I=2.2\times {10}^{-7}$and a total power of $6.2{P}_{\text{cr}}$ for the object considered (total optical power transmitted through the object). We assume that the power and the nonlinear phase shift scales linearly with the pulse energy measured with the photodiode.

The output of the cell was imaged by a 4f system onto a plane 22cm away from the camera plane and was recombined with the reference through beam splitter BS3. The reference was set at an angle of 1.3° in order to record an off-axis hologram of the field. The extra 22cm propagation distance allowed for the image of the filaments to spatially spread in order to avoid saturation of the detector. This defocus was compensated numerically using digital propagation of the recorded fields. The reference pulse was 10 times shorter than the signal pulse and therefore interference only occurs when the two pulses overlap. This allowed us to record the hologram of a selected time slice of the signal pulse.

#### 3.2 Results

We recorded series of images with increasing pulse power ${E}_{p}$, from the linear regime (${E}_{p}<0.1$μJ, $\Delta {n}_{\text{NL}}<2.2\times {10}^{-8}$, $P<0.55{P}_{\text{cr}}$) to the filamentation regime (up to ${E}_{p}$ = 5.6μJ, $\Delta {n}_{\text{NL}}<1.2\times {10}^{-6}$,$P=31{P}_{\text{cr}}$). Examples of recorded field and reconstructions are shown in Fig. 4, for a linear experiment (Fig. 4(a) and 4(b)) and a nonlinear experiment (Fig. 4(c) and 4(d)) performed with a pulse energy of 5μJ. The linear (low power) experiment is used as a control experiment to highlight the effects of the nonlinearity. It is assumed that only the beam power changes through the use of the tunable neutral density filter, the other parts of the experimental apparatus remaining unchanged. The reconstruction from the control experiment, shown in Fig. 4(b), serves as a pseudo-object${u}_{0}$that is used to measure the reconstruction error in the nonlinear experiments. In order to measure the quality of the reconstruction with the averaging method, we define the following errors: ${\epsilon}_{\text{amp}}(\Delta )\equiv \left|\Vert {u}_{0}+\Delta \right|-\left|{u}_{0}\right|\Vert /\Vert {u}_{0}\Vert $, where${u}_{0}$is the pseudo-object and${u}_{0}+\Delta $the reconstruction.

As illustrated in Fig. 4, the distortion in the digital reconstruction is clearly related to a filament that has formed on one of the bright stripe. This distortion has consistently the same appearance in all experiment, i.e. a dark area surrounded by a brighter ring as is clearly visible in Fig. 4(d). This defect has a characteristic spatial frequency given by the maximum of the spatial modulation instability spectrum [18, 7], the maximum gain frequency being ${k}_{\text{max}}=\sqrt{2{n}_{0}{\pi}^{2}{n}_{2}I/{\lambda}_{0}^{2}}$. Diffracted fields affected by filamentation thus display a characteristic feature size of $2\pi /{k}_{\text{max}}$. In the experiment shown in Fig. 4, the average intensity of the object (bright stripe) is $I=8.3\times {10}^{12}$ W/m^{2}. Thus ${k}_{\text{max}}=2.3\times {10}^{4}$rad/m and the characteristic size of the pattern is 0.27mm. By taking the Fourier transform of the distorted reconstructions (e.g. images in Fig. 4(d), 5(c) and 5(e)), we see that there is a corresponding peak at a frequency of $2\pi /0.28$mm, which is in good agreement with the theoretical prediction above (see Fig. 7 in appendix, Section 5.3).

In Fig. 5, we show the effect of the averaging method. The object, presented in Fig. 5(a) produces an output with one filament that was formed in the brightest region (Fig. 5(b)). In Fig. 5(c) we show the digital reconstruction from the measured field. It is strongly distorted due to the filament as described above. We added a random spatially band-limited noise (Fig. 5(d)) to the measurement, rescaled the total field to the initial power and performed the digital reconstruction again. The filaments in this reconstruction moved spatially (Fig. 5(e)). As expected, the reconstruction error becomes larger after adding the perturbation. However, by averaging the reconstructed field over 50 realizations of the noise, we obtain a smoother reconstruction (Fig. 5(f)) that is closer to the original object. According to the definition of the error, we have:${\epsilon}_{\text{amp}}({\delta}_{0}+<{\eta}_{0}>)<{\epsilon}_{\text{amp}}({\delta}_{0})$.

In Fig. 6(a) we show the evolution of the amplitude error as a function of the number of iterations. The error for the field (i.e. including the phase information), gets reduced by a lesser amount than the error on the amplitude because fluctuations in the reconstructed phase have a significant contribution to the error. If we only consider absorbing objects, then the phase becomes irrelevant and we only care about reconstructing the amplitude of the object. There is an optimum in the strength of the perturbation${\eta}_{\text{L}}$that yields the best improvement in the reconstruction of the amplitude (the optimum being around $\Vert {\eta}_{\text{L}}\Vert =0.25\Vert {u}_{\text{L}}+\delta {u}_{\text{L}}\Vert $). The existence of an optimum seems natural in this case because both too small or too large of a perturbation would fail to reduce error. However, this points to the fact that the statistical and spectral properties of${\eta}_{\text{L}}$could be engineered to reduce the error further. In these reconstructions, the bandwidth of${\eta}_{\text{L}}$was kept constant and only its relative power was changed.

In Fig. 6(b), we show the evolution of the condition number of the average field. The condition number decreases and stabilizes very soon as we iterate, even though the condition number of individual instances of propagation may get large because of parasitic filaments.

## 4. Conclusion

In conclusion, we presented a method based on the statistical averaging of numerical reconstructions to reduce the detrimental effect of filaments on images reconstructed after having propagated through focusing Kerr media. The error was quantitatively reduced by averaging, both for the total field and the field amplitude. The existence of an optimum for the strength of the perturbation suggests that the perturbation can be engineered to reduce the reconstruction error further. A measure of the invertibility of optical filaments (the condition number) in nonlinear Kerr media was also derived and used to interpret the experimental results. A large condition number implies that a small perturbation has a large potential to grow. In focusing media, steep increases in the value of the condition number can be directly related to small-scale filaments. The actual growth is eventually limited by the finite power contained in the field.

## 5. Appendix

In this appendix, we provide a detailed derivation of the condition number for the nonlinear propagation based on the split-step solution to the nonlinear Schrödinger equation. We provide a more formal description of the averaging method and relate it to the nonlinear Huygens-Principle. Finally, we show figures of the spectra of the nonlinear reconstructions that display the type of defect characteristic of filament-altered reconstructions.

#### 5.1 Derivation of the condition number

We consider the nonlinear propagation operator${F}_{\text{L}}$that transforms a scalar optical field ${u}_{0}(x,y)$ (the object field) into another optical field${u}_{\text{L}}(x,y)$(the output field) at a distance L, along the propagation axis *z*. By definition, we have: ${u}_{\text{L}}={F}_{\text{L}}({u}_{0})$. Similarly, we define the perturbations $\delta {u}_{0}$ that lead to the perturbations $\delta {u}_{\text{L}}$through the relationship: ${u}_{L}+\delta {u}_{L}={F}_{L}({u}_{0}+\delta {u}_{0})$. The *condition number* ${K}_{\text{L}}$ of operator ${F}_{\text{L}}$ is defined as follows:

Operator ${F}_{\text{L}}$ is nonlinear but it can be linearized locally by using the Taylor expansion (by ignoring high order terms) to get:

where$J$is the Jacobian of ${F}_{\text{L}}$. By definition of $\delta {u}_{\text{L}}$ and, we obtain:An expression for the Jacobian operator$J$can be derived from the SSF-BPM algorithm used to solve Eq. (1) [15]. In particular the nonlinear operator${F}_{\text{L}}$can be written as${F}_{\text{L}}={e}^{M\text{L}}$, where$A\equiv i{(2k)}^{-1}\left({\partial}^{2}/\partial {x}^{2}+{\partial}^{2}/\partial {y}^{2}\right)$,$B\equiv i{k}_{0}{\tilde{n}}_{2}|u{|}^{2}$ and $M=A+B$. The formal solution of Eq. (1) can be written now in an operator form as: $u(x,y,z)={e}^{Mz}u(x,y,0)$. Campbell-Baker-Hausdorff relation allows us to split the action of the two operators on the input field to two separate actions, by assuming small propagation distances. On other words, we discretize to *N* (*m* = 1,…, *N*) steps of $\Delta {z}_{m}$ distances and in such small interval we can approximate the solution as ${u}_{m+1}(x,y,{z}_{m}+\Delta z)\cong {\text{e}}^{A\Delta {z}_{m}}{\text{e}}^{{B}_{m}\Delta {z}_{m}}{u}_{m}(x,y,{z}_{m})$. A single step of the BPM can be viewed as a nonlinear operator ${F}_{\Delta z}$ that transforms the field over a propagation distance $\Delta z$: $u({z}_{0}+\Delta z)={F}_{\Delta z}[u({z}_{0})]$. The propagation step can be decomposed into a cascade of a linear$A$ and a nonlinear operator$B$. As a result the Jacobian operator$J$over many discrete steps can be expressed as:

We can now derive an upper bound ${\widehat{K}}_{\text{L}}$ for the condition number ${K}_{\text{L}}$ defined in Eq. (4). We have the following:

We can finally derive the following expression for the upper bound of the condition number in self-focusing nonlinear media:

*z*axis. As the discretization step$\Delta z$tends to zero, the condition number tends asymptotically to a finite value. We note finally that the condition number increases monotonically during propagation.

#### 5.2 Complementary information about the averaging method

To mathematically quantify the results of such an averaging method, we define the relative reconstruction error as, ${\epsilon}_{\text{field}}(\Delta )\equiv \frac{\Vert \Delta \Vert}{\Vert {u}_{0}\Vert}$ where${u}_{0}$is the object,$\Delta $the distortion and the Euclidean norm is defined as $\Vert f\Vert =\sqrt{{\displaystyle {\int}_{-\infty}^{+\infty}|f(x,y){|}^{2}}dxdy}$. Similarly, if we are only interested in reconstructing the field amplitude, we define the following error as: ${\epsilon}_{\text{amp}}(\Delta )\equiv \frac{\Vert |{u}_{0}+\Delta |-\left|{u}_{0}\right|\Vert}{\Vert {u}_{0}\Vert}$. Formally, if a perturbation ${\eta}_{\text{L}}$ is added to a measurement ${u}_{\text{L}}$ affected by a noise $\delta {u}_{\text{L}}$, we have: $<{F}_{\text{L}}^{-1}({u}_{\text{L}}+\delta {u}_{\text{L}}+{\eta}_{\text{L}})>={u}_{0}+\delta {u}_{0}+<{\eta}_{0}>$, where the brackets represent the statistical average over several realizations of ${\eta}_{\text{L}}$, with $<{\eta}_{\text{L}}>=0$. If the statistical characteristics of ${\eta}_{\text{L}}$ are properly chosen, essentially the amplitude and the spatial bandwidth, the following inequality holds on average $\Vert \delta {u}_{0}+<{\eta}_{0}>\Vert <\Vert \delta {u}_{0}\Vert $. A similar inequality holds on average for the error on the amplitude, meaning that

$\Vert |{u}_{0}+\delta {u}_{0}+<{\eta}_{0}>|-\left|{u}_{0}\right|\Vert <\Vert |{u}_{0}+\delta {u}_{0}|-\left|{u}_{0}\right|\Vert $. Note that these two inequalities may not be satisfied simultaneously and that they depend on ${\eta}_{\text{L}}$. At this point we have to mention that a more systematic way of analyzing and understanding the effect of statistical averaging in the image quality, is to perform a stochastic analysis of NLSE [20–22]. In particular, it is possible by using perturbation theory on a nonlinear type of Huygens-Fresnel principle (integral representation of NLSE based on Duhamel’s principle) to derive an infinite set of nonlinear evolution equations for the statistical average of the field [21]. This method has been applied to model the noise propagation in optical fibers [22]. We instead follow a direct numerical approach, which is less demanding computationally.

#### 5.3 Spatial spectra of the nonlinear reconstructions

As light propagates through focusing Kerr media, energy is transferred between spatial frequencies through four-wave mixing. The growth factor during propagation of the different spatial frequencies can be calculated [23] and reaches a maximum for a characteristic frequency${k}_{\text{max}}$. Field that have propagated through focusing Kerr media thus display characteristic features the size of which is given by $2\pi /{k}_{\text{max}}$. The distribution of the filaments and the spatial characteristic of the defects induced in reconstructions from noisy measurement follow similar spatial spectra. In Fig. 7, we show two examples of nonlinear reconstructions affected by filaments and the corresponding Fourier transforms. One example is the reconstruction from the raw detected field. The other example is the reconstruction including an artificial perturbation on top of the detected field, as we do when we apply the averaging method. The amplitude of the spatial the spectra of the reconstructions display a ring the radius of which corresponds to${k}_{\text{max}}$. The theoretical values are represented by the dashed circles on both spectra. Adding the perturbation does not alter the spatial characteristic of the reconstruction but it spatially shifts the filament enough to allow them to average out over several realization of the perturbation.

## References and links

**1. **A. Bramati, W. Chinaglia, S. Minardi, and P. Di Trapani, “Reconstruction of blurred images by controlled formation of spatial solitons,” Opt. Lett. **26**(18), 1409–1411 (2001). [CrossRef] [PubMed]

**2. **C. Barsi, W. Wan, and J. W. Fleischer, “Imaging through nonlinear media using digital holography,” Nat. Photonics **3**(4), 211–215 (2009). [CrossRef]

**3. **D. V. Dylov and J. W. Fleischer, “Nonlinear self-filtering of noisy images via dynamical stochastic resonance,” Nat. Photonics **4**(5), 323–328 (2010). [CrossRef]

**4. **A. Goy and D. Psaltis, “Digital reverse propagation in focusing Kerr media,” Phys. Rev. A **83**(3), 031802 (2011). [CrossRef]

**5. **C. Barsi and J. W. Fleischer, “Nonlinear Abbe theory,” Nat. Photonics **7**(8), 639–643 (2013). [CrossRef]

**6. **A. Goy and D. Psaltis, “Imaging in focusing Kerr media using reverse propagation,” Phot. Res. **1**(2), 96–101 (2013). [CrossRef]

**7. **R. W. Boyd, *Nonlinear Optics* (Academic, 2008).

**8. **G. A. Askar’yan, “Self-trapping of optical beams,” Sov. Phys. JETP **15**, 1088–1090 (1962).

**9. **R. Y. Chiao, E. Garmire, and C. H. Townes, “Self-trapping of optical beams,” Phys. Rev. Lett. **13**(15), 479–482 (1964). [CrossRef]

**10. **R. W. Boyd, S. G. Lukishova, G. Svetlana, and Y. R. Shen, eds., *Self-Focusing: Past and Present* (Springer, 2009).

**11. **B. Shim, S. E. Schrauth, A. L. Gaeta, M. Klein, and G. Fibich, “Loss of phase of collapsing beams,” Phys. Rev. Lett. **108**(4), 043902 (2012). [CrossRef] [PubMed]

**12. **N. Berti, W. Ettoumi, J. Kasparian, and J.-P. Wolf, “Reversibility of laser filamentation,” Opt. Express **22**(17), 21061–21068 (2014). [CrossRef] [PubMed]

**13. **M. Tsang, D. Psaltis, and F. G. Omenetto, “Reverse propagation of femtosecond pulses in optical fibers,” Opt. Lett. **28**(20), 1873–1875 (2003). [CrossRef] [PubMed]

**14. **M. D. Feit and J. A. Fleck Jr., “Beam nonparaxiality, filament formation, and beam breakup in the self-focusing of optical beams,” J. Opt. Soc. Am. B **5**(3), 633–640 (1988). [CrossRef]

**15. **A. Quarteroni, R. Sacco, and F. Saleri, *Numerical Mathematics* (Springer, 2000).

**16. **M. Centurion, Y. Pu, M. Tsang, and D. Psaltis, “Dynamics of filament formation in a Kerr medium,” Phys. Rev. A **71**(6), 063811 (2005). [CrossRef]

**17. **S. Courisa, M. Renarda, O. Fauchera, B. Lavorela, R. Chauxa, E. Koudoumasb, and X. Michautb, “An experimental investigation of the nonlinear refractive index (n_{2}) of carbon disulfide and toluene by spectral shearing interferometry and z-scan techniques,” Chem. Phys. Lett. **369**(3-4), 318–324 (2003). [CrossRef]

**18. **S. Kedenburg, M. Vieweg, T. Gissibl, and H. Giessen, “Linear refractive index and absorption measurements of nonlinear optical liquids in the visible and near-infrared spectral region,” Opt. Mater. Express **2**(11), 1588–1611 (2012). [CrossRef]

**19. **C. Ziehmann, L. A. Smith, and J. Kurths, “The bootstrap and Lyapunov exponents in deterministic chaos,” Physica D **126**(1-2), 49–59 (1999). [CrossRef]

**20. **J. Garnier, “Statistical analysis of noise-induced multiple filamentation,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **73**(4), 046611 (2006). [CrossRef] [PubMed]

**21. **C. Cecchi-Pestellini, L. Barletti, S. Aiello, and A. Belleni-Morande, “Mathematical methods for photon transport in random media,” J. of Quant. Spec. and Rad. Trans. **65**(6), 835–851 (2000). [CrossRef]

**22. **L. Barletti and G. Busoni, “Deterministic effective equations for the propagation of expectation in noisy nonlinear optical fibers,” Math. Methods Appl. Sci. **33**(10), 1221–1227 (2010). [CrossRef]

**23. **V. I. Bespalov and V. I. Talanov, “Filamentary structure of light beams in nonlinear liquids,” ZhETF Pis’ma **3**(11), 471–476 (1966).