Soliton Content of Fiber-Optic Light Pulses

Abstract

This is a review of fiber-optic soliton propagation and of methods to determine the soliton content in a pulse, group of pulses or a similar structure. Of central importance is the nonlinear Schrödinger equation, an integrable equation that possesses soliton solutions, among others. Several extensions and generalizations of this equation are customary to better approximate real-world systems, but this comes at the expense of losing integrability. Depending on the experimental situation under discussion, a variety of pulse shapes or pulse groups can arise. In each case, the structure will contain one or several solitons plus small amplitude radiation. Direct scattering transform, also known as nonlinear Fourier transform, serves to quantify the soliton content in a given pulse structure, but it relies on integrability. Soliton radiation beat analysis does not suffer from this restriction, but has other limitations. The relative advantages and disadvantages of the methods are compared.

1. Introduction

Light pulses in optical fibers carry the bulk of all telecommunication today. In comparison to electrical pulses in cables, fiber-optic transmission is vastly superior, mostly due to two fundamental advantages: light in fibers suffers extremely low power loss, and fibers provide an extremely wide bandwidth. For the loss, a rough number is $0.2\mathrm{dB}/\mathrm{km}$. This implies that even after 100 km, still, 1% of the launch power is left, which is unrivaled. The available bandwidth can be estimated as $30\mathrm{THz}$ (maybe $50\mathrm{THz}$ if one accepts slightly higher loss), which is several orders of magnitude better than anything that can be achieved with electronics.

Further comparison of optical and electrical transmission brings us to an imperfection shared by both: in either type of conduit, group velocity depends on frequency; as any signal requires a certain bandwidth, transmitted signals suffer from distortion due to group velocity dispersion. This is a linear distortion, which means that it can be perfectly compensated by adding an element of opposite dispersion.

In one respect, however, optical fibers are very different from electrical cables. The response of glass to light is nonlinear, whereas that of copper to current is not (within reasonable limits). Therefore, the physics of data transmission with light pulses through optical fibers is fundamentally different in that it involves a nonlinear response of the material to the signal [1,2]. Nonlinear distortions are not as easily compensated as linear ones. This has led to a widely-held conception that optical power in fibers must always be kept low enough so that nonlinear distortions are avoided.

To this argument, one can give two responses: One, nonlinearity can actually be put to good use, even be used to an advantage, when the concept of solitons is adopted. It has already been demonstrated even in a commercial setting. Two, the ever-growing demand for data-carrying capacity brings us close to a capacity crunch, unless means are found in just a few years to improve technology. It may become necessary to look at solitonic formats again.

2. The Nonlinear Schrödinger Equation and Some of Its Solutions

The glass of optical fiber is a dispersive material. With the index being a function of wavelength or optical frequency, i.e., $n=n\left(\omega \right)$, the group velocity is also frequency dependent. The glass is also a nonlinear material. This chiefly means that the refractive index follows the instantaneous intensity, i.e., $n=n\left(I\right)$. Across the duration of a short pulse of light, the index experiences a slight modulation, which gives rise to intensity-dependent phase shifts known as self-phase modulation [3].

A propagation equation for light pulses in fibers must contain these two influences, at least in leading order. The equation that fits this description is known as the nonlinear Schrödinger equation (NLSE). It plays an absolutely central role in all of the nonlinear fiber optics. It describes the evolution of the pulse amplitude envelope (the fast oscillation at some central optical frequency ${\omega }_0$ is removed) in a frame of reference comoving with the traveling light, i.e., with the group velocity at ${\omega }_0$. Written in physical quantities, it is given by:

$$i\frac{\partial }{\partial z}A-\frac{{\beta }_2}{2}\frac{{\partial }^2}{\partial t^2}A+\gamma {\left|A\right|}^{2}A=0.$$

(1)

Here, $A=A(z,t)$ is the complex amplitude with z the position along the fiber and t the retarded time. ${\beta }_2$ is the coefficient of group velocity dispersion, and $\gamma$ is the coefficient of nonlinearity. Depending on the situation, it may become necessary to add corrective terms for loss, higher-order dispersion or other nonlinear effects, as we shall see in Section 3.1. The equation received its name due to its formal similarity to the Schrödinger equation of quantum mechanics fame: The quantum mechanical version typically describes how a wave function spreads out in space as time goes by; here, the equation describes how a short pulse gets broadened temporally as it propagates down the fiber. Therefore, time and space coordinates switch roles. The potential, in the nonlinear case, comes from the self-phase modulation, which is captured in the nonlinear term.

The effort to find solutions to this equation benefited from the fact that the NLSE is integrable. Integrable equations have very special mathematical properties like an infinite number of preserved quantities. A first solution was found by Vladimir E. Zakharov and Aleksei B. Shabat in 1972 [4] using the inverse scattering technique, developed by Gardner et al. in 1967 [5] (see Section 4.1). It exists for a particular combination of algebraic signs of the two parameters; with $\gamma$ always positive, one would have the correct signs for ${\beta }_2<0$, i.e., for what is known as anomalous dispersion. Two years later, Junkichi Satsuma and Nubuo Yajima [6] treated the pertaining initial-value problem. Between them, these two papers already form a solid base of the understanding of solitons. In the time between those two, Akira Hasegawa and Frederick Tappert suggested that the NLSE was applicable to pulse propagation in optical fibers, so that fiber-optic solitons would exist [7]. At that time, optical fibers were a novelty and not yet technologically matured. Attempts to transmit light pulses for data transmission suffered from dispersive distortions. It was therefore a bold proposal that the use of solitons as signaling pulses might overcome dispersive broadening, in other words that nonlinearity could be exploited to cancel distortions arising from a linear mechanism.

Of course, the proposal required as a very first step that the existence of fiber-optic solitons be corroborated experimentally. In the anomalously dispersive regime, i.e., at wavelengths longer than ca. $1.3\mathsf{\mu }\mathrm{m}$, they were still too lossy for any meaningful test. By 1980, this situation had improved, and Linn F. Mollenauer, Roger H. Stolen and James P. Gordon could experimentally verify the existence of solitons [8].

2.1. Soliton Solutions

2.1.1. The Fundamental Soliton

In the anomalous dispersion regime, the NLSE Equation (1) supports a stable solution known as the fundamental soliton:

$$\begin{array}{ccc}\hfill A(z,t)& =& \sqrt{P_0}\mathrm{sech}\left(\frac{(t-t_{\mathrm{s}})-\mathrm{\Omega }{\beta }_2(z-z_{\mathrm{s}})}{T_0}\right)\hfill \\ & \times & \exp \left[i\frac{\gamma P_0}{2}(z-z_{\mathrm{s}})-i\mathrm{\Omega }\left((t-t_{\mathrm{s}})-\frac{1}{2}\mathrm{\Omega }{\beta }_2(z-z_{\mathrm{s}})\right)+i{\phi }_{\mathrm{s}}\right].\hfill \end{array}$$

(2)

Here, $P_0$ is the peak power; $T_0$ is the pulse duration; $\mathrm{\Omega }$ is the deviation of the soliton’s center optical frequency from that of the frame of reference, ${\omega }_0$. $t_{\mathrm{s}}$, $z_{\mathrm{s}}$ and ${\phi }_{\mathrm{s}}$ are the initial values of center time, start position and phase offset, respectively. As long as a single soliton is considered, frequency offset and initial values may be set to zero without loss of generality so that a considerably simplified version:

$$A(z,t)=\sqrt{P_0}\mathrm{sech}\left(\frac{t}{T_0}\right)\exp \left[i\frac{\gamma P_0}{2}z\right]$$

(3)

remains. In either case, amplitude and duration are coupled and must fulfil:

$$P_0{T}_0^2=\frac{|{\beta }_2|}{\gamma }.$$

(4)

The propagation of a fundamental soliton with $\mathrm{\Omega }=t_{\mathrm{s}}=z_{\mathrm{s}}={\phi }_{\mathrm{s}}=0$ is shown in Figure 1. As the temporal evolution (sech term) does not contain z, the temporal shape does not change during propagation. Without nonlinearity, the shape would change appreciably after one dispersion length $L_{\mathrm{D}}=T_0^2/\left|{\beta }_2\right|$. The other relevant length scale is the nonlinearity length $L_{\mathrm{NL}}=1/\left(\gamma P_0\right)$ after which self-phase modulation becomes appreciable. As Equation (4) shows, for the soliton $L_{\mathrm{D}}=L_{\mathrm{NL}}$. This equality expresses the balance between both.

Among the conserved quantities of the integrable equation are the soliton energy:

$$E_{\mathrm{sol}}=2P_0{T}_0$$

(5)

and the soliton center frequency $\mathrm{\Omega }$. As will be described below (Section 4.1), these two values are found from inverse scattering in the form of imaginary and real parts of a complex eigenvalue. We point out that in a dispersive medium, a nonzero frequency translates to a relative motion, so that the frequency is usually referred to as the velocity. The number of solitons is also preserved, even in cases when there is more than one.

2.1.2. Radiation and Higher-Order Solitons

To create a soliton, one can launch a sech-shaped pulse as in Equation (3) for $z=0$ into a fiber with parameters conforming to Equation (4). If, however, the parameters are only approximately right, we can do as follows: In the spirit of [6] and as also detailed in [9] (but they call A what we call N as in [1,2,10]), we write the initial condition as:

$$A(0,t)=N\sqrt{P_0}\mathrm{sech}\left(\frac{t}{T_0}\right).$$

(6)

where the pulse is scaled by a factor $N\ge 0$, $N\in \mathbb{R}$ called the soliton order. The energy is then:

$$E=N^2{E}_{\mathrm{sol}}.$$

(7)

Of course, at $N=1$, there is the fundamental soliton, but at non-integer N, there is also a part of the pulse energy that is not part of the soliton energy. This balance is a linear wave; in the soliton context, it is called radiation. Linear waves are subject to dispersive broadening, so that after a sufficiently long propagation, radiation will be dispersed away from the soliton, and the soliton itself emerges. Until that happens, the different phase evolution of soliton and radiation creates a beat note, visible as a slowly-decaying beat pattern. This is shown for $N=1.2$ in Figure 2 (left), where lengths are scaled to $z_0=(MM\pi /2)L_{\mathrm{D}}$. Evaluation of the beat pattern can yield soliton parameters, see Section 4.2.

There is a threshold value of N to generate a soliton. As long as $N<\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, all energy goes into radiation. At $N=1$, the radiative part dips to zero. At $N=\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, it reaches the same energy as at $N=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ again, and that suffices to generate a second soliton. At $N=2$, there is a combination of two solitons without any radiative part. This structure, for which an explicit expression was already given in [6], is called an $N=2$ soliton (Figure 2 right). More generally speaking, for any $N\in \mathbb{N}$, $N>1$, one has a so-called higher-order soliton of order N in the pure form, i.e., without radiation. An explicit expression for the $N=3$ case solitons was found in [11]. The power profile of all higher-order solitons oscillates with period $z_0$. For all non-integer values of N, there is a radiative contribution.

Figure 3 gives an overview of this situation. In its upper part, the pulse energy according to Equation (7), shown in units of $E_{\mathrm{sol}}$, is represented by the dashed parabola. Beginning at $N=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, the first soliton appears. Its energy then rises linearly and passes through unity at $N=1$. Similarly, at other half-integer values of N, more solitons begin, each of them contributing energy to a cumulative value shown in blue. The latter is tangent to the parabola at integer N values; here, the individual soliton energies are also integer in $E_{\mathrm{sol}}$ units. They take the first N values of the sequence $1,3,5,7,\cdots$, which add up to $N^2$. For all non-integer N, there is a gap between the cumulative solitonic and the total energy, which represents the radiative part of the energy. Radiation energy is shown in the figure’s lower part on an inverted, magnified scale (green curve). Whenever, at some half-integer N, it reaches $E_{\mathrm{rad}}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.E_{\mathrm{sol}}$, the next soliton is created. Always the integer number closest to N determines the number of solitons involved, while the fractional part determines the amount of radiation created along with them.

2.1.3. Soliton Interaction

When two or more pulses are launched in rapid succession, each of them is affected by the presence of the others due to the index nonlinearity. In an integrable system, they cannot merge (the number is preserved), but they can be set in relative motion as if there were interaction forces. Gordon investigated the soliton-soliton interaction [12] and found that an effective force would decay exponentially with growing separation between pulses and vary sinusoidally with their relative phase. This means that there can be both attraction (in phase) or repulsion (opposite phase). This prediction was tested experimentally in [13] and was found to be fully correct, with the only caveat that in the attractive case, the first close encounter (collision) leads to processes of higher order not captured in the NLSE.

Note that this view of interaction forces between particles is metaphorical and justifiable only if both are well separated. In a nonlinear system, the superposition principle does not hold; therefore, two sech-shaped pulses in close proximity to each other are not really two identical, but overlapping solitons, but form a somewhat more involved two-soliton compound [14]. For that, energy and velocity are preserved quantities.

2.2. Breather Solutions

There is a different family of solutions of the NLSE, which is best explained by starting with the continuous wave (cw) case. The NLSE is solved by the cw ansatz:

$$A=\sqrt{P_0}{\mathrm{e}}^{i\gamma P_{0}z}$$

(8)

with constant power $P_0$. It turns out that this solution is stable for normal (${\beta }_2>0$) and unstable for anomalous dispersion (${\beta }_2<0$). The instability implies growth of small perturbations of the cw; this is known as modulation instability [1]. The frequency of maximum gain is offset by $\pm {\omega }_{\max }$ from the carrier frequency, with:

$${\omega }_{\max }=\sqrt{\frac{2\gamma P_0}{|{\beta }_2|}}.$$

(9)

If the cw is perturbed by white noise, there is always energy in these Fourier components, and a ripple with frequency ${\omega }_{\max }$ will begin to grow with gain:

$$g_{\max }=2\gamma P_0.$$

(10)

The further fate was widely appreciated only a couple of years ago even though the mathematical results had all been worked out as early as 1986 [15]. There is a family of solutions of the NLSE characterized by an infinitely extended cw background and some modulation on top. Its general form can be written as [16]:

$$A(a,z,t)=\sqrt{P_0}\left[1+M(a,z,t)\right]\exp \left(i\gamma P_{0}z\right).$$

(11)

where $M(a,z,t)$ describes the modulation term. In the case of what is now known as the Akhmediev breather, it takes the form:

$$M(a,z,t)=\frac{2(1-2a)\cosh \left[b\left(a\right)\gamma P_{0}z\right]+ib\left(a\right)\sinh \left[b\left(a\right)\gamma P_{0}z\right]}{\sqrt{2a}\cos \left[{\omega }_{\mathrm{mod}}t\right]-\cosh \left[b\left(a\right)\gamma P_{0}z\right]}.$$

(12)

with parameters $0<a<\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ and $b\left(a\right)=\sqrt{8a-16a^2}$. The modulation frequency is ${\omega }_{\mathrm{mod}}={\omega }_{\mathrm{c}}\sqrt{1-2a}$, and:

$${\omega }_{\mathrm{c}}=\sqrt{\frac{4\gamma P_0}{|{\beta }_2|}}\ge \left|{\omega }_{\mathrm{mod}}\right|$$

(13)

is the frequency range where gain is possible. Gain maximum occurs at $a=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.$ with ${\omega }_{\max }={\omega }_{\mathrm{mod}}/\sqrt{2}$, and Equations (9) and (10) are recovered. This particular case is illustrated in Figure 4. From a cw background, the modulation grows, until at $z=0$, it reaches a maximum amplitude. Here, one has a periodic train of pulses. During further evolution, the modulation decays until for $z\to \infty$, the cw is recovered.

The first experiments to demonstrate the Akhmediev breather in optical fibers were described in [17] and in more detail in [18]. The process was initiated either from random perturbation [17] or from a suitable weak modulation as suggested in [19] and performed experimentally in [18].

If we let $a\to \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, this implies that ${\omega }_{\mathrm{mod}}\to 0$, and the temporal separation between the pulses in the train diverges. The modulation then takes the form:

$$M(z,t)=-\frac{4(1+i2\gamma P_{0}z)}{1+{\omega }_{\mathrm{c}}^2{t}^2+4{\left(\gamma P_{0}z\right)}^2}.$$

(14)

This is a structure that is localized in both time and space, known as the Peregrine soliton [20]. It is shown in Figure 5 (left); a first experimental demonstration was given in [21]. A further member of this family of NLSE solutions was studied in [22] and is known as the Kuznetsov–Ma soliton. Mathematically, that solution is described by the modulation term:

$$M(a,z,t)=\frac{2(1-2a)\cos \left[\left|b\left(a\right)\right|\gamma P_{0}z\right]-i\left|b\left(a\right)\right|\sin \left[\left|b\left(a\right)\right|\gamma P_{0}z\right]}{\sqrt{2a}\cosh \left[\left|{\omega }_{\mathrm{m}}\right|t\right]-\cos \left[\left|b\left(a\right)\right|\gamma P_{0}z\right]}$$

(15)

with a parameter range of $a>\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$. It consists of a periodically-recurring peak; see Figure 5 (right).

3. Extensions to the Nonlinear Schrödinger Equation

3.1. Generalized Nonlinear Schrödinger Equation

The NLSE neglects linear loss, approximates dispersion by the leading group velocity dispersion term and ignores nonlinear effects other than the intensity dependence of the refractive index. Depending on the situation and/or the required accuracy of the calculated results, terms must be added to include what is not represented. This leads to a generalized NLSE, which may take the form:

$$\begin{array}{ccc}\hfill i\frac{\partial }{\partial z}A-\frac{{\beta }_2}{2}\frac{{\partial }^2}{\partial t^2}A& -& \underset{\mathrm{higher-order}\mathrm{dispersion}}{\underbrace{i\frac{{\beta }_3}{6}\frac{{\partial }^3}{\partial t^3}A+\frac{{\beta }_4}{24}\frac{{\partial }^4}{\partial t^4}A\cdots }}+\gamma {\left|A\right|}^{2}A\hfill \\ & +& \underset{\mathrm{self-steepening}}{\underbrace{i\frac{\gamma }{{\omega }_0}\frac{\partial }{\partial t}\left({\left|A\right|}^{2}A\right)}}-\underset{\mathrm{Raman}}{\underbrace{T_{\mathrm{R}}\gamma A\frac{\partial }{\partial t}{\left|A\right|}^2}}+\underset{\mathrm{loss}}{\underbrace{i\frac{\mathsf{\alpha }}{2}A}}=0.\hfill \end{array}$$

(16)

As the generalized NLSE is not integrable, it usually must be treated numerically. The impact of the new terms is as follows:

The dispersion curve can be written as a Taylor expansion around the carrier frequency ${\omega }_0$; the NLSE is written so that the ${\beta }_2$ term appears. All higher-order terms (shown here up to the fourth order) provide corrections that become relevant when either the pulse spectrum gets quite wide or when the carrier frequency is close to a zero-dispersion point (wavelength where ${\beta }_2$ passes through zero) of the fiber.

The ‘self-steepening’ term [1,23] results from the intensity dependence of the group velocity and leads to an asymmetry in the pulse shape. It can produce shifts in spectral and temporal positions even in the absence of the Raman term.

The Raman effect leads to a redistribution of spectral power to the advantage of the low-frequency slope of the pulse spectrum. A first study of the Raman gain in fibers was done by R. H. Stolen et al. [24]. The gain curve peaks at a frequency offset of ≈13 $\mathrm{THz}$, but is nonzero all the way down to zero offset. As a result, the spectral ‘center of mass’ of a pulse is shifted to lower frequencies. This ‘redshift’ scales with the inverse fourth power of the pulse width [25]. The term expressing this in the NLSE is shown here in a simplified version, which contains the Raman response time $T_{\mathrm{R}}$. This model is valid for pulses being longer than ≈100 $\mathrm{fs}$; more accurate models are described in [26,27,28].

The loss term contains Beer’s loss coefficient $\mathsf{\alpha }$; the impact of the loss can be characterized by the characteristic loss length $L_{\mathsf{\alpha }}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\alpha }$}\right.$. As long as the loss is weak, i.e., $L_{\mathsf{\alpha }}\gg L_{\mathrm{D}},L_{\mathrm{NL}}$, the equilibrium of dispersive and nonlinear effects that characterizes the soliton is only mildly perturbed. In that situation, the soliton can rearrange its shape: With the peak power drooping and the width increasing slightly, Equation (4) can still be approximately fulfilled even when the energy is reduced [29,30].

3.2. Dispersion-Managed Fibers

As a further complication in a realistic assessment of pulse propagation in fibers, one must acknowledge that in the telecommunications field, a special type of fiber has been preferred for 20 years now in which ${\beta }_2$ is not a constant. Originally, an attempt was made to cancel dispersive effects by concatenating segments of fiber with alternatingly positive and negative group velocity dispersion. As it turned out, a perfect cancellation (path average ${\overline{\beta }}_2=0$)—possible at a single wavelength only anyway—did not work as well as if a small negative value were left. This is on account of the fiber’s nonlinearity.

Figure 6 illustrates the structure of a dispersion-managed fiber (DM fiber). The dispersion map is described by the spatial period $L_{\mathrm{map}}=L^{+}+L^{-}$ with which the alternating pattern is repeated [31]. Strictly speaking, the value of the nonlinearity coefficient $\gamma$ will also alternate from one fiber segment to the other.

It is not immediately clear whether in a dispersion-managed (DM) fiber something like solitons can exist, as they are solutions of the NLSE for anomalous dispersion only. However, as it turned out, a type of pulse exists that is stabilized by nonlinearity [32,33,34,35,36,37]. These pulses were then called dispersion-managed (DM) solitons. Their shape is neither of the usual unchirped sech type, nor constant. Rather, it breathes over a full dispersion period so that after a full period, the shape in both the amplitude and phase profile is restored. An example of a DM soliton is shown in Figure 7.

Mathematically, dispersion management is captured by making fiber parameters z-dependent:

$$i\frac{\partial }{\partial z}A-\frac{{\beta }_2\left(z\right)}{2}\frac{{\partial }^2}{\partial t^2}A+\gamma \left(z\right){\left|A\right|}^{2}A=0.$$

(17)

This modifies the NLSE into a non-integrable form so that it does not support solitons in the usual sense. Depending on the dispersion modulation strength (essentially, ${\beta }_2^{+}$–${\beta }_2^{-}$), the DM soliton shape may be nearly sech-shaped or closer to a Gaussian [38]. In the strong modulation limit DM, solitons exhibit oscillating tails [38,39]. There is no closed analytical expression for it, but good approximations to its shape can be found by numerical procedures like that given in [40]. During propagation, a DM soliton suffers from continuous radiative loss, which may be weak, but nevertheless must lead to its eventual decay. Fortunately, for applications, the decay distance typically exceeds any practically relevant fiber length.

It was eventually understood that dispersion-management provides advantages over a constant-dispersion fiber. It suppresses four-wave mixing, a process occurring in a transmission system with many wavelength channels filled simultaneously with independent data streams, where it leads to cross-talk and, thus, errors. An extensive review of dispersion-management is given in [31].

3.3. Cases Involving Several Solitons

We mention a few more cases in which solitons or their compounds play a role and which have received attention in the research literature.

3.3.1. Soliton Molecules

In DM fibers, stable compounds of solitons exist, which have been called soliton molecules [41]. The alternating dispersion creates a rapid to-and-fro of the chirp of DM solitons and, thus, a varying mutual force. If two DM solitons are in a particular mutual separation, the net force is zero, and an equilibrium is obtained [42,43]. For larger separations, there is attraction, for smaller, repulsion: it is a stable equilibrium. Soliton molecules also exist for more than two DM solitons and have been discussed as information carriers in data transmission [44,45]. As it turns out, there can be several equilibrium positions [46].

3.3.2. Soliton Gas and Crystal

Another situation in which soliton propagation is subject to periodic perturbation was described in a series of experiments about a synchronously-driven fiber resonator [47] (using non-DM fiber). Free propagation in the fiber alternated with interference at an input coupler. As a result, relatively long input pulses organized themselves into a pattern that contained a multitude of short pulses; each of those fulfilled the soliton condition Equation (4). The pattern could be either chaotic (soliton gas) [48] or periodic (called soliton crystal) [49], depending on the drive power.

Similar phenomena have also been reported from fiber lasers [50,51,52], which, however, are not described by the NLSE, as that does not consider nonlinear gain and loss. An equation including these effects is the cubic-quintic complex Ginzburg–Landau equation. One finds that there are soliton-like solutions of this equation that are called dissipative solitons [53,54,55].

3.3.3. Supercontinuum Generation

For some metrological purposes like coherence tomography, wide bandwidth light with good focusability is required. Thermal light sources cannot deliver adequate power due to both their insufficient spatial coherence and the limitations set by Planck’s radiation law. One can, however, generate broadband light in a fiber by exploiting the nonlinearity to create new frequency components. Experiments involve, in most cases, a powerful mode-locked laser and a piece of highly nonlinear fiber. Either a modulation instability according to an Akhmediev breather scenario evolves, but the structure then breaks apart, or the input pulse represents a high-order ($N\gg 1$) soliton, which undergoes fission [56]. In either event, the initial light signal breaks up into an arrangement of several pulses, many of which have soliton characteristics, and radiation. This arrangement has similarities to the soliton gas described above. Through the action of the Raman effect, four-wave mixing and cross-phase modulation combined with complex phase matching properties in the fiber’s dispersion curve (see, e.g., [57]), the light quickly develops a spectrum that may easily be one octave wide and often much more. Such a structure is referred to as an optical supercontinuum and has found much research interest; commercial supercontinuum sources are being offered. A good review is [58].

A corollary to this is the phenomenon of optical rogue waves. Rogue waves of the ocean have been identified as a real phenomenon, not a myth. Unusually, large single waves can occasionally appear without warning, possibly wreaking havoc on ships, and disappear again just as suddenly. It was first discussed in [59] that in the context of optical supercontinuum generation, single spikes of extreme power can occur every once in a while. There is an ongoing discussion about the mechanisms [60].

4. Methods to Verify Soliton Content

The inverse scattering technique allows one to find solutions to integrable nonlinear wave equations. However, in practice, one quite often is faced with the reversed task of determining the soliton content of a given signal. The task is relatively straightforward when the signal consists of isolated bell-shaped pulses. One can check whether the power profile approximates a ${\mathrm{sech}}^2\left(t\right)$ shape and whether the chirp is close to zero. Then, duration and peak power can be checked with Equation (4). This way, an $N=1$ soliton can be directly identified, and for pulses with $N\ne 1$, the local value of N can be found. However, the applicability of this direct method is limited to simple cases. It is of no use for waveforms, as they may be found in more complex situations like soliton gas or optical supercontinuum; see Section 4.1.3.

We will therefore briefly describe more advanced techniques. The direct scattering method is often considered as the benchmark for the determination of soliton parameters. However, it has its limitations, chief among which is that it relies on the integrability of the system. An alternative is soliton radiation beat analysis, which does not have that same restriction, but the price to be paid for that is that (i) more input data are required and that (ii) in complex cases, the obtained pattern may get too complicated to interpret.

4.1. The Inverse and Direct Scattering Transform

4.1.1. The Method

The inverse scattering transform (IST) is a technique that can be used to solve certain nonlinear partial differential equations (PDE) like the Korteweg–de Vries equation, the sine-Gordon equation or the nonlinear Schrödinger equation. The basic concepts of this method were presented in 1967 in a groundbreaking paper by Gardner and coworkers dealing with the Korteweg–de Vries equation [5]. Shortly after, Lax showed how to apply the method also to other nonlinear PDEs when certain conditions are met [61]. The IST was successfully adapted to the NLSE by Zakharov and Shabat in 1972 [4].

The method of IST is based on transformations between the time-space domain and a nonlinear spectral domain, in which the nonlinear evolution of some given input field reduces to a linear problem. In the spectral domain, a given field is represented by its scattering spectrum, which is calculated by the so-called direct scattering transform (DST). The inverse operation is called the inverse scattering transform. As these transformations can be regarded as extensions to the well-known Fourier transform, they are also known as nonlinear Fourier transform and inverse nonlinear Fourier transform, respectively. Using the nonlinear Fourier transform, one can get insight into the linear and nonlinear components of a given field by analyzing its scattering spectrum. Generally, it consists of a continuous part, which represents the small amplitude radiation, and a discrete part, which can be attributed to the solitons contained.

In recent years, several groups took up the idea of using multi-soliton pulses for nonlinear optical data transmission, which originally was suggested by Hasegawa and Nyu [62]. In the course of this renewed attention, several nonlinear transmission schemes and techniques for the nonlinear Fourier transform and its inverse transform have been (and still are) developed to overcome the linear capacity limits of fiber-optical transmission lines. An overview about the latest developments can be found in [63].

DST is a great tool to analyze given pulses, but its application is somewhat involved. The DST has been applied analytically to only a few simple pulse shapes like sech-shaped [6] or rectangular pulses [64]. Even for such a standard shape as a Gaussian,

$$A\left(t\right)=A_0\sqrt{P_0}\exp \left[-\frac{1}{2}{\left(\frac{t}{T_0}\right)}^2\right],$$

(18)

no analytical result is known to the authors’ best knowledge. Instead, approximate solutions can be found numerically. As an example, when setting $A_0=3$ and choosing $P_0$ and $T_0$ according to Equation (4), DST reveals that the pulse consists of two solitons and a small amplitude radiation part of nontrivial spectral shape (see Figure 8).

Various schemes for numerical calculations of the DST are available [63,65,66,67,68]. They all solve the Zakharov–Shabat eigenvalue problem of a discretized input field with M sample points and determine the nonlinear spectrum. One can distinguish between matrix methods that calculate the eigenvalues directly from a $M\times M$ matrix and search methods in which the eigenvalues are found in an iterative manner. For a recent discussion, see, e.g., [63,66,67,68].

In the following, we will sketch briefly how to calculate the scattering spectrum using a search method. Using the formulation of Ablowitz et al. [69], one finds that solving the dimensionless NLSE:

$$\frac{\partial q}{\partial z}=\frac{i}{2}\frac{{\partial }^{2}q}{\partial t^2}+i{\left|q\right|}^{2}q$$

(19)

is equivalent to the integration of the system of equations:

$$\begin{array}{cc}\hfill \frac{\partial \mathrm{\Phi }}{\partial t}+T\mathrm{\Phi }& =0,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \frac{\partial \mathrm{\Phi }}{\partial z}+Z\mathrm{\Phi }& =0.\hfill \end{array}$$

(21)

T and Z are operators (matrices) that determine the temporal and spatial evolution of some auxiliary function:

$$\mathrm{\Phi }=\left(\begin{array}{c}{\varphi }_1\\ {\varphi }_2\end{array}\right)$$

(22)

where ${\varphi }_1$ and ${\varphi }_2$ are called Jost functions. Explicitly, the operator T can be written as:

$$T=\left(\begin{array}{cc}-i\lambda & q\\ -q^{*}& i\lambda \end{array}\right).$$

(23)

Calculation of the DST can be done by integrating Equation (20) with respect to time with different test values of the eigenvalue candidate $\lambda =\zeta$. The initial condition can be chosen as:

$$\mathrm{\Phi }(z_0,t,\zeta )=\left(\begin{array}{c}1\\ 0\end{array}\right)\exp (-i\zeta t)\mathrm{with}t\to -\infty ,$$

(24)

when q is localized in the sense that for $t\to \pm \infty$, $q\to 0$ at an exponential rate [65,66]. As a result, one obtains the scattering coefficients $a\left(\zeta \right)$ and $b\left(\zeta \right)$ from $\mathrm{\Phi }(z_0,t,\zeta )$:

$$\begin{array}{cc}\hfill a(z_0,\zeta )& ={\varphi }_1(z_0,t,\zeta )\exp \left(i\zeta t\right)\mathrm{with}t\to \infty ,\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill b(z_0,\zeta )& ={\varphi }_2(z_0,t,\zeta )\exp (-i\zeta t)\mathrm{with}t\to \infty .\hfill \end{array}$$

(26)

Additionally, the derivatives of a and b with respect to $\zeta$ can also be calculated by extending the integration scheme to include ${\mathrm{\Phi }}^{\prime }=\partial \mathrm{\Phi }/\partial \zeta$.

In the scattering spectrum for anomalous dispersion, solitons are present in the form of discrete eigenvalues ${\zeta }_k$ with $a\left({\zeta }_k\right)=0$ and $\mathrm{Im}\left({\zeta }_k\right)>0$. $\mathrm{Re}\left({\zeta }_k\right)$ represents the frequency, and $2\mathrm{Im}\left({\zeta }_k\right)$ corresponds to the amplitude of the soliton; both values are preserved during propagation. The continuous part of the spectrum (radiation) can be calculated from the values of $a\left(\zeta \right)$, with $\mathrm{Re}\left(\zeta \right)=0$. The spectral power density ${\left|\mathcal{F}\right|}^2$ is then obtained as:

$${\left|\mathcal{F}\right|}^2=-\frac{1}{\pi }\ln \left(\right|a(\zeta ,z)\left|\right).$$

(27)

4.1.2. Numerical Restrictions

We have applied numerical DST to a case for which exact analytical results are known, i.e., the fundamental soliton Equation (3), to assess inaccuracies of the numerical procedure. They arise from various sources; we identify three main error sources.

Truncation error: At the edges of the computational time window of finite width $T_{\mathrm{W}}$, the pulse wings are truncated. A fraction $E_{\mathrm{trunc}}$ of the energy is not represented in the numerical field:

$${\displaystyle E_{\mathrm{trunc}}=2\underset{T_{\mathrm{W}}/2}{\overset{\infty }{\int }}\widehat{P}{\mathrm{sech}}^2\left(\frac{T}{T_0}\right)dT=E_{\mathrm{tot}}\left(1-\tanh \left(\frac{T_{\mathrm{W}}/2}{T_0}\right)\right).}$$

(28)

Sampling error: Sampling error arises from coarse sampling when the temporal discretization step is too wide for the pulse duration.

Floating point error: Computers store numbers in some format with a limited number of digits; ‘double precision’ usually offers 15 digits plus the exponent. After a large number of floating point operations, numerical errors accumulate.

To disentangle these contributions, we varied the number of sampling points and the width of the time window for a given pulse width; the results are shown in Figure 9. To the right, the computational window is sufficiently wide ($T_{\mathrm{W}}\gg T_0$). Then, the three curves shown are for discretization step numbers as indicated by the labels. For $2^{12}$ points, sampling is increasingly coarse when these points are spread over a wider window, and the error rises. Larger numbers of points improve the situation quadratically because the error scales with the square of the sampling point separation. To the left, truncation error dominates; the prediction of Equation (28), plotted as a solid curve, fits well. Floating point error becomes relevant only where these two contributions are small and where the number of calculation steps is large; in our example, it appears as ‘noise’ in the central section of the curve for $2^{22}$ points. We convinced ourselves that the transfer-matrix (TM) method presented in [65] yields the same results.

4.1.3. Applications and Limitations of DST

The DST method as described so far is quite generally applicable to arbitrary pulse profiles provided two conditions are met: One is integrability, which is not fulfilled for many structures in the generalized NLSE in Section 3.1 and for DM systems in Section 3.2. The other is that for time to infinity, the structure must decay to zero sufficiently fast [65,66]. This is not the case for the structures described in Section 2.2.

If integrability is violated only weakly either by loss or through the Raman effect, the concept of solitons remains valid in an approximate sense as pulses that maintain an (near-) equilibrium between dispersive and nonlinear influences. It might be more accurate to speak then of solitary pulses. Strictly speaking, though, all preserved quantities are no more. Perturbation theory has shown [29,30] that solitons rearrange their shape in terms of width and peak power; in the process, they acquire a new energy value, as well. In the sense that this may still be interpreted as solitons, some researchers have adopted the following strategy: In a mildly non-integrable system with loss [70] or Raman effect [71,72,73], the pulse structure is propagated numerically up to some position $z_i$, where DST is then applied to find what one may call ‘local eigenvalues’. Then, one lets $z_i$ slide and repeats DST, in order to find the position dependence, or spatial evolution, of the eigenvalues. An example of such a case is shown in Figure 10 where the Raman effect, third-order dispersion and loss are treated separately. The most conspicuous effects are the energy reduction by the loss term and the frequency downshift by the Raman effect.

When it comes to complex situations like in the context of optical supercontinuum generation where many interactions and collisions of several pulses take place, energy can be transferred between pulses, and more changes to local eigenvalues will occur. In any event, one has to interpret the results produced in this way with caution.

If one considers infinitely-extended structures, one first needs to acknowledge that they have infinite energy. One possible approach is to start with a finite piece (‘window’) of the infinite structure and determine eigenvalues, then increment the window width and follow the eigenvalues; finally, let the window width tend to infinity [74]. In the limit, one obtains an infinite number of eigenvalues, but with a finite density (number of eigenvalues per unit time). These eigenvalues cover a continuum of energies from zero up to some maximum; if one considers the simplest case of an unmodulated continuous wave, this maximum depends only on its power $P_0$ and is given by $E_{\max }=2\sqrt{\left|{\beta }_2\right|P_0/\gamma }$. With respect to the interpretation of these eigenvalues as signatures of solitons, the following was shown in [74]: Consider a finite window, which contains m eigenvalues, and allow it to propagate. On account of integrability, it contains a constant number of m solitons. If one then applies a perturbation, e.g., by Raman effect, precisely m separate pulses emerge; once they are sufficiently separated from each other after some propagation, it can be checked that each fulfills the soliton criterion Equation (4).

The other approach to infinitely extended fields is to consider periodic problems, like the Akhmediev breather. Adaptations of IST to periodic boundary conditions (periodic IST, also called finite gap theory) have been given in [75,76,77] where the obtained eigenvalues are not necessarily interpreted as solitons, but rather as the representation of linear and nonlinear modes [78,79,80]. Like the nonperiodic IST, this periodic case gained renewed interest in recent years. Different theoretical and numerical schemes have been developed with a view toward the application of nonlinear data transmission [63,66,67,81,82,83]. A variant of the periodic DST based on the artificial periodization of field components was used to analyze rogue events and numerically to find the eigenvalue spectra of the breather solutions of Section 2.2 [84]. Another approach to characterize rogue events via DST was taken in [85,86].

We have here concentrated on pulse-like initial conditions. Recently, the scope was extended to integrable turbulence [87] with a view towards understanding rogue waves. Initial conditions for the NLSE were generated by superimposing random perturbations in both amplitude and phase on a continuous wave; the resulting field was characterized by its eigenvalue spectrum.

4.2. Soliton Radiation Beat Analysis

Soliton radiation beat analysis (SRBA) is a numerical method to determine the soliton content of a pulse in an optical fiber even if the system is not integrable. In contrast to DST, which requires the pulse profile at a single position in the fiber, SRBA requires knowledge of the evolution over a certain length section of the fiber. This is easily understood: In integrable systems, soliton parameters are invariant so that if they are known in one place, they are known everywhere. SRBA can follow the evolution of soliton parameters in a non-integrable setting; the price to pay is that representative information about it must be introduced into the calculation.

As early as in 1974, it was noted [6] that the envelope of a pulse not exactly conforming to a soliton undergoes oscillations. Such a pulse contains a soliton and some radiation; the oscillation arises as a beating between both. The radiation evolves linearly; the soliton contains a power-dependent phase term as shown in Equation (3). The frequency of the beat note therefore allows a direct conclusion about the soliton energy.

The SRBA method was first introduced in [88]. The propagation of a sech pulse with the initial form $A(z=0,t)=N\sqrt{P_0}\mathrm{sech}(t/T_0)$ was simulated numerically. This was repeated while incrementing the soliton parameter N. As the oscillation of the peak power is damped due to dispersive spreading of the radiative part, it was found advantageous to evaluate the peak spectral power, which undergoes undamped oscillation. A windowed Fourier transform was applied to the spatial evolution of the power at the center frequency $|\tilde{A}{(\mathrm{\Omega }=0,z)|}^2$. This produces a spectrogram of spatial frequencies called the soliton radiation beat pattern. An example is shown in Figure 11 (left).

As N is increased, a soliton radiation beat first appears at its expected threshold near $N=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$. The figure shows that this beat note also has overtones. Its amplitude (grey scale) vanishes at integer N (apparent ‘white stripes’), which is of course because there, the radiation amplitude vanishes; the beat amplitude is proportional to the product of both amplitudes.

If there is more than one soliton involved, the beat pattern will contain the frequencies of beats between all pairs of solitons, as well as that of each soliton beating with radiation. In Figure 11 (left), traces pertaining to the second soliton begin at $N=\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ as expected and have amplitude nulls at $N=2$. Difference frequencies are visible as downward-bending curves. In the right panel of Figure 11, analytical results for beat frequencies are plotted for comparison. Shown are all calculated frequencies, up to four overtones, and all combination frequencies of this multitude. It is obvious that all SRBA traces of the left panel can be identified in this way. SRBA provides the additional information that there are notches in the amplitudes of all traces involving radiation at integer N values.

This shows that SRBA reproduces all known features of the soliton energies as known analytically from IST. In [89], this method was shown to work for chirped Gaussian pulses in the NLSE. Going beyond the NLSE to non-integrable situations is also possible as SRBA makes no assumption about integrability. A lossy fiber was treated in [90]. It was shown that a mild, ‘adiabatic’ loss, if it persists for a long fiber length, eventually turns diabatic because the relevant length scales $L_{\mathrm{D}}$ and $L_{\mathrm{NL}}$ also change in the process. Once the loss is diabatic, the soliton decays. In that situation, however, the applicability of DST even in the approximate sense as described above in Section 4.1.3 is ultimately lost, and interpretations based on that approach [70] are doubtful and at variance to the SRBA results.

A further important example of the application of SRBA is the determination of the dispersion-managed soliton content. This case is just not accessible for DST, not even in the sense of ‘local’ analysis as in Section 4.1.3. Other than SRBA, no other numerical method to analyze DM structures is known. In [88], a suitably-scaled Gaussian pulse was used as a good approximation to an ideal DM soliton, and an SRBA pattern was obtained. Both the modification of the soliton energy and the characteristic minima of soliton radiation beat traces were found.

In subsequent research, SRBA was also applied to analyze Akhmediev breathers and a sinusoidally modulated cw [91], as well as soliton crystals [92]. One can assume that in principle, the SRBA method can be applied to all non-integrable systems where a nonlinear phase shift appears. Potential applications are soliton molecules (see Section 3.3.1) or dissipative solitons (see Section 3.3.2).

SRBA as described so far relies on the evaluation of the central spectral power ${\tilde{A}}^2(\mathrm{\Omega }=0,z)$. In a generalization of the technique, the consideration is extended to the full spectrum, ${\tilde{A}}^2(\mathrm{\Omega },z)$ [93]. In this variant, the technique also yields information about soliton velocities. Figure 12 shows an example for the case of an $N=2$ soliton.

The pattern was generated using a value of $N=2.01$, so that there is always a little bit of radiation copropagating with the second-order soliton, and a beat note is generated. Parabola-like curves are due to radiation; horizontal flat lines touching the parabolas are due to solitons. Once the traces are identified, the soliton energies can be obtained from the spatial frequency of the nonlinear phase rotation. To illustrate the interpretation of such data, we point out that the two solitons appear with spatial frequencies in the ratio of 9:1. According to Equation (3), this means that the peak powers are $P_0\left(S_1\right)=9P_0\left(S_2\right)$. Equation (4) then dictates that $T_0\left(S_1\right)=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.T_0\left(S_2\right)$. The ratio of energies obtained from Equation (5) is then 3:1, precisely as expected from Figure 3. Any frequency shifts would appear as horizontal displacements; in this particular example, there are none. There is no scan of N here; this extended procedure has the advantage that only a single propagation simulation is necessary.

A particular difficulty with DM solitons is that their shape undergoes continuous changes. At various positions within $L_{\mathrm{map}}$, the shape varies considerably, but the pulse shape repeats after $L_{\mathrm{map}}$ (it is ‘stroboscopically’ stable). Even then, when the power is modified, the shape varies from ‘nearly sech’ to ‘nearly Gaussian’, and depending on the dispersion modulation depth, oscillating tails may appear. Therefore, extra measures need to be taken when the energy is scanned for an SRBA pattern calculation. For the data in the left panel of Figure 13, at 0.5 pJ increments of energy, the best approximation to the DM soliton shape was found using Nijhof’s method [40]. If the pulse shape were exact, there would be no beat note with radiation. This is why in 0.5 pJ increments, there seem to be white lines in the figure. Above and below each such step, the energy was scanned by simple scaling of the pulse shape to $\pm 0.25\mathrm{pJ}$. The slight deviation from the best pulse shape suffices to produce beat notes, and a full survey of their positions is obtained. For better clarity, the data from this calculation are transferred to the right panel of Figure 13, and different line thicknesses are standardized into three classes in the process. The data show that a DM soliton consists of several constituent solitons; the pertaining traces are labeled $S_1$–$S_3$.

A comparison of DST and SRBA reveals the advantages and weaknesses of both. SRBA requires more input data. These can only be available from numerical simulations; it is hard to imagine how suitable data could be obtained from an experiment. On the other hand, SRBA can treat cases that are inaccessible to DST. Its limit is reached when for a growing number of solitons involved, the number of traces becomes too large to disentangle them all. DST can be applied to many-soliton cases [94], but only in strictly integrable systems, or with the limitations in ‘nearly integrable’ systems as pointed out in Section 4.1.3.

5. Conclusions

Solitons in optical fibers are pulses of light that maintain a balance between linear and nonlinear distortions. Real-world systems are rarely, if ever, of the integrable kind. It is therefore appropriate to apply the word ‘soliton’ also to pulses for which the balance is maintained at least approximately; to do so is now common usage. However, one occasionally sees research papers in which the word is used in a very loose sense, with no effort being made to actually check whether the balance is fulfilled.

In this review, we outline the concept of fiber-optic solitons, with a particular view toward methods of verification of the solitonic nature of a given pulse, or group of pulses, or a similar structure. Soliton content has been discussed in the context of structures like soliton molecules, soliton gas and crystal, Akhmediev breather, optical supercontinuum and rogue waves. The simplest method is to check pulse parameters with Equation (4), which directly expresses an equilibrium condition. Unfortunately, this straightforward test is only applicable in very simple situations when there are clean, well-isolated pulses, and it does not fully quantify the soliton content. Grown from inverse scattering theory, the first method to analytically find the soliton solution of the NLSE is direct scattering transform. It basically tests whether the structure under test contains pulses that hang together as entities even when they collide, the property that gave solitons their ‘particle-like’ name [95]. An alternative method is soliton radiation beat analysis, which basically tracks the nonlinear phase rotation, which is the signature of solitons. We discuss some issues of practicality and accuracy and compare the ranges of applicability of both DST and SRBA.

There is no single method that covers all cases; indeed, in some relevant cases, no available method is fully satisfactory. In the complex evolution of supercontinuum, for example, the method is often adapted ad hoc to the situation, for want of alternatives. In some interesting cases, soliton content analysis has not been pursued yet, to our knowledge. Dissipative solitons are subject to yet another condition beyond Equation (4), but should be accessible for SRBA. This article is meant to provide the reader with a survey to what is currently known, to foster future research.