Numerical solver for supercontinuum generation in multimode optical fibers

Roman Khakimov; Igor Shavrin; Steffen Novotny; Matti Kaivola; Hanne Ludvigsen

doi:10.1364/OE.21.014388

1. Introduction

When short narrow-band laser pulses enter a nonlinear optical medium, they may experience extreme spectral broadening and convert into a broadband spectrally continuous output, a process known as supercontinuum generation (SCG). Since their first introduction, microstructured optical fibers (MOFs) have been particularly favored for generating supercontinua as broad as a few optical octaves. Their unprecedented flexibility in the design of the dispersion characteristics together with the achievable high nonlinearity have boosted the interest in SCG research, both experimentally and theoretically [1].

By now, comparisons between experiments and numerical simulations based on the generalized nonlinear Schrödinger equation (GNLSE) have opened the way towards a detailed physical understanding of the underlying complex dynamics of SCG [2, 3]. As a result, for selected cases, such as the single-mode pulse propagation in silica fibers, a good agreement between experiments and simulations can meanwhile be observed, which, however, can in general not yet be reached for many other situations, among them for instance soft glass fibers [4]. The high nonlinearity of MOFs results from the strong field confinement achieved when these fibers are designed with a small core and a high core-cladding index contrast. Although typically only the fundamental mode is effectively excited, these MOFs can in principle support numerous transverse spatial modes [5, 6]. Their coupling to each other depends on the difference of the propagation constants and can get enhanced when the fiber core diameter is increased. This can give a number of advantages: For instance, pulse propagation along higher-order modes will experience different dispersion properties, which will give access to novel nonlinear optical phenomena in the SCG process [7–11]. Furthermore, the steady demand for an increase in the achievable maximum power spectral density in SCG requires fibers with an increased core size, and thus possible multiple guiding modes, to push the fiber damage threshold beyond current limits.

Despite the fact that experiments have pointed out the importance of intermodal nonlinear effects [8, 12, 13], the complexity of SCG itself and the influence of intermodal phenomena in particular have mostly limited the theoretical studies of the pulse propagation in multimode MOFs to simple cases only. This has left GNLSE modeling constrained by either the number of modes, mostly to two birefringent [14–16] or spatially distinct modes [7, 17], or by the extent to which intermodal nonlinear phenomena have been taken into account [10, 18].

Only recently, the theoretical framework has successfully been extended to enable the study of pulse propagation in multimode fibers taking into account an, in principle, unlimited number of modes [19]. The derived multimode GNLSE (MM-GNLSE) comprises apart from the common nonlinearities, wavelength dependent nonlinear mode couplings and was first applied to study the dynamics of intermodal nonlinear effects in femtosecond SCG [20]. Since then, the importance of this theoretical model has been recognized by a number of researchers investigating into effects, ranging from pulse propagation in non-silica fibers [21] and multicore MOFs [22] to the role of multimode fibers in extending the fiber-capacity limit of modern telecommunications systems [23]. But also the goal to further extend the supercontinuum spectrum into either the UV or infrared, by, for instance, tapering MOFs with high air-fill fraction [24] or using very high NA soft glass fibers [1], has raised the interest in multimode SCG.

The aim of this paper is to present a straightforward approach to solving the MM-GNLSE so that modeling of pulse propagation in multimode fibers may become easily accessible to a broad research community. We transform the MM-GNLSE in the frequency domain to a system of standard ordinary differential equations (ODEs) enabling the use of readily available ODE solvers and thus giving straightforward access to the complex theoretical model. The potential of this numerical approach is demonstrated by simulating the nonlinear propagation of femtosecond pulses in a non-birefringent multimode MOF for a set of linearly polarized eigenmodes.

2. Modeling

We have adopted the MM-GNLSE proposed by Poletti and Horak [19] for numerically describing the ultrashort pulse propagation in the case when multiple transversal modes are excited in an optical fiber. This model extends the GNLSE typically used for describing single-mode propagation [1,2,25] by taking into account, besides Kerr and Raman nonlinearities, self-steepening, polarization and higher-order dispersion effects, also wavelength-dependent intermodal effects for an, in principle, unlimited number of modes. The mathematical formulation consists of a set of equations with their number being equal to the number of modes being considered. Selecting one of the modes (p), the evolution of the pulse envelope A_p(z, T) along the z-axis is described in the co-moving time frame by [19]

\begin{array}{l} \frac{\partial A_{p} (z, T)}{\partial z} = [i (β_{0}^{(p)} - β_{0}) A_{p} (z, T) - (β_{1}^{(p)} - β_{1}) \frac{\partial A_{p} (z, T)}{\partial T} + i \sum_{n ⩾ 2} \frac{β_{n}^{(p)}}{n!} {(i \frac{\partial}{\partial T})}^{n} A_{p} (z, T)] \\ + [i \frac{n_{2} ω_{0}}{c} \sum_{l, m, n} {(1 + i τ_{p l m n}^{(1)} \frac{\partial}{\partial T}) Q_{p l m n}^{(1)} 2 A_{l} (z, T) \int R (T^{'}) A_{m} (z, T - T^{'}) A_{n}^{*} (z, T - T^{'}) d T^{'} \\ + (1 + i τ_{p l m n}^{(2)} \frac{\partial}{\partial T}) Q_{p l m n}^{(2)} A_{l}^{*} (z, T) \int R (T^{'}) A_{m} (z, T - T^{'}) A_{n} (z, T - T^{'}) e^{2 i ω_{0} T^{'}} d T^{'}}] \\ = [𝒟^{(p)} (z, T)] + [𝒩^{(p)} (z, T)] . \end{array}

Generally, Eq. (1) can be divided into a dispersive 𝒟^(p) (first line term of Eq. (1) in square brackets) and a nonlinear part 𝒩^(p). The nonlinear part 𝒩^(p) comprises the couplings between the modes, characterized by the coupling coefficients $Q_{p l m n}^{(1, 2)}$ as defined by [19]

\begin{array}{l} Q_{p l m n}^{(1)} (ω) = \frac{ε_{0}^{2} n_{0}^{2} c^{2}}{12} \frac{\iint [E_{p}^{*} (ω) \cdot E_{l} (ω)] [E_{m} (ω) \cdot E_{n}^{*} (ω)] d S}{N_{p} (ω) N_{l} (ω) N_{m} (ω) N_{n} (ω)} \\ Q_{p l m n}^{(2)} (ω) = \frac{ε_{0}^{2} n_{0}^{2} c^{2}}{12} \frac{\iint [E_{p}^{*} (ω) \cdot E_{l}^{*} (ω)] [E_{m} (ω) \cdot E_{n} (ω)] d S}{N_{p} (ω) N_{l} (ω) N_{m} (ω) N_{n} (ω)} \end{array}

with the normalization coefficients

N_{p}^{2} (ω) = \frac{1}{4} \iint [E_{p}^{*} (ω) \times H_{p} (ω) + E_{p} (ω) \times H_{p}^{*} (ω)] \cdot e_{z} d S .

In Eq. (1) we have freely chosen ω₀, β₀ = β(ω₀) and β₁ = (∂β/∂ω)_ω=ω₀ to be the wavenumber, the propagation constant and the inverse group velocity of the fundamental mode (p = 1), while for the higher-order modes we have

β_{n}^{(p)} = {(\partial^{n} β^{(p)} / \partial ω^{n})}_{ω = ω_{0}}

. Fiber losses are neglected due to the short fiber lengths that will be considered here. The parameters n₀, n₂, ε₀ and c are the linear and the nonlinear refractive indices of the fiber material, vacuum permittivity and speed of light in vacuum, respectively. The delayed Raman response R(T) is defined as in Eq. (9) of [19] and

τ_{p l m n}^{(1, 2)}

is the shock term [19]. The required propagation constants β^(p) and the transverse electric and magnetic mode-field profiles E_p and H_p of each mode are here obtained independently from finite-element method calculations.

Note that through $Q_{p l m n}^{(1, 2)}$ in Eq. (2) this description allows us to account for the frequency dependence of the mode profiles across the SC bandwidth, analogously to the spectrally dependent effective area in the single-mode model for SCG [25]. Also, arbitrary profiles of the mode fields can be taken into account so that both circularly and linearly polarized modes can be studied with this model.

2.1. Numerical solving

In general, the GNLSE has no analytic solution and thus has to be solved numerically. The most widely employed numerical method is the split-step Fourier method as it provides a good compromise between efficiency, accuracy and coding simplicity. However, in order to guarantee the accuracy and convergence of the results, the numerical step sizes for the position coordinate have to be selected carefully as otherwise multiple simulations with decreasing step sizes or other means of error control and adaptive step-size algorithms will be required [26]. These limitations were removed by transforming the equation into frequency domain and integrating it by the well established and reliable fourth-order Runge-Kutta method, as demonstrated earlier for solving the single-mode GNLSE [1, 27, 28].

Here, we have adopted this idea to solve the multi-mode GNLSE. Similarly, the MM-GNLSE is transformed in frequency domain to a system of first-order ordinary differential equations (ODEs). We then utilized the fourth-order Runge-Kutta ODE solver provided by the MATLAB^® package to solve the system of ODE by direct integration. Implementation of the frequency domain MM-GNLSE solver thus becomes fast and customizable.

In order to rewrite the MM-GNLSE in the form of a system of first-order ODEs,

\frac{d y (z)}{d z} = f (y, z),

we start by first expressing the temporal pulse shape A_p(z, T) by the inverse Fourier transform ℱ⁻¹ of its complex spectral envelope Ã_p(z, ω − ω₀),

A_{p} (z, T) = ℱ^{- 1} {{\tilde{A}}_{p} (z, ω - ω_{0})} = \frac{1}{2 π} \int {\tilde{A}}_{p} (z, ω - ω_{0}) exp [- i (ω - ω_{0}) T] d ω .

Substituting A_p(z, T) by Eq. (5) only in the dispersive part of Eq. (1), and subsequently applying the Fourier transform to both sides of Eq. (1) yields

\frac{\partial {\tilde{A}}_{p} (z, ω - ω_{0})}{\partial z} - i [β^{(p)} (ω) - β_{0} - (ω - ω_{0}) β_{1}] {\tilde{A}}_{p} (z, ω - ω_{0}) = {\tilde{𝒩}}^{(p)} (z, ω - ω_{0}),

where 𝒩̃^(p) = ℱ𝒩^(p). Next we use the relation

(1 + i τ_{p l m n} \frac{\partial}{\partial T}) = ℱ^{- 1} (1 + (ω - ω_{0}) τ_{p l m n}) ℱ

and approximate the shock term by τ_plmn ≈ 1/ω₀ as discussed in [19]. Following [1], we perform a change of variables

{\tilde{A}}_{p}^{'} = {\tilde{A}}_{p} exp [- L^{(p)} (ω) z], L^{(p)} (ω) = i [β^{(p)} (ω) - β_{0} - (ω - ω_{0}) β_{1}] .

and finally arrive at the ODE describing the pulse propagation for mode p

\begin{array}{l} \frac{\partial {\tilde{A}}_{p}^{'} (z, ω - ω_{0})}{\partial z} = exp [- L^{(p)} (ω) z] i \frac{n_{2} ω}{c} \times \\ \times \sum_{l, m, n}^{M} ℱ {Q_{p l m n}^{(1)} 2 A_{l} (z, T) [R * (A_{m} A_{n}^{*})] + Q_{p l m n}^{(2)} A_{l}^{*} (z, T) [(R \cdot e^{2 i ω_{0} T}) * (A_{m} A_{n})]} . \end{array}

with * describing the convolution integral and superscript ^* the complex conjugate. For a total number of M modes and a frequency grid with N frequency points, the system of equations forms a first-order vector ODE (Eq. (4)). Comparing Eq. (8) with Eq. (4), we find, that the right hand side of Eq. (8) provides the entries of f(y, z) and

y (z) = (\begin{matrix} {\tilde{A}}^{'}_{1} \\ ⋮ \\ {\tilde{A}}^{'}_{M} \end{matrix}), {\tilde{A}}^{'}_{p} = (\begin{matrix} {\tilde{A}}_{p}^{'} (z, ω_{1} - ω_{0}) \\ ⋮ \\ {\tilde{A}}_{p}^{'} (z, ω_{N} - ω_{0}) \end{matrix})

with initial conditions

y (0) = (\begin{matrix} {\tilde{A}}_{1}^{'} (0, ω_{1} - ω_{0}) \\ ⋮ \\ {\tilde{A}}_{M}^{'} (0, ω_{N} - ω_{0}) \end{matrix}) .

Solving the first-order vector ODE, consisting of MN equations in total, yields Ã′_p. The spectral and temporal evolution of each mode is then obtained by performing a backward change of variables according to Eq. (7).

The MM-GNLSE solver presented here has been verified to reproduce the previously published results of femtosecond SCG simulations in multimode fibers [20], which were obtained by numerically solving the MM-GNLSEs using the split-step Fourier method.

3. Simulation results

In this section we demonstrate the potential of the theoretical model and the MM-GNLSE solver by simulating the SCG of ultrashort pulses launched into a multimode, non-birefringent MOF.

The design of the chosen silica MOF is described by a triangular grid of holes with a hole-to-hole pitch of Λ = 2.1μm and a hole-pitch ratio of d/Λ = 0.89 (see inset of Fig. 1). A full vectorial finite-element method (FEM) analysis utilizing COMSOL Multiphysics^® software yields as a result that this fiber geometry is non-birefringent and supports more than 6 guided modes at wavelengths below 1.7 μm. From the analysis we extract for the six lowest order modes (labeled M1 - M6) the electric E_p(ω) and magnetic H_p(ω) fields as well as the propagation constants β_p(ω). The respective group-velocity dispersion (GVD) parameters, obtained from β_p(ω), are plotted in Fig. 1(a) in the wavelength range from 0.4 μm to 1.5 μm. For modes M1 and M2 the zero GVD wavelengths (ZDW) are located near 0.8 μm, whereas the ZDWs of the modes M3 to M6 are shifted to around 0.6 μm. Within the numerical accuracy the field distributions of the modes are found to be transversal only. The mode fields and polarizations of modes M1 and M2 are shown in Fig. 1(b). We have calculated $Q_{p l m n}^{(1, 2)}$ only at the pump wavelength (ω = ω₀) in Eq. (8) to significantly shorten the computation time for $Q_{p l m n}^{(1, 2)}$ . As the rate of intermodal energy transfer is lower outside the narrow-band region of the pump, we will neglect by this simplification only small variations in the coupling efficiency and mode intensity profiles over the whole spectrum.

Figure 1 (a) Calculated GVD parameter for the six lowest-order modes of the MOF. Inset shows the fiber geometry used for the eigenmode analysis. (b) Mode fields (upper) and polarization (lower) of the two lowest-order modes M1 and M2.

Download Full Size | PDF

The generation of the supercontinuum is first simulated along a 15 cm long piece of the described MOF for 150 fs-long input pulses with a peak power of 10 kW. The chirp-free pulses are of secant-hyperbolic shape with a central wavelength at 0.8 μm. Initially, the light is linearly polarized and assumed to be fully coupled into mode M1. However, all modes require the inclusion of noise as otherwise no interaction between the modes is possible. This is an essential difference compared to single-mode studies and we thus have to take into account not only the pump laser noise but also the intrinsic quantum fluctuations present in each mode. This is realized here by exciting all six modes through quantum noise by including one photon per mode with a random phase in each spectral component, following [14]. Alternatively, depending on the pump pulse properties, a more elaborate model for the pump laser noise can be easily included to enable more detailed comparisons with experiments [29].

Solving the system of equations (Eq. (8)), we obtain the spectral contents of each mode at the end of the fiber as displayed in Fig. 2. We observe that energy has effectively only been transferred from the initially excited mode M1 to mode M2, with negligible amounts of energy in modes M3 to M6. Consequently, under these conditions the simulation reduces to the most simple multimode case when light is effectively only coupled between the two orthogonal polarization states of the fundamental mode (see Fig. 1(b)).

Figure 2 Simulated spectral contents of modes M1 to M6 after 15 cm of propagation with the initial pump pulse fully coupled to mode M1.

Download Full Size | PDF

Since the effective indices and mode intensity profiles of modes M1 and M2 match each other within the accuracy of the FEM analysis, the spectral output of the fiber is expected to be the same independent of in which of the two modes the light is propagating. Indeed, the spectra presented in Fig. 3 for an ensemble average of 100 simulations show good agreement between the combined spectra of modes M1 and M2 and simulations with light propagation in mode M1 only, equivalent to the scalar-GNLSE result. This confirms the correctness of the MM-GNLSE simulation. However, while the total spectral information can readily be obtained from solving the scalar-GNLSE considering mode M1 only, the actual complex interaction between modes M1 and M2 remains hidden in this simplified analysis.

Figure 3 Ensemble average of supercontinuum spectra obtained from 100 simulations with pulse propagation either restricted to mode M1 only (red line) or including coupling between modes M1 and M2 (blue line) for input pulses with 10 kW peak power fully coupled to mode M1.

Download Full Size | PDF

Figure 4 presents the relative mode energy ratios for different input pulse peak powers up to a fiber length of 30 cm obtained from MM-GNLSE simulations. The results originate from single-shot simulations by integrating the pulse envelope for each mode at position z and normalizing it to the input pulse energy (z = 0). We verified that in the absence of the Raman effect the total energy remains conserved during the propagation of the light pulses. The presence of the Raman effect, however, results in an energy loss during propagation, as represented by ∑ in Fig. 4. Further we confirmed that also the photon number $N_{p h} \propto \sum_{p = 1}^{M} \int {| A_{p} (ω) |}^{2} / \bar{h} ω d ω$ remains conserved. The energy transfer among the two modes strongly depends on the initial pulse peak power. While at peak powers of 1 kW only a negligible amount of energy is transferred to mode M2 by the end of the fiber, peak powers above 10 kW lead already to an energy transfer of more than ≈ 35 %. Also the onset of the energy transfer starts earlier with higher peak powers, for instance already at ≈ 1 cm for peak powers of 50 kW.

Figure 4 Evolution of energy distribution among modes M1 and M2 for 1 kW (solid line), 10 kW (dotted line) and 50 kW (dash-dotted line) peak pump powers and up to a propagation length of 30 cm. The respective gray lines represent the total energy ∑ of the two modes for each launched pulse. Input power is originally coupled to mode M1 only.

Download Full Size | PDF

Selecting 10 kW peak pump power and considering only the first 15 cm of the fiber, the spectral and temporal evolution of modes M1 and M2 are presented in Fig. 5. Initially, the pulse launched into mode M1 forms a higher order soliton. The soliton is then temporally compressed until soliton fission sets in at about 3 cm. As a result, a broad spectrum is formed which propagates in mode M1 with an almost unchanged width to the end of the fiber. The onset of light propagation in mode M2 can be observed after about 6 cm. By this point cross-phase modulation becomes very efficient due to a complete phase synchronization between modes M1 and M2, and results in significant energy transfer to mode M2. The power spectral density in mode M2 increases towards the end of the fiber, even reaching within a narrow spectral range a power level comparable to mode M1.

Figure 5 Spectral and temporal evolution of the supercontinuum in modes M1 and M2 generated from 10 kW pulses originally launched into mode M1.

Download Full Size | PDF

Although in the example presented here SCG is initiated from a fully linearly polarized pulse launched into mode M1, coupling to the orthogonally polarized mode M2 determines the overall polarization properties of the supercontinuum at the fiber output. In the past, the polarization properties of supercontinua generated in weakly and highly birefringent MOFs have been studied both experimentally and numerically [14–16,30–32]. Thereby, the numerical models were exclusively based on solving the coupled GNLSE with two circularly polarized eigenmodes [14, 16, 32]. Using the MM-GNLSE allows us to study the nonlinear polarization dynamics in SCG independent of the eigenmode basis.

To analyze the polarization properties further we utilize the partial polarization theory for pulsed optical beams introduced in [33] for the ensemble of simulated supercontinua presented earlier in Fig. 3. Following [33] we determine the spectral degree of polarization P_s(z, ω),

P_{s} (z, ω) = \sqrt{2 \frac{tr Φ^{2} (z, ω)}{{tr}^{2} Φ (z, ω)} - 1},

where at propagation distance z and frequency ω the spectral polarization matrix Φ(z, ω) has the following elements

Φ_{i j} (z, ω) = 〈 {\tilde{A}}_{i}^{*} (z, ω) {\tilde{A}}_{j} (z, ω) 〉,

which are obtained from the calculated complex spectral envelopes Ã_i(z, ω) with i, j = x, y. Here we assume x to refer to M1 and y to M2. By definition the spectral degree of polarization is limited between 0 ≤ P_s(z, ω) ≤ 1, indicating at its upper bound a spectrally fully polarized field for the ensemble average, whereas the lower bound corresponds to an unpolarized field. For each spectral component, the evolution of P_s along the fiber length is shown in Fig. 6(a). We find that the generated supercontinuum remains, as it should, spectrally fully polarized until the onset of mode M2. From this point on the spectral degree of polarization rapidly decreases at those wavelengths where mode M2 contributes to the total supercontinuum. On the other hand, a comparison of the degree of spectral polarization at z = 15cm (Fig. 6(b)) with the ensemble average of the spectral contents in each mode (Fig. 6(d)) reveals that only those parts of the spectrum where the contribution by mode M2 is negligible will remain fully polarized up to the end of the fiber.

Figure 6 (a) Evolution of the spectral degree of polarization P_s (Eq. (11)) along the fiber. (b) Spectral degree of polarization as a function of wavelength at the end of the fiber. (c) ( Media 1) Graphical illustration of the evolution of the polarization state at 0.8 μm along the fiber. The tip of the Poincaré vector is traced by the purple line. The green line on the Poincaré sphere indicates the corresponding polarization state of the polarized part of the light. (d) Ensemble average of the spectral contents in mode M1 (red line) and M2 (blue line).

Download Full Size | PDF

Furthermore, from the spectral polarization matrix Φ(z, ω) we can also determine the state of polarization for each spectral component by using the spectral Stokes parameters as given in [33]. Normalized to the spectral density of the field, the Stokes parameters describe the Poincaré vector which allows to geometrically illustrate the polarization state on the Poincaré sphere with unit radius. The length of the Poincaré vector is thereby identical to the degree of polarization [33], i.e., the state at the origin is unpolarized whereas states at the surface are fully polarized. By extending the Poincaré vector of a partially polarized state to the surface of the sphere, we may indicate the polarization state of its polarized part.

We studied the evolution of the polarization state along the fiber exemplarily at three different wavelengths: 0.7 μm, 0.8 μm and 0.9 μm. The results are illustrated in Fig. 6(c) ( Media 1) and Fig. 7 (( Media 2)and ( Media 3)), where the trace of the tip of the Poincaré vector is represented by the purple line and the corresponding state of polarization of the polarized part is represented by the green line on the Poincaré sphere. At 0.8 μm (Fig. 6(c), ( Media 1)) the light is initially fully linearly polarized so that the tip of the Poincaré vector is located on the sphere at the coordinates (S_s₁ = 1, S_s₂ = 0, S_s₃ = 0). Once also the orthogonally polarized mode M2 starts to develop, not only does the degree of polarization decrease, as already seen above, but also the polarization state changes. Light at 0.8 μm becomes elliptically polarized with the orientation and shape of the polarization ellipse varying along the fiber.

Figure 7 Graphical illustration of the evolution of the polarization state at wavelengths (a) ( Media 2) 0.7 μm and (b) ( Media 3) 0.9 μm along the fiber. For further explanations see Fig. 6(c).

Download Full Size | PDF

At 0.7 μm (Fig. 7(a), ( Media 2) and 0.9 μm (Fig. 7(b), ( Media 3), where the spectrum is initiated from noise, the initial degree of polarization is essentially zero. Consequently, the initial state of polarization is random but quickly evolves to a fully linearly polarized state with the onset of spectral evolution in mode M1 at 0.7 μm and 0.9 μm. The degree of spectral polarization reduces then again when the spectral components appear in mode M2, with the state of polarization changing to elliptical.

4. Conclusions

We have demonstrated that the set of equations forming the MM-GNLSE described in [19] can readily be solved in the frequency domain by transformation to a system of first-order ordinary differential equations. The complete temporal and spectral information of ultra-short pulses propagating in a multimode fiber can then be extracted in a straightforward manner by employing standard ODE solvers. We are convinced that this approach will make the accessibility easier and stimulate modeling of pulse propagation in multimode fibers within a broader research community.

We have verified the solver both by reproducing previously published results on SCG in multimode fibers [20] and by studying the simplest multimode case when energy is effectively only transferred between the two orthogonal states of linear polarization of the fundamental mode in a non-birefringent microstructured optical fiber. As one example among many applications, we have demonstrated for the latter case how modeling based on the MM-GNLSE may enable gaining deep insights into the complex nonlinear polarization dynamics in SCG.

Acknowledgments

This work has been financially supported by the Academy of Finland as part of the “Photonics and Modern Imaging Techniques” research programme (project 134857). RK and IS wish to thank the Graduate School of Modern Optics and Photonics. The computational part of this work has been made possible by the use of computational resources provided by the Aalto Science-IT project. Authors appreciate discussions with Dr. P. Horak.

References

1. J. M. Dudley and J. R. Taylor, eds., Supercontinuum Generation in Optical Fibers (Cambridge University Press, New York, USA, 2010) [CrossRef] .

2. J. M. Dudley, G. Genty, and S. Coen, “Supercontinuum generation in photonic crystal fiber,” Rev. Mod. Phys. 78, 1135–1184 (2006) [CrossRef] .

3. J. C. Travers, “Blue extension of optical fibre supercontinuum generation,” J. Opt. 12, 113001 (2010) [CrossRef] .

4. C. Agger, C. Petersen, S. Dupont, H. Steffensen, J. K. Lyngsø, C. L. Thomsen, J. Thøgersen, S. R. Keiding, and O. Bang, “Supercontinuum generation in ZBLAN fibers—detailed comparison between measurement and simulation,” J. Opt. Soc. Am. B 29, 635–645 (2012) [CrossRef] .

5. J. K. Ranka, R. S. Windeler, and A. J. Stentz, “Optical properties of high-delta air-silica microstructure optical fibers,” Opt. Lett. 25, 796–798 (2000) [CrossRef] .

6. R. Guobin, W. Zhi, L. Shuqin, and J. Shuisheng, “Mode classification and degeneracy in photonic crystal fibers,” Opt. Express 11, 1310–1321 (2003) [CrossRef] [PubMed] .

7. J. M. Dudley, L. Provino, N. Grossard, H. Maillotte, R. S. Windeler, B. J. Eggleton, and S. Coen, “Supercontinuum generation in air-silica microstructured fibers with nanosecond and femtosecond pulse pumping,” J. Opt. Soc. Am. B 19, 765–771 (2002) [CrossRef] .

8. A. Efimov, A. J. Taylor, F. G. Omenetto, J. C. Knight, W. J. Wadsworth, and P. St. J. Russell, “Nonlinear generation of very high-order UV modes in microstructured fibers,” Opt. Express 11, 910–918 (2003) [CrossRef] [PubMed] .

9. S. O. Konorov, E. E. Serebryannikov, A. M. Zheltikov, P. Zhou, A. P. Tarasevitch, and D. von der Linde, “Mode-controlled colors from microstructure fibers,” Opt. Express 12, 730–735 (2004) [CrossRef] [PubMed] .

10. R. Cherif, M. Zghal, L. Tartara, and V. Degiorgio, “Supercontinuum generation by higher-order mode excitation in a photonic crystal fiber,” Opt. Express 16, 2147–2152 (2008) [CrossRef] [PubMed] .

11. B. Zwan, S. Legge, J. Holdsworth, and B. King, “Spatio-spectral analysis of supercontinuum generation in higher order electromagnetic modes of photonic crystal fiber,” Opt. Express 21, 834–839 (2013) [CrossRef] [PubMed] .

12. F. G. Omenetto, A. J. Taylor, M. D. Moores, J. Arriaga, J. C. Knight, W. J. Wadsworth, and P. St. J. Russell, “Simultaneous generation of spectrally distinct third harmonics in a photonic crystal fiber,” Opt. Lett. 26, 1158–1160 (2001) [CrossRef] .

13. C. Lesvigne, V. Couderc, A. Tonello, P. Leproux, A. Barthélémy, S. Lacroix, F. Druon, P. Blandin, M. Hanna, and P. Georges, “Visible supercontinuum generation controlled by intermodal four-wave mixing in microstructured fiber,” Opt. Lett. 32, 2173–2175 (2007) [CrossRef] [PubMed] .

14. S. Coen, A. H. L. Chau, R. Leonhardt, J. D. Harvey, J. C. Knight, W. J. Wadsworth, and P. St. J. Russell, “Supercontinuum generation by stimulated Raman scattering and parametric four-wave mixing in photonic crystal fibers,” J. Opt. Soc. Am. B 19, 753–764 (2002) [CrossRef] .

15. M. Lehtonen, G. Genty, H. Ludvigsen, and M. Kaivola, “Supercontinuum generation in a highly birefringent microstructured fiber,” Appl. Phys. Lett. 82, 2197–2199 (2003) [CrossRef] .

16. Z. Zhu and T. G. Brown, “Polarization properties of supercontinuum spectra generated in birefringent photonic crystal fibers,” J. Opt. Soc. Am. B 21, 249–257 (2004) [CrossRef] .

17. D. A. Akimov, E. E. Serebryannikov, A. M. Zheltikov, M. Schmitt, R. Maksimenka, W. Kiefer, K. V. Dukel’skii, V. S. Shevandin, and Y. N. Kondrat’ev, “Efficient anti-Stokes generation through phase-matched four-wave mixing in higher-order modes of a microstructure fiber,” Opt. Lett. 28, 1948–1950 (2003) [CrossRef] [PubMed] .

18. T. Chaipiboonwong, P. Horak, J. D. Mills, and W. S. Brocklesby, “Numerical study of nonlinear interactions in a multimode waveguide,” Opt. Express 15, 9040–9047 (2007) [CrossRef] [PubMed] .

19. F. Poletti and P. Horak, “Description of ultrashort pulse propagation in multimode optical fibers,” J. Opt. Soc. Am. B 25, 1645–1654 (2008) [CrossRef] .

20. F. Poletti and P. Horak, “Dynamics of femtosecond supercontinuum generation in multimode fibers,” Opt. Express 17, 6134–6147 (2009) [CrossRef] [PubMed] .

21. J. H. V. Price, X. Feng, A. M. Heidt, G. Brambilla, P. Horak, F. Poletti, G. Ponzo, P. Petropoulos, M. Petrovich, J. Shi, M. Ibsen, W. H. Loh, H. N. Rutt, and D. J. Richardson, “Supercontinuum generation in non-silica fibers,” Opt. Fiber Technol. 18, 327–344 (2012) [CrossRef] .

22. X.-h. Fang, M.-l. Hu, L.-l. Huang, L. Chai, N.-l. Dai, J.-y. Li, A. Y. Tashchilina, A. M. Zheltikov, and C.-y. Wang, “Multiwatt octave-spanning supercontinuum generation in multicore photonic-crystal fiber,” Opt. Lett. 37, 2292–2294 (2012) [CrossRef] [PubMed] .

23. S. Mumtaz, R.-J. Essiambre, and G. P. Agrawal, “Nonlinear propagation in multimode and multicore fibers: Generalization of the Manakov equations,” J. Lightwave Technol. 31, 398–406 (2013) [CrossRef] .

24. U. Møller, S. T. Sørensen, C. Larsen, P. M. Moselund, C. Jakobsen, J. Johansen, C. L. Thomsen, and O. Bang, “Optimum PCF tapers for blue-enhanced supercontinuum sources,” Opt. Fiber Technology 18, 304–314 (2012). Fiber Supercontinuum sources and their applications [CrossRef] .

25. J. Laegsgaard, “Mode profile dispersion in the generalised nonlinear Schrödinger equation,” Opt. Express 15, 16110–16123 (2007) [CrossRef] [PubMed] .

26. O. V. Sinkin, R. Holzlohner, J. Zweck, and C. R. Menyuk, “Optimization of the split-step Fourier method in modeling optical-fiber communications systems,” J. Lightwave Technol. 21, 61–68 (2003) [CrossRef] .

27. J. Hult, “A fourth-order Runge–Kutta in the interaction picture method for simulating supercontinuum generation in optical fibers,” J. Lightwave Technol. 25, 3770–3775 (2007) [CrossRef] .

28. A. Heidt, “Efficient adaptive step size method for the simulation of supercontinuum generation in optical fibers,” J. Lightwave Technol. 27, 3984–3991 (2009) [CrossRef] .

29. M. H. Frosz, “Validation of input-noise model for simulations of supercontinuum generation and rogue waves,” Opt. Express 18, 14778–14787 (2010) [CrossRef] [PubMed] .

30. A. Apolonski, B. Povazay, A. Unterhuber, W. Drexler, W. J. Wadsworth, J. C. Knight, and P. St. J. Russell, “Spectral shaping of supercontinuum in a cobweb photonic-crystal fiber with sub-20-fs pulses,” J. Opt. Soc. Am. B 19, 2165–2170 (2002) [CrossRef] .

31. Z. Zhu and T. G. Brown, “Experimental studies of polarization properties of supercontinua generated in a birefringent photonic crystal fiber,” Opt. Express 12, 791–796 (2004) [CrossRef] [PubMed] .

32. M. Tianprateep, J. Tada, and F. Kannari, “Influence of polarization and pulse shape of femtosecond initial laser pulses on spectral broadening in microstructure fibers,” Opt. Rev. 12, 179–189 (2005) [CrossRef] .

33. T. Voipio, T. Setälä, and A. T. Friberg, “Partial polarization theory of pulsed optical beams,” J. Opt. Soc. Am. A 30, 71–81 (2013) [CrossRef] .

Numerical solver for supercontinuum generation in multimode optical fibers

Abstract

1. Introduction

2. Modeling

2.1. Numerical solving

3. Simulation results

4. Conclusions

Acknowledgments

References

Supplementary Material (3)

Cited By

Figures (7)

Equations (12)

Optics Express