Spatial and temporal pulse propagation for dispersive paraxial optical systems

The formalism for pulse propagation through dispersive paraxial optical systems first presented by Kostenbauder (IEEE J. Quant. Elec. 26 1148-1157 (1990)) using 4× 4 ray-pulse matrices is extended to 6 × 6 matrices and includes non-separable spatial-temporal couplings in both transverse dimensions as well as temporal dispersive effects up to a quadratic phase. The eikonal in a modified Huygens integral in the Fresnell approximation is derived and can be used to propagate pulses through complicated dispersive optical systems within the paraxial approximation. In addition, a simple formula for the propagation of ultrashort pulses having a Gaussian profile both spatially and temporally is presented. © 2016 Optical Society of America OCIS codes: (320.0320) Ultrafast optics; (320.5550) Pulses; (080.0080) Geometric optics; (080.1510) Propagation methods; (080.2730) Matrix methods in paraxial optics. References and links 1. A. G. Kostenbauder, “Ray-pulse matrices: a rational treatment for dispersive optical systems,” IEEE J. Quantum Electron. 26, 1148-1157 (1990). 2. S. Akturk, X. Gu, P. Gabolde, and R. Trebino,“The general theory of first-order spatio-temporal distortions of Gaussian pulses and beams,” Opt. Express 13, 8642-8661 (2005). 3. S. Akturk, X. Gu, P. Bowlan, and R. Trebino,“Spatio-temporal couplings in ultrashort laser pulses,” J. Opt. 12, 093001 (2010). 4. R. Bonifacio, C. Pellegrini, and L. Narducci, “Collective instabilities and high-gain regime in a free electron laser,” Opt. Commun. 50(6), 373-378 (1984). 5. R. Bonifacio, L. De Salvo, P. Pierini, N. Piovella, and C. Pellegrini, “Spectrum, temporal structure, and fluctuations in a high-gain free-electron laser starting from noise,” Phys. Rev. Lett. 73(70), 70-73 (1994). 6. K. J. Kim, “An analysis of self-amplified spontaneous emission,” Nucl. Instrum. Methods A 250(396), 396-403 (1986). 7. P. Emma, R. Akre, J. Arthur, R. Bionta, C. Bostedt, J. Bozek, A. Brachmann, P. Bucksbaum, R. Coffee, F.J. Decker, Y. Ding, D. Dowell, S. Edstrom, A. Fisher, J. Frisch, S. Gilevich, J. Hastings, G. Hays, Ph. Hering, Z. Huang, R. Iverson, H. Loos, M. Messerschmidt, A. Miahnahri, S. Moeller, H.-D. Nuhn, G. Pile, D. Ratner, J. Rzepiela, D. Schultz, T. Smith, P. Stefan, H. Tompkins, J. Turner, J. Welch, W. White, J. Wu, G. Yocky, and J. Galayda,“First lasing and operation of an angstrom-wavelength free-electron laser,” Nat. Photonics 4(9), 641-647 (2010). 8. L. Giannessi, D. Alesini, P. Antici, A. Bacci, M. Bellaveglia, R. Boni, M. Boscolo, F. Briquez, M. Castellano, L. Catani, E. Chiadroni, A. Cianchi, F. Ciocci, A. Clozza, M. E. Couprie, L. Cultrera, G. Dattoli, M. Del Franco, A. Dipace, G. Di Pirro, A. Doria, A. Drago, W. M. Fawley, M. Ferrario, L. Ficcadenti, D. Filippetto, F. Frassetto, H. P. Freund, V. Fusco, G. Gallerano, A. Gallo, G. Gatti, A. Ghigo, E. Giovenale, A. Marinelli, M. Labat, B. Marchetti, G. Marcus, C. Marrelli, M. Mattioli, M. Migliorati, M. Moreno, A. Mostacci, G. Orlandi, E. Pace, L. Palumbo, A. Petralia, M. Petrarca, V. Petrillo, L. Poletto, M. Quattromini, J. V. Rau, S. Reiche, C. Ronsivalle, J. Rosenzweig, A. R. Rossi, V. R. Albertini, E. Sabia, L. Serafini, M. Serluca, I. Spassovsky, B. Spataro, V. Surrenti, C. Vaccarezza, M. Vescovi, and C. Vicario, “Self-amplified spontaneous emission for a single pass free-electron laser,” Phys. Rev. ST. Accel. Beams 14, 060712 (2011). 9. E. Hemsing, A. Knyazik, M. Dunning, D. Xiang, A .Marinelli, C. Hast, and J. B. Rosenzweig, “Coherent optical vortices from relativistic electron beams,” Nat. Phys. 9(9), 549-553 (2013). #259327 Received 11 Feb 2016; revised 25 Mar 2016; accepted 25 Mar 2016; published 1 Apr 2016 © 2016 OSA 4 Apr 2016 | Vol. 24, No. 7 | DOI:10.1364/OE.24.007752 | OPTICS EXPRESS 7752 SLAC-PUB-16500 This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-AC02-76SF00515. 10. M. Xie, “Exact and variational solutions of 3D eigenmodes in high gain FELs,” Nucl. Instrum. Methods A 445(59), 59-66 (2000). 11. J. Feldhaus, E. Saldin, J. Schneider, E. Schneidmiller, and M. Yurkov, “Possible application of x-ray optical elements for reducing the spectral bandwidth of an x-ray sase fel,” Opt. Commun. 140(4), 341-352 (1997). 12. D. Ratner, R. Abela, J. Amann, C. Behrens, D. Bohler, G. Bouchard, C. Bostedt, M. Boyes, K. Chow, D. Cocco, F. J. Decker, Y. Ding, C. Eckman, P. Emma, D. Fairley, Y. Feng, C. Field, U. Flechsig, G. Gassner, J. Hastings, P. Heimann, Z. Huang, N. Kelez, J. Krzywinski, H. Loos, A. Lutman, A. Marinelli, G. Marcus, T. Maxwell, P. Montanez, S. Moeller, D. Morton, H.?D. Nuhn, N. Rodes, W. Schlotter, S. Serkez, T. Stevens, J. Turner, D. Walz, J. Welch, and J. Wu, “Experimental demonstration of a soft x-ray self-seeded free-electron laser,” Phys. Rev. Lett. 114, 054801 (2015). 13. Y. Feng, J. Hastings, P. Heimann, M. Rowen, J. Krzywinski, and J. Wu, “System design for self-seeding the LCLS at soft x-ray energies” in Proceedings of the 34th International Free-Electron Laser Conference, Nara, Japan 26-31 August 2012. 14. Q. Lin, S. Wang, J. Alda, E. Bernabeu, “Spatial-temporal coupling in a grating-pair pulse compression system analysed by matrix optics,” Opt. Quantum Electron. 27(9), 785-798 (1995). 15. S. A. Collins, Jr., “Lens-system diffraction integral written in terms of matrix optics,” J. Opt. Soc. Am. 60, 1168-1177 (1970) 16. A. E. Siegman, Lasers (University Science Books, 1986). 17. Wolfram Research, Inc., Mathematica, Version 10.0, Champaign, IL (2014). 18. R. Haberman, Applied Partial Differential Equations with Fourier Series and Boundary Value Problems (Pearson Education Inc., 2004).


Introduction
The simple and elegant method of Kostenbauder for the description of spatial-temporal couplings has been a powerful tool for modeling the propagation of ultrashort pulses through dispersive optical systems [1].The technique is based on the formulation of ray-pulse matrices that linearly describe the spatial and dispersive properties of an individual optical element or a cascaded series of elements that compose an optical beamline.The vectors that describe the rays take into account both the spatial and temporal variations in the propagating signal independent of the particular beamline components.As such, the matrix that describes the optical system properties has no dependence on the field and can therefore be used to independently evaluate the beamline.This powerful fact allows for the identification of spatiotemporal couplings inherent in a given optical system [2,3].Additionally, these matrices can be used to geometrically propagate rays that carry the spatial and temporal information through an entire beamline in a single step without having to model each individual element independently.A modified Huygens integral can also be constructed to propagate complex fields with arbitrary shape, both spatially and temporally and in both amplitude and phase, through the same beamline within the paraxial approximation in one step.Finally, pulses that have a Gaussian spatial and temporal shape can be propagated using a simple analytic expression.
The Kostenbauder description of optical elements and systems is limited to one spatial dimension (x) along with the temporal dimension (t).Pulses and optical systems with ever increasing complexity, however, will require a complete three-dimensional and nonseparable description that includes both transverse dimensions, (x, y).For example, selfamplified spontaneous emission free-electron laser (SASE FEL) light has significant spatial and temporal structure (see, for example, [4][5][6][7][8][9][10]).These narrow bandwidth pulses are composed of phase uncorrelated but coherent spikes of radiation that are randomly distributed along the longitudinal profile in both the temporal and spectral domains.The transverse mode content may also vary along the temporal profile of the pulse because the gain medium, which is a highly relativistic particle beam composed of free electrons, often has significant longitudinal slice variation of the properties that determine the growth rate of the dominant optical transverse mode.Furthermore, state-of-the-art schemes that are used to create temporal coherence from a SASE FEL often use optical systems that may introduce significant additional spatio-temporal couplings [11][12][13].
An extension of the ray-pulse formalism to both spatial dimensions for the specific case of linearly chirped systems is given in [14].There, however, the optical systems under consideration were only "ideal chirpers" and the pulses under investigation were assumed to be linearly chirped, i.e., the instantaneous frequency of the pulse only had a linear dependence on time.In this paper, the formalism used in [1] is extended to capture generalized threedimensional spatio-temporal couplings up to a quadratic phase.6 × 6 linear ray-pulse matrices are defined from which the eikonal in a modified and time-dependent Huygens integral can be derived and used to propagate pulses with complicated (in both amplitude and phase) spatial and temporal structure in the paraxial approximation.Finally, an analytic and simple propagation formula for three-dimensional Gaussian pulses is presented, which can also be used to analytically estimate the couplings introduced by a given optical beamline.

Ray-pulse matrix formalism
The description of the ray-pulse matrix formalism follows the outline in [1], which is repeated here to draw a direct connection to the original work.It is useful to begin with the definition of an 'optical axis' when using geometric optics to specify a dispersive system that can be described using propagation matrices.Here, the optical axis is defined as the path a reference ray with a given frequency, f 0 , called the reference frequency, takes through the system.For the purpose of this study, all optical elements, such as gratings, lenses, etc., are considered to be aligned to this reference path.A transverse plane with cartesian coordinates (x, y) is constructed orthogonal to the reference ray at each location along the reference path, which is defined to be located at the origin, (0, 0).In addition, clocks are located at each plane and are defined to read zero when the reference ray passes.For any particular section of the beamline, additional rays that propagate through the system can now be studied as deviations from the reference path using the relative coordinates (x, θ x , y, θ y ,t, f ).Here, the θ x,y refer to the slope normalized to the index of refraction a particular ray has relative to the reference path in each cartesian direction and f is the frequency deviation of the ray from the reference frequency.
At the beginning of an optical beamline, z = z i , a ray may be defined relative to the reference path using the six-dimensional vector A similar expression exists for the ray after propagation through the various optical elements that compose the beamline to the location z = z o , If the reference path is a monotonically increasing function in time, meaning the reference ray does not backtrack, the six-dimensional location of the output ray may be parameterized by a vector function of the input coordinates, This transfer map can potentially be a non-linear function of the input coordinates.Considering only small deviations, the function F can be Taylor expanded about the reference path, (0) = (0, 0, 0, 0, 0, 0), F (0) = 0, since this is just the path of the reference ray.Keeping only linear terms as a first approximation produces where are the Jacobian matrix elements of the the 6 × 6 matrix M for the transformation from initial to final coordinates.Equation ( 5) can be explicitly written as The determinant of this matrix gives the ratio of the trace space volume element of the output coordinates to the input coordinates.For the time-invariant (no time-dependent optical elements) and linear (non-linear optical effects such as second harmonic generation are not considered) systems studied here, and for the qualitative reasons given in [1], det (M) = 1.
The matrix in Eq. ( 7) has 36 terms, many of which can be simplified based on the following assumptions and constraints.Only passive optical elements are treated here.This implies that the output coordinates, Y, cannot depend on the time a ray with a particular frequency is launched into the system; only the output time, t o , depends on the input time, t i .This, in turn, requires that the elements in the fifth column of matrix M must be zero except for ∂t o /∂t i = 1.Also, since only linear propagation through dispersive optical systems is treated here the frequency of a particular ray is invariant, implying that the elements in the last row of matrix M must be zero except for where the partial derivatives in M have been replaced in order to represent a generalized propagation matrix.The A, B,C, D matrix elements have been written suggestively, indicating they correspond to the usual A, B,C, D matrix elements of single frequency but non-orthogonal systems (this will be proved later) [15].The E and F elements represent the output spatial coordinate and the output slope dependence on the ray frequency, respectively, while the G and H elements represent the output temporal dependence on the ray spatial and angular input coordinates, respectively.The subscripts indicate the coupling a particular matrix elements exhibits between two coordinates.For example, I t f represents the system on-axis dispersion and relates the output temporal dependence on the ray frequency.The number of terms needed to describe a system has already been reduced to 25.

Construction of the eikonal in a generalized time-dependent Huygens integral
It has been well established that monochromatic (in our notation f = 0 at the reference frequency) paraxial optical waves with a generalized transverse amplitude and phase can be transported through optical systems that are characterized by ABCD matrices using a Huygens integral in the Fresnell approximation that takes the form (see [16] and references therein) where L 0 is the on-axis optical path length through the system and the integral kernel is given by Here, the ray vector coordinates have been organized such that the output ray functional dependence on the input ray is described by or in shorthand notation, r where r = [x, y] and r = [θ x , θ y ].It should be noted that the eikonal function (the argument of the exponential in Eq. ( 10)) is only dependent on the input and output coordinates and not on the angles.Equation ( 9) is valid for both orthogonal (having rotational symmetry where each transverse coordinate can be analyzed independently -the cross terms in the transfer matrix in Eq. ( 11) are zero) and non-orthogonal optical systems.The eikonal function can be constructed by invoking Fermat's principle and equating the general optical path length a ray takes between object and image conjugate points that includes an ABCD optical system to the optical path length a ray takes along the optical axis.Equations ( 9) and (10) indicate that the optical path length a ray takes through the ABCD system includes the on-axis distance, L 0 , and an additional contribution that is quadratic in the displacement r and involves only the ABCD matrix elements of the system.Some issues in applying Eq. ( 9) to time-dependent fields in dispersive optical systems are evident.The path length along the optical axis, L 0 , can be frequency dependent if the lattice contains dispersive elements.The relative phase accumulation between the various frequency components in the field, therefore, would have to be calculated separately from the integral and carried through an entire optical system to obtain the correct longitudinal frequency distribution.Furthermore, the matrix used to propagate each frequency component would be slightly different, requiring an independent integral kernel for each frequency.One goal, therefore, is to find a single integral kernel that incorporates this time and frequency dependence.
To this end, the spatio-temporal analog of the single-frequency Huygen's integral is proposed, based on Eqs. ( 9) and (10), to account for these disparities in a simple way.This integral can be used to propagate time-dependent fields having a narrow bandwidth around the reference frequency, f 0 , through dispersive optical systems that can be described by the matrices given in Eq. ( 8) and include non-orthogonal spatio-temporal couplings (x, y,t) within the paraxial approximation: The integrals in Eq. ( 13) all range from [−∞, ∞].The foundation for this integral relation is discussed in more detail in Appendix A in section 7. We make the following ansatz for the spatio-temporal kernel functional form based on Eq. ( 10): where λ o = c/ f o is a reference wavelength, c is the speed of light, and the eikonal, L, is only dependent on the input and output spatial and temporal coordinates.Expanding the eikonal to second order in these coordinates produces a constant term, L o , which is equivalent to the constant phase factor in Eq. ( 9), terms that are linear in the coordinates, which are zero because of the way the optical axis and reference clocks have been defined, and terms that are quadratic in the coordinates just as in Eq. ( 10).The spatio-temporal Huygens integral can therefore be written as where Δ is needed to ensure energy conservation and T = t i − t o for time-invariant and linear systems.It should be noted that not all 25 real 6 × 6 matrix elements that describe the optical system in Eq. ( 8) can be independent since there are only 15 real coefficients that define the eikonal (e. g.AD − BC = 1 for orthogonal and monochromatic systems).
The correspondence between (α, β , γ,...) and A xx , H ty , I t f , . . .can be established by the procedure outlined in [1].Consider an incident wave with a well defined average position, slope, and frequency and perform the integration in Eq. (15).To this end, we consider the Gaussian input field The output wave involves an exponential of a very complicated second degree polynomial in x o , y o and t o with coefficients that depend on (α, β , γ,...).This polynomial can be separated into an envelope and a phase based on the real and imaginary parts of the argument of the exponential.
The position and time of the output pulse can be found by setting the derivative of the envelope with respect to (x o , y o ,t o ) to zero and solving for (x o , y o ,t o ).The frequency and slope of the output pulse can be found by evaluating the derivative of the phase at the position and time found above.Finally, the coefficients of the averages defined in Eq. ( 16) above correspond to A xx , H ty , I t f , . . . .Intermediate steps in the computation were performed with MATHEMATICA [17] and are omitted for brevity.The results of the involved and rather lengthy computation are summarized in section 8.1, Appendix B. In addition, the derivation of the coefficient Δ is given in section 8.2, Appendix B.
For the matrix elements A xx = A xx (α, β , γ,...), B xx = B xx (α, β , γ,...), etc., the matrix M as defined in Eq. ( 8) satisfies det (M) = 1.In addition, matrix M satisfies the constraint where Equation ( 17) is the symplectic condition (with an appropriate change of coordinates) and can be used to define the constraints on the overdetermined system of equations: The expressions for the matrix elements in section 8.1 along with the constraints listed above can be inverted to write the eikonal in Eq. ( 15) in terms of the matrix elements of the linear transport matrix, M. The expressions for α = α (A xx , B xx , G ty , . ..), β = β (A xx , B xx , G ty , . ..), etc., can be found in section 8.3, Appendix B. Many of these expressions contain the term and others in the denominator that can be problematic when evaluating the integrals in Eq. ( 13) if a particular beamline makes them zero.Special care must be taken when evaluating these systems.This is analogous to the care that must be taken in representing the monochromatic integral kernal, Eq. ( 10), for systems where |B| = 0.The importance of Eq. ( 15) should not be underestimated.It allows for fully three-dimensional fields to be numerically propagated through entire optical systems in a single step rather than resorting to more complicated wave optics propagation methods where single frequency components are propagated individually.

Limiting forms of the spatio-temporal Huygens integral
At this point it is useful to consider various simplified forms of Eq. ( 15) written in terms of (A xx , B xx , G ty , . ..) to ensure the spatio-temporal Huygens integral reduces to the appropriate expressions for specific limiting cases and to gain valuable insight and intuition into the propagation method.Consider, first, single-frequency propagation through a dispersive optical system at a frequency, f , defined relative to the reference frequency.The input field can be written as Plugging this expression into Eq.( 15), evaluating the time integral, and substituting the expressions α = α The eikonal of this kernel has components that are proportional to f and f 2 .The elements that are proportional to f 2 are independent of the spatial input variables.As such, these elements can be pulled out in front of the remaining spatial integrals.Furthermore, some of these terms are proportional to I t f , which is the linear transport matrix element that contains the information regarding the on-axis dispersion and is unique to the spatio-temporal formalism.This formalism, therefore, automatically captures potential frequency dependent path length and phase accumulation dependencies.The terms that are proportional to f (off-frequency) enter the eikonal linearly and are analogous to terms that enter the eikonal linearly in the spatial coordinates in monochromatic but misaligned (off-axis) optical systems.Frequencies that are detuned from the reference frequency, therefore, can be thought of as simple misalignments where the effect on the eikonal is captured with the help of the E, F, G, and H elements.Finally, the reduction of Eq. ( 15) from three to two integrals greatly reduces the computational resources that are needed to propagate fields that often require high spatial and temporal resolution as long as the incident field can be written in the form E i (x i , y i , f ).
Assuming that the optical beamline is spatially symmetric so that transverse components are separable (A xx = A yy = A, A xy = A yx = 0, etc.), the integral kernel in Eq. ( 21) can be further This is the two-dimensional and symmetric analogue of the integral kernel given in [1] for monochromatic spatio-temporal wave propagation, with the exception that the last two terms of the eikonal, 2 f Hλ o (x o + y o ) + 2 f E(x i + y i ), include a missing factor of f , which is needed simply based on dimensional arguments.Finally, analyzing the kernel in Eqs. ( 21) or ( 22) at the reference frequency ( f = 0) reduces the eikonal to the familiar two-dimensional non-separable or separable and symmetric monochromatic case respectively.Therefore, as mentioned previously, the generalized A, B, C, and D elements of the 6 × 6 ray-pulse matrix are identical to those of the usual monochromatic matrix.

Simple propagation formula for generalized three-dimensional Gaussian pulses
A simple propagation formula can be found using the 6 × 6 spatio-temporal ray-pulse matrices for generalized three-dimensional Gaussian pulses that is equivalent in form to both the usual one-dimensional and monochromatic formula (q o = (Aq i + B)/(Cq i + D)) as well as the spatiotemporal two-dimensional result obtained by Kostenbauder ( . This is done by expressing the input field in Eq. ( 15) as a generalized three-dimensional Gaussian pulse and performing the integration.The input field takes the form: where Q i is defined as above so that the field has a quadratic form.After performing the integration it can be shown that the output field is transformed according to which is identical to the transformation found by Kostenbauder except now A, B, C, and D are defined as Equation ( 25) is a very powerful tool.It allows Gaussian fields to be quickly propagated through a potentially complicated optical beamline in a single algebraic step.Furthermore, it allows for a more rigorous definition and identification of spatio-temporal distortions or couplings that can be introduced by a given beamline using the formalism of [2] rather than resorting to computationally costly wave optics numerical simulations.Analytic expressions that explicitly show the functional dependence of spatio-temporal couplings on an optical beamline's parameters should be useful for optical lattice design and optimization.

Conclusion
We have extended the Konstenbauder ray-pulse matrix formalism to include non-separable spatio-temporal couplings with both transverse dimensions.This formalism allows for the three-dimensional (x, y,t) propagation of narrow bandwidth pulses within the paraxial approximation through a time-invariant and linear optical beamline in a single step both geometrically, using 6 × 6 transport matrices, or physically using a modified Huygens integral.The extension of the formalism is important for fields that contain significant and non-separable spatio-temporal couplings or for beamlines that can introduce them.A simple propagation formula for generalized three-dimensional Gaussian pulses was also presented that should be useful in optical beamline design and optimization.
7. Appendix A: Discussion of Eq. ( 13); the spatio-temporal Huygens integral Consider the input-output relationship between a general monochromatic, free-space, sourcefree, linear system (for instance, see [18]) where G describes the impulsive response of the system to a point-source disturbance and L is a linear differential operator (e. g. the linear operator describing the scalar Helmholtz equation: L = ∇ 2 + k 2 ).Assuming the Fourier transform of the input signal, output signal, and impulsive response exist, Eq. ( 27) can be recast as For time-invariant systems this can be re-written as It is possible, therefore, to build up a generalized time-dependent input-output relationship (Eq.( 30)) from a monochromatic relationship (Eq.( 27)) by Fourier superposition, supporting the functional representation given in Eq. ( 13).
(A xx , B xx , G ty , . ..), β = β (A xx , B xx , G ty , . ..), etc. produces the following integral kernal for the output pulse exp { B yy (A yy B 2 xy − A xx B xy B yx + A xx B xx B yy − A xy B xy B yy )x 2 i + B xx B yy (A yy B xx − A xy B yx )y 2 i 2B xx B yy (A xy B yy − A yy B xy )x i y i + B xx B yy (B yy D xx − B yx D xy )x 2 o + B xx (B 2 xy D xx + B xx B yy D yy − B xy (B xx D xy + B yx D yy ))y 2 o + 2B xx B yy (B xx D xy − B xy D xx )x o y o − 2B xx B 2 yy x i x o + 2B xx B xy B yy x i y o + 2B xx B yx B yy x o y i − 2B 2 xx B yy y i y o + κψ − δ μϑ − κσε + μωε H ty = δ κ − β ω −β ψω + β σϑ + δ κψ − δ μϑ − κσε + μωε I t f = λ o (δ μ − β σ) −β ψω + β σϑ + δ κψ − δ μϑ − κσε + μωε 8.2.Derivation of ΔThe field energy is a conserved quantity and must be equal at the beginning and end of a particular beamline.This condition can be expressed as˚|E i (x i , y i ,t i )| 2 dx i dy i dt i = ˚|E o (x o , y o ,t o )| 2 dx o dy o dt o .−xiyi+2εxiti− x i t i + 2ψ y i t i − y i t i + 2x o β x i − x i + μ y i − y i + κ t i − t i + 2y o δ x i − x i + σ y i − y i + ω t i − t i + 2t o −ε x i − x i − ϑ y i − y i − ψ t i − t i dV i dV i dV o , (33)where dV i = dx i dy i dt i , dV i = dx i dy i dt i , and dV o = dx o dy o dt o .The terms in the argument of the exponential that are proportional to x o , y o , and t o result in delta functions (with factors of 2π) when evaluating the integral over dV o .Changing variables to u= 2π λ o β x i − x i + μ y i − y i + κ t i − t i , v = 2π λ o δ x i − x i + σ y i − y i + ω t i − t i , − x i − ϑ y i − y i − ψ t i − t i , (34)will reduce the number of volume integrals by leveraging the delta functions.Solving for x i , y i , and t i ,x i = x i + λ o (u [ψω − σϑ ] + v [−κσ + μω] + z [−κψ + μϑ ]) 2π (β ψω − β σϑ − δ κψ + δ μϑ + κσε − μωε) , y i = y i + λ o (u [δ ϑ − ωε] + v [−β ω + δ κ] + z [κε − β ϑ ]) 2π (β ψω − β σϑ − δ κψ + δ μϑ + κσε − μωε) , t i = t i + λ o (u [−δ ψ + σε] + v [β σ − δ μ] + z [β ψ − με]) 2π (β ψω − β σϑ − δ κψ + δ μϑ + κσε − μωε) , (35)and relating the differential volume elements dx i dy i dt i = |J| dudvdz, where the Jacobian determinant is given by i (x i , y i ,t i ) E * i x i , y i ,t i exp −