The broad spectral bandwidths, high peak powers, and (3+1)-d nature of intense femtosecond pulses result in a variety of complex and interesting propagation effects. These include self-channeling[1,2], temporal break-up of the pulse[3,4,5] and extreme spectral broadening-commonly called continuum generation[6,7]. To a large extent, these phenomena rely on the basic process of self-focusing, due to an index of refraction proportional to the intensity. However, other physical mechanisms turn out to be of equal importance when a broadband ultrashort pulse is considered. Material dispersion is typically one of the most important additional effects influencing the dynamics of an ultrashort pulse. Indeed, at moderate powers normal group velocity dispersion acts to arrest catastrophic collapse with the result of temporal splitting of the pulse into two[3,9,10,11]. However, recent observations of multiple splittings, reported in Refs. [4,5,12], indicate that sharp temporal features, small beam diameters, and broad bandwidths further require that space-time coupling, nonlinear shock, higher-order dispersion, and ionization effects be included in theoretical analysis aimed at quantitative comparison with experiments. The importance of these higher order terms is elucidated with “full-field” (amplitude and phase) measurements as provided by frequency-resolved optical gating (FROG) [13]. It also becomes evident that quantitative comparison between theory and experiment requires careful theoretical description of the measurement process, including propagation of the pulse from the nonlinear medium to a measuring device. Due to the broad spectrum and small spatial dimensions, the vacuum propagation modifies the spatio-temporal characteristics of the pulse. In the present communication we will illustrate this point by presenting quantitative analysis of our results regarding propagation of femtosecond pulses in a nonlinear, dispersive medium at powers several above the critical power for self focusing.

1. Model

We model the evolution of the complex envelope E(r, z, t) of the radially-symmetric field

with the following modified nonlinear SchrÖdinger equation [12,14]:

In Eq. 2, the time t is measured in the frame moving at the group velocity of the pulse. The transverse Laplacian in cylindrical coordinates ∇² = ∂ ²/∂r ² + (1/r)∂ / ∂r accounts for diffraction, while the second and third time derivatives describe group velocity dispersion (GVD) and third order dispersion (TOD). The temporal, longitudinal, and transverse coordinates are normalized to the characteristic pulse duration τ, the dispersion length l_D = 2τ ² /|k′′|, and the characteristic transverse length $l_{\perp }=\sqrt{l_D/2k,}$ respectively. In addition, ∊ ₃ = k′′′/(3k′′τ), and k = 2πn/λ, with n being the linear index of refraction at the central wavelength λ. The dispersion coefficients k′′ and k′′′ are the second and third derivatives of k with respect to frequency, evaluated at the central frequency ω ₀. Space-time focusing [15] is given by the time derivative of the Laplacian, while nonlinear shock [17] is given by the time derivative of the nonlinearity [g(|E|²)]. As shown in Eq. (2), both of these terms are proportional to ∊_ω = 1/ω ₀ τ. These two terms act together to shift energy towards the trailing edge of the pulse [14,15,16]. Due to the short duration of the pulse, it is important to account for both instantaneous and time-delayed Raman nonlinearities [18], such that

with

In Eq. (4), n ₂ is the nonlinear index of refraction, and α denotes the fractional amount of the nonlinearity due to the Raman effect.

Unless noted, in the following numerical solutions of Eq. (2) we use parameters typical of our recent experiments in fused silica [12]. The initial field is taken to be a hyperbolic secant in time and a real Gaussian in space, having intensity full-width at half maxiumum (FWHM) of 90 fs and 70 μm, respectively. Furthermore, the beam waist is located at the entrance face of the sample. The linear index of refraction is n = 1.45 at the center wavelength of λ = 0.8 μm. The nonlinear index of refraction is n ₂ = 2.5 × 10^-16 GW/cm², resulting in a self-focusing critical power of P _crit = (0.61λ)² π/(8nn ₂) = 2.6 MW. The GVD and TOD coefficients are k′′ = 360 fs²/cm, and k′′′ = 275 fs³/cm. For the Raman response of fused silica we use α = 0.15, τ_r = 50 fs, and ω_rτ_r = 4.2[18]. Equation (2) is solved using a symmetric split-step technique. During the linear part of each step, the time-dependent terms are solved in the frequency domain, while the transverse Laplacian is computed using a finite differences scheme. The time-dependent nonlinear part of the equation is also solved using finite differences.

For accurate comparison with our experiments, it is important to account for the propagation of the field from the output face of the nonlinear medium (fused silica) to the measurement apparatus. In our experiments after propagation in the fused silica, the entire field is allowed to diffract in air over 1.5 m. At this point an aperture selects the on-axis portion of the field for characterization by the second-harmonic generation (SHG) FROG apparatus[5,12]. Mathematically, we describe the linear propagation of the pulse from the exit face of the nonlinear medium (near field) to the measuring device (far field) by the equation:

where we are using longitudinal, transverse, and time scales introduced above for the nonlinear medium. The appearance of the linear refractive index n in the first term of Eq. (5) is due to the fact that the wave vector of the field in air (assumed to be vacuum) is smaller than in the medium by a factor of n. The amplitude of the radially-symmetric field at distance L from the nonlinear medium follows from the solution of Eq. (5) and is given by the equation

where J ₀ is the zeroth-order Bessel function. Equation (6) is evaluated numerically using the frequency domain representation E(r′,0,ω) of the field at the output of the nonlinear medium. Finally we note that the on-axis in the far field (r = 0, L → ∞), the amplitude of the pulse simplifies to the expression

This simplified expression tells us that the on-axis far field at time t can be approximated by the integral over the cross-section of the near field at the same local time. This result will be used below when discussing the physics of pulse splitting.

2. Results

 

Figure 1. (a) Animated surface plot of I(r, t) of a self-focusing femtosecond pulse in fused silica for an input power of 4 MW (583kB QuickTime movie). (b) Calculated intensity profile I(r, t) in the far field (z = 1.5 m) for the field shown in (a).

Download Full Size | PDF

The importance of detailed description of vacuum propagation is illustrated in Figure 1. The animated surface plot of Figure 1(a) shows the basic dynamics of pulse splitting at relatively small input powers as predicted by Eq. (2). The intensity I(r, t) is shown as a function of the propagation distance z (z runs from 0 to 30 mm). The spatial dimension is radially-symmetric, with position r = 0 at the center of the axis. The orientation of the figure is such that early times are at the back of the figure. The peak intensity at z = 0 cm is 72 GW/cm², corresponding to a peak power of 4 MW.

This input peak power is greater than the critical power, such that moderate self-focusing occurs about t = 0; however, it is not sufficient to cause pulse splitting inside the nonlinear medium. Nevertheless, splitting occurs as the pulse propagates in vacuum to the far field [Figure 1(b)]. As discussed above, the free propagation of the field from the output of the fused silica to the measuring apparatus involves the evaluation of Eq. (6). Figure 1(b) shows the result of this calculation of the intensity distribution I(r, t) at z = 1500 mm beyond the output of the 30 mm fused silica media. This far-field spatio-temporal shape of the pulse clearly demonstrates splitting. Numerical simulations show that for the above input intensity the pulse splitting becomes noticeable several centimeters after the exit face of the nonlinear medium.

We can understand the far-field splitting with reference to Eq. (7), which shows that the axial far-field amplitude can be approximated by the integral over the near-field spatial cross section at the same local time. As such, a maximum in the on-axis field can come either from a cross section with the largest values of the field, or from a cross section with possibly smaller values of the field, but nearly constant phase. This second possibility is quite generally created during the process of self-focusing, where the spatial phase written onto the field varies rapidly about the peak intensity (t = 0), but is more constant in the wings of the field. As a result, the differently-phased portions of the field about t = 0 in the near field may destructively interfere in the far field leading to a local minimum and the observed temporal pulse splitting of Figure 1(b).

 

Figure 2. (a) Animated surface plot of I(r, t) for input power of 4.7 MW (583kB QuickTime movie). (b) Calculated intensity profile I(r, t) in the far field for the field shown in (a). Early times are at the back of the figure.

Download Full Size | PDF

The animated surface plot of Figure 2(a) illustrates the basic dynamics of pulse splitting at moderate input powers. The peak intensity at z = 0 cm is 85 GW/cm², corresponding to a peak power of 4.7 MW. In Figure 2(a), we observe that an initially uniform Gaussian input self-focuses in both space and time before splitting into two separate pulses. The pulse splitting results when self-focusing moves off-axis energy towards the peak of the pulse, while positive dispersion acts to pull this energy away from t = 0. As this process continues, the peak intensity drops, stopping the collapse at t = 0. However, off-axis energy continues to focus at t ≠ 0 such that two pulses are resolved.

The temporal asymmetries seen in Figure 2(a) are predominantly the result of the space-time focusing and nonlinear shock terms of Eq. (2)[16,14]. Both of these terms act to move energy towards the back edge of the pulse, creating a trailing shock front that is seen around z = 19 mm in Figure 2(a). The resulting shock front causes extensive broadening on the blue side of the spectrum via self phase modulation. Although this shock formation initially results in a trailing pulse with higher intensity and shorter duration, this pulse spreads much faster than the longer and lower intensity leading pulse as seen for z > 20 mm in Figure 2(a). We also note that the Raman nonlinearity, which is included Eq. (2), counteracts the nonlinear shock and space-time focusing terms by transferring energy to the leading (red-shifted) pulse. However, this is a smaller effect.

 

Figure 3. (a) Axial intensity (black) and phase (red) measured with SHG-FROG. (b) Calculated axial intensity and phase (c) Measured axial spectrum of field shown in (a). The blue line is the square modulus of the Fourier transform the data of (a), while the red points are measured with a spectrometer. (d) Calculated axial spectrum. Frame (a) of this figure is linked to an animated graphic with sound (330 kB QuickTime movie with sound) that allows the reader to see and hear the frequency variations associated field.

Download Full Size | PDF

In the case of Figure 2(a), the transformation from the near to the far field can result in multiple splitting of the field as shown in Figure 2(b). While there were two pulses at the output of the nonlinear medium [z = 30 mm frame of Figure 2(a)], there are now three pulses near r = 0 in the far field. This is shown in more detail by the data of Figure 3. In this figure, we present a comparison of the measured and calculated axial field for the same conditions shown in Figure 2. Figure 3(a) is the measured axial (r = 0) intensity and phase of the complex envelope E(r = 0, z = 1500 mm,t). This data was acquired using the SHG-FROG technique under experimental conditions corresponding to those of the calculation of Figure 2. For comparison, Figure 3(b) is the calculated axial intensity and phase, demonstrating good agreement with the measurement. The measured and calculated frequency domain representation of the field is shown in Figure 3(c) and (d), respectively. In Figure 3(c) the blue line is the square modulus of the Fourier transform of the measured field of Figure 3(a), while the red points are an independently measured spectrum (using a 0.27 m spectrometer). The discrepancy seen on the short wavelength side of the spectrum is due to bandwidth limitations in the SHG-FROG measurement. Nonetheless, this data verifies the fidelity of the FROG measurement over several orders of magnitude.

 

Figure 4. (a) Intensity distribution I(r,t) after 30 mm of propagation in fused silica. The input power in this case is 5.5 MW, corresponding to an intensity of 100GW/cm². (b) Axial intensity corresponding to (a). (c) I(r,t) in the far field for z = 1.5 m beyond the exit face of the fused silica. (d) Axial intensity corresponding to (c)

Download Full Size | PDF

The broadening of the pulse spectrum (to ~ 150 nm) as seen in Figure 3(c) and (d) is the result of intensity induced modulation of the phase of the field. For reference, the near-transform-limited input field has a spectral width of ~ 45 nm at the 10^-3 level. As means of further illustrating the self-phase modulation that results in the spectral broadening, we present animated video with sound that is linked to Figure 3(a). In this multimedia clip the intensity (blue dots) and phase (black line) of Figure 3(a) are again shown; however, the pulse phase is additionally presented as a “sonogram.” The pitch of the sound one hears when playing the movie is proportional to the local frequency (indicated by the moving red point). This local, or instantaneous, frequency is defined as ω_inst = ω ₀ - ∂ϕ/∂t with ϕ(t) being the phase of the complex envelope. One will note that sharp negative slopes of the phase correspond to higher pitch and therefore correspond to blue shifts. On the other hand, positive slopes in the pulse phase correspond to red shifts as indicated by the lower pitch.

At yet higher input power the spatio-temporal evolution of the pulse is characterized by further broadening of its spectrum, especially on the blue side. The resulting temporal dynamics can be described as a far field coalescence of the pulse toward a single broad pulse as illustrated by Figure 4. Figure 4(a) is the calculated intensity profile I(r, t) at the output of 30 mm of fused silica for an input power if 5.5 MW. The axial intensity in this case is shown in Figure 4(b). Of significance here is the absence of repeated multiple splittings of the form seen in Figure 2. This implies that the original two pulses of Figure 2(a) do not simply create two additional daughter pulses of their own. In the far field, shown in Figure 4(c), we see the field beginning to coalesce towards a single dominant peak. Our measurements at this and higher powers show good agreement with these calculations[12].

The fact that the axial intensity is larger in the far field at times where minimal intensity is seen in the near field can be qualitatively explained by Eq. (7). It is from the smaller intensity region about +75 fs in the near field that the large trailing peak arises in the far field. As already described above, smaller nonlinear phase variations are written onto the less intense regions of the near field, resulting in reduced destructive interference in the far field at the same local time.