1. Introduction

The dynamics of high-power, femtosecond (fs) light pulses in nonlinear media is a fascinating area of current optical research, and fs pulses continue to play a major role in the expanding area of ultrafast optical physics on ever shorter timescales. Moreover, the underlying nonlinear processes are prototypes of effects that occur in many other systems, far beyond optics. It is therefore of great importance to understand the fs pulse dynamics in various media in detail. Since virtually all measurements are greatly complicated by both the high intensities and the extreme temporal and spatial scales, it is desirable to extract all available information from each type of experimental data.

In most cases, experiments with fs pulses concentrate on “global” and “integrated” quantities, such as fluence (i.e., time-integrated intensity), total plasma production per unit of propagation length, fluorescence intensity from the produced plasma column, burn patterns in various targets, etc. Unfortunately, there exists no unique way to relate this type of data to the details of the complex dynamics in the fs pulse. There are methods that provide a direct “window” into the temporal dimension of the pulse dynamics [3]. Variants of FROG, for example, can measure the waveform of the pulse as a function of the local pulse time even for quite complicated pulses, and can be applied to self-focusing pulses [4]. However, in this case the transverse dimensions are integrated out. To mitigate this drawback, methods were introduced that can map the whole profile of the pulse intensity both in time and transverse coordinates (see e.g., [5, 6, 7, 8]). Clearly, these approaches are invaluable experimental tools, but they are relatively complicated at the same time.

Recently, far-field spectra were measured for fs pulses propagating through bulk condensed media [1, 9, 10]. In gaseous media, where intensities are a few orders of magnitude higher and the measurement correspondingly even more difficult, one could use new techniques that allow one to effectively switch off the nonlinearity in the measurement phase [11] , and thus to obtain spectrally resolved far fields. It has been recognized that relatively simple measurement of angle-resolved spectra can provide quite detailed information about the processes that control the fs pulse propagation [1].

The goal of the present paper is to extend, and describe in detail, an effective three-wave mixing (TWM) picture, previously employed to understand the dynamics of X-waves in water [2, 12], that allows one to interpret observed far-field spectra. The TWM picture not only allows one to understand the physical origin of all the characteristic features of the spectra, but also provides one with a means of relating features of the spectra with process occurring within the propagating pulse, that is, intra-pulse dynamics such as pulse splitting. Following the derivation of the effective TWM picture we discuss a detailed procedure that should aid in the experimental interpretation of far-field spectra. It is also demonstrated that the method applies to both normal and anomalous GVD regimes.

 

Fig. 1. Tabulated complex susceptibility of water is used in the definition of the linear spectral propagator of the UPPE simulator. The vertical lines mark the wavelengths at which simulations were performed.

Download Full Size | PDF

2. Effective three-wave mixing picture

In this Section we present the effective TWM picture for the spectral dynamics of fs pulses propagating in dispersive, bulk nonlinear media. Physically, the basic idea is to view the nonlinear optical response of the system as due to scattering of the incident field from a “material wave” representing the nonlinear change in the optical properties which is characterized solely by its propagation velocity. Together with the linear dispersion relation of the medium, this allows one to determine the region in the three-dimensional far-field spectral space (spanned by frequency and two transverse wavenumbers) where the energy of the scattered waves becomes concentrated due to phase-matching constraints on the TWM processes.

Much of the development of the TWM picture has been motivated by observations of numerical simulations. To reveal the rationale behind certain steps, we include illustrative simulation results in this Section. Specifically, we simulated a fs pulse propagating through a sample of water. At the input face of the water cell, the pulse is Gaussian both in space and time. The beam waist and focus were chosen 0.1 mm and 5 cm, respectively. The duration of the pulse was 94 fs (FWHM of intensity), and we concentrated on a range of total energies of up to few μJ. Below we present data obtained at the pulse energy of 0.6 μJ. The central wavelength was chosen 527 nm, i.e., in the normal GVD region. We used the Unidirectional Pulse Propagation Equation simulator [13] that can faithfully capture linear chromatic properties and frequency-dependent losses in the medium in an arbitrary wide frequency region. This capability is crucial for obtaining correct supercontinuum (SC) spectra. Figure 1 illustrates how our solver captures both linear absorption and chromatic dispersion of water. The left panel shows how the absorption of water is implemented, and the right panel compares the “simulated” GVD to the experimental data [14] shown as solid dots. The computational domain radius was 0.5 mm and the temporal length was 2 ps. Radial and temporal dimensions were sampled by 512 and 8192 points, respectively. We refer the reader to Ref. [15] for details of our model for water.

To start, the electric field of the propagating pulse is expressed in the scalar approximation as a superposition of plane-waves with appropriate weight coefficients that depend on the propagation distance

where the sum runs over all discrete transverse wavenumbers and frequencies corresponding to a large, finite normalization volume V = X × Y × T. In Eq. (1) and in the following, the plane-wave “propagation constant” in the z direction

depends on transverse wavenumbers k_x and k_y , on the angular frequency ω and on the frequency-dependent linear dielectric permittivity ε(ω).

The expansion coefficients A_u,v,ω (z) satisfy the Unidirectional Pulse Propagating Equation [13]:

In this equation, Δχ(x,y,z,t) represents the local modification of the material susceptibility due to all nonlinear processes combined. It will be simply called the ‘nonlinear response’ in the following. As a function of spatial and temporal coordinates, it usually exhibits several peaks that move together with the optical pulse that generates the nonlinear response.

 

Fig. 2. On-axis intensity as a function of propagation distance and local time in the frame moving with the pulse. This picture shows a pulse splitting event that generates two split-off daughter pulses. The latter propagate with slightly different group velocities that manifest themselves in different angles with respect to the horizontal axis of the bright, high-intensity locations.

Download Full Size | PDF

The relation between the temporal profile of the propagating fs pulse and the temporal profile of the nonlinear response is illustrated in Figs. 2 and 3. Figure 2 shows the on-axis intensity “map” I(0,0,z,t-z/v_g ) as a function of the propagation distance z and of local time τ=t -z/v_g in the frame moving with the group velocity of the initial pulse, and shows a pulse-splitting event that creates two daughter pulses that subsequently propagate with slightly different velocities until they decay. For longer propagation distances, similar pulse splitting can repeat several times. Figure 3 shows the corresponding on-axis nonlinear response Δχ(0,0,z,t -z/v_g ). There is a close similarity to the profile of the intensity in Fig. 2, but the nonlinear response also exhibits a region of negative change in susceptibility due to plasma defocusing that is represented by blue vertical band in Fig. 3.

In the following derivation of the TWM picture we use the idea from the numerical solutions that the nonlinear response may be viewed as a sum of localized peaks that move at different velocities. In particular, we shall analyze the case of two well localized peaks, these being shown as the bright pink regions in Fig. 3. The angle between the peak and the horizontal axis reflects the difference between the group velocity v_g for the incident pulse and the propagation velocity of the peak. The animation in the right panel of Fig. 3 illustrates the relative motion of the response peaks in the frame moving with the group velocity of the initial pulse.

 

Fig. 3. On-axis nonlinear response as a function of propagation distance and local time in the frame moving with the pulse. Bright areas represent two nonlinear response peaks generated by the two daughter pulses created in the pulse-splitting event. The angles with respect to horizontal of the bright areas correspond to the “propagation velocities” of the response peaks. The blue region represents the negative susceptibility modification caused by free electrons generated by the pulse. The right panel is an animation (1.5MB) illustrating the evolution of the two response peaks.

Download Full Size | PDF

To proceed with our derivation of the TWM picture we denote by v_r the velocity of the r-th peak in the nonlinear response, and approximate the total nonlinear change in susceptibility Δχ as a sum of contributions from the (two) main peaks:

This decomposition doesn’t require any approximations, but is of greatest utility when Δχ_r (x,y,z,t -z/v_r ) depends weakly on its third argument z. To determine the loci of the maximal spectral energy concentration, we decompose each response peak in Eq.(4) into Fourier components via

After inserting the field expansion (1) and the response representation (5) into the propagation equation (3), integrating over time and transverse coordinates, and grouping the exponential phase terms, one arrives at

The first two terms under the sum change much more slowly with z than the remaining phase terms collected in the exponential term on the second line. We denote by l_r the length-scale along the propagation distance over which the r-th response peak maintains an elevated value, and we neglect variation of Δχ _{kx-u,ky -v,ω-Ω}(z)A _u,v,Ω(z) over this length scale. The length l_r determines the length-scale over which nonlinear interactions are significant. The most significant contribution to the SC spectrum will come from the processes for which the fast oscillating phase terms in the exponential term on the second line cancel each other, and are thus phase-matched. However, due to the finite interaction length ≃ l_r , processes mismatched by ≈ 2π/l_r will contribute as well. Thus, the spectral energy should be concentrated in the phase-matched region determined by

This equation can be interpreted as a phase matching condition for the process in which a “material wave” with frequency ω-Ω and transverse wavenumbers k_x - u,k_y - v scatters an incident optical wave with frequency Ω and the wave-vector {u,v, k_z (Ω,u,v)} to produce a scattered or output wave of frequency ω and transverse wavenumbers k_x,y ,y. The right-hand-side of (7) then specifies with what accuracy this phase matching condition needs to be satisfied.

By taking account of the fact that even during the SC generation process most of the spectral energy still remains located around the incident central frequency Ω = ω ₀ and transverse wavenumbers u,v = 0, and this can be considered as a “source” of in-waves that take part in the scattering process, we obtain the simpler phase-matching condition

which defines a region in the space (ω, k_x ,k_y ) into which the response peak r transfers the optical energy of the original pulse. The shape of this region depends on the linear chromatic dispersion of the medium through the frequency dependence of the k_z terms on the left-hand-side. The central line(s) of the phase matched region is further determined by the deviation of the velocities of the peaks in the nonlinear response 1/v_r from 1/v_g , as we shall see in the next Section. The extent of the region from the central line is then given by how far the peak propagates (via the term 1/l_r ). Of course, since we have not used any information on the spatio-temporal spectrum of Δχ _{kx-u,ky -v,ω-Ω}(z) we can only describe qualitatively the shape of the low-value contours of the far-field spectrum.

In the next section we present the application of the TWM picture to SC generation in water in both the normal- and anomalous-GVD regimes. In particular, as a first step we treat v_r and l_r in Eq. (8) as parameters to be determined by fitting the phase-matched region(s) to the shape of the numerically generated far-field spectra. We then compare how the extracted parameter values compare with their counterparts obtained directly from the simulation data of Δχ. Consistency between the fitted parameters and those extracted from the simulation data then provides validation of the TWM picture, and the utility of fitting the spectra to extract information about the intra-pulse dynamics.

3. Far field SC spectra in water

3.1. Normal GVD regime

As a generic example of our simulations Fig. 4 shows the far-field spectrum of the pulse for the water parameters described in the previous section, and after a propagation distance of 2 cm. There are several typical features evident in the spectrum. First, there is a well developed central feature in the form of an X. This is the X-wave spectral feature that was shown to be responsible for the long-distance propagation of light filaments. In addition, there are two other branches of the spectral concentration at very low and very high frequencies and high transverse wavenumbers. These structures are usually not visible in the experimental spectra, because they encompass an extreme frequency region and exhibit low intensities. Finally, there is also significant energy concentration propagating axially, with low transverse wavenumbers (small off-axis angles). Next we show that all these features can be understood in detail using the TWM picture developed in the previous Section.

 

Fig. 4. Logarithmic far-field spectral power density of a fs pulse after propagating for 2 cm in water.

Download Full Size | PDF

We used Eq. (8) and a trial and error approach to obtain a visual best fit between the predicted phase matched spectral regions and the simulated far-field spectrum. The result, shown in Fig. 5, was obtained in two steps. First, the central lines of the regions are obtained, and subsequently the breath of each region is adjusted. Thus, the first parameters to fit are the velocities of the response peaks v _1,2, because they determine the central lines of the phase matched regions shown as blue lines. The peak velocities control the gap between the two disjoint components of central lines in each region, and the precise direction of the arms. Thus, one can use Eq. (8) with zero right-hand-side to fit the two component of the central line (shown in blue in Fig. 5) to match the position and directions of the spectral “arms.” To specify the central line location, a convenient quantity is the difference of the inverse peak and pulse group velocities, Δ$v_r^{-1}$ = $v_g^{-1}$ - $v_r^{-1}$. Our visual fit gave as the value Δ$v_1^{-1}$ = -1.5 × 10^-11s/m, for the leading response peak and Δ$v_2^{-1}$ = +2.5 × 10^-11s/m for the trailing response peak.

Having fixed the peak velocities, the propagation lengths l_r can be easily adjusted to match lines of equal intensity for the low spectral powers. To this end, we use Eq. (8) again, this time with v _1,2 determined as just described, and treat l_r as a free parameter. We chose to optimize the match of the resulting phase-matched region to the outer high-frequency boundary of the far-field spectrum, corresponding to the relative power density of ≃ 10^-10. Of course the value of l_r obtained will be sensitive to what is chosen as the “outer boundary” of the far-field spectrum. This is of course consistent with the fact that Eq. (8) is only determined up to a factor that depends on the shape of the spatio-temporal spectrum of the nonlinear response. For the contours shown in Fig. 5 we obtained values of l ₁ = 6.3 mm and l ₂ = 1.6 mm for the leading and trailing response peaks, respectively.

In the two panels in Fig. 5 we have plotted separately the phase-matched regions associated with the two different peaks r = 1,2 in the nonlinear response, and taken together we see that they match very well the shape and characteristic dimensions of the far field spectrum. Animation (Fig. 6) of the far-field spectrum evolution shows how with increasing propagation distance the spectrum fills up the union of the phase matched regions. In particular, the combined phase matched regions exhibit the characteristic X-shape around the origin that is associated with X-waves. Of course, the extremities in each region are not “populated” because of the finite bandwidth of the spatio-temporal spectrum of the nonlinear response that is not reflected in the TWM picture. In fact, the transverse extent of the spectrum is controlled by the transverse “dimension” of the nonlinear response peak as can be seen from Eq. (3). Thus these transverse dimensions can also be estimated from the far-field spectrum.

 

Fig. 5. Phase matched regions induced by the trailing (upper panel) and leading (lower panel) peaks of the nonlinear response. The central line of each region (shown in blue) marks the loci of best phase-matching and consists of two distinct components. Each response peak generates one half of the central X-shaped feature.

Download Full Size | PDF

An important point to note is that one can qualitatively predict which frequencies and propagating directions are predominantly created using the more conventional picture of four-wave mixing between monochromatic plane waves. However, since it doesn’t account for anything except the linear chromatic properties of the medium, it only works for the central X-like portion of the spectrum. Further, the shape and direction of the spectral X-arms can only be obtained approximately. Namely, the four-wave mixing predicts that the arms of X meet at a finite angle. One can see easily from Eq. (8) as well as from Fig. 5 that the arms are actually parabolas that touch at their vertices. Moreover the four-wave mixing can’t explain the extreme high- and low-frequency far-field spectral structures.

Naturally, the question arises, how reliable are the parameters that we have just extracted from the numerically generated far-field spectrum. We finish this sub-section by extracting the response peak properties directly from Δχ measured in the simulation for the sake of comparison. Figure 7 shows the locations and values of the response-peak on-axis maxima. The first panel shows that the propagation lengths of the response peak are indeed in the millimeter range corresponding to the values of l_r obtained from the three-wave mixing argument. However, it is clear that l_r is not a precise measure of the propagation length of the response peak. The reason for this is twofold. First, the phase matched region equation is determined up to a factor of the order of unity that re-scales l_r as already pointed out. This factor is given by the spectral power level contour we consider to be the “outside boundary” of the spectrum, as well as by the details of the nonlinear response spatial and temporal profile. Second, the response-peak velocities are not really constant. The second panel of Fig. 7 shows the peak arrival time in the local frame moving with the initial pulse group velocity which indicate that the velocities are largest just after the pulse splitting and then decrease. The parabolic fits of the form

 

Fig. 6. Animation (2.4MB) illustrating the evolution of the far-field spectrum filling-up the union of the phase-matched regions.

Download Full Size | PDF

 

Fig. 7. Maximum (left) and arrival time (right) in the local pulse frame of two nonlinear response peaks. Symbols represent maxima located in the data depicted in Fig. 3. Full lines show parabolic fits used to obtain the peak velocities and decelerations (see Eq. (9)).

Download Full Size | PDF

where z ₀ is the distance at which the response peak reaches its maximum, provide values Δ$v_1^{-1}$ = -2.0 × 10^-11s/m and Δ$v_2^{-1}$ = +2.5 × 10^-11s/m for the leading and trailing pulses, respectively. These values compare very well to those estimated from the TWM picture.

Furthermore, inspection of the nonlinear response in Fig. 7 reveals that the two peaks persist for distances l ₁ ≃ 7 mm and l ₂ ≃ 2.5 mm, again in reasonable agreement with the values inferred from the TWM picture. The above procedure for fitting the motion of the peaks to parabolas also provide us with “deceleration” parameter values of the order α ≈ 10^-9. It can then be shown (see Appendix A) that nonzero α effectively increases the phase-matching tolerance such that l_r is calculated from

This explains why our estimate of l_r from the TWM picture is actually smaller than the actual propagation distance of the response peaks. In fact, this frequency-dependent correction further improves the shape of the phase matched regions which are somewhat wider at high and low frequencies. This is because the corrected required phase-matching accuracy decreases with frequency difference from the central pulse frequency. However, it is likely to be difficult to measure the experimental spectra with an accuracy that would make it possible to use this correction to estimate also the response peak accelerations.

To summarize, we have shown that the TWM picture provides an excellent fit to the shape of the SC far-field spectrum using the velocities of the peaks v _1,2 and the scale length l_r as fitting parameters. Furthermore, using the numerical data for the nonlinear response we have shown that the inferred values for v _1,2 and l_r are consistent with the actual values. This both validates the TWM approach and shows that fitting far-field spectra, experimental or theoretical, using the TWM phase-matching conditions is a bona fide way of identifying these important pulse parameters.

3.2. Anomalous GVD regime

 

Fig. 8. Far field spectrum generated in water by a femtosecond pulse centered at 1100 nm, i.e., in the anomalous GVD region, and for a propagation distance of 1.2 cm. Elliptic structures in the long-wavelength portion of the spectrum around the central pulse frequency are characteristic of the anomalous chromatic dispersion. Note that the spectrum extends well into the normal GVD region.

Download Full Size | PDF

The effective three-wave mixing picture outlined in Section 2 should be applicable to any chromatic dispersion “landscape.” In the previous subsections the input pulses were centered around wavelength at which the medium exhibits normal GVD. Next we consider input wavelength in the anomalous GVD region.

We choose our new initial wavelength to be 1100 nm. This value is a compromise between the need to reach as deep as possible into the anomalous GVD region and the necessity to keep the linear absorption, which increases rapidly towards the long wavelengths (see Fig. 1), low. In order to compensate for increased absorption, we have increased the input energy of the pulse to 2.3 μJ, all other parameters being kept the same as before. It should be noted that inclusion of frequency-dependent linear absorption doesn’t qualitatively influence the structure of the far field spectrum in the nominally transparent wavelength region investigated. This is because the main effect of absorption in this region is to gradually decrease the pulse energy in the propagation length preceding the nonlinear focus, whereas the far-field spectral features mainly develop during the relatively short propagation length during which high intensities are maintained in the nonlinear focus. In particular, the effects of absorption are minimal over the length scale of the nonlinear focus, and hence have little influence on the qualitative structure of the far-field spectrum.

 

Fig. 9. On-axis nonlinear response as a function of propagation distance and local time in the frame moving with the pulse for the pulse with central frequency in the anomalous GVD region.

Download Full Size | PDF

Figure 8 shows the far-field spectrum after 1.2 cm propagation through the water cell. There are two dominant features: At high frequencies, we have a conical emission component similar to that observed in the normal GVD case. In fact, this is the part of the spectrum which does extend deep into the normal GVD region. The second characteristic feature of the far-field spectrum is the elliptic-shaped texture around the input frequency. This is the analogue of the X-shaped feature of the normal GVD supercontinua from the previous subsection, and a manifestation of the anomalous GVD. The elliptic shaped texture of the spectrum around the fundamental frequency reflects the nature of the so-called linear O-waves [16].

 

Fig. 10. Phase matched region induced by the leading peak of the nonlinear response. The central line shown in blue has two components. The elliptic component passes through the point corresponding to the carrier wave of the fundamental frequency.

Download Full Size | PDF

The nonlinear response behavior is similar to that observed in the normal GVD case, and is illustrated in Fig. 9. The self-focusing collapse is arrested by weak plasma generation and pulse splitting. In this case the pulse splitting is rather asymmetric, with the trailing response peak being split-off from the strong leading peak.

As in the normal GVD case, we can identify nonlinear material response peaks that are responsible for the spectral structure by fitting the phase matched regions obtained from Eq. (8) for the central frequency ω ₀ corresponding to the wavelength of 1100 nm. At the chosen propagation distance, the spectrum can be fitted by a single phase-matched region shown in Fig. 10.

The central line of the phase-matched region shown in blue has two components. The component close to the central pulse frequency is elliptic, as expected for the anomalous GVD regime. The high-frequency component represents the conical emission that actually propagates in the normal GVD regime. The parameter Δv ^-1 obtained from the visual fit is -7×10^-12s/m, and thus represents the leading nonlinear response peak. This is in keeping with Fig. 9 which shows that at z = 1.2cm the strong leading peak already decays while the weaker trailing peak continues to propagate. The value of Δv ^-1 obtained directly from the recorded nonlinear response data is -1×10^-11 s/m which is not far from our value extracted from the far-field spectrum. At later propagation distances, signatures from the trailing peak appear in the spectrum but are rather weak and overlap with the spectrum generated by the leading response peak.

4. Conclusion

We have presented an effective three-wave mixing picture for spectral dynamics of femtosecond pulses. Physically, this TWM picture is built on the idea that the far-field spectrum generated in the interaction of a fs optical pulse with a nonlinear dispersive medium results from scattering of the input radiation off the “material waves” that constitute the nonlinear medium response. The material-waves are generated mainly by the optical Kerr effect but also reflect the presence of higher-order nonlinearities due to generated plasma. Simulations also reveal that certain experimental parameters, such as transverse beam shape and focusing can affect the properties of the induced nonlinear response profiles and through that also the supercontinuum spectrum.

We have demonstrated that the TWM picture can explain all characteristic features that occur in the far-field spectra, including conical and axial emissions as well as conical radiation occurring at extreme low and high frequencies. Moreover, using an equation for effective phase matching, it is possible to obtain quantitative estimates for the relative velocities of the nonlinear response peaks by comparing the phase-matched region to the shape of the far-field spectrum. The good agreement between the simulated shapes of far-field spectra and the phase-matched regions obtained from the effective three-wave mixing argument strongly suggests that the spectra could be well approximated within the first Born approximation with the nonlinear response in the role of the scattering potential. Thus, a natural question to ask is whether or not it is possible to solve the corresponding inverse problem, namely to obtain from the measured far-field spectrum the whole two-dimensional map of the nonlinear response, similar to those shown in Figs. 3 and 9. This would be extremely useful, as the nonlinear response map actually provides more information than the spatio-temporal pulse intensity profile, because it reveals directly the presence of plasma.