2.1

Nonlinear Schrödinger Equation model

The nonlinear Schrödinger equation is widely recognized as the workhorse model for numerical pulse propagation in transparent media, especially in the extensively used spectral region spanning from UV to near-IR wavelengths⁷.

The model equation which describes the evolution of the slowly varying electric field envelope E_Env along z follows:

Eq. (1) takes into account both linear (first term on the RHS) and nonlinear effects (second term on the RHS) that are acting on the electric field envelope E_Env. is the transverse Laplacian operator which is used to model diffraction, while D is the dispersive operator used to model the material chromatic dispersion. The operators , and account for space-time focusing and self-steepening of the laser pulse. The nonlinear source terms N_Kerr, N_NLL, N_Plasma on the RHS of Eq. (1) are associated with the optical Kerr effect, nonlinear losses from MPA and plasma related effects like plasma defocusing and absorption.

While the nonlinear Schrödinger equation is more than adequate for modeling the propagation of pulses of wavelengths spanning from the UV to the near-IR, we will show that it unable to properly describe the propagation of intense pulses in the mid-IR and long-IR regime where carrier-wave related effects become important.

2.2

Unidirectional Pulse Propagation Equation

The unidirectional pulse propagation equation (UPPE) model^{9, 10}, is essentially a unidirectional Maxwell propagator that fully resolves the underlying optical carrier-wave in the forward direction. Despite the higher computational requirements of such an approach, the resulting spectral solver shows close ideal parallel scaling on our in-house SGI UV2000 machine enabling propagation of radial symmetric simulations over multiple kilometers in a matter of a few hours. Results presented throughout this work are generated using the unidirectional pulse propagation equation (UPPE) model in both radial symmetry and xyzt resolved numerical grids.

Since UPPE is a fully spectral electric field resolving propagator, the native representation of the laser field is expressed in arrays representing spectral amplitudes which are functions of angular frequency ω, and transverse wavenumbers k_x, k_y. The spectral nature of the solver makes it possible to solve the linear propagation problem exactly. For a given medium permittivity ε(ω), the linear propagation properties are encoded in the propagation constants of the modal fields. In the case of a homogeneous medium these are plane waves, and the propagation constant depends on the angular frequency and the transverse wavenumbers:

The permittivity is complex, consequently both chromatic dispersion and absorption are accounted both for exactly. This is major advantage of the UPPE model since in our case studies the spectral windows are enormous and need to be accurately resolved.

The polarization and current density are given as functions of the electric field , thus the medium response is calculated from the history of the electric field independently at each spatial point of the numerical grid.

The exact system of equations describing the evolution of modal amplitudes along the z-axis for the forward and backward propagating field is given below. Eq.(3) describes the forward propagating wave

And Eqn.(4) describes the backward propagating wave

Equations (3) and (4) are the coupled unidirectional pulse propagation equations (UPPE)¹⁰, while and are the unit vectors that are transverse to the propagation direction.

Equations (3) and (4) are mutually coupled through the polarization and current density terms. In order to solve the system, both the forward and backward propagating electric field are required. At this point we basic approximation of the UPPE model is made, i.e., that an accurately enough medium response can be calculated by the forward portion only:

Hence the canonical form of the scalar now UPPE can be written as follows:

Where the exact linear propagator is

And the nonlinear response term is given by

This scalar UPPE model has been shown to be remarkably robust¹¹ against spontaneously generated backward generated waves. Additionally the UPPE model is ideal in describing the MWIR and LWIR USP pulse propagation regime in various gaseous media. This is best reflected by comparing the UPPE model to the Nonlinear Envelope equation (NEE) model described in the previous section. The failure of NEE is illustrated in Figure 1 where a dramatic departure between UPPE and NEE becomes evident after initial self-focusing has set in. The simulation is for a 100 fs duration 3.9 μm pulse in dry air. The combined action of the optical carrier shock and critical self-focusing act in tandem to limit the intensity growth in the UPPE case. On the other hand, the envelope model is, by default, unable to capture the basic physics driving the propagation in this regime, and will start to blow-up. There simply is no physics based justification to assume that NEE or any nonlinear envelope variant of NEE can properly capture the physics at these longer wavelengths where carrier-wave driven effects are.

Figure 1.

Comparison of the evolution of an identical initial 4 μm pulses via the UPPE and the NEE. Lower inset on the left contrasts the strong optical carrier shock with a sinusoidal waveform after ~ 5 m. Upper inset on the top right shows the eventual splitting of the NEE envelope due to weak group velocity dispersion with the pulse dispersing subsequently.

Figure 2 shows a “phase diagram” reflecting the operative physics when the wavelength increases from 800 nm to 8 μm. The vertical intensity axis indicates achievable peak intensities in typical filaments. There are two regions that exhibit fundamentally different physical behavior and consequently greatly impact the propagation characteristics of USPs. For combinations of intensity-wavelengths that are in the black shaded region we predict the classical critical self-focusing (critical collapse) behavior characteristic of 800 nm filaments. In this region, linear dispersion of the relevant species dominates over Kerr nonlinearity – In other words, the dispersion length is shorter than the characteristic nonlinear length scale which depends on the strength of the Kerr coefficient. Here physics of USP propagation is effectively captured by the classical Nonlinear Schrödinger equation (NLSE) and its many nonlinear envelope variants.

Figure 2.

Phase diagram showcasing the separation between nonlinear envelope (black region) and UPPE (blue region) behavior.

On the other hand, in the blue shaded region to the right which spans all wavelengths longer than about 3.5 μm for intensities typically seen in filaments, the physics changes fundamentally. Here, the characteristic nonlinear length scale is shorter than the dispersion length scale. Now a new singularity emerges in addition to the critical self-focusing collapse. The underlying optical carrier-wave no longer remains sinusoidal and starts to self-steepen generating an electromagnetic shock wave.

3.1

Simple model for atmospheric optical properties

In this work we employ the scalar UPPE for all simulations of USPs in the MWIR and LWIR. For this first approach we use a simple model for air in an attempt to showcase the basic physical mechanism driving USP in the mid-IR optical regime.

We model dry air taking into account chromatic dispersion based on the Sellmeier formula as parameterized by Peck and Reeder^{12, 13}. The real part of the linear susceptibility is depicted in Figure 3. The nonlinear refractive index is taking into account the instantaneous electronic response of n₂ = 9 x 10^-24 m²/W which was extrapolated by measurements¹⁴ from Milchberg’s group at 800 nm. χ⁵ and other higher order nonlinearities are not taken into account. Ionization was modeled using the Keldysh formulation with an effective power law of 6.7 x 10^-138 x I^8.5. We assume 2.5 x 10²⁵ particles per cubic meter, only accounting for single electron ionization of Oxygen with a gap of 12.03 eV.

Figure 3.

Real part of the linear susceptibility of air based on the classical Sellmeier equation^{12, 13}.

The initial wavepacket used in the simulation has a duration of 24 fs at FWHM, carrying 177 mJ of energy at a central wavelength of 4 μm. The initial beam has a waist of 1.5 cm and is collimated. Our choice of wavelength is motivated by the flatness of the dispersive curve and the relative high transparency of the atmosphere in this specific spectral region. In addition to that, the ionization losses are expected to be much lower than at shorter wavelengths, making this wavelength ideal to explore long distance mid-IR filamentation.

Our simulations predict a new filamentation paradigm that is in stark contrast to classical 800 nm filaments. Our results for our 4 μm – 24 fs - 1.5 cm pulse are depicted in detail in Figure 4. In the classical picture of 800 nm filamentation, the non-saturated nonlinearity of the optical Kerr effect leading to a spatio-temporal collapse is arrested by plasma absorption and defocusing. However in the mid-IR regime studied here, a carrier shock is formed that leads to collapse-accelerated self-steepening of the optical carrier-wave, well before ionization becomes relevant. The steep shock front generates an enormous supercontinuum that encompasses all even and odd harmonics. These being at shorter wavelengths, propagate as wavepackets away from the fundamental at a different group velocity. Group velocity dispersion, although weak, is critically important in regularizing the collapse in his way. The result is an intense light bullet that is sustained for extended propagation distances driven by multiple recurrent carrier wave shock steepening events, which in turn initiate bursts of spectrally broadened emissions at the back of the main pulse that act to temporally soften the shock front and arrest the collapse. Notably, this collapse arresting mechanism is spectral/temporal in contrast to the classical Kerr vs plasma filamentation paradigm that is mainly acting in the spatial domain.

Figure 4.

Showcase of the new filamentation paradigm found in the MWIR and LWIR regime in the atmosphere. Spatio-temporal representation of the electric field of the 4 μm pulse at (a) z = 10 m and (b) z = 38 m. Top insets show the peak intensity and beam waist vs propagation distance z, and the corresponding spectrum at each distance. The red vertical lines indicate the position in z that is depicted in the surface on the bottom.

The resulting mid-IR filament is much thicker (~1mm) than its 800 nm counterparts (~100μm), and is able to extend over tens of meters in the atmosphere in a nearly lossless fashion (~15% losses of total energy over 100 m). As can be seen in Figure 4., the leading edge of the electric field is remaining almost unchanged over 30 m of propagation, while higher harmonic burst are constantly emitted from the tail of the wavepacket.

4.1

Extended model for atmospheric optical properties: HITRAN dispersion and absorption

Moving from 4 μm pulse to 10 μm pulses, we need to extend and refine our material model for air, since dispersion and absorption lines are now close to the 3^rd and higher harmonics of 10 μm. In our effort to refine our model we now include many quantities routinely ignored in theoretical filamentation studies, by utilizing the extensive HITRAN¹⁵ database. We assume the composition of air as 78.1% N₂, 20.9% O₂, 0.9% Argon, 0.037% CO₂, 1.8 x 10^-4 % CH₄ and 10^-5 % CO. Only the most common isotope for each species is accounted for. In addition, water vapor is also included in our model for a modest 10% relative humidity (good optical atmospheric conditions). Temperature is 25 C. The linear optical properties of this system are calculated based on the formalism detailed in previous works¹⁶. The resulting real and imaginary parts of the complex linear susceptibility are depicted in Figure 5. (a) and (b) respectively for a wide spectral range. Notably, the black vs the other two lines only differ in the amount of water vapor concentration in the atmosphere which showcases the importance and impact even 10% of water vapor has on the linear optical properties of air.

Figure 5.

(a) Real and (b) imaginary parts of χ¹ of the atmosphere as modeled by the HITRAN database for various humidity values at a temperature of 25C and 1 atm. The black, blue, and red curves correspond to 0%, 10%, and 50% relative humidity. The case of 10% relative humidity represents a good compromise of ideal, but still realistic, atmospheric conditions for long range laser propagation.

4.2

Extended model for atmospheric optical properties: Many-body interactions and excitation induced dephasing

In addition to refining the atmospheric model for linear dispersion and absorption in the previous section, we have to also take into account many-body effects and excitation induced dephasing in this wavelength regime. As can be seen in Figure 6. the impact of many-body excitation induced dephasing is much more pronounced at 10 μm than at 4 μm. The overall impact on beam curvature (black curve) is greatly modified when compared to the ionization-free case (black dashed curve). At 10 μm we predict a strong suppression of the optical Kerr effect self-focusing lens due to many-body excitation induced dephasing, which should in principle smooth out the abrupt intensity spikes commonly observed in optical filaments.

Figure 6.

Impact of various physical effects on beam focusing for (a) 10 μm, (b) 4 μm pulses vs time. Red curve: diffraction. Blue curve: optical Kerr effect. Magenta curve: ionization due to many body excitation induced dephasing. Black curve: sum of the above effects. Black dashed curve: sum of the above without EID. Calculation is based on Argon gas.

While the full study of the many-body excitation induced dephasing for the all species found in the atmosphere is still an ongoing work, we estimate that that Argon is a close qualitative approximation of air in this regard^{17, 18}. The simplified many-body effects model can be described by the following equation

Where f is the spatio-temporal plasma density due to many-body interactions, and E the electric field of the laser pulse. The parameters c = 2.84 x 10^-7 (mV^-4s^-1) and s = 4.6 x 10¹⁸ (V²m^-2) are fitted to the full microscopic results for Argon gas.

4.3

Simulation result of 10 μm pulses in realistically modeled atmosphere

For our study of 10 μm pulses we also included the delay rotational Raman response of Nitrogen and Oxygen as modeled by Palastro et. al.¹⁹. We split the optical Kerr effect into an instantaneous electronic response with nonlinear indices of n_2,Oxygen = 8.3 x 10^-24 m²/W and n_2,Nitrogen = 1.67 x 10^-23 m²/W, while the rotational response is modeled by using a damped harmonic oscillator with the natural frequency damping and steady state index chosen to fit the quantum mechanical density-matrix calculation¹⁹. Ionization is modeled using many-body effects as described in the previous section, while multi-photon ionization is considered negligible for our working intensities. Linear dispersion and absorption is modeled by using the full HITRAN dataset as described previously.

We use a 10 μm Gaussian spatio-temporal wavepacket with beam waist of 3 cm at 1/e² radius and a varying pulse duration from 2.4 ps to 3 ps. The wavepacket starts out collimated and without any input chirp, having a diffraction length of 283 m and peak intensity of 8.2 x 10¹⁴ W/m², which corresponds to a peak power of 1.16 TW (2.3 critical powers).

Figure 7. shows the result of our simulations over a 5 km propagation distance. We can see that for the optimal pulse duration of 3 ps, depicted as the black curves in Figure (a) and (b), an extremely long and uniform filament is created. This soliton-like optical structure is persisting over multiple Rayleigh ranges with an almost invariant beam waist and peak intensity. For shorter pulse durations, the filament becomes more intense and consequently less stable and overall shorter. This behavior is in line with the cumulative nature of the many-body ionization and excitation induced dephasing which is a cumulative effect in time. For shorter durations, the cumulative effect is not enough to suppress the Kerr lens and intensity rises to higher values, effective cutting the overall propagation distance short.

Figure 7.

(a) Peak intensity and (b) beam waist, versus propagation distance z for 10 μm pulses with initial intensities of 8.2 x 10¹⁴ W/m², 3 cm beam waist and 3 ps (black continuous curve), 2.7 ps (blue dashed curve), and 2.4 ps (red dotted curve) durations. The thin black dashed curve corresponds to the 3 ps pulse where many-body effects have not been taken into account. The black-dotted curve corresponds to the linear propagation of the 3 ps pulse. Figure taken from our recent published work²⁰.

4.4

Nonlinear HITRAN

Beyond the linear contribution from rovibrational response, nonlinear contributions are to be expected for high enough intensities. In order to include these into simulations, we applied the optical Bloch equations to transition data from the HITRAN database²¹. Since the amount of data in HITRAN is rather prohibitive, a selection procedure based on transition strength and detuning from the carrier frequency was applied, reducing the number of transitions from more than 140,000 to 12,709. The linear rovibrational response is dominated by rotational transitions below the carrier frequency, leading to a negative response overall. This is enforced by the nonlinear response to short pulses at increased intensities, with a reversal at very high intensities. The negative contribution to the polarization response is weaker for longer pulses. The overall properties of the polarization response indicate a weakened focusing due to rovibrational response.

For the purpose of propagation simulations, a simplified approached is applied within the UPPE framework. Here, correlations between transitions are neglected in order to speed up the response calculation. This approximation cannot describe coherent excitations that become relevant at high intensities and limits the applicability of this approach to intensities below 4×10¹⁶ Wm^-2. The propagation results for pulses with a power of 1.6 TW (corresponding to 1×10¹⁶ Wm^-2 initial peak intensity and slightly more than one critical power) and 3.2 TW (2×10¹⁶ Wm^-2) are shown in Fig. N for a model atmosphere of argon with 1% water vapor. A comparison to a pure Argon atmosphere is given as well. It can be seen, that the presence of water and the resulting nonlinear rovibrational response lowers the focusing of the pulse compared to the water free case.

Figure 8.

Peak intensity for a pulse with a power of 1.6 TW (black) and 3.2 TW (red) in an Argon-water model atmosphere and pure Argon for comparison. The 3.2 TW pulse is cut off when the peak intensity goes above the range of applicability for the simplified model.