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:

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}}$.

 

Fig. 1 (a) An object field u₀ is projected onto the input window of a nonlinear medium and propagates through it. A field$u_{\text{L}}$is formed at the output window. The operator${\widehat{F}}_{\text{L}}^{}$ represents the physical nonlinear operator that transforms $u_{\text{0}}$into$u_{\text{L}}$. (b) The object is digitally reconstructed using operator$F_{\text{L}}^{-1}$, which represents the numerical model of${\widehat{F}}_{\text{L}}^{-1}$. Operator $F_{\text{L}}^{-1}$ is applied to the measurement$u_{\text{L}}+\delta u_{\text{L}}$and produces a reconstruction$u_{\text{0}}+\delta u_{\text{0}}$, where$\delta u_{\text{L}}$represents the experimental noise and$\delta u_{\text{0}}$the induced distortion in the reconstruction. Any other perturbation${\eta }_{\text{L}}$added on top of the experimental measurement gives rise to an additional term${\eta }_{\text{0}}$in the reconstruction.

Download Full Size | PDF

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]:

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¹⁰ 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.

 

Fig. 2 (a) Example of linear and nonlinear imaging experiments. The nonlinear experiment illustrates the effect of filaments on the reconstructed image. (b) Simulation illustrating the relationship between the occurrence of filaments in a focusing Kerr medium and the growth of perturbations with propagation. The image labeled “z = 0” shows the intensity of the object at z = 0mm. The other images show the intensity with increasing propagation distances in millimeters. Graph at the bottom: The dotted line shows the maximum intensity of the field in the corresponding z = d plane. The solid line shows the maximum (over 100 realizations of the noise) perturbation amplification factor $\gamma$ of the reconstruction from the field recorded at distance z = d. The vertical bars on the graph indicate the positions of the corresponding images. The surges in the curves are correlated with the occurrence of filaments.

Download Full Size | PDF

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²/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²/V² and the third order electrical susceptibility is${\chi }^{(3)}=1.02\times {10}^{-21}$m²/V² for toluene. The critical power for self-focusing, as given in [7] is: $P_{\text{cr}}=\pi {(0.61)}^2{\lambda }_0^2/(8n_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². 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.2P_{\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.

 

Fig. 3 Diagram of the experimental apparatus. The laser source was a regenerative amplifier the output compressor of which has been detuned in order to produced a temporally chirped pulse of 1.5ps. A variable neutral density filter (NDF) is placed at the laser output to adjust the power. Beam splitter BS1 sampled the laser beam and sends part of it to a photodiode (PD) to monitor the pulse energy. Beam splitter BS2 reflects 10% of the beam power to the compressor on the left (BG: blazed grating, VRR: vertical retroreflector). The distance between the gratings was adjusted to produce a transform-limited pulse of 150fs. This pulse served as a reference for the interferometric measurement. The reference is passed through a spatial filter (SF) to clean up the spatial profile. The pulse that was transmitted through BS2 is spatially reduced in order to reach the required intensity range for the nonlinear experiment. An absorptive object was placed at the input window of a silica cell filled with a nonlinear liquid (toluene). The output window of the cell was imaged by a 4f lens system in a plane 22cm far from the CCD detector plane. The image and reference beams are recombined through beam splitter BS3 before the CCD detector.

Download Full Size | PDF

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.55P_{\text{cr}}$) to the filamentation regime (up to $E_p$ = 5.6μJ, $\Delta n_{\text{NL}}<1.2\times {10}^{-6}$,$P=31P_{\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.

 

Fig. 4 Example of experimental image reconstructions. The gray levels represent the intensity of the field and are rescaled for better visualisation. (a) Intensity of the optical field measured at the output of the cell for the linear (low intensity) experiment. (b) Linear reconstruction from the field shown in (a) using linear back-propagation (Fresnel diffraction). (c) Output intensity of the nonlinear experiment. The arrow indicates the position of the filament. (d) Reconstruction form the field shown in (c) using the nonlinear reverse propagation.

Download Full Size | PDF

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{2n_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². 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).

 

Fig. 5 (a) Pseudo-object (obtained from the linear reconstruction of the low intensity control experiment). The gray levels represent the intensity of the field and are all mapped to the same scale, except the image in (a) which is a linearly scaled version of the linear image of the object. (b) Optical field intensity after 100mm propagation in a positive Kerr medium. The position of the filament is indicated by the arrow. (c) Nonlinear reconstruction from the output field (b). (d) One instance of the random perturbation${\eta }_{\text{L}}$used for averaging. In this example, the power of the perturbation is 40% of the power of the signal. (e) Digital reconstruction from one realization of the perturbation on the data shown in (b). (f) Intensity of the average reconstruction over 50 realizations of the perturbation.

Download Full Size | PDF

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.

 

Fig. 6 (a) Error in the reconstruction for the complex field as a function of the number of realizations. The horizontal solid line is the error of the unperturbed nonlinear reconstruction. The horizontal dashed line is the error of the linear reconstruction. The thick dotted line is the error for the average. The thin dashed-dotted line is the error (resp. condition number) of the current reconstruction. (b) Natural logarithm of the condition number. The horizontal line is the condition number of the unperturbed nonlinear reverse propagation. The dashed line is the condition number of the reverse propagation for a particular realization of the perturbation. The dotted line is the condition number of the average reverse propagation.

Download Full Size | PDF

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:

(4)$$K_{\text{L}}=\underset{\delta u_0\ne 0}{\sup }\gamma (\delta u_0)\text{,}\text{with}\gamma (\delta u_0)=\frac{\Vert \delta u_{\text{L}}\Vert /\Vert u_{\text{L}}\Vert }{\Vert \delta u_0\Vert /\Vert u_0\Vert },$$

where $\Vert \cdot \Vert$ denotes the Euclidean norm. The factor$\gamma$is a measure of the relative growth of the initial perturbation $\delta u_0$ through the system. The condition number is thus the largest value that can be reached by $\gamma$ as a function of $\delta u_0$. Since the power is a conserved quantity, we have that $\Vert u_0\Vert =\Vert u_{\text{L}}\Vert$, therefore the above equation takes the simpler form:

(5)$$K_{\text{L}}=\underset{\delta u_0\ne 0}{\sup }\frac{\Vert \delta u_{\text{L}}\Vert }{\Vert \delta u_0\Vert }.$$

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

(6)$$F_{\text{L}}(u_0+\delta u_0)\cong F_{\text{L}}(u_0)+J(u_0)\delta u_0,$$

where$J$is the Jacobian of $F_{\text{L}}$. By definition of $\delta u_{\text{L}}$ and, we obtain:

(7)$$\delta u_{\text{L}}=F_{\text{L}}(u_0+\delta u_0)-F_{\text{L}}(u_0)\approx J\delta u_0.$$

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:

(8)$$J=\underset{N\to \infty }{\lim }\prod _{m=1}^N{e}^{A\Delta z_m}{e}^{B_m\Delta z_m}.$$

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:

(9)$$\begin{array}{l}{K}_{\text{L}}=\underset{\delta u_0\ne 0}{\sup }\frac{\Vert \delta u_{\text{L}}\Vert }{\Vert \delta u_0\Vert }=\underset{\delta u_0\ne 0}{\sup }\frac{\Vert J\delta u_{\text{0}}\Vert }{\Vert \delta u_0\Vert }=\Vert J\Vert =\\ =\underset{N\to \infty }{\lim }\Vert \prod _{m=1}^N{e}^{M_m\Delta z_m}\Vert \le \underset{N\to \infty }{\lim }\prod _{m=1}^N\left(\Vert e^{M_m\Delta z_m}\Vert \right)=\prod _{m=1}^{\infty }\left(\underset{n}{\max }{\sigma }_n\right)={\widehat{K}}_{\text{L}},\end{array}$$

where the ${\sigma }_n$ are the singular values of the propagator $\exp (M_m\Delta z_m)$, and the norm of an operator $M$ is defined $\Vert M\Vert =\underset{x}{\sup }\Vert Mx\Vert /\Vert x\Vert$. The product of all the maximum singular values in each discrete step, gives us a first upper bound of the condition number. Of course, since singular values are associated to the order of magnitude of the growth of small perturbations and the calculation of the corresponding Lyapunov exponents [19], it would be more useful to derive another bound directly related to physical quantities of our problem. Therefore, by using the inequality relation of the norm of the product of two operators $\Vert M_1{M}_2\Vert \le \Vert M_1\Vert \stackrel{.}{\Vert M_2\Vert }$, the fact that $\Vert \exp (A_m\Delta z_m)\Vert =1$, and that$B_m$is diagonal because the nonlinear phase factor applied to $u_m$ depends only on the amplitude of$u_m$itself, we get:

(10)$$\begin{array}{l}{\widehat{K}}_{\text{L}}=\underset{N\to \infty }{\lim }\Vert {\displaystyle \prod _{m=1}^N{e}^{M_m\Delta z_m}}\Vert =\underset{N\to \infty }{\lim }\Vert {\displaystyle \prod _{m=1}^N{e}^{A\Delta z_m}{e}^{B_m\Delta z_m}}\Vert \le \underset{N\to \infty }{\lim }{\displaystyle \prod _{m=1}^N\left(\Vert e^{A\Delta z_m}\Vert \cdot \Vert e^{B_m\Delta z_m}\Vert \right)}=\\ =\underset{N\to \infty }{\lim }{\displaystyle \prod _{m=1}^N\left(\Vert e^{B_m\Delta z_m}\Vert \right)}\underset{N\to \infty }{=\lim }{\displaystyle \prod _{m=1}^N\left(\Vert e^{k_0{\tilde{n}}_2{\left|u_m\right|}^2\Delta z_m}\Vert \right)},\end{array}$$

and, as a result:

(11)$${\widehat{K}}_{\text{L}}=\exp \left[\underset{N\to \infty }{\lim }{\displaystyle {\sum }_{m=1}^N\Delta z_m{(k_0\Delta n_{\text{NL,}\text{max}})}_m}\right]=\exp \left[k_0{\displaystyle {\int }_0^{\text{L}}\Delta }{n}_{\text{NL,}\text{max}}dz\right].$$

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

(12)$$\ln ({\widehat{K}}_{\text{L}})=k_0{{\displaystyle \int }}_0^{\text{L}}\left|{\tilde{n}}_2\right|I_{\text{max}}dz=k_0{{\displaystyle \int }}_0^{\text{L}}\Delta n_{\text{NL,}\text{max}}dz,$$

where we have replaced $|u{|}^2$ by$I$and expressed the condition number in term of $\Delta n_{\text{NL},\text{max}}(z)$, the maximum nonlinearly-induced change in the refractive index undergone by the wave for each position along the 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.

 

Fig. 7 (a) Nonlinear reconstruction with no artificial perturbation. (b) Example of nonlinear reconstruction with a specific realization of the perturbation $\eta$. The arrows indicate the theoretical characteristic size of the filamentation pattern of 0.27mm for the experiment presented in the main article. (c) Amplitude of the Fourier transform of image (a). The dashed circle shows the spatial frequency corresponding to the characteristic defect size. (d) Amplitude of the Fourier transform of image (b).

Download Full Size | PDF

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₂) 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).