ABSTRACT Title of Dissertation: HYDRODYNAMIC AND ELECTRODYNAMIC IMPLICATIONS OF OPTICAL FEMTOSECOND FILAMENTATION Nihal Jhajj, Doctor of Philosophy, 2017 Dissertation directed by: Professor Howard Milchberg Department of Physics The propagation of a high peak power femtosecond laser pulse through a dielectric medium results in filamentation, a strongly nonlinear regime characterized by a narrow, high intensity core surrounded by a lower intensity energy “reservoir” region. The structure can propagate over many core diameter-based Rayleigh ranges. When a pulse of sufficiently high power propagates through a medium, the medium response creates an intensity dependent lens, and the pulse begins to focus in a runaway process known as optical collapse. Collapse is invariably mitigated by an arrest mechanism, which becomes relevant as the pulse becomes increasingly intense. In air, collapse is arrested through plasma refraction when the pulse becomes intense enough to ionize the medium. Following arrest, the pulse begins to “filament” or self- guide. In gaseous media, energy deposited in the wake of filamentation eventually thermalizes prompting a neutral gas hydrodynamic response. The gas responds to a sudden localized pressure spike by launching a single cycle acoustic wave, leaving behind a heated, low density channel which gradually dissipates through thermal diffusion. This dissertation presents a fundamental advance in the theory of optical collapse arrest, which is how a pulse transitions from the optical collapse regime to the filamentation regime. We provide experimental evidence, supported by theory and numerical simulation that pulses undergoing collapse arrest in air generate spatiotemporal optical vortices (STOVs), a new and previously unobserved type of optical vortex with phase and energy circulation in a spatiotemporal plane. We argue that STOV generation is universal to filamentation, applicable to all collapsing beams, independent of the initial conditions of the pulse or the details of the collapse arrest mechanism. Once formed, STOVs are essential for mediating intrapulse energy flows. We also study the hydrodynamic response following filamentation, with the intent of engineering the response to construct a variety of neutral gas waveguides. In a proof-of-concept experiment, we demonstrate that a transverse array of filamenting pulses can be used to inscribe two distinct types of waveguides into the air: acoustic and thermal waveguides. These waveguides can be used to guide very high average power laser beams or as remote atmospheric collection lenses. HYDRODYNAMIC AND ELECTRODYNAMIC IMPLICATIONS OF OPTICAL FEMTOSECOND FILAMENTATION by Nihal Jhajj Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park, in partial fulfillment of the requirements for the degree of Doctor of Philosophy 2017 Advisory Committee: Professor Howard Milchberg, Chair Professor James Drake Professor Ki-Yong Kim Professor Daniel Lathrop Professor Phillip Sprangle © Copyright by Nihal Jhajj 2017 ii Acknowledgements I am deeply grateful for the opportunity provided to me by my advisor, Prof. Howard Milchberg, to work and learn in his group. Howard leads by example with his passion for physics and indefatiguable work ethic. I owe much to his guidance, and the confidence he has placed in me as a scientist. My good friend and colleague, Eric Rosenthal, was there with me from the beginning. He was often the first person I turned to when I needed help, staying up with me for long nights of data taking, and sharing in the ups and downs of being a grad student - his companionship made the everyday reality of working in the lab worthwhile. I must also thank Jared Wahlstrand for being the senior guy in the lab who showed me the ropes. Jared was always generous with his time when I came to him with questions. More than a few ideas in this dissertation are the direct result of brainstorming sessions I had with the “small lab crew”, including Eric and Jared, but also Sina Zahedpour, Arman Fallahkair, and our newest addition, Ilia Larkin. I’ve had the pleasure of working with Ilia over the past few years, and I’m happy to know that the future of the experiment is in capable hands. I’d also like to thank my other colleagues: John Palastro, Yu-Hsin Chen, Brian Layer , Sung Yoon, Andrew Goers, Jennifer Elle, George Hine, Linus Feder, Bo Miao, Fatholah Salehi, Daniel Woodbury, and Robert Schwartz for all the assistance and useful discussion. I’m grateful to the IREAP administrative and facilities staff; in particular Nolan Ballew and Jay Pyle for helping me out with the occasional machining project, and to Bryan Quinn for all the assistance in maintaining and upgrading the facilities in the lab. iii Finally, I’m deeply indebted to my girlfriend, Dina Genkina, for keeping me sane over these years, and to my mother, Inga Jhajj, father, Inderjit Jhajj, and sisters Tess and Tara Jhajj, without your love and support I would not have been able to start this process, much less finish it. I dedicate my thesis to you. iv Table of contents Acknowledgements ....................................................................................................... ii Table of contents .......................................................................................................... iv List of figures ............................................................................................................... vi Chapter 1: Introduction ................................................................................................. 1 1.1 Dissertation outline ............................................................................................. 1 1.2 Linear optics........................................................................................................ 3 1.2.1 Maxwell’s equations .................................................................................... 3 1.2.2 Linear polarization response ........................................................................ 4 1.2.3 Wave propagation in linear media ............................................................... 6 1.3 Nonlinear optics ................................................................................................ 14 1.3.1 The anharmonic oscillator .......................................................................... 14 1.3.2 Formalism for nonlinear optics .................................................................. 16 1.3.3 Optical Kerr effect ..................................................................................... 18 1.3.4 Ionization ................................................................................................... 21 1.3.5 Nonlinear polarization ............................................................................... 25 1.4 Filamentation of high power optical pulses ...................................................... 27 1.4.1 Filamentation summary ............................................................................. 27 1.4.2 Optical collapse and collapse arrest ........................................................... 28 1.4.3 Core/reservoir model of filamentation ....................................................... 31 1.4.4 Spectral broadening and self-steepening ................................................... 33 1.4.5 Modulational instability and multiple filamentation.................................. 34 1.4.6 Long timescale gas response to filamentation ........................................... 35 Chapter 2: Optical beam dynamics in a gas repetitively heated by femtosecond filaments ...................................................................................................................... 38 2.1 Introduction ....................................................................................................... 38 2.2 Experimental setup............................................................................................ 39 2.3 Results and discussion ...................................................................................... 41 2.3.1 Beam deflection ......................................................................................... 41 2.3.2 Density hole evolution ............................................................................... 46 2.4 Conclusion ........................................................................................................ 49 Chapter 3: Demonstration and analysis of long-lived high power optical waveguiding in air ............................................................................................................................ 50 3.1 Overview ........................................................................................................... 50 3.2 Direct imaging of the acoustic waves generated by femtosecond filaments in air ................................................................................................................................. 51 3.2.1 Introduction ................................................................................................ 51 3.2.2 Single mode filament acoustic response .................................................... 51 3.2.3 Multi-mode filament acoustic response ..................................................... 54 3.2.4 Identifying the guiding mechanism behind the single filament acoustic guide .................................................................................................................... 57 3.2.5 Conclusions ................................................................................................ 60 3.3 Generation of long-lived optical waveguides in air .......................................... 61 v 3.3.1 Background ................................................................................................ 61 3.3.2 Gas hydrodynamics initiated by femtosecond filaments ........................... 62 3.3.3 Experimental setup..................................................................................... 64 3.3.4 Multi-filament-induced guiding structure .................................................. 67 3.3.5 Fiber analysis of air waveguides ................................................................ 69 3.3.6 Injection and guiding experiments ............................................................. 70 3.3.7 Simulations of waveguide development and guiding ................................ 73 3.3.8 Conclusions and potential future application to guide high average power beams .................................................................................................................. 75 3.4 Modal analysis of air waveguides ..................................................................... 77 3.4.1 Introduction ................................................................................................ 77 3.4.2 Modelling ................................................................................................... 78 3.4.3 Modal structure of acoustic air waveguides ............................................... 80 3.4.4 Fundamental mode and leakage rate for thermal air waveguides .............. 82 3.4.5 Coupling efficiency .................................................................................... 86 3.4.6 Conclusions ................................................................................................ 87 Chapter 4: Spatiotemporal Optical Vortices ............................................................... 89 4.1 Overview ........................................................................................................... 89 4.2 Introduction ....................................................................................................... 89 4.3 STOV generation in a collapsing pulse ............................................................ 92 4.3.1 Modeling STOV generation ....................................................................... 92 4.3.2 Experiment ................................................................................................. 99 4.3.3 Experimental results................................................................................. 110 4.3.4 STOV dynamics and energy flow ............................................................ 116 4.3.5 Conclusions .............................................................................................. 120 4.4 Linear propagation and angular momentum of STOV-carrying pulses in material media ....................................................................................................... 120 4.4.1 Angular momentum in STOV-carrying pulses ........................................ 121 4.4.2 Free space STOV mode propagation ....................................................... 128 4.4.3 Conclusion ............................................................................................... 132 Chapter 5: Summary and future directions ............................................................... 133 5.1 Summary ......................................................................................................... 133 5.2 Acoustic manipulation of the filament-induced air density hole .................... 134 5.3 Scaling the thermal air waveguide to a ~50 meter span ................................. 136 5.3.1 Intra-beam phase variation over long filamenting distances ................... 137 5.3.2 Microphone array for single-shot filament visualization ......................... 139 5.4 Linear generation of STOVs using optical elements ...................................... 140 Appendices ................................................................................................................ 144 List of publications by the candidate ........................................................................ 165 Bibliography ............................................................................................................. 167 vi List of figures Figure 1.1. Susceptibility of a dielectric near a resonance. .......................................... 8 Figure 1.2. Schematic depiction of multiphoton ionization (left) and tunneling ionization (right). ........................................................................................................ 23 Figure 1.3. Left: schematic representation of dynamic self-channeling of air filament due to competition between self-focusing nonlinearity and plasma refraction. Right: end mode image of air filament. ................................................................................. 28 Figure 1.4. Simulation of the time evolution of a gas density hole created in the wake of a filamenting pulse taken from Ref. [9]. At early times < 0.3 µs, the high pressure volume originally occupied by the filamenting core expels mass and launches a single cycle acoustic wave, reaching pressure equilibrium by ~1 µs. ................................... 36 Figure 1.5. Log-log plot of density hole FWHM vs time. Measured FWHM vs time is shown in circles, with slope of ½ fitted to the experimentally derived values. Data is taken from earlier work published in Ref. [9]............................................................. 37 Figure 2.1. Experimental setup: the pump and probe travel collinearly through the gas cell. The probe beam is relay imaged at the axial center of the plasma where the density hole is produced.............................................................................................. 40 Figure 2.2. Two-dimensional plots of gas density in 1 atm xenon and air created by 1 kHz pulse trains with 0.12 mJ and 1.5 mJ per pulse respectively. The image was taken in between pulses in the pulse train. The deep center region was created largely by the pulse prior to the probe, while the larger, shallower background is due to the accumulation of heat from earlier pulses in the train. ................................................. 42 Figure 2.3. Top panel: Image of deflected white light supercontinuum spot as a function of pulse number in the 1 kHz pulse train. Lower panels: Hole parameters as a function of pulse number in the pulse train. The fast shutter is opened just before pulse 1. Equilibrium is approached by about the 250 th pulse (at 250 ms).The inset in the second panel is a typical relative density profile extracted from the interferogram. ..................................................................................................................................... 43 Figure 2.4. Far field beam deflection angle versus pulse number in 1 kHz plus train. Blue squares: measurement; red squares: ray optics calculation using measured gas density profiles. ........................................................................................................... 45 Figure 3.1. Acoustic gas dynamics following an ionizing pulse in air. The refractive index change n is shown at time T after the passage of the pulse. Blue: T = 80 ns; green: T = 280 ns; red: T = 480 ns. (a) Hydrocode simulation. The calculated n is plotted as a function of transverse coordinate. (b) Sequence of n profile lineouts obtained by interferometry of a 2 mm long plasma at 10 Hz. .................................... 53 vii Figure 3.2. Measured refractive index profiles induced by 60 J pump pulses at 1 kHz. The delays shown are with respect to passage of the pump pulse. The thermal gas density hole left between pulses is seen in the 60 s delay panel. ....................... 54 Figure 3.3. Measured and simulated air refractive index profiles induced by a TEM01 focus (vacuum peak lobe intensity 14 22.5 10 W cm  ) ((a)-(d)) and an * 04LG focus (vacuum peak lobe intensity 14 23 10 W cm  ) ((e)-(h)). The inset between (e) and (f) shows an expanded view of the focused * 04LG beam mode at low intensity. The simulations in (g) and (h) assume a continuous ring heat source at t=0 with a Gaussian cross section at the ring location. ................................................................ 56 Figure 3.4. Measurement of the axial dependence of energy deposition in a single filament and simulation of the propagation of a probe beam in the post-filament gas density profile. (a) Measured acoustic amplitude along the filament. The filamenting pulse propagates right-to-left. Error bars indicate the standard deviation of the signal amplitude over 100 shots at each position. The inset shows a sample acoustic trace measured by the microphone. (b-d) Beam propagation method simulations of a probe beam propagating left-to-right using hydrocode-generated gas index profiles weighted by the acoustic data of (a). The leftmost panels show index profiles associated with the maximum acoustic amplitude (near z = 42 cm) for the post-filament time delays in (b-d)............................................................................................................................. 58 Figure 3.5. Gas dynamics following a single filament in air. (a) Interferometric measurement of the refractive index change following a short pulse as a function of the time delay of the probe pulse. (b) Hydrodynamic simulation, assuming a 60 m FWHM Gaussian heat source of peak initial density 32 mJ/cm 3 . ............................... 63 Figure 3.6. Generation of a filament array using half pellicles. (a) A 55 fs, 800 nm, 10 Hz pulsed laser is used to generate an array of four filaments. A pulse propagates through two orthogonal half-pellicles, inducing π phase shifts on neighboring quadrants of the beam, and then are focused to produce a 4-filament with a TEM11 mode (actual low intensity image shown). A 7 ns, 532 nm 10 Hz pulsed laser counter- propagates through the filament and is imaged either directly onto a CCD for guiding experiments or through a folded wavefront interferometer and onto a CCD for interferometry. (b) Rayleigh scattering as a function of z with a bi-filament produced by a single half pellicle (the bi-filament far-field mode is shown in inset). The bottom row shows burn patterns produced by a 4-filament produced by two orthogonal half pellicles. ...................................................................................................................... 65 Figure 3.7. Interferometric measurement of the air density evolution induced by a 4- filament. (a) The acoustic waves generated by each filament cross in the middle, generating a positive index shift, producing the acoustic guide. (b) The acoustic waves propagate outward, leaving behind a density depression at the location of each filament. (c) The density depressions produce the thermal guide, with a higher central density surrounded by a moat of lower density. (d,e) The density depressions viii gradually fill in as the thermal energy dissipates. A movie of the 4-filament-induced gas evolution is provided in the supplementary material [81]. .................................. 69 Figure 3.8. Demonstration of guiding of 7 ns, =532 nm pulses in 70 cm long acoustic and thermal air waveguides produced by a 4-filament. The panel in the upper left shows the probe beam, which is imaged after the filamentation region, with and without the filament. The time delay of the probe was 200ns, which is in the acoustic guiding regime. The effect of the thermal waveguide, the shadow of which can be seen in the image in the top center (with a red dashed circle showing the position of the lower density moat), is shown in the bottom row, where the probe beam is imaged after the exit location of the air waveguide with and without the filamenting beam. The coupling efficiency vs. injected pulse delay is shown in the upper right. Peak energy guided was ~110 mJ. ....................................................................................... 71 Figure 3.9. Simulation of the evolution of and guiding in a thermal air waveguide. The top row shows the index of refraction shift produced by the 4-filament-induced temperature profile as a function of time. The bottom row shows a BPM simulation of the guided laser beam profile at the end of a 70 cm waveguide produced by the 4- filament-induced refractive index change. .................................................................. 74 Figure 3.10. Mode properties at =532 nm of single and octo-filament induced acoustic guides. Measured air refractive index profiles for single (a) and octo-filament excitation (c), along with corresponding hydrocode-simulations ((b) and (d)). Panel (e) shows an index profile simulation for the larger transverse scale typical of axially extended air waveguides. Calculated modes for the single filament-induced acoustic annular guide are shown in (f) p=0, m=0 and (g) p=0, m=5. The lowest order mode for the guide of (e) is shown in (h). Panel (i) shows the maximum azimuthal mode number mmax vs delay for the single filament-induced annular guide. Points: simulation, red curve: best fit...................................................................................... 81 Figure 3.11. Typical refractive index profiles and modes of thermal air waveguides. (a) Interferometrically measured refractive index profiles, (b) simulated time evolution of eight-lobe filament with 500 μm ring radius, (c) simulated time evolution of ring heat source with 500 μm ring radius and initial energy deposition equal to the eight-lobe filament case, (d) calculated fundamental leaky modes at =532 nm of the octo-lobe and ring air waveguides shown in (b) and (c). ............................................ 83 Figure 3.12. (a) Energy loss at =532 nm in dB/100 m for the octo-lobe index and ring index structures of Fig. 2 (blue and red respectively) and from the expression for  in Eq. (2) (black). (b) Loss in dB/100 m vs. time for a range of wavelengths (color scale in log base 10). ................................................................................................... 86 Figure 3.13. Coupling efficiency of a Gaussian beam at focus as a function of time and waist size for the (a) single filament acoustic guide, (b) octo-filament acoustic guide, and (c) octo-filament thermal guide. ................................................................ 87 ix Figure 4.1. (a) Phase and intensity projections of simulated pulse, in local coordinates, showing STOV generation. Propagation medium: atmospheric pressure air. Black dashed lines encircle vortices v1, v2 and v3, and arrows point in the direction of increasing phase. Our phase winding convention considers a (ξ ξ0) +i(rr0) winding about a null at (ξ0, r0) to denote a +1 STOV, giving the v1 STOV a +1 charge. The white dots on the intensity projection are centered on the locations of vortices v1, v2 and v3. (b) Simulation of an air filament crossing the air-helium boundary for the conditions of (a), showing the pulse fluence and plasma density. Nonlinear propagation terminates as the pulse transitions from air to helium, whereupon the pulse and a reference pulse (not shown, see Appendix A) are directed to an interferometer. .................................................................................................... 95 Figure 4.2. Toy model showing birth of a vortex pair via spatiotemporal phase front shear. The white curve and arrow depicts the axial (temporal) intensity profile and propagation direction, while the “core” and “periphery” labels denote the spatial intensity step. The z=0 panel shows the initial condition where phases are aligned, the z=zv panel shows the birth of the null (vortices overlap) and the z=2zv panel shows continued shear carrying the vortices apart, with the +1 vortex moving to the temporal front, and the -1 vortex moving backward. Our phase winding convention considers a (ξ ξ0) +i(xx0) winding about a null at (ξ0, x0) to denote a +1 STOV. In the third panel, red and black arrows indicate the direction of increasing phase. The two arrows connect the same lines of constant phase; the arrows show that the phase difference between head and tail is ill-defined due to the vorticity of the -1 STOV. .................. 97 Figure 4.3. Simulations of beam phase (top) and intensity (bottom) at the axial/temporal slice v where the STOV pair first appears. From left to right, the pulse is advancing along z, with the input shown at z=0, a collapsing beam at z=156 cm, ionization onset at 160 cm (not shown), just before the vortices spawn at 166 cm, just after the vortices spawn at 167 cm, and an image where the vortex pair has moved out of the =v plane. The linear intensity images are overlaid with centered lineouts of 10 0log ( )I I , where 13 2 0 10 W cmI   . Experimental parameters were used as code inputs: w0 =1.3 mm, pulse energy=2.8 mJ (P/ Pcr=6.4), pulse FWHM 45fs. ............. 99 Figure 4.4. Experimental setup for measuring the in-flight intensity and spatial phase profiles of collapsing and filamenting femtosecond pulses over a ~2m range. ........ 102 Figure 4.5. (a) Line emission (in arbitrary units) of the neutral helium 3 3 D – 2 3P , =587.6 nm transition induced by a tightly focused 800 nm pulse as a function of helium cell position. (b) Simulated filament intensity and (c) spatial phase at 800 nm (black) just before the ~4 mm air-helium transition, as well as reconstructed intensity and phase (red). The reconstruction is performed by propagating the solution 4 mm through the transition region followed by 50 cm of helium. This simulated far field is then back-propagated through 50.4 cm of vacuum to simulate reconstruction from imaging optics. .......................................................................................................... 105 x Figure 4.6. Spectral phase (radians, top) and spectral intensity (arbitrary units, bottom) taken from simulation at z=180 cm as the vortex core in the spatiospectral domain crosses 2.35 PHz (λ=800 nm). Simulation parameters are the same as Fig. 3 of the main text: w0 = 1.3 mm, pulse energy = 2.8 mJ (P/ Pcr = 6.4), pulse FWHM 45 fs. From left to right, the top row (bottom row) shows the full spatial-spectral phase (intensity), as well as spatial phase (intensity) slices spectrally fore and aft of 800 nm. The dashed black line in the top left panel intersects the vortex core. ..................... 107 Figure 4.7. Comparison of the spatiospectral phase of the pulse at 800 nm, to the weighted spatial phase measured about 800 nm using the same simulation parameters considered in Fig. 3: w0=1.3 mm, pulse energy=2.8 mJ (P/ Pcr=6.4), pulse FWHM 45 fs. Flipping of the spatial phase occurs at z=179 cm at 800 nm and z=185 cm for the weighted spatial phase. ............................................................................................. 109 Figure 4.8. Beam on-axis phase shift (with respect to flat-phase reference arm) as a function of pulse power at zh = 150 cm. The phase jumps abruptly by ~2 at P/Pcr ~ 5, providing a clear signature of the transition from the pre-collapse to post-collapse beam. The red points are averages over 2600 shots (blue points) in 150 energy bins. ................................................................................................................................... 112 Figure 4.9. Top row: retrieved spatial phase (in radians) and intensity (arbitrary units) images at z = 165cm, P/Pcr = 4.4 for i) a pre-collapse beam, ii) and iii) beams where a vortex ring is on either side of the reference central wavelength of 800 nm. The bottom row shows that as the maximal phase gradient in the images increases, the maximal phase shift saturates at π for all cases of P/Pcr leading to beam collapse. . 114 Figure 4.10. Demonstration of energy flow about a +1 “r-vortex” STOV. White arrows correspond to size and direction of j in Equation (4.7). There are three distinct regimes: the left panel shows the saddle regime which exists in regularly dispersive media (2>0), the middle panel shows the degenerate regime which can exist in regular or anomalously dispersive media and has a dominant axis for energy flow, and the right panel shows the spiral regime which exists in anomalously dispersive media (2<0). ............................................................................................................ 117 Figure 4.11. Simple model of the propagation effect of the envelope phase of a pulse with an embedded STOV. Color gradient indicates change in phase and central circular black arrow indicates the direction of phase winding. In regions a) and c) respectively, the phase winding blue and red shifts the pulse as indicated by the blue and red arrows (indicating momentum) and vertical lines (wavefronts). In regions b) and d) respectively, the phase winding deflects the local momentum density of the beam downward and upward as shown by the solid black arrows/wavefronts. Dashed black arrows indicate the momentum density and unperturbed wavefronts for the case of a constant envelope phase..................................................................................... 124 Figure 4.12. Evolution of pulse intensity of a quasi-mode (top) in normally dispersive media, and mode (bottom) in anomalously dispersive media. The pulse envelopes have been rescaled so that they have equal size at all z locations. Pulse intensity is xi shown in the y=0 plane. Phase winding due to vorticity is indicated with red arrows. ................................................................................................................................... 131 Figure 5.1. Filament deflection by 25 Hz acoustic source. Figure a), on the left, shows deflections with a speaker current of 4 A RMS, with figure b) showing deflections with a 10 A speaker current. ..................................................................................... 135 Figure 5.2. Third harmonic generation scheme employing a speaker and grating to periodically modulate the location of the density hole. ............................................ 136 Figure 5.3. Phase interferograms collected from interference of the two lobes of a TEM01 filamenting mode after 20 meters of nonlinear propagation. a) Plots of different input power (per lobe) are shown, where lineouts of single shot interferograms are shown in different colors. b) Blue dots plot the transverse spatial phase difference between the two lobes, and red dots show the standard deviation in the transverse spatial phase difference. ..................................................................... 139 Figure 5.4. Left: close up of individual microphones, each with local A/D converter. Right: full mounted array of 8 microphones alongside short filament (which is on the other side of the mounting board), showing data hub on the bottom right. .............. 140 Figure 5.5. Schematic for using a 4f pulse shaper and SLM or static masks to impose a line STOV on a beam. ............................................................................................ 143 Figure A.0.1. Schematic of light launched from end of single cycle annular acoustic guide. D is the size of the annulus, Δr is the thickness of the compressive section of the wave, Δθ is the spreading angle for light launched from the guide, and Δzint is the scale length of the peak on-axis interference from light launched by the acoustic guide. ......................................................................................................................... 150 Figure A.0.2. Depiction of approximation for moving from smoothly varying radial index profile (left) to stepwise profile (right). .......................................................... 155 1 Chapter 1: Introduction 1.1 Dissertation outline Filamentation is a dynamic nonlinear propagation regime for high power optical pulses, where competing focusing and defocusing nonlinearities allow for the extended propagation of a high intensity filamenting core over many Rayleigh lengths [1]. Interest in filamentation stems from the rich nonlinear physics of the system, and the many applications that this affords. The large optical field present in filaments enables efficient nonlinear conversion for electromagnetic sources spanning from THz to x-rays: optical filaments can be used to generate broadband coherent optical fields [2] which can be compressed down to nearly single cycle duration [3], they can be used as a source for high peak power single cycle THz pulses [4,5], and for high harmonic generation [6]. Filamentation in air is of particular interest for remote applications, where filaments are unique in their ability to deliver high field intensities (up to 14 2~10 W cm ) at kilometer distances, allowing for ranged applications such as LIDAR [7] and laser induced breakdown spectroscopy [8]. This dissertation presents experiment, simulation and theory with two central thrusts: i) new applications for filaments made possible by exploiting the long timescale gas evolution occurring in the wake of a filamenting pulse [9], in particular the construction of long-lived optical waveguides [10], and ii) the discovery of the spatiotemporal optical vortex (STOV), a new type of optical vortex self-generated in all filamenting beams [11]. We establish that STOV generation is intrinsic to the initiation of filamentation, and central to the phase evolution and intra-envelope energy transport occurring in the subsequent propagation. 2 The remainder of Chapter 1 provides a background for filamentation. Starting from Maxwell’s equations, linear optics is briefly reviewed with an emphasis on dispersive pulse propagation. Nonlinear processes central to air filamentation are discussed, followed by a short review of filamentary phenomena. In Chapter 2 we explore the effects of the density hole left in the wake of a filamenting beam on a 1 kHz pulse train of filaments, where it is seen that the cumulative effects of the pulse train lead to convective motion of the gas and a sustained deflection of the beam pointing. Chapter 3 develops and demonstrates the idea of acoustic and thermal air waveguides, where we show it is possible to guide high power (both peak and average) and high energy laser pulses in the wake of structured transverse arrays of filaments. Proof-of-concept ~1 meter waveguides are generated, the gas hydrodynamics leading to guide formation are time-resolved using interferometry, and modal analyses of the various guides are also performed. Chapter 4 presents experiment, simulation and theory establishing the universality of STOV generation in the optical collapse process. We interferometrically reconstruct the transverse intensity and phase profiles of an air filament in midflight to detect the presence of the STOV, and develop a model of optical collapse arrest generic to all filamentation processes, which predicts the generation of STOVs. Finally, we theoretically investigate properties of linear pulses with embedded STOVs, identifying a STOV- carrying mode as well as assessing the orbital angular momentum of a STOV- carrying pulse in a dispersive medium. Chapter 5 presents a summary of the thesis, reviews ongoing work, and suggests future directions. 3 1.2 Linear optics 1.2.1 Maxwell’s equations Maxwell’s Equations, together with the Lorentz force law govern classical light-matter interactions [12]. The macroscopic formulation, in SI units, of Maxwell’s Equations useful for studying electromagnetic fields in material media is  D (1.1a) 0 B (1.1b) t  E B 0 (1.1c) t  H D j (1.1d) while the Lorentz force on a point charge q moving at velocity v is  q  F E v B (1.2) Here, E is the electric field, D is the displacement field, B is the magnetic field, H is the magnetizing field,  is the free charge density, j is the free current density, and q is a charge. The electric and displacement fields are related by the polarization response of the medium 0 D E P , while the magnetic and magnetization fields are related by the magnetization response of the medium 0 1   H B M , where P is the polarization, M is the magnetization, and 0 and 0 are the permittivity and permeability of free space. In this work, we consider propagation in dielectric materials where magnetization is negligible M 0 . Field ionization will lead to the presence of free charge, however, once created, this negative charge density is locally balanced by the positive charge density from the ions so that 0  . This situation can change with 4 large electromagnetic fields (far greater than those found in air filaments considered in this dissertation), for which the ponderomotive force can induce electron drift with respect to the ions, resulting in charge separation. In what follows, we include the bound and free electron response in the material’s polarization P . 1.2.2 Linear polarization response For sufficiently small field amplitudes, bound electrons make small excursions from their equilibrium positions resulting in a polarization response of the medium. The motion of the massive nucleus ( 3electron proton 10 m m  ) is negligible by comparison. Much insight can be obtained using a classical model for the polarization response of a single bound electron, where the electron motion is modelled using a damped, driven spring (Lorentz model of the atom). The damping force dampF m x  models momentum changes (primarily through radiation and collisions), the atomic binding force on the electron is 2 atom 0F kx m x    where 0 k m   is the resonant frequency of the atomic oscillator of spring constant k, and the driving force is driveF eE  . Here we neglect magnetic contributions to the electron motion, as v E E c v B for nonrelativistic electrons. At optical frequencies, the magnetic force term (and relativistic electron dynamics) becomes significant at intensities 17 210 W cm  , where the normalized vector potential 20 0 / ~ 0.1a eA mc  . The classical model becomes 2 0F mx m x m x eE      (1.3) 5 Assuming a time harmonic field  0( ) Re i tE t E e  the dipole response of a single oscillator in steady state is   2 02 2 0 ( ) ( ) Re i t e p t ex t E e m i                  (1.4) When the driving frequency is far slower than the resonant frequency 0  , the polarization response is nearly in phase with the electric field. In air, the nearest resonance to the optical frequency band is 0 0 2 240nm c     [13], resulting in a characteristic response time of 0.8fs  . For the titanium sapphire laser used in this thesis, which generates pulses with central wavelength 800nm  , the response can be considered instantaneous. Generalizing to a medium of N molecules per volume, where each molecule has jf electrons with natural frequencies 0 j and damping coefficients j , the medium polarization is given by 2 2 2 0 j j j j fNe m i            P E (1.5) Equation (1.5) specifies a susceptibility in the frequency domain 0( ) ( ) ( )    P E . In general, the susceptibility is a 2 nd rank tensor, as the structure of a molecule allows for an orientation dependent polarization response. This effect is important in crystalline materials, where the molecular constituents are all oriented in the same manner, but in gases comprised of randomly oriented constituents, the polarization response is isotropic, and the susceptibility can be treated as a scalar. 6 2 2 2 0 0 ( ) j j j j fNe m i                 (1.6) Note that one option is to include the steady state free electron response in Equation (1.6) by using index 0j  to identify 0 ef N N as the free electron density, 00 0  as the free electron resonance frequency and 0  as the electron collision frequency, giving the free electron contribution to ( )  as   2 ( ) p e i          where 1 2 2 0 e p N e m          is the plasma frequency. The electric and displacement fields are related by 0( ) ( ) ( ) ( ) ( )        D E P E , where  0( ) 1 ( )      is the frequency dependent permittivity. This yields the frequency dependent refractive index     1 2 1 2 0( ) ( ) 1 ( )n c        , where   1 2 0 0c     is the speed of light in vacuum. 1.2.3 Wave propagation in linear media Deriving a wave equation from the time-domain Maxwell’s equations (1.1) yields,   2 2 02 2 1 c t t t                   E P E E M j (1.7a) Here, the free current density j and the polarization bound current density t   P are shown as separate terms, but can be combined as discussed in Section 1.2.1. 7 above. Note that in general the dielectric function ( , ) x may both depend on frequency and position (the medium may be nonuniform) so that    E 0 . However, for media with sufficiently mild nonuninformity such that E E     ,    E 0 is an excellent approximation. This gives 2 2 2 02 2 2 1 c t t         E P E (1.7b) where we have combined the free electron response with the bound response on the right side j P P t t       , and taken the magnetization to be zero M 0 . Equation (1.7b) can be Fourier transformed ( , ) ( , )t E x E k and ( , ) ( , )t P x P k , where we may use the frequency domain constitutive relation  0( , ) ( ) ( , )     P k E k ,   2 2 0 02 ( , ) ( , ) ( ) (k, )k c          E k + E k = E (1.8) Which yields the dispersion relation ( )k n c   , where n is the index of refraction as discussed above, and k  k . Because N , 1 is the case for all gases up to several atmospheres in pressure and we can approximate 1 1 2 k c c              , which gives the real and imaginary part of the wavenumber r ik k ik  ,  Re 1 2 rk c         and Im 2 ik c         . Thus, the real and imaginary parts of the susceptibility govern the phase evolution and attenuation of waves travelling through the medium. 8         2 22 0 2 2 2 2 0 0 Re j j j j j fNe m              (1.9a)       2 2 2 2 2 0 0 Im j j j j j fNe m             (1.9a) Figure 1.1. Susceptibility of a dielectric near a resonance. Figure 1.1, shown above, plots the susceptibility as a function of frequency for a dielectric near a resonance. In the plot, we use 0 5PHz  and 1PHz  . When the driving frequency is close to the resonance frequency 0  , bound electrons oscillate in resonance with large amplitude, and considerable energy is lost to damping, leading to a larger attenuation coefficient ik . There is also a sharp dip in the real part of the susceptibility, called anomalous dispersion, characterized by phase 9 velocity Re 1 2 c c n         increasing as a function of frequency. Away from resonances, attenuation is suppressed   2 2 2 0 ik      , where 0 is the nearest resonance, and we are in the normally dispersive regime where phase velocity is a decreasing function of frequency. A medium is said to be optically transparent when resonances lie outside of the optical frequency span  400 1000nm . The situation is typified by low losses (often approximated as lossless) and normal dispersion. A general wave solution can be constructed by summing the Fourier components of Equation (1.18) as follows:   ( )3 0( , ) ( , ( )) i k t t d k k e       k x E x E k (1.10) where ( )k  is a constraint imposed by the dispersion relation. Of particular interest are paraxial electric fields, where the angular wavenumber distribution is narrowly concentrated about a particular value ( ) 1 ( ) k k     . Paraxial fields, which are directional and characterized by low divergence, are often used in optical experiments where long paths containing common optical elements such as lenses and mirrors place constraints on the transverse size and divergence of optical fields traversing the beam line. Consider a paraxial electric field travelling along the z-direction, which, for simplicity, we take to be time harmonic and propagating in vacuum (dispersionless):   ( , ) ( ) c i k z t t e  E x A x . Here we have factored the field into a planar term 10 corresponding to the central wavenumber c ˆckk z , and a field envelope A , which satisfies  2 22 c z zik    A A (1.11) Here, 2 2 2 x y    is the transverse Laplacian. Because the paraxial field is comprised of plane waves with a narrow distribution centered about c ˆckk z , the change in E along the propagation direction is dominated by the planar term so that 22 c z zik  A A which is known as the paraxial approximation. More rigorously, note that    3 ( )c i k z t i t c ze ik i d k e k z z            k xE A E E k , where we express the derivative on E using the envelope function (first equality) and the Fourier decomposition (second equality). Writing the Fourier component as z ck k k  , and assuming E is paraxial, so that c 1 k k  we establish that        3 3( ) ( )ci k z t i t i tc ce i d k e k ik d k e ik z              k x k xA E k E k E , from which can we recover the paraxial approximation 22 c z zik  A A . For a beam comprised of linearly polarized plane waves, the spread in angle of the polarization unit vectors cen 1 k k     is similarly small, resulting in a beam with polarization that is nearly spatially invariant 0 ˆ ˆ( ) A x A . Thus, Equation (1.11) simplifies to the time harmonic, paraxial, scalar envelope equation 22 zik A A   (1.12) 11 Of special relevance is the 00TEM fundamental laser mode solution to Equation (1.12), where 22 2 ( ) 2 ( )( )0 00 0( , ) ( ) krr i z R zw zwE z E e e w z            r (1.13) Here,  ,x y r is the transverse vector, 0E is the field amplitude, 2 0( ) 1 R z w z w z         is the z-dependent spot size, 2 0 R w z    , the characteristic length scale for the field envelope called the Rayleigh range, 0 ( 0)w w z  is the beam waist,  is the wavelength, 2 R( ) 1 z R z z z           is the beam radius of curvature, and (z) arctan R z z         is the Guoy phase. The transverse amplitude distribution 2 2 ( ) r w ze  dilates with propagation, but is otherwise invariant, and exhibits the lowest possible divergence of any linearly propagating beam with equal or lesser beam waist. The fundamental mode is the lowest order solution in a variety of modal decompositions of paraxial waves, the most popular of which are the Hermite-Gauss or TEMmn modes which have rectilinear symmetry, and the Laguerre-Gauss or LGlp modes which have cylindrical symmetry. We make use of both modal decompositions in the work presented in this thesis as a means of generating filament arrays (see Chapter 3).  2 22 2 2 ( ) 2 ( ) ( )0 0 2 2 TEM ( , , ) H H ( ) ( ) ( ) k x y x y i z R z w z mn mn m n w x y E x y z E e e w z w z w z                               (1.14a) 12 22 2 2 ( ) 2 ( )( )0 2 2 2 LG ( , , ) L ( ) ( ) ( ) l l krr i z l R zw z lp lp p E r r E r z e e w z w z w z                            (1.14b) Here, Hm is the m th order Hermite polynomial and Llp are the generalized Laguerre polynomials, where m, n, and p, are all positive integers, and l is an integer. In both cases, the higher-order modes also have transverse amplitude distributions which are invariant under propagation up to dilation, the same phase curvature, and the same transverse expansion rate as the fundamental mode, but have different transverse distributions, and accumulate more Guoy phase through propagation   R ( ) 1 arctan z z N z          where N m n  for TEM modes and 2N l p  for LG modes. Next, we consider the propagation of polychromatic beams in dispersive media. Fourier transforming Equation (1.7b) with respect to time gives: 2 2 2 ( , ) ( ) ( , ) 0 c       E x E x (1.15) We expand the permittivity about a central frequency c  ,     2 c2 c 0 ( ) ( ) ( ) ! n n n n                , where  2 c( ) ( ) n       is the n th derivative of 2 ( )   evaluated at c  , and exchange powers of  c  for temporal derivatives. After inverse Fourier transforming, we have     c c2 2 c2 0 ( , t) ( ) ( ) ( , )e 0 ! i t n n i tn t n e i t c n              E x E x (1.16a)  2 2 c2 0 1 ( , t) ( ) ( ) ( , ) 0 ! n n n t n i t c n            E x E x (1.16b) 13 Where in Equation (1.16b) we have decomposed the field into a temporal envelope and a carrier frequency ( , ) ( , ) c i t t t e E x E x , a form which is useful for describing a pulsed field with a relatively narrow spectral distribution surrounding a central frequency c . Equation (1.16b) is a useful form of Equation (1.7) suitable for handling polychromatic fields. Equation (1.16) can be simplified considerably by truncating the series for ( )  at some finite order. Generally speaking, more dispersive media, broader spectral bandwidth, and longer propagation lengths will result in increased sensitivity to higher order corrections to ( )  . The first order correction is 2 02 1 2 cc c k k c               , which is proportional to c g 1k v      . Here, gv is the group velocity, which alters the pulse propagation speed. Higher moments distort the shape of the pulse envelope by accounting for a nonlinear spread in the phase velocity of neighboring frequencies, with the first correction being 2 2 2 02 2 2 1 2 c c c k k c                 , where c 2 2 k       is the group velocity dispersion. As a final step, we combine the various approximations made in this section: paraxial, scalar, and truncated dispersion to create a pulse envelope evolution equation useful for modelling linear ultrafast pulse propagation. For the purposes of numerical simulation, it is useful to apply a change of coordinates: gv t z   and z z  (where we neglect the prime going forward). The transformation specifies a longitudinal coordinate  which co-moves with the pulse at the group velocity (in the original coordinates, the pulse will move along the longitudinal direction to the edge 14 of the numerical window, but in the transformed/co-moving the pulse remains stationary along the longitudinal axis). Applying the approximations results in   2 222 z ik A A         (1.17) where A is an electromagnetic field envelope, such as E . Equation (1.17) is used throughout the thesis as the linear component of a numerically solved model for filament propagation (see Section 1.3.4 for the entire model). In many cases, temporal variations in the pulse envelope are slow relative to an optical cycle [14], and Equation (1.17) may be simplified by approximating ik ik  , which is known as the temporal slowly varying envelope approximation (the temporal analogy of the paraxial approximation). 2 2 22 zik A A       (1.18) Equation (1.18) is more symmetrical and easily manipulated for the purposes of making analytical arguments. 1.3 Nonlinear optics 1.3.1 The anharmonic oscillator The linear polarization model developed in Section 1.2.2 was predicated on the assumption of small field amplitudes driving small excursions of bound and free electrons. Considering the bound response only, larger electron excursion amplitudes can be modelled by small nonlinear corrections to the restorative force.  20 NLF m x x x eE F       (1.19) Provided the atom or molecule in question is rotationally symmetric, or we are considering an ensemble average response x , NLF over a gas of randomly 15 oriented molecules, the requirement that NL NL( ) ( )F x F x   implies that NLF contains only odd terms 3 5 NL 3 5 ...F k x k x    , with lower order terms typically dominating over higher order terms 2 2 i i i ik x k x   . Consider the first nonlinear term 3 NL 3F k x  . Equation (1.19) can be solved perturbatively using (1) (3) ...x x x   where the terms ( )ix are ordered in decreasing magnitude: ( ) ( 2)i ix x  .       st st nd nd 3 (1) 2 (1) (1) (3) 2 (3) (3) (1) 0 0 3 higher order1 order 1 order 2 order 2 order ...m x x x m x x x eE k x            (1.20) Equation (1.20) groups terms from Equation (1.19) into different orders. By inspection, we see that the first order term (1)x satisfies the linear model of Equation (1.3), while the third order term (3)x can be viewed as being driven nonlinearly by (1)x . Critically, even if the incident field is harmonic,  0 1 . . 2 i tE E e c c  , making the first order solution also harmonic  (1) (1)0 1 . . 2 i tx x e c c  , the third order term, driven by   3 (1)x is driven at frequencies  and 3 . Assuming a harmonic incident field  0 1 . . 2 i tE E e c c  , a polarization (1) (3)P P P Nex    can be derived in a manner similar to Section 1.2.2. For convenience, we define the polynomial 2 2 0( )D i      . 244 3 3 0 0(3) 33 0 4 3 4 3 31 1 ( ) . . 2 4 (3 ) ( ) 2 4 ( ) ( ) i t i t k Ne Ek Ne E P t e e c c m D D m D D                        (1.21) 16 Nonlinear corrections to the restorative force have led to a third order susceptibility (3) entailing corrections in the medium response at the fundamental frequency, as well as the generation of a new frequency at the third harmonic. 1.3.2 Formalism for nonlinear optics The example considered above, the anharmonic oscillator, can be generalized to explain a host of phenomena in nonlinear optics. In materials which violate inversion symmetry, even order terms allow for second harmonic generation, as well as effects sensitive to the orientation of the field polarization. By considering a driving field composed of three distinct frequencies: 1 , 2 , and 3 the third order medium response allows for general four-wave mixing processes          where  ,  , and  can each be any of 1 , 2 , and 3 , allowing for a total of 22 distinct output frequencies. Like the Lorentz model considered in Section 1.2.2, the anharmonic oscillator should be viewed as a phenomenological model for the bound electronic response. In practice, physicists do not seek to directly measure spring constants ik , nor is it always possible to perform ab initio quantum mechanical calculations on physically realistic models of molecules, instead, the medium is usually experimentally probed with incident electric fields to determine the polarization response. In general, the polarization response to a set of discrete optical frequencies is (1) (2) (3) (1) (2) (3) 0 0 0... ...              P P P P χ E χ E E χ E E E (1.22) 17 where ( )NP denotes the N th order nonlinear polarization, ( )Nχ , the Nth order nonlinear susceptibility, is a tensor of rank 1N  , and E is a monochromatic electric field. Elements of the N th order nonlinear polarization can be computed for a set of monochromatic frequencies as follows:         0 0 1 1 1 2 ( ) ( ) 0 1 2 ... 1 2 ... 1 2 ; , ,..., ; , ,..., ( ) ( )... ( ) N N N N N n N n n n N n n n N P K E E E                       (1.23) Here, the nonlinear polarization oscillates at 1 N a a     , and we sum over all possible sets of N frequencies  1,..., N  summing to  . Note that a real, monochromatic field   1 ( ) ( ) ( ) 2 i t i tE t E e E e     provides equal spectral power at both positive and negative frequencies, as *( ) ( )E E   . K is a degeneracy factor which accounts for different permutations of a fixed set of input frequencies (see ref. [15]), and there is an implicit Einstein summation over the electric field polarization directions  , ,in x y z for 1,2,...,i N . The time domain polarization response corresponding to the term in equation (1.23) is    0 0 ( ) ( ) , 1 ( ) . . 2 i tN N n nP t P e c c        . As an example, for the anharmonic oscillator (considered above), the third order polarization contribution for the fundamental frequency can be expressed as (3) 0 3 ( ) ( ; , , ) ( ) ( ) ( ) 4 x xxxx x x xP E E E            , where 3 4 K  , we have taken the 18 scalar electric field to be oriented along the x direction, and 4 3 4 3 0 ( ; , , ) ( ) ( ) xxxx k Ne m D D              . Equation (1.2.3) can also be used to model pulsed electric fields in the adiabatic limit. Consider a pulsed field   1 ( ) ( ) . . 2 i tE t E t e c c   , where ( )E t is a temporal envelope with bandwidth  . Provided the susceptibility tensor ( )N is nearly constant (dispersionless) over the frequency range 2     , we are in the adiabatic limit and Equation (1.2.3) may be used for the pulsed field as well. More generally, Equation (1.2.3) can be integrated over the spectral domain to determine the polarization response from pulsed fields:     0 0 1 1 1 2 ( ) ( ) 0 1 ... 1 1 2 1 1 ( ) ... ; ,..., ; ,..., 2 ( ) ( )... ( ) . . N N a a N N N n N n n n N Ni t n n n N a a P t K E E E e d c c                          (1.24) Using the Convolution theorem, the right side of Equation (1.2.4) can also be expressed in the time domain: 0 0 1 1 ( ) ( ) 0 ... 1 1 1 ( ) ... ( ,..., ) ( )... ( ) N N N N N n n n n N n n N a a P t t t E E d             (1.25) For a general review of nonlinear optical phenomenon see ref. [15] In the remaining sections we focus on nonlinear processes relevant to filamentation, especially in air. 1.3.3 Optical Kerr effect 19 The third order susceptibility tensor (3) ( ; , , )xxxx      , responsible for the optical Kerr effect, is of special relevance for filamentation, as it causes the field to modulate its own phase. Assuming a linearly polarized electric field, an effective susceptibility can be defined by a combination of first and third order susceptibilities 2(1) (3) eff 3 ( ; ) ( ; , , ) ( ) 4 xx xxxx E              which can be used to construct an intensity dependent refractive index eff eff 0 21 ( )n n n I     , where 2 (1) 0 1 ( ; )xxn      , 2 0 0 1 ( ) ( ) 2 I n c E   is the spectral intensity, and (3) 2 2 0 0 3 ( ; , , ) 4 xxxxn n c          is the Kerr coefficient. The optical Kerr effect produces frequency shifting along the temporal axis referred to as self-phase modulation. Assuming an instantaneous response, the frequency shift after propagation over a distance L is  2( ) ( ) d d t kn I t L dt dt       . For a pulse with a bell-shaped temporal intensity profile, red-shifting occurs at the temporal front of the pulse, and blue-shifting occurs in the rear. The spatial analog to self-phase modulation is self-focusing, and occurs when a bell shaped transverse profile establishes its own intensity dependent lens. As discussed in Section 1.2.2 on the linear polarization response, the bound electronic response time 0.8fs  can be viewed as effectively instantaneous for the 800nm  pulses considered in this thesis, which have period 800nm 2.7fs  and bandwidths satisfying c 1    . With respect to the electronic Kerr effect, the 20 response is in the adiabatic regime, and we may thus approximate the time domain response as instantaneous 0 2,elec( ) ( )n t n n I t  . An additional intensity dependent effect exists for diatomic molecules, such as the air constituents 2N and 2O . The long axis (symmetry axis) of a diatomic molecule exhibits greater bound electron mobility and thus a stronger polarization response relative to the directions orthogonal to the symmetry axis. If a diatomic molecule is not aligned with the incident field polarization, the anisotropic polarization response exhibits a torque on the molecule, where the torque depends on the intensity of the field envelope of the pulse. In the presence of an intense electric field, the electron cloud distorts and thus “drags” the nuclei to align with the electric field polarization, producing a transient birefringence of the medium. The rotational response time of the molecule is determined by the moment of inertia of the molecule, and is naturally much slower than the bound electron response because of the large difference in mass. In air, the response time is 100fs [16], meaning that a sufficiently short pulse will induce a transient birefringence in its wake without experiencing the effects itself. Classical and quantum mechanical derivations for the rotational response of the molecules are both possible. The two approaches have good agreement provided the initial thermal distribution populates a sufficient number of rotational states, a condition which is satisfied when the temperature and moment of inertia of the molecules are sufficiently large [17]. One important fundamental difference is that the quantum mechanical model exhibits periodic realignments, or “quantum echoes”, which have been measured in 21 detail for air constituents [18] and can have significant effect on filamentation [19]. Echoes are seen because the rotational states are all periodic over the fundamental periodicity set by the lowest energy rotational state. In the classical model, where molecular rotation occurs for a continuous distribution of speeds, no such realignment is possible. For the symmetric molecules N2 and O2, full realignments occur at the fundamental periodicity, with partial realignments also occurring at 1 4 , 1 2 and 3 4 of the fundamental periodicity. Air exhibits its largest realignment at 8.4ps , where the full realignment of nitrogen roughly coincides with the 3 4 revival of oxygen. Realignments can be used to guide subsequent pulses travelling in the wake of the first pulse [19,20]. It is possible to model the rotational response as an effective (3) susceptibility, where the slow response of the molecules necessitates a temporally nonlocal response as in Equation (1.2.5). The time domain index response is given by: 2,rot( ) ( ) ( ) t n t n t I d       [21]. 1.3.4 Ionization The ionization potential of air constituents 2N and 2O are 2N 15.2eVU  and 2O 12.1eVU  , while single photon energies in the optical spectrum are 1-2eV , meaning that multiple photons must be absorbed simultaneously to liberate a bound electron. For low field intensities, the multiphoton ionization rate can be modelled as K KW I , where W is the ionization rate, K is the K-photon ionization cross section, I is the field intensity, and K is the minimum number of photons required to 22 overcome the ionization potential K U  . The rate can be derived from time- dependent perturbation theory by considering the lowest order relevant term. As the field intensity or optical period increase, photoionization enters the tunneling regime, where the electron tunneling time is now less than an optical period, and the lowest order perturbation term is no longer sufficient for explaining ionization. In the tunneling regime, one considers the combined potential of the coulomb field and incident field 0 ( ) 4 qZ q t r  E r , where 1Z  is the charge of the atomic core. The potential dips below the energy level of the bound electron at a finite distance, allowing the electron to tunnel out of its bound state. Figure 1.2, shown below, diagrammatically displays the difference between multiphoton ionization and tunneling ionization. Multiphoton ionization, displayed on the left of Figure 1.2, shows an unperturbed coulomb potential in black, and the simultaneous absorption of several photons (red arrows) to raise the electron energy above the barrier. Tunneling ionization, on the right of Figure 1.2, shows the combined incident field and coulomb potential distorting to levels where there is significant probability for the electron to tunnel beyond the resulting potential barrier. 23 Figure 1.2. Schematic depiction of multiphoton ionization (left) and tunneling ionization (right). A measure of the degree to which ionization is in the multiphoton or tunneling regime is given by the Keldysh parameter i p2 K U U   , where iU is the ionization potentiation and 22 p 2 e 04 e E U m   is the ponderomotive potential of the incident field. For 1K , ionization is dominated by multiphoton absorption, while for 1K , tunneling is the dominant mechanism. In air, where the lowest ionization potential comes from oxygen 2O 12.1eVU  , and filaments with pulses at 800nm  reach peak intensities of 13 25 10 W cmI   (see Section 1.4.3), 1.6K  , indicating that multiphoton absorption is favored, but we are in an intermediate regime. Keldysh was able to unify the two regimes into a single model [22], which was further refined by Perelmov, Popov and Terent’ev into what is now referred to as the PPT model [23]. As already discussed in Section 1.2.2, once electrons are born, their susceptibility can be modelled using the Lorentz model with the atomic restoring e- e- ħω 24 force set to zero, and the damping force replaced by a collision rate) e e m    E x x , which gives a susceptibility   2 e 2 e 0 1n e m i        (1.26) Where en is the free electron density, and  is a collision rate. For the low density plasmas generated by femtosecond filaments in air (ionization fraction ~0.001 [24]), the electron-neutral collision rate en 2THz  so that en 1   [25]. The free electron susceptibility can then be simplified to e c n n    , where 2 0 c 2 emn e    is the critical plasma density. Avalanche ionization due to electron-neutral collision is similarly negligible for the 100fs pulses considered in this dissertation, since 1avalanche [26]. Plasmas generated by filamentation in air are quasi-neutral: electron-ion charge separation can be shown to be negligible by comparing the ponderomotive force 2 2 pond 2 e4 e E m   F to the electrostatic force field established due to space charge separation between the electrons and ions. The electron and ion charge distributions can be modelled as cylinders of constant charge, which have corresponding electric fields inside 02    r E inside the cylinder, and 2 outside 0 ˆ 2 R r   E r outside the cylinder, where  is the charge density, r is the distance from the axis of the cylinder, and R is the radius of the cylinder. Displacing the electron cylinder by a small distance x R 25 results in an electric field magnitude 02 x E     in the sliver where the cylinders are non-overlapping, from which the space charge felt by a single particle may be estimated s.c. 02 e x F     . Equating the two forces gives the ponderomotive force- induced displacement as 2 0 2 e2 e E x m       . Estimating 2 2 core E E r  , where the filament core radius is core 50μmr (see Section 1.4.3), and the peak intensity achieved in the filamenting core as 2 13 2 0 1 10 W cm 2 c E  (also see Section 1.4.3),and 16 310 cme   [24], gives a displacement core0.2nm rx  . As a result, to a very good approximation, the plasma can be treated as quasi-neutral. 1.3.5 Nonlinear polarization In Section 1.2.2 we modelled the linear polarization response derived from the Lorentz model of the atom, neglecting all nonlinearities. Now we incorporate nonlinearities stemming from the bound electronic response, molecular rotational response, and the free electronic response of the medium. L NL 2 (1) (3) 0 pl ( , ) ( , ) 3 ( , ) ( , ) ( , ) ( , ) ( , ) ( , ) 4                         P r P r P r r r E r r E r (1.27) Here, L ( , )P r and NL ( , )P r are, respectively, the linear and nonlinear polarization. The linear susceptibility, as discussed in Section 1.2.2, is (1) ( , ) r , (3) incorporates bound electron and rotational susceptibilities, and pl models the free electron 26 response. The nonlinear polarization can be used to extend the linear pulse propagation Equation (1.17) to a nonlinear propagation equation suitable for modelling filamentation. To derive the result, start at the frequency domain wave Equation (1.15), splitting the polarization into linear and nonlinear components as in Equation (1.27), and then follow the derivation of Section 1.2.3 (inverse transform back to the temporal domain, transform to the moving frame gv t z   , apply paraxial approximation, and truncate linear dispersion at a finite order), this time keeping the nonlinear polarization on the RHS as a driving term.     2 2 2 2 NL 0 1 ˆ2 z ik A ik P               (1.28) Here, A and NLPˆ are envelope functions for the scalar electric field ( , , ) ( , , ) ikE z A z e    r r and the time domain nonlinear polarization NL NL ˆ( , , ) ( , , ) ikP z P z e    r r , which can be separated into NL,elec 0 2,elec ˆ ( , , ) 2 ( , , ) ( , , )P z n I z A z     r r r , NL,rot 0 2,rot ˆ ( , , ) 2 ( ) ( ) ( , , )P z n I d A z              r r , and NL,plas 0 ( , , )ˆ ( , , ) ( , , )e c n z P z A z n       r r r , with rate equations governing the populations of the free electrons and neutrals 2 2 2 2 e N N O O dn n n d      , 2 2 2 O O O dn n d     , and 2 2 2 N N N dn n d     , where i are ionization rates given by the PPT model, and in are number densities for the free electrons and neutrals. Note that we have used the 27 adiabatic approximation (see Section 1.3.2) in expressing the time-domain polarization response from plasma as well as the electronic Kerr effect 1.4 Filamentation of high power optical pulses 1.4.1 Filamentation summary When a high intensity pulse propagates through a dielectric medium it generates its own co-propagating lens due to the optical Kerr effect. If the pulse power is above a critical power threshold crP , the self-lens is able to overcome diffraction and the pulse begins to shrink in size. As the pulse shrinks in size and the peak intensity grows, the self-lens becomes increasingly dominant over diffraction in a self-reinforcing process known as optical collapse. The collapse is invariably mitigated by a “collapse arrest mechanism” which only becomes relevant to propagation as the pulse becomes increasingly singular. Following collapse arrest the pulse begins to “filament”, a dynamical competition ensues between focusing and defocusing nonlinearities that support extended propagation of a high intensity filamenting core over many Rayleigh lengths. Figure 1.3, below, shows a schematic representation of filament self-channeling in air, where self-focusing from the Kerr effect competes with refraction from plasma generated from field ionization. An end mode image of an air filament is also displayed, in which the intense and spectrally broadened “core” contrasts with the lower intensity “reservoir” (see Section 1.4.3). 28 Figure 1.3. Left: schematic representation of dynamic self-channeling of air filament due to competition between self-focusing nonlinearity and plasma refraction. Right: end mode image of air filament. 1.4.2 Optical collapse and collapse arrest Optical collapse, wave propagation where self-focusing dominates diffraction, is the first phase of filamentation. Self-focusing due to the optical Kerr effect is ubiquitous, as the process is automatically phase matched, produced from a single frequency, present regardless of the symmetry of the material, and often the most pronounced effect at lower intensities. Optical collapse was first experimentally observed by Hercher in 1964 [27], in a series of experiments where high power lasers were seen to drill transparent glass plates, leaving cylindrical damage patterns up to 2cm long and only several wavelengths wide. The large aspect ratio of the damage patterns cast doubt on potential explanations using linear optics. Chiao [28] and Kelley [29] were the first to identify self-focusing nonlinearities as responsible for lensing the beam. Following Kelley’s original arguments to demonstrate optical collapse and estimate the critical power, we consider a monochromatic beam of sufficient intensity as to be lensed by the optical Kerr effect. plasma self-focusing defocusing 29 2 2 2 0 2 2 ( ,z) ( , z) ( , z)z k n ik A I A n               r r r (1.29) The monochromatic propagation equation (1.29) can be derived from the filament model considered in Equation (1.28) by taking the long pulse limit, where ik  , dispersion is negligible 2 2 0   , the medium response is in the adiabatic limit, and we focus on the optical collapse phase of filamentation where the Kerr effect is the only relevant nonlinearity. Equation (1.29) is known as the Nonlinear Schrodinger Equation, or NLS. Suppose the beam is of sufficient intensity, that initial propagation is dominated by the nonlinear term 2 2 22A k n IA , so that the phase curvature of the beam evolves due to nonlinear lensing much more quickly than beam changes due to diffraction. The evolution of the beam at this early stage is then approximated by 2 ( ,z 0)( ,z) ( ,z 0)eikn I zA A     r r r (1.30) where ( , 0)A z r is the beam envelope at that the start of propagation. We may now estimate a length scale over which diffraction becomes relevant by estimating the size of the diffractive term using the solution in Equation (1.30). For definiteness, consider the initial condition 2 2 02 0( , z 0) r A E e     r r . By imposing 2 ( ,z 0)2 22 ( ,z) ( ,z) ( ,z 0) e ikn I z zik A A A              r r r r , we arrive at a self-focusing scale length 0 0SF 22 r n L n I  . By considering a diffractive scale length in the absence of nonlinearities 2 0 0 diff 2 1.22 r n L   , we obtain Kelley’s original estimate for the onset of 30 optical collapse: SF diffL L which can be rearranged to create a condition on the pulse power   2 2 cr 0 0 2 1.22 16 P P r I n n     . If self-focusing can overcome diffraction, can self-focusing be tuned to balance it? Solitonic solutions to the NLS (Equation (1.29)), can be found by assuming a propagation invariant solution where diffraction exactly balances self- focusing 2 2 2 0 2 ( ) ( ) ( ) k n A I A n      r r r , from which an entire family of solitonic solutions may be derived, with the lowest order solution known as the Townes profile. All solutions to the stationary NLS are unstable; infinitely small perturbations to the solitons either lead to catastrophic collapse of the wave, or eventual diffractive spreading, making physical realizations of solitonic propagation in this regime impossible [30]. Each transverse beam distribution has its own critical power, with the critical power for the Townes profile providing an absolute lower bound on the power necessary to initiate optical collapse 2 cr,Townes 0 2 3.72 8 P n n    . The Gaussian initial distribution relevant in many laser experiments results in a slightly larger critical power 2 cr,Gauss 0 2 3.77 8 P n n    . According to the NLS, Equation (1.29), a beam propagating above its critical power leads to catastrophic collapse, where, in a finite distance, a finite amount of energy in the beam aggregates to a point. The collapse length cL is generically a function of the linear diffractive length diffL and the ratio of the power to the critical 31 power cr P P . For the case of a Gaussian beam, Dawes and Marburger determined a formula for approximating the collapse length [31] by numerically solving the NLS for a variety of initial conditions, and fitting the obtained collapse length to an analytical formula diff c 2 1/2 0.367 0.852 0.0219 cr L L P P             (1.31) Of course, catastrophic collapse is not a real physical scenario. Invariably, an additional physical mechanism termed the collapse arrest mechanism becomes relevant as the pulse becomes increasingly intense. In air and other gaseous media, the collapse arrest mechanism is refraction from plasma which is generated by field ionization from the intense beam [1]. In water and many other solids, spectral broadening from increased self-phase modulation during collapse results in a dispersive arrest where temporal stretching of the pulse results in a lowered intensity [32]. Even in a plasma, which exhibits a self-focusing nonlinearity resulting from relativistic corrections to the quiver trajectory of free electrons, electron cavitation from the ponderomotive force of the beam leads to collapse arrest [33]. In the absence of an additional material effect, the NLS (Equation (1.29)) eventually becomes non-physical as the beam becomes non-paraxial, at which point it can be shown that vectorial effects will also arrest collapse [34]. 1.4.3 Core/reservoir model of filamentation 32 In the optical collapse scenario detailed in Section 1.4.2 an initially collimated beam with transverse scale length 0r begins to self-focus, transferring energy inward until the on-axis intensity diverges and a finite amount of energy has aggregated to a point [30]. Even in this unphysical collapse scenario absent an arrest mechanism, a substantial portion of the pulse energy remains distributed over the original transverse scale length 0r . Collapse arrested filamenting beams possess two transverse scale lengths, one for the small, high intensity filamenting “core”, and another for the larger, low intensity filament “reservoir” which surrounds the core and contains the majority of the energy in the pulse [35]. Typically, the arrest mechanism is active in the core, with Kerr lensing guiding additional energy inward from the reservoir, fueling continued propagation of the core. In an air filament, the core may be loosely defined by the maximal radius where the pulse is intense enough to produce ionization. The relatively small filamenting core produces an even narrower plasma column (due to the highly nonlinear dependence of ionization on intensity) which then acts to refract light out of the core. Courvoisier et al., conducted experiments where water droplets impeded the propagation of the filamenting core [36]. Beyond the droplet, the reservoir was seen to form a new core, indicating the “self-healing” property of filamenting beams. Provided the filamenting reservoir contains more than the critical power for self- focusing, optical collapse should continue, forming a new core. The intensity in the core of an air filament can be estimated by balancing the refractive index contributions of Kerr self-focusing and plasma refraction 33 (defocusing) e 2 c ( ) 2 n I n I n  , where the multiphoton ionization rate may be used to estimate the electron density (see Section 1.3.3) e ( ) ~ K K pn I I  . Here, K is the K- photon absorption cross section and p is the pulse length. Slight increases in the intensity lead to much higher electron density and stronger defocusing from plasma refraction. As a result, the intensity is effectively clamped somewhere near the balance set by Kerr self-focusing and plasma refraction 13 2 clamp ~ 2 10 W cmI   . This estimate is borne out by measurements of air filament intensity by Lange et al. [37], who found that an air filament that abruptly transitions to an argon gas cell has a harmonic cutoff of 23, which can be linked via simple theory [38] to a peak intensity of the laser of 13 25 10 W cm  . Additionally, assuming that each filamenting core contains roughly a critical power, the transverse scale of the core can be estimated as cr core clamp ~ 40μm P r I  . 1.4.4 Spectral broadening and self-steepening Phase shifts due to the optical Kerr effect and the plasma response can lead to spectral broadening. Ignoring diffraction, total phase accumulation of the pulse may be modelled as ( , ) ( , ) n t z t kz t t c             x x , with the instantaneous frequency given by einst 2 c 1 2 nd z I n dt c t n t                where we use an instantaneous Kerr response for simplicity. The Kerr response causes red shifting of the leading (rising) 34 edging of the pulse, and blue shifting of the trailing (falling) edge. Ionization, which is strongly nonlinear in the intensity, is negligible up to an intensity threshold, at which point the free electron density rapidly increases and saturates as the intensity clamps (see Section 1.4.3 for a discussion of intensity clamping). The sharp rise in free electron density contributes to a blue shifting of the pulse. Temporal variations in the index are also responsible for a reshaping of the amplitude of the temporal envelope of the pulse in a process known as self- steepening. The index enhancement from the Kerr effect results in an intensity dependent group velocity gv k    , where the high intensity peak of the pulse travels slower than the low intensity wings. Self-steepening results in a stretching of the rising edge of the pulse, and a compression of the falling edge. The combined effects of spectral broadening and self-steepening results in self-compression of filamenting beams. Compression of a 45 fs, 800 nm, 5mJ pulse down to an 8 fs, 3.8 mJ pulse has been demonstrated by Stibenz et al. [3], and is competitive with the pulse compression capability of modern hollow core fibers [39] without the additional cost and complexity of dealing with several meters long rigid glass fibers. 1.4.5 Modulational instability and multiple filamentation For powers well above the critical power crP P , a filamenting beam will produce multiple filamenting cores. The cores are nucleated as the beam undergoes collapse, as well as throughout propagation as individual cores expire and new cores are created. During collapse the Kerr effect creates an instability in the transverse 35 intensity distribution, producing multiple high intensity peaks or “hot spots”. Campillo et al. [40,41], considered a planar monochromatic wave, and found that the growth rate for transverse spatial modulations was 2 cr 8 ( ) 2 k I k k k P      which has maximal growth for opt cr 2 I k P    . This leads to a characteristic distance between filaments fil opt 2 d k    , from which Couairon and Berge [42] estimated a typical power per filament of 2 cr 4 P  , consistent with an experimental study by Ting et al. [43], who found cr1P per filament for a 20 mJ, 50 fs beam filamenting at 800 nm in air. 1.4.6 Long timescale gas response to filamentation All of the discussion regarding filamentation has thus far focused on ultrafast phenomena, that is, phenomena with timescales similar to or faster than the pulse itself. We now consider the behavior of gaseous media following propagation of a filamenting pulse. A long, thin plasma is generated in the wake of a filamenting pulse [35]. Over a timescale depending on the gas type and pressure (~1 ns for air at 1 atm [24]), the plasma then recombines and repartitions most of its initial thermal energy into translational and rotational degrees of freedom of the neutral gas. If the gas is diatomic, the ultrashort pulse torques the ensemble of molecules, depositing energy by coherently exciting rotational states of the gas. The rotational excitation decoheres 36 through collisions (~100 ps collisional timescale for air at 1 atm [44]) where the energy is also repartitioned to translational and rotational modes of the neutral gas. In a filament, energy is primarily deposited by the high intensity filamenting core (~50 µm radius for an air filament), and, owing to the air low thermal conductivity, the deposited energy does not diffuse appreciably before it is converted to heated air in approximately the same volume. Seeking pressure balance, the heated channel of neutral gas launches a single cycle acoustic wave and expels mass out of the volume. Figure 1.4, taken from Cheng et al. [9], is a hydrodynamic simulation of the early time evolution of the gas as it relaxes into pressure equilibrium. Figure 1.4. Simulation of the time evolution of a gas density hole created in the wake of a filamenting pulse taken from Ref. [9]. At early times < 0.3 µs, the high pressure volume originally occupied by the filamenting core expels mass and launches a single cycle acoustic wave, reaching pressure equilibrium by ~1 µs. 37 After ~1 µs the gas reaches pressure equilibrium, and a hot, low density channel occupies roughly the same volume as the filamenting core [9]. The gas begins to cool through diffusion, and the density depression expands and fills in over a timescale of milliseconds. To verify diffusively driven evolution, Figure 1.5, also taken from Cheng at. al, shows the relaxation of the density hole created in a variety of gases. The log-log plot of time vs. hole FWHM displays a slope of 1 2 , characteristic of diffusion into the 2 transverse dimensions. Further, the expansion rate of the density holes was used to infer a thermal diffusivity for the medium, which was found to be within ~10% of values found in the literature [9]. Figure 1.5. Log-log plot of density hole FWHM vs time. Measured FWHM vs time is shown in circles, with slope of ½ fitted to the experimentally derived values. Data is taken from earlier work published in Ref. [9]. 100 200 300 400 500 600 700 delay (s) h o le s iz e ( F W H M ) ( m ) 60 100 300 500 1500 N 2 1 atm O 2 1 atm air 1 atm Ar 1 atm N 2 0.5 atm O 2 0.5 atm Kr 100 torr 38 Chapter 2: Optical beam dynamics in a gas repetitively heated by femtosecond filaments 2.1 Introduction In this chapter we consider the effects of long timescale neutral gas evolution on the beam pointing dynamics of high repetition rate (≥1 kHz) filaments. Most effects pertinent to the filamentation process itself occur on a timescale on the same order or shorter than the pulse length (typically ~100 fs). For air filamentation [35], this includes self-phase modulation induced electronically and rotationally through the optical Kerr effect, as well as ionization from multiphoton absorption and tunneling. Filamentation has undergone intensive investigation for applications stemming from these ultrafast effects, such as harmonic generation [45,46], supercontinuum generation [2], and generation of few-cycle pulses through self- compression [3,47]. Missing from prior analyses of filamentation are longer timescale effects such as thermal blooming, which has been found to distort the propagation of high power CW lasers in the atmosphere [48,49]. Despite the high peak power of pulses in femtosecond filaments, the average power in the beam is relatively small, so thermal effects in filamentation are subtle and were only reported recently [9]. In general, models of filamentation had assumed that each laser pulse interacts with a uniform medium and that any perturbations caused by the previous pulse have vanished by the time the next pulse arrives. Recently, it was shown that filamentation in gases at kilohertz repetition rates can be affected by the density depression that accumulates due to laser heating of the gas [9]. Lensing by the density depression 39 alters the propagation of the laser pulse, affecting the onset of filamentation and the spectrum of generated supercontinuum [9]. In this chapter, we show that this long-lived density depression, coupled with convective motion of the heated gas, leads to a reproducible deflection of the filamenting beam. We present time-resolved measurements of the accumulating kHz pulse-train-driven gas density depression inside a gas cell and the associated beam deflection. We find that the deflection is well-described by a simple model of beam refraction associated with the evolution of the gas density hole, whose time and space scales are well described by a simple fluid analysis. The results presented here provide quantitative understanding of thermal effects on beam propagation of femtosecond filaments, which will point the way to both stabilization and long range control of the filamentation process. Results are shown for Xe and air, but are broadly applicable to all configurations employing high repetition pulse rate femtosecond laser propagation in gases. For a review of the hydrodynamic response following filamentation, see Section 1.3.6 or the paper by Cheng et al. [9]. 2.2 Experimental setup Here, we consider a situation in which a shutter is opened suddenly and a 1 kHz pulse train is focused into a gas cell. The first pulse experiences a uniform gas density. The time separation between successive pulses is short enough that each pulse, which itself heats the gas, is affected by the cumulative density hole, which acts as a negative lens. As the hole evolves, the pulses in the sequence are steered, eventually reaching a steady state deflection. Here we measure and quantify the hole evolution and beam deflection. 40 Figure 2.1. Experimental setup: the pump and probe travel collinearly through the gas cell. The probe beam is relay imaged at the axial center of the plasma where the density hole is produced. The experimental setup is shown in Figure 2.1. A 1 kHz pump pulse train, apertured to 4mm with 45 fs, 120 µJ pulses was focused by a 60 cm lens into a gas cell filled with air or xenon at 1 atm. An additional setup was employed where a 40 cm lens focused the beam into a 2.7 atm gas cell also filled with xenon. The latter is the same setup and gas cell we use for supercontinuum generation for our spectral interferometry experiments [50]. A fast shutter synchronized to the laser system was employed to create a finite pulse train. To capture the transverse location of the pump beam in the far field we routed the output beam, attenuated using neutral density filters, from the gas cell onto a quadrant detector. The beam spot size on the 7.8 mm diameter quadrant detector was ~4 mm. The recorded beam deflections were downward. Gas density evolution was measured using a folded wavefront interferometer with a CW HeNe laser probe beam [9]. The probe propagated collinearly with the pump pulse train, and was relay imaged from the axial center of the plasma, through a folded wavefront interferometer and onto a CCD. The phase shift from the density gas cell HeNe laser (probe) CCD objective Ti:Sapphire regenerative amplifier (pump) relay imaging lenses folded wavefront interferometer quad detector fast shutter 41 perturbation was then extracted from the resulting interferogram using standard interferometric techniques [51]. The phase shift was converted to gas density using the known linear refractive index of Xe [52]. In the case of the 2.7 atm cell, the probe beam propagated at a 5 degree angle with respect to the pump, and Abel inversion was used to extract the size and shape of the density perturbation. Temporal resolution was achieved by time-gating the CCD and triggering the electronic exposure in order to sample the perturbation after a given pulse in the pulse train. The temporal resolution was limited by the minimum CCD exposure duration of ~40 μs, which was short compared to the density profile evolution diffusive timescale as previously measured [9]. 2.3 Results and discussion 2.3.1 Beam deflection A crucial aspect of the gas dynamics not investigated previously is that buoyant forces move the density depression upward. This is seen in Figure 2.2, which shows the 2D gas density profile measured in 1 atm xenon and air in response to filament generation at 1 kHz. Here, the filaments are 2.5 cm and 4.75 cm long for xenon and air from visual inspection of the plasma fluorescence. The profile shows the gas density 100 µs before the next pulse in the train. The accumulated density hole from earlier pulses, which has expanded to hundreds of microns in width, is seen to be displaced upward from the pump beam center, centered at (0,0) on the image, while the pulse train deflects downward. 42 Figure 2.2. Two-dimensional plots of gas density in 1 atm xenon and air created by 1 kHz pulse trains with 0.12 mJ and 1.5 mJ per pulse respectively. The image was taken in between pulses in the pulse train. The deep center region was created largely by the pulse prior to the probe, while the larger, shallower background is due to the accumulation of heat from earlier pulses in the train. Results from the interferometry experiment particularly useful for understanding beam deflection are shown in Figure 2.3. Here we use data taken from the 2.7 atm xenon gas cell setup, where the filament is ~3 cm long. Plotted are the hole vertical displacement relative to the pulse train, the hole depth, and the hole half- width-at-half-maximum (HWHM). Each point is an average of 5 measurements. Shot to shot variations were small. Y ( m m ) X (mm) xenon / -0.6 0 0.6 2 1.5 1 0.5 0 -0.5 -0.02 -0.01 0 air -0.6 0 0.6 -0.06 -0.04 -0.02 0 43 Figure 2.3. Top panel: Image of deflected white light supercontinuum spot as a function of pulse number in the 1 kHz pulse train. Lower panels: Hole parameters as a function of pulse number in the pulse train. The fast shutter is opened just before pulse 1. Equilibrium is approached by about the 250 th pulse (at 250 ms).The inset in the second panel is a typical relative density profile extracted from the interferogram. The hole displacement was measured by sampling the cumulative density hole just before and after a pulse arrived. The location of the newest pulse is marked by a small, sharp dip in the density superimposed on the larger cumulative density hole, as seen in Figure 2.2. The hole depth and HWHM were extracted from an Abel inverted lineout of the interferogram and were measured 100 µs before a pulse arrives. As can be seen from Figure 2.2, the hole is significantly more extended above than below its deepest point. Since the hole drifts upwards and the beam deflects downwards, it is the bottom portion of the hole that is relevant to the deflection dynamics, and it is the d e fl e c ti o n ( m ra d ) #1 #50 #100 #200 intensity (a.u.) -1.75 0 1.75 3.5 0 1 0 100 200 d is p la c e m e n t (  m ) 0 500 1000 0 0.01 0.02   /  distance (m) 0.05 0.1 0.15 d e p th (   /  ) 0 200 400 600 800 100 150 200 pulse number or time in milliseconds H W H M (  m ) 44 HWHM of the bottom side of the hole shown plotted in Figure 2.3. Also displayed in the top panel of Figure 2.3 is a sequence of filament white light beam images in the far field as a function of pulse number in the sequence, showing progressive deflection with negligible beam distortion. Furthermore, we use these beams in our spectral interferometry experiments [50] and detect no phase front distortion from the density hole-induced steering. We model the beam deflection by the density hole using the ray propagation equation, d d n n ds ds        r , applied in the vertical plane, where s is the optical path along the 2D ray trajectory, ˆ ˆz x r z x points to the ray tip, with zˆ along and xˆ vertically perpendicular to the propagation direction, and where 0n n n  is the refractive index. For our situation of a small refractive index perturbation 0 1 n n  , 0 1n and small ray deflection angles ( ) 1 dx z dz at large z , this equation simplifies to 2 2 d x d n dz dx   (2.1) The heated density depression relaxes through thermal diffusion [9], and so the index perturbation is approximately Gaussian: 2 2 0 exp( )n n x z      (2.2) Here  and  characterize the length scales of the perturbation in the x and z directions respectively, and 0n (<0 for a density hole) gives the maximum amplitude of the perturbation. In the experiment, the axial extent of the density hole  1 2  of 45 ~ 3 cm was measured by translating the object plane of the imaging system along the axis of the hole, and 0n and  were obtained by fitting the measured density hole to a Gaussian. Light ray trajectories are linear outside of the density hole region, and so the problem can be treated as classical scattering, where the amplitude and dimensions of the index perturbation, impact parameter, and deflection (scattering) angle are the relevant variables. The deflection angle  is given by integrating equation (2.1), using equation (2.2) as the index profile:    2 00 04 x z z d n x z ,z xdx dx dz n e dz dz dx             Δ Δ (2.3) The incident ray is assumed to be initially travelling parallel to the optical axis with impact parameter 0x , where in the experiment 0x corresponds to the density hole displacement. Figure 2.4. Far field beam deflection angle versus pulse number in 1 kHz plus train. Blue squares: measurement; red squares: ray optics calculation using measured gas density profiles. 0 200 400 600 800 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 pulse number or time in milliseconds d e fl e c ti o n ( m R a d ) y = A*[exp(-t/) -1]  = 108 ms quad model 46 Figure 2.4 shows good agreement between the measured deflection and deflection from the ray optics model. The deviation at earliest times is due to uncertainty in the density hole position for small displacements. The ray optics model contains no free parameters and uses the density hole parameters from Figure 2.3 measured simultaneously with the deflection. That this model works so well for filaments is consistent with density hole lensing of the filament ‘reservoir,’ the lower intensity beam surrounding the filament that exchanges power with the filament core [53]. The beam deflection reaches steady state in ~250 ms with amplitude of ~3.5 mrad. Note that the magnitude of the steady state deflection is reproduced by the model, while the time to steady state agrees with the hole parameter plots of Figure 2.3. As can be seen from Figure 2.3, the density hole width only increases by a factor of 2 or 3, while the depth and displacement increase by very large factors. The widening of the hole quickly stagnates on the bottom side as convection begins balancing thermal diffusion, preventing heat from spreading downward, as seen in Figure 2.2. As a result, the growth in the deflection angle with time is primarily driven by the density hole depth and displacement, and follows their relaxation timescale. 2.3.2 Density hole evolution We can estimate the relaxation timescale to steady state using fluid equations in steady state for conservation of momentum and energy. The physical picture is that the gas heated by the laser pulses begins to rise by buoyancy. The gas below the heated region then also rises, setting up a velocity flow field. One can show using the momentum equation that the vertical flow velocity is ultimately limited by viscous 47 drag to 2gL u    , where  (>0) is the mean decrement in gas density averaged transverse to the flow, L is the transverse scale length of the density hole, g is the gravitational acceleration, and  is the dynamic viscosity. We next consider the energy conservation equation in steady state and apply it to the region just below the localized laser-filament thermal source. Given the vertical flow velocity u , the density depression scale length in the vertical direction, L , adjusts itself until there is a balance between thermal transport and diffusion, giving B 0 5 2 Nk T uT L   , where 0 and 0T are the background gas density and pressure, N is the gas number density, T (>0) is the transverse average temperature increment, 0      is the relative density depression,  is the gas thermal conductivity, and Bk is Boltzmann’s constant. Combining these equations and using pressure equilibrium, 0 0 T T      , leads to an effective length scale   1 3 1 3 2 eff 2 B 2 5 L L L mN k g           for the system. Substituting effL for L in the expression for u gives 1 3 2 2 B 0.16 mg u N k          for the limiting flow velocity. For a given gas species, note the weak cube root dependence of u and effL on the experimentally controllable parameters of gas number density N and hole depth  , as well as on the relatively weak temperature dependence of  48 and  [54]. This makes the scale of these estimates robust over a wide range of conditions. For our experiment in Xe, 3 -1 -15.65 10 W m K     [54], 23μPa s   [54], 252.2 10 kgm   , 19 -37 10 cmN  at 2.7 atm, and we take a hole depth of 0.1 from Figure 2.3. This gives a limiting gas flow velocity of 115μm msu  and ff 150μmeL . The approximate time to reach this steady state is given by u t a  , where a g  is the acceleration from the buoyancy force. Examination of Figure 2.3 shows that in the early phase of hole displacement, 0.01 0.05  , giving 20 100mst  , a range of timescales in reasonable agreement with the approach to steady state shown in Figure 2.3 and Figure 2.4. In an alternative, essentially kinematic approach for estimating these space and time scales, we take an estimate of the initial height rise of the hole due to buoyancy, 2 rise 1 2 d gt and set it equal to the approximate hole radius as determined by thermal diffusion, thermal 2d t [9], the idea being that the thermal source location cannot move below the expanding hole width as the hole rises. Here, t is the elapsed time after the laser shutter is opened and B 2 5 Nk    is the thermal diffusivity. The result is 2 3 1 3 rise 4 t t g         and 1 3 2 rise 32 d g       . Using the above parameters gives rise 100mst and 1mmd , in reasonable agreement with the experiment and the prior scale estimates. 49 2.4 Conclusion In summary, we have examined the propagation of a high repetition rate filament pulse train in the transiently evolving gas that it heats. When a filament pulse train is suddenly initiated in a gas, the local density is reduced, driving buoyant motion of the gas. The upward drift of the gas density hole steers subsequent pulses downward. Eventually the buoyant motion is damped by viscous forces, establishing a steady gas flow field and hole density profile, and the downward beam steering stabilizes. Unlike conventional thermal blooming with high power CW lasers, the beam mode is preserved. Further, a simple ray model explains the transient deflection of the beam by the accumulated density hole, and simple fluid analysis explains the experimentally observed space and time scales for the density hole dynamics and beam steering. 50 Chapter 3: Demonstration and analysis of long-lived high power optical waveguiding in air 3.1 Overview In this chapter we explore the possibility of creating waveguides using the long timescale neutral gas hydrodynamic response following filamentation. In Chapter 2, we demonstrated that a high repetition rate (>1 kHz) filamenting beam deflects itself due to the cumulative thermal density hole generated by the heated gas from previous pulses in the pulse train. Such simple density holes will only defocus or deflect optical beams. However, by producing a ring of filaments, we show in this chapter that optical guiding structures can be generated in air in two ways: through collisions of single cycle acoustic waves launched by each filament in the array, or by the long timescale annular density hole that is left over after the acoustic waves have propagated away. Section 3.2 presents time resolved interferometry of the acoustic response following filamentation, which we use to characterize acoustic guides and resolve an open question in the literature regarding a waveguide created in the wake of a single filament. In Section 3.3, we demonstrate the generation of air waveguides and their guiding properties in a proof-of-concept experiment. In Section 3.4 we analyze the optical properties of thermal and acoustic air waveguides by performing a modal analysis. 51 3.2 Direct imaging of the acoustic waves generated by femtosecond filaments in air 3.2.1 Introduction In this section we investigate the time evolution of the acoustic response following femtosecond filamentation. Long-lived thermal depressions produced by short pulses were recently measured interferometrically [9], but those measurements did not have the sub-microsecond time resolution required to measure the full evolution of the acoustic response. Levi et al. recently used a delayed 150 ns pulse to probe the axially extended gas density perturbation produced by an ultrashort filamenting pulse in air [55]. They claimed a positive gas density perturbation on axis with a microsecond lifetime that supports optical guiding. Here, we present interferometric measurements of the acoustic response induced by filaments from single- and multi-mode beams, at repetition rates of 10 Hz and 1 kHz, and we compare our measurements to hydrodynamic simulations. We also perform simulations of the propagation of light in post-filament gas density profiles. For single filaments, at no delay do we find on-axis enhanced gas refractive index profiles that can support guided modes; such profiles are produced only by multi- filaments. Gas dynamics after deposition of energy by an ultrashort optical pulse has been discussed previously [9,10,56–58] and is reviewed in Section 1.4.6. We simulate the gas dynamics using the 1-dimensional radial Lagrangian hydrodynamical code detailed in Appendix A.2. 3.2.2 Single mode filament acoustic response 52 Simulation results assuming a 30 m FWHM Gaussian heat source of energy density 330mJ cm (consistent with the focal spot, typical plasma density and temperature of ~2x10 16 cm 3 and ~5 eV, and rotational heating of 4 meV/molecule [9,59]) are shown in Figure 3.1(a). The refractive index of air at standard atmospheric temperature and pressure, 4 0 1 2.7 10n    [60], was used to calculate the refractive index change from the gas density. The single cycle acoustic wave propagates outward at the speed of sound, 1340m s in air [61], and the amplitude of the density perturbation falls as 1 2r . The central density depression rapidly reaches its maximum depth by ~100 ns and then changes very little over microsecond timescales after the acoustic wave is launched. Short filaments were generated with 50 fs, 800 nm Ti:Sapphire pulses at 10 Hz to investigate single shot gas response and at 1 kHz to probe for a possible dependence on repetition rate [9,58]. A 532 nm, 7 ns laser probe pulse timed to the Ti:Sapphire laser is used in a counter-propagating geometry. The optical setup is similar to that in [9] and [55], with the important exception that here and in [9] we perform imaging interferometry for direct and unequivocal measurement of the evolving refractive index profile. We use folded wavefront interferometry [9,10,58], which measures the 2D phase shift ( , )x y of the probe beam, from which the refractive index perturbation n is extracted. Timing jitter between the two laser systems is <10 ns, which is unimportant for the longer timescales explored in this experiment. For a reliable measurement of ( , )x y , it is critically important to keep the interaction length of the probe short enough that refraction from the gas perturbation negligibly distorts the probe [62]. This is achieved by limiting the 53 filament length to ~2 mm by focusing the pump laser beam relatively tightly at f/30, where filament formation is dominated by the lens focusing. In the absence of refractive distortion, the change in probe phase due to the index perturbation n is     eff( , ) , , , ,0x y k n x y z dz k n x y L     , where k is the vacuum wavenumber of the probe, z=0 is the pump beam waist location, and Leff =2 mm. Figure 3.1. Acoustic gas dynamics following an ionizing pulse in air. The refractive index change n is shown at time T after the passage of the pulse. Blue: T = 80 ns; green: T = 280 ns; red: T = 480 ns. (a) Hydrocode simulation. The calculated n is plotted as a function of transverse coordinate. (b) Sequence of n profile lineouts obtained by interferometry of a 2 mm long plasma at 10 Hz. The measured refractive index change induced by the 10 Hz laser with a 96 J pulse (vacuum peak intensity 14 24 10 W cm  ) is shown in Figure 3.1(b) for a few delays following excitation. The density depression along the optical axis is clearly seen, as is the outwardly propagating single cycle acoustic wave. Good agreement is found between the experiment (Figure 3.1(b)) and the calculation (Figure 3.1(a)). 54 We repeated our measurements with the 1 kHz laser with pulse energy 60 J (vacuum peak intensity 14 22.5 10 W cm  ) and the same focusing geometry. The measured refractive index profile is shown in Figure 3.2. We found no qualitative differences between the 10 Hz and 1 kHz cases apart from the broad and shallow density hole left over from the previous pulse, where the depth of the hole at 1 ms is a few percent of ambient density [9]. Such a hole can negatively lens an optical pulse over an extended region [9,58], but has negligible effect on a tightly focused beam. We also observe that the hydrodynamic evolution of the gas initiated by individual pulses in the 1 kHz train is essentially unaffected by the presence of the preexisting density hole. Figure 3.2. Measured refractive index profiles induced by 60 J pump pulses at 1 kHz. The delays shown are with respect to passage of the pump pulse. The thermal gas density hole left between pulses is seen in the 60 s delay panel. 3.2.3 Multi-mode filament acoustic response 500 mm 0.2 ms 0.9 ms 1.4 ms 60.0 ms Dn −7 −6 −5 −4 −3 −2 −1 0 1 x 10 −5 55 Next, we generated double and octuple filaments using TEM01 and Laguerre- Gaussian * 04 0,4 0,-4LG LG LGi  modes at the focus. For a TEM01 focal mode, the pre-focused beam was passed through a half-pellicle angled so as to phase shift one half the beam by  with respect to the other half. For the * 04LG focus, with eight azimuthal intensity lobes, we use normal incidence reflection of the pre-focused beam from an eight segment mirror with alternating triangular segments recessed by /4~200 nm or /2. Using comparable methods, other groups have produced multi- lobed beams for filament generation [63,64]. Images of the acoustic response to multi-lobed modes are shown in Figure 3.3. Figure 3.3(a) and (b) show the measured response at 0.2 s and 0.5 s for a TEM01 mode. Figure 3.3(c) and (d) show simulations assuming that the gas dynamics is linear, so that two radially offset single filament results can be added. This is reasonable since the peak relative amplitude of the acoustic wave, 0.05    , is small. The approximation is borne out by the good agreement with the measurements. The measured acoustic response from the * 04LG mode is shown in Figure 3.3(e) and (f). The low intensity 8-lobe mode is shown in the Figure inset. The gas response to the * 04LG focus was simulated assuming a continuous ring heat source at t=0 with a Gaussian cross section at the ring location. The results, shown in Figure 3.3(g) and (h), are in excellent agreement with the experimental results. The ring does a very effective job of simulating the merged acoustic response to closely spaced azimuthal lobes. As seen in the experiment and simulation, at longer delays (Figure 3.3(f) and (h)), two single cycle sound waves propagate away from the beam axis. 56 This is because the lobes (or ring) launch both inward- and outward-directed acoustic pulses. The inward-directed annular acoustic wave collides with itself on axis and produces a very strong density increase of ~30% (Figure 3.3(e) and (g)). It then passes through the axis and re-emerges as a wave trailing the originally outward- directed wave. Movies of the gas evolution induced by single- and multi-filaments – both experiment and simulation- can be viewed at [65]. Figure 3.3. Measured and simulated air refractive index profiles induced by a TEM01 focus (vacuum peak lobe intensity 14 22.5 10 W cm  ) ((a)-(d)) and an * 04LG focus (vacuum peak lobe intensity 14 23 10 W cm  ) ((e)-(h)). The inset between (e) and (f) shows an expanded view of the focused * 04LG beam mode at low intensity. The 57 simulations in (g) and (h) assume a continuous ring heat source at t=0 with a Gaussian cross section at the ring location. 3.2.4 Identifying the guiding mechanism behind the single filament acoustic guide We have generated local density enhancements on-axis only in the case of gas dynamics initiated by a multi-lobed focus, as seen in Figure 3.3 for * 04LG and TEM01 modes, and with TEM11 modes (see [65]). As described above, in those cases, acoustic waves are launched from the lobe locations and superpose near the axis. However, for single filaments generated by a wide range of pump intensities and focal spot sizes, we always find that the central density (and refractive index) change is negative with respect to the surrounding gas. Such a profile would not be expected to guide an on-axis beam. But the results of Levi et al. [55] show apparent on-axis guiding of a 150 ns probe pulse at ~1 s delay in their post-single-filament gas density profile. What can explain this? In Levi et al., the filament length is ~300 mm, while for our interferometry measurements it is intentionally kept short at Leff ≈ 2 mm. There is no reason to believe that the radial hydrodynamics are significantly different in the two cases. Our simulation assumes an infinitely long heat source and it matches our short filament measurements extremely well. However, what remains to be considered are the details of probe pulse propagation in the extended gas density profile produced by a long filament. 58 Figure 3.4. Measurement of the axial dependence of energy deposition in a single filament and simulation of the propagation of a probe beam in the post-filament gas density profile. (a) Measured acoustic amplitude along the filament. The filamenting pulse propagates right-to-left. Error bars indicate the standard deviation of the signal amplitude over 100 shots at each position. The inset shows a sample acoustic trace measured by the microphone. (b-d) Beam propagation method simulations of a probe beam propagating left-to-right using hydrocode-generated gas index profiles weighted by the acoustic data of (a). The leftmost panels show index profiles associated with the maximum acoustic amplitude (near z = 42 cm) for the post-filament time delays in (b-d). To investigate this issue, we performed experiments and simulations for an extended interaction length; results are shown in Figure 3.4. A ~20 cm filament, here dominated by Kerr self-focusing and plasma defocusing, was generated by a 10 Hz, 2.8 mJ, 50 fs, 800 nm pulse focused at f/200. For a better understanding of a long filament’s effect on a probe beam, one must know the axial dependence of the induced refractive index perturbation. To characterize the axial profile of the energy 59 deposited by the filament, we used a sonographic technique [66]. A rail-mounted microphone is positioned ~3 mm away from the filament and scanned along it. The amplitude of the leading feature of the sound trace (see inset of Figure 3.4(a)) is plotted as a function of axial position in Figure 3.4. Our hydrodynamic simulations verify that the sound amplitude is proportional to the laser energy deposited. This information was fed into a paraxial beam propagation method simulation [67], the results of which are shown in Figure 3.4(b),(c), and (d) for probe injection delays of 0.4 s, 0.8 s, and 1.2 s, respectively. Here, the simulated gas density structure was built by weighting the results of the hydrocode by the axial slice-by-slice relative laser energy deposition measured in Figure 3.4(a). The structure is stationary on the timescale of the 7 ns probe pulse. Note the enhanced light intensity at the location of the expanding annular acoustic wave, where there is a positive refractive index perturbation. The annulus serves to trap probe light either when it overfills the entrance of the structure or when it is more tightly focused into the entrance. Either way, the density depression on axis repels probe light from the center of the structure, and a portion of the diverging beam is trapped by the annulus. The results of Figure 3.4 are for the probe beam overfilling the entrance, but simulation of a more tightly focused probe gives qualitatively similar results. Past the end of the filament (z=48 cm), the annular beam generates a strong interference maximum on axis within a few centimeters, producing the appearance of on-axis guiding if the plane containing the maximum is imaged onto a camera. The well- known Bessel beam [68] can be considered as an example of a self-interfering annular beam. The range of time delays for which an apparent on-axis mode is seen 60 depends on the position of the object plane and the diameter of the acoustic annulus, which is set by the sound speed and time delay. Beyond ~1 s, the annulus has moved sufficiently outward that the interference maximum shifts beyond the imaging system’s object plane. The acoustic wave’s amplitude, and thus its ability to trap light, also rapidly falls as it propagates outward. See Appendix A.3 for an explanation of the on-axis interference which occurs after light is launched from the annular acoustic ring. These effects create the appearance of an optimum temporal window for on-axis light trapping. Experimental probe images for our ~20 cm filaments show on-axis maxima whose behavior is consistent with this picture. For a 150 ns probe pulse [55], these effects should be present, though temporally smeared out. 3.2.5 Conclusions In conclusion, we have directly measured the sub-microsecond gas dynamics following ultrashort pulse single- and multi-mode filamentation in air, obtaining excellent agreement with hydrodynamic simulations. Accurate measurements depend on minimizing refractive distortion of the probe, which demands a short probe interaction length. For the post-single filament gas response, measurements and simulations show that while there is no positive on-axis refractive index enhancement at any delay, the annular sound wave can guide injected light, leading to the appearance of an on-axis interference mode of limited axial extent beyond the end of the filament. An array of filaments, however, can produce an on-axis refractive index enhancement owing to superposition of acoustic waves, and this structure may serve as an optical waveguide, as recently verified [10]. We have also recently shown that 61 an even more effective and much longer-lived waveguide can be generated using the relaxation of a filament array in the thermal diffusion regime (see Section 3.3) [10]. 3.3 Generation of long-lived optical waveguides in air 3.3.1 Background Femtosecond filaments are notable for their ability to deliver high intensity and high peak power at distances ranging up to several kilometers [35], enabling applications such as stand-off laser induced breakdown spectroscopy, backwards lasing through the atmosphere, and LIDAR [7,8,69]. Despite these applications, it remains a significant limitation that femtosecond filamentation cannot deliver high average power over long distances in a single tight spatial mode. This is due to the fact that for laser pulses with P ~ several Pcr , the beam will collapse into multiple filaments [70] with shot-to-shot variation in their transverse locations. For Pcr ~ 5-10 GW, this means that single filament formation requires pulses of order ~1 mJ. For a 1 kHz pulse repetition rate laser, this represents only 1 W of average power. Here, we demonstrate a method employing filaments that can easily supersede this limitation by setting up a robust, long range optical guiding structure lasting milliseconds. It opens the possibility for optical guiding of megawatt levels of average power over long distances in the atmosphere. The guiding structures demonstrated here have substantial potential for directed energy applications [71]. The generation of long-lived thermal guiding structures in air using filament arrays also has the potential to enhance other photonics applications in the atmosphere. For instance, they could be used to concentrate heater beams for remote atmospheric lasing schemes [72] or for inducing characteristic emission for standoff 62 detection of chemical compounds. Many remote detection applications rely on the collection of fluorescence [7,8,69], for these remote-sensing schemes, where detection over large distances may be desired, very little of the isotropically emitted fluorescence reaches the detector at a distance. The long-lived guiding structures demonstrated here could be used as an effective collection lens, enhancing the signal. They may also find use in atmospheric laser communication [73]. Finally, they might also be used to enhance and control the propagation of an injected ultrashort filamenting pulse [74], similar to what has been done with a permanent refractive index structure in glass [75] and recently with the plasma from an array of filaments in air [76]. We note that there has been much recent work on using the refractive index of the plasma generated by an array of filaments to form guides for microwaves [77,78] and nanosecond optical pulses [79]. We emphasize that the guiding we demonstrate here does not use the optical response of the plasma – rather, it uses the ~106 times longer duration hydrodynamic response of the gas after heating by the filaments. 3.3.2 Gas hydrodynamics initiated by femtosecond filaments Recently we found that a femtosecond filament, starting at an electron temperature and density of a few eV and a few times 10 16 cm -3 [24], acts as a thermal source to generate long-lived gas density hole structures that can last milliseconds and dissipate by thermal diffusion [9]. In air, additional heating can occur from molecular excitation [9], with peak deposited energy density in the plasma and molecules of as much as ~100 mJ/cm 3 . Owing to the finite thermal conductivity of the gas, the initial energy invested in the filament is still contained in a small radial zone, but it is 63 repartitioned into the translational and rotational degrees of freedom of the neutral gas. The result is an extended and narrow high pressure region at temperatures up to a few hundred K above ambient. In air, this pressure source launches a radial sound wave ~100 ns after the filament is formed. By ~1 s, the gas reaches pressure equilibrium with an elevated temperature and reduced gas density in the volume originally occupied by the filament, after which the ‘density hole’ decays by thermal diffusion on a few millisecond timescale [23]. Figure 3.5. Gas dynamics following a single filament in air. (a) Interferometric measurement of the refractive index change following a short pulse as a function of the time delay of the probe pulse. (b) Hydrodynamic simulation, assuming a 60 m FWHM Gaussian heat source of peak initial density 32 mJ/cm 3 . The full dynamics are clearly seen in Figure 3.5(a), which presents a time resolved measurement of the 2D density hole evolution (expressed as air refractive index shift) of a short air filament from nanoseconds through microseconds after filament formation. A 1D radial fluid code simulation, described in Appendix A.2, is shown in Fig. 3.5(b) for comparison and the results are in excellent agreement with the measurements. The experimental results verify that the density hole first deepens over tens of nanoseconds, and launches a sound wave which propagates beyond the 64 ~200m frame by ~300ns. By ~ 1-2 s, pressure equilibrium is reached and the hole decays by thermal diffusion out to millisecond timescales. 3.3.3 Experimental setup 65 (a) (b) Figure 3.6. Generation of a filament array using half pellicles. (a) A 55 fs, 800 nm, 10 Hz pulsed laser is used to generate an array of four filaments. A pulse propagates through two orthogonal half-pellicles, inducing π phase shifts on neighboring quadrants of the beam, and then are focused to produce a 4-filament with a TEM11 mode (actual low intensity image shown). A 7 ns, 532 nm 10 Hz pulsed laser counter- propagates through the filament and is imaged either directly onto a CCD for guiding 66 experiments or through a folded wavefront interferometer and onto a CCD for interferometry. (b) Rayleigh scattering as a function of z with a bi-filament produced by a single half pellicle (the bi-filament far-field mode is shown in inset). The bottom row shows burn patterns produced by a 4-filament produced by two orthogonal half pellicles. A =532 nm, 7 ns duration beam counter-propagates along a femtosecond filament structure generated by a 10 Hz Ti:Sapphire laser system producing =800 nm, 50 fs pulses up to 100 mJ. The optical arrangement is similar to our earlier experiment of ref. [9]. Here, the 532 nm pulse serves as either a low energy interferometric probe of the evolving gas density profile, using a folded wavefront interferometer, or as an injection source for optical guiding in the gas density structure. The transverse gas density profiles shown in Fig. 1(a) were obtained using the 532 nm pulse as an interferometric probe of a single short ~2 mm filament. The short filament length is essential for minimizing refractive distortion of the interferometric probe pulse [62]. The 2D density profiles were extracted from the interferograms as described in ref. [9]. The delay of the 532 nm probe/injection pulse is controlled with respect to the Ti:Sapphire pulse with a digital delay generator. The pulse timing jitter of <10 ns is negligible given the very long timescale gas evolution we focus on. For the injection experiments, up to 110 mJ is available at 532 nm. The Rayleigh side-scattering image of Figure 3.6(b) was obtained by concatenating multiple images from a low noise CCD camera translated on a rail parallel to the filament. The images were taken through an 800 nm interference filter. We note that at no probe delay do we see an on-axis refractive index enhancement that might act as a waveguiding structure and explain a recent report of filament guiding [55], an issue further discussed in ref. [80] and Section 3.2. At the 67 longer delays of tens of microseconds and beyond, the thermal gas density hole acts as a negative lens, as seen in our earlier experiments [9]. 3.3.4 Multi-filament-induced guiding structure Although a single filament results in a beam-defocusing gas density hole, a question arises as to whether a guiding structure can be built using the judicious placement of more than one filament. We tested this idea with a 4-lobed focal beam structure using two orthogonal ‘half-pellicles’. As seen in Figure 3.6(a), the pellicles are oriented to phase-shift the laser electric field as shown in each near-field beam quadrant. Below the filamentation threshold, the resulting focused beam at its waist has a 4-lobed intensity profile as shown, corresponding to a Hermite-Gaussian TEM11 mode, where the electric fields in adjacent lobes are  phase shifted with respect to each other. Above the threshold, the lobes collapse into filaments whose optical cores still maintain this phase relationship and thus 4 parallel filaments are formed. As a demonstration of this, the top panel of Figure 3.6(b) shows an image of the Rayleigh side-scattering at 800 nm from a 2-lobed filament produced by a single half pellicle, indicating that the  phase shift is preserved along the full length of the filament. This image was obtained by concatenating multiple images from a low noise CCD camera translated on a rail parallel to the filament. The images were taken through an 800 nm interference filter. The bottom panel shows burn patterns taken at multiple locations along the path of a ~ 70 cm long 4-lobe filament used later. For the 70 cm 4-filament, the filament core spacing is roughly constant at ~300 m over a L~ 50 cm region with divergence to ~ 1 mm at the ends. 68 The effect of a 4-filament structure on the gas dynamics is shown in Figure 3.7, a sequence of gas density profiles measured for a short ~2 mm filament (produced at f/35) to minimize refractive distortion of the probe beam. The peak intensity was 14 2<10 W cm , typical of the refraction-limited intensity in more extended filaments, so we expect these images to be descriptive of the gas dynamics inside much longer filaments. Inspection of the density profiles shows that there are two regimes in the gas dynamical evolution which are promising for supporting the guiding of a separate injected laser pulse. A shorter duration, more transient acoustic regime occurs when the sound waves originating from each of the four filaments superpose at the array’s geometric center, as seen in panel (a) of Figure 3.7, causing a local density enhancement of approximately a factor of two larger than the sound wave amplitude, peaking ~80 ns after filament initiation and lasting approximately ~50 ns. A far longer lasting and significantly more robust profile suitable for guiding is achieved tens of microseconds later, well after the sound waves have propagated far from the filaments. In this thermal regime, the gas is in pressure equilibrium [9]. As seen in panels (c) and (d) of Figure 3.7, thermal diffusion has smoothed the profile in such a way that the gas at center is surrounded by a ‘moat’ of lower density. The central density can be very slightly lower than the far background because its temperature is slightly elevated, yet it is still higher than the surrounding moat. The lifetime of this structure can be several milliseconds. In both the acoustic and thermal cases, the diameter of the air waveguide “core” is approximately half the filament lobe spacing. A movie of the 4-filament-induced gas evolution is provided in the supplementary material [81]. 69 Figure 3.7. Interferometric measurement of the air density evolution induced by a 4- filament. (a) The acoustic waves generated by each filament cross in the middle, generating a positive index shift, producing the acoustic guide. (b) The acoustic waves propagate outward, leaving behind a density depression at the location of each filament. (c) The density depressions produce the thermal guide, with a higher central density surrounded by a moat of lower density. (d,e) The density depressions gradually fill in as the thermal energy dissipates. A movie of the 4-filament-induced gas evolution is provided in the supplementary material [81]. 3.3.5 Fiber analysis of air waveguides Having identified two potential regimes for optical guiding, a short duration acoustic regime, and a much longer duration thermal regime, it is first worth assessing the coupling and guiding conditions for an injected pulse. We apply the fiber parameter V for a step index guide [82] to the air waveguide,     1 2 1 2 2 2 co clco cl 2 2 2a a V n nn n        , where the effective core and cladding regions have refractive indices co,cl 0 co,cln n n  , 0n is the unperturbed background air index ( 4 0 1 2.77 10n    at room temperature and pressure [60]), con and cln are the (small) index shifts from background at the core and cladding, and the core diameter is 2a, taken conservatively at the tightest spacing of the filament array. The numerical 70 aperture of the guide is 2 V NA a    . Because accurate density profile measurements are restricted to short filaments, we use the results of Figure 3.7 and apply them to much longer and wider-lobe-separated filaments that are inaccessible to longitudinal interferometry owing to probe refraction. As typical filament core intensities are restricted by refraction (‘intensity clamping’) to levels 14 -210 W cm  [35], we expect that the measurements of Figure 3.7 apply reasonably well to longer filaments and different lobe spacings. For the acoustic guide, we used a filamenting beam with lobe spacing of 150 m, so 2a ~75 m. Using co 0 0.05 1 n n   , and cl 0 0.02 1 n n   from Figure 3.7 then gives V~ 2.8 (> 2.405) and NA~ 6.310-3, indicating a near-single mode guide with an optimum coupling f-number of 0.5 # 80f NA  . For the long thermal guide, we used a filamenting beam with lobe spacing of ~300 m, so 2a~150 m. From Figure 3.7, the core index shift is nco ~0 and the cladding shift is the index decrement at the moat, cl 0 0.02 1 n n    , giving V~2.9, corresponding to a near-single mode guide with NA~3.210-3, corresponding to f/# ~160. For a more complete analysis of the modal properties of the thermal and acoustic guides, see Section 3.4 3.3.6 Injection and guiding experiments The experimental setup is shown in Figure 3.6(a). An end mode image from injection and guiding of a low energy =532 nm pulse in the acoustic waveguide produced from a 10 cm long 4-filament is shown in Figure 3.8. In order to 71 differentiate between guiding and the propagation of the unguided beam through the fully dissipated guide at later times (>2 ms) we define the guiding efficiency as g ug tot ug E E E E   where Eg is the guided energy within the central mode, Etot is the total beam energy and Eug is the fraction of energy of the unguided mode occupying the same transverse area as the guided mode. Figure 3.8. Demonstration of guiding of 7 ns, =532 nm pulses in 70 cm long acoustic and thermal air waveguides produced by a 4-filament. The panel in the upper left shows the probe beam, which is imaged after the filamentation region, with and without the filament. The time delay of the probe was 200ns, which is in the acoustic guiding regime. The effect of the thermal waveguide, the shadow of which can be seen in the image in the top center (with a red dashed circle showing the position of the lower density moat), is shown in the bottom row, where the probe beam is imaged after the exit location of the air waveguide with and without the filamenting beam. The coupling efficiency vs. injected pulse delay is shown in the upper right. Peak energy guided was ~110 mJ. Best coupling occurred at an injection delay of ~200 ns and f/# >100, with a peak guided efficiency of 13%, although the guides were not stable on a shot-to-shot basis. Efficient guiding in the acoustic regime takes place over an injection delay 72 interval of only ~100 ns, consistent with the time for a sound wave to cross the waveguide core region, s 100ns a c , where 2a=75 m and 4 -1 s 3.4 10 cm sc   is the air sound speed [54]. We found that for longer filaments with wider lobe spacings, the acoustic guides were even less stable. Unless the 4-filament lobes were well balanced in energy and transverse position, the sound wave superposition would not form a well-defined air waveguide core. This is why a shorter 10 cm filament was used for the acoustic guide experiment. While the acoustic superposition guide is a promising approach, future experiments will need filaments generated by very well-balanced multi-lobe beam profiles, an example of which is seen in ref. [80]. By comparison, the thermal guides were far more robust, stable, and long- lived. Results from the thermal guide produced by a 70 cm long 4-filament are also shown in Figure 3.8, where optimal coupling was found for f/#=200, in rough agreement with the earlier fiber-based estimate. An out of focus end mode image (not to scale) is shown to verify the presence of the thermal guide’s lower density moat. Here we note that owing to the much greater lobe spacing of its long 4-filament, the thermal guide of Figure 3.8 lasts much longer (milliseconds) than that from the short 4-filament of Figure 3.7 (~10 s). Guided output modes as function of injection delay are shown imaged from a plane past the end of the guide, in order to minimize guide distortion of the imaging. These mode sizes are larger than upstream in the guide where 4-filament lobe spacing is tighter, but where we are unable to image reliably. We injected up to 110 mJ of 532 nm light, the maximum output of our laser, with 90% energy throughput in a single guided mode. This corresponds to a peak guiding efficiency of 70%. Guiding efficiency vs. injected pulse delay is plotted in Figure 3.8. 73 As seen in that plot, peak guiding occurs at ~600 s and persists out to ~2 ms where the guiding efficiency drops to ~15%. Based on the guide core diameter of 2a~150 m and the portion of the filament length with constant lobe spacing, L~50 cm, the guided beam propagates approximately 2 15 L a   Rayleigh ranges. A movie of the thermal waveguide output beam, during real time rastering (at 10 Hz) of the injected beam across the guide entrance, is shown in the supplementary material of ref. [10]. 3.3.7 Simulations of waveguide development and guiding Owing to the linearity of the heat flow problem, the evolution of the 4-filament- induced density structure in the thermal regime can be calculated by finding the solution ( , , )T x y t to the 2D heat flow equation, 2tT T   , for a single filament source located at ( , ) (0,0)x y  and then forming 4 4 1 ( , , ) ( , , )j j j T x y t T x x y y t     , where ( , )j jx y are the thermal source locations in the 4-filament. Here pc    , where  and pc are the thermal conductivity and specific heat capacity of air. To excellent approximation, as shown in ref. [9], ( , , )T x y t is Gaussian in space. Invoking pressure balance, the 2D density evolution is then given by       2 2 2 4 0 0 4 b 2 2 1b 0 0 , exp 4 4 j j j x x y yT R N r t N T R t R t                     , where 0R is the initial 1/e radius of the temperature profile of a single filament and 0T is its peak value above bT , the background (room) temperature. Using 0 50μmR  , 0 15KT  , 74 2 -10.21cm s   for air [9], and source locations separated by 500 m, approximating our 70 cm 4-filament, gives the sequence of gas density plots shown in the upper panels of Figure 3.9, clearly illustrating the development and persistence of the guiding structure over milliseconds. Figure 3.9. Simulation of the evolution of and guiding in a thermal air waveguide. The top row shows the index of refraction shift produced by the 4-filament-induced temperature profile as a function of time. The bottom row shows a BPM simulation of the guided laser beam profile at the end of a 70 cm waveguide produced by the 4- filament-induced refractive index change. The propagation of the 532 nm beam in the waveguide was simulated in the paraxial approximation using the beam propagation method (BPM) [67]. The calculated intensity at the output of the waveguide is shown in the lower panels of Figure 3.9. At early delays <100 s, characteristics of a multimode waveguide are observed in the simulation, including mode beating. At later times, as the refractive index contrast decreases, the propagation is smoother, indicating single mode behavior, consistent with the estimates in Section 3.3.5 using the fiber parameter. The simulation is in reasonable agreement with the experimental results. Axial nonuniformity in the waveguiding structure could explain the absence of four fold 75 symmetry in the experimental data, whereas it is pronounced in the simulations. For a more complete assessment of the modal properties of the thermal waveguide as well as the acoustic guide, see Section 3.4. 3.3.8 Conclusions and potential future application to guide high average power beams We have demonstrated the generation of very long-lived and robust optical waveguides in air, their extent limited only by the propagation distance of the initiating femtosecond filament array and the axial uniformity of its energy deposition. Assuming a sufficiently uniform filament, this is ultimately determined by the femtosecond pulse energy absorbed to heat the gas. Based on a single filament diameter of ~100 m, an electron density of ~31016 cm-3 [24], ionization energy of ~10 eV per electron, and 5 meV of heating per air molecule [9], approximately 0.5 mJ is needed per meter of each filament. With a femtosecond laser system of a few hundred millijoules pulse energy, waveguides hundreds of meters long are possible. What is the optical-power-carrying capacity of these guides? For high peak power pulses, the peak guided intensity,  -2p W cmI  , is limited by self-focusing in the waveguide core. For 1μm  guided light, the approximate B-integral condition is 14 2 p g ~ 2 10 W cmI L    where gL is the guide length (cm) and the nonlinear index of air is obtained from Ref. [16]. For our ∼100 cm waveguide, this gives 12 2 p g ~ 2 10 W cmI L    , corresponding to a pulse energy of ∼4 J for the 150 μm diameter core, assuming a 10 ns pulse. At higher repetition rates (> 1 kHz) and over longer ranges (hundreds of meters), heating from non-linear absorption mechanisms 76 such as stimulated Raman scattering (SRS) may also become relevant for guided energetic nanosecond pulses. However, the real utility of these air waveguides, in the thermal formation regime, derives from their extremely long millisecond-scale lifetime. This opens the possibility of guiding very-high-average powers that are well below the self-focusing or SRS thresholds. We now consider the robustness of our filament-induced waveguides to thermal blooming [49,71] from molecular and aerosol absorption in the atmosphere. For thermal blooming, we consider the deposited laser energy which can raise the local gas temperature by a fraction  of ambient, g 1.5 P t p A     , where Pg is the guided laser power, t is the pulse duration,  is the absorption coefficient, A is the waveguide core cross sectional area, and p is the ambient pressure. Thermal blooming competes with guiding when  is approximately equal to the relative gas density difference between the core and cladding. In our measurements of the thermal air waveguide, the typical index (and density) difference between the core and cladding is of the order of ~2% at millisecond timescales. Taking =0.02, p=1 atm, and =2×10-8 cm-1 [71], gives g 5 21.5 10 J cm P t A      as the energy flux limit for thermal blooming. For example, for a 1.5 mm diameter air waveguide core formed from an azimuthal array of filaments [80], the limiting energy is g 2.7kJP t . Note that we use a conservative value for  at ~1 m which includes both molecular and aerosol absorption for maritime environments [71], which contain significantly higher aerosol concentrations than dry air. If a high power laser is pulsed for t~2 ms, consistent with the lifetime of our 10 Hz-generated thermal waveguides, the peak 77 average power can be 1.3 MW. It is possible that in such environments, air heating by the filament array itself could help dissipate the aerosols before the high power beam is injected, raising the thermal blooming threshold and also reducing aerosol scattering. An air waveguide even more robust against thermal blooming and capable of quasi-continuous operation may be possible using a kHz repetition rate filamenting laser. We have already shown that the cumulative effect of filamenting pulses arriving faster than the density hole can dissipate leads to steady state hole depths of order ~10% [9]. 3.4 Modal analysis of air waveguides 3.4.1 Introduction In this section, we explore the modal and light transport properties of air waveguide structures using analytical and numerical methods. The optical properties of air waveguides depend on the evolving air density profile, which is determined by the axial and transverse distribution of energy deposited in the gas by the filament or filament array. For typical pulse durations of ~40-100 fs, approximately 325mJ cm is deposited over the ~50 µm radius filament core [9,59,80] by plasma generation and molecular rotational excitation. The hydrodynamic response of the gas is again simulated using our 1D Lagrangian hydrocode where we take the initial energy deposition to be 325mJ cm . Details of the code can be found in Appendix A.2, and in earlier work [9,10,80]. We consider the two regimes of guiding demonstrated in Section 3.3.6, acoustic and thermal guiding. Acoustic guiding, which typically occurs over a microsecond timescale interval [10,80,83], works by confining light in the positive 78 density (refractive index) crest of the single cycle annular acoustic wave launched from a single filament [55,80], or in the enhanced density peak resulting from collision of acoustic waves from multiple filaments [10,80] (see Figure 3.10). In the much longer lived, millisecond scale thermal guiding regime [9,10,83], a ring of  4 filaments leaves a cladding ‘moat’ of diffusively merged density holes surrounding a near ambient air density core (see Figure 3.11). 3.4.2 Modelling In order to discuss mode structure, we consider idealized air waveguides whose transverse refractive index profile is axially invariant. While real air waveguides have axial variation [10,80,83] the scale is over many Rayleigh ranges of the guided beam. We also neglect any absorption or scattering losses, which can be added with an exponential attenuation factor. To find the modes ( , ) i zA z e r , where A is a slowly varying field amplitude and r and z are transverse and axial coordinates, we consider the paraxial scalar wave equation 2 2 2 2 (r , ) 2 [ ( (r ) )] (r , ) 0 A z i k n A z z             (3.1) where k is the vacuum wavenumber, 0( ) ( )n n n  r r is the air refractive index (assumed dispersionless for our range of considered wavelengths), ( )n  r is the filament-induced refractive index profile change, and  is the waveguide wavenumber. Use of the scalar wave equation is well justified, as the air waveguides created in the wake of filamentation are all well within the weak guiding 79 regime [82,84]; index shifts are of order 5 0 10 n n  . For later use in our analysis, we define an effective mode index n by  0k n n   . The type of mode supported depends on the profile details of the guide. Profiles with a region or regions of positive refractive index with respect to the surrounding air can support strictly bound modes [82], which propagate with no transverse leakage loss. This applies to both the single and multifilament acoustic guides, where there are regions of increased air density with respect to the ambient density. In contrast, thermal guides support only leaky modes [82] because the index perturbation from the density holes is strictly negative. Leaky modes have an oscillatory component at large r that corresponds to the leakage of energy out of the guide as the light propagates [85]. However, we will show that the propagation losses per unit length for these guides can be quite small. To find the bound modes of air waveguides that can support them, we set 0 z    in Equation (3.1), leaving the transverse Helmholtz equation which we solve numerically. For azimuthally symmetric guides, the solutions are indexed by radial (p) and azimuthal (m) mode numbers. To find the fundamental leaky mode of an arbitrary leaky structure, we apply the beam propagation method (BPM) [67] to Equation (3.1) (for 0 z    ) using a suitable initial guess ( ,0)A r and sufficient propagation along z that the solution, up to an exponential decay along z, approaches 80 steady state for small A z   . This method is equivalent to shedding any higher order leaky modes. 3.4.3 Modal structure of acoustic air waveguides Figures 3.10(a) and (b) show experimental and simulated index profiles for the single filament-induced acoustic guide at time t=500 µs after the filament. The amplitude and radial scale (or wavelength) of the outgoing density perturbation is imposed by the strength and diameter a of the initiating filament energy deposition. Figures 3.10(f) and (g) show calculated bound modes at =532 nm for (p,m) = (0,0) and (0,5), confirming that the mode conforms to the acoustic wave peak [55,80]. At typical levels of energy deposition in filaments, the density modulation is not deep or wide enough to support higher order (p >0) radial modes. As the acoustic wave expands outward, higher order azimuthal modes with cosm dependence are allowed, with the highest azimuthal mode number scaling as 0.6 max 0 R m R , as seen in Figure 3.10(i), where 0R is the size of the filamenting core, sR c t is the ring radius, cs is the acoustic speed and t is time after the filament energy deposition. Prior experiments and simulations of trapping in single filament acoustic guides have shown a useful window of ~1 µs over which optimally positioned Gaussian beams can be efficiently coupled [55,80]. 81 Figure 3.10. Mode properties at =532 nm of single and octo-filament induced acoustic guides. Measured air refractive index profiles for single (a) and octo-filament excitation (c), along with corresponding hydrocode-simulations ((b) and (d)). Panel (e) shows an index profile simulation for the larger transverse scale typical of axially extended air waveguides. Calculated modes for the single filament-induced acoustic annular guide are shown in (f) p=0, m=0 and (g) p=0, m=5. The lowest order mode for the guide of (e) is shown in (h). Panel (i) shows the maximum azimuthal mode number mmax vs delay for the single filament-induced annular guide. Points: simulation, red curve: best fit. In the multifilament acoustic guide, each of the filaments, equidistant from a common center, launches a single cycle acoustic wave. The interference maximum produced when the waves meet on axis typically lasts s 1μs a c   , a timescale a) single b) nx10-5 500 m -6 -4 -2 0 c) in d e x s tr u ct u re octo d) e) 500 m -2 0 2 4 6 f) E -f ie ld m o d e s single g) single h) a.u. octo500 m -1 0 1 0 0.5 1 1.5 0 2 4 6 8 time in s m a x a zi m u th a l n u m b e r i) single 82 scale set by the sound speed and acoustic wavelength (here ~300 ns). Figures 3.10(c) and (d) show an interferometric measurement of an octo-filament-induced acoustic guide at the moment of peak central index enhancement at delay ~200 ns and its hydrocode simulation using a ring-shaped pressure source [80]. Guides as in Figure 3.10(c), amenable to longitudinal interferometry, are shorter with tighter transverse spatial features [10,80]. A hydrocode simulation more representative of our axially extended structures used for guiding [10,83] is shown in Figure 3.10(e) for an initial ring pressure source of radius 200 m with the same total energy deposition as for 8 filaments. The p=0, m=0 mode for this guide is shown in Fig 3.10(h), peaking on axis, in contrast to the fundamental ring-shaped mode of the single filament-induced acoustic guide (Figure 3.10(c)). For typical parameters, such multifilament acoustic guides only support the lowest order mode. 3.4.4 Fundamental mode and leakage rate for thermal air waveguides 83 Figure 3.11. Typical refractive index profiles and modes of thermal air waveguides. (a) Interferometrically measured refractive index profiles, (b) simulated time evolution of eight-lobe filament with 500 μm ring radius, (c) simulated time evolution of ring heat source with 500 μm ring radius and initial energy deposition equal to the eight-lobe filament case, (d) calculated fundamental leaky modes at =532 nm of the octo-lobe and ring air waveguides shown in (b) and (c). As discussed earlier, for thermal guides generated by a ring array of filaments, leaky guiding occurs, with the field energy contained in the ambient density core bounded by the lower density moat where the filaments were located. Row (a) of Figure 3.11 shows interferometric measurements of short octo-filament arrays for several times during the evolution of the thermal guide. The three panels resolve distinct stages of the evolution: At t = 2 μs the density holes from individual filaments are still distinct, at t = 10 μs the holes have merged through thermal diffusion to form t = 2 s a) expt. t = 10 s t = 50 s 150 m n -6 -4 -2 0 x 10 -5 t = 500 s b) octo sim. t = 1000 s t = 1500 s 1 mm n -2 -1 0 x 10 -6 t = 500 s c) ring sim. t = 1000 s t = 1500 s 1 mm n -2 -1 0 x 10 -6 d) octo guide t = 500 s a.u. ring guide 1 mm t = 500 s 0 0.5 1 84 a continuous ring, and at t = 50 μs the holes have diffused to the center and have washed-out the guide. Figures 3.11(b) and (c) show larger transverse scale simulations consistent with longer multi-filaments used in guiding experiments [10,83]. As can be seen, the relevant timescales increase with the transverse scale size, with guide washout now occurring at > 1.5 ms. In Fig. 3.10(b) the linearity of gas hydrodynamics for small density perturbations allows the summation of single filament results from our 1D hydrocode (see ref. [80] or Appendix A.2) to obtain an excellent approximation of the 2D response to the octo- filament. Figure 3.11(c) shows the efficacy of treating the octo-filament array as a continuous ring heat source with the same total initial energy. In contrast with the acoustic regime, the thermal regime produces a much longer lasting guide whose lifetime 2 thermal 1ms 4 R    [9] is set by the thermal diffusivity of the gas ( 2 -119μm μs   for air at standard conditions [61]) and the transverse length scale of the guide (here R=500 μm). The fundamental leaky modes for the eight-lobe and idealized ring thermal guides, shown in Figure 3.11(d), are close to Gaussian in shape and nearly identical. Even at the early time t=500 μs, where there is a marked difference between the eight-lobed and ring index structures (Figures 3.11(b) and (c)), the overlap integral between the two modes is >0.995. We find that the propagation leakage loss from the ring guide of Figure 3.11(c) closely approximates that for the octo-thermal guide of Figure 3.11(b). Computed loss rates shown in Figure 3.12(a) for the ring guide (in red) and the octo- guide (in blue) are nearly identical except at early times, when the gaps between the 85 octo-lobes cause higher loss. Approximate expressions for the axial loss rate Im( )  and effective index n of the lowest order mode can be derived assuming a step index decrement n (<0) at the moat of width r, cl cl 2 co cl 2 cl co 1 r r e k r e              (3.2a)   2 2 2 0 0 co2 2 co 1 1 1 1 ln(2 ) 2 2 n n n n k r k r        (3.2b) Where cor is the core radius, and   1 2 co 02k n n   and    1 2 cl 02k n n n   are propagation constants in the core and cladding. Since the lowest order leaky mode is nearly Gaussian (Figure 3.11(d)) we used a variational equation with Gaussian trial function to derive Equations (3.2). For details regarding derivation of the loss rate  and effective index n see Appendix A.4. At early times up to ~500 μs, propagation loss decreases as there is less transverse leakage through a broad and shallow cladding region than through a narrow and deep one. Beyond ~500 μs the trend reverses as the density hole diffusion reduces the central density, decreasing the index contrast between the core and cladding. Loss rates vs. delay for a range of wavelengths are plotted in Figure 3.12(b). 86 Figure 3.12. (a) Energy loss at =532 nm in dB/100 m for the octo-lobe index and ring index structures of Fig. 2 (blue and red respectively) and from the expression for  in Eq. (2) (black). (b) Loss in dB/100 m vs. time for a range of wavelengths (color scale in log base 10). 3.4.5 Coupling efficiency Figure 3.13 shows coupling efficiency into various guides for a varying Gaussian beam waist diameter, with the waist located at the waveguide entrance and centered on the guide axis. Coupling efficiency is defined as the overlap integral between the Gaussian beam and the modes of the guide. For the single filament acoustic guide, coupling only occurs into the lowest order mode, as azimuthal symmetry of the Gaussian beam prevents coupling into any higher order modes. Highest coupling efficiency is achieved at early times, but the presence of the deep density hole on axis requires the coupled mode to be large. At longer times, the optimal waist size increases linearly in time, following the acoustic wave location sr c t . However, coupling efficiency falls as 1t because the beam area increases as 2r to overlap an acoustic ring area increasing as r . Maximum coupling efficiency to the multifilament acoustic guide occurs when the inward propagating acoustic waves en e rg y lo ss ( d B /1 0 0 m ) w av el en gt h in μ m time in ms 0 1 2 101 10-1 10-3 0 1 2 1.5 2 0.5 1 2 0 -2 -4 a) b) 87 collide on axis, which occurs just beyond 0.5 μs in Figure 3.13(b). For a waist of 70 μm, the coupling efficiency is >85% over a duration of ~ 250 ns. For the octo- filament-induced thermal guide (Figure 3.13(c)), the lowest order mode is nearly Gaussian; we find that optimal coupling efficiency is >95% for approximately the first millisecond. As one can see, there is a very broad range of nearly optimal waist sizes with a full width at half maximum spanning ~400 μm. Figure 3.13. Coupling efficiency of a Gaussian beam at focus as a function of time and waist size for the (a) single filament acoustic guide, (b) octo-filament acoustic guide, and (c) octo-filament thermal guide. 3.4.6 Conclusions In conclusion, we have found and characterized optical modes of the distinct types of air waveguides that are formed in the wake of femtosecond filaments or filament arrays. The acoustic wave from a single filament supports only ring-shaped optical modes, with an increasing number of azimuthal modes with time. However, the coupling efficiency for an injected Gaussian beam falls as 1t , fading away by ~1 s. The multifilament acoustic guide primarily supports lowest order modes peaked on axis, but likewise has an efficient coupling lifetime of ~1 s. The thermal air waveguides studied here support near-Gaussian leaky modes. We examined their w a is t in m m a) 0.07 1 1.5 0.7 0 time in s b) 0.33 0.8 c) 0 2000 0 0.5 1 88 leakage loss as a function of delay, wavelength, and guide characteristics. Gaussian beam coupling of >95% efficiency is possible over the millisecond scale lifetime of the guide. 89 Chapter 4: Spatiotemporal Optical Vortices 4.1 Overview In this chapter we study spatiotemporal optical vortices (STOVs), and their connection to optical collapse arrest, which is a universal process in all self-focusing scenarios. STOVs are a new kind of optical vortex, which exhibit phase winding and the flow of electromagnetic energy in a spatiotemporal plane. Section 4.2 introduces the concept of an optical vortex, briefly reviews existing literature on optical vortices, and then provides an example construction of an electric field with an embedded STOV. Section 4.3 presents experiments, supported by theory and numerical simulation, constituting the first experimental detection of a STOV. We establish that STOV generation is fundamental to the optical collapse process, occurring in all optical pulses undergoing optical collapse. Finally, Section 4.4 considers the pulse evolution and angular momentum of linearly propagating STOV-carrying pulses. 4.2 Introduction Vortices – localized regions in which the flow of some quantity such as mass or electromagnetic energy circulates about a local axis – are a common and fundamental element of classical [86] and quantum fluids [87–89] as well as optics [90]. In quantum fluids, the circulating quantity is the spatial atomic probability density; in optics it is the electromagnetic energy density. Both densities are expressed as the magnitude squared of a complex scalar field iue   derivable from a Schrödinger-like equation (SE), where the vortex circulation is l c d   , u and  are real scalar fields, and the integral is on a closed contour about the local 90 axis. A non-zero value of  implies a discontinuity or ‘defect’ in . Furthermore, the single-valuedness of fields demands that the circulation be quantized, 2m  , where m, an integer, is called the ‘topological charge’. Because remains constant as one shrinks the contour around the local axis, or vortex core, the phase  becomes undefined at the core center, (‘phase singularity’) where the field magnitude is necessarily zero. In optics, vortices have been heavily investigated for decades, but as far as we know, the types that have been studied experimentally are entirely those that can be supported by monochromatic waves, that is, they have phase and energy circulation purely in the spatial domain. In this chapter, we present experimental evidence, backed by theory and simulation, for a dynamic phase vortex that mediates energy flow spatially and temporally. The spatiotemporal character of the vortex dictates that energy flow near the vortex core is saddle or spiral depending on the sign of dispersion of the propagation medium. Further, we show that vortex formation is an inherent feature of nonlinear collapse arrest, and that vortex dynamics are associated with changes to the envelope of a nonlinearly propagating pulse. In an early, influential publication, Nye and Berry introduced the concept of dislocations (field nulls) to wave theory and gave many examples [91]. A well-known example of a vortex in linear optics is the linearly polarized Laguerre Gaussian (LG) mode lpLG ( , ,z)r  of integer radial index p0 and non-zero integer azimuthal index l, with azimuthal dependence exp( )il ). This mode has an on-axis field null, a 91 topological charge 2 l    , an orbital angular momentum (OAM) of Nlħ, where N is the number of photons , and a spiral energy density flux about the propagation axis [92]. The phase singularity at beam center is called a screw dislocation [91]. Such modes can be generated by phase plates [93]. LG modes are solutions to the paraxial wave equation 2 0i z       for propagation along z, which is of the SE form with zero potential [92]. An example of an LG field at z=0 with a null on axis is     2 2/1 0, 1E (r ) rw rx iy w e      r , where ˆ ˆx y  r x y , and wr is the beam waist. Now consider, as an example, a construction similar mathematically to the LG field but physically quite different: the pulsed Gaussian field   2 2 22 E( , ) r ww r x i e e w w                r r , where gv t z   is the local position in a frame moving with the pulse group velocity gv , and w is the axial pulse length. This field has a moving null at  0, 0x   along the y-axis, perpendicular to the direction of propagation, and the circulation  on closed spatiotemporal contours in  ,x  around the null is quantized. We refer to optical vortices with phase winding along spatial and temporal axes as spatiotemporal optical vortices (STOVs). We contrast the STOV construction with other recent work, such as by Eilenberger et al. [94] or Mihalache et al. [95], where beams with non-trivial spatiotemporal coupling are created with embedded vorticity, but the vorticity is fundamentally similar to the LG construction given above, where phase and energy circulation from the vortex occur only in the transverse spatial dimensions. We refer to optical vortices 92 with phase circulation solely in spatial dimensions as spatial optical vortices, and note that, as far as we know, STOVs have only been considered theoretically, with the first mention by Yangirova and Sukhorukhov in the context of solitons [96–98]. Outside of optics there has been theoretical and experimental work on systems with analogs of the STOV construction, such as work by Nicholls et al. in acoustics [99], and Kartashov et al. in quantum fluids [100]. In this chapter, we demonstrate that STOVs are a fundamental and universal feature of optical pulse collapse and arrest in self-focusing media. Their existence in nonlinear ultrafast pulse propagation appears to be ubiquitous, and their creation, motion, and destruction is intrinsically linked to the complex spatiotemporal evolution of the pulse. 4.3 STOV generation in a collapsing pulse 4.3.1 Modeling STOV generation Optical collapse is a fundamental phenomenon in nonlinear optics [30] (see Section 1.4.2 of the introduction). It occurs when the laser pulse-induced change in the medium’s refractive index generates a self-lens whose focusing strength increases with intensity. Above a critical power level (Pcr), self-lensing exceeds diffraction, and the pulse experiences runaway self-focusing. In the absence of “arrest” mechanisms terminating self-focusing, the pulse would collapse to a singularity. In reality, additional physical effects intervene. For example, in the case of femtosecond filamentation in ionizing media [35], plasma generation acts to defocus the pulse when the peak self-focused pulse intensity reaches the ionization threshold. Other collapse arrest mechanisms include dispersion-induced pulse lengthening [32], 93 cascaded third order nonlinearities [101], vectorial effects from pulse non- paraxiality [34], and, in the case of relativistic self-focusing, electron cavitation [33]. In both calculations and simulations [14,102–104], it has been found that the following modified paraxial equation for wave evolution (see Section 1.3.4), in the form of a time-dependent nonlinear Schrödinger equation (NLSE), is well-suited to describe optical pulse collapse and collapse arrest for pulse propagation along z:   2 2 2 2 2 2 V 0ik k z                     (4.1) Here,  cE ( , , ) i k z t z e      r is the dimensionless electric field component,  is the pulse envelope, c 2 2 2 c 2 k c k             is the dimensionless group velocity dispersion (GVD), r and  are as before, and the axial position z can be viewed as a time-like coordinate. The physics of self-focusing and arrest is contained in the functional  V  . To demonstrate toroidal STOV generation in the arrest of self-focusing collapse, we perform propagation simulations using Eq. (1), imposing azimuthal symmetry, and include electronic, rotational, and ionization nonlinearities in  V  (see Section 1.3.4 or ref. [102,105]). The Gaussian input pulse is 3 mJ, 45 fs (Gaussian FWHM) with an input waist w0=1 mm, and the propagation medium is atmospheric pressure air. As our experiments (see Section 4.2.3) are intentionally operated in the single filamentation regime, where multi-filamentary modulational instabilities are precluded, azimuthal symmetry is a good approximation [35]. Figure 4.1(a) is a post-collapse plot of the pulse phase at z=120 cm, which shows the emergence of two oppositely wound and oppositely propagating STOVs (v1 (+1 94 charge, forward) and v2 (-1 charge, backward)) entrained between the higher intensity core of the beam and the beam periphery. The delayed rotational nonlinearity from the N2 and O2 air constituents [16] forms an additional vortex pair ~100 fs behind the main pulse, v3 (+1, forward) and v4 (-1, backward), where v3 is shown bisected in Fig. 1(a) and v4 has exited the simulation window. Further evolution of the simulation shows that v2 and v3 collide and annihilate, while v1 continues propagating with the most intense part of the pulse (to be discussed later). We note that the delayed generation of v3 and v4, where the pulse intensity is many orders of magnitude weaker than at its peak, shows that STOVs can also be generated linearly by an imposed spatiotemporal index transient. STOVs are not merely mathematical free-riders on intense propagating pulses: as we will see, a real energy flux j circulates either as a saddle (for 2>0) or as a circle (for 2<0) in the  ,r plane surrounding the STOV core. 95 Figure 4.1. (a) Phase and intensity projections of simulated pulse, in local coordinates, showing STOV generation. Propagation medium: atmospheric pressure air. Black dashed lines encircle vortices v1, v2 and v3, and arrows point in the direction of increasing phase. Our phase winding convention considers a (ξ ξ0) +i(rr0) winding about a null at (ξ0, r0) to denote a +1 STOV, giving the v1 STOV a +1 charge. The white dots on the intensity projection are centered on the locations of vortices v1, v2 and v3. (b) Simulation of an air filament crossing the air-helium boundary for the conditions of (a), showing the pulse fluence and plasma density. Nonlinear propagation terminates as the pulse transitions from air to helium, whereupon the pulse and a reference pulse (not shown, see Appendix A) are directed to an interferometer. To understand the generation of these spatiotemporal toroidal structures, recall that phase vortices in fields are closely associated with localized field nulls [91]. In (b) pulse to interferometer nonlinear propagation terminates 500 μm (a) v1v2v3  r rad TW/cm2 50 fs 96 arrested self-focusing, field nulls occur as a natural part of the dynamics and spawn toroidal vortices of opposite charge. This can be illustrated by the toy model of Figure 4.2, which shows the effect on a beam of abruptly spatially varying self-induced change in refractive index, which we model here as a sharp transverse step. We consider left-to-right propagation of the ‘half plane-wave’ pulse       2 g g0 NL/v / ( , , ) 0E( , , ) E δ h( ) e e e i z k v i x z x z x          , neglecting dispersion and diffraction, where  and z are as before, <<1, h(x)=1 for x<0 and 0 otherwise, 2 NL 2( , , ) E( , , )x z kn x z z   is the sharply stepped nonlinear Kerr phase, and n2 is the nonlinear index of refraction. At z=0, the phase fronts are aligned for all x. As the pulse propagates, the Gaussian intensity distribution causes the front of the pulse to red shift and the back of the pulse to blue shift, with spatiotemporal phase front shear developing between the pulse ‘core’ at x<0 and periphery at x>0. When the propagation reaches 2 1 v 2 0( )z z kn E   , the peak of the pulse in the core is  out of phase with the periphery, and the electric field magnitude nulls out at the single point v( 0, 0, )x z z    , marked by a circle, forming an ‘edge dislocation’ [91]. Upon further propagation, continued phase shearing spawns two null points of opposing (1) phase winding, marked by circles (z=2zv). These two STOVs, whose axes are along y (perpendicular to page), immediately begin moving apart; one vortex advances in time towards the front of the pulse while the other vortex moves to the back. Similar dynamics are theorized to exist in monochromatic breather-solitons, where ring shaped vortices are formed quasi-periodically throughout propagation [106]. 97 Figure 4.2. Toy model showing birth of a vortex pair via spatiotemporal phase front shear. The white curve and arrow depicts the axial (temporal) intensity profile and propagation direction, while the “core” and “periphery” labels denote the spatial intensity step. The z=0 panel shows the initial condition where phases are aligned, the z=zv panel shows the birth of the null (vortices overlap) and the z=2zv panel shows continued shear carrying the vortices apart, with the +1 vortex moving to the temporal front, and the -1 vortex moving backward. Our phase winding convention considers a (ξ ξ0) +i(xx0) winding about a null at (ξ0, x0) to denote a +1 STOV. In the third panel, red and black arrows indicate the direction of increasing phase. The two arrows connect the same lines of constant phase; the arrows show that the phase difference between head and tail is ill-defined due to the vorticity of the -1 STOV. These general features occur in simulations of self-focusing collapse arrest. In Figure 4.3, we show 2D profiles of the pulse’s phase v( , , )r z   and intensity 2 0 v 1 ( , , ) 2 nc E r z   for a simulation using Equation (4.1), tracking the moving plane =v where the phase singularity and field null first appear. Log lineouts of the intensity (normalized to 13 210 W cm ) are overlaid on the 2D profiles. As the initial (z=0) Gaussian input beam self-focuses, a strongly peaked high intensity central core (similar to the Townes profile [107]) develops, surrounded by a lower intensity periphery, with a sharp transition knee between them (z=156 cm). The associated 98 phase plot shows the central core having accumulated a much larger nonlinear phase shift than the periphery. During collapse, the knee moves radially inward, preventing phase shear from accumulating substantially at any particular radial location. However, as the core peak intensity rises to the point where ionization begins (z =160 cm, not shown), the location of the intensity knee stabilizes and highly localized phase shear builds up, greatly steepening the transition between core and periphery, with the field beginning to dip towards a null (z=166 cm). Finally, only a short propagation distance later (z=167 cm), sufficient shear develops that the core and periphery are  out of phase at the slice =v, with the dip in the simulation becoming orders of magnitude deeper, forming a ring-shaped null surrounding a core of relatively flat intensity and phase. Note the core-periphery phase difference cp jumps by 2 between z=166 and 167 cm. (As seen in the computation of Figure 4.1(a), the null then spawns two oppositely charged (1) toroidal vortices that propagate forward and backward.) The plane =v is still shown at z=178 cm, but the vortices have now migrated out of that plane. Later in the paper, we derive and discuss equations of motion describing vortex dynamics within the moving reference frame. 99 Figure 4.3. Simulations of beam phase (top) and intensity (bottom) at the axial/temporal slice v where the STOV pair first appears. From left to right, the pulse is advancing along z, with the input shown at z=0, a collapsing beam at z=156 cm, ionization onset at 160 cm (not shown), just before the vortices spawn at 166 cm, just after the vortices spawn at 167 cm, and an image where the vortex pair has moved out of the =v plane. The linear intensity images are overlaid with centered lineouts of 10 0log ( )I I , where 13 2 0 10 W cmI   . Experimental parameters were used as code inputs: w0 =1.3 mm, pulse energy=2.8 mJ (P/ Pcr=6.4), pulse FWHM 45fs. 4.3.2 Experiment 4.3.2(a) Experimental overview In order to experimentally confirm the existence of STOVs, we image the spatiospectral phase and intensity profiles of femtosecond laser pulses mid-flight during their pre- and post-collapse evolution in air. Until now, we have been discussing STOVs as a spatiotemporal phenomenon, but they are also vortices in their spatiospectral representation, which has enabled us to unambiguously observe them. Why should a pulse with a STOV have a vortex in the spatiospectral domain? For a small temporal chirp of the electric field where the vortex is embedded, there is a simple linear mapping between time and frequency. This leads to the vortex appearing in the spatiospectral as well as the spatiotemporal domain. Section 4.2.3(c) 100 discusses a Gaussian pulse with a temporally centered toroidal STOV and shows that the spatiospectral representation of the pulse possesses a vortex with the same spatial radius. We note that the relationship between vorticity in the spatiotemporal and spatiospectral domains is, in general, complex, with the full field distribution (including vortices) in one domain contributing to an individual vortex in the conjugate domain. Because our experimental images identify STOVs after optical collapse arrest has occurred, our results avoid this complexity, as simulations show that by this point there is only one STOV propagating with the filamenting pulse. Direct measurement of the spatial phase and intensity profile of a filament in mid- flight is not amenable to standard techniques; the use of relay imaging or beam splitters is subject to the severe distorting effects of nonlinear propagation and material damage. However, by interrupting nonlinear pulse propagation by an air- helium interface, the in-flight beam intensity and phase profile can be linearly imaged through helium, taking advantage of the very large difference in instantaneous nonlinear response between helium and air ( 2,He 2,air 0.04 n n  [16,50]). This helium cell technique was first employed by Ting et al. [43] to measure the in-flight intensity profile of a femtosecond filament. Here we extend the technique to also enable measurement of the pulse transverse phase profile. Section 4.2.3(b) discusses the details of how we produce the air-helium interface and how the <4 mm thick air- helium transition layer is sufficiently thin to enable distortion-free imaging of the air filament cross section. The interface is movable along the laser propagation axis to allow mid-flight filament intensity and phase imaging over the full propagation path. 101 A conceptual view of the experiment is presented in Figure 4.1(b), where we show the post-collapse pulse from the simulation of Figure 4.1(a) encountering the air-helium interface and terminating nonlinear propagation, after which it is relay imaged through a folded wavefront interferometer and combined there with a weak femtosecond reference beam with flat spatial and spectral phase. The resulting interferogram encodes the 2D spatial phase and intensity of the in-flight filamenting (signal) pulse at the air-helium interface, averaged over the ~20 nm bandwidth of the λ=800 nm reference arm. In effect, the reference arm spectrally gates the STOV, resulting in a spatiospectral interferogram centered at the reference pulse bandwidth. Spectral gating is crucial to our measurements. If the reference arm had the same wide bandwidth as nonlinearly-generated in the signal (filamenting) arm, the signature of the STOV, which is present over a smaller spectral window, would be washed out. The remaining sections of 4.3.2 provide detailed descriptions and justifications for the experimental method: a discussion of the experimental setup is given in Section 4.3.2(b), the spatiospectral representation of STOVs is addressed in Section 4.3.2(c), and the interferometric analysis is explained in Section 4.3.2(d). 4.3.2(b) Experimental setup The experimental apparatus makes possible the reconstruction of the transverse spatial phase and intensity profiles of a femtosecond optical air filament in mid-flight. To do this, we use an abrupt air-helium transition to halt nonlinear propagation, as ionization yield and self-focusing are both negligible in helium. The beam is then relayed linearly from the transition zone through the helium and 102 interfered with a reference pulse in an interferometer, enabling extraction of the transverse amplitude and phase profiles of the filament. Figure 4.4. Experimental setup for measuring the in-flight intensity and spatial phase profiles of collapsing and filamenting femtosecond pulses over a ~2m range. The experimental setup is shown in Figure 4.4. Our filamentation source is a chirped pulse amplification Ti:Sapphire amplifier ( λ=800 nm, 45 fs, 0-5 mJ). The beam from the laser is spatially filtered using a pinhole to produce a Gaussian mode with flat phase fronts – this is important for the reference arm in the experiment, which requires a flat spatial phase. After spatial filtering, the pump/filamenting (signal) arm and reference arm are generated using uncoated wedges in a Mach- Zehnder (MZ) configuration to create a large difference in power between the two pulses (~10 4 : 1). Here, the low power reference arm reflects off the front faces of the wedges, while the high power signal arm transmits, creating a dispersion imbalance that is corrected further downstream. The pulses are then compressed, with the signal arm now at 45 fs full width at half maximum (FWHM) intensity, rotated 90 degrees 103 using a periscope (converting P-polarization to S), and down collimated to a waist of 1.3 mm using a reflective off-axis telescope. After down-collimation, the pulses are launched a variable distance spanning 50-225 cm beyond the last optic of the telescope before nonlinear propagation terminates inside the nozzle of the translatable helium cell. Past the air-helium transition, both pulses propagate linearly in helium 50 cm, with the intense filament core of the high power beam expanding transversely in size. Both pulses are then wedge-attenuated before being directed out of the cell through a 200 µm thick BK7 window. Outside the cell, the high power pulse is attenuated to match the power of the reference arm via reflections from a second set of wedges in MZ configuration. In order to maintain polarization purity, polarization rotation from the upstream periscope was necessary, as wedge reflections preferentially select for s-polarization. Wedge transmissions by the reference arm fix the dispersion mismatch created in the pre-compressor MZ interferometer. The pulses are recombined at the output of the interferometer and sent through a lens which images an upstream plane, just before the nozzle’s gas transition region, to a CCD camera. The air-helium interface is formed by the non-turbulent flow of helium, at slightly positive pressure, into the ambient air through a 1/4” diameter nozzle on a translatable rail-mounted cell. The filament propagates from air into the nozzle and nonlinear propagation terminates over the sharp 4 mm transition from air to helium. The helium-air transition was measured, as in Ting et al. [43], by monitoring the strength of the 3 3 D – 2 3P , =587.6 nm helium line as the helium cell nozzle is moved through a tightly focused 800 nm, 45 fs ionizing pulse. As seen in Figure 104 4.5(a), the rapid drop-off of the helium line indicates that there is negligible helium beyond a 4 mm 10%-to-90% transition layer at the nozzle. To confirm the fidelity of imaging and phase reconstruction through the helium cell, we used a time domain propagation code [102] in 2+1 dimensions (r, z, ξ) to model the propagation of a filamenting pulse through the 4 mm air-helium transition into the far field in the bulk helium at atmospheric pressure. The accuracy of phase and intensity reconstruction was verified by reverse-propagating the beam via phase conjugation through vacuum back to the air region just before the transition. The results are displayed in Figure 4.5(b) and (c), which show that the reconstructed spatial intensity and phase at 800 nm (red) closely track that of the input electric field just before the transition region (black). In addition, we verified that small deviations from the correct imaging plane (+/- 1 cm) did not affect the results. 105 Figure 4.5. (a) Line emission (in arbitrary units) of the neutral helium 3 3 D – 2 3P , =587.6 nm transition induced by a tightly focused 800 nm pulse as a function of helium cell position. (b) Simulated filament intensity and (c) spatial phase at 800 nm (black) just before the ~4 mm air-helium transition, as well as reconstructed intensity and phase (red). The reconstruction is performed by propagating the solution 4 mm through the transition region followed by 50 cm of helium. This simulated far field is then back-propagated through 50.4 cm of vacuum to simulate reconstruction from imaging optics. 4.3.2(c) STOVs in spatiospectral space Although the STOV is a vortex in the spatiotemporal domain, our experiment measures the spatial phase of the filamenting pulse by interference with a reference pulse that is centered at λ=800 nm (see Section 4.2.3(d) below). What is the signature of a STOV in spatiospectral space? -5 0 5 10 15 0 0.5 1 position in mm lin e s tr e n g th ( a .u .) a) 0 0.5 1 0.2 0.4 0.6 0.8 1 transverse space in mm in te n s it y ( a .u .) b) original reconstructed 0 0.5 1 -4 -3 -2 -1 0 p h a s e i n r a dc) 106 Consider an ‘R’ vortex ring centered at 0 0( , )r  embedded in a Gaussian background field with a temporal chirp, 2 2 2 2 2w 0 0 0( , ) [ i( )] rr w iaE r E e e r r           (4.2) where rw and w are transverse and longitudinal widths, respectively, and a is a chirp parameter. The Fourier transform along the  axis is, up to a complex coefficient,    22 2 2 2g 2 /w / ( ) 2 2 0 0( , ) (1 ) ( )( ) 2 r cr v ia w c g iw E r e iaw r r aw i v                              (4.3) Where gk v  , the variable conjugate to  is k , and c is the central frequency of the pulse. In spatiospectral space the vortex is centered at 0 0 0 0 0 g2 4 0 1 1 (1 ), 2 v cr r a r aw aw                 (4.4) It is seen that the STOV is manifested at a definite spectral location depending on its axial position 0 and pulse chirp a . For a STOV perfectly centered in the pulse ( 0 0  ), its spatiospectral representation is located at the central frequency ( c  ) and at the same spatial radius ( 0 0r r  ). For an air filament, the STOV is, indeed, positioned temporally about the peak intensity of the filamenting pulse, as this is where the Kerr effect-driven inward flow of optical energy yields to the plasma refraction-induced outflow, as will be discussed in Section 4.2.5. 107 Figure 4.6. Spectral phase (radians, top) and spectral intensity (arbitrary units, bottom) taken from simulation at z=180 cm as the vortex core in the spatiospectral domain crosses 2.35 PHz (λ=800 nm). Simulation parameters are the same as Fig. 3 of the main text: w0 = 1.3 mm, pulse energy = 2.8 mJ (P/ Pcr = 6.4), pulse FWHM 45 fs. From left to right, the top row (bottom row) shows the full spatial-spectral phase (intensity), as well as spatial phase (intensity) slices spectrally fore and aft of 800 nm. The dashed black line in the top left panel intersects the vortex core. Figure 4.6 displays the spatiospectral representation of a STOV using simulation output. We use the same simulation parameters presented in Figure 4.3: w0 = 1.3 mm, pulse energy = 2.8 mJ (P/Pcr = 6.4), pulse FWHM 45 fs. The top row shows spatiospectral phase, while the bottom row displays the spatiospectral intensity. The leftmost column shows the complete spatiospectral phase and intensity, and makes clear that there is phase vorticity in the spatiospectral domain. The middle (rightmost) column shows the phase and intensity spectrally fore (aft) of the vortex where we present the data in a manner similar to the experimental images of Figure 108 4.9 shown in Section 4.2.4. Here, ‘fore’ and ‘aft’ refer to left and right of the dashed black line in the upper left panel. The diameter of the core is well reproduced, as are the abrupt phase jumps of ~π from core to periphery on both sides of the vortex. The spectral intensity images show the dark ring of the vortex core. The location of the vortex core changes as the pulse propagates. The phase flip seen in Figure 4.6 by sampling spectrally to either side of the core would also be seen by observing a fixed frequency while scanning the input power or by scanning along the propagation axis of the filament, which is how we observe the phase flip experimentally in Figures 4.9 and 4.10. 4.3.2(d) Interferometric reconstruction Since the collimated low power reference arm has a flat spatial phase, the spatial phase difference, extracted from the interferogram using standard techniques [51,108], is just the spatial phase accumulated by the nonlinearly propagating pulse. Such a pulse can develop a complicated time (and therefore frequency) dependence [35,109,110]. However, our interferograms are spectrally gated by the narrower reference pulse spectrum. What is actually measured, therefore, is a weighted average of the spatial phase as a function of frequency, where the weighting factor is the product of the (narrower) spectral amplitude of the low power reference arm and the (broader) spectral amplitude of the filamenting arm. The oscillatory portion of the signal on the CCD is given by sin * ref sigint( , ) 2Re ( , , ) ( , , ) ikxx y e d A x y A x y               (4.5) 109 where x and y are transverse coordinates in the pulse, refA and sigA are Fourier transforms of the field envelopes of the reference and signal pulses, k is the central wavenumber of the reference pulse, and θ is the crossing angle of the two beams. As the spatially filtered reference pulse propagates linearly, refA is fully known, and is well approximated by a Gaussian with flat spectral phase,  2 2 2 2 20 ref 0( , , ) x y w A x y A e        and the weighted spatiospectral phase of the signal pulse is then extracted. Figure 4.7. Comparison of the spatiospectral phase of the pulse at 800 nm, to the weighted spatial phase measured about 800 nm using the same simulation parameters considered in Fig. 3: w0=1.3 mm, pulse energy=2.8 mJ (P/ Pcr=6.4), pulse FWHM 45 fs. Flipping of the spatial phase occurs at z=179 cm at 800 nm and z=185 cm for the weighted spatial phase. 110 What is the difference between the weighted spatial phase which can be extracted from Equation 4.5 and the spatial phase at a single frequency (such as at the spectral center =800 nm)? Figure 4.7 uses simulation output to directly compare the spatiospectral phase at =800 nm and the weighted spatiospectral phase. We simulate using the same input parameters as Figure 4.3: w0 = 1.3 mm, pulse energy = 2.8 mJ (P/ Pcr = 6.4), pulse FWHM 45 fs. The Figure tracks the evolution of the two different spatial phase quantities as the STOV in spatiospectral space crosses =800 nm. It is apparent from Figure 4.7 that the two quantities follow each other closely, and critically, they both exhibit the flip in phase we use to identify the vortex. The two main differences in the phase quantities are that the 800 nm spatial phase flips first at z=179 cm, followed by the weighted phase at z=185 cm (likely due to asymmetric dispersion about 800 nm), and that the weighted phase gives a core region ~20% larger than the exact spatial phase at 800 nm. Figure 4.7 establishes that the weighted spatial phase is a good proxy for the spatial phase at =800 nm, and can be used to detect the STOV in spatiospectral space as outlined in section 4.2.3(c). 4.3.3 Experimental results The experiment was performed by scanning the helium cell axial position over a range covering both pre- and post-collapse propagation for all energies used. The data consists of a densely spaced collection of beam intensity and phase images at various cell positions ( hz ) and input pulse energies ( i ). For all measurements, collimated beams were launched with w0=1.3 mm and FWHM intensity pulse width =45 fs. 111 Figure 4.8 shows the beam on-axis phase shift  with respect to the interferometric reference pulse as a function of P/Pcr at a fixed position of zh =150 cm after launch, where 2 cr 0 2 3.77 8 P n n    for our Gaussian input beam profile and iP    is the input power. The red points are averages over 2600 shots (blue points) in 150 energy bins. It is important to note that the scatter in  of roughly 1rad at any given laser power is constant across all powers measured, including cr 1 P P , where we could not detect any nonlinear phase. Therefore, the scatter is due to the shot-to-shot interferometric instability of the measurement and is not intrinsic to the filamentation process. The most striking aspect of the plot is the abrupt jump in beam central phase of approximately ~2 at cr 5 P P . The phase goes from positive and rapidly increasing (increasing self-focusing) to abruptly negative (defocusing), providing a clear signature of the transition from the pre-collapse to post-collapse beam. For nominally constant laser power right at the jump, the phase fluctuates in the range ~, showing the extreme sensitivity of the phase flip. 112 Figure 4.8. Beam on-axis phase shift (with respect to flat-phase reference arm) as a function of pulse power at zh = 150 cm. The phase jumps abruptly by ~2 at P/Pcr ~ 5, providing a clear signature of the transition from the pre-collapse to post-collapse beam. The red points are averages over 2600 shots (blue points) in 150 energy bins. A more revealing way to display what is happening at the collapse is shown in Figure 4.9. Here, for given zh, the phase images ( , )h iz  were searched for i or cr P P r where the central phase appeared to randomly flip sign from shot to shot. These are the powers at which pulse collapse was observed for each position, just as cr 5 P P was for zh=150 cm in Figure 4.7. The top row of Figure 4.9 shows beam phase and intensity images for input power cr 4.4 P P  (at zh= 165 cm). Because the onset of collapse arrest is extremely sensitive to fluctuations in the pulse energy (as seen in Figure 4.8), these images span the possibilities of pre-arrest through post- arrest of the 113 collapse, and typically three types of images appear. Panel (i) shows strongly-peaked intensity and positive phase; the beam is collapsing, but arrest has not yet begun. Panel (ii) shows radically different images, as does panel (iii): the intensity images show narrow ring-shaped nulls embedded in a relatively smooth background, and the phase images show a sharp yet smoothly transitioning jump close to π or -π across the rings, with the phase jumps flipped between (ii) and (iii). We note that the smooth phase transition from the periphery to the core rules out 2π phase ambiguities in interferometric phase extraction. The bottom row of Figure 4.9 plots, for a range of cr P P , the phase difference cp between the core and periphery of the beam. To do this, for each phase image we compute the difference between the maximum and minimum values of the phase within a 60 micron box centered about the largest spatial phase gradient, the radial location of which defines ‘core’ and ‘periphery’. For each nominal value of cr P P , it is clear that as the phase gradient becomes large, the phase difference saturates at π. Near saturation, roughly 50% of the shots have the core phase advanced from the periphery while the others show the reverse. 114 Figure 4.9. Top row: retrieved spatial phase (in radians) and intensity (arbitrary units) images at z = 165cm, P/Pcr = 4.4 for i) a pre-collapse beam, ii) and iii) beams where a vortex ring is on either side of the reference central wavelength of 800 nm. The bottom row shows that as the maximal phase gradient in the images increases, the maximal phase shift saturates at π for all cases of P/Pcr leading to beam collapse. The evidence from Figures 4.8 and 4.9, and comparison to the simulation of Figure 4.3, strongly suggest that we are imaging spatiotemporal vortex rings: the abrupt appearance of ring-shaped nulls in the field magnitude accompanied by 2π phase jumps in cp across the null– these are exactly the signatures of a vortex. Because the circulation around a general singly charged vortex is 2, examining our vortex in the spatiospectral domain ( , )r , one would expect, depending on the sign of vortex winding, that the core-periphery phase difference cp jumps by 2π from  slices just before the vortex ( cp    ) to  slices just after ( cp   ). For example, before a vortex of charge +1 the core is phase-advanced with respect to the periphery; after the vortex it is phase-retarded. This is exactly what we observe experimentally and what is predicted in the simulation of Fig. 3 and its spatiospectral 115 counterpart, where even the ~400m diameter of the vortex ring is accurately determined. Of the four STOVs simulations show are generated at collapse arrest in air (see discussion of Fig. 1(a)), only the temporally foremost +1 STOV does not annihilate or separate from the bulk of the pulse . Using Fig. 1(a) as a guide, we interpret our results as a spectral “fly-by” of a +1 STOV from the blue to the red side of our reference pulse spectrum centered at 800 nm. We note that a similar fly-by of a 1 STOV from red to blue would present itself in an identical manner. How are STOVs born in real collapsing pulses? In our ( , , )r z simulations, vortices are immediately born as tori around the beam owing to azimuthally symmetric ( -independent) phase shear. In real beam collapse, where there is  variation in the laser field, topological considerations lead us to expect that shear in higher E-field locations will first lead to a point null, followed by an expanding and distorted torus on one side of the beam that progressively wraps to the other side of the beam and then, meeting itself, splits into two toroidal STOVs of opposite phase winding. The onset of these STOVs, aligned with planes of constant  , has a beam- regularizing influence, as seen in the images of Figure 4.9, which show remarkably flat phase and intensity profiles inside the ring. This could be the reason for the notably high quality spatial modes and supercontinuum beams (so-called ‘spatial cleaning’ [111]) observed in filamentation. We are performing 3D propagation simulations to verify this scenario. We also note that the ring null forms a natural and well-defined boundary between what had been qualitatively labelled the ‘core’ and ‘reservoir’ regions [35] in femtosecond filaments. 116 4.3.4 STOV dynamics and energy flow Once STOVs are generated, it is important to understand how they propagate. Following the method of ref. [112] as applied to spatial vortices, we approximate the local form of the STOV as a spacetime “R-vortex” vortex 0 0( ) (r r )i      of charge 1 with a linear phase winding about 0 0( , )r , embedded in a background field envelope bg such that bg vortex   . If we take bg ie   , where  and  are the real amplitude and phase of the background field, then as the pulse propagates along z , the nonlinear Schrödinger equation in ( , , )r z moves the vortex location vortex 0 0( , )r r according to: 0 0 0 vortex 2 2 ( , , ) ( , , ) 1 1ˆ ˆ ˆˆˆ ˆ 2 r z r z k z r r r                                         r r ξ r ξ ξ (4.6) where the derivatives are evaluated at the present vortex core location 0 0 0( , , )r z (See Appendix A.5). Equation (4.6) demonstrates interesting analogies with fluid vortices. The term ˆ 2r  ξ propels the vortex forward or backward depending on its charge and radius (curvature), and strongly resembles the speed 4 r  of a toroidal fluid vortex (such as a smoke ring) [113]. Identifying 2 ˆ k r          j r as the local effective velocity associated with the background electromagnetic flux (see Appendix A.6), we interpret it as a charge-independent drag-like term, expanding (contracting) the STOV for power outflow (inflow) for 0( 0) r     . The term 1ˆˆ r       r  is a Magnus-like 117 motion [86], here propelling the STOV along ˆξ , perpendicular to the vortex circulation vector ˆ and the ring expansion/contraction direction rˆ . The terms 2 ˆ      ξ and 2 1ˆ ˆ       ξ  are their spatiotemporal analogues. In gases, small 2 (~10 -5 ) makes these terms negligible; they contribute much more significantly in solid media. We note that bg is not a fixed field independent of vortex motion. Equation (2) should be understood as a stepwise predictor of vortex motion based on an updated bg . Figure 4.10. Demonstration of energy flow about a +1 “r-vortex” STOV. White arrows correspond to size and direction of j in Equation (4.7). There are three distinct regimes: the left panel shows the saddle regime which exists in regularly dispersive media (2>0), the middle panel shows the degenerate regime which can exist in regular or anomalously dispersive media and has a dominant axis for energy flow, and the right panel shows the spiral regime which exists in anomalously dispersive media (2<0). 118 To understand energy flow near STOVs, it is useful to consider the electromagnetic energy flux associated with the full field envelope iue   in the moving frame of the pulse (see Appendix A.6) 2 2 1 ( )u k         j ξ (4.7) where one can see that the sign of 2 determines whether the energy flow near the STOV is spiral ( 2 0  ) or saddle ( 2 0  ). What are the relative contributions of longitudinal and transverse energy flow about a STOV? For filamentation in air, 5 2 10  , the characteristic axial pulse length and width are ~c×10 fs and ~100 µm (filament core), and one finds 410 1 r j j   . Here, the distinction between saddle and spiral is not important, as we are in the degenerate case. For solids, however, where 2 2 10  , we expect the distinction between saddle and spiral energy flow to be very important. Figure 6 provides intuitive visualization of the energy flow pattern for the saddle, degenerate and spiral cases. Equations (4.6) and (4.7) are useful for an intuitive picture of the governing dynamics of STOVs, especially when viewed together with propagation simulations. For example, for the four STOVs seen in the simulations of Figure 4.1(a), their dominant early movement is governed by 1ˆ r      ξ in Equation (4.6), which propels the +1 (-1) STOV temporally forward (backward), with the forward motion initially being superluminal. (We will explore the detailed implications of superluminal STOV motion in a future publication.) A consequence of the opposing directions for the  119 charges is collision and annihilation of v2 and v3 of Figure 4.1(a), as discussed earlier. Remarkably, the v3 STOV superluminally climbs from a region of negligible intensity through many orders of magnitude of increasing intensity to reach and annihilate v2, whereupon a local depression is left in the field that more gradually dissipates. The v1 STOV eventually settles in the temporal middle of the highest intensity portion of the pulse, propagating at nearly vg. Evidently,  self-consistently evolves to balance the ξˆ terms in Equation (4.6) and  flattens along the radial dimension (as indicated in the experiment), preventing expansion or contraction of the STOV. Our simulations show that a surviving +1 STOV is always coupled to the filamenting pulse. This is no coincidence, as the energy flow for a +1 STOV is toward (away from) the pulse axis temporally in front of (behind) the STOV. This is exactly as expected, where the front of the pulse draws energy in by Kerr self-focusing, and energy at the back of the pulse is directed outward by plasma refraction. A link and descriptions of STOV movies are at ref. [65]. There appears to be a deep connection between the STOV picture of filamentation—spontaneous generation of STOVs followed by STOV- governed energy flow in the pulse—and the spontaneous formation of conical nonlinear waves (X and O waves), which have been used to explain propagation dynamics of a filamenting pulse [114–116]. In real beams without  symmetry, we expect collisions of oppositely charged STOVs to be much more complex, although the beam regularization observed in experiments may conspire to promote collisions. Auxiliary 3D+1 linear propagation simulations in which we imposed STOVs as initial conditions on Gaussian pulses show repulsion of like-charged STOVs, which pass around each 120 other, and splitting of higher charge STOVs into multiple STOVs of single charge. We note that our measured air-based STOVs are not solitons, as diffraction does not balance self-focusing for a dark object. STOV solitons could exist, however, in an anomalously dispersive, self-defocusing medium [96,97]. 4.3.5 Conclusions A STOV is an optical vortex with phase circulation in a spatiotemporal plane. STOVs form naturally as a consequence of arrested self-focusing collapse and their dynamics influences subsequent pulse propagation. STOVs can also be imposed linearly via prescribed spatiotemporal or spatiospectral phase shifts, making possible their engineering for applications. While evidence for STOV generation was demonstrated in experiments and simulations of short pulse filamentation in air, we expect that STOVs, whose dynamics are subject to topological constraints, are a fundamental and ubiquitous element of nonlinear propagation of intense pulses. STOV-STOV interactions should prove to be a fundamental mediator of intra- and inter-beam dynamics. 4.4 Linear propagation and angular momentum of STOV- carrying pulses in material media Section 4.3 presents the first experimental detection of a STOV, which we show are self-generated in pulses undergoing optical collapse arrest and are intrinsic to the subsequent self-guided evolution of the pulse. While the presence of STOVs in such nonlinearly propagating pulses provides a compelling reason to study their fundamental properties, STOVs are not an inherently nonlinear phenomena. In this 121 section, we study the evolution of STOV-carrying pulses in the linear dispersive propagation regime, with particular emphasis on their angular momentum. We propose a simplified model to calculate the orbital angular momentum in a pulse with a STOV. The model provides physical intuition for the origins of the angular momentum. For STOV-carrying pulses in dispersive media, we identify a modal solution for the case of anomalous dispersion, and a “quasi” modal solution for normal dispersion. 4.4.1 Angular momentum in STOV-carrying pulses Here, we wish to calculate the angular momentum present in a linearly polarized STOV-carrying optical pulse. We consider the case where the pulse is propagating through a dispersive dielectric material, with spectral content far from any resonances, so that the material can be approximated as lossless. When speaking about the angular momentum of the pulse, we include contributions from both the electromagnetic field as well as the material response. The “correct” partitioning of momentum and angular momentum into electromagnetic and material subsystems is the topic of a century long investigation: the Abraham- Minkowski controversy [117–121], which we do not engage with as we are considering the total angular momentum, and do not discriminate between electromagnetic and material contributions. We consider Maxwell’s equations (in SI units) for the fields in a uniform, linear, lossless dielectric with zero free charge density and zero free current density: 0D  , 0B  , 0t  E B , 0t  H D . The material linear response 122 is expressed in the frequency domain as 0( , ) ( ) ( , )r    D r E r (tilde indicates Fourier transform) while 1( , ) ( , )H r B rt t , with 0  . If the time domain fields are taken to be real, then *( , ) ( , )E r E r   , *( , ) ( , )D r D r   , and *( ) ( ) ( )r r r         since the medium is lossless. The last relation allows expansion of the permittivity ( )r  as an even function of  about zero: 2 r r,2 0 ( ) nn n        (4.8) The conserved linear and angular momentum densities within a lossless, dispersionless dielectric are well known [12], and can be derived through direct manipulation of Maxwell’s equations. In the presence of dispersion, momentum and angular momentum densities were only recently derived by Philbin et al., from an action principle using Noether’s Theorem [122,123], 2 2 10 r,2 1 1 ( 1) ( ) ( ) 2 n n m n m m n t t n m                p D B E E (4.9a) 2 1 20 r,2 1 1 ( 1) ( ) ( ) 2 n n m m n m n t t n m               r p E El (4.9b) where the momentum density p is seen to be the Minkowski momentum D B plus an additional dispersive term, and the angular momentum density l takes on the familiar form r p plus an additional dispersive term. The additional dispersive terms are, in general, difficult to compute exactly, but are typically small compared to leading terms provided the pulse is far from any resonances, and the spectral bandwidth is small relative to the central frequency. 123 To calculate total angular momentum, we assume a form for an optical pulse with an embedded line STOV. We take the electric field at time t=0 to be  0 c( , 0) ( , , )exp( )t E z w im x w x w y w z w ik z    E r A (4.10) Here, x, y, and z are spatial variables, with the pulse travelling in the z-direction. The prefactor E0 is the electric field amplitude, m=±1 determines the direction of the phase winding, w|| and w⊥ are the longitudinal and transverse scale lengths of the pulse, kc is the central wavelength of the pulse, and A is the pulse envelope. We impose several restrictions on A so that the angular momentum in the pulse represents the intrinsic, or coordinate independent, orbital angular momentum [124]. We take the beam to be approximately linearly polarized, ˆA A x , so that there is no angular momentum contribution from the polarization state [124], and we assume A is symmetric about the x and z axes so that a lop-sided field distribution does not contribute to extrinsic orbital angular momentum, orbital angular momentum that depends on choice of coordinates [124]. Additionally, we assume that the spectral bandwidth of A is small compared to the central frequency c . 124 Figure 4.11. Simple model of the propagation effect of the envelope phase of a pulse with an embedded STOV. Color gradient indicates change in phase and central circular black arrow indicates the direction of phase winding. In regions a) and c) respectively, the phase winding blue and red shifts the pulse as indicated by the blue and red arrows (indicating momentum) and vertical lines (wavefronts). In regions b) and d) respectively, the phase winding deflects the local momentum density of the beam downward and upward as shown by the solid black arrows/wavefronts. Dashed black arrows indicate the momentum density and unperturbed wavefronts for the case of a constant envelope phase. To obtain physical insight and to straightforwardly derive expressions for angular momentum, we develop a simple four-quadrant model of the local spacetime environment of a STOV, as depected in Figure 4.11. The diagram shows the top and bottom portions of the vortex region respectively blue-shifted and red-shifted, while the temporal front and back of the pulse experience downward and upward tilt of the a c bd z x momentum density 125 phase fronts. Because the medium is dispersive, red and blue shifting sections of the pulse result in a differential in the momentum density from top to bottom, while the upward and downward deflections of the phase fronts cause the direction of the momentum density to deflect upward and downward as indicated on the diagram. Deviations from the momentum density of a plane wave of equal field amplitude travelling in the z-direction reveal the circular and saddle character of the momentum flux for the beam in anomalous and normally dispersive media. The simple model shown in Figure 4.11 suggests an intuitive method for estimating the angular momentum contained in the fields: approximate each quadrant as a section of a plane wave, with the top and bottom quadrants approximated by red and blue shifted plane waves travelling in the z-direction, and the left and right quadrants approximated by plane waves at the central frequency which are deflected upward and downward into the x-axis from the z-axis. Maxwell’s equations admit plane wave solutions of the form:   0 1 . 2 i t e c c     k r E E ,   0 1 . 2 i t e c c     k r D D , and   0 1 . 2 i t e c c     k r B B , where 0 0 0( )  D E and 0 0    k B E . For plane waves, the exact form of the angular momentum density (4.9b) simplifies to     c 0 r g* 20 0 0 0 ˆIm 4 2 r t t t nd E d                      l r D B E E r k , where ( )g k n c     is the group index, and t denotes a cycle average. For the top and bottom quadrants, labelled a) and c) in Figure 4.11, we take the plane wave as directed along the z-direction with frequencies shifted from the central 126 frequency of the pulse due to the winding phase of the vortical term. Upon integration the angular momentum is   r g2,ac 0 0 r g r g 1 2 8 y t nw wm u l V E m n n c n w                 (4.11) Here, ,acy t l is the cycle-averaged angular momentum in the y-direction (vortex axis) resulting from integrating regions a) and c). V is the volume of a quadrant, ( )r r c      are the blue and red shifted dielectric functions, g g c( )n n     , n is the refractive index at c , 2 0 02u V E is the vacuum energy of the pulse, and 4 c n w    is the frequency shift inherent to the phase winding. Calculating the angular momentum contribution from the left and right quadrants, sections b) and d) in Figure 4.11, we find a cycle-averaged angular momentum r g ,bd 8 y t n wm u l n w     (4.12) Where ,bdy t l is the cycle-averaged angular momentum in the y-direction from integrating regions b) and d) of Figure 4.11, and we have estimated the plane wave deflection into the x direction as c4 x z k k w k    , which is valid for paraxial pulses c 1k w . The total, cycle-averaged angular momentum is: r g r g ,ac ,bd 8 y y yt t t n n wwm u l l l n w w                (4.13) 127 The overall angular momentum is proportional to the vacuum pulse energy u , with an orientation which depends on the sign (charge) of the vortex (such as 1m   ), and a numerical factor of 8  specific to our simple model and resulting from our choice of how the phase winds about the pulse z x im w w  . The ratios involving the pulse transverse size and longitudinal length w w  and w w result from the spatial dependence of the angular momentum density l r . The contribution to the angular momentum from regions a) and c) of Figure 4.11 (first term of Equation 4.13) can be split into two parts, proportional to r     and gn    . The first term, r     , is attributable to a differential in the magnitude of the material response between regions a) and c), while the second term, proportional to the group velocity dispersion gn c GVD      , occurs due to a differential in the group velocity between regions a) and c). Contributions to the angular momentum from regions b) and d), attributable to wavefront tilt, complement results already derived by Bliokh and Nori [98]. The result is also similar to calculations done by Allen et al. [92] who considered orbital angular momentum in vortex carrying monochromatic Laguerre-Gauss modes. Although Allen et al. consider spatial optical vortices, the reason the beam contains orbital angular momentum remains fundamentally the same: wavefront tilt due to the presence of a vortex. Allen et al. showed that LG modes with azimuthal index l have 128 angular momentum equal to Nl , where N is the number of photons. To compare to Allen et al., we consider propagation in vacuum, and impose w w , which simplifies Equation (4.12) to ,bd 8 y t l Nm   , where we express results in terms of the number of photons u N   . The results are identical up to the numerical factor of 8  . A factor of 2 difference is expected because our result only considers half of the pulse (regions b) and d)) while the Allen result considers the entire beam – the remainder of the difference in their ratio 0.8 is attributable to differences in the distribution of the electric field. While the plane wave model provides an intuitive explanation for the origin of angular momentum in a STOV carrying pulse, and clearly delineates between contributions due to wavefront tilt and contributions due to dispersive effects, a more rigorous derivation is required to confirm the validity of the model. We are presently in the process of performing this derivation, which involves deriving an approximate form to the momentum and angular momentum densities of Equations (4.9) suitable for slowly varying envelopes, which can then be applied to a fully consistent solution of Maxwell’s equations containing a STOV. 4.4.2 Free space STOV mode propagation How do pulses with STOVs evolve as they propagate? To address this we consider the scalar wave equation under paraxial and slowly varying envelope approximations, subject to 2 nd order dispersion: 129  2 222 zik u u     (4.14) Here x and y are transverse dimensions, gv t z   is the longitudinal moving frame coordinate, gv is the group velocity, z is the longitudinal direction, k is the central wavenumber, 2 2 2 x y    is the transverse Laplacian, u is the electric field envelope which specifies the scalar electric field: ( , , , ) ( , , , ) ikE x y z u x y z e    , and 2 2 2 2 c c k c k             is the dimensionless group velocity dispersion. Separation of variables leads to the familiar Hermite-Gauss solutions to the Helmholtz equation: ,a y,b ,cxu u u u (4.15a) 22 2 ( ) 2 ( )( )0 ,a 2 ( ) ( ) ( ) a xx kxx i z R zw zx x a x x w x u x H e e w z w z                  (4.15b) Here, 0xw is the beam waist, 2 0( ) 1x x Rx z w z w z         is the z-dependent spot size, 2 0 2 x Rx kw z  is the Rayleigh range, aH is the a th order Hermite polynomial, 2 ( ) 1 Rxx z R z z z           is the curvature, and 1( ) ( 1) tana Rx z z a z           is the Guoy phase. Solutions for ,y bu are identical in form to ,x au , while solutions for ,cu are obtained with the substitutions: 2 2 , sgn( ) k k z z     . Scaling the effective wavenumber to 2 k  accounts for the different rate of spreading present in dispersion relative to diffraction, while flipping the sign of z for normal dispersion results in 130 phase accumulation identical to a left-going wave  i kz t e    , opposite the right- going wave convention  i kz t e   used for the transverse directions in Equation (4.14). A pulse containing a STOV oriented along the y-axis, similar in form to the pulse considered in Equation (4.10) can be constructed out of the modal solutions (4.15). Here we wish to construct a pulse that exhibits modal evolution, and so impose Rx Ry R Rz z z z   .  y,0 ,1 ,0 ,0 ,1x xu u u u imu u   ,a y,b ,cxu u u u (4.16a)    2 2 2 2 2 23 2 22 2 2 0 0 2 0 2 ( )2 sgn( ) ( ) ( ) ( ) sgn( ) ( )0 2 2 ( ) ( ) x x ik x y x y R zi z w z i z i zx x x w im u e e e xe e w z w z                                (4.16b) As before, 1m   determines the winding of the vortex, and we have expressed the pulse using pulse parameters from the xu family of solutions provided in (4.15b). Adding ,1 ,0 ,0 ,1x xu u imu u  is analogous to using 0 th and 1 st order TEM modes to produce a first order Laguerre-Gaussian pulse, but with the mode now lying in the  ,x  plane instead of the transverse spatial plane. Note that phase accumulation in the pulse envelope in Equations (4.16b) is equal for both terms when 2 0  , resulting in an amplitude which simply dilates with propagation (a modal solution), but when 2 0  the Guoy phase present in the vortical term (bracketed) exhibit right-going and left-going signatures respectively, resulting in what we call a quasi-mode. 131 Figure 4.12. Evolution of pulse intensity of a quasi-mode (top) in normally dispersive media, and mode (bottom) in anomalously dispersive media. The pulse envelopes have been rescaled so that they have equal size at all z locations. Pulse intensity is shown in the y=0 plane. Phase winding due to vorticity is indicated with red arrows. How do the solutions (4.16b) behave as they propagate? Figure 4.12 displays plots of the intensity as a function of propagation along the z-axis in the y=0 plane. Intensity plots are shown for five different locations: 0, ,Rz z   . For ease of viewing we have normalized the spatial scales so that all pulse envelopes have the same size. We have also chosen to plot solutions which are identical at the waist, z=0. For anomalous dispersion, all separable solutions evolve according to the “right- propagating-wave” convention, resulting in a free-space mode where the pulse amplitude simply dilates with propagation, in a manner analogous to a transverse LG01 beam, but with the pulse lengths scaled to reflect different magnitudes of diffraction and dispersion. In the case of normal dispersion, the transverse dimensions still obey the right-going convention, but the moving frame coordinate  now obeys the left-going convention, and the result is the quasi-mode evolution shown in the first row of Figure 4.12. The sign of the phase winding, indicated with a red arrow, flips while propagating from z   to the waist at 0z  , and flips once again while propagating from the waist to the positive far field at z  , with intermediate regions at Rz z  where the pulse is transitioning between opposing windings. 132 Interestingly, the “quasi” modal evolution of the STOV pulse in a normally dispersive medium changes the direction of its phase winding while preserving total orbital angular momentum. We are presently in the process of examining this situation, and deriving equations for angular momentum conservation for scalar pulses in dispersive media. 4.4.3 Conclusion We have studied the angular momentum and pulse evolution for linearly propagating pulses with STOVs. We identified two distinct sources for net orbital angular momentum due to the presence of a STOV: wavefront tilt and dispersion from the polarization response of the medium. Contributions to STOV beam angular momentum from wavefront tilt are fundamentally similar to wavefront tilt contributions to the orbital angular momentum of spatial optical vortices [92], while the dispersive contributions are physically distinct and specific to STOVs. We have also considered the linear propagation of STOV-carrying scalar optical pulses, identifying a free space mode in anomalously dispersive media, and a quasi-mode for normally dispersive media. We are presently in the process of rigorously deriving expressions for the total orbital angular momentum in STOV-carrying pulses and related equations for angular momentum conservation in dispersive media. 133 Chapter 5: Summary and future directions 5.1 Summary The research in this dissertation explores two topics: i) how to use the long timescale, neutral gas hydrodynamic response following filamentation to guide and deflect optical beams, and ii), the discovery of spatiotemporal optical vortices (STOVs), for which we provide the first experimental detection, and establish that STOV generation is a fundamental part of the optical collapse process, occurring in all filamenting beams. Chapters 2 and 3 focus on the long timescale, neutral gas hydrodynamic response following a filamenting pulse. In Chapter 2 we demonstrate that a pulse train of filamenting beams undergoes a sustained deflection in its trajectory due to the cumulative density hole generated from the pulse train. This effect is present in all high repetition rate (>1 kHz) filamenting pulse trains in gaseous media. In Chapter 3 we show how the neutral gas hydrodynamic response may be used for the formation of acoustic and thermal waveguides. These guides may be useful for ranged atmospheric applications such as remote sensing and the guiding of high average power beams. In Chapter 4 we introduce STOVs, and provide the first experimental detection, supported by theory and simulation, of a STOV. We develop a simple model of optical collapse arrest generic to all short pulse optical filamentation and self-guiding scenarios. We argue that STOV formation is universal, applying to all beams undergoing optical collapse in any system, independent of the details of the collapse arrest mechanism. Propagation simulations of air filaments show that in the 134 self-guided phase of propagation one of the STOVs has settled down to form a torus surrounding the filament core and co-propagates with it, mediating energy flow to and from the reservoir. Finally, through studying the linear propagation of STOV beams, we identify free-space mode solutions and calculate the orbital angular momentum of STOV-carrying beams. The remainder of this chapter suggests extensions to existing work. 5.2 Acoustic manipulation of the filament-induced air density hole The density depressions produced in the wake of a filament were shown in Chapter 2 to lead to a persistent deflection in the pulse train of a filamenting beam. By acoustically manipulating the position of the density hole, a filamenting pulse train can be dynamically steered. Figure 5.1, below, shows preliminary data where we demonstrate deflections in the beam pointing which result from passing a beam above an acoustic source (in this case a car subwoofer oscillating at 25 Hz). With sound waves driven with speaker currents of 4 A RMS, we stimulate deflections in a 1 kHz pulse train which vary sinusoidally in time, tracking the waveform of the driver. At 10 A RMS speaker current, we observe a rapid decline in the deflection magnitude at the peak of the oscillation. We interpret this result as indicating that the density hole displacement has become large enough that the deflection effect saturates. The presence of the density hole creates a transverse gradient in the refractive index, which acts to deflect pulses (see Chapter 2). For small displacements from the center of the density hole, increasing the displacement acts to increase the index gradient sampled by a pulse, which results in a larger deflection. However, this effect 135 eventually saturates and reverses as the density hole displacement becomes sufficiently large (in the limit of very large displacements, a pulse samples unperturbed air and does not deflect at all). Figure 5.1. Filament deflection by 25 Hz acoustic source. Figure a), on the left, shows deflections with a speaker current of 4 A RMS, with figure b) showing deflections with a 10 A speaker current. In addition to dynamic steering of a filamenting pulse train, acoustically modulating the location of the density hole may be useful for quasi-phase matching third harmonic generation in filamentation. Due to dispersion, the phase walk-off length for third harmonic generation is given by 33 L k k       . For example, phase matching in air with an 800 nm pulse has a walk-off length 7L mm  [13]. Figure 5.2, below, shows a scheme where we propose using a solid grate to obstruct sound waves emitted from a speaker. A filamenting pulse train passes over the grate, creating an elongated density hole, local sections of which (those above the sound- transmitting slots of the grate) are driven out of the path of subsequent filaments by 136 the speaker. The combined effect of the density hole, speaker and grate is to create a periodically modulated gas density for the purpose of quasi phase matching third harmonic generation. We should be able to achieve ~1-10% gas density modulation in air at 1 kHz [9,80], with substantially larger modulations achievable by using a more readily ionizable gas such as Xenon [58] or higher gas pressure (which leads to more energy absorbed per pulse), or a higher repetition rate laser (leading to a deeper cumulative density hole). A similar scheme involving a cavity with a longitudinal acoustic mode to produce density variations has been proposed [125], but it is very difficult to generate transducer-driven acoustic waves with amplitudes approaching ~1% [126]. Figure 5.2. Third harmonic generation scheme employing a speaker and grating to periodically modulate the location of the density hole. 5.3 Scaling the thermal air waveguide to a ~50 meter span speaker sound waves grate density hole 137 One of the primary use cases for filament generated air waveguides is guiding secondary beams over large distances in the atmosphere. To this end, an important next step is to test the scalability of the guide over larger distances. Our laser lab is situated at one end of a ~50 meter hallway in our institute, and we have already begun experimenting with ranged filamentation and waveguide generation over the ~50 meter span of the hallway. 5.3.1 Intra-beam phase variation over long filamenting distances Multi-lobe filamenting arrays, like the TEM11 mode used in the 1 meter thermal guide demonstration [10] rely on a relative phase shift of  between neighboring lobes to maintain lobe separation. Recent theoretical studies have suggested that nonlinear phase accumulation in filamenting beams is rapid and highly sensitive to input conditions [127,128]. As a result, it is claimed that the spatial phase of a filament is effectively stochastic, strongly constraining applications which may rely upon maintaining phase coherence across a complicated multi-filamenting beam over many meters of propagation. Results from our initial demonstration of quad filaments [10] already cast some doubt on these claims, as the four-lobe beam maintains lobe separation throughout the initial linear propagation, optical collapse, ~1m of filamentary propagation, and post-filament propagation. This strongly suggests that the inter-lobe phase relationship is maintained over the full pulse evolution. Interferometric measurements of a two-lobe TEM01 filamenting beam propagating ~20 meters down the hallway are shown in Figure 5.3. Images were 138 obtained by interferometrically combining the two lobes of the post-filament pulse on a CCD camera, positioning the two lobes so that their centroids (determined by eye) overlap. This allows for a measurement of the transverse spatial phase difference between the lobes: lobe1 01 01 lobe2 02 02( , ) ( , ) ( , )x y x x y y x x y y       , where 0 0( , )m mx y is the location of the centroid of the thm lobe. The fringe patterns in Figure 5.3(a) are lineouts of the interferograms. The coherence of the fringe patterns indicates that the two neighboring lobes maintain their phase relationship all the way up to 67.5 critical powers ( cr67.5 P ) per lobe (maximum power available in our laser system at the time of the experiment). Given the small but not insignificant shot-to- shot intensity differences between the lobes, references [127,128] would have predicted random variation in their nonlinear phase difference from shot to shot. Figure 5.3(b) plots ,x y N  (blue dots), phase difference between the two lobes averaged over the full transverse beam profile ,x y and a sequence of 100 shots N . It is seen that ,x y N  for the full range of lobe powers. Figure 5.3(b) also quantifies shot-to-shot differences in  , by computing the standard deviation of ,x y  over all shots at a given power:  ,STD x y . The “loss of phase” references [127,128], which argue that the phase in a filamenting beam is extremely sensitive to input peak laser power, would have predicted a dramatic increase in standard deviation in the phase owing to the fact that absolute power fluctuations increase with nominal power. However, what is measured is a relatively flat standard deviation with no evident increase with lobe power. 139 These measurements reinforce the notion that intra-beam phase differences engineered into multi-filamenting beams will be maintained over long propagation distances. Figure 5.3. Phase interferograms collected from interference of the two lobes of a TEM01 filamenting mode after 20 meters of nonlinear propagation. a) Plots of different input power (per lobe) are shown, where lineouts of single shot interferograms are shown in different colors. b) Blue dots plot the transverse spatial phase difference between the two lobes, and red dots show the standard deviation in the transverse spatial phase difference. 5.3.2 Microphone array for single-shot filament visualization 0.125 P cr 2.5 P cr 8.75 P cr 18.75 P cr 31.25 P cr 43.75 P cr 56.25 P cr 67.5 P cr 0 20 40 60 0 2 4 P/P cr p h as e in r ad 2 mm a) b) in te n si ty ( a. u .) 140 Scaling up the thermal guide entails filamenting with a larger mode, and sufficient power so that a multifilamentation pattern can diffusively evolve into a contiguous index cladding. Preliminary results of a multifilamenting beam nucleated from a TEM00 input profile suggest contiguous energy deposition, but have thus far only been “measured” subjectively using the human ear while walking along the beam propagation path. To improve upon this, we have purchased a set of 64 microphones which can be operated synchronously, allowing for a single shot reconstruction of the energy deposition from ~50 meter filament. Further progress on the thermal air waveguide awaits a laser upgrade of our laser from 2 TW to 10 TW, to provide additional power required for creating a contiguous ring-shaped energy deposition on a ~50 meter scale. Figure 5.4. Left: close up of individual microphones, each with local A/D converter. Right: full mounted array of 8 microphones alongside short filament (which is on the other side of the mounting board), showing data hub on the bottom right. 5.4 Linear generation of STOVs using optical elements In Section 4.2 we demonstrated the nonlinear generation of ring STOVs through optical collapse arrest, but also alluded to the possibility of linear generation 141 of STOVs using passive optical elements. A collapsing beam with a preexisting STOV may provide a smoother transition to post collapse filamentation, shedding less energy in the collapse arrest, and perhaps producing a more stable and longer lived filament. Spatial optical vortices have found use in optical tweezing [129,130], where the beam can be used for trapping and transferring angular momentum to micro and nano sized particles, as well as stimulated emission depletion microscopy [131], where the amplitude null is exploited to create images with resolution exceeding the diffraction limit. The same characteristic phase vortex features, which manifest as rapid spatiotemporal variations in a beam with a STOV, may find use in manipulation of relativistic particles as well as for providing rapid temporal signatures which could be useful in time resolved measurement. Here we present a scheme which employs a 4f spatiotemporal pulse shaper [131] to embed a line STOV onto a pulse. To motivate the scheme, consider a pulse with an embedded line STOV BG V( , , )E x y E E  , where we partition the field into a “background” and “vortex” term, mirroring the approach (and symbols) used in Section 4.3.4. Transforming the pulse into spatiospectral space yields an electric field BG V ˆ ˆ ˆ( , y,k )E x E E   , where the “hat” denotes Fourier transform, and the  denotes convolution. We consider a line STOV oriented along the y-axis, so that    V 0 0( , )E x im x x      , and  V 0 0ˆ ( ,k ) ( ) ( )E x im x x k i k            . Thus, in spatiospectral space   BGBG 0 0 ˆ ˆ ˆ( , y,k ) dE E x E im x x i dk          . 142 Considering a Gaussian background field 2 2 2 c2 2 BG x y ik w w E e          , where   2 22 2 c 2 40 BG ˆ 2 k k wx y w E e         , the spatiospectral representation is given by     2BG 0 0 cˆ ˆ( , y,k ) 2 i E x E im x x k k w              (5.1) Starting with a beam without a line STOV, we may impose a STOV by transforming to spatiospectral space, and modifying the amplitude and phase appropriately to impose the bracketed term in Equation (5.1). When centering the line STOV about the temporal axis of the pulse, so that 0 0  , we see that the bracketed term of Equation (5.1) can be interpreted as having a dividing line where the field in spatiospectral space passes through a zero given by   20 c 1 2 x x k k w m     , across which the phase of the field jumps discontinuously by  . A 4f pulse shaper [131] can be used to transform the pulse into spatiospectral space, where it is straight forward to use a programmable spatial light modulator, or static amplitude and phase masks to impose the desired pattern on the beam. Figure 5.5, below, shows a schematic for the 4f pulse shaper: a grating is used to disperse the pulse along one axis, after which a cylindrical lens at a distance f away can be used to spectrally resolve the pulse a distance f beyond the lens. A mask is inserted in the Fourier plane where the pulse is transformed into spatiospectral space (note that the grating and cylindrical lens only act on one transverse dimension of the pulse, leaving 143 the other unchanged). Beyond the mask the sequence is reversed to transform back from the spectral domain. Figure 5.5. Schematic for using a 4f pulse shaper and SLM or static masks to impose a line STOV on a beam. SLM or static mask cylindrical lenses gratings 144 Appendices A.1 Fluid scalings for convective motion of gas driven by pulse train of filaments Starting from momentum and energy balance: 2P t                v v v g v (A.1.1) 2 2 1 1 2 2 P S t                           v v v v q (A.1.2) Here  is the mass density, v is the fluid velocity, P is the pressure, g is gravitational force,  is the coefficient of viscosity,  is the internal energy density, q is the heat flux and S is a heat source. Recall that the heating source comes from the long, thin cylindrical filamenting core. We take the y-direction to point along the axis of the cylinder, making the problem translationally invariant in this direction, let z point upwards (opposing gravity) and x be the remaining transverse direction. Assume a steady state flow 0 t    v with low Reynold’s number 2 1     v v v (i.e. a laminar flow where viscosity dominates convection). Momentum balance in the z simplifies to: 2 0 P g v z          (A.1.3) Where we denote zv v . Assuming hydrostatic pressure equilibrium, 0 P g z      where 0    results in a balance between viscosity and buoyancy from density perturbations. 145 2g v     (A.1.4) The vertical scale length will be much larger than the transverse scale length, as the flow field is oriented along the z-direction zL L which results in an estimate for the terminal fluid velocity based off of scaling arguments: 2gL v    (A.1.4) Where the quantities v and  denote characteristic velocity and density perturbations. To obtain a second scaling we consider the energy equation in steady state and away from the source: 2 1 0 2 P               v v v q (A.1.4) For an ideal gas, the internal energy density is proportional to the pressure . . . 2 D O F P P   , with 3 2   for monatomic gases, and 5 2 for diatomic molecules in the experiment (where two additional D.O.F. are present for rotations, but vibrational modes are still frozen out). Also note that the fluid kinetic energy density can be neglected as it is much smaller than the internal energy density: 2 2 2B 1 s k T P c v m            where sc is the sound speed. We use Fourier’s Law to model conductivity, T  q , where  is the thermal conductivity of the gas. The energy balance simplifies to:   1 0P T     v (A.1.5) 146 Applying a transverse average over the x-direction: L L dx    , notice that ( ) 0xv L  and 0 x L T x     , to give:  1 0 T Pv z z             (A.1.6) The overbar denotes transverse average, and we assume no transverse variation in pressure (as it is mostly varying along the vertical direction). Integrating the equation along z, starting from below the disturbance, which we take to be at z=0 gives:     01 1 T Pv Pv z          (A.1.7) Here we assume that the thermal gradient below the disturbance is zero and that the velocity is 0v . By continuity we can relate the velocities:  0 0 0 0 0 1v v v v                which gives:   0 1 T P v z          (A.1.8) From here, we can arrive at a second scaling relation for the vertical velocity: v L   (A.1.9) Here, 0 pc     is the thermal diffusivity,   B1p k c m   is the specific heat capacity, m is the molar mass, and L is the vertical length scale. We can combine the two relations to get an effective length scale: 147     1 32 eff B1 L L L g mn nk g            (A.1.10) this can in turn be used to estimate the vertical velocity 22 effgLgLv       (A.1.11) 148 A.2 Simulation of filament-induced gas dynamics Simulations of the gas hydrodynamic evolution are performed in cylindrical geometry using a one-dimensional Lagrangian one-fluid hydrocode, in which the conservation equations for mass, momentum and energy,  i i iS t        v  , were solved numerically. For the mass equation, 1  and 1  0 , for the momentum equation, 2  v and 2  I (where I is the unit tensor), and for the energy equation, 2 3 1 2     v and 3 P v q . Here,  is the volume density of the conserved quantity,  is the flux of that quantity, and S refers to sources or sinks, while  is mass density,  is fluid internal energy density, v is fluid velocity, P is gas pressure, and T  q is the heat flux, where  and T are the gas thermal conductivity and temperature. The code includes an artificial viscosity term proportional to 2 v which prevents shock formation, but is not relevant for the neutral gas hydrodynamical response following filamentation, where shock formation does not occur. The radiation of the heated gas contributes negligibly to the energy balance, as can be verified by assuming the maximum emission of a black body and finding it to be tiny. At all times, 1 2 0S S  by mass and momentum conservation, but without approximations, 3 0S  because the thermal part of the energy density is changed by laser heating and by ionization/recombination of all the relevant species in the gas. However, we recognize that at times 10ns after laser filament excitation, all of the 149 energy initially stored in free electron thermal energy and in the ionization and excitation distribution is repartitioned into a fully recombined gas in its ground electronic state. The ‘initial’ radial pressure distribution driving the gas hydrodynamics at times 10ns is set by the initial plasma conditions    e0 e B e g ( ) f P r N r k T r f  , where Bk is Boltzmann’s constant, e ( )N r and e ( )T r are the initial electron density and electron temperature profiles immediately after femtosecond filamentation in the gas, and ef and gf are the number of thermodynamic degrees of freedom of the free electrons and gas molecules. Here, e 3f  , and g 5f  for air at the temperatures of this experiment ( 0 100KT  ). To simulate the neutral gas response at long timescales we solve the fluid equations for the i , using 3 0S  and the initial pressure profile given by 0 ( )P r above. 150 A.3 Scale length for on-axis interference launched from an annular acoustic guide Here we approximate the scale length of the peak on-axis interference launched from light guided in the compressive section of the annular single cycle acoustic wave. Figure A.0.1. Schematic of light launched from end of single cycle annular acoustic guide. D is the size of the annulus, Δr is the thickness of the compressive section of the wave, Δθ is the spreading angle for light launched from the guide, and Δzint is the scale length of the peak on-axis interference from light launched by the acoustic guide. The radial wavenumber spread emanating from the acoustic wave can be estimated from the thickness of the wave 1k r  which can in turn be used to estimate the angular spread 1k k k r      . From here, a simple geometric argument can be used to estimate the region of strongest axial interference intz : tot 2 4 D z    and int tot int 2 2 D z k rD z z        . D Δr Δθ Δθ/4 Δz int Δθ/2 Δz tot 151 Typical values from the annular wave are 500μmD and 25μmr , which for 0.5μm  yields int 20cmz . 152 A.4 Estimating the leakage rate and mode number of a leaky thermal guide Here we will derive the leakage rate and estimate the effective index of the fundamental leaky mode travelling through a thermal air waveguide. To do this, we will first derive the transverse energy flux for the paraxial equation, which we will use to find a general form for the leakage rate of a lossy mode, and then finally apply the general form to the index structure of a thermal air waveguide. The mode index will be estimated using a variational method. We start with the paraxial wave equation valid for modal solutions ( , ) i zA z e r , 2 2 2 2(r , )2 [ ( (r ) )] (r , ) 0 A z i k n A z z             (A.4.1) Where we use the same conventions as established in Section 3.3.2: A is the envelope to the scalar electric field,  is the wavenumber of the guided mode, k is the vacuum wavelength, n is the index, r lies in the transverse plane, and z is the propagation direction. Note that kn is the medium wavenumber,  is the z- component of the medium wavenumber, making 2 2 2k n   the transverse component of the medium wavenumber. Energy continuity equation We wish to derive a continuity equation for the energy, which will take the form of 153 2 z A    S (A.4.2) Where 2 A is proportional to energy, and S is proportional to the transverse energy flux. As the coefficients of proportionality will be irrelevant to the argument, we will neglect them, and refer to the aforementioned quantities as energy and energy flux. To derive S , add * *(A.4.1) (A.4.1)A A   * *1 2 A A A A i        S (A.4.3) Loss rate for leaky mode Assume the envelope has a leaky form: 1 ( , , ) ( , ) zA x y z f x y e N  (A.4.4) Here f is the transverse distribution of the wave, which decays away exponentially at the loss rate  , and 2 N f    normalizes f over the core and cladding of the guide, denoted  . Substituting the form of A given in Equation (A.4.3) into the LHS of the continuity equation (A.4.2) and integrating over  , we have 2 22 zz A e       , while the RHS can be evaluated by using the divergence theorem: cl2S S dA r S         , here clr is the cladding radius, and we have assumed an azimuthally symmetric index profile. Combining these results gives a general form for the loss rate of energy flowing out of a leaky guide: cl * *cl 2 r r r f f f f i N           (A.4.5) 154 Equation (A.4.5) can be simplified further by noting that outside the cladding region (1) 0 ( )f H r , where (1) 0H is a Hankel function physically representing radiation from the cladding of the guide. For strongly bound guides, 1clr , so we may use an asymptotic form for the Hankel function to estimate its derivatives (1) (1) 0 0 ˆ( ) ( ) rH r i H r    , and thus 2* * 2f f f f i f     , resulting in a simplified form for the loss rate cl 2 cl r r f r N       (A.4.6) Solving for loss rate of a thermal guide To obtain a closed form expression to approximate the loss rate of the thermal guide, we approximate the index structure as a radially stepwise function with azimuthal symmetry. As discussed in Chapter 3, the gas evolves via thermal diffusion, enforcing a pressure balance that maintains equal and opposite fractional changes in the gas density and temperature 0 0 T T       , with the change in index proportional to the change in density 0 01 n n       . The thermal energy initially deposited takes the form of a ring array with each lobe depositing energy from individual filaments. As the gas evolves diffusively, the lobes blend together (as seen in Figure 3.2 of Section 3.2.4) becoming indistinguishable from an azimuthally symmetric ring source. The radial dependence of the index perturbation is smoothly varying, but core and cladding radii cor and clr can be estimated by taking moments of the radial distribution, and an effective constant cladding index effn can be defined 155 by enforcing that the integrated index decrement be equal 2 2 eff cladding cladding nd x n d x    . Figure A.0.2. Depiction of approximation for moving from smoothly varying radial index profile (left) to stepwise profile (right). Exact solutions to the stepwise profile shown on the right of Figure A.2 can be obtained analytically by solving the paraxial wave equation over the three domains: clr r , co clr r r  , and clr r , and imposing matching conditions (continuity of first derivative) on the solutions at the boundaries of the three domains. 0 co 1 0 2 0 co cl (1) (2) 1 0 2 0 cl ( ) ( ) ( ) ( ) ( ) J r r r f a K r a I r r r r b H r b H r r r                (A.4.7) Here, 0J is a Bessel function of the first kind, 0K and 0I are modified Bessel functions, and ( ) 0 mH are Hankel functions. While straightforward, this approach necessitates solving a system of transcendental equations, so we will make further simplifications to obtain a solution in terms of elementary functions. Since the loss r Δn r co r cl r Δn r co r cl 156 rate is proportional to the energy density at the edge of the cladding over the total energy in the cladding cl 2 r r f N    we will aim to crudely estimate these quantities by approximating the mode as constant in the core (zeroth order expansion of 0J ) and exponentially decaying in the cladding (asymptotic form of 0K ), while imposing continuity.  co co co co cl 1 r r r r f r e r r r r          (A.4.8) From here it is straightforward to compute  cl cl co 2 2co clcl r r r r r f e r     and   cl cl co22 coco cl 1 r rr N r e         resulting in     cl cl co cl cl co 2 co cl 2 cl co 1 r r r r e k r e              (A.4.9) Here, co co( )r r   and cl co cl( )r r r    are the magnitudes of the transverse component of the medium wavenumber in the core and cladding respectively. Finding the eigenvalue of the leaky mode Equation (A.4.9) specifies the loss rate in terms of the index profile of the guide, but does not specify the eigenvalue  (1)leak 0k n n   , corresponding to the fundamental leaky mode (which appears in the loss rate implicitly through  ). One can view the continuum mode to which the leaky mode corresponds as a resonance condition in the transverse wavevector  allowing for efficient coupling of energy flowing radially inward to be trapped in the leaky. In principle, the fundamental leaky 157 eigenvalue can be found by solving the system of equations (A.4.7) for every value of  , and then searching for resonances in the ratio of the energy density inside and outside of the guide: cl 2 0 2 r r r f f     , with (1) leak corresponding to the largest resonance. As discussed, this approach requires solving systems of transcendental equations, so to circumvent this issue we approximate (1) leak by employing a variational method instead. So far we have been considering modal solutions of the waveguide, which obey equation (A.4.1). More generally, a non-modal solution   ( , ) ( , , ) i kz t E t A x y z e  x obeys the paraxial wave equation  2 2 22 1zik A k n A        (A.4.10) Additionally, non-modal solutions can be constructed by summing over modal solutions   mod , y, z ( ) ( , )eikz i z es A x e A x y d    where  ,A x y are a complete set of modes satisfying orthonormality 2 * 2 ( )A A d x         ,  is the Dirac delta function and   2 * 2A Ad x    is the projection of A onto A . Applying 2 2 *d xA  to both sides of Equation (A.4.10), it is straightforward to establish that     2 2 2 22 2 2 modes 2 ( ) 1k k d k n A A d x          (A.4.11) 158 We can find the eigenvalue (1) leak by considering an appropriate set of trial functions   22 , rx y e      . Note that there are continuum modes of lower order (larger  ) than the lowest order leaky mode. However, because the trial functions decay for large r , they have large projection onto leaky modes, and very small projection onto other continuum modes where energy density is not suppressed outside the guiding structure. More rigorously, one can consider the limit clr  . As the cladding becomes infinite in extent, leaky modes become purely bound, while the continuum solutions interspersed between the bound modes no longer exist. In this limit, the lowest order solution is the lowest order bound mode, which will be very similar in form to the lowest order leaky mode of the corresponding leaky guide in the limit of thick cladding cl 1r . Thus, we will estimate (1) leak by extremizing   2 2 22 2 21k n d x      for the corresponding infinite cladding problem.                2 2 co 2 22 2 2 22 2 2 2 2 co cl co 2 2 2 2 min co co cl2 co (1) 2 2 2 2 2 leak co co co cl2 co min 1 min 1 2 1 ln 2 1 1 1 ln 2 2 r k n d x k n k n n e k r n n r k n k r n n kr                                         (A.4.12) In summary, we have derived an approximation for the loss rate  of the fundamental leaky mode in terms of the index structure and propagation constant (1) leak . 159     cl cl co cl cl co 2 co cl 2 cl co 1 r r r r e k r e              (A.4.13a)    2 2 20 0 eff co2 2 co 1 1 1 1 ln 2 2 2 n n n n k r k r        (A.4.13b) In this last expression, we have expressed the eigenvalue in terms of an index  (1)leak 0k n n   , and have approximated 2 2 co cl 0 eff2n n n n   since 0 1 n n  for thermal air waveguides. 160 A.5 Deriving the equation of motion for the core of a STOV We assume azimuthal symmetry and consider the evolution of the complex envelope  associated with the scalar electric field E ,      , , , , ikz tE r z r z e     (A.5.1) Where 2 k    is the wavenumber, gv t z   is the moving frame longitudinal coordinate, and gv is the group velocity. In the paraxial and slowly-varying envelope approximation, the propagation equation for  is 2 2 2 2 2 2 ( ) 0ik k V z               (A.5.2) Where 2 2 2/ (1/ )( / )r r r       is the transverse Laplacian (assuming azimuthal symmetry) and 2 2 2 2 ( ) cen k c k         is the dimensionless group velocity dispersion (GVD) evaluated at the central frequency cen , and  V  is a nonlinear term. For the case 2 0  , Equation (A.5.2) drives transverse displacement of purely spatial optical vortices through the term 2  , as shown in refs. [38,39]. Here, the term containing 2 additionally drives vortex motion along the local space ( ) direction. Therefore, as the beam propagates, the vortex moves temporally as well. Suppose at 0z the position of the vortex ring is  0 0,vortex rr , and at 0z dz the position of the vortex ring is  0 0,vortex vortexd d r dr    r r . Then 161 0 0 0 0 0 0( , , ) ( , , )r dr d z dz r z dr d dz r z                     . But since 0  at the vortex, this leads to ST vortexd dz z       r (A.5.3) where the spacetime gradient is defined as ˆˆ ST r        r ξ . Following the method of ref. [38,39] as applied to spatial vortices, we approximate the local form of the STOV as a spacetime “R-vortex” 0 0( ) ( )vortex i r r      of charge 1 with a linear phase winding about  0 0, r , embedded in a background field envelope bg such that bg vortex   . If we take i bg e   , where  and  are the real amplitude and phase of the background field, then substitution of bg vortex   into Eq. (1.3) yields 2 2 1 1ˆ ˆ ˆˆˆ ˆ 2 vortexk z r r r                                   r r ξ r ξ ξ (A.5.4) This is Equation (4.6) in the main text. As the pulse propagates along z , Equation (A.5.4) gives the next move of the STOV based on the current self-consistent background field bg . Arriving at equation (A.5.4) by substituting bg vortex   into equation (A.5.3) is straightforward but algebraically tedious. It useful to first express the background field using a Cartesian representation bg u iv   and define a coordinate system centered about the vortex: 0     , 0r r r   . Thus, 162      bg vortex u iv ir u vr i v ur a ib                    (A.5.5) From here, factor the complex pulse evolution equation (A.5.2) into two coupled real equations using a and b . 2 2 2 1 2 za b k        (A.5.6a) 2 2 2 1 2 zb a k        (A.5.6b) Equations (A.5.6) are then combined with equation (A.5.3) to solve for  ,vortexd d drr 2 2 2 2 ST vortex at vortex dz a d b k               r (A.5.7a) 2 2 2 2 ST vortex at vortex dz b d a k               r (A.5.7b) Computing terms individually, one finds ˆ ˆST at vortexa u vr  , ˆ ˆ ST at vortex b v ur   , 2 0 2 rat vortex v a v r           , 2 0 2 rat vortex u b u r            , 2 2 at vortex a u    , and 2 2 at vortex b v    . Which, upon substitution into (1.7) yields 2 1 2 2 2 r at vortex u ud vdr u v dz k r                    (A.5.8a) 2 1 2 2 2 r at vortex v vd udr v v dz k r                     (A.5.8b) Equations (A.5.8) can be used to solve for d and dr : 163   2 2 2 22 2 / 2 r r at vortex dz k u v d u u u v v v v u ru v                     (A.5.9a)   2 22 2 / r r at vortex dz k dr v v v u u u u v u v              (A.5.9b) Finally, we return to the polar form for the background field i bg e   , where cosu   and sinv   . Upon substitution, and after applying various derivatives and trigonometric identities we arrive at 2 1 2 r at vortex dz d k r                 (A.5.10a) 2 r at vortex dz dr k             (A.5.10b) Which can be combined into the single vector equation shown in Section (4.2.5) 2 2 1 1ˆ ˆ ˆˆˆ ˆ 2 vortexk z r r r                                   r r ξ r ξ ξ (A.5.11) 164 A.6 Energy flow in an ultrashort pulse Without loss of generality, we take iue   , where u and  are the real amplitude and phase of the field. Based on conservation of energy we can derive an equation of the form 2u z     j where j is the energy current density in the laser pulse frame and 22u   is the normalized energy density. Following an argument similar to [38], in order to derive j we start with 2 * *( ) . .u c c z z z            and use the pulse envelope equation of motion 2 2 2 2 2 2 ( ) 0ik k V z               to replace z    . Near the vortex, the amplitude of the field is close to zero, making nonlinearities negligible, away from the vortex nonlinearities are only strictly negligible if they are real, resulting in * 1 . . 0 2 V c c ik         . After cancelling out the purely imaginary terms with their complex conjugates, we obtain * 2 2 2 2 1 1 . ( ) 2 c c u u ik k               and 2 2 * 2 22 2 2 2 . 2 c c u u ik k                       , resulting in an energy current density of 2 2 1 ( )u k         j ξ (A.6.1) 165 List of publications by the candidate 1. N. Jhajj, Y.-H. Cheng, J. K. Wahlstrand, and H. M. Milchberg, "Optical beam dynamics in a gas repetitively heated by femtosecond filaments," (2013). 2. N. Jhajj, E. W. Rosenthal, R. Birnbaum, J. K. Wahlstrand, and H. M. Milchberg, "Demonstration of Long-Lived High-Power Optical Waveguides in Air," Phys. Rev. X 4, 11027 (2014). 3. N. Jhajj, J. K. Wahlstrand, and H. M. Milchberg, "Optical mode structure of the air waveguide," Opt. Lett. 39, 6312 (2014). 4. N. Jhajj, I. Larkin, E. W. Rosenthal, S. Zahedpour, J. K. Wahlstrand, and H. M. Milchberg, "Spatiotemporal Optical Vortices," Phys. Rev. X 6, 31037 (2016). 5. Y.-H. Cheng, J. K. Wahlstrand, N. Jhajj, and H. M. Milchberg, "The effect of long timescale gas dynamics on femtosecond filamentation," Opt. Express 21, 4740 (2013). 6. T. I. Oh, Y. S. You, N. Jhajj, E. W. Rosenthal, H. M. Milchberg, and K. Y. Kim, "Scaling and saturation of high-power terahertz radiation generation in two-color laser filamentation," Appl. Phys. Lett. 102, 201113 (2013). 7. T. I. Oh, Y. S. You, N. Jhajj, E. W. Rosenthal, H. M. Milchberg, and K. Y. Kim, "Intense terahertz generation in two-color laser filamentation: energy scaling with terawatt laser systems," New J. Phys. 15, 75002 (2013). 8. J. K. Wahlstrand, N. Jhajj, E. W. Rosenthal, S. Zahedpour, and H. M. Milchberg, "Direct imaging of the acoustic waves generated by femtosecond filaments in air," Opt. Lett. 39, 1290 (2014). 9. E. W. Rosenthal, N. Jhajj, J. K. Wahlstrand, and H. M. Milchberg, "Collection of remote optical signals by air waveguides," (2014). 10. H. M. Milchberg, Y.-H. Chen, Y.-H. Cheng, N. Jhajj, J. P. Palastro, E. W. Rosenthal, S. Varma, J. K. Wahlstrand, and S. Zahedpour, "The extreme nonlinear optics of gases and femtosecond optical filamentation," Phys. Plasmas 21, 100901 (2014). 11. E. W. Rosenthal, J. P. Palastro, N. Jhajj, S. Zahedpour, J. K. Wahlstrand, and H. M. Milchberg, "Sensitivity of propagation and energy deposition in femtosecond filamentation to the nonlinear refractive index," J. Phys. B At. Mol. Opt. Phys. 48, 94011 (2015). 12. D. Kuk, Y. J. Yoo, E. W. Rosenthal, N. Jhajj, H. M. Milchberg, and K. Y. 166 Kim, "Generation of scalable terahertz radiation from cylindrically focused two-color laser pulses in air," Appl. Phys. Lett. 108, 121106 (2016). 13. E. W. Rosenthal, N. Jhajj, I. Larkin, S. Zahedpour, J. K. Wahlstrand, and H. M. Milchberg, "Energy deposition of single femtosecond filaments in the atmosphere," Opt. Lett. 41, 3908 (2016). 167 Bibliography 1. A. Braun, G. Korn, X. Liu, D. Du, J. Squier, and G. Mourou, "Self-channeling of high-peak-power femtosecond laser pulses in air," Opt. Lett. 20, 73 (1995). 2. P. B. Corkum, C. Rolland, and T. Srinivasan-Rao, "Supercontinuum Generation in Gases," Phys. Rev. Lett. 57, 2268–2271 (1986). 3. G. Stibenz, N. Zhavoronkov, and G. Steinmeyer, "Self-compression of millijoule pulses to 78 fs duration in a white-light filament," Opt. Lett. 31, 274 (2006). 4. K.-Y. Kim, J. H. Glownia, A. J. Taylor, and G. Rodriguez, "Terahertz emission from ultrafast ionizing air in symmetry-broken laser fields," Opt. Express 15, 4577 (2007). 5. T. I. Oh, Y. S. You, N. Jhajj, E. W. Rosenthal, H. M. Milchberg, and K. Y. Kim, "Scaling and saturation of high-power terahertz radiation generation in two-color laser filamentation," Appl. Phys. Lett. 102, 201113 (2013). 6. T. Vockerodt, D. S. Steingrube, E. Schulz, M. Kretschmar, U. Morgner, and M. Kovačev, "Low- and high-order harmonic generation inside an air filament," Appl. Phys. B 106, 529–532 (2012). 7. J. Kasparian, M. Rodriguez, G. Méjean, J. Yu, E. Salmon, H. Wille, R. Bourayou, S. Frey, Y.-B. André, A. Mysyrowicz, R. Sauerbrey, J.-P. Wolf, and L. Wöste, "White-Light Filaments for Atmospheric Analysis," Science (80-. ). 301, (2003). 8. P. Rairoux, H. Schillinger, S. Niedermeier, M. Rodriguez, F. Ronneberger, R. Sauerbrey, B. Stein, D. Waite, C. Wedekind, H. Wille, L. Wöste, and C. Ziener, "Remote sensing of the atmosphere using ultrashort laser pulses," Appl. Phys. B Lasers Opt. 71, 573–580 (2000). 9. Y.-H. Cheng, J. K. Wahlstrand, N. Jhajj, and H. M. Milchberg, "The effect of long timescale gas dynamics on femtosecond filamentation," Opt. Express 21, 4740 (2013). 10. N. Jhajj, E. W. Rosenthal, R. Birnbaum, J. K. Wahlstrand, and H. M. Milchberg, "Demonstration of Long-Lived High-Power Optical Waveguides in Air," Phys. Rev. X 4, 11027 (2014). 11. N. Jhajj, I. Larkin, E. W. Rosenthal, S. Zahedpour, J. K. Wahlstrand, and H. M. Milchberg, "Spatiotemporal Optical Vortices," Phys. Rev. X 6, 31037 (2016). 12. J. D. Jackson, Classical Electrodynamics (Wiley, 1999). 13. P. E. Ciddor, "Refractive index of air: new equations for the visible and near 168 infrared," Appl. Opt. 35, 1566 (1996). 14. T. Brabec and F. Krausz, "Nonlinear Optical Pulse Propagation in the Single- Cycle Regime," Phys. Rev. Lett. 78, 3282–3285 (1997). 15. R. W. Boyd, Nonlinear Optics (Academic Press, 2008). 16. J. K. Wahlstrand, Y.-H. Cheng, and H. M. Milchberg, "Absolute measurement of the transient optical nonlinearity in N_2, O_2, N_2_O, and Ar," Phys. Rev. A 85, 43820 (2012). 17. J. P. Palastro, T. M. Antonsen, and A. Pearson, "Models of the delayed nonlinear Raman response in diatomic gases," Phys. Rev. A 84, 13829 (2011). 18. Y.-H. Chen, S. Varma, A. York, and H. M. Milchberg, "Single-shot, space- and time-resolved measurement of rotational wavepacket revivals in H_2, D_2, N_2, O_2, and N_2O," Opt. Express 15, 11341 (2007). 19. S. Varma, Y.-H. Chen, and H. M. Milchberg, "Trapping and Destruction of Long-Range High-Intensity Optical Filaments by Molecular Quantum Wakes in Air," Phys. Rev. Lett. 101, 205001 (2008). 20. S. Varma, Y.-H. Chen, J. P. Palastro, A. B. Fallahkair, E. W. Rosenthal, T. Antonsen, and H. M. Milchberg, "Molecular quantum wake-induced pulse shaping and extension of femtosecond air filaments," Phys. Rev. A 86, 23850 (2012). 21. J. P. Palastro, T. M. Antonsen, S. Varma, Y.-H. Chen, and H. M. Milchberg, "Simulations of femtosecond atmospheric filaments enhanced by dual pulse molecular alignment," Phys. Rev. A 85, 43843 (2012). 22. L. V. Keldysh, "Ionization in the field of a strong electromagnetic wave," JETP 20, (1965). 23. M. V. T. A.M. Perelomov, V.S. Popov, "Ionization of atoms in an alternating electric fiedl," JETP 23, (1966). 24. Y.-H. Chen, S. Varma, T. M. Antonsen, and H. M. Milchberg, "Direct Measurement of the Electron Density of Extended Femtosecond Laser Pulse- Induced Filaments," Phys. Rev. Lett. 105, 215005 (2010). 25. K. Akhtar, J. E. Scharer, S. M. Tysk, and E. Kho, "Plasma interferometry at high pressures," Rev. Sci. Instrum. 74, 996–1001 (2003). 26. A. Houard, V. Jukna, G. Point, Y.-B. Andre, S. Klingebiel, M. Schultze, K. Michel, T. Metzger, and A. Mysyrowicz, "Study of filamentation with a high power high repetition rate ps laser at 1.03 um," Opt. Express 24, 7437 (2016). 169 27. M. Hercher, "Laser-induced damage in transparent media," J. Opt. Soc. Am. 54, (1964). 28. R. Y. Chiao, E. Garmire, and C. H. Townes, "Self-Trapping of Optical Beams," Phys. Rev. Lett. 13, 479–482 (1964). 29. P. L. Kelley, "Self-Focusing of Optical Beams," Phys. Rev. Lett. 15, 1005– 1008 (1965). 30. G. Fibich, The Nonlinear Schrödinger Equation : Singular Solutions and Optical Collapse (n.d.). 31. E. L. Dawes and J. H. Marburger, "Computer Studies in Self-Focusing," Phys. Rev. 179, 862–868 (1969). 32. D. Strickland and P. B. Corkum, "Resistance of short pulses to self-focusing," J. Opt. Soc. Am. B 11, 492 (1994). 33. A. J. Goers, G. A. Hine, L. Feder, B. Miao, F. Salehi, J. K. Wahlstrand, and H. M. Milchberg, "Multi-MeV Electron Acceleration by Subterawatt Laser Pulses," Phys. Rev. Lett. 115, 194802 (2015). 34. G. Fibich, G. Fibich, and B. Ilan, "Vectorial and Random Effects in Self- Focusing and in Multiple Filamentation," Phys. D 157, (2001). 35. A. Couairon and A. Mysyrowicz, "Femtosecond filamentation in transparent media," Phys. Rep. 441, 47–189 (2007). 36. F. Courvoisier, V. Boutou, J. Kasparian, E. Salmon, G. Méjean, J. Yu, and J.- P. Wolf, "Ultraintense light filaments transmitted through clouds," Appl. Phys. Lett. 83, 213–215 (2003). 37. H. R. Lange, A. Chiron, J.-F. Ripoche, A. Mysyrowicz, P. Breger, and P. Agostini, "High-Order Harmonic Generation and Quasiphase Matching in Xenon Using Self-Guided Femtosecond Pulses," Phys. Rev. Lett. 81, 1611– 1613 (1998). 38. P. B. Corkum, "Plasma perspective on strong field multiphoton ionization," Phys. Rev. Lett. 71, 1994–1997 (1993). 39. F. Böhle, M. Kretschmar, A. Jullien, M. Kovacs, M. Miranda, R. Romero, H. Crespo, U. Morgner, P. Simon, R. Lopez-Martens, and T. Nagy, "Compression of CEP-stable multi-mJ laser pulses down to 4 fs in long hollow fibers," Laser Phys. Lett. 11, 95401 (2014). 40. A. J. Campillo, S. L. Shapiro, and B. R. Suydam, "Relationship of self‐ focusing to spatial instability modes," Appl. Phys. Lett. 24, 178–180 (1974). 170 41. A. J. Campillo, S. L. Shapiro, and B. R. Suydam, "Periodic breakup of optical beams due to self‐focusing," Appl. Phys. Lett. 23, 628–630 (1973). 42. A. Couairon and L. Bergé, "Modeling the filamentation of ultra-short pulses in ionizing media," Phys. Plasmas 7, (1999). 43. A. Ting, D. F. Gordon, E. Briscoe, J. R. Peñano, and P. Sprangle, "Direct characterization of self-guided femtosecond laser filaments in air," Appl. Opt. 44, 1474 (2005). 44. Y.-H. Chen, S. Varma, A. York, and H. M. Milchberg, "Single-shot, space- and time-resolved measurement of rotational wavepacket revivals in H_2, D_2, N_2, O_2, and N_2O," Opt. Express 15, 11341 (2007). 45. N. Aközbek, A. Iwasaki, A. Becker, M. Scalora, S. L. Chin, and C. M. Bowden, "Third-Harmonic Generation and Self-Channeling in Air Using High- Power Femtosecond Laser Pulses," Phys. Rev. Lett. 89, 143901 (2002). 46. N. Kortsalioudakis, M. Tatarakis, N. Vakakis, S. D. Moustaizis, M. Franco, B. Prade, A. Mysyrowicz, N. A. Papadogiannis, A. Couairon, and S. Tzortzakis, "Enhanced harmonic conversion efficiency in the self-guided propagation of femtosecond ultraviolet laser pulses in argon," Appl. Phys. B 80, 211–214 (2005). 47. C. P. Hauri, W. Kornelis, F. W. Helbing, A. Heinrich, A. Couairon, A. Mysyrowicz, J. Biegert, and U. Keller, "Generation of intense, carrier-envelope phase-locked few-cycle laser pulses through filamentation," Appl. Phys. B 79, 673–677 (2004). 48. D. C. Smith and D. C., "High-power laser propagation - Thermal blooming," IEEE, Proceedings, vol. 65, Dec. 1977, p. 1679-1714. 65, 1679–1714 (1977). 49. V. V. Vorob’ev, "Thermal blooming of laser beams in the atmosphere," Prog. Quantum Electron. 15, 1–152 (1991). 50. J. K. Wahlstrand, Y.-H. Cheng, and H. M. Milchberg, "High Field Optical Nonlinearity and the Kramers-Kronig Relations," Phys. Rev. Lett. 109, 113904 (2012). 51. M. Takeda, H. Ina, and S. Kobayashi, "Fourier-transform method of fringe- pattern analysis for computer-based topography and interferometry," J. Opt. Soc. Am. 72, 156 (1982). 52. A. Dalgarno and A. E. Kingston, "The Refractive Indices and Verdet Constants of the Inert Gases," Proc. R. Soc. London A Math. Phys. Eng. Sci. 259, (1960). 53. M. Mlejnek, E. M. Wright, and J. V. Moloney, "Dynamic spatial replenishment of femtosecond pulses propagating in air," Opt. Lett. 23, 382 (1998). 171 54. "No Title," http://webbook.nist.gov/chemistry/fluid/. 55. O. Lahav, L. Levi, I. Orr, R. A. Nemirovsky, J. Nemirovsky, I. Kaminer, M. Segev, and O. Cohen, "Long-lived waveguides and sound-wave generation by laser filamentation," Phys. Rev. A 90, 21801 (2014). 56. F. Vidal, D. Comtois, Ching-Yuan Chien, A. Desparois, B. La Fontaine, T. W. Johnston, J.-C. Kieffer, H. P. Mercure, H. Pepin, and F. A. Rizk, "Modeling the triggering of streamers in air by ultrashort laser pulses," IEEE Trans. Plasma Sci. 28, 418–433 (2000). 57. S. Tzortzakis, B. Prade, Y.-B. Andre, M. Franco, A. Mysyrowicz, S. Huller, and P. Mora, "Femtosecond laser-guided electric discharge in air," in Conference Digest. 2000 International Quantum Electronics Conference (Cat. No.00TH8504) (IEEE, n.d.), p. 1. 58. N. Jhajj, Y.-H. Cheng, J. K. Wahlstrand, and H. M. Milchberg, "Optical beam dynamics in a gas repetitively heated by femtosecond filaments," (2013). 59. S. Zahedpour, J. K. Wahlstrand, and H. M. Milchberg, "Quantum Control of Molecular Gas Hydrodynamics," Phys. Rev. Lett. 112, 143601 (2014). 60. K. P. Birch, "Precise determination of refractometric parameters for atmospheric gases," J. Opt. Soc. Am. A 8, 647 (1991). 61. D. R. Lide, Handbook of Chemistry and Physics, 77th ed. (CRC Press, 1996). 62. K. Kim, I. Alexeev, and H. Milchberg, "Single-shot measurement of laser- induced double step ionization of helium," Opt. Express 10, 1563 (2002). 63. H. Gao, W. Chu, G. Yu, B. Zeng, J. Zhao, Z. Wang, W. Liu, Y. Cheng, and Z. Xu, "Femtosecond laser filament array generated with step phase plate in air," Opt. Express 21, 4612 (2013). 64. L. Liu, C. Wang, Y. Cheng, H. Gao, and W. Liu, "Fine control of multiple femtosecond filamentation using a combination of phase plates," J. Phys. B At. Mol. Opt. Phys. 44, 215404 (2011). 65. "No T," http://lasermatter.umd.edu/publications.html. 66. J. Yu, D. Mondelain, J. Kasparian, E. Salmon, S. Geffroy, C. Favre, V. Boutou, and J.-P. Wolf, "Sonographic probing of laser filaments in air.," Appl. Opt. 42, 7117–20 (2003). 67. M. D. Feit and J. A. Fleck, "Light propagation in graded-index optical fibers," Appl. Opt. 17, 3990 (1978). 68. J. Fan, E. Parra, and H. M. Milchberg, "Resonant Self-Trapping and 172 Absorption of Intense Bessel Beams," Phys. Rev. Lett. 84, 3085–3088 (2000). 69. J. Liu, J. Dai, S. L. Chin, and X.-C. Zhang, "Broadband terahertz wave remote sensing using coherent manipulation of fluorescence from asymmetrically ionized gases," Nat. Photonics 4, 627–631 (2010). 70. M. Rodriguez, R. Bourayou, G. Méjean, J. Kasparian, J. Yu, E. Salmon, A. Scholz, B. Stecklum, J. Eislöffel, U. Laux, A. P. Hatzes, R. Sauerbrey, L. Wöste, and J.-P. Wolf, "Kilometer-range nonlinear propagation of femtosecond laser pulses," Phys. Rev. E 69, 36607 (2004). 71. P. Sprangle, J. Penano, and B. Hafizi, "Optimum Wavelength and Power for Efficient Laser Propagation in Various Atmospheric Environments," (2005). 72. P. Sprangle, J. Peñano, B. Hafizi, D. Gordon, and M. Scully, "Remotely induced atmospheric lasing," Appl. Phys. Lett. 98, 211102 (2011). 73. K. W. Fischer*, M. R. Witiw, J. A. Baars+, and T. R. Oke, "Atmospheric Laser Communication," Bull. Am. Meteorol. Soc. 85, 725–732 (2004). 74. P. Panagiotopoulos, N. K. Efremidis, D. G. Papazoglou, A. Couairon, and S. Tzortzakis, "Tailoring the filamentation of intense femtosecond laser pulses with periodic lattices," Phys. Rev. A 82, 61803 (2010). 75. M. Bellec, P. Panagiotopoulos, D. G. Papazoglou, N. K. Efremidis, A. Couairon, and S. Tzortzakis, "Observation and Optical Tailoring of Photonic Lattice Filaments," Phys. Rev. Lett. 109, 113905 (2012). 76. S. Suntsov, D. Abdollahpour, D. G. Papazoglou, P. Panagiotopoulos, A. Couairon, and S. Tzortzakis, "Tailoring femtosecond laser pulse filamentation using plasma photonic lattices," Appl. Phys. Lett. 103, 21106 (2013). 77. R. R. Musin, M. N. Shneider, A. M. Zheltikov, and R. B. Miles, "Guiding radar signals by arrays of laser-induced filaments: finite-difference analysis.," Appl. Opt. 46, 5593–7 (2007). 78. M. Châteauneuf, S. Payeur, J. Dubois, and J.-C. Kieffer, "Microwave guiding in air by a cylindrical filament array waveguide," Appl. Phys. Lett. 92, 91104 (2008). 79. S.-B. Wen, C.-F. Chen, X. Mao, and R. E. Russo, "Guiding and focusing of a nanosecond infrared laser within transient hollow plasma femtosecond filament channels," J. Phys. D. Appl. Phys. 45, 355203 (2012). 80. J. K. Wahlstrand, N. Jhajj, E. W. Rosenthal, S. Zahedpour, and H. M. Milchberg, "Direct imaging of the acoustic waves generated by femtosecond filaments in air," Opt. Lett. 39, 1290 (2014). 173 81. "Supplemental Material," (n.d.). 82. A. W. Snyder and J. D. Love, Optical Waveguide Theory (Springer US, 1984). 83. E. W. Rosenthal, N. Jhajj, J. K. Wahlstrand, and H. M. Milchberg, "Collection of remote optical signals by air waveguides," (2014). 84. T. R. Clark and H. M. Milchberg, "Optical mode structure of the plasma waveguide," Phys. Rev. E 61, 1954–1965 (2000). 85. T. R. Clark and H. M. Milchberg, "Time-evolution and guiding regimes of the laser-produced plasma waveguide," Phys. Plasmas 7, (2000). 86. "An introduction to fluid dynamics. G. K. Batchelor, F.R.S., London (Cambridge University Press), 1967. Pp. xvii, 615; Plates 24; Numerous Figures. 75s. in U.K., $13.50 in U.S.A," Q. J. R. Meteorol. Soc. 94, 435–435 (1968). 87. M. J. H. Ku, W. Ji, B. Mukherjee, E. Guardado-Sanchez, L. W. Cheuk, T. Yefsah, and M. W. Zwierlein, "Motion of a Solitonic Vortex in the BEC-BCS Crossover," Phys. Rev. Lett. 113, 65301 (2014). 88. S. Donadello, S. Serafini, M. Tylutki, L. P. Pitaevskii, F. Dalfovo, G. Lamporesi, and G. Ferrari, "Observation of Solitonic Vortices in Bose-Einstein Condensates," Phys. Rev. Lett. 113, 65302 (2014). 89. G. P. Bewley, D. P. Lathrop, and K. R. Sreenivasan, "Superfluid Helium: Visualization of quantized vortices," Nature 441, 588–588 (2006). 90. M. R. Dennis, Y. S. Kivshar, M. S. Soskin, and G. A. Swartzlander Jr, "Singular Optics: more ado about nothing," J. Opt. A Pure Appl. Opt. 11, 90201 (2009). 91. J. F. Nye and M. V. Berry, "Dislocations in Wave Trains," Proc. R. Soc. London A Math. Phys. Eng. Sci. 336, (1974). 92. L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, "Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes," Phys. Rev. A 45, 8185–8189 (1992). 93. M. W. Beijersbergen, R. P. C. Coerwinkel, M. Kristensen, and J. P. Woerdman, "Helical-wavefront laser beams produced with a spiral phaseplate," Opt. Commun. 112, 321–327 (1994). 94. F. Eilenberger, K. Prater, S. Minardi, R. Geiss, U. Röpke, J. Kobelke, K. Schuster, H. Bartelt, S. Nolte, A. Tünnermann, and T. Pertsch, "Observation of Discrete, Vortex Light Bullets," Phys. Rev. X 3, 41031 (2013). 174 95. D. Mihalache, D. Mazilu, F. Lederer, Y. V. Kartashov, L.-C. Crasovan, L. Torner, and B. A. Malomed, "Stable Vortex Tori in the Three-Dimensional Cubic-Quintic Ginzburg-Landau Equation," Phys. Rev. Lett. 97, 73904 (2006). 96. A. P. Sukhorukov and V. V. Yangirova, "Spatio-temporal vortices: properties, generation and recording," in M. A. Karpierz, A. D. Boardman, and G. I. Stegeman, eds. (International Society for Optics and Photonics, 2005), p. 594906. 97. N. Dror and B. A. Malomed, "Symmetric and asymmetric solitons and vortices in linearly coupled two-dimensional waveguides with the cubic-quintic nonlinearity," Phys. D Nonlinear Phenom. 240, 526–541 (2011). 98. K. Y. Bliokh and F. Nori, "Spatiotemporal vortex beams and angular momentum," Phys. Rev. A 86, 33824 (2012). 99. K. W. Nicholls and J. F. Nye, "The paths of dislocations in wave pulses: an experimental test," J. Phys. A. Math. Gen. 19, 375–383 (1986). 100. Y. V. Kartashov, B. A. Malomed, Y. Shnir, and L. Torner, "Twisted Toroidal Vortex Solitons in Inhomogeneous Media with Repulsive Nonlinearity," Phys. Rev. Lett. 113, 264101 (2014). 101. P. Panagiotopoulos, P. Whalen, M. Kolesik, and J. V. Moloney, "Super high power mid-infrared femtosecond light bullet," Nat. Photonics 9, 543–548 (2015). 102. J. P. Palastro, T. M. Antonsen, and H. M. Milchberg, "Compression, spectral broadening, and collimation in multiple, femtosecond pulse filamentation in atmosphere," Phys. Rev. A 86, 33834 (2012). 103. M. Kolesik and J. V. Moloney, "Nonlinear optical pulse propagation simulation: From Maxwell’s to unidirectional equations," Phys. Rev. E 70, 36604 (2004). 104. J. Andreasen and M. Kolesik, "Nonlinear propagation of light in structured media: Generalized unidirectional pulse propagation equations," Phys. Rev. E 86, 36706 (2012). 105. E. W. Rosenthal, J. P. Palastro, N. Jhajj, S. Zahedpour, J. K. Wahlstrand, and H. M. Milchberg, "Sensitivity of propagation and energy deposition in femtosecond filamentation to the nonlinear refractive index," J. Phys. B At. Mol. Opt. Phys. 48, 94011 (2015). 106. A. S. Desyatnikov, D. Buccoliero, M. R. Dennis, Y. S. Kivshar, M. J. Padgett, and E. A. Cornell, "Spontaneous knotting of self-trapped waves," Sci. Rep. 2, 43001 (2012). 175 107. K. D. Moll, A. L. Gaeta, and G. Fibich, "Self-Similar Optical Wave Collapse: Observation of the Townes Profile," Phys. Rev. Lett. 90, 203902 (2003). 108. T. R. Clark and H. M. Milchberg, "Time- and Space-Resolved Density Evolution of the Plasma Waveguide," Phys. Rev. Lett. 78, 2373–2376 (1997). 109. A. Couairon, J. Biegert, C. P. Hauri, W. Kornelis, F. W. Helbing, U. Keller, and A. Mysyrowicz, "Self-compression of ultra-short laser pulses down to one optical cycle by filamentation," J. Mod. Opt. 53, 75–85 (2006). 110. S. Skupin, G. Stibenz, L. Bergé, F. Lederer, T. Sokollik, M. Schnürer, N. Zhavoronkov, and G. Steinmeyer, "Self-compression by femtosecond pulse filamentation: Experiments versus numerical simulations," Phys. Rev. E 74, 56604 (2006). 111. B. Prade, M. Franco, A. Mysyrowicz, A. Couairon, H. Buersing, B. Eberle, M. Krenz, D. Seiffer, and O. Vasseur, "Spatial mode cleaning by femtosecond filamentation in air," Opt. Lett. 31, 2601 (2006). 112. D. Rozas, C. T. Law, and G. A. Swartzlander, Jr., "Propagation dynamics of optical vortices," J. Opt. Soc. Am. B 14, 3054 (1997). 113. I. S. Sullivan, J. J. Niemela, R. E. Hersherberger, D. Bolster, and R. J. Donnelly, "Dynamics of thin vortex rings," J. Fluid Mech. 609, (2008). 114. M. Kolesik, E. M. Wright, and J. V. Moloney, "Dynamic Nonlinear X Waves for Femtosecond Pulse Propagation in Water," Phys. Rev. Lett. 92, 253901 (2004). 115. D. Faccio, A. Matijosius, A. Dubietis, R. Piskarskas, A. Varanavičius, E. Gaizauskas, A. Piskarskas, A. Couairon, and P. Di Trapani, "Near- and far- field evolution of laser pulse filaments in Kerr media," Phys. Rev. E 72, 37601 (2005). 116. A. Averchi, A. Lotti, A. Couairon, D. Faccio, D. Papazoglou, P. Di Trapani, and S. Tzortzakis, "Ultrashort laser pulse filamentation from spontaneous X Wave formation in air," Opt. Express, Vol. 16, Issue 3, pp. 1565-1570 16, 1565–1570 (2008). 117. S. M. Barnett, "Resolution of the Abraham-Minkowski Dilemma," Phys. Rev. Lett. 104, 70401 (2010). 118. S. M. Barnett and R. Loudon, "The enigma of optical momentum in a medium," Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 368, 927–939 (2010). 119. I. Brevik and I., "Experiments in phenomenological electrodynamics and the electromagnetic energy-momentum tensor," Phys. Rep. 52, 133–201 (1979). 176 120. U. Leonhardt, "Optics: Momentum in an uncertain light," Nature 444, 823–824 (2006). 121. R. N. C. Pfeifer, T. A. Nieminen, N. R. Heckenberg, and H. Rubinsztein- Dunlop, "Colloquium: Momentum of an electromagnetic wave in dielectric media," Rev. Mod. Phys. 79, 1197–1216 (2007). 122. T. G. Philbin and O. Allanson, "Optical angular momentum in dispersive media," Phys. Rev. A 86, 55802 (2012). 123. T. G. Philbin, "Electromagnetic energy momentum in dispersive media," Phys. Rev. A 83, 13823 (2011). 124. K. Y. Bliokh and F. Nori, "Transverse and longtiduinal angular momenta of light," Phys. Rep. 592, 1–38 (2015). 125. U. K. Sapaev, I. Babushkin, and J. Herrmann, "Quasi-phase-matching for third harmonic generation in noble gases employing ultrasound," Opt. Express 20, 22753 (2012). 126. D. J. A. Gallego-Juarez and L. Gaete-Garreton, "Experimental study of nonlinearity in free progressive acoustic waves in air at 20 kHz," Le J. Phys. Colloq. 40, C8-336 (n.d.). 127. B. Shim, S. E. Schrauth, A. L. Gaeta, M. Klein, and G. Fibich, "Loss of Phase of Collapsing Beams," Phys. Rev. Lett. 108, 43902 (2012). 128. A. Sagiv, A. Ditkowski, and G. Fibich, "Loss of phase and universality of stochastic interactions between laser beams," (2017). 129. A. T. O’Neil and M. J. Padgett, Axial and Lateral Trapping Efficiency of Laguerre–Gaussian Modes in Inverted Optical Tweezers (2001), Vol. 193. 130. N. B. Simpson, L. Allen, and M. J. Padgett, "Optical tweezers and optical spanners with Laguerre–Gaussian modes," J. Mod. Opt. 43, 2485–2491 (1996). 131. A. M. Weiner, "Femtosecond pulse shaping using spatial light modulators," Rev. Sci. Instrum. 71, 1929 (2000).