## Abstract

We present the generation of optical pulses with a spectral range of 500–2400 nm and energies up to 10 µJ at 1 kHz repetition rate by cascaded second-order nonlinear interaction of few-cycle pulses in beta-barium borate (BBO). Numerical simulations with a 1D+time split-step model are performed to explain the experimental findings. The large bandwidth and smooth spectral amplitude of the resulting pulses make them an ideal seed for ultra-broadband optical parametric chirped pulse amplification and an attractive source for spectroscopic applications.

© 2016 Optical Society of America

## 1. Introduction

Few-cycle light pulses with mJ-scale pulse energies are of high interest to produce isolated attosecond pulses in gases and on surfaces [1,2]. Classical laser systems have proved incapable of directly generating or amplifying such pulses due to the lack of broadband gain materials that would allow to enter the few-cycle regime. Therefore these systems – e.g. Ti:Sapphire systems – are used in combination with nonlinear broadening techniques like hollow core fiber (HCF) compression to overcome this limitation [3]. Such schemes have been in use for several years and are able to generate pulses with durations of less than 4 fs and energies up to 3.4 mJ [4].

For high-energy attosecond pulse generation, however, input energies of a few mJ are not sufficient. A common method to further amplify few-cycle pulses is optical parametric chirped-pulse amplification (OPCPA), where a strong narrowband pump beam amplifies a weak broadband signal inside a crystal possessing a second-order nonlinearity. [5–8]. It has been demonstrated that with this concept pulse energies of 80 mJ together with pulse durations of less than 5 fs are achievable [9].

On the way to even higher energies with OPCPA systems the damage threshold and the available aperture size of commonly used nonlinear crystals like BBO and Lithium Triborate (LBO) limit the upscaling. The Petawatt-Field-Synthesizer (PFS) project that aims for multi-Joule pulse energies tries to work around this issue by employing DKDP-crystals which can be grown to very large apertures [10, 11]. With these crystals and an Yb:YAG-based pump system [12,13] (frequency-doubled central wavelength of 515 nm) the phase-matchable spectral range for the seed is 700–1400 nm – a rather exotic wavelength band for high-energy OPA systems.

The generation of appropriate seed pulses in this spectral range is a non-trivial task [14, 15] and was the main motivation for the work described in this publication. The scheme we present in the following was adapted from Fattahi et al. [16] who used difference frequency generation (DFG) in BBO to generate pulses in the range of 1000–2500 nm for seeding a short-wavelength infrared (SWIR) OPCPA system. In contrast to their setup where the DFG was optimized for maximal energy around 2000 nm and long wavelength components were separated spectrally from the input pulses we took advantage of cascaded nonlinear interactions to produce pulses with an ultra-broad spectrum of 500–2400 nm and separated the output from the input by polarization selection.

## 2. Experimental setup

The experiment is set up as follows: The compressed output of a commercial multi-pass Ti:Sapphire amplifier (1.5 mJ, *λ*_{0} = 790 nm, 25 fs, 1 kHz) is focused into a gas-filled hollow core fiber and spectrally broadened to 500–950 nm. The pulses are compressed by chirped mirrors and a pair of fused-silica wedges to about 4.7 fs (Fourier limit 4.1 fs) and are linearly polarized by a thin nano-particle polarizer (Thorlabs Inc.). A reflective filter wheel is used to fine control the power. The beam is then focused with a spherical mirror into a BBO crystal (Type I, *ϑ*=20°, *φ*=90°, *d*=500 µm) and afterwards recollimated by another spherical mirror. A second polarizer is used to filter out ordinary polarized components which are then analyzed with a power meter and a set of three calibrated broadband spectrometers (Ocean Optics).

A crucial point for the experiment is the rotation angle *ψ* of the BBO crystal (the roll-angle in the “yaw-pitch-roll” terminology): If the crystal is oriented such that the polarization of the input beam is along the extra-ordinary axis the output beam is also exclusively extra-ordinary polarized (see Fig. 1(a)) and no light passes the ordinary-oriented polarizer. However, if one rotates the BBO by *ψ* around the optical axis the input beam splits when it enters the crystal into ordinary and extra-ordinary components whose electric field amplitudes *A*_{o} and *A*_{e} depend on the total input field amplitude *A*_{in} and on the angle *ψ*:

The energy that is contained in the respective polarization directions is then *E*_{e/o} ~ |*A*_{e/o}|^{2}. This scaling was confirmed experimentally with weak, temporally stretched pulses to suppress nonlinear contributions (see blue curve in Fig. 1(b)). If one increases the intensity to several TW/cm^{2} in the focus (*E*_{in} = 60µJ) and spectrally filters the output beam with a long-pass filter (RG1000) the nonlinear generation of new wavelengths >1000 nm can be detected (red curve in Fig. 1(b)). A strong dependence of the energy in this spectral region on *ψ* is visible as well as the existence of an angle where no new frequencies are generated. Throughout the experiments *ψ* was set to about 5° as this angle proved to give best results in terms of energy and spectral shape.

The measurement of the spectrum of the o-polarized output pulses reveals a continuum that spans more than two octaves from 500–2400 nm (see Fig. 2). This continuum contains spectral bands of different origins:

- In the region 500–950 nm we first of all find the ordinary-polarized fraction of the input beam (“seed” or “signal”). For
*ψ*= 5° its energy equals less than 1% of the total input energy while more than 99% is contained on the e-axis (“pump”). - By optical parametric amplification the long wavelength part (> 700nm) of the signal is amplified by the short wavelength part of the pump, for example:$${\omega}_{\text{pump},\lambda =600\text{nm}}\phantom{\rule{0.2em}{0ex}}\left(\mathrm{e}\right)\phantom{\rule{1em}{0ex}}\to \phantom{\rule{1em}{0ex}}{\omega}_{\text{signal},\lambda =825\text{nm}}\phantom{\rule{0.2em}{0ex}}\left(\mathrm{o}\right)\phantom{\rule{1em}{0ex}}+\phantom{\rule{1em}{0ex}}{\omega}_{\text{idler},\lambda =2200\text{nm}}\phantom{\rule{0.2em}{0ex}}\left(\mathrm{o}\right)$$This amplification of the signal boosts the output energy in the range of 700–950 nm significantly. By the same process the SWIR components above ~ 1400nm are generated as the difference frequency (“idler”) of pump and signal – represented by the 2200 nm product in Eq. (2). This interaction has been described in [16].
It should be pointed out that the broad bandwidth of the interaction would usually (for a narrow-band pump) require a non-collinear geometry for phase-matching. Because of the broadband pump employed here this requirement is relaxed, permitting a collinear geometry.

The most interesting spectral region is that between ~ 950–1400 nm which contains a large fraction of the total energy. The generation of photons in this band is at first glance surprising since these spectral components cannot be produced efficiently by a single

*χ*^{(2)}-process given the range of the input spectra.

A strong influence of *χ*^{(3)}-processes on the other hand can be excluded: if cross-phase-modulation (XPM) by the pump would cause a broadening of the signal into the infrared we should also see self-phase-modulation (SPM) in the pump spectrum (which we did not) as both effects generally appear at a similar intensity level [17]. Additionally the spectral broadening via XPM would strongly depend on intensity, whereas we can see in Fig. 2 that the spectral *shape* between 950–1400 nm is hardly changing with increasing input energy. Cross-polarized-wave (XPW) generation [18] does also not serve as a proper explanation since the XPW would not vanish for *ψ* = 0 which is the case in our measurements as already shown in Fig. 1(b).

Our explanation for the findings is therefore the cascading of two *χ*^{(2)}-processes, an effect that has been reported before in other experiments [19–21]. An exemplary interaction would be:

So the idler photons generated in a first DFG process by one pump wavelength (600 nm) interact subsequently with pump photons of a different wavelength (700 nm) in a second process and generate photons in the sought-after spectral region. Both processes are well phasematched by the employed BBO crystal leading to an efficient cascading.

## 3. Simulation

To verify our explanation that cascading processes contribute significantly to the generation of the measured output pulses, we performed simulations with a 1D+time split-step method that will be described in the following.

Many analytical and numerical simulations of nonlinear interactions of light waves limit themselves to interactions of discrete wavelengths or specific, pre-selected wavelength bands for the reason of simplicity, of faster computation or to keep problems analytically solvable [22, 23]. Usually this limitation is justified due to the small bandwidths involved or because the products of the interaction are known and parasitic effects can be neglected. However, for few-cycle pulses as described in this work very broad bandwidths and possible cascaded interactions have to be taken into account. Likewise, the common distinction of pump, signal and idler beam by specific wavelength bands becomes problematic in degenerated cases when their spectral regions overlap.

Therefore our simulation follows a more general approach similar to [24] and distinguishes pulses (or rather their electric fields) not by spectral regions but only by their polarization direction. Hence the model describes the propagation and nonlinear interaction of extra-ordinary and ordinary polarized pulses that both may span the whole spectral range of interest. Propagation and interaction are calculated in the frequency- and time-domain respectively by subsequently Fourier-transforming between the two domains at each calculation step (split-step). Because of the small thickness of the employed crystal and the long Rayleigh-length of the beam we could neglect effects like walk-off, diffraction and self-focusing.

For the calculations we start with the equation for the propagation of ultra-short pulses in the slowly-varying-envelope approximation from Appendix B of [25] omitting the diffraction term for simplicity:

*n*and the nonlinear polarization

_{i}The indices *i, j, k, l* denote the *X*-, *Y*- and *Z*-coordinates and
$\mathfrak{F}$ is the Fourier transformation. Throughout our calculations we use the Einstein sum convention.

Taking into account only second-order effects this yields a set of nonlinear equations for the pulse development with *Z*-propagation:

This set of equations describes the propagation and interaction of fields in the crystal coordinate system *X,Y,Z*. However, to accomplish phase-matching the beam propagation vector (defined by the *z*-axis of a beam coordinate system *x,y,z*) has to be rotated with respect to the crystal system. For easier comparison we follow the convention [26] to rotate the beam system first by *φ* around the *z*-axis, then by *ϑ* around the *y*-axis and finally by *ψ* around (again) the *z*-axis (see Fig. 3). The more practical view of this transformation is from a fixed beam coordinate system (i.e. the lab) where consequently the crystal is rotated. The respective transformation of crystal coordinates into lab coordinates is given by the matrix

With this matrix the effective nonlinear tensor ${\chi}_{lmn}^{(2)}{}^{\prime}$ in lab coordinates (indicated by the dash) can be calculated from the tensor ${\chi}_{ijk}^{(2)}$ in crystal coordinates by performing three matrix multiplications:

where*T*

^{−1}is the inverse matrix of

*T*.

To illustrate the idea we choose *ψ* = 0 and take a look at the component
${\chi}_{xyy}^{(2)}{}^{\prime}$. For a crystal with a 3m point group like BBO the
${\chi}_{ijk}^{(2)}\text{-tensor}$ has eleven non-zero components with three independent values, namely *d*_{22}, *d*_{31} and *d*_{33} [26]. Performing the matrix multiplications of Eq. (8) the value of
${\chi}_{xyy}^{(2)}{}^{\prime}$ is 2(*d*_{31} sin*ϑ − d*_{22} sin3*φ* cos*ϑ*) reproducing the well-known *d*_{eff} value from literature. In the typical case of *φ* = 90° the *x*-axis is the extra-ordinary axis and the *y*-axis the ordinary, hence
${\chi}_{xyy}^{(2)}{}^{\prime}$ (like *d*_{eff}) defines the coupling constant for the interaction *o* + *o → e*.

The advantage of the matrix formalism over a single *d*_{eff}-value is that it contains the effective coefficients for *all* possible interactions and that the tensor-dimension can be easily extended to include *χ*^{(3)}-, *χ*^{(4)}-,… interactions as well.

In Eq. (6) one can now replace
${\chi}_{ijk}^{(2)}$ by the obtained effective tensor
${\chi}_{lmn}^{(2)}{}^{\prime}$ of the rotated crystal. To solve the resulting differential equation numerically we divide the *z*-propagation distance (i.e. the crystal length) into *S* discrete slices of thickness ∆*z* and apply the Exponential Euler Method [27] to calculate the electric fields at a slice *s* + 1 from the fields at the previous slice *s*:

Note that also the electric fields are now in lab coordinates. For a better numerical stability reduced wavenumbers
${k}_{l}^{\ast}={k}_{l}-{k}_{0}$ are used, i.e. the wavenumber *k*_{0} of the central frequency is subtracted from the wavenumbers *k _{l}* (“moving frame”). Furthermore absorption is taken into account. Implementing Eq. (9) into MATLAB

^{®}with

*S*= 1000 slices and 2

^{14}equidistant

*ω*-points the computation takes less than two seconds on one core of an Intel

^{®}Xeon E3-1240.

For simplified calculation we choose *ψ* = 0 to have *x* as the extra-ordinary and *y* as the ordinary axis. The de facto polarization rotation of *ψ* = 5° is accounted for by applying Eq. (1) and distribute the input energy onto *x*- and *y*-axis accordingly. For the respective input fields
${{\tilde{E}}^{\prime}}_{x,0}\left(\omega \right)$ and
${{\tilde{E}}^{\prime}}_{y,0}\left(\omega \right)$ the experimentally measured spectra as well as the retrieved spectral phase obtained with the transient-grating frequency-resolved optical gating (TG-FROG) technique were used. The beam profile at the focus was taken to be Gaussian with 380 µm full width at half maximum (FWHM) and circular symmetry in good agreement with a CCD-measurement of the real beam. In order to take the spatial dependence of the intensity into account the beam profile was segmented into 200×200 pixels with corresponding intensity values. Simulations were performed for a range of 60 different input intensities to create a lookup-table of generated spectra. Subsequently each pixel can be linked to a spectrum in the lookup-table according to its intensity value. The summed up spectra of all pixels yield then the global spectrum (quasi-2D+time).

The results of these calculations – i.e. the generated global spectra – are shown in Fig. 4. To be able to compare them with the experimental data the same total energies were used as input parameter. As one can see, the forming of the multi-octave-spanning spectrum is qualitatively reproduced. Especially the efficient generation of new frequency components at 950–1400 nm is well confirmed. Quantitatively the order of magnitude of the output energy is correctly calculated but a higher conversion efficiency compared to the experimental results is predicted. This might originate from overestimating the peak intensity of the input pulse due to temporal and spatial artifacts from the hollow-core fiber broadening and the chirped mirror compression or from imperfections in the BBO crystal.

To visualize the cascading nature of the nonlinear interactions we also performed a simulation with narrow spectral bands but otherwise unchanged parameters (see Fig. 5). The small bandwidths allow us to identify the underlying frequency-mixing processes and the result confirms the contribution of cascaded interactions to the formation of the multi-octave spectra of Fig. 4.

It should be noted at this point that if carrier-envelope-phase (CEP) stability is of interest for an application, only (cascaded) DFG-products generated by an even number of interacting wavelengths are passively phase-stable, while the other products inherit the phase of the input pulses. This is also important if the entire spectral bandwidth is planned to be used since in the case of a non-CEP-stable input beam a variable phase jump at the border (~ 1400nm) between these two spectral regions would yield different temporal compression for each shot. However, simulations show that for actively stabilized systems CEP fluctuations of ±300mrad around CEP=0 and intensity fluctuations of up to ±20% are tolerable to maintain a stable waveform of the output pulses. Under these practically feasible conditions the achievable pulse duration of ~ 2.6fs (sub-cycle) is lengthened by not more than 10%. The spectral phase of these pulses is – dominated by dispersion – smooth and can be well approximated with a polynomial function (dispersion coefficients: GDD=28 fs^{2}, TOD=37 fs^{3}, FOD=−71 fs^{4}, *λ*_{central}=1000 nm). Hence good compression with standard tools is expected. A test with chirped mirrors designed for 800–1300 nm experimentally confirmed the compressibility of the pulses in at least this spectral range.

## 4. Conclusions

We have demonstrated that by making use of (cascaded) second-order nonlinear effects in BBO one can generate pulses with a spectral range of 500–2400 nm from 4.7 fs (500–950 nm) input pulses. The experimentally achieved overall conversion efficiency was up to 15%. With numerical simulations we could reconstruct the experimental findings and identify the processes contributing to the output spectrum.

The scheme described in this publication provides a simple, efficient, scalable and therefore attractive way of generating broadband seed pulses for OPA systems in the near- to shortwave-IR region. If a high contrast of the seed pulses is required the scheme even allows for the integration of an XPW stage: by placing an appropriate barium fluoride crystal directly in front of the BBO and setting the roll-angle *ψ* of the BBO to zero the ordinary-polarized components (that take part in the subsequent nonlinear interactions) originate solely from the XPW generation process. The temporal cleaning by this process ensures that no artifacts of the few-cycle input pulses are transmitted. A conducted proof-of-principle experiment with this setup demonstrated already its feasibility.

Simulations suggest that provided a CEP-stable input and an ultra-broadband compression scheme pulse durations of about 2.6 fs could be achieved rendering the generated pulses potentially useful for ultrafast spectroscopy.

## Acknowledgments

The authors would like to thank Ferenc Krausz for his support. This work was funded through the PFS grant of the Max-Planck Society, by EURATOM-IPP and by DFG through the Cluster of Excellence Munich-Centre for Advanced Photonics (MAP EXC 158).

## References and links

**1. **W. Schweinberger, A. Sommer, E. Bothschafter, J. Li, F. Krausz, R. Kienberger, and M. Schultze, “Waveform-controlled near-single-cycle milli-joule laser pulses generate sub-10 nm extreme ultraviolet continua,” Opt. Lett. **37**, 3573–3575 (2012). [CrossRef] [PubMed]

**2. **P. Heissler, R. Hörlein, J. M. Mikhailova, L. Waldecker, P. Tzallas, A. Buck, K. Schmid, C. M. S. Sears, F. Krausz, L. Veisz, M. Zepf, and G. D. Tsakiris, “Few-cycle driven relativistically oscillating plasma mirrors: A source of intense isolated attosecond pulses,” Phys. Rev. Lett. **108**, 235003 (2012). [CrossRef] [PubMed]

**3. **M. Nisoli, S. D. Silvestri, O. Svelto, R. Szipöcs, K. Ferencz, C. Spielmann, S. Sartania, and F. Krausz, “Compression of high-energy laser pulses below 5 fs,” Opt. Lett. **22**, 522–524 (1997). [CrossRef] [PubMed]

**4. **F. Böhle, M. Kretschmar, A. Jullien, M. Kovacs, M. Miranda, R. Romero, H. Crespo, P. Simon, R. Lopez-Martens, and T. Nagy, “Generation of 3-mJ, 4-fs CEP-stable pulses by long stretched flexible hollow fibers,” in “Research in Optical Sciences: Postdeadline Papers,” (Opt. Soc. Am., 2014), p. HW5C.2.

**5. **I. Ross, P. Matousek, M. Towrie, A. Langley, and J. Collier, “The prospects for ultrashort pulse duration and ultrahigh intensity using optical parametric chirped pulse amplifiers,” Opt. Commun. **144**, 125–133 (1997). [CrossRef]

**6. **A. Baltuška, T. Fuji, and T. Kobayashi, “Visible pulse compression to 4 fs by optical parametric amplification and programmable dispersion control,” Opt. Lett. **27**, 306–308 (2002). [CrossRef]

**7. **D. Brida, C. Manzoni, G. Cirmi, M. Marangoni, S. Bonora, P. Villoresi, S. D. Silvestri, and G. Cerullo, “Few-optical-cycle pulses tunable from the visible to the mid-infrared by optical parametric amplifiers,” J. Opt. **12**, 013001 (2010). [CrossRef]

**8. **O. Mücke, S. Fang, G. Cirmi, G. Rossi, S.-H. Chia, H. Ye, Y. Yang, R. Mainz, C. Manzoni, P. Farinello, G. Cerullo, and F. Kärtner, “Toward waveform nonlinear optics using multimillijoule sub-cycle waveform synthesizers,” IEEE J. Sel. Top. Quantum Electron. **21**, 1–12 (2015). [CrossRef]

**9. **L. Veisz, D. Rivas, G. Marcus, X. Gu, D. Cardenas, J. Mikhailova, A. Buck, T. Wittmann, C. M. S. Sears, S.-W. Chou, J. Xu, G. Ma, D. Herrmann, O. Razskazovskaya, V. Pervak, and F. Krausz, “Generation and applications of sub-5-fs multi- 10-TW light pulses,” in “2013 Conference on Lasers and Electro-Optics Pacific Rim,” (Opt. Soc. Am., 2013), p. TuD2_3.

**10. **Z. Major, S. A. Trushin, I. Ahmad, M. Siebold, C. Wandt, S. Kliengebiel, T.-J. Wang, J. A. Fülop, A. Henig, S. Kruber, R. Weingartner, A. Popp, J. Osterhoff, R. Hörlein, J. Hein, V. Pervak, A. Apolonski, F. Krausz, and S. Karsch, “Basic concepts and current status of the petawatt field synthesizer - a new approach to ultrahigh field generation,” The Review of Laser Engineering **37**, 431–436 (2009). [CrossRef]

**11. **C. Skrobol, I. Ahmad, S. Klingebiel, C. Wandt, S. A. Trushin, Z. Major, F. Krausz, and S. Karsch, “Broadband amplification by picosecond OPCPA in DKDP pumped at 515 nm,” Opt. Express **20**, 4619–4629 (2012). [CrossRef] [PubMed]

**12. **S. Klingebiel, C. Wandt, C. Skrobol, I. Ahmad, S. A. Trushin, Z. Major, F. Krausz, and S. Karsch, “High energy picosecond Yb:YAG CPA system at 10 Hz repetition rate forpumping optical parametric amplifiers,” Opt. Express **19**, 5357–5363 (2011). [CrossRef] [PubMed]

**13. **C. Wandt, S. Klingebiel, S. Keppler, M. Hornung, M. Loeser, M. Siebold, C. Skrobol, A. Kessel, S. A. Trushin, Z. Major, J. Hein, M. C. Kaluza, F. Krausz, and S. Karsch, “Development of a joule-class Yb:YAG amplifier and its implementation in a CPA system generating 1 TW pulses,” Laser Photon. Rev. **8**, 875–881 (2014). [CrossRef]

**14. **I. Ahmad, S. Trushin, Z. Major, C. Wandt, S. Klingebiel, T.-J. Wang, V. Pervak, A. Popp, M. Siebold, F. Krausz, and S. Karsch, “Frontend light source for short-pulse pumped OPCPA system,” Appl. Phys. B **97**, 529–536 (2009). [CrossRef]

**15. **A. Kessel, C. Skrobol, S. Klingebiel, C. Wandt, I. Ahmad, S. A. Trushin, Z. Major, F. Krausz, and S. Karsch, “Generation and optical parametric amplification of near-IR, few-cycle light pulses,” in “CLEO: 2014,” (Opt. Soc. Am., 2014), p. SM3I.3.

**16. **H. Fattahi, A. Schwarz, S. Keiber, and N. Karpowicz, “Efficient, octave-spanning difference-frequency generation using few-cycle pulses in simple collinear geometry,” Opt. Lett. **38**, 4216–4219 (2013). [CrossRef] [PubMed]

**17. **G. P. Agrawal, *Nonlinear Fiber Optics*, 3rd ed. (Academic, 2001).

**18. **N. Minkovski, G. I. Petrov, S. M. Saltiel, O. Albert, and J. Etchepare, “Nonlinear polarization rotation and orthogonal polarization generation experienced in a single-beam configuration,” J. Opt. Soc. Am. B **21**, 1659–1664 (2004). [CrossRef]

**19. **H. Tan, G. P. Banfi, and A. Tomaselli, “Optical frequency mixing through cascaded second-order processes in β-barium borate,” Appl. Phys. Lett. **63**, 2472–2474 (1993). [CrossRef]

**20. **G. I. Petrov, O. Albert, N. Minkovski, J. Etchepare, and S. M. Saltiel, “Experimental and theoretical investigation of generation of a cross-polarized wave by cascading of two different second-order processes,” J. Opt. Soc. Am. B **19**, 268 (2002). [CrossRef]

**21. **J. Matyschok, T. Lang, T. Binhammer, O. Prochnow, S. Rausch, M. Schultze, A. Harth, P. Rudawski, C. L. Arnold, A. L’Huillier, and U. Morgner, “Temporal and spatial effects inside a compact and CEP stabilized, few-cycle OPCPA system at high repetition rates,” Opt. Express **21**, 29656–29665 (2013). [CrossRef]

**22. **R. Boyd, *Nonlinear Optics*, 3rd ed. (Academic, 2008).

**23. **G. Arisholm, “General numerical methods for simulating second-order nonlinear interactions in birefringent media,” J. Opt. Soc. Am. B **14**, 2543–2549 (1997). [CrossRef]

**24. **T. Lang, A. Harth, J. Matyschok, T. Binhammer, M. Schultze, and U. Morgner, “Impact of temporal, spatial and cascaded effects on the pulse formation in ultra-broadband parametric amplifiers,” Opt. Express **21**, 949–959 (2013). [CrossRef] [PubMed]

**25. **T. Brabec and F. Krausz, “Intense few-cycle laser fields: Frontiers of nonlinear optics,” Rev. Mod. Phys. **72**, 545–591 (2000). [CrossRef]

**26. **R. L. Sutherland, *Handbook of Nonlinear Optics*, 2nd ed. (Marcel Decker, Inc., 2003).

**27. **M. Hochbruck and A. Ostermann, “Exponential integrators,” Acta Numerica **19**, 209–286 (2010). [CrossRef]