{'fraction_non_alphanumeric': 0.03512205350434221, 'fraction_numerical': 0.015411766761606877, 'mean_word_length': 4.2922137969299055, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the numerical accuracy of the method of multiple scales for nonlinear dispersive wave equations

In this paper we study dispersive wave equation using the method of multiple scales (MMS) and perform several numerical tests to investigate its accuracy. The key feature of our MMS solution is the linearity of the amplitude equation and the complex nature of the time-frequency. The MMS is tested as an initial value problem using three choices of the dispersion model, one toy and two Lorentz models. Depending on the parameters of the problem, the amplitude equation can be both well- or ill-posed. Despite the ill-posedness, the MMS solution remains a valid approximation of the solution to the original nonlinear model.


Introduction
For mathematical models of waves propagating in material media, the phenomenon of dispersion is frequently caused by the requirement of causality.
For example, for the case of a light pulse propagating through an isotropic and homogeneous media, the propagating pulse induces a local dipole density, P, which for the simplest cases takes the form meaning that the polarization at a time t only depends on the electric field at times previous to t. This memory effect, which is the embodiment of causality, is in optics called temporal dispersion, or just dispersion. In this paper, in order to make our discussion specific, we will focus on this particular context in our work, but our methods and results apply quite widely to dispersive wave propagation. The presence of dispersion evidently spells trouble for the integration of the governing equations for the waves. In general, they cannot be solved as an initial value problem.
This phenomenon is of course well known in the theory of wave propagation, and in optics in particular, and various more or less ingenious methods has been invented for getting around the problem.
In optics one frequently tries to get around this problem by solving Maxwell's equations for optical pulse propagation as a boundary value problem, rather as an initial value problem.
In fact, one could argue that the boundary value problem is more closely aligned with the way experiments are done, than the initial value problem is. Waves in a material slab is launched by shining a source laser at the interface of the slab, and therefore one could say that the incoming field at the boundary of the slab is fully controlled by using the laser. Thus, we have the data necessary for solving a boundary value problem for optical pulse propagation. This is the basic, and in fact the only premise, which underlie the Unidirectional Pulse Propagation Equation(UPPE) [1] [2] [3] approach to optical pulse propagation. Of course, if there is significant back scatter of light from the interface and/or the material, the boundary condition is not fully controlled, and UPPE, and also other, less general, boundary solvers, which all rely on unidirectionality, are in trouble. For such cases, one can still solve Maxwell's equations for optical pulse propagation as a boundary value problem by using a more general approach than UPPE. This approach is called the Bidirectional Pulse Propagation Equation(BPPE) [4] [5].
The BPPE approach is an exact method in the sense that no solutions has been lost when transitioning from Maxwell to BPPE. It is however also a purely numerical method and does not give any analytical insight into the pulse propagation problem. The UPPE is also an exact method, but only if one restricts to solutions of Maxwell than satisfy the condition of unidirectionality, and it is also a purely numerical method.
A much older approach to optical pulse propagation solves Maxwell's equations as an initial value problem by restricting to a class of solution that are spectrally narrow, also called narrow band solutions. This restriction makes it possible to derive equations, in general called amplitude equations, that, for a limited time, give a good approximation to Maxwell's equations, for solutions that are spectrally narrow. These amplitude equations can be solved as initial value problems. The systematic approach for deriving these equations is the method of multiple scales(MMS) [6] [7]. In optics, the best known such amplitude equation is the nonlinear Schrødinger equation(NLSE). Another well known amplitude equation, which is particularly useful for driven optical systems like a Laser, is the Complex Ginsburg-Landau equation [8].
From a purely numerical point of view, the great thing about amplitude equations is that they are much faster to solve numerically than the original Maxwell's equations. The root cause of this is that for the narrow spectrum solutions represented by the amplitude equations, the fast frequency at the center of the spectrum, the so-called carrier wave, needs not be temporally resolved, only deviations from the center frequency needs to be resolved, and the range of these deviations is by assumption small.
In addition to being fast to solve numerically as an initial value problem, amplitude equations tend to have a universal form, at least to leading order, and this universal form, quite frequently, makes the equation amenable to analytical investigations. In the best cases, a complete analytical solution can be found. This is true for the NLSE equation [9]. However, this analytic solvability is not robust. If we want to extend the amplitude equation beyond the leading order, which we must, if we want an equation that approximate the narrow spectrum solutions to Maxwell for a longer interval of time, the analytic solvability is typically lost, but the fast numerical solvability is not. So, one can say that the important feature of amplitude equations in optics, is that they form the basis for a fast numerical approach for obtaining certain types of solutions to Maxwell's equations.
The aim of this paper is two-fold. Firstly, we want to go beyond the use of stationary modes as a basis for the MMS expansion, and consider the more general case of modes that are decaying, usually because of loss of some sort in the material, or gain for that matter. In this paper we will focus on the situation of lossy materials. The case of weak losses can treated by using stationary modes, assuming that the decaying terms are small compared to the leading part of Maxwell's equations. For situations where the loss is too large to be included as a perturbation, one must use a MMS based on decaying modes. This is a situation that occur if one is investigating optical pulse propagation close to a material resonance. It is also the case if one want to derive amplitude equations in near-zero index situations. In fact, this last case is the major motivation for the work we do in this paper. Two interesting features of the resulting amplitude equations, features that we explore in detail later in this paper, is that the amplitude equations we derive are necessarily linear and that they are frequently ill posed . The interesting thing is that they, nevertheless, accurately represents both the linear and the nonlinear dynamics of the narrow band solutions of Maxwell's equations, for time intervals for which the amplitude equations according to the MMS procedure should approximate such solutions well.
Secondly, we want to investigate the validity of the amplitude equations, as a fast numerical scheme for narrow band solutions to Maxwell's equations, by comparing the numerical solution of the amplitude equation to the corresponding numerical solution of Maxwell's equations. Since we cannot solve Maxwell's equation numerically as an initial value problem, this validity check has an obvious problem that needs to be handled. Handling this problem is the second major focus of this paper. It relies on the fact that in optics one almost always approximate the linear dispersion of materials using Sellmeier equations [10].

A model wave equation, including nonlinearity and general temporal dispersion
The basic model equation we will use to illustrate our methods, and for which we will state our main conclusions, is the simplest nontrivial wave equation from nonlinear optics. It is scalar, and includes general linear dispersion and a Kerr, cubic nonlinear material response. Our methods and conclusions apply much more widely than this, but for the sake of clarity and because numerical methods plays an important role in this paper, it is necessary to work within a specific class of equations.
In order to put our model equation into a real physical context, we start this section by deriving the equation from Maxwell's equations under some reasonable physical assumptions on the material response.
Maxwell's equations for a situation where there are no free charges or currents, are given by Most materials show no magnetic response at optical frequencies, thus we assume that The polarization is in general a sum of a term that is linear in E and ones that are nonlinear in E. We thus have We assume that the linear material response is isotropic, homogeneous and causal In the explicit calculations that we do in this paper, we will assume that the nonlinear polarization is restricted to the Kerr effect. Thus we will assume that P N L = ε 0 ηE · EE, (2.5) where η is the Kerr coefficient. This is a choice we make just to be specific, the applicability of our methods, and the validity of our conclusions, derived in this paper, in no way depend on this particular choice for the nonlinear response. Inserting (2.2)-(2.5) into (2.1), we can rewrite Maxwell's equations into the form where we have used an alternative form of (2.4) which is derived in Appendix A. The factor √ 2π is a consequence of our conventions for the Fourier transform.
We will now restrict ourselves to solutions of the form where for the assumed Kerr effect we have By taking cross derivatives it is easy to eliminate the magnetic field components and arrive at the equation which is the basic model equation we will be using in the rest of this paper.

Scaling of the model equation
If the aim is to solve Maxwell's equations numerically, in some specific physical context, it is not really necessary to scale the equation, and frequently this is not done, even if one can argue that it could still be a useful thing to do. However, if one is going to derive an amplitude equation for the same physical situation, it might not be essential to scale the equation, but it certainly is extremely useful. After all, the essence of MMS is the ordering of terms in certain expansions according to size, and ensuring that this ordering, according to size, persists, up to some time of our choosing. Equation (2.11) will be the staring point for our multiple scale approach. Let us start by picking scales Z 0 , T 0 and E 0 , for space, time and electric field amplitude, so that we have where the primed symbols are the scaled quantities. With these choices of scales the model equation (2.11) takes the form We are at this point free to choose the time scale T 0 , so let us choose it so that the factor before the z-derivative becomes one We next set the scale for the electric field to be the initial peak electric field amplitude. Thus With this we can write the model equation in the form where we now have dropped primes on all quantities, since from this point on, only scaled quantities will appear. In this equation we have introduced the dimensionless parameter ε = √ ηE 0 . Typically, the Kerr parameter is fixed for any given material, whereas E 0 is at our disposal to vary over many order of magnitude, depending on the strength of the laser used to generate the initial electric field. This means that ε can be made to vary over many orders of magnitude, but for realistic field intensities it is always smaller than one, usually much smaller than one. ε is the small perturbation parameter which we need for the MMS expansion.
In this paper, Fourier transforms and plane waves play a prominent role, and whenever that is the case, it is convenient to pick the units for wave number and frequency in such a way that the phase of plane waves, and the Fourier transform, and its inverse, retain the same symbolic form in scaled and unscaled quantities. It is easy to verify that this is the case if we make the choice ) Thus we measure frequency in the well known unit Hertz, or cycles per unit time. Wavelength is in a similar way measured in periods per unit length.
Note that (2.14), which fixes the time scale in terms of the length scale, also, because of (2.18), fix the frequency scale in terms of the wave number scale Thus, the only scale that remains to pick is the one for wave number. The initial field is in a lab situation generated using a laser. These days most optical labs have lasers that produce pulses of femto second duration [11], and labs with more specialized equipment can produce pulses down to atto second durations [12]. Such short pulses have a very broad spectrum and their dynamics are therefore hard to approximate using amplitude equations, which require narrow band pulses. There are versions of MMS that can handle such broad band pulses, but only in the weak dispersion limit.
In this paper we apply MMS in a form that is tailored to the opposite limit of strong dispersion. Thus, in this paper we must assume that the initial field has a narrow wave number spectrum centred on a wave length determined by the lasing wave length of the laser generating the initial pulse. What the word "narrow" in the previous sentence means will be clarified in the MMS expansion in the next section.
In this paper we are focused on validating amplitude equations derived using decaying modes in the vicinity of a material resonance, and it thus makes sense to pick the scale, K 0 , for wave number, in such a way that the position of the resonance is centred on scales wavelengths that are of order one. The initial laser pulse will then have a narrow band spectrum centred on a scaled wave length that is also of order one. We will in the rest of the paper assume that this has been done, and from now on only refer to scaled quantities, both in the model equation and in the specific material models that we will introduce in the sections to come.

Decaying mode, amplitude equation, for the model equation
We proceed with the multiple scale method by introducing the expansions ∂ t = ∂ t 0 + ε∂ t 1 + ε 2 ∂ t 2 + . . . , ∂ z = ∂ z 0 + ε∂ z 1 + ε 2 ∂ z 2 + . . . , e = e 0 + εe 1 + ε 2 e 2 + . . . , (2.20) where the connection between the multiple scale field amplitude e and the electric field amplitude E is given by The functionχ(i∂ t ) is expressed as a Taylor series as follows. (2.25) For the order ε 0 equation we choose the wave packet solution and where ω = ω(k) is a complex solution to the dispersion relation Here, the complex refractive index, n(ω), is defined by Observe that our multiple scale expansion is based on a complex, decaying mode, not a complex, stationary mode, which is usual when one apply MMS far from any resonances of the material. We will see here and also in later sections that this fact will change many aspects of the resulting amplitude equations. We must now calculate the right-hand side of the order ε 1 equation. Inserting (2.26) into (2.25) and (2.26), we get In order to remove secular terms we must postulate that Observe that from the dispersion relation (2.28) we have (2.32) Thus (2.31) can be written in the form For the case of stationary modes this quantity is real and by definition equal to the group velocity for an initial light pulse with a narrow spectrum centred on the wave number k.
The order ε 1 equation simplifies into We now must compute the right-hand side of the order ε 2 equation. Inserting (2.35) into the right-hand side of the order ε 2 equation we get where ω i = Im ω. At this point it is worth observing that none of the nonlinear terms in (2.36) are secular. In addition to the usual nonsecular term e i3θ 0 , we also have the term e iθ 0 e 2t 0 ω i , which would have been included into the secular terms, and thus, in the final amplitude equation, if it wasn't for the complex nature of ω. This makes MMS based on decaying modes essentially different from the case of stationary modes.
Continuing the calculation, we observe that in order to remove secular terms from (2.36), we must postulate that The factor multiplying the term ∂ t 2 A 0 can be expressed using (2.32), and equation (2.37) therefore simplifies as follows By removing the secular terms from equation (2.36), the order ε 2 equation turns into This equation is not homogeneous, and we chose at this point to solve for e 2 , using only a particular solution. One such particular solution is evidently where Defining an amplitude A(z, t) by and proceeding in the usual way, using (2.33) and (2.38) we finally get the following amplitude equation The amplitude A is related to the electric field amplitude E through the formula By design, for (2.46), (2.47) to be an approximate solution to (2.16), we must have where we recall that ε is a number much smaller than 1. Given these circumstances, we observe that and thus, the, second order in time, amplitude equation (2.46), is asymptotically equivalent to the more convenient, first order in time, equation Note that if we want to have an amplitude equation which is a good approximation to the model equation beyond a time of order 1 ε 2 , we have to extend the MMS procedure to higher order in ε. It is evident from what we have said about secular terms in the paragraph following equation (2.1), that these extended amplitude equations will all be linear, no matter to which order the MMS procedure is extended.
Since the amplitude equation is linear, the solution space, S A , is of course a linear space. The relation (2.47) amounts to a map, M , from S A into the solution space, S E , of the exact equation (2.16). The map is certainly not surjective, and neither is it injective. The lack of injectivity means that the map cannot be used to induce a nonlinear superposition principle on its image set, M (S E ) ⊂ S E , using the usual pullback/pushforward approach.
In order to get the approximate solution to our model equation defined by the amplitude equation (2.50), and the relation (2.47), we decided to pick particular solutions at order ε and ε 2 , at both orders disregarding the general solution to the homogeneous equation. The consequence of adding a solution to the homogeneous equation, in the form of wave packets, at one or both orders, would be to add one or two new independent amplitudes to the problem. Each of these amplitudes would, through the removal of secular terms at order ε, ε 2 and ε 3 , satisfy their own amplitude equations. Both these extra amplitude equations would also be linear, and uncoupled from each other and the one for the amplitude A. The relation defining the electric field in terms of the amplitudes would now be much more complicated and involve sums of products of all three amplitudes.
The deciding factor for whether we include these extra amplitudes, or not, is what kind of solutions of the model equation we are trying to approximate. This comes down to which kind of initial conditions for the model equation we are able to represent faithfully using our amplitude equation.
Any choice of a narrow-band initial amplitude for the amplitude equation (2.50), will lead to a narrow band wave packet solution of the model equation, with dispersive properties determined by the choice of a complex solution ω = ω(k) to the dispersion relation (2.28), at order ε 0 of our MMS expansion. This is certainly a valid choice of initial condition, if our aim is to validate the amplitude equation (2.50) together with its corresponding defining relation (2.47) for E.
However, such a choice of initial condition for E is awkward from a physical point of view. Usually the initial field is determined by a laser, whose output, in most cases, can be approximated spectrally by a narrow Gaussian centred at the operating wave length of the laser. The initial condition for E, defined above using a narrow-band initial condition for A, consists of two Gaussians centred at k and 3k, and with a special relationship between the amplitudes and phases of the two Gaussians.
From the relation (2.47), it is evident that in order for E to be a Gaussian centred on some wave number k, the amplitude A must be a Gaussian centred on k = 0, and we have to introduce a solution to the homogeneous solution at order ε 2 of the form leading, in the end, to a relation determining the electric field in terms of the amplitudes A and By fixing the initial value of the amplitude B to be we ensure that the initial spectrum for E is a Gaussian, whose center and width, is determined by the initial spectrum for the amplitude A. The initial value for the amplitude A is found by solving for A using for example Newton's method. Since our aim is to get an approximation to the solution of the model equation to order ε 2 for times t ≤ ε −2 , we need to remove the secular terms generated by the amplitude B, at order ε 3 , since such terms would grow linearly and potentially disturb the spectrum at 3k, already for times of order ε −1 .
From detailed calculations done while deriving the amplitude equation (2.50), it is not hard to see that the resulting amplitude equation for the amplitude B must take the form For the validity tests we will discussing later in the paper, we will always introduce such an extra amplitude B, whenever it is necessary for faithfully representing a Gaussian initial spectrum for E.

Testing the validity of the amplitude equations
The first goal of this paper was to use the MMS approach to derive the amplitude equation close to a material resonance for nonlinear wave equations with arbitrary, but of course causal, material response. Our main model equation and motivation comes from the field of nonlinear optics, but our methods and results can evidently be applied much wider than this. Our second goal, which is the focus in this section and the rest of the paper, is to test the validity of the derived amplitude equations using high precision numerical simulations. If one is far from a material resonance, the dispersion relation of the material can to a good approximation be modelled by a polynomial function of the frequency, usually a second order polynomial is sufficient like in the Lorentz model [13]. In this situation, a direct simulation of the model equation can fairly easily be achieved, and the validity of the amplitude equation can thus be tested. However, even in this situation, numerical validation is not all that common. One reason for this is that the narrow band solutions, which are well approximated by the amplitude equations, are also the type of solutions to Maxwell that are most challenging to simulate numerically. This is because such solutions, which basically are wave packets, have a slowly varying, and consequently very wide, envelope, and at the same time contain a very large number of oscillation under the envelope. Thus one need both a large computational domain and a very high resolution of that domain. This makes for a large number discretization points and thus long running times.
When one is close to a material resonance, low order polynomial approximations does not work as well as when one is far from a resonance. In this situation different and more complicated approximations must be used. This typically turns the model equation, which is a differential equation in time, into an equation that is a pseudo-differential equation in time. Solving such a thing as an initial value problem is not an easy proposition.
A frequently used class of such, more general approximations, are the rational functions. This type of approximations is in particular very much used for approximating the electric susceptibility in optics. For this situation, one can exactly transform the model equation into a equation that is a differential equation in time, and whose initial value problem can be solved numerically. The transformation is based on what is called the Sellmeier formulas in optics, and it thus make sense to denote the associated transformation for the Sellmeier transformation.

The Sellmeier Transformation
The Sellmeier formulas in optics are approximations of the electric susceptibility in terms of sums rational functions of the simple type.
where a, b and c are complex constants. Causality, and questions of loss or gain in the material, put restrictions on the constants that we leave aside for now. Later, when we do our numerical simulation in order to validate the amplitude equations, these restrictions come to the fore. Finite sums of functions of the type (3.1.1) will produce general rational functions. Thus the Sellmeier formulas are simply rational approximations to the electric susceptibilitŷ In order to describe the Sellmeier transformation, observe first that (2.16) can be written in the form L(E, E 3 ) = 0, where L is a suitable operator that produces our equation (2.16) and which includes the integral operatorχ(i∂ t ). Secondly, observe that the Fourier transform of our model equation can be rewritten into the formL Assuming the form of the susceptibility to be described by the very simplest Sellmeier formula (3.1.1),χ(ω) = R(ω), the operatorL is the differential operator It is evident from (3.1.3) that any solution to the differential equatioñ is also a solution to the original pseudo differential equation ( Our idea is now to restrict the numerical validation to solutions in the smaller solution space S D . This of course only make sense if the MMS procedure applied to the Sellmeier transformed equation (3.1.5) produce the exact same amplitude equation as the one we got earlier, starting from the original pseudo differential equation (2.16). From our short discussion of the relation between S M and S D in the previous paragraph, we certainly expect to get the same amplitude equation, but we still feel that it is prudent to directly verify that this is the case for each of the two explicit examples of Sellmeier formulas discussed in the following sections.

A toy model for dispersion
In this first validation calculation for our amplitude equation (2.50) and associated reconstructed electric field amplitude (2.47) we chose a material response function of the simple form for some real positive constants u and v. The corresponding electric susceptibility, which is the Fourier transform of (3.2.1), is given by the formulâ Defining parameters γ = v u , a = 1 u , the formula for the susceptibility can be written in the more convenient formχ We are not aware of any material response that in the optical regime is described by this electric susceptibility, but it does satisfy the, all important Kramer-Kronig relations, and thus describe a causal response. For us, this choice is a device for testing the amplitude equation in the simplest setting possible. Later in the paper we will investigate the validity of the amplitude equation when the electric susceptibility describe a very common atomic model, the Lorentz oscillator. For this case the validation is conceptually the same as for the simple model above, but it is technically much harder.

The Sellmeier transformation
The model for electric susceptibility described above, leads to a Sellmeier transformed equation that, in our scheme of things, is as simple as possible to handle numerically. Following the procedure described in the previous section, we find that the Sellmeier transform of our model equation is This equation is, as expected, a differential equation and not a pseudo differential equation. It can be solved numerically using, for example a pseudo spectral method, which is what we do in this paper. As noted earlier, for the validation test to make sense, we must ensure that the MMS procedure applied to (3.2.4) gives us the same amplitude equation (2.50), which we got from the original model equation (2.16). We do this is Appendix B and observe there that the derived amplitude equation is indeed the same as the original one (2.50).

Numerical results
Actually solving the Sellmeier transformed equation (3.2.4), is awkward, because it is not solved explicitly with respect to the highest derivative. However, the fact that ε is small means that, by iteration, we can approximate (3.2.4) by an equation that is solved explicitly with respect to highest derivative. We achieve this by expanding out the offending term a∂ ttt E 3 as follows and substituting for ∂ ttt E using (3.2.4). Dropping terms of order ε 4 , which we must do in order to be consistent with our order ε 2 MMS expansion, we get the following explicit equation Since the explicit equation (3.2.6) agree with the, exact, implicit equation (3.2.4) to order ε 2 , the amplitude equation to order ε 2 for these two equations must be the same. Thus we can use (3.2.6) to validate our original amplitude equation (2.50).
Recall that the amplitude equation is derived using one specific mode of the linearization. Thus, to do the validation we must select one of these modes. Given this mode, all the parameters in the amplitude equation are determined in terms of the parameters of (3.2.6). The detailed expressions are derived in Appendix B. This linear mode, which is an explicit function of z and t, can now be used to determine all the three initial conditions needed to solve (3.2.6) numerically, and thus to complete the validation.
The equation (3.2.6) has three independent modes. The solution to the linear part of (3.2.6) can be expressed as inverse Fourier integrals of all three modes (3.2.7) The electric field E(z, t) must be real so (3.2.7) must be equal to its complex conjugate Let us first assume that k > 0. From the form of the dispersion relation (B.6) we observe that if ω(k) is a solution to p(k, ω) = 0, then we have This implies that if ω(k) is a solution, then −ω * (k) is a solution as well. Thus, assuming that the three solutions are distinct, we can number them in such a way that Next observe that the dispersion relation (B.6) is even in k. Using this fact, we can number the solutions for negative k is such a way that ω 1,2,3 (k) = ω 1,2,3 (−k). Using these relations, together with (3.2.11) and (3.2.12), formula (3.2.8) turns to For (3.2.7) and (3.2.13) to be the same, the following relations between the amplitudes A 1,2,3 (k) must hold (3.2.14) From these relations, we see that the amplitudes for k < 0 are determined from their values for k > 0 and vice versa. While deriving the asymptotic solution in Appendix B, we assumed the solution to the 0-th order equation in the perturbation hierarchy (B.1) to be one of the three possible modes. The amplitude equation (2.50) is the correct equation only for narrow band solutions consisting of one such mode, centred around wave number k 0 which we, without loss of generality, can assume is positive. To be specific, assume that this mode is A 1 (k). The narrow band property implies that A 1 (k) = 0 for k < 0, and from this the relations (3.2.14) implies that we can consistently choose A 3 (k) = 0 for all k and A 2 (k) = 0 for k > 0.
Using these assumptions and (3.2.14), we get from (3.2.13) where we let ω 1 (k) = ω(k) to match the parameters in our asymptotic solution (2.47) with (3.2.15). We now turn back to equation (3.2.6). Upon transforming this PDE into its spectral domain k, we get a third order ODE's with a nonlinear right-hand side. We want to solve this ODE as an initial value problem, and therefore need three initial conditions, here denoted by f (k), g(k) and h(k). These three initial conditions we obtain from (3.2.15) f (k) =Ê(k, 0), These conditions can be also expressed in terms of the amplitude A 1 (k) using (3.2.15)

2.19)
The only thing left to do now, is to compute the initial condition for the amplitude equation (2.50) in terms of our chosen amplitude A 1 (k).
As discussed at the end of section 2.3, the initial condition for E, which is natural from a physical point of view, is one whose spectrum is a narrow Gaussian. This is taken care of by letting the spectral amplitude A 1 (k) be a Gaussian centered at some wave number k 0 where D, δ > 0.
As also discussed at the end of section 2.3, we might have to introduce an extra amplitude at order ε 2 in order to faithfully represent the initial condition on E in terms of an initial condition for the amplitude equation.
The way to determine if any extra amplitude has to be introduced at order ε 2 , is to assume the opposite. Given this, the relation between the amplitude A and the electric field E is determined by (2.47).
We now take the inverse Fourier transform of equation (3.2.17) and equate it to the right-hand side of (2.47) evaluated at t = 0. Matching separately the first part and the second part, which is the complex conjugates of the first part, on both sides, we get which is a nonlinear system of algebraic equations, consisting of two equations and two unknowns A(z, 0) and A * (z, 0). This system can easily be solved numerically, for example using Newtons method.
Note that the parameter δ controls the width of A 1 (k). Given that the amplitude A(z, t) should be slowly varying in z, ∂ z A ∼ O(ε), the parameter δ should be chosen such that δ ∼ 1/ε.
In order to do a numerical comparison, we fix the refractive index by choosing the parameter values γ = 5 and a = 20. In figure 1 we see the real and imaginary part of the resulting refractive index of our material.
At this point, the only remaining parameter to set is the, all important, nonlinearity parameter ε. In all our derivation, we assumed that this parameter was small, only then will there be separation of scales, which is the key assumption underlying our derivations of amplitude equations using MMS. Usually, one is deriving amplitude equations in order to gain some analytical insight into the wave propagation problem at hand, and the actual numerical value of ε does not need to be fixed.
Here however, where we are doing a numerical validation of the amplitude equations, the situation is different, here we need to give a specific value to the nonlinearity parameter. And, with respect to this, we face a tradeoff.
On the one hand, we should pick a value for ε that is as small as possible, in order for the separation of scales to be as large as possible. Only for such values can we expect our amplitude equations to accurately approximate the original model equation.
On the other hand, a very small value for the nonlinearity parameter, means that we have to propagate the waves for very long distances for the nonlinearity in the model equation to influence the spectrum of the waves. Very long distances means very long times which translates into very long running time for the simulations of the model equation.
In this paper we choose the value for the nonlinearity parameter to be ε = 10 −1 . This certainly does not seem to be a very small number, but as we will see, even for a nonlinearity parameter as large as this, our amplitude equation does actually approximate the exact model equation well. This is just another example of what one could call the phenomenon of the unreasonable accuracy of asymptotic methods.
As it turns out, our chosen value for the nonlinearity parameter is also physically reasonable. For example, for a light pulse in a visible part of the spectrum, at 586 nm, propagating through Argon gas at atmospheric pressure, of an intensity equal to half the critical ionization threshold, the value for the nonlinearity parameter is ε = 0.13 [14].
Let us now point our attention to calculating the correct initial condition for A. In order for the amplitude equation (2.50) to be valid, the amplitude A must be spectrally narrow. In fact, the whole MMS expansion is based on the assumption that the width of the spectrum of A is of order epsilon.
The spectrum for A which we find assuming that no new amplitude is required at order ε 2 is displayed in figure 2. The parameters used to calculate the initial condition for A were a = 20 and γ = 5. Evidently, the width of this spectrum is not of order ε. The main part of the spectrum is narrow, but there is an additional peak in the spectrum centred at 2k 0 whose distance from the origin is of order one, not order ε. Thus, assuming that no new amplitude is needed at order ε 2 leads to a contradiction.
From relation (2.47) we observe that the peak in the spectrum for A corresponds to a peak at 3k 0 for E. Furthermore, from figure 2 we observe that the height of the peak is of order ε 2 . From these two observations it is evident that the peak in the spectrum for A can be taken into account by introducing an extra amplitude B in the MMS expansion at order ε 2 . As we outlined at the end of section 2.3, the amplitude B will come equipped with its own linear amplitude equation, decoupled from the one for A.
The initial values for A and B are now found by solving equation (2.54) for A, using the approach from (3.2.21), and then using the identity (2.53) to determine the initial value for B in terms of the one for A. In figure 3 we display an example of the initial condition A(k, 0). In this example we use the value k 0 = 2π. Note that this means that our initial pulse spectrum is located to the right of the frequency defining the material resonance, in a region of anomalous dispersion, which means that the real part of the refractive index decrease for increasing frequency. Now we have everything we need in order to compare the numerical solution of the model equation   For the values of the parameters used in this comparison, the complex frequency of the chosen mode is given by ω 0 = ω(k 0 ) = 6.28 − i2.5 × 10 −2 .
On order to get more insight into the accuracy, let us compare the wave number spectrum of these two solutions. This is what is displayed in figure figure 5, which includes both the major bumps in the k spectrum. The larger one has its support at the frequency k = k 0 = 2π as expected. This represents the linear part of the solution. The smaller one is located around the frequency k = 3k 0 = 6π. It makes sense that it sits precisely at 3k 0 , since the Kerr nonlinearity has the form E 3 , which means that when inserting a plane wave for E into E 3 , the wave numbers get multiplied by a factor of 3.
From figure 5b it does appear that there is a small deviation between the exact numerical solution and the solution derived from the amplitude equation.This deviation is however one or more orders of magnitude smaller than ε 2 , and such deviations are to be expected for for a MMS expansion truncated at order ε 2 .
Based on the experience from testing the MMS solution, it became clear that in order to stay in the correct asymptotic regime, where the MMS solution is valid, the parameter values fixing the problem are subject to certain constraints.
First of all, the constants α and β occurring in the amplitude equation, defined in (2.39), (2.40), include the derivatives of the susceptibilityχ(ω). These derivatives come from the Taylor expansion of this function. We expect this Taylor series to converge, so by assumption, these derivatives must not break the order of the preceding terms in the expansion. This assumption depends mainly on the  Figure 5: The spectrum of the solutions of (3.2.6) using a system of first order ODE's (blue graph) compared with the spectrum of the solution (2.47) using the amplitude A(z, t) computed from (2.50) (red graph) for the toy model dispersion.
parameters γ and a in (3.2.2) and the frequency ω around which the Taylor series is expanded. To preserve the order of these terms is closely tied to the possibility of making the assumptions in (2.48) and consequently arriving at the amplitude equation. The assumption made in (2.48) are related also to the choice of the parameter δ in the initial condition (3.2.20). The amplitude A(z, t) is by assumption slowly varying in z compared with the exponential e ik 0 z because of (2.45). Therefore the initial condition (3.2.20) should be fast varying in k and the parameter δ needs to be chosen accordingly, for example δ = 1/ε.

Stability and well-posedness
From the numerical solution of the amplitude equation in the previous section, and its ability to accurately approximate the exact, narrow band, solutions for Maxwell, one might think that our amplitude equation is just fine. However, nothing could be further from the truth. In this section we will show that the amplitude equation is in fact ill posed as a PDE.
Since the amplitude equation is linear, this fact can be easily proven by merely calculating the growth curve for the equation. This curve we find by inserting A(z, t) = e λ(k)t e ikz into the equation, cancelling the common factor e λ(k)t e ikz , and extracting the real part of the resulting algebraic equation.
where a 1 and a 2 are real parameters that depend on the complex parameters α, β and ω (k) is the derivative of the dispersion relation (2.28). This function is a parabola which pass through the origin, and is displayed in figure 6, using the parameter values from our numerical test. From figure 6 it is evident that the amplitude equation is ill posed by definition; arbitrary high spatial frequencies will grow exponentially in time, with no upper bound for the growth rate.
Being ill posed is usually, for good reasons, thought of as ringing the death knell for any proposed mathematical model. But, still, in the situation investigated in the previous section, the ill posed amplitude equation is an excellent tool for simulating narrow band solutions to Maxwell's equation. We regard this as another example of the way asymptotic methods makes use of, and gives meaning to, otherwise meaningless mathematical expressions. Anyone that is a user of asymptotic methods can not help noticing this fact. The classical case of this is Eulers example, where the useless power series n (−1) n n! x n , with zero radius of convergence, is, nevertheless, an accurate approximation to a certain exponential integral. Similarly, the expansion for the energy of the anharmonic quantum oscillator is an asymptotic series with zero radius of convergence [15]. In fact, essentially all series used in quantum theory, in both the particle and the field incarnations of the theory, is known to, or believed to, have zero radius of convergence [16] [17]. But, nevertheless, their usefulness is indisputable, their predictions give some of the most accurate correspondences between theory and experiment in all of science.
The reason why the ill posed amplitude equation nevertheless is an accurate numerical model for the narrow band solutions of Maxwell, for which it was designed, is simple. From the form of the stability curve (3.2.22) it is clear that it passes through 0 for every parameter. Combining this with that fact that the initial condition for A will always be centered at k = 0, we can conclude that in the case of ill-posedness, the growth rate around k = 0 will not be large. For any solution that satisfies the assumptions used to derive the amplitude equation, the fastest growing wave number component in the spectrum of the amplitude A(z, t), will, during the time for which the amplitude equation is valid, t ≤ ε −2 , not grow large enough to affect the solution to order ε 2 or greater.

Lorentzian model of dispersion
In the previous section we used the simplest possible rational approximation to the electric susceptibility. As we noted there, this model is, as far as we know, not a realistic physical model, even though it is causal and thus does satisfy the optical Kramer-Kronig relations. In this subsection we consider the next simplest one, which is of the form This dispersion is a realistic physical model of the electric susceptibility. It can be derived from a purely classical model of an atom consisting of a positive charge, representing the nucleus, together with all the electrons, save one [18]. The remaining electron is singled out by being the one that resonantly respond to an imposed oscillatory electric field. Since the single electron is much lighter than the rest of the atom, which after all contains the nucleus, we in effect describing a simple oscillator, which is a harmonic oscillator, unless the field is strong enough to pull the electron too far from the atom. The model is called the Lorentz oscillator, and have been applied to a vast range of materials in the gaseous, liquid and solid phase. The parameters ω r , ω p and γ are interpreted as the resonance frequency of the oscillator, the plasma frequency and the absorptive loss. The factor of 2π is, as noted before, a consequence of our Fourier transform conventions, which are introduced in Appendix A. The Lorentz oscillator model is causal and thus satisfy the Kramer-Kronig relations. We now repeat the calculations from the previous section for the Lorentz oscillator. We start by finding a differential equationL (E, E 3 ) = 0, corresponding to our model pseudo differential equation (2.16), using the Sellmeier transformation. After that, we use MMS to verify that we get the same amplitude equation for the Sellmeier transformed equation as the one we got from the original equation (2.16). We then move on to doing a numerical comparison of the accuracy of the amplitude equation with respect to the Sellmeier transformed equation. As noted earlier in this paper, this amounts to a direct comparison of the numerical accuracy of the amplitude equation with respect to the original model equation (2.16), which is the goal of this paper. In our calculations we rewrite the function (3.3.1) in a for us more convenient form where we have defined a = −1/ω 2 p , b = −γ/ω 2 p and c = ω 2 r /ω 2 p . In our case, the refractive index is defined as n 2 (ω) = 1 + √ 2πχ(ω). In Figure (7) we display an example of such a refractive index, corresponding to the numerical values a = −0.01, b = −7 × 10 −3 and c = 0.16.
In the next subsections we derive the equation (3.1.3) for the Lorentzian model of dispersion and obtain the MMS solution for it. Next we test the MMS solution using two choices of parameters for the Lorentz model.

The Sellmeier transformation
The Lorentzian model of dispersion leads to the following Sellmeier transformed equation Again, the derived equation (3.3.3) does not include any pseudo differential operator. We proceed in this case in the same way as in chapter 3.2.1 and obtain MMS solutions to (3.3.3) and see if they match with the solutions in the asymptotic regime for (2.22). A detailed calculation of this MMS solution is done in Appendix C. The derived amplitude equation is found to be the same as (2.50).

Numerical results
We proceed by testing the results based on the amplitude equation (2.50) and a solution to (3.3.3) for the Lorentz model of dispersion. Two tests will be conducted, each of them with a different choice of parameters a, b, c for the Lorentz model. The Sellmeier transformed equation (3.3.3) is implisit in its highest time derivative, but can be approximated by an equation explicit in its highest time derivative. Following the same procedure as in chapter 3.2.2 we have As before, the amplitude (2.50) equation is an equation for one of the four independent modes equation (3.3.4) has. In order to test the MMS solution we need to choose one mode for which we have the amplitude equation (2.50). In order to arrive at the point where we can chose the mode, we proceed like we did in subsection 3.2.2. The electric field E(z, t) is expressed as the inverse Fourier transform of the sum of all 4 modes and is then equated with its complex conjugate. In order to proceed we need to establish relations between the mode frequencies found by solving the dispersion equation (C.6).
We observe that Thus like in subsection 3.2.2, solution space of p(k, ω) = 0 is the same as for p(k, −ω * ) = 0. We can therefore conclude that if ω(k) is a solution, then −ω * (k) is as well. Assuming that all solutions are distinct, we can enumerate the four the solutions in such a way that the following relations between ω 1,2,3,4 (k) hold: The equation p(k, ω) = 0 is also even in k, and thus we can enumerate the four solutions corresponding to negative k in such a way that ω i (−k) = ω j (k), for some i and j. Using this fact, together the relations (3.3.7), (3.3.8), the reality of the electric field implies the following relations between the amplitudes A 1,2,3,4 (k):

3.9)
The amplitudes for negative argument are defined from their values for k > 0. As before, we want the electric field consisting of only one amplitude, for which we obtained the amplitude equation (2.50), therefore we set A 2 (k) = A 3 (k) = A 4 (k) = 0 for k > 0. This implies A 3,4 = 0, ∀k. We thus arrive at The amplitude A 1 (k) can be chosen arbitrary and it will define the initial condition E(z, 0). The equation (3.3.4) is solved as a 4-th order ODE in the spectral domain k. For this to work, we need four initial conditions, denoted byf j (k) for j = 0, 1, 2, 3. These conditions are obtained from (3.3.11) and can be expressed in terms of the given amplitude A 1 (k) as

3.12)
By choosing A 1 (k), k > 0, all four initial conditions are defined. From the relation (2.47) between E(z, t) and A(z, t), we then calculate the initial condition for the amplitude equation (2.40) in exactly the same way as in (3.2.21). The same procedure also applies in terms of the extra mode B(z, t) whether or not it should be included. As we will see, this extra mode will be included in the following numerical tests. We proceed to the implementation part and choose the initial amplitude A 1 (k) to be the same as for the toy model. A Gaussian centered at the wave number k 0 .
for some D, δ > 0. We will use different values for k 0 in the next two numerical tests.
Lorentz test for ultraviolet resonance The model in chapter 3.1.2 resulted in an ill-posed amplitude equation that turned out to give us a correct solution within the asymptotic regime for our perturbation scheme. In this test we picked the parameters in such a way that it results in a well-posed amplitude equation. In figure 8 we see the real and imaginary part of the refractive index. The choice of parameters in this figure is 3.14) The parameter for the nonlinear term is again chosen to be ε = 10 −1 . When comparing the two solutions from both equations (3.3.4) and the MMS solution in figure 9, we see that they overlap sufficiently at all the presented times. However, comparing the solutions in the spectral domain can reveal the more subtle differences. On figure 10 we see that they overlap quite nicely. The part of the spectrum we are interested in is the smaller gaussian hump centered at 3k 0 = 24, where the nonlinearity manifests. In order to see the details how well these two functions overlap, we can look at the figure 10 for both the gaussian humps. Figure 10b is rather convincing us of the accuracy of the MMS solution. On figure 10a is the main gaussian that represents the linear part of the solution. As we can see, the MMS solution is indeed accurate up to the order of ε 2 for both the toy model and the Lorentz model of dispersion. The parameters used in figures 9 and 10 are k 0 = 8, ω 0 = ω(k 0 ) = 7.9 − i1.99 × 10 −2 . In the next paragraph we run one more numerical test with a different set of parameters.
Lorentz test for infrared resonance Let us now choose the parameters a, b, c for the Lorentz test for infrared resonance. In this model we are using the following parameters in the refracting index:

3.15)
We can see the plot of this refraction index in figure (11). We also notice that the graph of this function is somewhat unusual. The resonance for this function is much closer to 0 than for the index in figure 8. However, the resonance area in the non-scaled version of figure 11 is located at the frequencies that are around 1.5 × 10 14 (with the scaling factor Ω 0 = 1.5 × 10 15 ) which is in the infrared range. The initial condition for the amplitude remains the same as in (3.3.13) as well as the parameters δ, D, η and E 0 . The refractive index parameters give us the scaled frequency of the pulse obtained from the dispersion relation ω 0 = 6.29 − i4.91 × 10 −2 .
As seen from the results in figures (12) and (13) we can conclude that the MMS solution proved itself and it is sufficiently accurate compared to the high precision numerical solution up to the error of order ε 2 = 10 −2 . The parameters used in figures 12 and 13 are k 0 = 2π, ω 0 = 6.29 − i4.91 × 10 −2 .
In the next subsection we investigate the stability of the amplitude equation for both the Lorentz tests for ultraviolet and infrared resonance.

Stability and well-posedness
In this subsection we look at the stability of the amplitude equation (2.50) with the parameters used in Lorentz tests for ultraviolet and infrared resonance. The formula for the stability is the same as in (3.2.22). On figure (14a) we can see the graph of the function Re λ(k) for both tests. As mention earlier in Lorentz test for ultraviolet resonance, we chose the parameters for this test such  We did three numerical tests; the first corresponds to a toy model of dispersion, and the other two the, more physical, Lorentz model of dispersion. The MMS solution performed very well in all three cases. For the two Lorentz cases, we chose the frequency of the initial pulse to be in front of and behind the resonance, where the first case gave us a well-posed amplitude equation while the second one turned out to be an ill-posed one.
The key idea behind MMS is to maintain the ordering in asymptotic expansions. This ordering also applies to the Taylor expansion of the susceptibility (2.22). It is clear, that this ordering depends on the parameters inχ(ω) itself, but also on the frequency ω 0 around which we are expanding. In our countless numerical tests we have tried many different combinations of the parameters. We observed that this ordering failed for certain values for the parameters. Especially when one wants to be near or at resonance. Certainly, further investigation is needed to pinpoint exactly where it can go wrong. The derivatives ofχ(ω) also appear in the constant α, which is required to be of a certain order since it appears in the ε 2 part of the amplitude equation. Alternatively, different methods may be used to expand the susceptibility, for example rational function expansion.
Another problem with our approach that also relates to the values of parameters inχ(ω), is the growth of the nonlinear term in the Maxwell's equation. As expected, if the height of the spectral peak at 3k 0 exceeds the order of ε 2 in the numerical solution, our MMS solution fails to maintain its accuracy. This problem, however, occurred less often than the aforementioned violation of ordering in theχ(ω) expansion. In order to resolve these issues, more extensive work is required.

Appendix A
In this appendix we introduce a more convenient representation of the linear polarization (2.4). First, let us mention that our convention for Fourier transform is Then using (A.1) we have according to (2.4) whereχ(ω) is the Fourier transform of χ(t).

Appendix B
In this appendix we derive the amplitude equation and the MMS solution to (3.2.4). The multiple scales were introduced in (2.20). Using these on the equation (3.2.4), we get the following perturbation hierarchy ε 0 : a∂ t 0 t 0 t 0 e 0 + (γ + 1)∂ t 0 t 0 e 0 − a∂ z 0 z 0 t 0 e 0 − γ∂ z 0 z 0 e 0 = 0, and where ω = ω(k) is a complex number and a solution to the dispersion relation Note that k is the scaled initial spatial frequency and is used to obtain the complex frequency ω.
We now proceed to the ε order equation. Inserting (B.4) into (B.2) we get In order to remove secular terms we postulate that Observe that by differentiating the dispersion relation (B.6) with respect to k we get Using (B.9), the equation (B.8) can be written in the form The ε order equation (B.7) now simplifies into a∂ t 0 t 0 t 0 e 1 + (γ + 1)∂ t 0 t 0 e 1 − a∂ z 0 z 0 t 0 e 1 − γ∂ z 0 z 0 e 1 = 0. We now compute the right-hand side of the order ε 2 equation. Inserting (B.12) into the right-hand side of the order ε 2 equation (B.3) we get are the non-secular terms and where ω i = Im ω. In order to remove secular terms we postulate that From the equation (B.8) we express the ∂ t 1 z 1 derivative in terms of ∂ t 1 t 1 in the following way Substituting (B.16) back to (B.15) and using the relation (B.9) by then term ∂ t 2 A 0 we get Let us deal with the factor by the term ∂ t 1 t 1 A 0 separately.
Now we turn once again for help to the dispersion relation (B.6) and find (B.20) Using this result in (B.17), our final amplitude equation is By removing the secular terms from the equation (B.13), the order ε 2 equation turns to which we solve for e 2 taking only the particular solution. where , (B.24) where ω i = Im ω.
Using the same argument as in (2.48) we can simplify the amplitude equation (C.33) as ∂ t A + ω (k)∂ z A − i∂ zz A β − α (ω (k)) 2 = 0, (C.38) which can be solved as an initial value problem.