3 i.e. d dz ′ ∫ + ∞ −∞ | A ( z $^{′}$, τ ) | $^{2}$dτ = 0. It can also be shown that the total energy is conserved even if the non ap- proximated nonlinear polarization envelope A$_{p}$ [Eq.(10)] is used. we set the reference frequency ω$_{0}$ to be equal to the sec- ond harmonic of the input pulse: in this way the second harmonic is stationary in the reference frame ( z $^{′}$, τ ). We solved Eq. (11) by split-step Fourier method ex- ploiting fourth order Runge-Kutta scheme for the non- linear step. Figure 1a) shows the evolution of the electric field enve- lope amplitude | A | from numerical solution of Eq. (11). We can see the typical scenario of the propagation of femtosecond pulses in highly group velocity mismatched (GVM) process: the fundamental frequency (FF) pulse generates its second harmonic (SH) during propagation, and the generated SH pulse has the typical shape of a ini- tial peak followed by a long tail, whose duration is fixed by the product between GVM and crystal length. Figure 1b) shows the electric field envelope amplitude at the end of the crystal. It can be seen a peak that corresponds to the faster frequency components located around FF, fol- lowed by a long tail that ends with a second lower peak. This long pulse corresponds to the generated SH compo- nents. The SH pulse is smooth, indicating that no beat- ing with eventual FF components is present. Whereas in the residual FF pulses centered around τ ≈ − 1 . 6 ps , there is a clear fast oscillation, indicating that FF and SH components are superimposed. To test the results, we simulated the same set-up with a standard coupled wave model [13], by inserting the val- ues of first and second order dispersion evaluated at FF and SH. To compare the results we filtered A around FF and SH. Figure 1c) shows the electric field amplitudes at the end of the crystal. The results of the two mod- els are practically indistinguishable. It is worth noting that this simulation shows the validity of Eq. (11) over a bandwidth of ω$_{0}$ . As a second example we consider the propagation of a femtosecond pulse into a highly mismatched period- ically poled lithium niobate (PPLN) sample, that was demonstrated experimentally to generate an octave span- ning supercontinuum spectral broadening [10]. To model the refractive index dispersion we employed a Sellmeier model fitted from experimental data [14] and nonlinear coefficient is d$_{33}$ = χ (2) $_{LN}$/ 2 = 27 pm/V . In the numerical code we inserted the exact dispersion relation k ( ω ). We assumed a QPM grating with a period Λ = 30 µm (phase matched for second harmonic generation at around 2 µm fundamental wavelength). We included higher order QPM terms, since the huge bandwidth can phase match different spatial harmonics. We injected a T = 50 fs FWHM long gaussian pulse, centered around 1580 nm , with I = 15 GW/cm 2 peak intensity. In the simulation we set the reference wavelength λ$_{0}$ = 700 nm . Figure 2a) shows the evolution of the spectrum during the propaga- tion into a L = 7 mm crystal. We can see a consistent broadening and redshift of the FF part of the spectrum that, at the end of the crystal, reaches an octave-spanning bandwidth from 1200 nm to 3000 nm . We can also see the generation of spectral components at the second and third harmonics. At the second harmonic the spectrum initially broadens and has an evolution ruled by highly In order to show the validity of our equation, we simulated the propagation of a femtosecond pulse in a L = 5 mm long periodically poled lithium tantalate sam- ple (PPLT). To model the refractive index dispersion we employed a Sellmeier model fitted from experimen- tal data [12] and nonlinear coefficient is d$_{33}$ = χ (2) $_{LT}$/ 2 = 10 . 6 pm/V . In the numerical code we inserted the ex- act dispersion relation k ( ω ). We assumed a first or- der quasi phase matching (QPM) grating, with a period Λ = 17 . 4 µm . We thus allowed a periodic variation of the nonlinear coefficent χ $^{(2)}$= χ $^{(2)}$( z ) = 2 /πχ (2) $_{LT}$e i 2 π/ Λ $^{z}$+c.c.. We injected a T = 60 fs FWHM long gaussian pulse, centered around 1400 nm , with I = 10 GW/cm 2 peak in- tensity. The corresponding residual phase mismatch is Δ k = 2 k ( ω$_{in}$ ) − k (2 ω$_{in}$ ) = 10000 m − $^{1}$, where ω$_{in}$ is the carrier frequency of the input pulse. In the simulation FIG. 1: Propagation of a femtosecond pulse in a PPLT crys- tal. a) Evolution of the field amplitude | A | from numerical solution of Eq. (11). b) Electric field amplitude at the crystal output. c) Comparison between coupled wave solution (dot- ted curves) and | A | filtered around fundamental and second harmonic (solid curves). The initial pulse has gaussian shape and the parameters are T = 60 fs , I = 10 GW/cm $^{2}$, λ in = 1400 nm , λ$_{0}$ = 2 πc/ω$_{0}$ = 700 nm , d$_{33}$ = χ (2) $_{LT}$/ 2 = 10 . 6 pm/V −2 −1.5 −1 −0.5 0 0.5 0 0.5 1 1.5 2 x 10 8  [ps] | A |  [V/m] b) −2 −1.5 −1 −0.5 0 0.5 0 0.5 1 1.5 2 x 10 8  [ps] | A |  [V/m] Input FF SH c)