Modulational instability in dispersion-kicked optical fibers
Abstract
We study, both theoretically and experimentally, modulational instability in optical fibers that have a longitudinal evolution of their dispersion in the form of a Dirac delta comb. By means of Floquet theory, we obtain an exact expression for the position of the gain bands, and we provide simple analytical estimates of the gain and of the bandwidths of those sidebands. An experimental validation of those results has been realized in several microstructured fibers specifically manufactured for that purpose. The dispersion landscape of those fibers is a comb of Gaussian pulses having widths much shorter than the period, which therefore approximate the ideal Dirac comb. Experimental spontaneous MI spectra recorded under quasi continuous wave excitation are in good agreement with the theory and with numerical simulations based on the generalized nonlinear Schrödinger equation.
pacs:
I Introduction
Modulational instability (MI) refers to a process where a weak periodic perturbation of an intense continuous wave (CW) grows exponentially as a result of the interplay between dispersion and nonlinearity. MI constitutes one of the most basic and widespread nonlinear phenomena in physics, and it has been studied extensively in several different physical systems like water waves, plasmas, and optical devices (1). For a cubic nonlinearity, as the one occurring in the nonlinear Schrödinger equation used to model optical fibers, the underlying physical mechanism can be understood in terms of four-wave mixing between the pump, signal and idler waves. However, the scalar four-wave interactions in a homogeneous fiber can be phase matched, and hence efficient, only in the anomalous group-velocity dispersion (GVD) regime. In the normal GVD regime, on the other hand, MI can occur in detuned cavities (2); (3), thanks to constructive interference between the external driving and the recirculating pulse. Alternatively MI with normal GVD can also arise in systems with built-in periodic dispersion (4); (5); (6); (7), among which dispersion oscillating fibers (DOFs) have recently attracted renewed attention (13); (8); (9); (10); (11); (12). In this case, phase matching relies on the additional momentum carried by the periodic dispersion grating (quasi-phase-matching). The occurrence of unstable frequency bands can then be explained using the theory of parametric resonance, a well-known instability phenomenon which occurs in linearized systems for which at least one parameter is varied periodically during the evolution (13). Up to now, most experimental investigations realised in optical fibers have been performed with basic sinusoidal (8); (9); (10); (11); (12) or amplitude modulated (14) modulation formats. In this work, on the other hand, we study a radically different periodic modulation of the GVD, in the form of a periodic train (or comb) of Dirac delta spikes. This is a fundamental and widespread modulation format, encountered in a variety of physical systems. In optics, delta combs have been exploited to model lumped amplification in long haul fiber optic transmission systems (18); (17), or to model the power extraction in soliton based fiber lasers (19). Moreover, comb-like dispersion-profiled fibers have been exploited to generate trains of solitons starting from a beat signal (20). At more fundamental level kicked systems are widely investigated as a paradigm for the emergence of chaos in perturbed Hamiltonian systems, with the delta-kicked rotor being the most renowned example (15). Its quantum version is described by a Schrödinger equation forced by a Dirac comb and has been extensively analyzed to study chaos in quantum systems (16). Recirculating fiber loops have been used to reproduce the quantum kicked rotor with an optical system, to study chaos and Anderson localization (21); (22), and to illustrate how an optical system can be used to mimic other physical systems that are more difficult to reproduce experimentally. In the same vein, we hope the experimental setup we propose in this paper could be used as an experimental platform to investigate such phenomena in the presence of nonlinearities, a topic of much current interest. Finally the approach that we propose to analyse MI in the fiber with delta-kicked GVD allows us, on one hand to enlighten the featuress of the parametric resonance that are not dependent on the specific format of the modulation, and on the other hand to compare and contrast the features of the ideal delta-kicked profile with other formats including non-ideal (physically realizable) kicking as well as widely employed profiles such as oscillating GVD.
The paper is organized as follows. In Section II we provide a simple argument allowing to determine the central frequencies of the unstable sidebands for general periodically modulated fibers. In Section III, we then use Floquet theory to analytically compute the width of the gain bands and as well as their maximum gain for dispersion-kicked fibers. In Section IV we investigate numerically the effect of the smoothing of the delta comb. In Section V, we describe the experimental set-up and we compare the experimental results with theory and numerical simulations based on the generalized nonlinear Schrödinger equation. We draw our conclusions in Section VI.
Ii Identifying the gain band central frequencies
Consider the NLSE
(1) |
We will assume the dispersion and the nonlinearity coefficient are of the form
(2) |
where and are periodic functions of period such that , and .
Let be a stationary solution of (1). We consider a perturbation of in the form , where the perturbation satisfies . Inserting this expression in (1), and retaining only the linear terms we find
(3) |
Writing , with and real functions, we obtain the following linear system:
Finally, taking the Fourier transform of this system in the time variable , leads to
(4) |
where we used the definiton . Note that this is a Hamiltonian dynamical system in a two-dimensional phase plane with canonical coordinates . Analyzing the linear (in)stability of the stationary solution therefore reduces to studying the solutions to (4) for each . Since the coefficients in the equation are -periodic with period , Floquet theory applies. This amounts to studying the linearized evolution over one period , to obtain the Floquet map which in the present situation is the two by two real matrix defined by . As a result . Note that necessarily has determinant one, since it is obtained by integrating a Hamiltonian dynamics, of which we know that it preserves phase space volume. As a consequence, if is one of its eigenvalues, then so are both its complex conjugate and its inverse . This constrains the two eigenvalues of considerably: they are either both real, or lie both on the unit circle. Now, the dynamics is unstable only if there is one eigenvalue satisfying , in which case both eigenvalues are real. We will write for the two eigenvalues of . We are interested in studying the gain, that is
(5) |
as a function of , and . It measures the growth of . The gain vanishes if the two eigenvalues lie on the unit circle. A contour plot of the gain in the plane, for the case of the delta comb dispersion modulation that is the main subject of this paper, can be found in Fig. 2. The regions where the gain does not vanish are commonly referred to as Arnold tongues. We will explain below that, whereas their precise form depends on the choice of , the position of their tips does not.
Since the system (4) is not autonomous, it cannot be solved analytically in general. Nevertheless, the above observations will allow us to obtain some information about its (in)stability for small , , and valid for all perturbations , whatever their specific form.
To see this, we first consider the case . It is then straightforward to integrate the system (4). The linearized Floquet map is then given by
(6) |
where
(7) |
Here (normal average dispersion), since we restrict our investigations to the defocusing NLS. Note that the matrix has determinant equal to , as expected. The eigenvalues of can be readily computed as
(8) |
What will happen if we now switch on the interaction terms and ? It is then no longer possible, in general, to give a simple closed form expression of the solution to (4), which is no longer autonomous, and hence of the linearized Floquet map . Nevertheless, we do know that, for small , the eigenvalues of must be close to the eigenvalues . We then have two cases to consider.
Case 1. .
Now , they are distinct, and they both lie on the unit circle, away from the real axis. They then must remain on the unit circle under perturbation since, for the reasons explained above, they cannot move into the complex plane away from the unit circle.
Consequently, in this case, the stationary solution is linearly stable under a sufficiently small perturbation by and , and this statement does not depend on the precise form of or of . In fact, with growing and/or , the two eigenvalues will move along the unit circle until they meet either at or at for some critical value of the perturbation parameters. Only for values of the latter above that critical value can the system become unstable. A pictorial description of this situation is shown in the left hand side of Fig. 1.
Case 2. , . Now is a doubly degenerate eigenvalue of . Under a small perturbation, the degeneracy can be lifted and two real eigenvalues can be created, one greater than one, one less than one in absolute value. The system has then become unstable! Of course, it will now depend on the type of perturbation whether the system becomes unstable, remains marginally stable (the two eigenvalues don’t move at all, but stay at or ), or becomes stable (the two eigenvalues move in opposite directions along the unit circle). A pictorial description of this situation is shown in the right hand side of Fig. 1. For the Dirac comb modulation of , which is our main object of study in this paper, the details are given in the next section.
In conclusion, examining (7), one sees that only if , where
(9) |
can an infinitely small Hamiltonian perturbation of lead to an unstable linearized dynamics near the fixed points considered. These values of therefore correspond to the tips of the Arnold tongues, that is, to the positions of the (centers of) the unstable sidebands of the defocusing NLS under a general periodic perturbation . This is illustrated for a Dirac comb modulation of the GVD in Fig. 2. One also observes in that figure that, for a value of close to some , the system becomes unstable only for a small but nonzero critical value of , that we shall compute below for the Dirac delta comb GVD.
Equation (9) was derived in (13) by appealing to the theory of parametric resonance and Poincaré-Lindstedt perturbation theory. Our argument above is elementary and shows in a simple manner that the resonant frequencies do not at all depend on the form of or . Note that, if and , a case considered in (9)-(10), the system (4) is equivalent to the equation of a harmonic oscillator of (spatial) frequency , sinusoidally modulated with period . In that case the system leads to a Mathieu equation for which it is known that resonance occurs when the period of the modulation is a integer multiple of the half (spatial) period of the oscillator, which is .
Additional physical insight can be obtained by expanding Equation (9) for small power, i.e. assuming . At zero order we recover the well known quasi-phase-matching relation (5); (8); (17)
(10) |
Equation (10) entails the conservation of the momentum, made possible thanks to the virtual momentum carried by the dispersion grating, of the four wave mixing interaction between two photons from the pump, going into two photons in the symmetric unstable bands at lower (Stokes) and higher (antiStokes) frequencies with respect to the pump.
Iii Calculation of the Modulational Instability gain bands: Dirac comb
We now turn our attention to the computation of the gain , in particular for values of close to the resonant frequencies. We concentrate on the special case where the GVD is a Dirac delta comb:
(11) |
Since in the rest of this paper, , we will drop it from the notation. To compute the gain, we need to compute the linearized dynamics and determine the behaviour of its eigenvalues in the neighbourhood of and in the -plane.
In this case the linearized Floquet map is easily seen to be explicitly given by
(12) |
where is defined by Equation (6), but now with , and
(13) |
The characteristic polynomial of is given by
so that the eigenvalues of (12) can be computed explicitly as:
(14) |
with
and .
A Taylor expansion of about yields
(15) |
where
(16) | ||||
(17) |
The dependence in (not in ) entails that the sign of the kick has no incidence in this regime, i.e. assuming .
Formula (15) shows that is a saddle point for . If is even, occurs close to , and if is odd, close to . More precisely,
from which we can find an estimate of the gain amplitude and of the bandwidth near the tips of the tongue at , as:
(18) |
(19) |
Note that the threshold value for above which instability occurs can be read off from the above by setting which corresponds to
This confirms again, as expected, that an arbitrary small will generate instability right at . In Fig. 2a we show an example of the analytically computed MI gain, showing the first two Arnold tongues. As can be seen, for a small enough strength of perturbation, let’s say , the approximation (19) gives a good estimate of the width of the parametric resonance (see red curves). This situation is detailed further in Figs. 2b,c, showing a section for and , respectively.
Finally, a straightforward calculation gives the asymptotic behaviour of the gain at for large and fixed, that is
(20) |
with whenever , and otherwise.
In Fig. 3, we show an example of the analytically computed MI gain at as a function of . We compare it to the approximation (20), which is very accurate, even for small (see red circles). Note in particular that the oscillating behaviour of the gain is well captured by (20) which, for large enough and small, can be approximated by
Summing up. It is clear from the above that, precisely at the values , which only depend on and on , but not on the precise form of , any small perturbation can create an instability and hence a gain. At frequencies near these particular values, a minimal threshold strength of is needed to create an instability. This minimal value, and even the fact that an instability is indeed generated, does depend on the precise form of . For the Dirac comb the explicit expression for the gain in this regime can be read off from (20).
Iv Approximations of the delta function
In order to shed light on the dependence of the gain on the shape of the periodic modulation, and also with an eye towards the experimental realization of the Dirac comb fiber, we now analyze what happens when the Dirac comb is approximated by a train of physically realizable “kicks”. We thus consider a smoothened Dirac comb described by
(21) |
where we normalize the positive function in order to have . For a rectangular pulse of width , we get
(22) |
For a gaussian function , the maximum amplitude of the kick can be calculated as
that in the limit gives
Note that, in these models, we have ()
Hence corresponds to a rather symmetric situation where is close to the midpoint between and , so that fluctuates symmetrically about its average value. Whereas corresponds to a very asymmetric situation where has a large abrupt peak. The parameter therefore controls the shape of the GVD modulation at fixed and (or ).
As shown in the previous section, by changing the shape of the kick, we do not change the frequency of the parametric resonances. The smoothing of the delta function nevertheless does modify the characteristics of the MI by changing the value of the gain, as we now illustrate by computing the gain numerically at the resonant frequencies . An example of how the changing shape of the modulation modifies the first parametric resonance is illustrated in Fig. 4, that shows the gain at fixed and , as a function of the peak amplitude (or, equivalently, the width ) of the kicks. We make the following observations. First, a good approximation of the gain given by the Dirac comb is obtained for , both for the rectangular and gaussian pulses. Second, for , the gain of the square pulse modulation is zero, as expected, since we are then in the limit case of a constant modulation (and normal GVD). Third, it is apparent that the Dirac comb gives the highest possible gain, for a fixed area of the kicks and fixed and . Finally, it is interesting to note that a sinusoidal modulation , with the same value of , gives a gain close to one half with respect to the delta case. Indeed for a sinusoidal modulation, it has been shown that (see Equation (7) from Ref. (10))
(23) |
By expanding Equation (23) for small , we get
(24) |
at first order in . In conclusion, a large concentrated perturbation of the GVD about its average enhances the MI gain.
V Experimental results
It is well known that in homogeneous fibers, the GVD depends on the diameter of the fiber. One can therefore modulate the GVD by modulating the diameter of the fibers as a function of , as in (8); (9); (10). We manufactured three different microstructured optical fibers modulated by a series of Gaussian pulses to approximate the ideal Dirac delta comb studied in Section III. The change of their outer diameters along the fiber is represented in Fig. 5(a). As can be seen in the inset, their diameters have a gaussian shape with a standard deviation , which is the same for all three fibers, and very small ( m) compared to the period of the comb (10 m). Hence we can write
where . The three fibers have a very similar minimum diameter , while their maximum values are different. We have 172 , 207 and 240 for fibers labelled A, B and C, respectively, corresponding to . To understand how the two experimental parameters and control the quality of the approximation of the delta function on the one hand, and the value of on the other hand, we proceed as follows. First, , so that . A first order Taylor expansion of about yields
Comparing this to (2) we find
and
Hence, with the notation of Section IV, . This corresponds to proving that these Gaussian pulses should induce a very similar parametric gain compared to ideal Dirac delta functions (see Fig. 4). Furthermore, the height of the Gaussian pulse controls , which will allow us to investigate the impact of this parameter on the first MI side lobe gain, as it was done in the theoretical study and illustrated in Fig.2.
For all three fibers, the ratio of the diameter of the holes over (the pitch of the periodic cladding) is assumed to be constant along the fiber and estimated to about 0.48 from scanning electron microscope images. The diameter variations of the fibers are proportional to those of the pitch, with corresponding to the minimum value of the diameter (, green line in Fig. 5(a)) and and , blue, black and red curves respectively, for the maximum values. As an example, the Group Velocity Dispersion (GVD) curve corresponding to the minimum pitch value has been calculated from Ref.(24) and is represented in Fig. 5(b) as a green curve. Its Zero Dispersion Wavelength (ZDW) is located at 1055 nm while those of the GVD curves corresponding to the maximum values of the diameters of fibers A, B and C are red-shifted to 1110 nm, 1136 nm and 1168 nm respectively (Fig. 5(b)). In order to give another illustration of the large dispersion variations induced in these fibers by varying their diameters, the maximum GVD values for fibers A, B and C have been calculated at a fixed wavelength (1052.5 nm) and compared to the background value. As can be seen in Fig. 5(b), an increase of the diameters by a factor of only 1.27, 1.45 and 1.75, leads to a one order of magnitude improvement on the GVD values: 21, 29 and 35 respectively. Under such large variation of the fiber diameter, the GVD can no longer be considered as proportional to the pitch value, as was the case in Ref.(9) for instance. This is illustrated in Fig.5(c), where the evolution of the GVD calculated at 1052.5 nm as a function of the fiber diameter is represented. It can be seen to be well approximated by an affine function in the range between 140m and 180 m, but not beyond. As a consequence, the shape of the GVD variations will be slightly different from the one of the diameter, specifically for fiber C. However, we checked numerically that this can be considered as relatively weak distortions that do not significantly impact the gain of the MI process. We can still consider that the key parameters remain the different heights in GVD of the Gaussian-like pulses in fibers A, B and C. In order to get a more complete picture of the impact of the fiber diameter variations on its guiding properties, the variation of the nonlinear coefficient is plotted as a function of the fiber diameter in Fig. 5(d). The most important feature to note here is that the amplitude of variation is much smaller, and only the same order as the one of the diameter itself. Hence, these variations are more than one order of magnitude lower than those of the GVD and we have checked numerically that their impact on the MI process is negligible. Consequently, we can infer that these fibers represent a good prototype to validate our theoretical investigation in the previous sections, where only longitudinal GVD variations have been taken into account.
The experimental setup is schematized in Fig. 6. The pump system is made of a continuous-wave tunable laser (TL) diode that is sent into an intensity modulator (MOD) in order to shape 2 ns square pulses at 1 MHz repetition rate. They are amplified by two Ytterbium-doped fiber amplifiers (YDFAs) at the output of which two successive tunable filters are inserted to remove the amplified spontaneous emission in excess around the pump. These quasi-CW laser pulses have been launched along the birefringent axis of the fibers. The pump peak power has been fixed to 6.5 W and the pump wavelength at 1052.5 nm for fiber A. The output spectrum recorded at its output is represented as a blue curve in Fig. 7 (a). Two MI side lobes, located at +/-4.8 THz appear on both sides of the pump. These experimental results have been compared with numerical simulations performed by integrating the generalized nonlinear Schrödinger equation. We used the GVD, and variations calculated from the measured diameter values (see Figs. 5). Other parameters are extracted from experiments and are listed in the caption of Fig. 7. Note that we checked that except the longitudinal GVD variations, all other parameters can be assumed to be constant and equal to the average values. As can be seen in the blue curve in Fig. 7(b), two symmetric MI side lobes also appear in the simulated spectrum, in a very good agreement with experiments. Their positions have been compared with the predictions of Equation 9, represented by green dashed lines in Fig. 7(b) (calculated with ). An excellent agreement is also obtained. In order to show that the MI gain is larger when the weight of the Dirac delta function is increased, we performed similar experiments in fibers B and C where the areas of the Gaussian pulses are larger than in fiber A. However, in experiments, due to the fact that the Dirac comb has been approximated with a series of Gaussian functions, changing their amplitudes also modifies the average value of the dispersion. As a consequence, MI side lobes would be generated at different frequency shifts. In order to keep constant the position of the MI side lobes, and then provide a correct comparison with the theoretical study, one have to take care to keep the average GVD value constant in all the fibers. To do so experimentally in fibers B and C, we slightly tuned the pump wavelength until the first MI side lobe is generated at 4.8 THz as in fiber A. As can be seen in Fig. 7(a) (red and black curves), the position of the first MI side lobe in fibers B and C is indeed located at that frequency by tuning the pump wavelength to 1061.8 nm and 1067 nm, respectively. We can therefore consider that the average GVD values are very similar in the three fibers and hence that only the areas of the Gaussian functions, i.e. the equivalent of the Dirac weights, vary. The amplitudes of the first MI side lobes generated in fibers B and C are indeed larger compared to fiber A, as predicted by the theory. This is in pretty good agreement with numerical simulations (Fig.7(b)), where the same procedure was used. We found that the average values in fibers B and C are and respectively. The small discrepancy between these values is attributed to spurious longitudinal fluctuations arising during the drawing process. Indeed, as can be seen in Fig. 5(a), the background over which Gaussian pulses are superimposed is not perfectly flat, and in fibre C, it is not horizontal. To counterbalance these imperfections, it was necessary to adjust the average dispersion values. We checked that with a series of perfect Gaussian pulses superimposed on a flat and horizontal background, the same average GVD value would be obtained. Furthermore, we can note that in fibers B and C, additional MI side lobes are generated due to the periodic modulation of the GVD (labelled MIi in 7(a), up to 5 in fiber B). Their positions are also well predicted by numerical simulations and by Equation 9 (green lines in 7(b)). This excellent agreement confirms that their positions indeed scale approximately as , being the side lobe order, that is the typical signature of the MI process occurring in dispersion oscillating fibers. It was already reported experimentally in Refs.(9); (14), with a sinusoidal variation of the GVD, modulated in amplitude or not and it is now illustrated in this paper with a Dirac delta comb. Moreover, in fibers B and C two symmetric side lobes that are not predicted by the theory appear around the pump at about 2.15 THz (labelled spurious side lobes in Fig.7(a)). They result from a non-phase matched four-wave mixing process involving the pump, the first and second MI side lobes. The energy conservation relation involving these waves predicts a frequency shift of 2.2 THz from the pump for the fourth wave, that is in good agreement with the shift of 2.15 THz measured experimentally.
Vi Conclusions
Modulation instability has been investigated theoretically and experimentally in dispersion kicked optical fibers. An analytical expression of the parametric gain has been obtained allowing to predict the behavior of the MI process in such fibers. Specifically, it was shown that increasing the weights of the Dirac functions leads to larger MI gains for the first MI side lobe. We exploit the fact that the Dirac delta comb can be well approximated by a series of short Gaussian pulses in order to perform an experimental investigation using microstructured optical fibers. We then experimentally report, for the first time to our knowledge, multiple MI side lobes at the output of these dispersion kicked optical fibers. We demonstrate that they originate from the periodic variations of the dispersion. We also validate experimentally that increasing the height of the modulation leads to a larger gain for the first MI side lobe. This illustrates that optical fibers constitute an interesting platform to realize experimental investigations of fundamental physical phenomena.
Acknowledgments
The present research was supported by the Agence Nationale de la Recherche in the framework of the Labex CEMPI (ANR-11-LABX-0007-01), Equipex FLUX (ANR-11-EQPX-0017), and by the projects TOPWAVE (ANR-13-JS04-0004), FOPAFE (ANR-12-JS09-0005) and NoAWE (ANR-14-ACHN-0014).
References
- V. E. Zakharov and L. Ostrovsky, Physica D 238, 540 (2009).
- M. Haelterman, S. Trillo, and S. Wabnitz, Opt. Lett. 17, 745 (1992).
- S. Coen and M. Haelterman, Phys. Rev. Lett. 21, 4139 (1997).
- J. Bronski and J. N. Kutz, Opt. Lett. 21, 937 (1996).
- N. J. Smith and N. J. Doran, Opt. Lett. 21, 570 (1996).
- F. Kh. Abdullaev, S.A. Darmanyan, A. Kobyakov, and F. Lederer, Phys. Lett. A 220, 213 (1996).
- F. Kh. Abdullaev and J. Garnier, Phys. Rev. E 60, 1042 (1999).
- M. Droques, A. Kudlinski, G. Bouwmans, G. Martinelli, and A. Mussot, Opt. Lett. 37, 4832 (2012).
- M. Droques, A. Kudlinski, G. Bouwmans, G. Martinelli, A. Mussot, A. Armaroli, and F. Biancalana, Opt. Lett. 38, 3464 (2013).
- M. Droques, A. Kudlinski, G. Bouwmans, G. Martinelli, and A. Mussot, Phys. Rev. A 87, 013813 (2013).
- C. Finot, J. Fatome, A. Sysoliatin, A. Kosolapov, and S. Wabnitz, Opt. Lett. 38, 5361 (2013).
- M. Conforti, A. Mussot, A. Kudlinski, and S. Trillo, Opt. Lett. 39, 4200 (2014).
- A. Armaroli and F. Biancalana, Opt. Express 20, 25096 (2012).
- Copie, FranÃ§ois and Kudlinski, Alexandre and Conforti, Matteo and Martinelli, Gilbert and Mussot, Arnaud, Opt. Express 4, 3869 (2015).
- B. V. Chirikov, Phys. Rep. 52, 263 (1979).
- G. Casati, B. V. Chirikov, F. M. Izraelev, and J. Ford, Lecture Notes in Physics 93, 334 (1979).
- F. Matera, A. Mecozzi, M. Romagnoli, and M. Settembre, Opt. Lett. 18, 1499 (1993).
- N. J. Smith, K. J. Blow, and I. Andonovic, J. Lightwave Techol. 10, 1329 (1992).
- S. M. J. Kelly, Electron. Lett. 28, 806 (1992).
- S. V. Chernikov, J. R. Taylor, and R. Kashyap, Opt. Lett. 19, 539 (1994).
- B. Fisher, A. Rosen, and S. Fishman, Opt. Lett. 24, 1463 (1999).
- S. Atkins, A. Rosen, A. Bekker, and B. Fisher, Opt. Lett. 28, 1463 (2003).
- A. H. Nayfeh and D. T. Mook, Nonlinear Oscillations (Wiley, New York, 1979).
- K.Saitoh and M. Koshiba, Opt. Express 13, 267 (2005).