Tunable Modulational Instability Sidebands via Parametric Resonance in Periodically Tapered Optical Fibers

We analyze the modulation instability induced by periodic variations of group velocity dispersion and nonlinearity in optical fibers, which may be interpreted as an analogue of the well-known parametric resonance in mechanics. We derive accurate analytical estimates of resonant detuning, maximum gain and instability margins, significantly improving on previous literature on the subject. We also design a periodically tapered photonic crystal fiber, in order to achieve narrow instability sidebands at a detuning of 35 THz, above the Raman maximum gain peak of fused silica. The wide tunability of the resonant peaks by variations of the tapering period and depth will allow to implement sources of correlated photon pairs which are far-detuned from the input pump wavelength, with important applications in quantum optics.


Introduction and motivations
Parametric resonance (PR) is a well-known instability phenomenon which occurs in systems the parameters of which are varied periodically during evolution, see the classical treatments in Refs. [1,2]. For example, a harmonic oscillator the frequency of which is forced to vary in time, will become unstable if its internal parameters and the amplitude of the frequency variation happen to be inside special regions, known as resonance tongues. The study of the properties of resonance tongues has a long history and relies on a variety of geometrical approaches [3,4].
It is natural that such a general phenomenon was associated to the equally important instability process that is ubiquitous in infinite dimensional dynamical systems: modulation instability (MI), also known as Benjamin-Feir instability [5,6]. MI is known to exist in different branches of physics such as fluid-dynamics [7], plasma physics [8,9,10], Bose-Einstein condensates [11] and solid-state physics [12]. In nonlinear optics [13], it manifests itself as an exponential growth of two spectrally symmetric sidebands on top of a plane wave initial condition, by virtue of the interplay between the cubic Kerr nonlinearity and the group velocity dispersion (GVD). In optical fibers it leads to the breakup of a plane wave into a train of normal modes of the system, i.e. solitons [14,15].
The link between PR and MI has been established in 1993, see [16], in relation to the periodic re-amplification of signals in long-haul telecommunication optical fiber cables. This was based on a nonlinear Schrödinger equation (NLS) where the coefficient of the nonlinear term is varied along the propagation direction. Importantly, this peculiar type of MI occurs in both normal and anomalous GVD. This prediction was later partially verified in experiments, see [17,18].
Moreover, in long-haul fibers, dispersion management is a commonly used technique which introduces periodic modulation of fiber characteristics. The possibility of instability phenomena disrupting adjacent communication channels has been thoroughly analyzed, see e.g. [19,20,21,22]. Specifically, in [19] the partial suppression of the conventional MI in anomalous GVD due to a large swing dispersion management is discussed, while in [20] the degenerate case of zero average dispersion was studied. The combination of both loss and dispersion compensation is studied in [21,22]. The main interest in those works was on step-like variations of the GVD coefficient.
At the same time the effects of smooth periodic or random variations of fiber parameters were studied in [23,24,25]. Also some work has been done on the effect of the perturbation of fiber parameters on soliton propagation [26,27].
It turns out that the variation of dispersion and nonlinearity can enhance or suppress the PR, while higher order nonlinear effects such as self-steepening proved less important. Quite surprisingly, experiments on micro-structured fibers have been reported only recently for the first time, see Refs. [28,29], where a photonic-crystal fiber (PCF, [30]) of varying diameter is used. In that experiment, the dispersion is periodically switched from normal to anomalous, but this feature is not required to achieve PR, while the effect of Raman scattering plays an important role in the relative magnitude of the PR peaks.
The conventional explanation is in term of a grating-assisted phase matching process [16,21,22,28], with the exception of Ref. [24], which relies on the theory of Mathieu equation -a paradigmatic example of the application of Floquet theory, see Refs. [2,3]-and analyzes the case of nonlinearity variations only.
In this work we present an improved analytical treatment based on the so-called Poincaré-Lindstedt perturbation method and averaging method, which provide excellent predictions on resonant frequency, gain and bandwidth of the unstable peaks appearing in the output spectra. The new feature here is the tunability of the instability bands, specifically the lowest order one, which possesses a gain larger than the other PR peaks. Finally we provide physical estimates for a PCF-based system which leads to instability bands with frequency detunings larger than 13 THz, the detuning of the Raman gain peak of fused silica. This system can be used for the generation of entangled photon pairs in a narrow bandwidth, far detuned from the Raman gain peak, which would reduce the impact of Raman-induced decorrelations, with important applications in quantum optics, especially in quantum computation and cryptography [31, 32].

NLS with varying parameters and derivation of Hill equation
Let us consider the dimensionless NLS equation for a slowly-varying electric field envelope A(z,t) of a guided mode at carrier frequency ω 0 , with both GVD and nonlinearity varying along the propagation direction, namely s and n are normalized coefficients, s(z) ≡ β 2 (z)/β 0 2 and n(z) ≡ γ(z)/γ 0 , where β 2 (z) and γ(z) are the GVD and nonlinear coefficients, respectively, and the 0 superscript denotes their mean values. n and s are assumed to be periodic functions of z. Finally z ≡ Z/Z nl is the dimensionless distance in units of the nonlinear length Z nl ≡ (γ 0 P t ) −1 , and t ≡ (T − v 0 g Z)/T s is the dimensionless retarded time in units of T s ≡ Z nl |β 0 2 |, and v 0 g is the mean group velocity. P t is the total input power injected in the mode, and A is the dimensionless modal intensity scaled by √ P t . We look for for a steady state solution of (1) in the form A = √ P 0 exp (iφ (z)): it can be verified that φ (z) = P 0 z −∞ n(z ′ )dz ′ . We then perturb this steady state by adding a small complex time dependent contribution a(z,t), i.e. A(z,t) = √ P 0 + εa(z,t) exp (iφ (z)), with ε ≪ 1. Inserting this ansatz in Eq. (1) and taking only the terms which are first order in ε, one finds that a obeys the following equation: Finally we substitute in (2) the ansatz a(z,t) = a A (z)e −iδt + a S (z)e iδt to obtain the following Schrödinger equation where the dot denotes z-derivative, H s (z) is a non-Hermitian Hamiltonian (which thus allow unstable modes to grow), |ψ ≡ (a A , a * S ) T , andσ i are the Pauli matrices. In the remaining part of this section we will assume that dispersion and nonlinearity exhibit the simplest possible periodic behavior s(z) = s 0 +s(z) = s 0 + hs 1 cos αz, n(z) = n 0 +ñ(z) = n 0 + hn 1 cos αz, where generally s 0 = ±1 for normal (anomalous) dispersion and n 0 = 1; α is the normalized spatial angular frequency for the parameter oscillation. The forcing amplitude is controlled by parameter h, which must be small to guarantee the validity of our perturbative expansions. However, we find below that our estimates are reliable even for h ∼ 0.5. The conventional way in which PR is approached in fiber optics is to split the dispersion coefficient into two parts: a constant average term and an oscillating term. By means of a change of variables this latter is included in a single varying nonlinear coefficient, see Refs. [19,21]. Then all quantities are expanded in Fourier series and it is assumed that only one Fourier coefficient resonates at each PR order. This leads to the nonlinear phase-matching condition for the m-th peak resonant frequency s 0 δ 2 m + 2n 0 P 0 = mα (5) where m is the PR order and δ m the corresponding resonance detuning. This condition is valid for m > 0 (m < 0) if GVD is normal (anomalous). Below, in our derivation, we assume instead that m is a positive integer. Eq. (5) can also be understood as a grating assisted phase matching condition, the grating being the periodic modulation of the fiber parameters. We performed detailed numerical simulations and verified that the result of (5) is quite coarse. The purpose of the present section is to prove that a more accurate estimate can be made, in a more general case than the variation of nonlinear coefficient only, which was already discussed in [24]. We define the parameter space as the plane (δ , h) and we investigate the different features of PR instability (resonant frequency, bandwidth and gain) on this plane, for a fixed value of |n 1 |/|s 1 |.
In order to compare to the theory of PR in a classical mechanical oscillators, it is instructive to derive a 2nd order ODE from the system of Eqs. (3). Let us first transform it in phase-quadrature variables, by applying the rotation We define the new state |φ = R|ψ which satisfies a system analogous to that in Eq. (3), with Hamiltonian matrix with c 1 (z) = − 1 2 s(z)δ 2 and c 2 (z) = c 1 (z) − 2n(z)P 0 . We can thus more easily derive a second order ODE for both φ 1,2 :φ 1,2 −ċ 1,2 which still contains a first derivative of φ 1,2 . This can be eliminated by the substitution We finally obtain the following ODË Notice that in contrast to a usual dissipative term, which corresponds to a threshold in the magnitude of the perturbation needed to achieve PR instability (see Ref. [2]), our variable transformation is simply an oscillating factor, which does not affect the instability. The left-hand side of (8) contains, even in our simple case of Eq. (4), several harmonic terms. This means that Eq. (8) is a Hill equation, which generalizes the Mathieu equation found in [24], whereċ 1 = 0. It is well known that Mathieu equation is an exceptionally simple case, [1,4], with a predictable structure of stability and instability regions in the parameter space. Moreover since we consider a Hill equation, we expect the instability regions to be irregular and appear in the form of instability islands separated by stable parts. This is consistent with the statements made in Ref.
[19], i.e. that a large switching of dispersion may suppress instability.
Despite Eq. (8) highlights all these important properties of the system, it is difficult to handle analytically. We thus turn back to (6) and derive our analytical estimates from it, in order to report the simplest possible treatment.
In the following part we present the main results of our calculations, which rely on the commonly used methods of averaging and on the Poincaré-Lindstedt perturbation method, see [33].

Estimate of Resonant Frequency and Parametric Gain
We now explore the properties of the system with H pq defined in (6). In the limit of vanishing perturbation we have a simple harmonic oscillator written in complex variables φ 1,2 . It is widely known that, in this limit, parametric resonance occurs if the natural frequency of the oscillator is a multiple of half the forcing frequency and this is the basic point behind any perturbative approach, see [2]. Thus, by setting , the natural frequency of the unforced oscillator is simply ω 0 = c 0 1 c 0 2 . The resonance condition becomes ω 0 = mα/2, with m = 1, 2, 3, . . . which expressed in terms of the NLS parameters corresponds to a detuning Note that, if we assume mα 2n 0 P 0 ≫ 1, we obtain the quasi-phase matching condition given in Eq. (5) (taking care of the different convention on m, which in our relation, Eq. 10 is only positive, while in Eq. (5) can also be negative), which is thus a coarse approximation if α ≈ n 0 P 0 .
In order to obtain the peak gain at resonance we use the method of averaging, which consists in posing substituting in (9) and averaging over a period of the forcing term T z ≡ 2π/α, we obtain where ρ 1 ≡ hs 1 /s 0 and ρ 2 ≡ h(δ 2 s 1 + 4n 1 P 0 )/(δ 2 s 0 + 4n 0 P 0 ). It is thus apparent that the m = 1 resonance occurs at 2ω 0 = α and in this case the matrix of the system (12) has purely real eigenvalues of opposite sign. The positive one corresponds to the peak gain of the first PR band and turns out to be the maximum achievable value. This reads as where δ 1 is calculated according to Eq. (10). We notice that there may exist a set of parameter values which suppresses or even forbids the occurrence of the instability, specifically s 1 /s 0 = n 1 /n 0 . This implies that, e.g., in normal GVD (s 0 = 1) if both nonlinearity and dispersion undergo parallel increase and decrease, the instability is suppressed. Instead the same condition maximizes the gain under anomalous GVD (s 0 = −1).
The higher order resonances are more difficult to characterize by this method, but their gain is generally of order h 2 , since in this case the first-order contribution vanishes, see Eq. (12).

Estimate of resonance bandwidth
Floquet theory predicts that our system (9) has quasi-periodic solutions (composed by a periodic function with period T z and a phase-factor, a complex function of unit absolute value). Moreover, stability margins correspond to a pair of periodic or anti-periodic solutions, defined by |φ (z + T z ) = ±|φ (z) , which possess periods equal to T z and 2T z respectively. Bearing this in mind, we apply the Poincaré-Lindstedt method, by making the following perturbation expansion in powers of h: (14) We choose d 0 such that it corresponds to the first resonant frequency, i.e. d 0 = −δ 2 1 /2. At zeroth order we obtainφ 10 The first resonant band occurs at ω 0 = α/2, thus Eq. (15) has solutions with period 2T z . For a fixed value of h, the stability margins correspond in general to different frequencies, which implies they must correspond to linearly independent eigenfunctions. The higher order corrections, which provide the stability margins, are obtained by solving for the successive terms of the perturbation series (|φ 1 (z) , |φ 2 (z) , . . .) by making use alternatively each of the independent solutions of Eq. (15). We thus first pose φ 10 = cos ω 0 z, and at first order in h we obtain: We then impose that the secular term vanishes and solve for d 1 , If we set φ 10 = sin ω 0 z we obtain the same value with opposite sign, so that the instability margins of the first band satisfy δ 2 The four solutions of Eq. (17) are denoted by ±δ ± 1 , and the bandwidth (i.e. the range of unstable detuning values) is given by Higher order corrections and bandwidth of overtone resonances can be found by solving for the higher order terms |φ 1 (z) , |φ 2 (z) , . . . of the perturbative expansion, but the calculation is too lengthy to report here.

Comparison of analytical and numerical results of the Hill equation
The Hill equation (8) can be solved numerically by means of an ODE solver, following the prescriptions of the Floquet theory. The monodromy matrix is calculated by setting two linearly independent initial conditions and evaluating the solution at z = T z , see [3]. The eigenvalues of the monodromy matrix provide the instability gain (Floquet exponents).
Since PR bands are typically narrow, this requires a fine grid in the (δ , h) parameter space. In order to speed-up the calculation we (i) compute the exact instability margins, then (ii) calculate the gain values of each instability tongue.
The calculation of the exact stability margins can be carried out by the standard Hill determinant method or harmonic balance based on the truncated Fourier expansion of the variables (φ 1,2 ) and forcing terms, see [34]. In contrast with more conventional spectral problems, where the eigenvalue appears explicitly in the equation, we need to find the values of detuning δ at a fixed h, while the coefficients depend on δ in a nontrivial way. This implicit dependence is solved by means of a root finding algorithm based on the minimization of the least singular value, see [35]. The second step involve the numerical evaluation of the Floquet exponents in the close proximity of each band. We compare the results of the numerical calculation of resonant tongues and analytical estimates for the case of maximal gain and bandwidth, i.e. n 1 /n 0 = −s 1 /s 0 , see (13) and (16).
First we show in Figs. 1 and 2 the instability regions for both normal and anomalous GVD at a fixed frequency of parameter variation α = 10 for the first two PR bands, m = 1 (a) and m = 2 (b). Each resonance tongue stems exactly from the detuning predicted by Eq. (10) and the maximum gain (shown in the insets) only slightly drifts away from that value. Our analytical estimates refer to the first resonance band and are shown in Figs. 1(a) and 2(a) to be quite accurate. Instead in Figs. 1(b) and 2(b) only numerical results are provided. We observe that the second order PR exhibits a threshold value for h below which the instability gain is virtually zero. This occurs also for higher-order PR and both in normal and anomalous GVD. As explained above this is not due to the first derivative in Eq. (7) which is not a damping term, but can be qualitatively ascribed to the the fact that in the Hill equation the forcing function contains several overtones of α and they could in special cases suppress completely a subset of resonances; see [33] and references therein.  Fig. 2. Same as in Fig. 1, but in the anomalous GVD regime (s 0 = −1). In this case n 1 = +s 1 = 1. We provide also the curves of resonant frequency, gain and bandwidth as a function of α for fixed amplitude of the variation of nonlinearity and dispersion, |s 1 | = n 1 = 1 and h = 0.5, Figs. 3 and 4. The forcing terms s 1 and n 1 are of equal amplitude as above, and the perturbation is quite large. Nevertheless our estimates for the first band prove quite reliable; moreover the resonant frequency hardly differs from its h → 0 value at any resonant peak.
We notice that in the normal GVD regime (Fig. 3) each resonant frequency converges to zero as α → 0 while the gain g → h m from below as α → ∞. Anomalous dispersion (Fig. 4) causes gain to approach the same limit but from above, while the resonant frequency is limited by the conventional MI band which cuts off at δ = 2.

Including higher order effects
To conclude this section we briefly discuss how to compute the resonant frequency in a generalized NLS model including higher order terms. As above the linearized equation can be cast as where we split as before the Hamiltonian in two parts, a constant and an oscillating one. We also assume that any common diagonal term of H gnls is eliminated by a suitable change of variable, which is always possible. At that point we have a traceless 2 × 2 matrix, thus if we define ±ω 0 to be the eigenvalues of H 0 , and assume they are both real (for a certain choice of parameters and detuning), the resonance condition becomes ω 0 = c 0 1 c 0 2 = α/2 and from the expression of ω 0 we derive the corresponding detuning.
We consider a generalized model which includes higher-order dispersion (HOD) terms up to the fourth order. It is well known that only even HOD terms contribute to MI [36]. We thus just need to modify c 1,2 in Eq. (6) by making the following substitution where β n 4 is the normalized fourth order dispersion (FOD) and is assumed to be constant. The third order dispersion (TOD) appears only as a common diagonal term and can be removed. The resonance condition is thus a quartic equation in δ 2 and if β n 4 < 0 there may exist more than just one pair of sidebands which undergo parametric resonance. This is analogous to the conventional MI in the presence of HOD, see e.g. [37].

Comparison with numerical solution of NLS
In order to assess the correctness of our analysis we solve Eq. (1) by means of the split-step Fourier method. We use the same parameters as above: s 0 = 1, α = 10, −s 1 = n 1 = 1 and h = 0.5. We use a large α in order to achieve widely separated PR resonances and large gain in the normal GVD regime, see Fig. 3. Finally we operate under normal GVD, because anomalous GVD gives rise to the conventional MI bands which have a larger gain (four times larger than the PR gain in the case of h = 0.5).
In Fig. 5(a) we show the output spectrum after a propagation distance z = 38. One can clearly distinguish the first three peaks (m = 1, 2, 3), plus two peaks corresponding to the four-wave mixing of the carrier and the m = 1 PR band. In table 1 we report the values of peak detuning and bandwidth extracted from the spectrum of Fig. 5(a), the numerical results of the linearized problem discussed in the previous section and predictions obtained by phase-matching considerations, Eq. (5). The peak positions are in very good agreement with the results of the previous section. The approximate phase-matching always underestimates the value of the actual detuning.
In Fig. 5(b) we show how the three PR peaks evolve during propagation: the instability leads to amplification at the rate predicted by the linearized problem, i.e. Eq. (13) (see the dashed lines), apart from saturation occurring towards the end of the evolution.
Next we discuss the small scale oscillations reported in [28, 29], i.e. an amplificationdeamplification cycle undergone by the resonant peaks. Since Fig. 5(b) does not allow us to visualize them, it is useful to introduce a new set of parameters which gives larger gain, so that fewer periods are sufficient. Thus we choose h = 0.9 and normal GVD conditions. Moreover, in order to accurately visualize the oscillations, we set a larger period, i.e. α = 5. In Fig. 5(c) we show the evolution of the first peak, together with the solution of the averaged problem, Eq. (11), with ω 0 = α/2 (a phase shift is applied so that the two curves overlap). The numerical and approximate solutions of the linear problem agree quite well, thus providing a good explanation of the process without explicitely resorting to a three wave dynamical system (as in [29]). The amplification-deamplification depends in practice exclusively on α, i.e. it stems from the forced oscillations impressed on resonance by the variation of parameters, superimposed to the unstable growth. This is completely analogous to a pendulum the pivot of which is displaced periodically: a small displacement from the position in which the bob points down initiates oscillations (at the natural frequency) which are further amplified at each cycle.
For the sake of completeness we compare, in table 2, the position of resonant peaks with the predictions of the grating phase matching and the Hill equation. Our estimate Eq. (10) is the closest to numerical NLS result despite it is of zeroth order in h, while the simple phasematching argument, Eq. (5)

The design of a periodically tapered PCF
In this last section we propose a feasible system which supports PR instability peaks exhibiting relatively large frequency detunings. This is important, e.g., in quantum optics where the implementation of new sources of entangled photons with a reduced Raman decorrelation is of great practical interest, and this can provide an alternative approach to what proposed in Refs. [31,32]. In our calculations we use a periodically tapered PCF whose index contrast, small core and design flexibility allow to obtain a large nonlinearity, together with regions of small GVD. A similar system was already used in [28,29], but here we predict the possibility of observing far detuned, tunable instability peaks, by employing short tapering periods.
At first we explore the design space (pitch and filling fraction) in order to operate in a region of small normal GVD near λ 0 = 1064 nm, and such that the zero-dispersion point (ZDP) can be approached by slightly varying the PCF geometry. The modal analysis is performed by means of COMSOL Multiphysics [38].
In Fig. 6 we show the properties of a PCF made of pure silica with a triangular lattice of air holes. We assume the air filling fraction d/Λ = 0.4 (d is the air hole diameter, Λ is the pitch) to be constant so that by varying Λ, d is adjusted accordingly. We are interested mainly in two quantities: the GVD (β 2 ), see We propose to operate between Λ min = 3.25 µm and Λ max = 3.6 µm. In this range, β 2 ∈ [0.4, 2.7] ps 2 /km and γ ∈ [7.8, 9.4] W −1 /km. Thus the GVD undergoes, when compared to its average value, very large oscillations, equivalent to h ≈ 0.85 in Eq. (4). Note that we do not cross into the negative GVD region in order to inhibit any spurious occurrence of the classical MI. Moreover the dependence of GVD on Λ is, in good approximation, linear. Thus a cosineshaped tapering of the PCF leads to a cosine variation of GVD. The value of γ is only slightly modified by the tapering and can be also approximated as a cosine. This allows to straightforwardly apply the theory developed in the previous sections. We thus use the following simple formulas for the variations of the parameters: Λ = Λ 0 + Λ 1 cosαz, β 2 = β 0 2 + β 1 2 cosαz, γ = γ 0 + γ 1 cosαz, where all the parameters are reported in table 3. We use as above the superscripts 0 and 1 to denote the average value and first Fourier coefficient, and we useα for the dimensional spatial frequency of tapering (α = α/Z NL ). As far as HOD terms are concerned, they are in a good approximation constant and are also reported in table 3. Finally in our simulations we include the self-steepening term and stimulated Raman scattering response of silica, see [36], in order to obtain a very realistic simulation.  We set our design target: the period of variation along the direction of propagation (T Z ≡ 2π/α) has to be such that the m = 1 PR band is located at ∆ f = ±35 THz. This quite large detuning allows to clearly distinguish PR from Raman gain, in contrast with [28] where PR peaks occur in the Raman gain spectral region and thus are further amplified by it (in the cited paper the spectrum exhibits a frequency asymmetry which is a clue to the enhancement of sidebands provided by Raman amplification).
We verified that self-steepening and Raman effects are of little importance in our case, nevertheless at large detuning we cannot neglect the effect of HOD. As discussed above, in sec. 2.5, only even order terms contribute to the instability, and amount to a simple modification of c 1,2 in Eq. (6). In order to complete our design, we thus set P t = 25 W, correct the PR condition including β 4 and finally obtainα = 47.7 m −1 , corresponding to a period of tapering of T Z = 13.2 cm. Neglecting β 4 , the estimated period would be less than 10 cm. The predicted gain at ∆ f = ±35 THz is approximately g ≈ 100 km −1 . We simulate a 250 m long fiber with a periodically tapered central part of about 225 m, which corresponds to 1708 periods. The power level and fiber tapering periods are only slightly more demanding than those used in [28]. The use of a highly nonlinear fiber (for example by using materials other than fused silica) can scale down power levels and thus nonlinear lengths conveniently, but this is beyond the scope of the present work.
The spectrum at the end of the propagation is shown in figure 7. We notice two main peak pairs appearing beyond Raman-gain bands (the two broad asymmetric peaks at ∆ f = ±13 THz). The first pair is located at ∆ f = ±34.8 THz, the second pair at ∆ f = ±56 THz. The bandwidth of each individual peak belonging to the first pair is 0.25 THz, and the gain agrees well with our theoretical prediction (g ≈ 100 km −1 ). The second PR peak pair exhibits smaller gain (g ′ ≈ 70  Fig. 7. Output spectrum after the propagation in a periodically tapered PCF. All fiber parameters are listed in Tab. 3. We detect the two Raman gain bands at ∆ f = ±13 THz, Stokes (red-detuned) and anti-Stokes (blue-detuned), the former exhibits as expected a larger gain. We also label the main two PR instability peak pairs which correspond respectively to the design requirement of ∆ f = 35 THz and to the additional phase matching allowed for by FOD. km −1 ) and narrower bandwidth, and can be ascribed to FOD. which allows for an additional solution of the nonlinear phase matching condition, see e.g. [37]. In summary, the presence of β 4 < 0 has the advantage of leading to a period T Z larger than that predicted neglecting all HOD terms. It has also a drawback that an additional PR band appears, but we verified this applies, for our specific choice of parameters, only at the first PR order. If we include losses in our simulations we obtain a smaller available power, which affects mainly the peak gain, see Eq. (13). Indeed the resonant detuning, see Eq. (10), is nearly independent on input power for α ≈ 200, which corresponds to ourα.

Conclusion
In this work we have presented a thorough analysis of parametric resonance instabilities occurring in a generalized nonlinear Schrödinger equations with varying dispersion and nonlinearity. We have shown that, in the case of GVD and nonlinear coefficients varying in a simple sinusoidal way, it is possible to predict the maximum gain and instability margins. The calculation of the resonant frequencies is possible even in more general cases, since it depends only on the average values of coefficients and the period of variation. Our calculations provide analytical estimates and are more accurate than those found in previous literature. We have validated our theory by means of numerical dynamical simulations and found a very good agreement. Finally we have designed a periodically tapered photonic crystal fiber which allows to achieve instability peaks at a large tunable frequency detuning from a given pump wavelength. Several higher-order resonant peaks are present but their gain and bandwidth are generally smaller than the one occurring at the smallest detuning. Such a system can be useful in quantum optical applications such as the efficient generation of entangled photon pairs in regions of frequencies far from the Raman gain peak.