Open Access
Add to BibSonomyAbstract
An economic, space- and time-resolved method to model ultra-short, intense pulse propagation in waveguides is described. Simulations of supercontinuum generation on a chip demonstrate the utility of the approach. Comparisons with the generalized nonlinear Schrödinger equation elucidate spatial effects, which influence pulse dynamics and the generation of new spectral components.
© 2014 Optical Society of America
1. Introduction
In recent years, there has been much research in nonlinear optics pertaining to ultrashort pulses propagating in what can be characterized as “confined geometries.” Beyond the traditional nonlinear fiber optics, this includes extreme nonlinear processes in hollow-core fibers [1–5], supercontinuum generation in micro-structure fibers [6–10], high-harmonic generation in capillaries [11], femtosecond filamentation in slab waveguides [12–15], and on-chip supercontinua and frequency-comb generation in integrated, highly nonlinear waveguides [16–21]. Naturally, the development of methods for computer-aided modeling has been crucially important for these areas as interpretation of experiments must often be supported by numerical simulations.
The methods most often used to simulate nonlinear pulse propagation in waveguides are based on the projection of the optical field onto the fundamental mode of the given waveguide. The result is the well-known Generalized Nonlinear Schrödinger Equation, which can be applied in many different situations. It is able to and produce simulated spectra closely comparable to the experimental ones [22, 23]. Another well established model option is represented by the Forward Maxwell Equation (FME), which made it possible to predict the soliton fission mechanism of supercontinuum generation in fibers [24–26].
However, the propagation models that rely on the single transverse mode approximations have drawbacks, which have been addressed in various ways (e.g., [27, 28]). One of the issues is that in some cases the shape of the mode changes strongly with the wavelength and that the effective nonlinearity varies significantly over the bandwidth of the pulse. Another issue arises in hollow waveguides where it may be necessary to account for higher-order modes. Multimode approaches can be used in capillaries [29], and in waveguides with complex geometries [30–32]. These methods straddle the space between single-mode projections and approaches with full spatial resolution. While going beyond single-mode methods in that they can capture the multimode dynamics as well as to some degree the dependence of the modal properties on on the wavelength, they are also computationally significantly faster than spatially resolved methods, including the one discussed here. On the other hand, coupled-mode methods require modal coupling coefficients, and these must be evaluated for the specific nonlinear interactions. In effect, the propagation equations depend on which type of light-matter interactions is taken into account.
A seemingly obvious answer to these challenges would be a fully space- and time-resolving simulation of optical pulses in waveguide structures. Indeed, the research community interested in beam-propagation methods have addressed this problem via both the classical beam-propagation techniques as well as time-domain pulsed-beam modeling [33–38]. Nevertheless, problems remain on a practical level because the fully resolved simulations inherit difficulties from two sides. On the one hand, they require short integration steps just as the Beam-Propagation Method (BPM) approaches suitable for waveguides with strong material contrast. On the other hand, they must propagate wide-bandwidth optical waveforms over significant distances. As a consequence, such an approach necessitates a truly large-scale numerical simulation, and this is just one of the reasons that, despite the increasing computing power, such methods are far from commonplace.
It is the aim of this paper to demonstrate the next step toward practical applications of fully resolved modeling of extremely nonlinear processes in waveguides. In particular:
- We describe application of the generalized unidirectional pulse propagation equations to a waveguide structure with sharp material interfaces and strong refractive index contrast. Our implementation is capable of modeling femtosecond pulse propagation over experimentally relevant propagation distances.
- We demonstrate that simulations of highly nonlinear pulse propagation in a waveguide with full space- and time-resolution, while being orders of magnitude more computationally expensive than modeling based on the generalized nonlinear Schrödinger equation (GNLS), become feasible with commonly available computing resources.
- We use comparative simulations, with the traditional GNLS method on one side, and our fully resolved simulation approach on the other, to elucidate higher-order-mode related physical effects that contribute to shaping of the generated supercontinuum and pulse dynamics in general.
- We find a natural description beyond the commonly used fundamental-mode GNLS in terms of what we call the “spatial envelope” that describes the nonlinear deformation of the propagating waveform as compared to the field shape of the fundamental mode.
What sets our approach apart from other methods designed specifically for optical pulse propagation in waveguides is that it captures arbitrary material chromatic dispersion properties in the spectral domain, and that the treatment of medium nonlinearities is completely separated from the beam propagation. These properties allow us to apply the method to situations with ultra-fast dynamics and extreme spectral broadening, without any restrictions as to what type of nonlinear interactions one can deal with. From the implementation point of view, it does not matter, for example, if the core is filled with water with weak ionization, and.or higher-order nonlinearity, or if the material exhibits two-photon absorption.
In what follows, we briefly recapitulate the theoretical background of the generalized Unidirectional Pulse Propagation Equations (gUPPE) [39] which serve as the implementation framework for the simulator used in this work. Then we describe the approximations, simulation parameters and our model waveguide structure, and continue with discussion of the results from our comparative simulations. Finally, in our discussion we point out how the insight from this work could help to build even faster space- and time-resolved simulators of pulse propagation in confined geometries.
2. Generalized unidirectional pulse propagation equations
Generalized Unidirectional Pulse Propagation Equation (gUPPE) is a mixed-representation pulse propagator. The optical field is described in terms of its spectral amplitudes, A⃗(ω, x, y, z) which are related to the temporal spectrum of the pulse, E⃗(ω, x, y, z), through
where the linear operator L accounts for the chromatic dispersive material properties, described by z-invariant permittivity ε(x, y, ω), and for the transverse structure of the waveguide. It operates on the transverse components of a given vector field Ψ⃗, In this expression, ∇⊥ and Δ⊥ operate in the transverse space. Of course, L is nothing but the operator that appears in the Helmholtz wave equation, and which is frequently used in beam-propagation problems. We will only consider piece-wise constant ε, so that the last term vanishes, but account must be taken of proper conditions for field values and its derivatives at material interfaces.The rationale for working with A instead of E is that spectral amplitudes are the natural “slow variables” which only change with z if nonlinear light-medium interactions are present. Then A obeys this generalized unidirectional pulse propagation equation
where ℱ stands for Fourier transform along temporal and/or frequency dimension, and the nonlinear medium response is represented by the operator N, acting on the real-space representation of the pulse optical field, E(t, x, y, z): Here P⃗(E⃗) encapsulates the model of the nonlinear medium response. It is usually implemented in terms of a function that accepts the history of the optical field, E(t) at a given point (x, y), and returns the corresponding nonlinear polarization (note that the medium response in terms of the induced current density can be treated in an analogous way). It is often convenient to use the analytic signal of the electric field, which is obtained from the spectral amplitudes through a Fourier transform in which the integration is restricted to positive frequencies In practice, the frequency integration interval is further reduced to the so-called active frequencies between ωmin and ωmax, as required by the maximal spectral extent attained.The above equations constitute the pulse propagation framework we use in this work. Note that the light-medium interactions are formally separated from the linear propagation properties of the given structure. This property makes it possible to marry ultra-fast pulse simulation techniques [41] with any of the many available beam-propagation methods.
2.1. Approximations
Derivations of pulse propagation equations naturally lead to pairs of coupled equations, one for the forward and one for the backward propagation direction. If and only if nonlinear interactions are weak enough so that they fail to generate strong backward propagating waves, the system can be restricted to the forward-propagating field alone — this is what we term the unidirectional approximation. In principle, it is the only approximation required. In practice, the implementation of the linear propagator for a structured medium represents a set of uncoupled beam-propagation problems, and as such it requires further approximations. Fortunately, the gUPPE formulation allows one to take advantage of the many methods developed for beam propagation modeling.
In this study, we compare a single-mode propagation method versus a spatially resolved treatment. Since the latter requires significantly more computational effort, we will consider a simple waveguide geometry, and we furthermore adopt the semi-vectorial approximation, with the vector field approximated by its dominant component aligned along and/or perpendicular to the material interfaces of the wave-guide core.
In many situations, including supercontinuum generation on chips [21], which we choose for our illustrations, the pulse propagation can be described quite accurately using solely the fundamental mode of the waveguide. With that in mind, we rewrite the linear propagator as follows:
and choose β(ω) equal to the propagation constant of the fundamental mode. Assumption of a paraxial regime is not necessary for approximating the square root. Rather, we require that the transverse shape of the propagating pulse will be close to that of the fundamental mode. Importantly for our comparison with the single-mode approximation, propagation properties of the fundamental mode are exactly preserved, because it is an eigenvector of L, for which the argument of the second exponential vanishes.2.1.1. GNLS: Single-mode approximation
The standard method for simulation of nonlinear pulse propagation in waveguides is the generalized nonlinear Schrödinger equation (GNLS). To mimic GNLS with gUPPE, one merely needs to choose the appropriate approximation of the linear propagator, namely
and to reduce the transverse computational domain to a single point. It is customary to normalize the optical field in GNLS simulations to the power, and express the nonlinearity in terms of the effective mode area and the waveguide nonlinear constant. In gUPPE, we work with the real electric field on the waveguide axis, which in turn requires an effective value of the nonlinear index. Notwithstanding different normalizations of the computational variables, the gUPPE and GNLS simulation approaches are equivalent provided the latter is properly implemented in the spectral domain. This will be our reference method against which we validate and compare the spatially resolved simulations.2.1.2. gUPPE: Space-resolved simulation
Computationally, gUPPE constitutes a large system of ordinary differential equations (ODEs) that control the evolution of the spectral amplitudes. An ODE solver requires implementation of the right-hand-side in Eq. (3), which in turn needs evaluation of the action of the linear propagator for an arbitrary step:
where A is a vector of amplitudes supplied by the ODE solver [41]. Because L is diagonal in the frequency domain, this operation can be executed independently for each active frequency. This is nothing but an application of a beam-propagation method, which can be realized in a number of ways. We have implemented the Crank-Nicolson (CN) scheme in two transverse dimensions, and the Operator Splitting (OS) scheme, both second-order accurate in Δz. In the 2D CN approach, a linear system of equations is set up for all variables comprised in the transverse cross-section of the optical field. The system connects “old” and “new” beam discretizations for each frequency separately. Each of these frequency slices is therefore updated in one beam-propagation step. In principle, various discretization schemes, and different solvers can be utilized at this stage. We have utilized an iterative sparse-matrix solver with the BICGSTAB algorithm.While we obtained accurate agreement between the two methods, the OS approach turned out to be significantly faster. More importantly, in some (rare) cases the iterative solver in the CN method failed to converge quickly. The reason that the iterative solver is less than optimal for the present application is that a lot of computing time is wasted on frequency slices that carry negligible power and, in the end, do not affect the solution. For these reasons, in the linear propagator we default to using the operator splitting method.
For each alternating direction, x or y, a Crank-Nicolson step is used to advance the solution, e.g.
where discretization of Δxx is done with a five-point stencil, giving fourth-order accuracy. This requires a penta-diagonal linear solver, but allows a coarser grid, which is crucially important for our application. In the vicinity of the core-cladding interface, the discrete Laplacian operator is modified such that the electric field satisfies the appropriate (dis)continuity condition depending on the TE/TM orientation of the interface. Moreover, the five-point, fourth-order scheme allows one to ensure that the normal derivative of the field is continuous (in our semi-vectorial approximation). The discretization naturally depends on where one assumes the interface to be positioned with respect to the closest grid points. In our implementation we assumed that it was located “just outside” of the last grid point within the core.3. Comparative simulation of supercontinuum generation on a chip
For our comparative simulations, we choose a rectangular waveguide made of silicon nitride, embedded in cladding of SiO2. The core dimensions are 1008 nm × 720 nm, and we assume quasi-TE polarization for our semi-vectorial simulations.
We utilize the gUPPE propagation in imaginary distance to obtain the fundamental mode propagation constant β(ω) of this waveguide over a range of frequencies. In order to exclude the influence of the numerical dispersion from our comparison, the same propagation constant is used both in the GNLS, and in the space-resolving approach. This ensures that any deviations observed between the two methods are indeed due to spatial, or higher-order mode effects present in the fully resolved approach.
The chromatic properties of the simulated waveguide are illustrated in Fig. 1. The most important features that control the properties of the generated supercontinua are the two zeros of the group velocity dispersion, located at ∼ 930 nm and ∼ 1420nm. In order to check that our grid spacings are sufficiently fine, calculations were repeated with doubled resolution. The difference in the effective waveguide susceptibility (top panel), which serves as the input for the gUPPE simulator, is barely discernible, while there is only a ∼ 20 nm shift of the zero GVD position (bottom panel). This comparison indicates that with the fourth-order accurate discretization, one can use a relatively coarse grid with Δx = Δy = 72 nm. Taking advantage of the symmetry of the solution, we restrict the computational domain to one quadrant of the waveguide structure. It is important to use transparent boundary conditions, because part of the pulse energy is radiated out, especially during the initial stage of propagation when the input beam adjusts to the waveguide. To keep our grid as small as possible, we have used a novel transparent boundary method that eliminates the smooth transition between the inner domain and the outer perfectly matched layer [40]. This allows us to reduce the number of grid points to 24×24 in the transverse domain. This grid is coupled with the spectral lattice, carrying 816 active frequencies spanning the interval between 8 × 1014 and 4 × 1015 s−1. As a result, our simulation has 959616 degrees of freedom to represent the optical pulse in a waveguide. Such a large number of variables, together with the short integration step and the necessity to propagate the pulsed waveform over a long distance results in a large-scale computation. The gUPPEcore framework has a built-in domain-decomposition parallelization, which we utilize to speed up our simulation runs. For the data sets shown in this paper, we used twelve parallel processes running on a desktop computer. A typical run required about six hours to complete.
Fig. 1 Simulated waveguide dispersion properties. Top-left panel shows the input for the gUPPE simulator, the effective fundamental mode susceptibility, calculated for two different grid resolutions. Bottom-left panel compares the respective group velocity dispersion results. The computational domain, indicated in the scheme, assumes a symmetric solution.
The simulations were initiated with a 200 fs duration pulse with a Gaussian spatial profile adjusted to maximize the coupling into the fundamental mode. Over the first few hundred microns of propagation, a portion of the initial waveform is radiated out and is absorbed by the transparent boundary conditions at the edges of the computational grid. During this initial period, the pulse acquires a slight chirp. For the sake of comparison with the GNLS method, we have adjusted the initial chirp in the single-mode simulations accordingly. Also, because real waveguides exhibit appreciable linear losses, we have added 0.8 dB/cm loss [21] to our waveguide model. While this loss effectively increases the pulse energy required for efficient supercontinuum generation, we did not find qualitative difference from the simulations without the account for losses.
As for the light-matter interactions included in our model, we parametrize the instantaneous Kerr effect with the nonlinear index of n2 = 2.4 × 10−19 m2/W for silicon nitride [42]. The Raman effect, together with two-photon (and higher) absorption has been neglected [21], as was the third-harmonic generation because it falls outside the region of the generated supercontinuum.
3.1. GNLS versus gUPPE
Figure 2 shows supercontinuum spectra simulated in a 40 mm long waveguide for a range of different (output) pulse energies. In both cases, major features correspond to those observed experimentally in similar waveguides [21]. In particular, two spectral peaks due to dispersive waves appear, one on the short wavelength side, and one on the long-wavelength side of the spectrum. With increasing pulse energy, the long-wavelength dispersive peak moves further away from the center of the spectrum, while the short-wavelength peak moves closer, but at a much lower rate. Also, at higher energies, the central portion of the spectrum significantly broadens and shifts toward higher frequencies. This blue shift is responsible for the movement of the dispersive peaks [43].
Fig. 2 Supercontinuum generation in a waveguide after 40 mm propagation. Left: GNLS results. Right: gUPPE simulation with spatial resolution. Labels indicate output pulse energies.
Comparison of the results from the two methods reveals that the space-resolved simulation produces these peaks at slightly different wavelengths, namely closer to the center. Figure 3 shows that the modification is especially evident on the long wavelength side of the spectrum where the difference between the methods is about 80 nm in the peak wavelength. Because both methods exhibit identical dispersion properties when restricted to the fundamental mode, the different locations of the dispersive-wave features in the supercontinuum spectra indicate that higher-order modes play a role in the pulse dynamics.
Fig. 3 Dispersive-wave spectral peak locations as functions of pulse energy for the short-wavelength peak (left) and the long-wavelength side of the spectrum (right).
Let us look closer at the quantitative differences between the compared methods. Figure 4 compares the evolution of the maximal pulse intensity with the propagation distance. One can see the same behavior in both methods, with the GNLS result achieving somewhat higher intensity during the peak events. On the right in Fig. 4, we compare the temporal shape of the pulses at the points along the waveguide where they achieve their highest intensities. In general, the waveform shape is dominated by the central maximum that is surrounded by a lower-intensity pedestal exhibiting a structure due to interference between different spectral components. Once again, the qualitative behavior of the pulses is the same, but with appreciable quantitative differences.
Fig. 4 Single-mode versus space-resolved simulation. Left: the on-axis maximal intensity as a function of the propagation distance for the pulse energy of 18 pJ. Right: The temporal profile of the pulse intensity taken at the point of the highest intensity, i.e. at the top of the second peak in the respective plots on left.
In order to better visualize the difference between the dynamics over the length of the waveguide, Fig. 5 shows the on-axis intensity temporal profile versus the propagation distance. It is evident that while the behavior is correctly captured in the single-mode approximation, the maxima propagate with slightly different group velocities, manifesting themselves in different slopes of the bright traces in the space-time maps of Fig. 5. The propagation velocities are important in affecting where exactly the dispersive-wave spectral peaks appear in the generated supercontinuum, because they arise thanks to group-velocity matching between the source (i.e. the most intense part of the pulse) and the dispersive radiation. The following condition controls the frequency ω for the dispersive wave generated from a “source” pulse [25,26] of central frequency ωs:
with Vs the wavepacket group velocity, and the intensity-dependent last term reflecting the nonlinear modification of the propagation constant. However, one can see that a difference between VS, as produced by GNLS and/or gUPPE simulation, alone can not explain the different location of the peaks. Indeed, having seen that intensities are similar, we can infer that the last term in Eq. (10) is similar in both approaches. Moreover, the difference in the VS is much smaller in the lower-energy pulse (see similar slopes in Fig. 5 left panel) than in the higher-energy one (Fig. 5 right panel), while Fig. 3 shows that the difference in the dispersive wave peak location already exists at the lower energy. This means that we have to look further for what modifies the dynamics of the dispersive-wave features in the fully resolved simulation.Fig. 5 On-axis intensity vs propagation distance and local time for fully resolved (F) simulation and single mode (S) approximation. These maps show the last 30 mm propagation in the waveguide for 18 pJ pJ (left) and 21.4 pJ (right) pulse energy.
Obviously, the differences in dynamics discussed above are caused by the higher-order modal components that develop during the nonlinear interaction of the pulse with the waveguide. Inspection of the transverse profile of the propagating pulse reveals that deviations from the fundamental mode shape are rather small. Based on the projection onto the fundamental mode, we estimate that in the pulse slices most affected by spatial reshaping, less than ten percent of the power propagates in higher-order modes. Therefore, it is most likely that only several lowest-order modes have significant contributions.
Instead of using a modal decomposition, here we wish to present an alternative way to understand the pulse waveform. We introduce a “spatial envelope,” S(x, y), to describe the transverse profile of the pulse as
where M represents the fundamental waveguide mode at the central wavelength. Figure 6 depicts the spatial envelope for the time slice t corresponding to the maximal on-axis intensity of the pulse. In the initial portion of the waveguide, S is practically constant as the pulse settles into a shape extremely close to that of the fundamental mode. As the intensity around the peak of the pulse increases in the second half of the waveguide, the spatial envelope increases in the on axis region, while it decreases around and outside of the core-cladding interface. This means that the most intense portion of the pulse is more localized than the fundamental mode is. Figure 6 shows that this deformation is almost constant when the intensity of the pulse is high. In other words, the spatial envelope is a natural “slow variable” that describes the trans-verse shape of the pulse. By definition, it is a continuous function across both the TE and TM interfaces, while it only exhibits discontinuity of the normal derivative across the TE interfaces. As a consequence, S could be sampled with a coarser discrete grid. This makes us believe that it can provide an efficient way to represent the pulse in numerical simulations.Fig. 6 Spatial envelope S, Eqn. (11), of the nonlinear pulse propagating along z in a waveguide. The initially constant value represents the propagation in the fundamental mode. In the second half of the waveguide length, narrowing of the profile around the center at x = 0 indicates that the propagating pulse experiences an effectively narrower waveguide due to nonlinear interactions.
Narrowing of the spatial envelope in the high-intensity regimes also suggests a qualitative explanation for the difference between the supercontinuum spectra of the GNLS approximation and the space-resolving simulation. The dispersive-wave spectral peak locations are determined by the group velocity dispersion profile. Roughly speaking, their creation can be viewed as a scattering from the most intense pulse component into a frequency region that is group-velocity matched to the main peak [43, 44]. This is why such a phase-matching condition usually connects wavelengths on the opposite sides of zero-GVD wavelength as indicated by Eq. (10). Qualitatively, the effective narrowing of the pulse profile due to nonlinearity has a similar effect as shrinking of the waveguide, which results in zero-GVD locations being closer to each other, which in turn brings the dispersive-wave peaks also closer to the central wavelength. This is a purely spatial effect not present in the GNLS method. We thus conclude that there are in fact two distinct effects that modify the location of the dispersive-wave peak in the spectrum. One is the change of the effective index due to the Kerr effect, and this is captured by GNLS, and often taken into account in the supercontinuum spectrum analysis. The second, only present in the spatially resolved treatment, is the modal narrowing that modifies the effective dispersion properties for the most intense portion of the pulse and thus contributes an additional shift of the spectral peak(s).
4. Conclusion
We have described a space- and time-resolved method to simulate nonlinear pulse propagation in waveguide structures with strong confinement of the light field. It utilizes the recently derived generalized unidirectional pulse propagation framework to implement an approach capable of capturing extremely wide-bandwidth pulsed waveforms propagating in waveguides over experimentally relevant scales. While the new method is still a few orders of magnitude more computationally demanding than the traditional single transverse mode approximation, we have demonstrated that fully resolved simulations are feasible with rather modest computational resources.
Our approach thus offers a novel alternative for situations in which it is important to capture the possibility of higher-order mode excitation due to nonlinear interactions, and also for situations in which the modal shapes change with the wavelength and thus affect the nonlinearity experienced by different spectral components of the pulse.
We have performed simulations of supercontinuum generation in a waveguide on a chip to validate and compare the new method to the single-mode model based on GNLS. In this context, the single-mode method has been previously shown to produce results in good agreement with experiments. Our comparative simulations reveal that there are subtle spatial effects that can further modify the resulting supercontinuum spectrum. In particular, we have identified a spatial effect that adds to the spectral shift of dispersive-wave peaks. Similar to the intensity-dependent modification of a phase-matching relation that relates the pulse central wavelength to that of the dispersive wave, this effect modifies (narrows) the transverse shape and through it also the effective dispersion of the nonlinearly propagating wavepacket.
It appears that the spatial narrowing of the pulse is “locked” to the most intense portion of the pulse. While the present method does not evaluate contributions of higher modes explicitly, this observation suggests that a portion of the energy carried by the higher-order modes that contribute to this deformation do not propagate with their natural phase or group velocities. Rather, they follow the most intense feature of the pulse. This kind of a behavior is likely to be analogous to the phase “locking” between the fundamental and third-harmonic radiation [45] which, mutually coupled through the Kerr nonlinearity, tend to propagate together. We propose that a natural way to describe this is in terms of what we call a spatial envelope. Our modeling indicates that except at the very beginning of the waveguide, the spatial envelope is a simple function that can be sampled with a relatively small number of grid points. We are currently investigating the possibility of using the spatial envelope as an extremely economic space-resolved representation of the pulsed waveform for simulations in strongly confining waveguides.
Acknowledgments
This work was funded by the United States Air Force Office for Scientific Research under grant number FA9550-11-1-0144. A. B. was supported by grant number FA9550-13-1-0228.
References and links
1. P. S. J. Russell, P. Hölzer, W. Chang, A. Abdolvand, and J. Travers, “Hollow-core photonic crystal fibres for gas-based nonlinear optics,” Nature Photonics 8, 278–286 (2014). [CrossRef]
2. J. C. Travers, W. Chang, J. Nold, N. Y. Joly, and P. St J Russell, “Ultrafast nonlinear optics in gas-filled hollow-core photonic crystal fibers [invited],” JOSA B 28, A11–A26 (2011). [CrossRef]
3. M. Nurhuda, A. Suda, K. Midorikawa, M. Hatayama, and K. Nagasaka, “Propagation dynamics of femtosecond laser pulses in a hollow fiber filled with argon: constant gas pressure versus differential gas pressure,” J. Opt. Soc. Am B 20, 2002–2011 (2003). [CrossRef]
4. P. Hölzer, W. Chang, J. Travers, A. Nazarkin, J. Nold, N. Joly, M. Saleh, F. Biancalana, and P. S. J. Russell, “Femtosecond nonlinear fiber optics in the ionization regime,” Phys. Rev. Lett. 107, 203901 (2011). [CrossRef] [PubMed]
5. M. Nisoli, S. De Silvestri, O. Svelto, R. Szipöcs, K. Ferencz, Ch. Spielmann, S. Sartania, and F. Krausz, “Compression of high-energy laser pulses below 5 fs,” Opt. Lett. 22, 522–524 (1997). [CrossRef] [PubMed]
6. M. A. Foster, A. C. Turner, M. Lipson, and A. L. Gaeta, “Nonlinear optics in photonic nanowires,” Opt. Exp. 16, 1300–1320 (2008). [CrossRef]
7. A. L. Gaeta, “Catastrophic collapse of ultrashort pulses,” Phys. Rev. Lett. 84, 3582–3585 (2000). [CrossRef] [PubMed]
8. A. L. Gaeta, “Nonlinear propagation and continuum generation in microstructured optical fibers,” Opt. Lett. 27, 924–926 (2002). [CrossRef]
9. X. Gai, D.-Y. Choi, S. Madden, Z. Yang, R. Wang, and B. Luther-Davies, “Supercontinuum generation in the mid-infrared from a dispersion-engineered As2S3 glass rib waveguide,” Opt. Lett. 37, 3870–3872 (2012). [CrossRef] [PubMed]
10. S. P. Stark, J. C. Travers, and P. S. J. Russell, “Extreme supercontinuum generation to the deep UV,” Opt. Lett. 37, 770–772 (2012). [CrossRef] [PubMed]
11. T. Popmintchev, M.-C. Chen, D. Popmintchev, P. Arpin, S. Brown, S. Ališauskas, G. Andriukaitis, T. Balčiunas, O. D. Mücke, A. Pugzyls, A. Baltuška, B. Shim, S. E. Schrauth, A. Gaeta, C. Hernández-García, L. Plaja, A. Becker, A. Jaron-Becker, M. M. Murnane, and H. C. Kapteyn, “Bright coherent ultrahigh harmonics in the keV X-ray regime from mid-infrared femtosecond lasers,” Science 336, 1287–1291 (2012). [CrossRef] [PubMed]
12. H. S. Eisenberg, R. Morandotti, Y. Silberberg, S. Bar-Ad, D. Ross, and J. S. Aitchison, “Kerr spatiotemporal self-focusing in a planar glass waveguide,” Phys. Rev. Lett. 87, 043902 (2001). [CrossRef] [PubMed]
13. S. Akturk, C. L. Arnold, B. Zhou, and A. Mysyrowicz, “High-energy ultrashort laser pulse compression in hollow planar waveguides,” Opt. Lett. 34, 1462–1464 (2009). [CrossRef] [PubMed]
14. S. Chen, A. Jarnac, A. Houard, Y. Liu, C. L. Arnold, B. Zhou, B. Forestier, B. Prade, and A. Mysyrowicz, “Compression of high-energy ultrashort laser pulses through an argon-filled tapered planar waveguide,” J. Opt. Soc. Am. B 28, 1009–1012 (2011). [CrossRef]
15. C. L. Arnold, S. Akturk, M. Franco, A. Couairon, and A. Mysyrowicz, “Compression of ultrashort laser pulses in planar hollow waveguides: a stability analysis,” Opt. Exp. 17, 11122–11129 (2009). [CrossRef]
16. L. Caspani, D. Duchesne, K. Dolgaleva, S. J. Wagner, M. Ferrera, L. Razzari, A. Pasquazi, M. Peccianti, D. J. Moss, J. S. Aitchison, and R. Morandotti, “Optical frequency conversion in integrated devices,” J. Opt. Soc. Am. B 28, A67–A82 (2011). [CrossRef]
17. Y. Okawachi, K. Saha, J. S. Levy, Y. H. Wen, M. Lipson, and A. L. Gaeta, “Octave-spanning frequency comb generation in a silicon nitride chip,” Opt. Lett. 36, 3398–3400 (2011). [CrossRef] [PubMed]
18. D. Duchesne, M. Peccianti, M. R. E. Lamont, M. Ferrera, L. Razzari, F. Lgar, R. Morandotti, S. Chu, B. E. Little, and D. J. Moss, “Supercontinuum generation in a high index doped silica glass spiral waveguide,” Opt. Exp. 18, 923–930 (2010). [CrossRef]
19. L. Zhang, Y. Yan, Y. Yue, Q. Lin, O. Painter, R. G. Beausoleil, and A. E. Willner, “On-chip two-octave super-continuum generation by enhancing self-steepening of optical pulses,” Opt. Exp. 19, 11584–11590 (2011). [CrossRef]
20. K. Saha, Y. Okawachi, B. Shim, J. S. Levy, R. Salem, A. R. Johnson, M. A. Foster, M. R. E. Lamont, M. Lipson, and A. L. Gaeta, “Modelocking and femtosecond pulse generation in chip-based frequency combs,” Opt. Exp. 21, 1335–1343 (2013). [CrossRef]
21. R. Halir, Y. Okawachi, J. S. Levy, M. A. Foster, M. Lipson, and A. L. Gaeta, “Ultrabroadband supercontinuum generation in a CMOS-compatible platform,” Opt. Lett. 37, 1685–1687 (2012). [CrossRef] [PubMed]
22. J. M. Dudley, G. Genty, and S. Coen, “Supercontinuum generation in photonic crystal fiber,” Rev. Mod. Phys. 78, 1135–1150 (2006). [CrossRef]
23. G. Genty, P. Kinsler, B. Kibler, and J. M. Dudley, “Nonlinear envelope equation modeling of sub-cycle dynamics and harmonic generation in nonlinear waveguides,” Opt. Exp. 15, 5382–5387 (2007). [CrossRef]
24. E. A. Golovchenko, E. M. Dianov, A. M. Prokhorov, and V. N. Serkin, “Decay of optical solitons,” JETP Lett. 42, 74–77 (1985).
25. A. V. Husakou and J. Herrmann, “Supercontinuum generation of higher-order solitons by fission in photonic crystal fibers,” Phys. Rev. Lett. 87, 203901 (2001). [CrossRef] [PubMed]
26. J. Herrmann, U. Griebner, N. Zhavoronkov, A. Husakou, D. Nickel, J. C. Knight, W. J. Wadsworth, P. S. J. Russell, and G. Korn, “Experimental evidence for supercontinuum generation by fission of higher-order solitons in photonic fibers,” Phys. Rev. Lett. 88, 173901 (2002). [CrossRef] [PubMed]
27. M. Kolesik, E. Wright, and J. Moloney, “Simulation of femtosecond pulse propagation in sub-micron diameter tapered fibers,” Appl. Phys. B 79, 293–300 (2004). [CrossRef]
28. S. Afshar V. and T. M. Monro, “A full vectorial model for pulse propagation in emerging waveguides with subwavelength structures part I: Kerr nonlinearity,” Opt. Exp. 17, 2298–2318 (2009). [CrossRef]
29. C. Courtois, A. Couairon, B. Cros, J. R. Marques, and G. Matthieussent, “Propagation of intense ultrashort laser pulses in a plasma filled capillary tube: Simulations and experiments,” Phys. Plasmas 8, 3445 (2001). [CrossRef]
30. F. Poletti and P. Horak, “Description of ultrashort pulse propagation in multimode optical fibers,” J. Opt. Soc. Am. B 25, 1645–1654 (2008). [CrossRef]
31. F. Poletti and P. Horak, “Dynamics of femtosecond supercontinuum generation in multimode fibers,” Opt. Exp. 17, 6134–6147 (2009). [CrossRef]
32. Yunhong Ding, Jing Xu, Haiyan Ou, and C. Peucheret, “Mode-selective wavelength conversion based on four-wave mixing in a multimode silicon waveguide,” Opt. Exp. 22, 127–135 (2014). [CrossRef]
33. J. Shibayama, M. Muraki, J. Yamauchi, and H. Nakano, “Comparative study of several time-domain methods for optical waveguide analyses,” J. Light. Technol. 23, 2285 (2005). [CrossRef]
34. K. Q. Le, T. Benson, and P. Bienstman, “Application of modified Padé approximant operators to time-domain beam propagation methods,” J. Opt. Soc. Am. B 26, 2285–2289 (2009). [CrossRef]
35. N.-N. Feng, G.-R. Zhou, and W.-P. Huang, “An efficient split-step time-domain beam-propagation method for modeling of optical waveguide devices,” J. Light. Technol. 23, 2186 (2005). [CrossRef]
36. H. M. Masoudi and M. S. Akond, “Time-domain BPM technique for modeling ultra short pulse propagation in dispersive optical structures: Analysis and assessment,” J. Light. Technol. 32, 1936–1943 (2014). [CrossRef]
37. H. M. Masoudi, “A novel nonparaxial time-domain beam-propagation method for modeling ultrashort pulses in optical structures,” J. Light. Technol. 25, 3175–3184 (2007). [CrossRef]
38. H. M. Masoudi and M. S. Akond, “Efficient iterative time-domain beam propagation methods for ultra short pulse propagation: Analysis and assessment,” J. Light. Technol. 29, 2475–2481 (2011). [CrossRef]
39. J. Andreasen and M. Kolesik, “Nonlinear propagation of light in structured media: Generalized unidirectional pulse propagation equations,” Phys. Rev. E 86, 036706 (2012). [CrossRef]
40. J. Andreasen, A. Bahl, and M. Kolesik, “Reflectionless beam propagation on a piecewise linear complex domain,” J. Light. Technol. submitted.
41. A. Couairon, E. Brambilla, T. Corti, D. Majus, O. d. J. Ramírez-Góngora, and M. Kolesik, “Practitioner’s guide to laser pulse propagation models and simulation,” Eur. Phys. J. Special Topics 199, 5–76 (2011). [CrossRef]
42. K. Ikeda, R. E. Saperstein, N. Alic, and Y. Fainman, “Thermal and Kerr nonlinear properties of plasma-deposited silicon nitride/ silicon dioxide waveguides,” Opt. Exp. 16, 12987–12994 (2008). [CrossRef]
43. M. Kolesik, L. Tartara, and J. V. Moloney, “Effective three-wave-mixing picture and first Born approximation for femtosecond supercontinua from microstructured fibers,” Phys. Rev. A 82, 045802 (2010). [CrossRef]
44. D. Faccio, A. Averchi, A. Couairon, M. Kolesik, J. V. Moloney, A. Dubietis, G. Tamosauskas, P. Polesana, A. Piskarskas, and P. DiTrapani, “Spatio-temporal reshaping and X-wave dynamics in optical filaments,” Opt. Exp. 15, 13077–13095 (2007). [CrossRef]
45. M. Kolesik, E.M. Wright, A. Becker, and J.V. Moloney, “Simulation of third-harmonic and supercontinuum generation for femtosecond pulses in air,” Appl. Phys B 85, 531–538 (2006). [CrossRef]
