Mode profile dispersion in the generalized nonlinear Schrödinger equation

J. Lægsgaard

doi:10.1364/OE.15.016110

1. Introduction

The generalized nonlinear Schrödinger equation (GNLSE) has become a standard tool for simulating the propagation of strong light pulses in optical fibers. The GNLSE has been found to be highly useful for making both qualitative and quantitative predictions of nonlinear experiments in standard optical fibers, and has also been succesfully used in the description of highly complex phenomena such as supercontinuum generation in small-core photonic crystal fibers (PCFs) [1]. Recently, a number of experiments involving nonlinear propagation in both solid-core [2, 3] and hollow-core [4, 5, 6, 7] photonic bandgap (PBG) fibers have been reported, and also in these cases simulations based on the GNLSE have been instrumental in interpreting the observed phenomena. Notwithstanding these successes, the advent of PBG fibers calls for a careful reexamination of the GNLSE derivation, since the PBG fibers differ from standard optical fibers in several respects. Most importantly, they guide light in finite frequency windows, and show a frequency dependence of the mode profiles which is much stronger than what is typically found in both standard fibers and PCFs not based on the photonic bandgap effect. Furthermore, for hollow-core PBG fibers one has complications arising from the hybrid nature of the nonlinearity, which has non-neglibible contributions from both silica and air, and from the existence of surface modes surrounding the core, with much higher nonlinear coefficients than the core modes [8, 6].

The purpose of this paper is to address the first of the questions raised above, namely the impact of mode profile dispersion on nonlinear propagation in single-mode PBG fibers. In the GNLSE, the transverse degrees of freedom have been integrated out, and the mode profile enters the calculation through the effective area parameter, A_eff. The effective area contributes to the nonlinear coefficient γ, which is the prefactor of the nonlinear term in the GNLSE, as follows:

γ = \frac{3 χ_{xxxx}^{(3)} ω^{2}}{8 c^{2} β (ω) A_{eff}}

where β (ω) is the propagation constant of the guided fiber mode, ω is the frequency, and χ⁽³⁾_xxxx is the diagonal element of the third-order susceptibility tensor of the fiber material, here assumed to be amorphous silica. The most common approach to treating mode profile dispersion is currently to include a frequency dependence of A_eff in Eq. (1) [9, 10, 11]. However, several authors have pointed out that, at least in principle, one needs to go beyond this approach [12, 13, 14, 15], and that a correct description of mode profile dispersion in fact requires a generalization of the effective-area concept. In this paper, the earlier results are rederived, a systematic approximation scheme is developed for numerical solution of the resulting equations, and the accuracy of a modified GNLSE relative to the full calculation is tested. The main conclusion is that a slight reformulation of the GNLSE [12], which can be formulated in terms of a frequency-dependent A_eff, and has the same numerical complexity as the traditional GNLSE, can greatly enhance the agreement with the more complete theory.

The rest of the paper is structured as follows: In section 2, the formal theory of nonlinear propagation in a single-mode fiber in the presence of mode profile dispersion is derived. In section 3, practical numerical approaches to the problem are developed, while in section 4, results of some numerical tests are presented. Section 5 summarizes the conclusions.

2. Derivation of a 1+1D propagation equation

The starting point for the GNLSE derivation is the Maxwell equations for the electromagnetic field, in the presence of nonlinear polarization:

\nabla \times H = ε (r_{⊥}) \frac{\partial E}{\partial t} + \frac{\partial P_{NL}}{\partial t}, \nabla \times E = - μ_{0} \frac{\partial H}{\partial t}

Here ε(r⊥) is the transverse distribution of the dielectric constant defining the waveguide. The fields are partitioned as:

E (r, t) = \sum_{m} \frac{1}{2 π} \int d ω {\tilde{G}}_{m} (z, ω) e_{m} (r_{⊥}, ω) \exp [i (β_{m} (ω) z - ω t)],

H (r, t) = \sum_{m} \frac{1}{2 π} \int d ω {\tilde{G}}_{m} (z, ω) h_{m} (r_{⊥}, ω) \exp [i (β_{m} (ω) z - ω t)],

where the modal fields e m and h _m are assumed to be eigenstates of the linear eigenvalue problem for a z-propagating mode, with propagation constants β_m(ω). In Ref. [13] this ansatz was shown to lead to the 1+1D propagation equation:

\frac{\partial {\tilde{G}}_{m} (z, ω)}{\partial z} = \frac{2 π i ω}{2 N_{m} (ω)} \int d r_{⊥} \exp (- i β_{m} (ω) z) e_{m}^{*} (r_{⊥}, ω) \cdot P_{NL} (r, ω)

P_{NL} (r, ω) = \int d t \exp (i ω t) P_{NL} (r, t)

where N_m(ω) is a normalization parameter given by:

\int d r_{⊥} [e_{m} \times h_{n}^{*} + e_{n}^{*} \times h_{m}] = 2 N_{m} (ω) δ_{mn}

The nonlinear polarization in the time domain can be written as:

P_{NL} (r, t) = ε_{0} \int d t' R (t - t') χ^{(3)} ⋮ E (r, t) E (r, t') E (r, t')

where R(t) is the Raman response function and χ⁽³⁾ is the full third-order susceptibility tensor. Both R and χ ⁽³⁾ will in this work be assumed independent of position. It is important to note that, at this stage, the field E(r,t) is purely real, which means that the frequency integrals in Eqs. (3), (4) extend over both positive and negative frequencies, with β_m(-ω)=-β_m(ω), e _m(r⊥,-ω)=e*_m(r⊥,ω) etc. Inserting the frequency-domain expansion of E(r,t), Eq. (3), one obtains:

P_{NL} (r, ω) = \frac{ε_{0}}{{(2 π)}^{2}} \sum_{nqp} \int d ω_{1} d ω_{2} G_{n} (z, ω_{1}) G_{q}^{*} (z, ω_{1} + ω_{2} - ω) G_{p} (z, ω_{2}) \times

where the definition

G_{m} (z, ω) = {\tilde{G}}_{m} (z, ω) \exp (- i β_{m} (ω) z)

was used to simplify the final result. Inserting in Eq. (5) yields:

\frac{\partial {\tilde{G}}_{m} (z, ω)}{\partial z} = \frac{i ω ε_{0}}{4 π N_{m} (ω)} \exp (- i β_{m} (ω) z) \sum_{nqp} \int d ω_{1} d ω_{2} G_{n} (z, ω_{1}) \times

K_{mnqp} (ω, ω_{1}, ω_{1} + ω_{2} - ω, ω_{2}) = \int d r_{⊥} e_{m}^{*} (r_{⊥}, ω) \cdot χ^{(3)} ⋮ e_{n} (r_{⊥}, ω_{1}) e_{q}^{*} (r_{⊥}, ω_{1} + ω_{2} - ω) e_{p} (r_{⊥}, ω_{2})

Eq. (11) is fairly general, and useful for describing both scalar and vector nonlinear effects (e.g. polarization rotation). Since the focus of this work is on the treatment of mode profile dispersion, Eq. (11) will in the following be simplified to avoid unnessecary complications. An all-silica PBG fiber with low index contrast will be studied, which allows use of the scalar approximation. It will further be assumed that the light is launched into a single mode with linear polarization. Under these assumptions, only a single term is retained in the sum over modal fields, and this field can be written as

e (r_{⊥}, ω) = F (r_{⊥}, ω) \hat{x}, h (r_{⊥}, ω) = n_{eff} (ω) ε_{0} c F (r_{⊥}, ω) \hat{y}

where n_eff=cβ/ω, x̂, ŷ are unit vectors, and F is a real function. In the example PBG fiber, full-vectorial control calculations showed that the fraction of field energy carried by the neglected vector components was at most on the order of 10^-3. For the following discussion it is useful to introduce the normalizations:

N (ω) = n_{eff} (ω) ε_{0} c, \tilde{G} (z, ω) \to \sqrt{2 n_{0} ε_{0} c} \tilde{G} (z, ω)

\int d r_{⊥} {∣ F (r_{⊥}, ω) ∣}^{2} = 1, E_{p} (z) = \frac{1}{2 π} \int_{0}^{\infty} d ω {∣ \tilde{G} (z, ω) ∣}^{2} \frac{n_{eff} (ω)}{n_{0}}

where E_p is the total pulse energy, and n₀ is some representative refractive index, here chosen to be 1.45.

The remaining steps first of all involve a separation of the positive and negative frequency components of Eq. (11), assuming that these are well separated, even in the nonlinear polarization term. This requires that the width of the spectrum does not exceed ω ₀/3, where ω0 is some suitably chosen base frequency [16]. This assumption is automatically fulfilled in the PBG fiber examples to be discussed here due to the finite spectral width of the photonic bandgap. Furthermore, a coordinate transformation into a comoving time frame, t′=t-β1z, is done, where β ₁=dβ/dω, evaluated at ω ₀. This is equivalent to transforming β(ω) into β(ω)-β ₁(ω₀)×(ω-ω₀). In addition, ω is shifted by the base frequency ω ₀ and β(ω) is correspondingly measured relative to β₀=β(ω₀). With the normalization conventions in Eq. (14), and the above changes, Eq. (11) becomes:

\frac{\partial \tilde{G} (z, ω)}{\partial z} = i \frac{n_{2} ω^{2} \exp (- i \tilde{β} (ω) z)}{2 π c^{2} β (ω)} \int d ω_{1} d ω_{2} G (z, ω_{1}) G^{*} (z, ω_{1} + ω_{2} - ω) \times

K (ω, ω_{1}, ω_{1} + ω_{2} - ω, ω_{2}) = \int d r_{⊥} F (ω) F (ω_{1}) F (ω_{1} + ω_{2} - ω) F (ω_{2})

\tilde{β} (ω) = β (ω) - β (ω_{0}) - β_{1} (ω_{0}) \times (ω - ω_{0})

n_{2} = \frac{3 χ_{xxxx}^{(3)}}{4 n_{0} ε_{0} c}

Note that this choice of n ₂ differs from the conventional definition of n ₂=3χ⁽³⁾_xxxx/4n²₀ε₀c, which follows from approximating ω/(β(ω)in the prefactor on the RHS of Eq. (16) by 1/n ₀, a convention not adopted in this work. It is from now on implicit that frequency integrations are done over positive frequencies only. Apart from some differences with respect to normalizations these results correspond to the earlier findings in Refs. [12, 14].

If the frequency variation of F(r⊥,ω) is neglected, Eq. (16) becomes:

\frac{\partial \tilde{G} (z, ω)}{\partial z} = i \frac{n_{2} ω^{2} \exp (- i \tilde{β} (ω) z)}{2 π c^{2} β (ω) A_{eff}} \int d ω_{1} d ω_{2} G (z, ω_{1}) G^{*} (z, ω_{1} + ω_{2} - ω) G (z, ω_{2}) R (ω - ω_{1})

A_{eff}^{- 1} = \frac{\int d r_{⊥} F^{4} (r_{⊥})}{{[\int d r_{⊥} F^{2} (r_{⊥})]}^{2}}

where Eq. (15) was used to show that K is in this case equal to the inverse of the well-known effective area parameter. As discussed further in the next section, numerical solution of Eq. (20) is a much easier task than solution of Eq. (16). Many authors have sought to include the effects of mode profile dispersion simply by including frequency dependence into the effective area parameter, to obtain the equation:

\frac{\partial \tilde{G} (z, ω)}{\partial z} = i \frac{n_{2} ω^{2} \exp (- i \tilde{β} (ω) z)}{2 π c^{2} β (ω) A_{eff} (ω)} \int d ω_{1} d ω_{2} G (z, ω_{1}) G^{*} (z, ω_{1} + ω_{2} - ω) G (z, ω_{2}) R (ω - ω_{1})

Since this procedure is obviously not well founded in light of Eq. (16), one may ask if it represents any improvement at all over Eq. (20). In fact, there are several cases in which one can argue for this. One example is a soliton, which shifts its frequency during propagation due to the Raman effect. As the frequency shifts with z, the soliton experiences a changing effective area parameter, even if its spectral width at any particular value of z is small enough to merit a neglect of the frequency variation of F. Another case where the use of Eq. (22) may be adequate is that in which the different spectral components of the propagating light are well separated in time (e.g. in a strongly chirped pulse). Since the Kerr interaction is strictly local in time, and the Raman response function has a width of ~100 fs, the different frequency components in such pulses will mainly experience nonlinear interactions with themselves. The difference between Eq. (16) and (22) becomes apparent when frequency components with a significant difference in the mode profile have a temporal overlap. This observation is quite interesting, since in PBG fibers a strong mode profile dispersion is often associated with a strong group velocity dispersion, meaning that different frequency components of a short pulse will quickly disperse away from each other.

On the other hand, if it is the case that different spectral components do interact significantly with each other, correct inclusion of mode profile dispersion can be important. One well-known example of such interaction is four-wave mixing. Indeed, to derive the four-wave mixing equations with the correct overlap integrals for nonlinear interactions [17], one must start out from Eq. (16) rather than Eq. (22). Another example is the injection of powerful pulses close to a zero-dispersion point of the fiber, a process which is important in supercontinuum generation [1].

3. Approximate forms of the propagation equation

The appearance of K in the double frequency integral in Eq. (16) implies that this integral cannot be evaluated as a series of convolutions, as is the case for the integrals in Eq. (20), and (22). Since convolutions can be done by the fast Fourier transform (FFT) with a numerical effort on the order of Nln(N) where N is the number of points on the frequency mesh, whereas a full 2D integral, evaluated at all values of ω, would require N ³ operations, this is a very significant drawback, severely limiting the practical usefulness of Eq. (16). Therefore, it is important to derive approximate approaches to this equation which are numerically feasible. Kolesik et al. [13] presented a working procedure based on separating radial and angular coordinates in a circular symmetric fiber, but such an approach is not useful in the more general case of microstructured fibers. In this section, two generally applicaple methods are described. One is a numerical procedure, which allows one to solve Eq. (16) with an accuracy that can be systematically improved, at the price of a strongly increased, though not prohibitive, computational load. The second is a modification of Eq. (22), originally proposed in a slightly different form by Mamyshev and Chernikov [12], which does not increase the computational burden at all, and is very easy to implement in GNLSE codes developed to solve Eq. (22).

3.1. Basis set expansion

One general way of making Eq. (16) tractable, is to expand the object K into one-dimensional basis functions:

K (ω, ω_{1}, ω_{1} + ω_{2} - ω, ω_{2}) = \sum_{μ ν γ δ} K_{μ ν γ δ} k_{μ} (ω) k_{v} (ω_{1}) k_{γ} (ω_{1} + ω_{2} - ω) k_{δ} (ω_{2})

The integral on the RHS of Eq. (16) then becomes:

\int d ω_{1} d ω_{2} G (z, ω_{1}) G^{*} (z, ω_{1} + ω_{2} - ω) G (z, ω_{2}) R (ω - ω_{1}) K (ω, ω_{1}, ω_{1} + ω_{2} - ω, ω_{2}) =

G_{μ} (z, ω) = G (z, ω) k_{μ} (ω)

The frequency integrals can now again be done in O[Nln(N)] operations using fast Fourier transforms. The price is, that the number of such integrals to be evaluated now increases with the number of terms in the expansion of K.

One obvious way to bring about an expansion of the form (24) is to expand the field distribution of the guided mode, F(r⊥,ω), into a set of frequency-independent basis functions:

F (r_{⊥}, ω) \approx \sum_{μ} a_{μ} (ω) F_{μ} (r_{⊥})

With such an expansion, the quantities introduced in Eq. (23) become:

k_{μ} (ω) = a_{μ} (ω), K_{μ ν γ δ} = \int d r_{⊥} F_{μ} (r_{⊥}) F_{ν} (r_{⊥}) F_{γ} (r_{⊥}) F_{δ} (r_{⊥})

A simple and systematic procedure for doing this is to use a Taylor-series like expansion:

F (r_{⊥}, ω) \approx F_{0} (r_{⊥}) + \sum_{μ = 1}^{M} F_{μ} (r_{⊥}) {(ω - ω_{c})}^{μ}

The functions F_µ could in principle be determined as the frequency derivatives of F at ω_c, but a more robust procedure is to use a polynomial fitting scheme as detailed below.

To set up the expansion in Eq. (28), a semivectorial finite-difference scheme, implemented in MATLAB, was employed to determine the eigenmodes F(r⊥,ω) of an all-silica PBG fiber. The structure chosen for study consists of a triangular array of Ge-doped rods in a background of pure silica, with a missing rod comprising the low-index core defect. The structure is depicted by the black circles in Fig. 1. Only a quarter of the structure is shown, since the finite-difference calculation was restricted to modes that are even under reflection in the x- and y-axes. Losses are not accounted for in the calculation, so periodic boundary conditions have been assumed for simplicity. Such fibers have been extensively studied in the last few years [18, 19, 20], and their linear propagation properties are well known [21, 22, 23]. It has been shown, that such fibers guide light in the core by the photonic bandgap effect in finite intervals separated by the cutoff frequencies for low-order guided modes in the high-index rods [23]. In the present case, the rods are assumed to have a step-index profile, with a difference in refractive index of 0.02 to the pure silica background. The dispersive properties of such fibers have earlier been found to be well described within the scalar approximation [18]. The rod diameter is 0.5Λ, where Λ is the pitch, or center-to-center distance between neighbouring rods. Propagation in the third bandgap, from V=3.95 to V=5.0 is considered here, where V is the classical V-parameter of the high-index rods:

V = \frac{π d}{λ} \sqrt{n_{rod}^{2} - n_{S i O_{2}}^{2}}

This V-range is well within the bandgap-guiding region, but still a significant mode profile dispersion occurs at the edges, as shown below.

Fig. 1. Guided-mode profiles in the all-silica PBG fiber described in the text. Black circles denote high-index rods. Profiles for V=3.95 (left), V=4.5 (middle) and V=5 (right) are shown. The computational domain shown in the figure is discretized in 250×250 gridpoints.
Download Full Size | PDF

In a first step, guided-mode profiles were calculated at 43 equidistant frequency points in the V-parameter range from 3.95 to 5.0, with Λ=12 µm. The fields in the center and at the edges of this range are shown in Fig. 1. Subsequently, at each point on the finite-difference grid, the 43 field values were fitted to a Mth-order polynomium in the normalized and centered frequency variable (ω-ω_c)/ω_c where ω_c was taken in the center of the frequency range. The set of M+1 polynomial coefficients from each gridpoint then constitute theM+1 basis functions in Eq. (28). Examples of these basis functions are shown in Fig. 2, for M=2.

Fig. 2. Basis functions, as in Eq. (28), derived by fitting the fields calculated by the finitedifference method. The plots show F ₀ (left), F ₁ (middle) and F ₂ (right) for an M=2 expansion.
Download Full Size | PDF

Due to the large number of summation terms in Eq. (24), it is extremely important for the practical use of this approach that a reasonable description of the mode profiles can be obtained with a small number of basis functions. This was checked by calculating the frequencydependent effective area, Eq. (21), using either the calculated mode profiles, F(r⊥,ω), or the approximate mode profiles given by the polynomial expansion in Eq. (28). In Fig. 3, the ‘exact’ effective-area curve is compared to curves calculated with expansions of M=2 (i.e. 3 basis functions in total), and M=3. It is seen, that already the M=2 expansion provides a useful approximation to the effective-area curve, and that the M=3 expansion gives a very satisfactory agreement. The M=2 expansion implies 81 terms in the sum occurring in Eq. (24), whereas M=3 implies 256 terms. The sums can be somewhat reduced by using the symmetry of K under permutation of the frequency arguments. In the implementation used here, a single evaluation of the nonlinear term using the M=3 expansion requires 164 complex and 320 real FFT operations, compared to just two complex and two real FFT’s when using Eq. (22). In the latter case, computational overheads arising from the evaluation of exponential functions and the use of a fairly advanced (Bulirsch-Stoer) stepper algorithm are not completely insignificant, so the total difference in runtime between the two approaches was about a factor of 50. While the use of Eqs. (16), (24) thus lead to a significantly increased computational load, they are no longer unfeasible, and can therefore be used to check the accuracy of simpler formulations. Still further reductions in the computational overhead may be achieved by casting away small terms in K_µνγδ, and/or do linear transformations within the determined basis set to minimize some of these coefficients, but such schemes were not studied in this work.

3.2. Modified GNLSE

The purpose of this subsection is to derive an equation which fully regains the numerical efficiency of the traditional GNLSE, Eq. (22), while providing an improved approximation to the more complete Eq. (16). Let us suppose that some basis-state expansion of the general form (26), but not necessarily of the particular form (28), exists, for which the K ₀₀₀₀ term in Eq. (24) is much more important than the other terms. If only the K ₀₀₀₀ term is retained, this equation becomes:

Fig. 3. Effective area for a PBG fiber as described in the text with Λ=12 µm. The black curve denotes the results of the semivectorial finite-difference calculation, whereas the red and green curves denote results with M=2 and M=3 polynomial expansions for the fields.
Download Full Size | PDF

\int d ω_{1} d ω_{2} G (z, ω_{1}) G^{*} (z, ω_{1} + ω_{2} - ω) G (z, ω_{2}) R (ω - ω_{1}) K (ω, ω_{1}, ω_{1} + ω_{2} - ω, ω_{2}) =

\bar{G} (z, ω) = G (z, ω) a_{0} (ω)

In the same approximation, the frequency-dependent effective area calculated from Eq. (21) becomes:

A_{eff}^{- 1} (ω) = K_{0000} a_{0}^{4} (ω)

With these results, Eq. (16) can be approximated by:

\frac{\partial \tilde{G} (z, ω)}{\partial z} = i \frac{n_{2} ω^{2} \exp (- i \tilde{β} (ω) z)}{2 π c^{2} β (ω) A_{eff}^{1 ⁄ 4} (ω)} \int d ω_{1} d ω_{2} \bar{G} (z, ω_{1}) {\bar{G}}^{*} (z, ω_{1} + ω_{2} - ω) \bar{G} (z, ω_{2}) \overset{}{R} (ω - ω_{1})

\bar{G} (z, ω) = \frac{G (z, ω)}{A_{eff}^{1 ⁄ 4} (ω)}

A similar equation was derived in Ref. [12]. This equation retains the simplicity of the traditional GNLSE and is in fact very easily implemented in codes designed to solve Eq. (22).

The approximation leading to Eq. (33) is correct if the basis states with µ>0 in Eq. (28) have effective areas much larger than F ₀. In the PBG fiber example discussed here, the F ₁ state, shown in the middle panel of Fig. 2, has an effective area ~2.5 times larger than that of F ₀ (left panel), and the area of the F ₂ state (right panel) is about 6 times larger than the F ₀ area. Thus, the validity of Eq. (33) is not obvious in this case. There are, however, two other arguments of a more general nature which favour the use of Eq. (33) over Eq. (22). The first is that both Eq. (16) and Eq. (33) conserve the classical photon number, which with the normalizations used here is proportional to:

N_{phot} \propto \int d ω n_{eff} (ω) \frac{{∣ G (z, ω) ∣}^{2}}{ω}

That this is a conserved quantity when solving Eq. (16) or Eq. (33) is readily verified following the procedure in Ref. [16]. Eq. (22), on the other hand, is found to conserve the quantity

\int d ω n_{eff} (ω) A_{eff} (ω) \frac{{∣ G (z, ω) ∣}^{2}}{ω}

In the presence of strong mode profile dispersion, the difference between these two quantities may be significant.

The second argument favouring Eq. (33) is that it has first-order accuracy in ω-ω _c as will now be demonstrated. If one assumes that the relation

F (r_{⊥}, ω) \approx F (r_{⊥}, ω_{c}) + F' (r_{⊥}, ω_{c}) (ω - ω_{c})

holds, and that all higher-order terms appearing in the evaluation of K can be discarded, the expression for K becomes:

K (ω, ω_{1}, ω_{1} + ω_{2} - ω, ω_{2}) \approx K_{0} + 2 K_{1} (ω_{1} + ω_{2} - 2 ω_{c})

K_{0} = \int d r_{⊥} F^{4} (r_{⊥}, ω_{c}), K_{1} = \int d r_{⊥} F' (r_{⊥}, ω_{c}) F^{3} (r_{⊥}, ω_{c})

In a similar approximation, the inverse effective area becomes:

A_{eff}^{- 1} (ω) \approx K_{0} + 4 K_{1} (ω - ω_{c})

Finally, a first-order expansion of A^1/4_eff in Eq. (33), and use of Eq. (40) yields:

\frac{1}{A_{eff}^{1 ⁄ 4} (ω)} \int d ω_{1} d ω_{2} \bar{G} (z, ω_{1}) {\bar{G}}^{*} (z, ω_{1} + ω_{2} - ω) \bar{G} (z, ω_{2}) R (ω - ω_{1}) \approx \frac{A_{eff}^{3 ⁄ 4} (ω_{c})}{4 A_{eff}^{3 ⁄ 4} (ω_{c})} \times

The last equation is seen to agree with Eqs. (16),(38), thus confirming that Eq. (33) will be accurate if the spectral weight of the pulse extends over a region where the mode profile varies linearly. A similar result cannot be obtained for Eq. (22).

4. Numerical results

This section investigates the agreement between Eq. (22), hereafter denoted GNLSE, Eq. (33), hereafter denoted M-GNLSE, and the full Eq. (16) solved with a polynomial expansion of M=3 as described in the previous section. For a fair comparison, the GNLSE and M-GNLSE are solved using the effective area calculated with the M=3 expansion, e.g. the green curve in Fig. 3. Material dispersion was included in the finite-difference calculations, assuming that the index difference between the high-index rods and the pure-silica background remained constant. The dispersion curve of the fiber with Λ=12 µm is shown in Fig. (4). The zero-dispersion point is at 1.025 µm. The calculations were done on a frequency grid with 2 ¹⁵ points, extending from λ ₀=0.91 µm to λ0=1.15 µm. The value for n ₂, as defined by Eq. (19), was put to 3.857·10^-8 W/µm², and the Raman response function was parametrized according to Ref. [17].

In the first example to be discussed, fundamental solitons are launched at a center wavelength of 1100 nm, with a temporal width T ₀ matched to the peak power P ₀, the effective area and the dispersion coefficient of the fiber according to the well-known relation [17]:

P_{0} = \frac{∣ β_{2} ∣}{γ T_{0}^{2}}

where β ₂ is the second derivative of β with respect to ω, and γ is defined in Eq. (1). The dispersion coefficient of the fiber at 1100 nm is 146 ps/nm/km, and a P ₀ of 100 kW corresponds to a T ₀ of about 32 fs. Due to the Raman effect, the soliton will shift towards longer wavelengths during propagation [24], and it will deviate from the ideal soliton form due to the presence of higher-order group velocity dispersion and mode profile dispersion.

Fig. 4. Dispersion curve for the fiber structure described in the text with Λ=12 µm, and material dispersion included in the calculation.
Download Full Size | PDF

Fig. 5. Average wavelength of a soliton with a peak power of 100 kW launched at 1100 nm as a function of propagation distance. Inset shows the shifts after 10 meters of propagation as a function of peak power, P ₀, calculated with the GNLSE and M-GNLSE.
Download Full Size | PDF

In Fig. 5, the mean wavelength of the soliton is shown as a function of propagation distance for an input P ₀ of 100 kW, as found with the three calculational approaches discussed above. The calculation using the M=3 expansion took about 28 hours on a 3GHz PC, whereas solution of the M-GNLSE required about 35 minutes. The M-GNLSE is seen to provide a near-perfect match to the results of the M=3 expansion, whereas the wavelength shift found with the GNLSE deviates by about 8.4% after 10 meters of propagation. The redshifting of the wavelength happens quickly at the launch end of the fiber, but gradually slows down during propagation. This is understandable, since the soliton adjusts its width to the increasing group velocity dispersion at longer wavelengths, thereby increasing T ₀ according to Eq. (42). The rate of redshifting decreases with increasing T ₀, and is roughly proportional to T ^-4 ₀ for long pulses [24]. The difference between the shifts predicted by the GNLSE and M-GNLSE is roughly constant when measured relative to the total shift at a given propagation distance. In the inset, the average wavelength after 10 meters of propagation is shown as a function of P ₀, calculated with either the GNLSE or M-GNLSE. At the lowest peak power of 10 kW, the difference between the predicted redshifts is about 2.5% of the total shift, indicating that the deviation between the GNLSE and M-GNLSE follows the spectral width of the soliton. It should also be noted that the relative deviations found are of the same order of magnitude as the relative variation in A _{e f f} over the width of the soliton.

The agreement between the redshift predictions of the M-GNLSE and the M=3 expansion is reflected in the calculated spectra after 10 meters of propagation, as shown in Fig. 6 for P ₀=100 kW. Apart from being shifted in frequency, the GNLSE spectrum is wider than those predicted by the other methods, because it resides at wavelengths with a lower group velocity dispersion.

Fig. 6. Spectral weight of a soliton launched at 1100 nm with a peak power of 100 kW after 10 meters of propagation. |G(z,ω)|² has been multiplied by c/λ ² to convert into units of energy per wavelength.
Download Full Size | PDF

As a second numerical example, the propagation of a high power pulse injected at the zero-dispersion point (1025 nm) is considered. The injected pulse has a Gaussian shape, with a full width at half maximum (FWHM) of 0.93 ps. The peak power is 50 kW, and the pulse is propagated over a distance of 60 cm. In the initial stages of propagation, self phase modulation shifts spectral weight into both normal- and anomalous-dispersion spectral regions. Pumping at a zero-dispersion point ensures that the generated sidebands can copropagate over a substantial length of fiber, making the problem of interest for the present investigation of mode profile dispersion. After 60 cm of propagation, a large spectral broadening has already occurred, as seen in Fig. 7. The GNLSE is seen to provide qualitatively correct predictions, however there are significant quantitative deviations from the full theory in the various spectral features. The results of the M-GNLSE, on the other hand, is seen to provide a perfect match to the predictions of the full M=3 expansion.

Fig. 7. Spectrum of a Gaussian pulse of peak power 50 kW, initial width of 0.93 ps and a center wavelength of 1025 nm after 60 cm of propagation.
Download Full Size | PDF

The above results indicate that the M-GNLSE provides results in perfect agreement with those of the full M=3 expansion, for a strongly reduced computational effort. It is, however, important to keep in mind that the M-GNLSE is not exact. For example, it does not correctly represent the overlap integrals controlling the strength of four-wave mixing processes between frequencies with substantially differing mode profiles [17]. In the example fiber studied here, pump pulses injected in the normal-dispersion regime slightly below the zero-dispersion point are phased-matched to signal/idler waves at well-separated wavelengths. Studies of such four-wave mixing processes showed that neither the GNLSE nor the M-GNLSE consistently reproduced the results of a full calculation. A related problem is the modulation instability of a ps pulse launched at the long-wavelength edge of the PBG transmission window, where the mode profile dispersion is strong. To show an example of this process, a Gaussian pulse of 11 ps FWHM, a peak power of 40 kW and an initial center wavelength of 1120 nm was propagated over a distance of 1.2 meters. Modulational instability creates sidebands to the main peak, which then evolve into a complex spectral structure, with significant parts of the spectral weight being shifted both up and down in frequency. Due to the fairly long length of the pulse, dispersive effects do not prevent the different frequency components from interacting with each other. The resulting spectra are shown in Fig. 8. It is clear that neither the GNLSE nor the M-GNLSE provide a perfect reproduction of the full calculation in this case. Considering the integrated spectral functions shown in the inset of Fig. 8, the M-GNLSE provides a slightly improved description of the spectral weight distribution at short and intermediate wavelengths, whereas the GNLSE is better at the long-wavelength part of the spectrum. Calculations with other parameters did not indicate this as a consistent trend. One advantage of the M-GNLSE in this respect is that it conserves the true photon number, Eq. (35), and therefore reproduces the pulse energy found in the full calculation with good accuracy. The spectrum calculated with the full M=3 expansion integrates to 471.4 nJ, the M-GNLSE spectrum to 471.3 nJ, and the GNLSE spectrum to 465.3 nJ. These numbers are not exactly equal to the actual pulse energies due to the neglect of the factor n_eff(ω)/n₀ in Eq. (15).

Fig. 8. Spectrum of a Gaussian pulse of peak power 40 kW, initial width of 11 ps and a center wavelength of 1120 nm after 1.2 m of propagation. Inset show the cumulative integral of the spectral functions, starting from the long wavelength edge.
Download Full Size | PDF

5. Conclusion

In conclusion, the correct treatment of mode profile dispersion in the generalized nonlinear Schrödinger equation has been investigated, using a solid-core all-silica PBG fiber as an example case. In agreement with previous works, it has been shown that mode profile dispersion in the general case calls for a generalization of the effective-area concept. A systematic numerical approach to the full equation has been developed, and used to test the performance of the conventional GNLSE as well as a different approximate equation, with the same degree of numerical complexity. The latter is found to have superior performance compared to the full calculation in almost all cases studied, and should therefore be generally preferable over the traditional GNLSE formulation.

Acknowledgements

It is a pleasure to acknowledge Peter John Roberts, Per Dalgaard Rasmussen, Michael Frosz and Ole Bang for their many helpful suggestions and critical readings of the manuscript. This work was financially supported by the Danish High Technology Foundation.

References and links

1. J. M. Dudley, G. Genty, and S. Coen, “Supercontinuum generation in photonic crystal fiber,” Rev. Mod. Phys.78(4) (2006).

2. A. Fuerbach, P. Steinvurzel, J. Bolger, A. Nulsen, and B. Eggleton, “Nonlinear propagation effects in antiresonant high-index inclusion photonic crystal fibers,” Opt. Lett. 30(8), 830–832 (2005). [CrossRef]

3. A. Fuerbach, P. Steinvurzel, J. Bolger, and B. Eggleton, “Nonlinear pulse propagation at zero dispersion wavelength in anti-resonant photonic crystal fibers,” Opt. Express 13(8), 2977–2987 (2005). [CrossRef]

4. D. Ouzounov, F. Ahmad, D. Muller, N. Venkataraman, M. Gallagher, M. Thomas, J. Silcox, K. Koch, and A. Gaeta, “Generation of megawatt optical solitons in hollow-core photonic band-gap fibers,” Science 301(5640), 1702–1704 (2003). [CrossRef]

5. D. G. Ouzounov, C. J. Hensley, A. L. Gaeta, N. Venkateraman, M. T. Gallagher, and K. W. Koch, “Soliton pulse compression in photonic band-gap fibers,” Opt. Express 13(16), 6153–6159 (2005). [CrossRef]

6. C. J. Hensley, D. G. Ouzounov, A. L. Gaeta, N. Venkataraman, M. T. Gallagher, and K. W. Koch, “Silica-glass contribution to the effective nonlinearity of hollow-core photonic band-gap fibers,” Opt. Express 15(6), 3507–3512 (2007). [CrossRef]

7. F. Gerome, K. Cook, A. George, W. Wadsworth, and J. Knight, “Delivery of sub-100fs pulses through 8m of hollow-core fiber using soliton compression,” Opt. Express 15(12), 7126–7131 (2007). [CrossRef]

8. J. Lægsgaard, N. A. Mortensen, J. Riishede, and A. Bjarklev, “Material effects in airguiding photonic bandgap fibers,” J. Opt. Soc. Am. B 20, 2046–51 (2003). [CrossRef]

9. N. Karasawa, S. Nakamura, and N. Nakagawa, “Comparison Between Theory and Experiment of Nonlinear Propagation for A-Few-Cycle and....” IEEE J. Quantum Electron. 37(3), 398–404 (2001). [CrossRef]

10. G. Chang, T. B. Norris, and H. G. Winful, “Optimization of supercontinuum generation in photonic crystal fibers for pulse compression,” Opt. Lett. 28(7), 546–548 (2003). [CrossRef]

11. B. Kibler, J. M. Dudley, and S. Coen, “Supercontinuum generation and nonlinear pulse propagation in photonic crystal fiber: influence of the frequency-dependent effective mode area,” Appl. Phys. B 81(2–3), 337–342 (2005). [CrossRef]

12. P. Mamyshev and S. Chernikov, “Ultrashort-pulse propagation in optical fibers,” Opt. Lett. 15(19), 1076–1078 (1990). [CrossRef]

13. M. Kolesik, E. Wright, and J. Moloney, “Simulation of femtosecond pulse propagation in sub-micron diameter tapered fibers.” Appl. Phys. B: Lasers & Optics 79(3), 293–300 (2004). [CrossRef]

14. A. Ferrando, M. Zacares, P. de Cordoba, D. Binosi, and A. Montero, “Forward-backward equations for nonlinear propagation in axially invariant optical systems,” Phys. Rev. E (Statistical, Nonlinear, and Soft Matter Physics) 71(1), 16,601 (2005).

15. Y. Mizuta, M. Nagasawa, M. Ohtani, and M. Yamashita, “Nonlinear propagation analysis of few-optical-cycle pulses for subfemtosecond compression and carrier envelope phase effect,” Phys. Rev. A (Atomic, Molecular, and Optical Physics) 72(6), 63,802 (2005).

16. K. Blow and D. Wood, “Theoretical description of transient stimulated Raman scattering in optical fibers,” IEEE J. Quantum Electron. 25(12), 2665–2673 (1989). [CrossRef]

17. G. P. Agrawal, Nonlinear Fiber Optics (Academic Press, San Diego, 2001).

18. J. Riishede, J. Lægsgaard, J. Broeng, and A. Bjarklev, “All-silica photonic bandgap fibre with zero dispersion and large mode area at 730 nm,” J. Opt. A: Pure and Applied Optics 6, 667–70 (2004). [CrossRef]

19. A. Argyros, T. Birks, S. G. Leon-Saval, C. M. B. Cordeiro, F. Luan, and P. S. J. Russell, “Photonic bandgap with an index step of one percent,” Opt. Express 13, 309–14 (2005). [CrossRef] [PubMed]

20. G. Bouwmans, L. Bigot, Y. Quiquempois, F. Lopez, L. Provino, and M. Douay, “Fabrication and characterization of an all-solid 2D photonic bandgap fiber with a low-loss region (<20 dB/km) around 1550 nm,” Opt. Express 13, 8452–9 (2005). [CrossRef] [PubMed]

21. A. K. Abeeluck, N. M. Litchinitser, C. Headley, and B. J. Eggleton, “Analysis of spectral characteristics of photonic bandgap waveguides,” Opt. Express 10, 1320–33 (2002). [PubMed]

22. N. M. Litchinitser, S. C. Dunn, B. Usner, B. J. Eggleton, T. P. White, R. C. McPhedran, and C. M. de Sterke, “Resonances in microstructured optical waveguides,” Opt. Express 11, 1243–51 (2003). [CrossRef] [PubMed]

23. J. Lægsgaard, “Gap formation and guided modes in photonic bandgap fibres with high-index rods,” J. Opt. A: Pure and Applied Optics 6, 798–804 (2004). [CrossRef]

24. J. P. Gordon, “Theory of the soliton self-frequency shift,” Opt. Lett. 11(10), 662–664 (1986). [CrossRef]

Mode profile dispersion in the generalized nonlinear Schrödinger equation

Abstract

1. Introduction

2. Derivation of a 1+1D propagation equation

3. Approximate forms of the propagation equation

3.1. Basis set expansion

3.2. Modified GNLSE

4. Numerical results

5. Conclusion

Acknowledgements

References and links

Cited By

Figures (8)

Equations (49)

Optics Express