Dispersion in time and space: what propagating optical pulses in time (¬ space) forces us to confront

I derive a temporally propagated uni-directional optical pulse equation valid in the few cycle limit. Temporal propagation is advantageous because it naturally preserves causality, unlike the competing spatially propagated models. The exact coupled bi-directional equations that this approach generates can be efficiently approximated down to a uni-directional form in cases where an optical pulse changes little over one optical cycle. They also permit a direct term-to-term comparison of the exact bi-directional theory with its corresponding approximate uni-directional theory. Notably, temporal propagation handles dispersion in a different way, and this difference serves to highlight existing approximations inherent in spatially propagated treatments of dispersion. Accordingly, I emphasise the need for future work in clarifying the limitations of the dispersion conversion required by these types of approaches; since the only alternative in the few cycle limit may be to resort to the much more computationally intensive full Maxwell equation solvers.


I. INTRODUCTION
The importance of having simple and robust methods for the propagation of optical pulses has undergone increasing attention in recent years. This is due to the multitude of applications [1] in which increasing short pulses are relied upon to act as a kind of strobe-lamp to image ultrafast processes [2,3], or even designer electric field profiles [4][5][6] to excite specific atomic or molecular responses. Alternatively, strong nonlinearity is used to construct equally wide-band but also temporally extended pulses -i.e. (white light) supercontinnua [7][8][9] -or even come full circle and use the strong nonlinearity to generate sub-structure that is again temporally confined, as in optical rogue waves [10]; or even the temporally and spatially localized filamentation processes [11][12][13].
It is clear, therefore, that progress towards diminishing pulse durations as well as their increasing spectral bandwidths, and higher pulse intensities are all either stretching existing pulse propagation models to their limits, or breaking them. In such regimes, we need to be sure that our numerical models still work, and have a clear idea of what has been neglected, and what the side-effects of those approximations are. In this paper I use a modern directional approach to show how a straightforward and relatively simple derivation allows a side-by-side comparison of exact and approximate propagation equations, whilst still providing the numerical and analytical convenience of a first-order wave equation.
However, unlike the typical approach in nonlinar optics, I consider propagation in time rather than along some chosen spatial axis [14,15]. The resulting evolution equations for the wave profiles in space no longer give us access to the frequency spectrum [16], but instead we have the wavevec- * Electronic address: Dr.Paul.Kinsler@physics.org tor spectrum, as well as a causally appropriate time evolution [17]. This requires an alternate approach to the temporal response of the propagation medium, which highlights the existing and somewhat poorly characterised approximations inherent in spatially propagated treatments of dispersion.
I start in section II by deriving a customized second order wave equation from Maxwell's equations, and reorganize it to define the material properties appropriately and set up the factorization stage. In section III, I use the method of factorization that allows us to construct an explicitly bidirectional model, and which is then reduced to the popular uni-directional limit in section IV. Section IV D discusses typical modifications that can be applied to the equations given in sections III and IV in order to and simplify them appropriately and compare them to existing models. Having derived this main result -the propagation equations and their various specialisations -I turn in V to the consequences of temporal propagation and the new perspective it provides on the handling of dispersive effects in this and standard (spatially propagated) models. This is then followed in section VI by a specific comparison between the ordinary spatially propagated nonlinear Schrödinger (NLS) equation -perhaps the most widely investigated equation in nonlinear optics -and its temporally propagated counterpart. Finally, in section VII, I present my conclusions.

II. SECOND ORDER WAVE EQUATION
Most optical pulse problems consider a uniform and source free dielectric medium. In such cases a good starting point is the second order wave equation for the electric field, which results from the substitution of the ∇ × H = ∂ t D + J Maxwell's equation into ∇ × E = −∂ t B (see e.g. [18]). Here, however, I want to follow the displacement field D (not E) since it naturally occurs in conjuction with a time derivative -likewise, Dr.Paul.Kinsler@physics.org http://www.kinsler.org/physics/ the other important field is B rather than H. Thus we rexpress the curl Maxwell's equations as Note that for both of these equations, we will also (as usual) want to specify the material response, whether dielectric or magnetic, in terms of a dielectric polarization P = D − ε 0 E and a magnetization M = B − µ 0 H. Either of these equations given above might be either in the spatial domain (with argument r = (r x , r y , r z )), or the wavevector domain (with argument k = (k x , k y , k z )). Now we define the total magnetization M in terms of a component M L linearly dependent on B, and another component M B containing the rest. As a result, the excess polarization M µ ≡ M B , which will typically include any nonlinearity, is also a function of B and not H. In traditional H-based spatially-propagated approaches, the temporal response ("dispersion") can be incorporated naturally, here it cannot. We have that which in the spatial frequency (wavevector k) domain is since after a Fourier transform from r into k, the products become spatial convolutions ("⋆"); which emphasises that anything wavevector dependent is necessarily a non-local concept. Note that this is the converse of the usual situation where a temporal response is described using time-domain convolutions, and these become products in the frequency domain. Now we define the total polarization P in terms of a component P L linearly dependent on D, and another component P D containing the rest. As a result, the excess polarization P ε ≡ P D , which will typically include any nonlinearity, is also a function of D. We have that P(r) = P L (r) + P D (r) = α ε (r)D(r) + P D (r), (2.5) which in the spatial frequency (wavevector k) domain is With these definitions for ∂ t D, ∂ t B, and for M and P set up, we can now proceed with the derivation of a second order wave equation for D - Now take the time derivative of this equation, use the fact that α µ is independent of time, define a µ = 1 − α µ and a ε = 1 − α ε ; and then substitute for ∂ t B, so that Remember that in this picture, all field properties are specified as functions of our chosen primary field D; notably, the magnetic field B which is usually the main argument for M B must be determined from the calculated D using the relevant Maxwell's equation (i.e. eqn. (2.2)).

A. Simple case: isotropic & homogeneous
Since trying to treat all the details of eqn. (2.15) correctly leads to excessively complicated expressions, I will first simplify it down into the case where (a) the material's reference properties a µ , a ε do not vary in space, and (b) the effective monopole current is zero (K = 0). Since (magnetic) monoples do not exist, and all non reference properties can still be encoded in P D , M B , or J, this simplification does not (of itself) impose any approximation. However, the first simplification means that our reference properties cannot incorporate any exact knowledge about the spatial structure of the propagation medium 1 . This restriction on obtaining the best possible match between the reference behaviour and the exact behaviour of the propagation medium will typically affect later approximations, where we assume deviations from the reference behaviour are small. One additional advantage of selecting constant reference parameters α µ and α ε is that we can replace them with matrices to allow for e.g. birefringence and cross-polarization couplings.
The simplified second order wave equation for changes in the displacement field D(r) as it propagates forward in time Note that the retention of the ∇ · D and ∇ · P D terms allows us to include charge effects, as is needed when discussing highpower nonlinear situations such as filamentation [19].
In the spatial Fourier "wavevector" regime, we replace ∇ → ık, so that k 2 = k · k, and reorganize slightly. The equation for changes in the displacement field D(k) as it propagates forward in time is where we have defined a reference frequency Ω for the propagation, based on the reference material properties, i.e.
We could -if we believed we understood the spatial properties of the propagation medium sufficiently well -simply alter this definition by incorporating an alternate k dependence that mimics an assumed spatial dispersion. However, it needs to be understood that such a step is distinctly ad hoc, since it lacks the convolutions over k necessary for any accurate representation of the material structure. But, that said, such an approach may still have considerable practical utility; and indeed is functionally equivalent to how other methods handle spatial dispersion, as a non-local (and hence, strictly speaking, a non-causal) process: e.g. those that import a waveguide dispersion calculated from the (spatial) modes given by the guide's transverse cross-section. To allow for this possibility, in the following I will replace the physically rigorous magnitude k argument for Ω with k, its vector counterpart 2 . This allows for not only more complicated assumed spatial dispersions, but also additional orientation dependence. However, these benefits (only) arise as a result of approximation. Strictly speaking, both on the grounds covered above, and ones taking stricter view of causality [16], any notion of spatial dispersion is incompatible with locality. Unlike for the spatially propagated scheme in [14], we have not included temporal response of the propagation medium in our definition of the directional fields, since it is incompatible with our temporal propagation. To include this we must model it explicitly as is done in Maxwell solvers such as FDTD [20,21]; or approximate it as if it were a spatial effect (see e.g. [22], or the discussions in the later section V.

III. FACTORIZATION
I now factorize the wave equation [14,23], a process which neatly avoids some of the approximations necessary in traditional approaches. Factorization takes its name from the fact that the LHS of eqn. (2.2) or (2.4) is a simple difference of squares which might be factorized, indeed this is what was done in a somewhat ad hoc fashion by Blow and Wood in 1989 [24]. Since the factors are just of the form (∂ t ∓ ık), each by itself looks like a forward (or backward) directed wave equation. It therefore allows us to define a pair of counterpropagating Greens functions, and so divide the second order wave equation into a bi-directional pair of coupled first order wave equations. In addition, it allows us to directly compare the exact bi-directional and approximate uni-directional theories term for term, whereas in other approaches the backward parts simply vanish and are not directly available for comparison. A short summary of the factorization method is given in appendix A.
To re-iterate an important point -the choice of a ε , a µ and therefore Ω(k) in eqn. (2.4) defines the specific Greens functions used; it therefore also defines the basis upon which we will then propagate the displacement field D.

A. Bi-directional wave equations
A pair of bi-directional wave equations suggests similarly bi-directional fields, so I split the electric displacement field into forward (D + ) and backward (D − ) directed parts, with D = D + + D − . In the following equations, remember that P D , M B , etc all depend on the full field D, and not only one partial field (e.g. only D + ). This is an important point since we see they drive both forward and backward evolution equations equally.
The coupled bi-directional first order wave equations for the directed fields D ± are propagated forward in time whilst Here I have replaced the partial time derivatives on the RHS with over-dots, this being to emphasise that if our expression is to remain strictly causal [17], then our models for current J and magnetization M B must return their time-derivatives as explicit functions of known quantities (here, typically D or t).
Note that factorized equations can be rebuilt into a single second-order equation by taking the sum and difference, then substituting one into another with the assistance of a derivative w.r.t. t (see section IV.B of [23]).
A remaining clarification required is to understand the meaning of "forward" and "backward" in the context of three dimensional space -this being less straightforward than in the spatially propagated case, where one spatial axis is selected 3 , and forward and backward refers to time t. Here, we should understand "forward" and "backward" as being with respect to a particular choice of k, and those individual vector components, that we are interested in.

B. Propagation, evolution, and directed fields
Note that since our solutions of the wave equations enforce propagation towards later times t, the fields D ± (t) are directed forwards and backwards in space; these fields then evolve forwards and/or backwards in space as t increases. Note that I use this terminology (propagated, directed, evolved) throughout this paper to mean these three specific things.
When examining the wave equation eqn. (3.2) which evolves the directed fields D ± as they propagate forward in t, we see that the RHS has two types of terms: which I label the "reference" and "residual" parts [25].
Reference evolution (or "underlying evolution") is that given by ±ıΩ(k)D ± term, and is determined by our chosen a ε and a µ , and any additional ad hoc refinements. By itself, it would describe an ordinary oscilliatory evolution where the 3 Typically a cartesian axis, although other choices are possible [15]. field oscillations in D ± (r) would move forward (+) or backwards (−) in space. This is analogous to the choice of reference when constructing directional fields [26], but here is done for temporal rather than spatial propagation [27].
Residual evolution accounts for the discrepancy between the true evolution and the underlying/reference evolution, and is every part of the material response not included in a ε , a µ (and hence Ω(k)). i.e. it is all the remaining terms on the RHS of eqn. (3.2). These typically include any non-linear polarization, or spatially or orientationally dependent linear terms; they are analogous to the correction terms used in directional fields models. In the language used by Ferrando et al. [23], these residuals are "source" terms. Although we might hope they will be only a weak perturbation, so that we might make the (desirable) uni-directional approximation discussed later, the factorization procedure is nevertheless valid for any strength.
C. Reference evolution: choice of Ω and the resulting D ± I now examine how the choice of Ω affects the relative sizes of the forward and backward directed D + and D − . To do this I consider the simple example of a medium for which the field is known to propagate with frequency ω; but for demonstration purposes we choose a reference evolution determined by a frequency Ω that is different from ω. E.g., for a linear isotropic medium we could exactly define ω 2 = Ω 2 + ∆ 2 ; but in general we would just have some residual (source) term Q. This means that our definitions of forward and backward directed fields do not exactly correspond to what the wave equation will actually evolve forward and backward as we propagate towards larger t.
The second order wave equation is (∂ 2 t + Ω 2 )D = −Q, which in the linear case has Q = ∆ 2 D, so that (∂ 2 t + ω 2 )D = 0. The factorization in terms of Ω is then Now if we select the case where our field D only evolves forward, we know that D = D 0 exp [ıωt]. Consequently D ± must have matching oscillations: i.e. D ± = D ± 0 exp[ıωt], even though D − is directed backwards. Substituting these into eqn. (3.1) gives which specifies how much D − we need to combine with D + so that our pulse evolves forward; since the D − will be dragged forward by its coupling to D + . This interdependent D ± behaviour is generic -no matter what the origin of the discrepancy between Ω and the true evolution of the field (i.e. the residual/source terms such as mismatched linear behaviour, nonlinearity, etc): some non-zero backward directed field D − must exist but still evolve forwards with D + . Analogous behaviour can be seen in the directional fields approach of Kinsler et al. [26].

D2OWE
Time propagated optics and dispersion Dr.Paul.Kinsler@physics.org http://www.kinsler.org/physics/ Usually we hope that this residual D − contribution is small enough so that it can be neglected. If we assume D − ≃ 0, then we find that ω ≃ Ω + ∆ 2 /2Ω, which is just the expansion of ω = (Ω 2 + ∆ 2 ) 1/2 to first order in ∆ 2 /Ω 2 . Following this, we find that eqn. (3.2) then says that D − 0 ≃ (∆ 2 /4Ω 2 )D + 0 , which has come full circle and provided us with the scale on which D − can be considered negligible. Outside the restricted (linear) case where we know ∆ 2 , the true frequency ω might be difficult to determine, and in nonlinear propagation any local estimate of the frequency can change as the pulse propagates.
There is a further important point to notice: if we choose Ω = Ω(k) with a wave-vector dependence, then we see that the source-like terms inherit that dispersion; and this is almost inevitable, since typically Ω ∼ ck. This means, for example, that even if we started with nonlinear polarization model without any k dependence, our factorized equations will have a k dependence on the nonlinear terms. This is the complementary process to what happens in the spatially propagated directional fields approach of Kinsler et al. [26], where choosing a dispersive reference has an equivalent effect on the correction terms.

IV. UNI-DIRECTIONAL WAVE EQUATIONS
Making only a single, well defined class of approximation I can now reduce the exact coupled bi-directional evolution of D down to a single uni-directional first order wave equation. As with the complementary spatially propagated theory [14], I do not require a moving frame, a smooth envelope, or to assume inconvenient second order derivatives are somehow negligible, as was necessary in traditional treatments. The approximation is that the residual terms are weak compared to the (reference) ±ıΩD term -e.g. weak nonlinearity or orientation dependence. This enables me to assert that if I start with D − = 0, then D − will remain negligible. In this context, "weak" means that no significant change in the backward field is generated in a distance shorter than one wavelength ("slow evolution"); and that small effects do not build up gradually over propagation times of many periods ("no accumulation").
Slow evolution is where the size of the residual terms is much smaller than that of the underlying linear evolutioni.e. smaller than ΩD. This allows us to write down straightforward inequalities which need to be satisfied. It is important to note the close relationship between these and a good choice of Ω, as discussed in subsection III C. If Ω is not a good enough match, there always be significant contributions from both forward and backward directed fields; and even if nothing ends up evolving backwards, an ignored backward directed field will result in (e.g.) miscalculated nonlinear effects, since the total field D = D + + D − will be wrong.
No accumulation occurs when the evolution of any backward directed field D − is dominated by its coupling via the residual terms to the forward directed field D + ; and not by its preferred underlying backward evolution. No accumulation means that forward evolving field components do not couple to field components that evolve backwards; this the typical behaviour since the phase mismatch between forward evolv-ing and backward evolving components is ∼ 2Ω; in essence it is comparable to the common rotating wave approximation (RWA). This rapid relative oscillation means that backward evolving components never accumulate, as each new addition will be out of phase with the previous one; it is not quite a "no reflection" approximation, but one that asserts that the many micro-reflections do not combine to produce something significant. An estimate of the terms required to break this approximation are given in appendix B; generally speaking this is a much more robust approximation than the slow evolution one. Of course, a periodic temporal modulation of the medium will give periodic residual terms, and these can be engineered to force phase matching. In spatially propagated analyses, where we talk of spatial (wavevector) phase matching and not this temporal (frequency) phase matching [28], this would be a periodicity based on a relatively small phase mismatch (see e.g. quasi phase matching in Boyd [29]); but might even go as far as matching the backward wave (see e.g. [30]).
It is also important to note that the same small size of perturbation from the residual terms can accumulate on the forward evolving field components (or, indeed, the backward perturbation on the backward evolving field components). Although the magnitude of the residual terms acting on the forward and backward field evolution are identical, forward evolving components of the residuals can accumulate on the forward evolving field because they are phase matched; whereas backward residuals are not, and rapidly average to zero.

A. Residual terms
The handling and criteria for different types of residual terms has been discussed already [14] for the spatially propagated case. While it would be straightforward to repeat that level of detail here, in fact the basic process is virtually identical except for the different roles played by the t and r arguments. Accordingly I will consider only the most commonly considered residuals -the dielectric scalar and vector terms as well as currents. Those interested in magnetic or magnetoelectric effects can adapt the process to their own systems in a similar manner.
The relevant total polarization P D can then be decomposed into pieces before seeing how each might satisfy the slow evolution criteria. Thus we write P D (D, r,t) = a ε (D, r,t)D(r,t) + U D (D, r,t) = a L (r,t) * D(r,t) + a N (D, r,t) * D(r,t) + U L (r,t) * D(r,t) + U N (D, r,t).
The part which is scalar in nature is represented by a ε , and might contain (non-reference) linear parts and time response (a L ), perhaps a convolution (" * ") over the past, or angle dependence; it might also contain nonlinear contributions a N such as the third order Kerr nonlinearity with a N D ∝ (D· D)D.
The vector part U D , if present, could be from e.g. a second order nonlinearity, which couples the ordinary and extraordinary field polarizations. Note that this description of the D2OWE Time propagated optics and dispersion Dr.Paul.Kinsler@physics.org http://www.kinsler.org/physics/ material parameters does not restrict them in any way, and can e.g. include any order of nonlinearity.

B. Slow evolution?
Now I will treat each possible residual term in order, where the oppositely directed field is negligible: i.e., for D ± , we have that D ∓ ≃ 0. where the scalar ε L contains the linear response of the material that is both isotropic and lossless (or gain-less); since here it is a time-response function, it is convolved with the electric displacement field D. Note that the field vector D, and indeed the material parameter ε L are all functions of time t and space r = (x, y, z); the polarization P ε and its components φ ε , V ε are a functions of time t, space r, and the field D.
Below I will refer to field components D i and D ± i , where D + + D − = D ≡ (D x , D y , D z ) and i ∈ {x, y, z}; also to wave vector components k i from k = (k x , k y , k z ), with k 2 ⊥ = k 2 x + k 2 y . However, note that in the constraints below, that k j is also used as a substitute symbol to represent any one of k x , k y , or k z .
Any corrections due to imperfect choice of Ω will appear in the P D polarization correction term. However, unlike in the spatially propagated case, there is no diffraction-like termhere, that is quite naturally part of the reference behaviour.
Firstly, scalar polarization terms φ , which can be either linear (φ L (k)) or nonlinear (φ N (D, k)). These might encode e.g. the time-response of the medium (and hence its dispersion), or be some nonlinearity response such as the third-order Kerr nonlinearity. Since these are locally correlated in space, in the wavevector picture they include convolutions over k, and so give us the criterion In the linear case, φ ≡ φ L is independent of D, so only the material parameters are constrained and the pulse properties play no role. In the nonlinear case, the criteria makes demands on the peak intensity of the pulse -but does not apply smoothness assumptions or bandwidth restrictions.
Secondly, linear and nonlinear vector terms from U. These will have a criterion broadly the same as the scalar cases in eqn. (4.1), but with U replacing φ D. Thus for i ∈ {x, y, z}, we can write down constraints for each component U i , which are In the linear case, U ≡ U L , and since U L and D have some linear relationship, this criterion only constrains the material parameters contained in U L , not the field profile or amplitude. The nonlinear case is similar to the scalar one, and the peak pulse intensity is restricted; e.g. for a χ (2) medium, |U N | ∼ χ (2) |D|. However, one complication of the vector cases is that a field consisting of only one field polarization component (e.g. D + x ) may induce a driving in the orthogonal (and initially zero) components (e.g. D ± y ). Hence both D ± y fields will be driven with the same strength, so that it is far from obvious that we can set D − y to zero, but still keep the D + y without being inconsistent. However, as described above, it is the phase matching which ensures that forward residuals accumulate, whilst the non-matched backward residuals are subject to the RWA, and become negligible: hence we can still rely on eqn. (4.2), albeit under caution.
Thirdly, we have the residual terms indicated by the presence of time-varying currents. These are Fourthly, we have the divergence terms kk · D, kk · P D , which amount to charge density contributions either from D or P D , with σ P (k) = a −1 ε k · P D and σ D (k) = k · D. Thus if the wavevector is oriented along a unit vector u, then is the relevant criterion.
To summarize, these various criteria assert that modulations away from the reference evolution must be weak. Notably, weak nonlinearity is invariably guaranteed by material damage thresholds [31]. However, and in contrast to the usually favoured (but strictly speaking less causal) spatially propagated methods, the lack of any rigorous way to incorporate dispersion into the reference frequency Ω(k) makes it harder to satisfy the criteria for the ordinary linear response.

C. Uni-directional equation for D +
In the case where all of the slow-evolution criteria listed above hold, we can be sure that the backward directed field D − is negligible, and if the no-accumulation condition also holds, then neither will there be any non-negligible backward evolving contributions to the field. Consequently, we can be sure that an initially negligible D − remains so, and eqn.(3.2) becomes where now the polarization P D , diffraction, and divergence are solely dependent on the forward directed field D + .

D. Modifications
Just as for the comparable spatially propagated derivation [14] we might also apply some of the strategies used in other D2OWE Time propagated optics and dispersion Dr.Paul.Kinsler@physics.org http://www.kinsler.org/physics/ approaches to get a more tractable and/or simpler evolution equation. Unlike in "traditional" derivations [32][33][34][35], none of these are required, but they nevertheless may be useful.
Briefly, a co-moving frame might be added, using z ′ = z + tv f . This is a simple linear process that causes no extra complications; the leading RHS ıcD + term is replaced by ı(c ∓ v)D ± . Note that setting c = v will freeze the phase velocity of the pulse, not the group velocity. Also, a carrierenvelope separation could be implemented using − k 1 z)] defining the envelope A(z) with respect to wave vector k 1 and carrier frequency ω 1 . Multiple envelopes centred at different wave vectors and carrier frequencies (k i , ω i ) could also be used to split the field up. This is typically done when the field principally comprises a number of distinct narrowband components (e.g. [29,36,37]). The wave equation can then be separated into one equation for each piece, coupled by the appropriate wavevector-matched polarization terms (cf [28]).

V. DISPERSION
Consider a dielectric slab waveguide of thickness d in the perpendicular x direction, propagation in the z direction, and with core and cladding permittivities ε 1 and ε 0 . In the CW limit, the spatial structure of the waveguide gives rise to frequency dependent properties that are usually called "spatial dispersion". Even for this simple slab design, the dispersion relation has a non trivial form; in this case it is given by the solution to a transcendental equation [38]. Typically, this is written in a way implying we want to calculate k x and β = k z from ω, i.e. for the transverse electric (TE) field modes we have and T (k x d) is either tan(k x d) or − cot(k x d).
However, we can now see that this is an implicitly spatially propagated representation, since in time propagated systems we do not (cannot) know ω, but we do know k. Therefore eqn. (5.2) should be rewritten accordingly, using k z for β (since wavevector is no longer a propagation reference), and Ω for ω (since frequency now is a propagation reference); thus Here we see that the natural (reference) propagation frequency is defined solely by the transverse wavevector k x , as is the matching evolution wavevector k z . While apparently simpler than the traditional form, since Ω is directly given by k x , we usually want to neglect any direct knowledge of the transverse properties. This means that what we really want is not Ω(k x ), but Ω(k z ), the inverse of the traditionally preferred β (ω). Despite these complications, which are in any case rather similar to those we see in the traditional spatially propagated picture, we can still find the same waveguide modes and their dispersions, but this time given as functions of k z , not ω. So for some n-th mode of the waveguide, as expanded about some central evolution wavelength k zn , which corrresponds to a reference frequency Ω n = Ω(k zn ), we might have where γ 1 = c +γ 1 has two contributions -that due to the vacuum (just c), and the estimated modifications due to spatial effects (γ 1 ). We may also wish to add the effects of whatever temporal response the material might have to our propagation calculation; but also decide that an approximate model where we modify the expression above in eqn. (5.5) is sufficent. This is the converse of the traditional procedure, where the waveguide dispersion is added to whatever β (ω) is appropriate for the medium's time response. Consequently we need only make the same assumptions: (a) that the two sources of dispersion are weak and so can simply be added, (b) that the temporal response is assumed to have negligible effect in the waveguide's propagating modes, and (c) that the temporal dispersion is the same at all points in the waveguide.
For simplicity we consider the usual situation where our dispersion has already been characterised and represented as a truncated power series in frequency. In line with the recommendations of [25], I exclude imaginary components (e.g. loss or gain terms) in either the wavevector reference β or the frequency reference Ω. Consequently, the reference wavevector considered here has the simple form where κ 1 = c −1 +κ 1 has two contributions -that due to the vacuum (just 1/c), and the modification due to the temporal response (κ 1 ). I have chosen this form for β (ω) because it is possible to calculate its inverse, although in general -if there are more terms in the expansion, or β has some specific and non-trivial analytic form -the problem is more complicated. For this we again replace ω with Ω to denote the new role as a reference frequency, and replace β with k z . With Ω n = ck zn , this gives which we could expand into a power series for small k zn − k z . Since κ 1 includes the vacuum, and dispersion κ 2 is assumed to be small, we should make the top sign choice since it gives the physically relevant small frequency corrections. The linear term in the series has the coeffient γ 1T = 1/κ 1 and the quadratic term γ 2T = κ 2 /κ 3 1 . We now combine the two dispersive effects -the spatial/geometric/waveguide ones based on parameters γ i , and converted temporal ones based on κ i , to give a total reference D2OWE Time propagated optics and dispersion Dr.Paul.Kinsler@physics.org http://www.kinsler.org/physics/ frequency based on k z of (5.9) Note that to avoid incorporating the vacuum contribution twice when summing the dispersive effects, I have replaced γ 1 withγ 1 .
In the complementary case, useful for spatial propagation, we can follow the converse -and indeed the usual procedure -to convert the CW effect of the waveguide's spatial structure Ω(k z ) into a frequency domain form [38,39]. This then can be combined with temporal response of β (ω) to get (5.10)

Discussion
We might have hoped to find a less approximate comparison between the spatially propagated and the temporally propagated approaches to dispersion what has just been discussed. Further, although qualitative comparisons can also be made, ideally we would like a specific and quantitative comparison as done for e.g. testing the unidirectional approximation [31].
That this is not as simple as it might sound can be demonstrated as follows: write two mathematically identical but physically and notationally distinct propagation equations, i.e.
∂ z E(ω) = +ıK(ω)E(ω) (5.11) ∂ t D(z) = +ıΩ(z)D(z). (5.12) For simplicity, this comparison is done in the 1+1D limit and for unidirectional equations; the residual ω or z behaviour has been merged with the reference behaviour to give single propagation parameters K(ω) or Ω(z). Superficially, the comparison between eqn. (5.11) and (5.12) looks promising. However, the vacuum part of K(ω) is ω/c, but the vaccum part of Ω(z) is in contrast a spatial derivative, i.e. c∂ z . Even if it were possible to match up the functional forms of the non-vacuum temporal response of K(ω) − ω/c with a spatial structure giving Ω(z) − c∂ z , the vacuum parts would remain unmatched. Alternatively, we could compare the z-propagated ω-domain to the t-propagated k-domain. However, although in this case the reference behaviours match, but the k-domain model now has a convolution that does not appear in the ω-domain one.
As a result, it is unclear how well approximating temporal dispersion as spatial properties works in practise; but then it is equally unclear how well approximating spatial properties as temporal dispersion works -even if this latter process is so ubiquitous in optics as to be effectively invisible. Removing either of the mismatches discussed above -either in reference behaviour or convolution -requires a narrowband (CW) limit, a process which necessarily obscures the parameter regime where interesting tests and comparisons can (and should) be made -i.e. the short pulse, large bandwidth limit. Indeed, the stringent nature of the CW approximation would defeat the entire point of the comparison. Nevertheless, it may be possible to design numerical simulations containing structures and temporal responses whose scales are carefully graduated so that the mundane bandwidth issues can be kept separate from the interesting propagation ones. Such a design would, as desired, allow the propagation approximations (only) to be compared and contrasted by numerical simulation. However, this is a non-trival task, and one which I leave for later work.

VI. PULSE PROPAGATION
Third order nonlinearities are common in many materials, e.g. in the silica used to make optical fibres [18,38,39]. Here I compare 1+1D propagation such a waveguide, in order to elucidate the differences between the standard spatially propagated picture and the temporally propagated picture.
The assumptions made are those of transverse fields, weak dispersive corrections, and weakly nonlinear response; these all allow us to decouple the forward and backward wave equations. This decoupling then allows us, without any extra approximation, to reduce our description to one of forward only pulse propagation. The specific example chosen here is for an instantaneous cubic nonlinearity, but it is easily generalized to non-instantaneous cases or even other scalar nonlinearities.

A. Spatially propagated NLSE
The spatially propagated NLS for E + x (ω) based on a directional factorization [14] can be written where F[...] is the Fourier transform that converts the nonlinear polarization term into its frequency domain form. Here we have chosen a fixed, non-dispersive reference wavevector β , the temporal response of the propagation medium is encoded in the coefficients κ j . For the temporal (time response) of the material, this is an in-principle exact re-representation of that response as a Taylor series expansion; although in practise the series is usually truncated after a few terms and the propagating fields restricted to within a bandwidth where only those terms are significant.
The important point to remember is that in this spatially propagated picture we do know the full time history, at least D2OWE Time propagated optics and dispersion Dr.Paul.Kinsler@physics.org http://www.kinsler.org/physics/ in a computational sense -i.e. as we compute a solution to the propagation equation. This "computational past" [16], comprising all previously computed values of the field properties and material state means that a full frequency spectrum and response is known. As a result, the time response of the material can be directly encoded in an accurate and useful dispersion relation, and then be approximated as a Taylor series in ω.
Since this propagation is assumed to be taking place in a waveguide, we also have its geometric -wavelength sensitive -properties to consider. However, as discussed in the previous section we typically just re-represent those as if they were instead generated by a temporal response (see e.g. [38,39]) Thus the coefficients κ j usually combine both the effect of time response and geometric properties into a single total dispersion; this is the very basis of dispersion compensation in optical fibres, since the material (temporal) dispersion is offset by the waveguide design (spatial dispersion) [38].

B. Temporally propagated NLSE
Using the approach derived and presented in this paper, we can develop a time propagated counterpart to eqn. (6.1). In a purely mathematical sense, it is, apart from notation, identical to eqn. (6.1). However, in physical terms, the change is not trivial, since the physical meaning and boundary conditions differ significantly. Starting from eqn. (4.1) we can get a propagation equation in wavevector space for D + x (k) which reads where G[...] is the Fourier transform that converts the spatially localized nonlinear polarization into its wavevector domain form. Here, all dispersive behaviour has been combined into the coefficients γ j , which mimic the common notion of spatial dispersion. Such terms can arise directly from a Taylor expansion in k of a relevant CW dispersion relation about some suitable reference wavevector; however, although we can use the expansion to represent the spatial structure of the propagation medium (e.g. a waveguide), the conversion from spatial structure to k expansion is only relatively straightforward in the CW limit.
If the propagation medium also has a temporal response (i.e. temporal dispersion), and if we do not want to model it explicitly, then we can use the approach described in the previous section, where the temporal response of the material was (approximately) represented as if it were instead due to spatial properties.
Finally, note that in ordinary temporally propagated simulations, the computational past and the physical (temporal) past are identical; this is not the case for spatially propagated simulations, as mentioned above.

VII. CONCLUSION
I have derived a general first order wave equation for unidirectional pulse propagation in time that allows for arbitrary dielectric polarization, diffraction, and free electric charge and currents. This propagation equation has a utility all of its own, allowing efficient unidirectional propagation in a strictly causal manner (i.e. time directed). However, it requires comparable approximations to those used in deriving spatially propagated unidirectional propagation equations, although the specific tradeoffs are different.
The propagation equation was derived by first factorizing the second order wave equation for D into an exact bidirectional model, it applies the same type of well defined type of approximation to all non-trivial effects (e.g. nonlinearity, diffraction), and so reduces the propagation equations to a first order uni-directional wave equation. One nice feature of this factorization method is that it provides for a termto-term comparison of the exact bi-directional theory with its approximate uni-directional counterpart.
In addition to its use where the causal simulation of a propagation problem is of strict concern, it also illuminates the rather shadowy process of handling material response and spatial structure in waveguides by constructing a single "combined dispersion" parameter. As discussed in section V, we can now see that this simple addition of temporal and spatial dispersions in an ad hoc manner stands on rather poor foundations, particularly in the few-cycle pulse limit. Further, it is hard to see any straightforward way of testing the limitations of this common procedure, which poses an interesting challenge for the future. That said, in most applications, we expect that a combined dispersion will remain a perfectly adequate approximation, since the pulses are rarely so short, and it seems reasonable to assume that the character of the pulse evolution should remain similar, even if details might differ. Lastly, the lessons learnt from this investigation could also be applied to other types of wave propagation, most notably acoustic systems (see e.g. [27]).