format: text

text:
# 1 Introduction

## 1 Introduction

The HERA collider, where 27.5 GeV positrons collide with 820/920 GeV protons, is traditionally considered as a natural extension of Rutherford’s experiment and the process of deep inelastic scattering ($`ep`$ DIS) is interpreted as a reaction in which a virtual photon, radiated by the incoming positron, probes the structure of the proton. However, there are regions of phase space for $`ep`$ collisions where it is other way round: where partons in the proton probe the photon structure.

The structure of the talk will be the following: it will start with posing the problem i.e. a description of our understanding of the structure of the photon in terms of Quantum Field Theory (QFT) and the motivations for its investigation. Next, we briefly review measurements of the photon structure in $`e^+e^{}`$ experiments and its general properties. In the main part, we review HERA measurements sensitive to the structure of the photon. We end with a summary and conclusions.

Various aspects of the subject discussed here have been covered recently in many excellent review articles e.g. , , , .

## 2 Structure or interaction dynamics?

The photon is the gauge boson of QED and, as far as we know, elementary. In any QFT the existence of the interaction means also that the quanta themselves have structure. Coupling between field quanta will make possible fluctuations, or splitting of any quantum over limited time, into two (or more in higher order) quanta. For example in QED a photon can fluctuate into an electron-positron ($`e^+e^{}`$) pair, or, as first noted by Ioffe , into a $`q\overline{q}`$ pair. If the fluctuation time defined in the probing objects rest frame as $`t_f2E_\gamma /m_{q\overline{q}}^2`$ is much larger than the interaction time $`t_{int}`$, the photon builds up structure in the interaction. Here, $`E_\gamma `$ is the energy of the fluctuating photon and $`m_{q\overline{q}}`$ is the mass into which it fluctuates. The fluctuation time of a photon with virtuality $`Q^2`$ is given by $`t_f2E_\gamma /(m_{q\overline{q}}^2+Q^2)`$, and thus at very high $`Q^2`$ one does not expect the condition $`t_ft_{int}`$, unless the probe has even higher virtuality.

It is clear that what we described here as the photon structure is a part of the interaction dynamics. In practice in any high energy interaction, the distinction between structure and interaction is determined by the factorisation scale. This is a momentum transfer squared scale below which any parton activity is considered to be part of the parent structure and above which any parton activity is considered to be a part of the interaction dynamics. It is exactly in this sense that in $`ep`$ DIS we measure the proton structure function $`F_2(x,Q^2)`$ at a factorisation (or resolution) scale $`Q^2`$: $`Q^2`$ is almost always the largest momentum transfer squared scale in DIS, so all parton activity at lower transverse momenta should be associated with the target i.e. proton. The structure function is interpreted as a density of quarks and antiquarks which carry a fraction $`x`$ of the proton momentum. Any structure function is thus dependent on the choice of the factorisation (resolution) scale. It was shown that structure functions for hadrons can be used universally in any high energy interaction at a given scale. In the same way it is also expected that photon structure $`F_2^\gamma (x_\gamma ,Q^2)`$ will be similarly universal and factorisable.

When interactions are switched off the photon remains structureless so we say that photon structure is driven by fluctuations. It is in many respects different from that of hadrons which have structure even after interactions are switched off (valence driven structure). In addition, the photon structure can be investigated as a function of its virtuality i.e. size, which is not possible for hadrons. We have thus several strong motivations for the photon structure studies:

* Measurements of the photon structure function provide interesting tests of QCD. Some pose even a challenge to theory e.g. there is no full theory of the parton density suppression due to increasing virtuality (only an asymptotic prediction exists)
* The parametrisation of the universal and factorisable photon structure function based on LEP and HERA measurements is a necessary component of the description of many high energy processes
* At future linear colliders it will be very important to understand the large number of events from photon interactions. Some photon interactions can be important probes to look for new physics, for example the production of Higgs bosons at a photon linear collider ($`\gamma \gamma H`$) .

## 3 Photon structure from $`e^+e^{}`$ and its general properties

The hadronic structure of the photon can be probed either by virtual photons in $`e^+e^{}`$ collisions or by virtual partons in $`ep`$ collisions at HERA. This article is devoted mainly to the recent measurements at the HERA collider, but for the sake of completeness some $`e^+e^{}`$ results are discussed in this section.

The principle of the measurement of the hadronic structure function of the photon $`F_2^\gamma (x_\gamma ,Q^2)`$ is depicted in Figure 1 (left). A highly virtual $`\gamma ^{}`$ with large $`Q^2=q^2`$ probes a quasi-real $`\gamma `$ with virtuality $`P^20`$. In this measurement the energy of the target is not known so $`x_\gamma `$ (the fraction of the photon momentum carried by the interacting parton) has to be reconstructed from the hadronic final state. This leads to systematic uncertainties and a limited reach for small values of $`x_\gamma `$. To underline the analogies and differences we depict the principle of the measurement of the proton structure function on the right of Figure 1.

Many measurements of the hadronic structure function of the photon have been performed at several $`e^+e^{}`$ colliders. The range of resolution scale covered by various experiments is $`0.24Q^2700\mathrm{GeV}^2`$.

In Figure 2 we present a compilation  (recently updated with OPAL data ) of photon structure function measurements in various $`e^+e^{}`$ experiments. $`F_2^\gamma `$ is plotted as a function of $`Q^2`$ for several ranges of $`x_\gamma `$. One can see positive scaling violation for all $`x_\gamma `$. This is the most striking feature of the photon structure in comparison to the structure of the proton, for which positive scaling violation is observed only at low $`x`$ values. This different behavior can be understood in terms of perturbative QCD as coming from the splitting $`\gamma q\overline{q}`$, having no analogue in the proton case. In addition, and again contrary to the proton case, it also causes the photon structure function to be large for high $`x`$ values.

A simple LO parametrisation of the photon structure function

$$F_2^\gamma =3\underset{q}{}e_q^4\frac{\alpha }{\pi }x[x^2+(1x)^2]\mathrm{ln}\frac{Q^2}{\mathrm{\Lambda }^2}+[{}_{}{}^{}\mathrm{hadron}_{}^{}\mathrm{VDM}]$$
(1)
can be written down in QCD. The expression is asymptotic i.e. all terms not multiplied by large $`\mathrm{ln}Q^2`$ are dropped. The second term in the equation – “hadron” part, corresponds to $`\gamma q\overline{q}`$ fluctuations with virtuality below a cutoff value $`\mathrm{\Lambda }^2`$. This noncalculable piece is traditionally represented by photon fluctuations into vector mesons, as proposed 30 years ago in terms of the Vector Dominance Model . Witten  and others showed that the above mentioned features of the LO formula (1) remain valid in calculations which include higher order perturbative terms. Thus a rising dependence of $`F_2^\gamma `$ with increasing probe scale $`Q^2`$ throughout most of the range of $`x`$ is thereby seen to be a characteristic feature of the “fluctuation driven” (rather than “valence driven”) structure functions. Measurements compiled in Figure 2 provide therefore an interesting and important test of pQCD.

Since the measurements of $`F_2^\gamma `$ span a wide range of $`Q^2`$, it is possible to apply the DGLAP pQCD formalism  as is done for the proton structure function. As a result parton density functions (pdfs) for the photon are available.

There exist several pdfs for real and also for virtual photons in leading and next-to-leading order, which are based on the full evolution equations. They are constructed in manner very similar to that for the proton pdfs. They differ in the assumptions made about the starting scale, the input distributions assumed at this scale and also in the amount of data used in fitting their parameters. Those most quoted in the literature are LAC , GRV , AGF  and SAS . These parametrisations were developed before the publication of the new LEP data. LEP has measured the photon structure function $`F_2^\gamma `$ in the range $`0.002x_\gamma 1.0`$ and $`1.86<Q^2<700\mathrm{GeV}^2`$. These data represent an important step in the reduction of statistical and systematic uncertainties and cover a larger kinematic region both in $`x_\gamma `$ and $`Q^2`$. From Figure 2 (right) it is clear that the LEP data provide new information on photon structure and the existing pdf parametrisations will have to be revisited.

## 4 Photon structure from HERA

At HERA the photon structure is investigated mainly by measurements of high $`E_T`$ jets and charged particles -. The partonic structure of quasi-real or virtual target photons is probed by a parton from the proton producing a pair of partons of large transverse momentum squared $`E_T^2`$ much larger than the photon virtuality $`Q^2`$. Due to the large cross section at HERA the photon can be probed at even larger factorisation scales than at LEP.

The photoproduction of jets through a hard parton differs from the jet production in hadron-hadron collisions in one respect: the photon interacts either by means of its direct coupling to high $`E_T`$ partons or by means of partons in its structure, which are resolved in the interaction.

This is depicted in Figure 3. The terminology “direct” and “resolved” processes is now established, though in reality the direct photon contains fluctuation driven structure not resolved at a given scale.

The fraction of the photon momentum carried by an interacting parton can be reconstructed at hadron level from jet kinematics

$$x_\gamma ^{jet}=\underset{jets}{}\frac{E_{Tjet}e^{\eta _{jet}}}{E_\gamma }$$
(2)
where $`E_{Tjet}`$ and $`\eta _{jet}`$ are the transverse jet energy and pseudorapidity of the jets. The sum extends over high $`E_T`$ jets (i.e. does not include the photon and proton remnant jets). The “true” (i.e. parton level) fraction $`x_\gamma `$ can be unfolded from the data using simulated $`x_\gamma x_\gamma ^{jet}`$ correlations. In order to reach low values of $`x_\gamma `$ the measurement should extend to low values of the jet transverse momentum and/or high values of their pseudorapidity (i.e. close to the proton remnant). Both requirements are experimentally challenging. In particular, the low cut on $`E_T`$ results in large corrections for hadronisation and soft underlying events i.e. secondary interactions of the partons belonging to the photon and proton remnants.

In the LO QCD approximation to di-jet photoproduction, the $`ep`$ cross section is written

$$\frac{\mathrm{d}^5\sigma }{\mathrm{d}y\mathrm{d}x_\gamma \mathrm{d}x_p\mathrm{d}\mathrm{cos}\theta ^{}\mathrm{d}Q^2}=\frac{1}{32\pi s_{ep}}\frac{f_{\gamma /\mathrm{e}}(y,Q^2)}{y}\underset{ij}{}\frac{f_{i/\gamma }(x_\gamma ,E_T^2,Q^2)}{x_\gamma }\frac{f_{j/p}(x_\mathrm{p},E_T^2)}{x_p}|M_{ij}(\mathrm{cos}\theta ^{})|^2$$
(3)
where $`f_{i/\gamma }`$ and $`f_{i/p}`$ are the pdfs for each parton species $`i`$ in a photon and a proton, respectively. They are evaluated at the factorisation scale $`E_T^2`$. The $`M_{ij}`$ are QCD matrix elements for $`22`$ parton-parton hard scattering processes. The quantity $`s_{ep}`$ is the square of the centre of mass energy in the $`ep`$ collision, and $`\theta ^{}`$ is the polar angle of the outgoing partons in the parton-parton centre of mass frame.

The measurement of the di-jet cross section can be exploited to gain information about the parton densities in the photon. There are two conceptually different approaches, which can be nicknamed “extraction” and “universality test”.

We can extract the parton density in the photon treating the measured di-jet cross section, the parton density in the proton (well constrained by inclusive $`ep`$ DIS data) and the QCD matrix element as input to formula (3). The result of the measurement is a fundamental, parton level quantity. This approach has one important limitation: it is the LO procedure, because it employs the LO notion of the effective parton density and the LO Monte Carlo which relates parton to hadron levels. The H1 collaboration extracted the parton density in photon in the region of the low $`x_\gamma `$ where gluon component of the parton density is of particular interest and poorly constrained by $`e^+e^{}`$ measurements.

In the “universality test” we have different objective: to check if the parametrisation of parton density in photon constrained by $`e^+e^{}`$ measurements is adequate for the description of the inclusive di-jet cross section within NLO QCD theory. In this approach one of the existing parametrisations of the parton density in photon - is used as input information in the NLO calculation  of jet photoproduction cross section at the parton level. In order to be consistent, we have to avoid the application of the LO Monte Carlo for hadronic and soft underlying event corrections and thus restrict the analysis to the region where these corrections are expected to be small: high transverse jet energy hence high $`x_\gamma `$ (see formula (2)). The calculated parton level cross sections are compared with the data at the hadron level, which means that the data are corrected for detector effects only.

It should be underlined that neither approach is satisfactory, especially at the level of the measurement precision recently reached by both HERA experiments. Tools for extraction of the NLO pdf in photons from di-jet measurements, similar to those existing for the proton pdf, are missing. Obviously, there are important technical differences between electromagnetic and strong probes: in the case of $`ep`$ or $`e\gamma `$ DIS the resolution scale does not have to be unfolded for hadronisation and other QCD effects, because it is known directly from the electron measurement. In addition, the strong probe comes in two species (quark or gluon), and this certainly complicates the fitting procedure. The LO solution to this problem, i.e. Single Effective Parton Density approach (see sec. 4.1), may be not applicable in NLO theory.

### 4.1 Measurement of the Effective Parton Density of quasi-real photons

At HERA we cannot distinguish experimentally which partons initiated the hard scattering process in Figure 3, therefore formula (3) is not suitable for the extraction of parton densities. To avoid this difficulty, the H1 Collaboration adopted the Single Effective Subprocess (SES) approximation , which exploits the fact that the dominant contributions to the cross section come from the parton-parton scattering matrix elements which have similar shapes. Introducing the effective pdf in the photon

$$f_\gamma ^{\mathrm{eff}}(x_\gamma ,E_T^2)\underset{n_f}{}(f_{q/\gamma }(x_\gamma ,E_T^2)+f_{\overline{q}/\gamma }(x_\gamma ,E_T^2))+\frac{9}{4}f_{g/\gamma }(x_\gamma ,E_T^2),$$
(4)
the effective proton pdf $`f_p^{\mathrm{eff}}`$ and the SES matrix element $`M_{\mathrm{SES}}`$, we can express the resolved photon contribution to the cross section by the product $`f_\gamma ^{\mathrm{eff}}f_p^{\mathrm{eff}}|M_{\mathrm{SES}}|^2`$. As $`f_p^{\mathrm{eff}}`$ is well constrained by the data and $`M_{\mathrm{SES}}`$ is calculable in perturbative QCD, $`f_\gamma ^{\mathrm{eff}}`$ can be directly determined from the measurement of the di-parton cross section.

In Figure 4 (left) we show the H1 measurement  of the di-jet cross section as a function of the parton momentum fraction $`x_\gamma ^{jet}`$.

The histograms represent a LO QCD MC calculation showing the contributions of the direct photon-proton interactions as well as the quark and gluon induced processes using the GRV pdfs for photon and proton. It should be noted that the relatively low cut on the jet transverse energy ($`E_T^{jets}>`$ 6 GeV) allows a precise measurement down to $`x_\gamma =0.05`$. Sensitivity to the gluon content of the photon is clearly seen. In the definition of $`d\sigma /dx_\gamma ^{jet}`$ pedestal of the jet transverse momentum due to the soft underlying event is subtracted. This is carefully checked by comparing the energy flow around the jets in data and simulation which agrees very well. The remaining differences are added to the systematic errors. In Figure 4 (right) we show the parton density in the photon extracted from the differential cross section. $`x_\gamma ^{jet}`$ has been unfolded to parton level $`x_\gamma `$ and the effective parton density is determined by the comparison of MC and the data

$$f_\gamma ^{DATA}=f_\gamma ^{MC}\frac{\sigma ^{DATA}}{\sigma ^{MC}}.$$
(5)
Figure 5 shows the gluon content of the photon obtained from the effective photon pdf by subtraction of the GRV quark contribution. A good agreement with the gluon density derived from the photoproduction of high transverse momentum charged particles  should be noted, as these two measurements have different systematics. Thus the large gluon component of the photon pdf at low values of $`x_\gamma `$ may be considered as established. The systematic errors dominate above $`x_\gamma 0.1`$, but for low $`x_\gamma `$ values the statistical errors are still important and in principle it is still possible to reach better precision at low $`x_\gamma `$.

Figures 4 and 5 show that at least some parametrisations of the photon pdf constrained by $`e^+e^{}`$ measurements are compatible with jet photoproduction in $`ep`$ as calculated in LO. We can consider this fact as a demonstration of $`e^+e^{}`$-$`ep`$ universality of the photon structure, but taking into account the large errors involved, at the semi-quantitative level only.

### 4.2 Measurement of the effective parton density of the virtual photon

One can study the structure of virtual photons in a similar way to that described in the previous section. Such a study  is presented in Figure 6 (left), based on sample of $`6\mathrm{pb}^1`$ of DIS events in the kinematic range $`1.6Q^280\mathrm{GeV}^2`$, $`0.1y0.7`$ with at least 2 jets with average squared transverse energy $`\overline{E}_t^2>30\mathrm{GeV}^2`$. Jets are restricted to be asymmetric in energy, and the constraints are such that neither jet has $`E_T4\mathrm{GeV}`$ and the sum is always $`11\mathrm{GeV}`$.

In Figure 6 (left) the di-jet cross section is plotted as a function of $`x_\gamma ^{jets}`$ for different regions of $`Q^2`$ and $`E_T^2`$. One sees a clear excess over the expectation of direct photon reactions, indicating that virtual photons also have a resolved part. The direct interactions are manifest in the region of $`x_\gamma ^{jets}1`$. With increasing $`Q^2`$ at fixed $`E_T^2`$ the direct contribution increases while the resolved is significant only for $`E_T^2>Q^2`$ i.e. when the spatial extent of the photon exceeds the resolution of the probe. This is a beautiful demonstration of QFT at work.

It should be noted that pQCD does not provide a full theoretical description of the virtuality dependence of the parton density in the photon. There is only an asymptotic prediction  in LO for $`F_2^\gamma ^{}`$ i.e. the structure function of the virtual photon $`\gamma ^{}`$

$$F_2^\gamma =3\underset{q}{}e_q^4\frac{\alpha }{\pi }x[x^2+(1x)^2]\mathrm{ln}\frac{Q^2}{P^2},$$
(6)
valid for $`Q^2P^2\mathrm{\Lambda }^2`$. Here $`P^2`$ is the photon virtuality while $`Q^2`$ is the virtuality of the probe (resolution scale). The formula (6) is analogous to (1) for a quasi-real photon. Now, kinematics provides its own cut-off $`P^2`$ as long as it exceeds $`\mathrm{\Lambda }^2`$. The “hadron piece” in formula (1) no longer counts as $`F_2^{VDM}1/P^4`$.<sup>2</sup><sup>2</sup>2Note that formulae (1) and (6) refer to the photon structure function as measured in $`e^+e^{}`$ experiments, and there are conflicting notations: in $`e^+e^{}`$ $`Q^2`$ always denotes the virtuality of the probe (i.e. resolution scale) while in $`ep`$ experiments at HERA $`Q^2`$ traditionally denotes the virtuality of the photon which here happens to be a target. Throughout the article we stick to traditional notations, leaving the reader to determine the correct recognition of the symbols.

For the description of the triple differential di-jet cross section we can employ phenomenological models of the dependence of the effective parton density in the photon on its virtuality $`Q^2`$ , , . The Drees-Godbole model (DG)  starts with parton densities in the real photon and suppresses them by a factor $``$

$$(Q^2,E_T^2,\omega ^2)=\frac{\mathrm{ln}\frac{E_T^2+\omega ^2}{Q^2+\omega ^2}}{\mathrm{ln}\frac{E_T^2+\omega ^2}{\omega ^2}},$$
(7)
where $`\omega `$ is a free parameter. The quark densities in the real photon are suppressed by $``$ and the gluon densities by $`^2`$. This ansatz, based on analysis in , is designed to interpolate smoothly between $`\mathrm{ln}(Q^2)`$ and $`\mathrm{ln}(Q^2/P^2)`$ in the LO asymptotic formulae for real (1) and virtual (6) photons. In the model of Schuler and Sjöstrand , the perturbative anomalous component of the parton density in a virtual photon is designed to approach the asymptotic result (6) as in the DG model, while the VDM part is rapidly suppressed by the factor $`m_V^2/(m_V^2+Q^2)`$. In this approach the parton density in a quasi-real photon has to be initially decomposed into perturbative and non-perturbative components, so that the description of the $`Q^2`$ suppression of the photon pdf is inseparable from the parametrisation of the quasi-real photon pdf. This is not the case for the DG model, which can be applied to any quasi-real photon pdf.

In Figure 6 (left) the predictions of the HERWIG MC with GRV photon pdf multiplied by the DG factor are compared with the data. Clearly the data are well described by the MC with suitable choice of the parameter $`\omega =0.2`$ which controls the onset of the $`Q^2`$ suppression.

In Figure 6 (right) we show the effective parton density of the virtual photon extracted from the triple differential di-jet cross section (on the left). In the whole measurement range the systematic errors dominate the total error. The $`x_\gamma `$ dependence of the effective parton density shows a tendency to rise with increasing $`x_\gamma `$ as predicted by formula (6). For quasi-real photons in a similar $`x_\gamma `$ range the pdf is flat. This fact can be interpreted as the effect of the VDM contribution, which preferentially fills out the low $`x_\gamma `$ region flattening the otherwise rising $`x_\gamma `$-dependence of the anomalous pQCD contribution. With rising $`Q^2`$, the VDM contribution rapidly decreases, leaving only the rising $`x_\gamma `$-dependence. The errors on the experimental data are clearly too large to differentiate between SAS and DG, although only SAS has a rising tendency for reasons explained before (the DG model has no explicit VDM component).

In Figure 7 (left) we show the extracted photon pdf as a function of $`Q^2`$. The data are compared with the pure VDM prediction and the DG and SAS models discussed above. The DG and SAS parametrisations describe the data up to $`Q^225\mathrm{GeV}^2`$. The last data point $`25Q^280\mathrm{GeV}^2`$ is above both parametrisations. This is however the region where $`Q^2E_T^2`$ and non-leading terms not accounted for in the models are expected to be important and may affect the extraction of the effective parton distribution from the data (breakdown of factorisation).

The ZEUS collaboration investigated the photon pdf suppression with increasing $`Q^2`$ in a different way. They use an operational definition of direct and resolved di-jet cross sections: $`\sigma ^{dijet}(\mathrm{direct})=\sigma ^{dijet}(x_\gamma ^{jet}>0.75),`$ $`\sigma ^{dijet}(\mathrm{resolved})=\sigma ^{dijet}(x_\gamma ^{jet}<0.75).`$ In Figure 7 (right) we show the ratio of the so-defined “resolved” to “direct” contributions as a function of $`Q^2`$. The SAS prediction is below the data points, but the shape of the dependence is correct.

At the end of this section we would like to underline that there is no basic difficulty in the theoretical description of the region $`E_T^2Q^2`$. On the contrary, in this region of phase space, terms with $`\mathrm{ln}(E_T/Q)`$ appearing in the perturbative expansion are small, so that even fixed order theory is applicable in this case. It is only the photon structure function approach which breaks down for $`QE_T`$. Let us also note that in the range $`E_T^2Q^2\mathrm{\Lambda }^2`$ the photon structure approach can, in principle, be replaced by resummed $`(\mathrm{ln}(E_T^2)/Q^2))^n`$ theory, because for the photon, unlike for hadrons, the purely perturbative (anomalous) part of pdf dominates completely at high values of $`E_T^2`$. In this sense the structure function of the virtual photon could and should be considered as a convenient parametrisation of higher order QCD corrections which are not available yet.

## 5 Comparison of the measured di-jet cross section with a NLO QCD calculation

The ZEUS Collaboration analysis of di-jet photoproduction , , based on a large sample of events, confronts a precision measurement with a NLO QCD calculation. The cut on the transverse energy of jets $`E_T^{jet1,2}>14,11\mathrm{GeV}`$ is so high that corrections for hadronisation and the soft underlying event are certainly small (less than $`10\%`$). Therefore the data, corrected only for detector effects, are compared directly to NLO QCD calculations  at the parton level. There is good general agreement of the di-jet cross section with the NLO calculation as seen e.g. in Figure 9 where it is plotted as a function of the transverse energy. While $`E_{Tleading}^{jet}`$ increases from 14 to 50 GeV, cross section falls by 3 orders of magnitude and no obvious deviations from theory are seen, except for events with forward jets $`1<\eta _{1,2}^{jet}<2`$ and $`E_{Tleading}^{jet}<25\mathrm{GeV}`$. Let us note that the theory–data discrepancy appears exactly in the region where the distance between “direct” and “direct + resolved” points is large, meaning that the influence of the resolved component is large. This tendency of the data becomes even more evident when the di-jet cross section is plotted as a function of the photon momentum fraction . In Figure 9 (left) the di-jet cross section is plotted as a function of $`x_\gamma ^{obs}x_\gamma ^{jet}`$ in several bins of the transverse energy of the leading jet. Except for the last bin (dominated by direct process) all data points are above the NLO cross section.

Figure 9 (right) shows the relative value of the theory-data discrepancy as a function of $`x_\gamma ^{obs}`$. The relative discrepancy reaches 60$`\%`$. Even taking into account the uncertainty of the theory due to the choice of the NLO scale (shown as shaded area) and the experimental uncertainty due to the hadronic energy scale, we have to admit that there is an indication that the parton densities in the photon are underestimated. Incidentally, this conclusion seems to be confirmed by recent measurements of the di-jet cross section in $`\gamma ^{}\gamma ^{}`$ collisions  . Even so it is difficult to accept without reservation that the pdf in photon is much underestimated. The argument against such interpetation of ZEUS and OPAL results is the general good agreement of the LEP measurements of photon structure function with GRV parametrisation. This problem certainly deserves further scrutiny.

## 6 Summary and Conclusion

Measurements of $`e\gamma `$ DIS and jet electroproduction at HERA confirm the simple picture of photon structure first identified in $`e^+e^{}`$ experiments at PETRA. There is general compatibility of $`e^+e^{}`$ and $`ep`$ measurements which demonstrates the validity of a factorisable structure which can be assigned to the photon, both real and virtual. Properties of the photon structure function are in agreement with theoretical expectations, in particular its dependence on the photon momentum fraction and resolution scale. Also the parton density suppression with increasing photon virtuality is in agreement with phenomenological models based on asymptotic QCD predictions.

The recent precision measurements of the di-jet cross section seem to indicate that the existing parametrisations (based on $`e^+e^{}`$ measurements) of the photon pdf are inadequate for a precise description of the data using NLO theory. This problem certainly deserves further scrutiny.

The accumulated LEP, LEP2 and HERA data with much reduced systematic uncertainties over a very wide kinematic range now wait for a combined analysis, which would result in a new NLO parametrisation of the photon pdf. The concepts and tools for such an analysis have yet to be developed.


text:
# EFFICIENT QUANTUM COMPUTING INSENSITIVE TO PHASE ERRORS

## 1 Introduction

The problem of quantum computation has attracted recently a great deal of attention $`^{\mathrm{?},\mathrm{?},\mathrm{?}}`$. This interest stems from the fact that the massive parallelism permitted by quantum mechanics enables to reach exponential efficiency of computation in certain quantum algorithms. The most famous example is the Shor algorithm which allows to factor large numbers exponentially faster than any known classical algorithm $`^\mathrm{?}`$. Recently, other types of exponentially efficient algorithms has been developed for the simulation of various physical systems $`^{\mathrm{?},\mathrm{?},\mathrm{?}}`$. The exponential gain in computation is related to the exponentially large Hilbert space of the quantum computer which is composed of multi-qubit states operating in parallel (each qubit is a two-level quantum system). Usually the algorithms are constructed for ideal quantum computers operating free of noise and imperfections. In reality, any physical realization of such a computer involves a certain level of imperfections, noise in gate operations and decoherence. First investigations showed that quantum computation can tolerate a sufficiently low level of errors $`^{\mathrm{?},\mathrm{?}}`$. More recently, it has been found that quantum computation is tolerant to quantum errors when simulating classical chaotic dynamics for which classical computer errors grow exponentially with time $`^\mathrm{?}`$. However in general the quantum errors grow with the number of gate operations and any realistic quantum computer is faced with this problem. To deal with this problem of fault-tolerant computation, quantum error-correcting codes were recently developed $`^{\mathrm{?},\mathrm{?},\mathrm{?}}`$. They allow to reduce the level of errors in a systematic way, but require the introduction of many redundant qubits and additional gates, which significantly complicates the computational process. The complexity of these codes depends strongly on the type of errors they should correct. Indeed, while simpler classical codes need to correct only bit errors, the quantum ones should in addition simultaneously correct the quantum phase errors. The quantum phase errors seem to be of primary importance since the massive parallelism of quantum computing is related to entanglement in the Hilbert space which is directly based on phase coherence. Therefore, according to this common lore it seems impossible to perform efficient and accurate quantum computations in presence of uncontrolled strong phase errors. In this paper we show on an explicit example that it is not always the case, and that there are algorithms insensitive to phase errors which perform accurate and efficient computation. Our example is based on the simulation of classical chaotic dynamics, which is very hard to simulate accurately on classical computers since this dynamics is unstable and round-off errors grow exponentially with time. In spite of this, our quantum algorithm, including measurement, remains insensitive to phase errors for arbitrary time.

## 2 The Arnold-Schrödinger cat quantum algorithm

To illustrate this phenomenon, we choose an algorithm which simulates the classical chaotic dynamics of the well-known Arnold cat map $`^{\mathrm{?},\mathrm{?}}`$. It was recently shown $`^\mathrm{?}`$ that quantum computers can simulate this dynamics with exponential efficiency. In addition it was shown that a small level of quantum errors in the gate operations of order $`ϵ`$ allows to simulate accurately this map on times of order $`O(1/ϵ^2)`$. Thus quantum computers can face classical exponential instability and chaos.

The dynamics of the map we are studying is given by:

$$\overline{y}=y+x\text{(mod}\text{1)},\overline{x}=y+2x\text{(mod}\text{1)},$$
(1)
where bars denote the new values of the variables after one iteration. This is an area-preserving map, in which $`x`$ can be considered as the space variable and $`y`$ as the momentum. A discretized version on a $`N\times N`$ square lattice is also described by this map. In $`^\mathrm{?}`$, a quantum algorithm called Arnold-Schrödinger cat map was introduced, and it was shown to simulate this dynamics on the lattice with exponential efficiency. In this paper, we modify this algorithm in such a way that exponential efficiency is preserved and in addition it becomes insensitive to phase errors. This is obtained by the introduction of a new measurement procedure.

The quantum algorithm introduced in $`^\mathrm{?}`$ simulates the discrete classical dynamics given by (1) and operates with $`3n_q1`$ qubits. These qubits are organized in three registers. Two of them describe the classical phase space with $`N^2`$ points and $`N=2^{n_q}`$. The third register with $`n_q1`$ qubits is used as workspace. In this way, an initial classical phase space density can be represented by a quantum state $`_{i,j}a_{ij}|x_i>|y_j>|0>`$, with $`x_i=i/N`$, $`i=0,\mathrm{},N1`$ and $`y_j=j/N`$, $`j=0,\mathrm{},N1`$, written in binary representation, and we choose initially $`a_{ij}=0`$ or $`1/\sqrt{N_d}`$ where $`N_d`$ is the number of points in the classical distribution. Then, iterations of the map (1) are performed with the help of additions of integers modulo (N) (modular additions). The quantum algorithm we use for this operation is similar to the one described in $`^\mathrm{?}`$ (see also $`^\mathrm{?}`$). First we compute all the carries of the addition, using two Toffoli gates and one controlled-not (CNOT) gate per qubit. Then we perform the addition starting from the last qubit and erasing the carries by running the inverse of the preceding step. This needs two CNOT gates per qubit addition and the same gates as above to erase the carries. The result is taken modulo (N) by eliminating the last carry. After these operations, the amplitudes $`|a_{ij}|`$ describe the classical phase space distribution function after iteration of (1). In total, one needs $`16n_q22`$ Toffoli and CNOT gates per map iteration. On the contrary, a classical computer will require $`O(2^{2n_q})`$ operations per iteration for $`N_d=O(N^2)`$ orbits.

It is important to stress that during the whole process the classical distribution function is determined only by the probabilities $`|a_{ij}|^2`$ of the quantum computer wavefunction expanded on the Hilbert space basis of register states $`|x_i>|y_j>`$ (after each map iteration the third register is always in the state $`|0>`$). Thus, the information about the classical distribution function is stored in these probabilities, and is not sensitive to the relative quantum phases of $`a_{ij}`$. This suggests that the phase errors accumulated during gate operations will not affect the quantum computer simulation of this classical chaotic dynamics.

## 3 Effect of phase errors

To study the effects of quantum errors, one usually uses the fidelity $`^\mathrm{?}`$, defined as: $`f(t)=`$

$`|<\psi _ϵ(t)|\psi _0(t)>|^2=|_{i,j}a_{ij}^{(ϵ)}(t)a_{ij}^{(0)}(t)|^2`$. Here $`|\psi _0(t)>=_{i,j}a_{ij}^{(0)}(t)|x_i>|y_j>`$ is the quantum state after $`t`$ perfect iterations, while $`|\psi _ϵ(t)>=_{i,j}a_{ij}^{(ϵ)}(t)|x_i>|y_j>`$ is the quantum state after $`t`$ imperfect iterations. Obviously, this quantity is very sensitive to the relative phases of $`a_{ij}`$. Since the classical phase space density is not sensitive to these phases, we introduce another characteristic which is related only to the amplitudes $`|a_{ij}|`$. We call it faithfulness and define it by: $`f_a(t)=(_{i,j}|a_{ij}^{(ϵ)}(t)a_{ij}^{(0)}(t)|)^2`$. This quantity can be considered as a generalization of the usual fidelity. As well as $`f(t)`$, the faithfulness $`f_a(t)`$ is always $`1`$, and it determines the deviation from the exact amplitudes (the value $`1`$ is reached only for $`|a_{ij}^{(ϵ)}(t)|=|a_{ij}^{(0)}(t)|`$ for all $`i,j`$). In addition, one has always $`f_a(t)f(t)`$. Contrary to the usual fidelity, the faithfulness measures only amplitude errors, being insensitive to the quantum phases. We note that its definition is related to the preferential basis chosen initially in the Hilbert space.

To study the dependence of fidelity and faithfulness on phase errors, the nondiagonal part of each Toffoli and CNOT gate used in the algorithm was multiplied by a diagonal matrix with elements $`\mathrm{exp}(i\theta _m)`$, with random phases $`\theta _m`$ homogeneously distributed in $`[ϵ_\varphi ,ϵ_\varphi ]`$. The results are shown in Fig. 1 (Left). Here the initial state was chosen in the form of a cat’s smile (see Fig. 1 of $`^\mathrm{?}`$ and the coarse-grained version in Fig. 2). In the presence of phase errors only, the fidelity $`f`$ decreases with the number of map iterations and drops almost to zero for sufficiently strong phase noise. At the same time, the faithfulness $`f_a`$ is not affected even by the maximal possible phase noise. We also checked that $`f_a=1`$ is not affected if each $`a_{ij}`$ is multiplied after each gate by a random phase $`\mathrm{exp}(i\theta _m)`$ with $`\theta _m`$ distributed in $`[\pi ,\pi ]`$, although in this case the fidelity is almost zero after one map iteration.

To show the difference between phase and amplitude (bit) errors, we computed the faithfulness in presence of a small unitary noise in the gates. For that, in addition to large phase errors, for each gate the nondiagonal part was diagonalized, and each eigenvalue was multiplied by a random phase $`\mathrm{exp}(i\eta )`$, with $`ϵ<\eta <ϵ`$. This introduces both phase and amplitude errors, and Fig. 1 (Left) shows that the faithfulness starts to drop slowly with time. Hence despite the presence of strong phase errors the faithfulness is sensitive only to the amplitude errors.

Thus the information stored in the amplitudes is not sensitive to phase decoherence. Still, one should find a way to retrieve a part of this exponentially large information. Usually one performs a quantum Fourier transform (QFT) and measures the maximal harmonics of the distribution, as was suggested in $`^\mathrm{?}`$. However, the QFT is extremely sensitive to the quantum phases of $`a_{ij}`$, as is illustrated in Fig. 1 (Right). The zero harmonic $`Q_0(t)=_{i,j}a_{ij}(t)/N=\sqrt{N_d}/N`$ is time-independent in the absence of errors, but drops rapidly with $`t`$ if phase errors are present. To avoid the effects of phase errors, one can measure only the first $`n_g`$ qubits from the $`n_q`$ qubits present in the register $`|x>`$ and the same for the register $`|y>`$. This procedure introduces a coarse-graining of the phase space ($`x,y`$), with the number of cells $`N_g=2^{2n_g}`$. The result of such a measurement is determined by the total probability inside each cell $`W_{i_gj_g}=_{<i,j>}|a_{ij}|^2`$ where the summation is performed over all ($`i,j`$) inside the cell ($`i_g,j_g`$). This probability is not sensitive to phase errors, and can be extracted by a number of measurements which is polynomial in $`N_g`$. Fig. 1 (Right) shows that the probability in a chosen cell ($`i_g^{},j_g^{}`$) is indeed insensitive to phase errors. We note that this coarse-grained probability is a very natural quantity for the dynamical system under investigation. One is not interested in the exponential amount of information present in all $`a_{ij}`$ since one cannot even store it classically, and therefore it is better to operate with coarse-grained characteristics as is usually done in chaotic dynamical systems. The number of cells $`N_g`$ can be kept constant while the number of iterated classical orbits increases exponentially with $`n_q`$. All these orbits are iterated accurately and with exponential efficiency during quantum computation, and the constant number of cells $`N_g`$ is used only to extract the essential information generated by this chaotic dynamics.

Another test of the effects induced by phase and amplitude errors can be performed on the basis of time-inversion. Indeed, the exact map (1) is exactly time-invertible. However, in presence of imperfections this reversibility can be destroyed. In the régime of classical chaos, the classical round-off errors grow exponentially with time and destroys time-reversibility in a logarithmically short time. For quantum simulations, it was shown in $`^\mathrm{?}`$ that quantum errors grow only polynomially with time. Due to that, time-reversibility is preserved in quantum computation of (1) for relatively small errors. However, it is naturally expected to be destroyed in the case of strong errors. Contrary to this expectation, Fig.2 (Right) shows that time-reversibility is exactly preserved in the presence of phase errors of maximal amplitude, and the classical distribution is exactly reproduced for all $`t`$. In the right column, the difference between the two images at $`t=0`$ and $`t=100`$ is on the level of classical computer precision. On the contrary, strong enough amplitude errors completely destroy time-reversibility, as is shown on Fig.2 (Left).

## 4 Discussion

Thus, all the data clearly show that our quantum algorithm simulating (1) is insensitive to phase errors. This result can be understood in the following way. All nondiagonal parts of the gates used in the algorithm are represented by the operator $`\sigma ^x`$, while the noncommuting part of phase errors is represented by $`\sigma ^z`$. Of course, $`\sigma ^x`$ and $`\sigma ^z`$ do not commute. However, the action of $`\sigma ^x`$, $`\sigma ^z\sigma ^x`$ and $`\sigma ^x\sigma ^z`$ on a two-component spinor gives the same amplitudes of the components (with different relative phases). Thus any quantity encoded in the amplitudes, in our case the classical distribution function, remains invariant in presence of $`\sigma ^z`$ (phase) errors. On the contrary, it is sensitive to $`\sigma ^x`$ (amplitude) errors. Another way of understanding this insensitivity to phase errors is to remark that all used gates belong to a very specific subgroup among unitary transformations of the Hilbert space, that is the group of permutations of the basis formed by the states where each qubit is polarized in the $`z`$ direction (each qubit is either $`|0>`$ or $`|1>`$). The amplitudes in this basis are insensitive to phase errors if only such transformations are present in the algorithm. Indeed, any permutation can be written as a product of transpositions which exchange only two states. By the same argument as for $`\sigma ^x\sigma ^z`$ given above, such a transformation is immune to phase errors, and hence any permutation. We stress again that phase errors do affect the final state through the relative phases, but do not affect the measurement which gives the cell probabilities $`W_{i_gj_g}`$.

The above mathematical argument explains the insensitivity to phase errors. In a more physical way, we can say that the map (1) describes the classical dynamics of the Arnold cat map, which naturally should not be sensitive to quantum phases. Of course, one can imagine other quantum algorithms which will simulate this classical dynamics using both phases and amplitudes of the wave function, and therefore will be sensitive to phase errors. However, on the basis of the Arnold-Schrödinger cat algorithm discussed in this paper, we make the conjecture that classical Hamiltonian dynamics of generic systems can always be simulated on a quantum computer in a way insensitive to phase errors. Indeed, for such a dynamics the classical information can be naturally encoded in the amplitudes only $`^\mathrm{?}`$. It is rather likely that such a situation can appear in quantum computations which are not connected with classical mechanics, for example probing the range of values of a function.

The implementation of such algorithms insensitive to phase errors can be enormously simpler than in the case of other algorithms sensitive to quantum phases. Indeed, the necessity to correct both phase and amplitude errors significantly complicates quantum error-correcting codes $`^{\mathrm{?},\mathrm{?},\mathrm{?}}`$. If only amplitude errors are to be corrected, one can use much simpler codes close to the classical ones. Also, in some physical systems phase errors can be naturally much stronger than amplitude ones. For example, recent studies of the emergence of quantum chaos in a quantum computer $`^\mathrm{?}`$ showed that for sufficiently strong residual inter-qubit interaction, exponentially many states are mixed and amplitude errors become enormously strong. On the contrary, below the quantum chaos border, amplitude errors are very small whereas phase errors are still important. In spite of that, a quantum computer in this régime can efficiently simulate algorithms of the type discussed here.

We note that while the algorithm presented above is exponentially faster than any deterministic classical algorithm iterating the map (1) nevertheless one can try to compete with it with the help of classical Monte Carlo simulation with a polynomial number of trajectories. Such an approach does not produce the exact density distribution with an exponential number of orbits which is hidden in the quantum wavefunction. However, the statistical accuracy of both methods can be comparable since one makes a polynomial number of measurements of the quantum final state. At the same time one should keep in mind that such a Monte Carlo simulation is based on the statistical assumption that a polynomial number of trajectories can correctly describe the fine structure of classical phase space. In contrast, the quantum simulation takes exactly into account the dynamics on all scales. We also stress that without large phase errors the QFT gives access to information which is unaccessible even for classical Monte Carlo algorithms.

In conclusion, we have shown the existence of quantum algorithms which can simulate efficiently certain computational problems and at the same time are insensitive to phase errors. Our explicit example is related to the simulation of classical motion and we make the conjecture that classical Hamiltonian dynamics can always be simulated in a way immune to phase decoherence. The existence of such efficient quantum algorithms insensitive to the relative phases shows that contrary to the common lore, the massive parallelism of quantum computing is not necessarily related to quantum interference. Actually, quantum mechanics allows to follow in parallel exponentially many computational paths in a way insensitive to phase decoherence.

## Acknowledgments

We thank C. W. J. Beenakker for discussions and the IDRIS in Orsay and the CalMiP in Toulouse for access to their supercomputers. This work was supported in part by the EC RTN contract HPRN-CT-2000-0156.

## References


text:
# DYNAMICS OF SOLITONS AND QUASISOLITONS OF CUBIC THIRD-ORDER NONLINEAR SCHRÖDINGER EQUATION

\[

## Abstract

The dynamics of soliton and quasisoliton solutions of cubic third order nonlinear Schrödinger equation is studied. The regular solitons exist due to a balance between the nonlinear terms and (linear) third order dispersion; they are not important at small $`\alpha _3`$( $`\alpha _3`$ is the coefficient in the third derivative term) and vanish at $`\alpha _30`$. The most essential, at small $`\alpha _3`$, is a quasisoliton emitting resonant radiation (resonantly radiating soliton). Its relationship with the other (steady) quasisoliton, called embedded soliton, is studied analytically and in numerical experiments. It is demonstrated that the resonantly radiating solitons emerge in the course of nonlinear evolution, which shows their physical significance.

PACS numbers: 42.65.Tg, 52.35.Mw

\]

I. INTRODUCTION

Equations describing soliton processes are usually obtained by certain approximating procedures affecting nonlinearity and dispersion. Now, in extensive studies of ultrafast processes, the classical approximations often appear to be insufficient and higher order effects become of importance. A typical example is high speed systems like nonlinear transmission lines in femtosecond regime for soliton communications, etc. For such systems soliton solutions in classical sense do not generally exist. Only for very specific choices of parameters one can find localized solutions. Instead of regular solitons there may appear nonlocal steady or (and) unsteady soliton like structures (which may be called quasisolitons). These may, however, have significant importance for the nonlinear dynamics.

In this paper we consider dynamics described by the extended third-order cubic nonlinear Schrödinger (NLS) equation

$$i_T\mathrm{\Psi }+\frac{1}{2}_X^2\mathrm{\Psi }+|\mathrm{\Psi }|^2\mathrm{\Psi }+i\alpha _1|\mathrm{\Psi }|^2_X\mathrm{\Psi }+i\alpha _2\mathrm{\Psi }_X|\mathrm{\Psi }|^2$$

$$+i\alpha _3_X^3\mathrm{\Psi }=0,$$
$`(1)`$
which plays an important role in many nonlinear problems, in partiqular, in the nonlinear fiber optics . The third derivative term describes higher order correction to the linear dispersion and the cubic terms with the first derivative are dispersive nonlinear corrections. Here, we assume that all $`\alpha _n`$ ( $`1\alpha _n3`$) are real parameters.

Eq.(1) has regular soliton solutions \[3-5\], vanishing at $`|X|0`$. In particular cases, Eq.(1) may have exact N \- soliton solutions  or even be integrable . An important feature of the solutions describing regular solitons is that they are degenerating at $`\alpha _30`$ (see below), which substantially reduces their physical importance. Their numerical investigation at $`\alpha _{1,2}1`$ where done e.g. in Ref. (see also references therein)

There also exist some solutions to Eq.(1), describing quasisolitons. One of them is a steady solution looking like a soliton type pulse embedded into a small amplitude plane wave (a soliton on plane wave pedestal). We call it embedded soliton (ES), using the word proposed in a different context in Ref.. Apart from that, Eq. (1) has other type of quasisoliton solutions which describe soliton-like pulses permanently emitting resonantly generated radiation (see, e.g., Refs. \[8-12\] for the particular case at $`\alpha _1=\alpha _2=0`$ and  for the full Eq.(1)). We will call them resonantly radiating solitons (RRS); they are unsteady because of losses caused by the radiation. The lifetime of radiating soliton is sufficiently large if $`\alpha _3`$ is small enough and, naturally, one can speak about the soliton only in this quasisteady case. (However, at long times, the losses caused by the radiation become essential for applications .) At $`\alpha _30`$, the RRS turns into the regular soliton of Eq.(1) without the third derivative.

There is an interesting and important connection between the two types of quasisolitons, ES and RRS, and we demonstrate it in numerical simulations. We will also study the role of solitons and quasisolitons in the nonlinear processes described by Eq.(1). In partiqular, the regular solitons and RRS compete between themselves in nonlinear processes. As far as the regular solitons disappear at $`\alpha _30`$ and the RRS have short lifetimes at large $`\alpha _3`$, it is clear that at large $`\alpha _3`$ the regular solitons are more important while the RRS may play a decisive role at sufficiently small $`\alpha _3`$. The most interesting case is, of course, small $`\alpha _3`$ , because the third derivative term emerges as the result of an expansion. (The effect of next, 4th derivative, term can be seen in Refs.; at certain relalionship between the coefficients before third and fourth derivatives there can be no radiation at all which, naturally, may happen also when other high-order terms are taken into account.)

The paper is organized as follows. In Sec.I we describe important properties of Eq.(1), which are used below. In particular, we discuss regular soliton solutions , Galilean transformations and conservation laws . The embedded solitons are studied in Sec.3. As we have already mentioned, the ES consists of a steady soliton-like pulse on the plane wave background. This structure is rater common in different nonlinear highly dispersive systems \[9-20\]. At $`\alpha _30`$ , the plane wave disappears and the pulse turns into a regular soliton of Eq. (1) at $`\alpha _3=0`$ . The plane wave amplitude, increasing with $`\alpha _3`$ , may become unstable at sufficiently large $`\alpha _3`$ due to the modulational instability. Considering small $`\alpha _3`$, we show that the pulse part of ES is rather close to the pulse in RRS and the wavenumber of the plane wave coincides with the wavenumber of resonantly emitted radiation (by the RRS) . Next, we introduce the cut off operation, cutting off both wave wings from the ES (Sec.IV). Then we see that the remained soliton-like pulse is transformed into the RRS, emitting radiation only in one direction , according to the direction of the group velocity of resonant radiation. The resonant radiation disappears at $`\alpha _1=6\alpha _3`$. In this case, considered in Sec.V, we show that Eq.(1) can be transformed, by means of the Galilean transformation, to the complex modified Korteweg-de Vries (MKdV) equation . At particular initial conditions, it is reduced to the real MKdV equation which is integrable and therefore has N \- soliton solutions . This does not mean the complete integrability of Eq.(1) at $`\alpha _1=6\alpha _3`$ , because the reduction to the real MKdV equation is possible only for particular initial conditions. However, the Painlevé analysis  of the above mentioned complex MKdV equation shows that at $`6\alpha _3=\alpha _1=2\alpha _2`$ it possesses the Painlevé property. This is just the case when Eq. (1) is integrable . If $`\alpha _2=0`$, we arrive at the case considered by Hirota  who has shown the existence of N \- soliton solutions in this case. Our approach, however, is more general and (in fact, as in Ref.  ) we can see that the N -soliton solutions are possible even for nonvanishing $`\alpha _2`$ (at $`\alpha _23\alpha _3`$ ). Then we present a numerical study of the initial value problem in the Hirota case:$`\alpha _1=6\alpha _3`$, $`\alpha _2=0`$ taking as initial conditions complex pulses and one can see how the robust Hirota solitons are emerging from these pulses. In Sec.VI, we describe numerical experiments at $`\alpha _16\alpha _3`$. Solving the initial value problem, we show that, at sufficiently small $`\alpha _3`$, the initial disturbances decay into the resonantly radiating solitons (RRS). On the other hand, at small $`\alpha _3`$, the Potasek- Tabor solitons  were not detected. This indicates that the RRS are important physical objects playing a significant role in nonlinear dynamics. In Sec. VII , a summary of obtained results is given. In Appendix A, the modulational instability of a plane wave, described by Eq.(1), is considered; it is helpful in the study of the stability of embedded solitons. In Appendix B we investigate analytically, by means of conservation laws, the evolution of RRS caused by the radiation. This analysis is in agreement with numerical results described in Secs.IV and VI.

II. IMPORTANT PROPERTIES OF EQ. (1)

a. Exact soliton solutions

First, we discuss the exact soliton solutions of Eq.(1). They can be written in the form

$$\mathrm{\Psi }_s=asech[b(XV_sT)]e^{ikXi\omega T},$$
$`(2a)`$
$$a^2=\frac{6\alpha _3}{\alpha _1+2\alpha _2}b^2(\alpha _12\alpha _2),$$
$`(2b)`$
$$V_s=\kappa 3\alpha _3\kappa ^2+\alpha _3b^2,$$
$`(2c)`$
$$\omega =\frac{1}{2}\kappa ^2\alpha _3\kappa ^3\frac{16\alpha _3\kappa }{2}b^2,$$
$`(2d)`$
$$\kappa =\frac{\alpha _1+2\alpha _26\alpha _3}{12\alpha _2\alpha _3}.$$
$`(2e)`$
In fact, these are the soliton solutions found by Potasek and Tabor  (with corrected misprints). We will call Eqs. (2a-e) Potasek-Tabor (PT) soliton solutions. At $`\alpha _3=0`$, PT solitons do not exist; this signifies that their existence is a result of a ballance between the nonlinear terms and the linear third order dispersive term . If $`\alpha _20`$, from Eq. (2e) it follows that $`\kappa `$ is a fixed number, uniquely determined by the coefficients of Eq.(1). However, the solitons (2) exist even at $`\alpha _2=0`$ provided that

$$\alpha _1=6\alpha _3.$$
$`(3)`$
In this case $`\kappa `$ can be arbitrary , because Eq. (2e) follows from the equation $`12\alpha _2\alpha _3\kappa =\alpha _16\alpha _3+2\alpha _2`$.

In case (3) and $`\alpha _2=0`$, Eq. (1) is the so called Hirota equation that can be transformed to the complex modified Korteweg-de Vries equation which has N -soliton solutions . For the Hirota soliton, from Eq.(2b) it follows

$$b=a.$$
$`(4)`$
In the other particular case, $`6\alpha _3=\alpha _1=2\alpha _2`$, Eq. (1) is integrable . (See also Sec.V.) Some results for solitons (2) with $`\alpha _20`$ were reported in Ref. (see also references therein) . In particular, it was shown numerically that they emerge in a solution of initiall value problem. However, this result was obtained at $`\alpha _1=\alpha _3=1`$; at small $`\alpha _3`$, as it will be shown below, they cannot compete with RRS.

Apart from the solution (2a-e), describing ”bright” solitons, which was called in Ref.  $`sech`$-family, Potasek and Tabor have also found $`tanh`$-family which will not be discussed here. (Some generalizations, with both $`sech`$ and $`tanh`$ terms, are derived in Ref.  by means of rather tedious algebra; the results of this paper are mostly contained in Ref. and the generalizations , which are valid for very specific choices of the parameters, can be obtained by the much simpler approach of Potasek and Tabor.) The results of Ref. are repeated also in some other papers (e.g.).

b. Galilean transformation

In the following we shall frequently use the Galilean transformation describing the transition to the reference frame moving with the velocity V . For Eq. (1) it reads

$$\mathrm{\Psi }(X,T)=\psi (XVT,T)\mathrm{exp}[i(KX\mathrm{\Omega }T)],$$
$`(5)`$
where $`K`$ and $`\mathrm{\Omega }`$ are defined by equations

$$V=K3\alpha _3K^2,$$
$`(6)`$
$$\mathrm{\Omega }=\frac{1}{2}K^2\alpha _3K^3.$$
$`(7)`$
The function $`\psi (x,t)`$ with

$$x=XVT,t=T$$
$`(8)`$
satisfies the equation

$$i_t\psi +\frac{1}{2}a_2_x^2\psi +q|\psi |^2\psi +i\alpha _1|\psi |^2_x\mathrm{\Psi }+i\alpha _2\psi _x|\psi |^2$$

$$+i\alpha _3_x^3\psi =0,$$
$`(9)`$
with

$$a_2=16\alpha _3K,$$
$`(10)`$
$$q=1\alpha _1K.$$
$`(11)`$
c. Conservation laws

One can check by means of straightforward calculations that Eq. (1) has the following integral of motion

$$N=_{\mathrm{}}^{\mathrm{}}|\mathrm{\Psi }(X,T)|^2𝑑X.$$
$`(12)`$
At

$$\alpha _2=0,$$
$`(13)`$
we have two other conserving integrals

$$P=\frac{1}{2}i_{\mathrm{}}^{\mathrm{}}\left(\mathrm{\Psi }_X\mathrm{\Psi }^{}\mathrm{\Psi }^{}_X\mathrm{\Psi }\right)𝑑X,$$
$`(14)`$
$$H=_{\mathrm{}}^{\mathrm{}}\{\frac{1}{2}|_X\mathrm{\Psi }|^2\frac{1}{2}|\mathrm{\Psi }|^4\frac{1}{4}i\alpha _1|\mathrm{\Psi }|^2(\mathrm{\Psi }^{}_X\mathrm{\Psi }\mathrm{\Psi }_X\mathrm{\Psi }^{})$$

$$\frac{1}{2}i\alpha _3(\mathrm{\Psi }^{}_X^3\mathrm{\Psi }\mathrm{\Psi }_X^3\mathrm{\Psi }^{})\}dX.$$
$`(15)`$
III. EMBEDDED SOLITONS

We start with the solutions to Eq. (1) of the form

$$\mathrm{\Psi }(X,T)=\chi (x)e^{i\mathrm{\Lambda }t},$$
$`(16)`$
where $`x`$ and $`t`$ are defined in (8). Substituting (16) into Eq.(1), we arrive at the ordinary differential equation

$$iV_x\chi +\frac{1}{2}_x^2\chi +|\chi |^2\chi +i\alpha _1|\chi |^2_x\chi +i\alpha _2\chi _x|\chi |^2$$

$$+i\alpha _3_x^3\chi =\mathrm{\Lambda }\chi .$$
$`(17)`$
Imposing periodic boundary conditions at the ends of a sufficiently broad interval and considering V and $`\chi (0)`$ as given parameters and $`\mathrm{\Lambda }`$ as the eigenvalue, we arrive at a nonlinear eigenvalue problem which can be numerically solved by a kind of shooting method. An example of such solution at $`\alpha _2=0`$ and

$$\alpha _1=0.3,\alpha _3=0.1,\chi (0)=1.5,V=0$$
$`(18)`$
is shown in Fig.1(a,b) . It is a nonlocal steady pulse with $`|\mathrm{\Psi }(x)|_{max}=|\mathrm{\Psi }(x_0)|`$ and small symmetric ”wings”. At large $`|xx_0|`$ we can linearize Eq. (17); then we see that

$$\chi (x)\text{const}e^{i\kappa x}(|xx_0|1),$$
$`(19)`$
where $`\kappa `$ is a real root of cubic algebraic equation

$$\mathrm{\Lambda }\kappa V+(1/2)\kappa ^2\alpha _3\kappa ^3=0.$$
$`(20)`$
Finding $`\mathrm{\Lambda }`$ from the numerical solution of the eigenvalue poblem and $`\kappa `$ from (20), we have for the case (18)

$$\mathrm{\Lambda }=0.979,\kappa =5.34.$$
$`(21)`$
On the other hand we can determine $`\kappa `$ directly from $`\chi (x)`$ which is found numerically from Eq. (17), together with $`\mathrm{\Lambda }`$. From Fig. 1b we see that the numerically found $`\chi (x)`$ has indeed asymptotic behavior (19) with $`\kappa 5.3`$ which is very close to the root of Eq. (20), written in (21). This agreement is an evidence of the correctness of the numerical solution of eigenvalue problem by the shooting method. In a similar way, the solution of eigenvalue problem (17) and Eq. (20) for $`\alpha _2=0`$ and

$$\alpha _1=0.8,\alpha _3=0.1,\chi (0)=1.5,V=0.225$$
$`(22a)`$
gives

$$\mathrm{\Lambda }=1.209,\kappa =5.03.$$
$`(22b)`$
The asymptotic behavior of numerically found $`\chi (x)`$ gives the same $`\kappa 5.03`$.

For the other case with $`\alpha _2=0`$ and

$$\alpha _1=1.1,\alpha _3=0.1,\chi (0)=1.5,V=0.35$$
$`(23a)`$
we have

$$\mathrm{\Lambda }=1.937,\kappa =6.12.$$
$`(23b)`$
while from the numerically found $`\chi (x)`$ we obtain $`\kappa =6.1`$.

In fact, similar quasisoliton solutions, looking like embedded solitons, were obtained, by different approaches, for other highly dispersive systems as well \[15-17,9,11,12,17-20\].

Now consider the embedded solitons, starting from Galilean transformation (5). Writing it in the form

$$\mathrm{\Psi }(X,T)=\psi (x,t)e^{iKx}\text{exp}[i(KV\mathrm{\Omega })t]$$
$`(24)`$
and assuming that

$$\psi (x,t)=\stackrel{~}{\psi }(x)\text{exp}\left(i\frac{1}{2}\lambda ^2t\right),$$
$`(25)`$
where $`\lambda `$ is a constant parameter, we compare (24) with (16). Then we have

$$\stackrel{~}{\psi }(x)=\chi (x)\text{exp}(iKx),$$
$`(26)`$
$$\lambda ^2=2(\mathrm{\Lambda }+\mathrm{\Omega }KV).$$
$`(27)`$
From (25) and (9) we arrive at the following equation for $`\stackrel{~}{\psi (x)}`$

$$\frac{1}{2}\lambda ^2\stackrel{~}{\psi }+\frac{1}{2}a_2_x^2\stackrel{~}{\psi }+q|\stackrel{~}{\psi }|^2\stackrel{~}{\psi }+i\alpha _1|\stackrel{~}{\psi }|^2_x\stackrel{~}{\psi }+i\alpha _2\stackrel{~}{\psi }_x|\stackrel{~}{\psi }|^2$$

$$+i\alpha _3_x^3\stackrel{~}{\psi }=0,$$
$`(28)`$
Substituting here (26), (27) and taking into account (6), (7) and (10), (11) we have Eq.(17) as one should expect. From Eq. (26) it is seen that $`|\stackrel{~}{\psi }|^21`$ at $`|x|1`$. Linearizing Eq. (28) we obtain

$$\stackrel{~}{\psi }e^{ikx}(|x|1),$$
$`(29)`$
where $`k`$ is a root of the equation

$$2\alpha _3k^3a_2k^2\lambda ^2=0.$$
$`(30)`$
Substituting (19) and (29) into Eq. (26), we see that

$$\kappa =k+K.$$
$`(31)`$
Then, using (6), (7) and (10), (27) we can easily prove that Eqs. (30) and (20) are equivalent, as it should be.

At small $`\alpha _3`$, the solution of Eq.(28) can be written as

$$\stackrel{~}{\psi }(x)=[u_s(x)+f(x)]\text{exp}[i\varphi _s(x)].$$
$`(32)`$
with small $`f(x)`$. Here $`u_s(x)`$ and $`\varphi _s(x)`$ are defined by the requirement that

$$F(x)=u_s(x)\text{exp}[i\varphi _s(x)]$$
$`(33)`$
is a soliton solution of Eq. (28) without the last term i.e.

$$\frac{1}{2}\lambda ^2F+\frac{1}{2}a_2_x^2F+q|F|^2F+i\alpha _1|F|^2_xF+i\alpha _2F_x|F|^2=0$$
$`(34)`$
and $`F(x)0`$ at $`x=\pm \mathrm{}`$. Solving Eq. (34), we have

$$u_s(x)=\lambda \sqrt{\frac{2p}{q}}\left[\text{cosh}\left(\frac{2\lambda }{\sqrt{a_2}}x\right)+p\right]^{1/2},$$
$`(35)`$
$$\varphi _s(x)=\frac{\alpha _1+2\alpha _2}{2A}\text{arctan}\left[\sqrt{\frac{1p}{1+p}}\text{tanh}\left(\frac{\lambda x}{\sqrt{a_2}}\right)\right],$$
$`(36)`$
$$p=\frac{\sqrt{a_2}q}{\sqrt{4A^2\lambda ^2+a_2q^2}},$$
$`(37)`$
$$A^2=\frac{4\alpha _1(\alpha _1+2\alpha _2)(\alpha _1+2\alpha _2)^2}{12}.$$
$`(38)`$
Eqs. (35)-(38) were obtained in Ref.  for a pulse part of RRS. From (35) it follows that the soliton amplitude is

$$u_0=\lambda \sqrt{\frac{2p}{(1+p)q}}$$
$`(39)`$
and its width is given by

$$\delta =\sqrt{a_2}/\lambda .$$
$`(40)`$
Thus one must require $`a_2>0`$ and from Eq. (10) it follows that at $`\alpha _3K>0`$,

$$|K|<\frac{1}{6}|\alpha _3|,$$
$`(41)`$
which is a restriction on the soliton velocity (6).

A small term $`f(x)`$ in Eq.(32) expresses the effect of the third order dispersion. At $`|x|\sqrt{(}a_2)/\lambda `$, or less, it describes a modification of the pulse, arising due to the last term in Eq. (28), while at large $`x`$

$$f(x)\stackrel{~}{\psi }(x)e^{ikx}$$
$`(42)`$
where $`k`$ is a root of Eq. (30). Note, that this equation was derived in Ref.  for the wavenumber of the resonantly generated radiation by the RRS in the reference frame where the RRS is at rest. From all that we conclude that (32) coincides with the asymptotic expression for the RRS and its radiation at large t.

Let us now compare the embedded solitons, obtained above numerically, with the solution (32) at condition (13) . From Eq.(38) we have

$$A=\frac{1}{2}|\alpha _1|.$$
$`(43)`$
Eqs. (6) and (10) give

$$K=\frac{1\sqrt{112\alpha _3V}}{6\alpha _3},$$
$`(44)`$
$$a_2=\sqrt{112\alpha _3V}.$$
$`(45)`$
Using (39) and (37), we obtain

$$u_0^2=\frac{\sqrt{a_2}}{2A^2}(\sqrt{4A^2\lambda ^2+q^2a_2}q\sqrt{a_2}),$$
$`(46)`$
where q is defined in (11). From (27) and (6),(7) we find

$$\lambda ^2=2\left(\mathrm{\Lambda }\frac{1}{2}K^2+2\alpha _3K^3\right).$$
$`(47)`$
Then we have for the case (18), (21), which is shown in Fig. 1: $`u_01.37`$. For the case (22) $`u_01.41`$, and for (23) $`u_01.42`$. So, in all three cases we approximately have $`u_01.4`$ which is rather close to $`\chi (0)=1.5`$, assumed in the nonlinear eigenvalue problem for all three cases. This supports the conclusion that expressions (35) and (36) approximately describe the pulse in the embedded soliton (at small $`\alpha _3`$ ).

Now consider $`f(x)`$. At small $`\alpha _3`$, the roots of Eq. (30) have simple analytical expressions . Neglecting the first term in (30), we have two smallest roots

$$k\pm i\frac{\lambda }{\sqrt{a_2}}.$$
$`(48)`$
Substituting this into (42), we have

$$\stackrel{~}{\psi }f(x)\text{exp}(\lambda /\sqrt{a_2}x).$$
$`(49)`$
This is in agreement with the asymptotic behavior of expression (33) at large $`x`$. The third root can be approximately obtained if one neglects the last term in (30). This gives

$$ka_2/2\alpha _3,$$
$`(50)`$
which approximately determines the wavenumber of the plane wave in the wings. Expressions (50) and (31) approximately give the roots of Eq. (20) in analytical form. From (50) it follows

$$\text{sgn}k=\text{sgn}\alpha _3.$$
$`(51)`$
And, finally, from the results of Ref.  it follows that at small $`\alpha _3`$

$$f(x)B\left(\frac{\sqrt{a_2}|k|}{A}\right)^{\frac{1}{2}}\text{exp}\left[\frac{\pi \sqrt{a_2}|k|}{4\lambda }\left(1+\frac{2}{\pi }\text{arcsin}p\right)\right]e^{ikx},$$
$`(52)`$
where B is a complex constant with $`|B|110`$. This expression is valid at

$$\frac{\sqrt{a_2}|k|}{\lambda }1,$$
$`(53)`$
i.e. when the wavenumber $`k`$ is much larger than the inverse width of the soliton (33). In expression (52), it is also assumed that $`2A\lambda \lambda /\sqrt{a_2}|k_{1,2}|`$ (this does not exclude $`A\lambda 1`$). If

$$2A\lambda \lambda /\sqrt{A_2}|k_{1,2}|(\text{or}\mathrm{\hspace{0.33em}\hspace{0.33em}2}A\lambda <\lambda /\sqrt{a_2}|k_{1,2}|),$$
$`(54)`$
which may be satisfied only at $`p1`$ and $`q1`$, we have

$$f(x)=B\sqrt{a_2}|k|\text{exp}\left(\frac{\pi \sqrt{a_2}|k|}{2\lambda }\right)e^{ikx}.$$
$`(55)`$
It is easy to check that we can arrive at Eq. (55) by substituting in Eq. (52) the first of conditions (54) and p=1 .

From the all foregoing, one can see a connection between the embedded and resonantly radiating solitons. In the next section we present numerical experiments disclosing this connection from another side.

IV. THE CUT OFF OPERATION

Let us define the cut off operation transforming the embedded soliton into an isolated pulse. Turning to the function $`\chi (x)`$ in (16), we write

$$\chi _{cut}(x)=\chi (x)r(x),$$
$`(56)`$
where $`r(x)`$ is a cutting factor which we take in the form

$$r(x)=\frac{1}{2}\left[\text{tanh}\left(\frac{xx_0+\mathrm{\Delta }x}{\gamma }\right)\text{tanh}\left(\frac{xx_0\mathrm{\Delta }x}{\gamma }\right)\right].$$
$`(57)`$
Here $`x_0`$ is the center of the pulse and $`\mathrm{\Delta }x>0`$ is the width of the cutted function $`\chi _{cut}(x)`$. According to (57), $`r(x)`$ vanishes at $`|xx_0|\mathrm{}`$ and the positive parameter $`\gamma `$ characterizes the ”sharpness” of vanishing. Assuming that $`\gamma `$ is small enough and $`\mathrm{\Delta }x`$ is such that the factor $`r(x)`$ cuts off only the wings without essential disturbing the pulse in $`\chi (x)`$ ( $`\gamma `$ and $`\mathrm{\Delta }x`$ can be properly chosen in numerical tests) we then take the pulse $`\chi _{cut}(x)`$ as the initial condition to Eq.(1) (see Fig.2a)

$$\mathrm{\Psi }(X,0)=\chi _{cut}(X).$$
$`(58)`$
In Fig.2b one can see that the cutted pulse emits radiation at $`t>0`$. However, on the left hand side of the pulse the radiation spreads out with time; this shows that on the left hand side there is a transient radiation, emitted at small t due to initial condition (58). On the contrary, on the right hand side we see a wave train with approximately constant amplitude, with the front propagating to the right; so, the length of the wavetrain increases with time. Therefore the cutted pulse permanently emits radiation to the right. The spectrum, shown in Fig.1c, has a peak at $`\kappa `$ approximately equal to the root of Eq. (20) with account of finite velocity of the pulse at $`t=138`$ (it is still rather close to $`\kappa `$ from Eq. (21), i.e. to the wavenumber of the wing waves in embedded soliton).

Analyzing the time behavior of $`\text{arg}\mathrm{\Psi }`$, we find that $`d(\text{arg}\mathrm{\Psi })/dt=\mathrm{\Lambda }(t)`$ is a slow function of t, with $`\mathrm{\Lambda }(0)0.98`$ and $`\mathrm{\Lambda }(130)0.85`$. Note that this $`\mathrm{\Lambda }(0)`$ coincides, with a good accuracy, with (21). Similar results were obtained for cutted pulses in cases (22) and (23).

We conclude that the cutted pulses are radiating solitons; and the radiation is permanently emitted only in one direction. As far as $`\kappa `$ is connected with $`k`$ from Eq. (31), which is the real root of Eq. (30), obtained from the resonant condition, we conclude that the cut off operation transforms ES ino RRS. The front of radiated wavetrain propagates with the group velocity $`U(k)`$ given by

$$U(k)=a_2k3\alpha _3k^2a_2/4\alpha _3.$$
$`(59)`$
From this it follows that

$$\text{sgn}U(k)=\text{sgn}k=\text{sgn}\alpha _3.$$
$`(60)`$
Therefore at $`\alpha _3<0`$ it should be $`k<0`$ and $`U(k)>0`$. This means that the soliton in Fig. 2a should permanently emit radiation to the right while the peak in the spectrum of the wavetrain should be at negative $`k`$. This is in agreement with the results presented in Figs. 2b and 2c.

To describe analytically the whole system after the cut off, one can use at small $`\alpha _3`$ the equation

$$\stackrel{~}{\psi }(x,t)=[u_s(x)+\eta (x,t)]\text{exp}[i\varphi _s(x)],$$
$`(61)`$
where $`u_s(x)`$ and $`\varphi _s(x)`$ are given (in adiabatic approximation) by Eqs. (35) and (36) and $`\eta (x,y)`$, at large $`x`$ and $`t`$, has the following asymptotic expression

$$\eta (x,t)f(x)\mathrm{\Theta }(Ux)\mathrm{\Theta }(|U|t|x|).$$
$`(62)`$
Here $`f(x)`$ is given by Eq. (52) and $`\mathrm{\Theta }(Y)`$ is the Heaviside function

$$\mathrm{\Theta }(Y)=1(Y>0),\mathrm{\Theta }(Y)=0(Y<0),.$$
$`(63)`$
Eq. (62) expresses that the soliton radiates in the direction of group velocity $`U`$ and the radiation front propagates with the velocity $`|U|`$. The mentioned above adiabatic approximation means that at sufficiently small $`\alpha _3`$, the radiation is so small that the soliton parameters in $`u_s(x)`$ and $`\psi _s(x)`$ (as well as the wavenumber $`k`$) can be considered as constant. However, the soliton losses may be essential at large times.

The variation of soliton parameters, caused by the radiation, can be estimated by means of the integrals of motion (Appendix B). A vast information about the soliton evolution, caused by the radiation, follows from the numerical solution of Eq.(1). For the parameters used in Figs. 2, the results following from this solution are presented in Fig.3. One finds that the mean value of $`|\mathrm{\Psi }|_{max}=u_0(t)`$, shown in Fig.3a, decreases logarithmically, similar to what was found for the case $`\alpha _1=\alpha _2=0`$ . Analyzing the soliton position (Fig.3b) one can find that the soliton velocity $`V(t)`$ increases from $`V(0)=0`$ to $`V(138)0.29`$, i.e. the soliton is accelerating and the function $`V(t)`$ increases (also logarithmically). Now, using $`\mathrm{\Lambda }(0)0.98`$, $`\mathrm{\Lambda }(138)0.85`$ and Eqs. (47),(44), we have $`\lambda (0)1.4`$, $`\lambda (138)1.28`$. Therefore, the soliton width $`\mathrm{\Delta }=\sqrt{a_2}/\lambda `$ increases from $`\mathrm{\Delta }(0)0.71`$ to $`\mathrm{\Delta }(138)0.99`$, i.e. the radiating soliton is widening with time . All that can be considered as a numerical confirmation of analytical results obtained in Appenix B from the conservation laws, in particular that

$$\text{sgn}\frac{dV}{dt}=\text{sgn}k=\text{sgn}U,$$
$`(64)`$
which means that the soliton is accelerating in the direction of the group velocity of resonant radiation.

V. SPECIAL CASE $`\alpha _1=6\alpha _3`$. HIROTA SOLITONS

The coefficient q in Eq. (9) disappears at

$$K=1/\alpha _1.$$
$`(65)`$
Substituting this in Eq. (10) , we have $`a_2=16\alpha _3/\alpha _1`$. Therefore, at condition (3) the coefficient $`a_2`$ also vanishes and Eq. (9) takes the form

$$_t\psi +6\alpha _3|\psi |^2_x\psi +\alpha _2\psi _x|\psi |^2+\alpha _3_x^3\psi =0.$$
$`(66)`$
which is related to Eq.(1) by a Galilean transformation (5),(8) with

$$V=\frac{1}{12\alpha _3},K=\frac{1}{6\alpha _3},\mathrm{\Omega }=\frac{1}{108\alpha _3^2}.$$
$`(67)`$
Consider a particular solution of Eq.(66)

$$\psi (x,t)=e^{i\theta }\zeta (x,t),$$
$`(68)`$
where $`\theta =\text{const}`$ and $`\theta `$, $`\zeta (x,t)`$ are real. Then Eq. (66) is reduced to the modified Korteweg-de Vries (MKdV) equation

$$_t\zeta +2(3\alpha _3+\alpha _2)\zeta ^2_x\zeta +\alpha _3_x^3\zeta =0.$$
$`(69)`$
It is completely integrable and, in particular, has exact N \- soliton solutions  if $`\alpha _23\alpha _3=(1/2)\alpha _1`$ \[cf. with Eq. (2b)\]. As far as Eq. (66) is reduced to Eq. (69) only for particular initial conditions, it may not be , generally, integrable; however, one can assert that it also possess N \- soliton solutions at any $`\alpha _2(1/2)\alpha _1`$, which is a generalization of the Hirota solution at $`\alpha _2=0`$. Below, we call Eq. (66) complex MKdV equation. The Painlevé test , applied to Eq. (66), shows that it has the Painlevé property if $`\alpha _2=3\alpha _3`$, i.e.

$$6\alpha _3=\alpha _1=2\alpha _2.$$
$`(70)`$
This is just the integrability condition for Eq.(1) found by Sasa and Satsuma . These conclusions are in agreement with those following from the Painlevé test of Eq.(1) .

If, in addition to condition (3), $`\alpha _2=0`$ (Hirota case ), we arrive at the complex MKdV equation

$$_t\psi +6\alpha _3|\psi |^2_x\psi +\alpha _3_x^3\psi =0.$$
$`(71)`$
The soliton solution to this equation is

$$\psi _s(x,t)=a\text{sech}[a(xct)]\text{exp}[i(px\sigma t)]$$
$`(72)`$
with arbitrary p and

$$c=3\alpha _3p^2+\alpha _3a^2,\sigma =\alpha _3p^3+3\alpha _3pa^2.$$
$`(73)`$
Then

$$\mathrm{\Psi }_s(X,T)=a\text{sech}[a(XV_sT)]\text{exp}[i(\kappa X\omega T)],$$
$`(74)`$
where

$$V_s=V+c,\kappa =K+p,\omega =\sigma +\mathrm{\Omega }+pV.(75)$$

Therefore $`\kappa `$ is arbitrary and $`V_s(\kappa )`$, $`\omega (\kappa )`$ coincide with Eqs. (2c,d) respectively. This means that Eqs. (74) and (75) indeed describe the Hirota solitons, mentioned in Sec.2a. It is also seen that though the Hirota conditions do not ensure the integrability ( Painlevé tests are negative !), Eq. (1) in Hirota case definitely has exact N -soluton solutions.

Let us now consider numerical solutions of the initial value problem to Eq.(1) in Hirota case

$$\alpha _1=6\alpha _3=0.6,\alpha _2=0.$$
$`(76)`$
As initial condition, we choose a complex pulse

$$\mathrm{\Psi }(X,0)=A\text{sech}\left[\frac{1}{2}(XX_0)\right]\text{exp}[iC(XX_0)]$$
$`(77)`$
with A=1 , C=0.5 and $`X_0=85.3`$. The structures of $`\mathrm{\Psi }(X,T)`$ at T=0 and T=168 are shown in Figs. 4(a,b). We see that $`\mathrm{\Psi }(X,T)`$ evolves into two pulses, well separated at T=168 . The spetrum of $`\mathrm{\Psi }(X,T)`$ at T=168 (Fig.5) gives an evidence that, in addition to these pulses, no significant wavetrains are produced at initial condition (77). The time variation of the amplitude of the pulse with largest amplitude, $`|\mathrm{\Psi }^{(1)}|`$, is shown in Fig.6. It approaches the constant limit, approximately equal to 1.5. The velocity of that pulse approaches the limit approximately equal to 0.35 . These results suggest that the leading pulse asymptotically approaches a regular soliton with $`a^{(1)}1.5`$ and $`V_s^{(1)}0.35`$. The shapes of $`\mathrm{\Psi }(X,T)`$, obtained numerically for $`T>150`$, confirm this conclusion. Then, using Eqs. (2c,d) and (4), we can calculate the shifted frequency $`\omega ^{}=\omega (k)kV_s^{(1)}`$, which is the soliton frequency in the reference frame where the soliton is at rest. The plot of $`\omega ^{}`$ versus $`V_s`$, shown in Fig.7, gives $`\omega ^{}(0.35)1.5`$. This must be equal to $`\text{arg}\mathrm{\Psi }^{(1)}`$. On the other hand, the numerical solution permits to obtain $`\text{arg}\mathrm{\Psi }^{(1)}`$ versus time in the frame where the pulse is at rest. This also gives $`\text{arg}\mathrm{\Psi }^{(1)}1.5`$, at sufficiently large T . Similar results were obtained for the second pulse. Thus we conclude that the initial complex pulse (77) decays into two Hirota solitons. Qualitatively similar results were obtained for the initial pulse

$$\mathrm{\Psi }(X,0)=\text{exp}\left[\frac{(XX_0)^2}{16}+i\frac{XX_0}{2}\right],X_0=85.3.$$
$`(78)`$
Other features of the dynamics of Hirota solitons, in particular their collisions , will be examined in another paper.

VI. PULSE EVOLUTION AT $`\alpha _16\alpha _3`$, $`\alpha _2=0`$

Here, we report on the numerical solutions of Eq.(1) in general case. First we take, as initial condition, the pulse (77) with A=1.9, C=0 and $`X_0=100`$ and assume that

$$\alpha _1=0.5,\alpha _3=0.1.$$
$`(79)`$
The behavior of solution at these parameters is shown in Fig.8, where the case $`\alpha _3=0`$ (with $`\alpha _1=0.5`$) is also presented for comparison. In the latter case, the initial pulse splits into four solitons of the form (33), (35)-(38), with $`a_2=1`$ and moving to the left (Fig.8a). On the other hand, in the case (79) the initial pulse splits into three radiating solitons. Two of them propagate to the right and the smallest one - to the left.

The radiating wavetrains propagate to the right with the group velocities larger than the soliton velocities (Figs. 8b and 8c); this is in agreement with Eq. (60). The time dependence of the soliton coordinates is shown in Fig.9 . At $`\alpha _3`$, the soliton velocities are constant. At $`\alpha _3=0.1`$, the solitons are accelerating in positive direction, which agrees with Eq. (64).

The Fourier spectra for both cases, $`\alpha _3=0`$ and $`\alpha _3=0.1`$, are shown in Fig. 10. The spectrum at $`\alpha _3=0`$ exhibits no resonant radiation (Fig.10a) while at $`\alpha _3=0.1`$ we see three distinct resonant peaks with negative $`\kappa `$, which are the wave numbers of the resonant radiation emitted by the solitons in the ”laboratory” frame (Fig.10b). The difference between the wave numbers follows from Eq. (20) which is valid, as we have seen, both for embedded and radiating solitons. As far as the solitons have different velocities and amplitudes (the latters are determined by the soliton parameter $`\mathrm{\Lambda }`$), they have different wavenumbers because $`\kappa `$, being a root of Eq. (20), depends on the soliton velocity and amplitude.

An important difference between the two cases, $`\alpha _3=0`$ and $`\alpha _30`$, is seen in Figs. 11 and 12. In the first case, the average soliton amplitudes have constant limits at $`t\mathrm{}`$ while in the second case they must slowly decrease because of the soliton radiation.

In the aforegoing, we studied the case when $`\alpha _1\alpha _30`$. Now let us consider $`\alpha _1\alpha _30`$. We take

$$\alpha _1=0.6,\alpha _3=0.1,$$
$`(80)`$
and initial condition (77) with $`C=1/2`$, $`X_0=850`$. The solution at $`T=370`$ is shown in Fig. 13. We see a pulse with radiation on both its sides. An analysis shows that on the left hand side there is the resonantly generated radiation (in agreement with Eq. (60), where now $`\alpha _3>0`$, $`k>0`$, $`U<0`$). On the right hand side, there is a radiation composed from the harmonics of continuous spectrum; it satisfies the linearized Eq.(1) and has been emitted in the transient period of time. The spectal distribution is shown in Fig.14. One can see two spectral maximums at positive wavenumbers. The narrow one corresponds to the resonant radiation, while the broader peak is composed of the continuous spectrum. Its structure can be understood from the dispersion equation

$$\omega =(1/2)k^2\alpha _3k^3$$
$`(81a)`$
and the correspondig expression for the group velocity of continuous spectrum

$$V_g(k)=k3\alpha _3k^2.$$
$`(81b)`$
From this we see that, at $`\alpha _3>0`$, the continuos radiation on the right hand side of the pulse is composed from

$$0<kk_0=1/6\alpha _3,$$
$`(81c)`$
where $`V_g(k_0)=\text{max}V_g(k)`$. As far as $`k_0`$ is less then the wavenumber of resonant radiation, approximately given by Eq.(50) (note that in our case , according to (45), $`a_2>1`$), one can see why the peak of continuous spectra is less than $`\kappa `$ of resonant radiation. All this shows that the pulse in Fig. 13 is nothing than RRS. The time behavior of its amplitude is shown in Fig. 15 ; qualitatively, it is similar to Fig. 3a.

It is reasonable to compare these results with the case $`\alpha _1=0.6`$, $`\alpha _3=0`$ , at the same initial condition. Then we have two regular solitons, propagating to the right. The soliton positions versus time for both cases are shown in Fig. 16. Note a very small acceleration of RRS in the direction of group velocity (which is now negative) .

VII. CONCLUSIONS

We have considered soliton and quasisoliton solutions of Eq. (1) and their relationship. Solutions (2), describing regular solitons, degenerate at $`\alpha _30`$; this shows that they exist due a balance between nonlinear terms and linear third order dispersion. On the other hand the quasisolitons (embedded and resonantly radiating)

turn at $`\alpha _30`$ into the regular soliton solutions of Eq. (1) without the third derivative term. Apart from that, at $`\alpha _2=0`$, the regular solitons exist only in the special case $`\alpha _1=6\alpha _3`$.

The resonantly radiating solitons (RRS) are nonsteady; the amplitudes of their radiation at small $`\alpha _3`$ are exponentially small and so their parameters change logarithmically slow in this case. Only such quasisteady solitons have sufficiently large lifetimes to be of physical significance. The embedded solitons (ES), which are steady structures consisting of pulses embedded into plane waves (wings), have close connection with RRS, similar to other systems. The corresponding analytical treatment in Sec.III is suplemented by numerical simulations with the cut off operation (Sec.IV) ; it is demonstrated that it leads to the transformation of ES into RRS that radiates in the direction of the radiation group velocity . Just after the cutt off, the amplitude, velocity and wavenumber of the radiated wavetrain coincide with those of wings; marching in time, we have seen that the amplitude decreases, the soliton velocity increases and the wavenumber changes according to Eqs.(30) and (31). This is in agreement with conservation laws (Appendix B).

In Sec. V we have investigated a special case $`\alpha _1=6\alpha _3`$, interesting from a theoretical point of view . Then the quasisolitons, ES and RRS, do not exist and Eq. (1) can be reduced to the complex MKdV equation which has N -soliton solutions (consisting of regular solitons). This generalizes the Hirota solutions obtained for $`\alpha _1=6\alpha _3`$, $`\alpha _2=0`$. The complex MKdV equation possesses the Painlevé property only at $`\alpha _2=3\alpha _3`$, which is the case when Eq.(1) is integrable.

Presumably, in the nonlinear processes only regular solitons and RRS may be of the physical significance. From Eq. (2b) it follows, however, that the amplitude of regular soliton vanishes at $`\alpha _30`$, while for the RRS it remains finite. Therefore one should expect that at small $`\alpha _3`$, the RRS are more important in the nonlinear evolution. This was confirmed by solving the initial value problem (Sec. VI), where the RRS and not regular solitons are seen emerging from the initial disturbances. (Small $`\alpha _3`$ is the most interesting case from physical point of view because the third derivative term is in fact the result of an expansion) .

This work was partly supported from INTAS grant No 99-1068.

Appendix A

Modulational instability of plane wave, according to Eq.(1)

Eq. (1) has exact plain wave solution

$$\mathrm{\Psi }=A\text{exp}(i\kappa Xi\omega T).$$
$`(A1)`$
Substituting (A1) into (1), we have

$$\omega =\frac{1}{2}\kappa ^2(1\alpha _1\kappa )A^2\alpha _3\kappa ^3.$$
$`(A2)`$
Substituting a slightly perturbed wave (A1)

$$\mathrm{\Psi }=A(1\chi )\text{exp}(i\kappa Xi\omega T),$$
$`(A3)`$
into Eq.(1), we have

$$i\chi _T+\frac{1}{2}(16\alpha _3K)\chi _{XX}+i(K+\alpha _1A^23\alpha _3K^2)\chi _X$$

$$+(1\alpha _1\kappa )(\chi +\chi ^{})A^2+i\alpha _2(\chi _X+\chi _X^{})A^2.$$
$`(A4)`$
Writing $`\chi =u+iw`$ and assuming that

$$(u,w)\text{exp}(ipXirT).$$

we obtain dispersion equation, which is convenient to write in the form

$$\mathrm{\Gamma }^22\alpha _2kA^2\mathrm{\Gamma }\frac{1}{4}(16\alpha _3K)k^2[(16\alpha _3K)k^2$$

$$4(1\alpha _1K)A^2]=0,$$
$`(A5)`$
where

$$\mathrm{\Gamma }=r+\alpha _3p^3(\kappa 3\alpha _3\kappa ^2+\alpha _1A^2)p.$$
$`(A6)`$
Thus, at real p , $`\text{Im}\gamma =\text{Im}r`$ and the stability condition is

$$(16\alpha _3\kappa )^2p^24A^2[(1\alpha _1\kappa )(16\alpha _3\kappa )\alpha _2^2A^2].$$
$`(A7)`$
Unlike the plain wave solution of the regular NLS, now the plain wave can be stable at any p , if

$$(1\alpha _1\kappa )(16\alpha _3\kappa )\alpha _2^2A^2.$$
$`(A8)`$
At $`\alpha _2=0`$ this is possible at

$$\frac{1}{\alpha _1}\kappa \frac{1}{6\alpha _3}(6\alpha _3<\alpha _1)$$
$`(A9a)`$
$$\frac{1}{6\alpha _3}\kappa \frac{1}{\alpha _1}(\alpha _1<6\alpha _3)$$
$`(A9b)`$
Applying the stability criterion to the wings of embedded solitons, where $`\kappa 1/2\alpha _3`$, we see that this may satisfy condition (A9b) ; then one can expect that the wings are stable. We should also take into account only $`p>L^1`$, where L is the period in the numerical scheme; this relaxes the limitatations following from the stability criterion.

Appendix B

Investigation of the soliton evolution by means of the integrals of motion.

First , we present a general analysis of conserving integrals N, P, H. Substituting (5) into expressions (12), (14), (15) we have after simple algebra

$$N=_{\mathrm{}}^{\mathrm{}}|\psi (x,t)|^2𝑑x,$$
$`(B1)`$
$$P=KN+\frac{1}{2}i_{\mathrm{}}^{\mathrm{}}(\psi _x\psi ^{}\psi ^{}_x\psi )𝑑x.$$
$`(B2)`$
where $`\psi (x,t)`$ satisfies Eq.(9).

To have a convenient expression for H, we substitute in Eq. (15) $`\alpha _3_x^3\mathrm{\Psi }`$ and $`\alpha _3_x^3\mathrm{\Psi }^{}`$ from Eq.(1), to obtain

$$H=\frac{i}{2}_{\mathrm{}}^{\mathrm{}}𝑑X(\mathrm{\Psi }^{}_T\mathrm{\Psi }\mathrm{\Psi }_T\mathrm{\Psi }^{})+\frac{1}{2}_{\mathrm{}}^{\mathrm{}}𝑑X|\mathrm{\Psi }|^4$$

$$+\frac{i}{4}\alpha _1_{\mathrm{}}^{\mathrm{}}𝑑X|\mathrm{\Psi }|^2(\mathrm{\Psi }^{}_X\mathrm{\Psi }\mathrm{\Psi }_X\mathrm{\Psi }^{}).$$

Then using the Galilei transformation in the form (24) and taking into account that at large $`t`$

$$_ti(\lambda ^2/2)\psi ,$$

we arrive at the following asymptotic expression

$$H\left(\mathrm{\Omega }KV\frac{1}{2}\lambda ^2\right)N+PV+\frac{1}{2}q_{\mathrm{}}^{\mathrm{}}|\psi |^4𝑑x$$

$$\frac{1}{2}i\alpha _1_{\mathrm{}}^{\mathrm{}}|\psi |^2(\psi _x\psi ^{}\psi ^{}_x\psi )𝑑x.$$
$`(B3)`$
Substituting (25) and (32) into (B1), we have at large t (cf. Refs. ),

$$NN_s+N_r,$$
$`(B4)`$
where

$$N_s=_{\mathrm{}}^{\mathrm{}}|u_s(x,t)|^2𝑑x,N_r=_{\mathrm{}}^{\mathrm{}}|\eta (x,t)|^2𝑑x$$
$`(B5)`$
are contributions from the soliton and the radiation. Using (35) we have

$$N_s=\frac{2a_2}{\alpha _1}\text{arctan}\left(\frac{\alpha _1\lambda }{\sqrt{a_2}q}\right)$$
$`(B6)`$
and, at large t,

$$N_r|f|^2|U|_t$$
$`(B7)`$
where f is given by (52) or (55). Evidently, $`|f|^2`$ does not depend on x. Then from the conservation of N it follows

$$\frac{d}{dt}\left[\frac{2a_2}{\alpha _1}\text{arctan}\left(\frac{\alpha _1\lambda }{\sqrt{a_2}q}\right)\right]|f|^2|U|.$$
$`(B8)`$
In a similar way, from (B2) and the conservation of P we have the following asymptotic equation at large t

$$N\frac{dK}{dt}+\frac{d}{dt}_{\mathrm{}}^{\mathrm{}}u_s^2(x)_x\psi _s(x)dx+k|f|^2|U|0.$$
$`(B9)`$
From (36) at $`\alpha _2=0`$, we have

$$_x\psi _s(x)=\frac{\alpha _1}{2a_2}u_s^2(x).$$
$`(B10)`$
Taking into account that

$$_{\mathrm{}}^{\mathrm{}}u_s^4(x)𝑑x=\frac{2a_2q}{\alpha _1}\left(\frac{2\sqrt{a_2}\lambda }{q}N_s\right),$$
$`(B11)`$
we transform Eq. (B9) to the form

$$N\frac{dK}{dt}\frac{d}{dt}\left[\frac{q}{\alpha _1}\left(\frac{2\sqrt{a_2}\lambda }{q}N_s\right)\right]+k|f|^2|U|0.$$
$`(B12)`$
Now we turn to Eq.(B3). After differentiation over t and simple calculations we obtain

$$P\frac{dV}{dt}\left(\lambda \frac{d\lambda }{dt}+a_2K\frac{dK}{dt}\right)N$$

$$=\frac{d}{dt}\left[\frac{1}{2}q_{\mathrm{}}^{\mathrm{}}|\psi |^4𝑑x+\frac{i}{4}\alpha _1_{\mathrm{}}^{\mathrm{}}|\psi |^2(\psi ^{}_x\psi \psi _x\psi ^{})𝑑x\right]$$

$$=\frac{d}{dt}[\frac{1}{2}q_{\mathrm{}}^{\mathrm{}}|\psi _s|^4dx+\frac{i}{4}\alpha _1$$

$$_{\mathrm{}}^{\mathrm{}}|\psi _s|^2(\psi _s^{}_x\psi _s\psi _s_x\psi _s^{})dx].$$
$`(B13)`$
Using

$$\frac{dV}{dt}=a_2\frac{dK}{dt},$$
$`(B14)`$
we obtain

$$a_2(PKN)\frac{dK}{dt}\lambda \frac{d\lambda }{dt}N$$

$$\frac{d}{dt}\left(\frac{1}{2}q_{\mathrm{}}^{\mathrm{}}u_s^4𝑑x+\frac{\alpha _1^2}{2a_2}_{\mathrm{}}^{\mathrm{}}u_s^6𝑑x\right).$$
$`(B15)`$
Let us apply these relations to small $`\alpha _1`$ (which is just the case shown in Figs.3). Assuming that

$$\frac{\alpha _1\lambda }{\sqrt{a_2}q}1,$$
$`(B16)`$
we have

$$N_s\frac{2\sqrt{a_2}\lambda }{q}\left(1\frac{\alpha _1^2\lambda ^2}{3a_2q^2}\right).$$
$`(B17)`$
Therefore Eq. (B8) is reduced to

$$\frac{d}{dt}\left[\frac{2\sqrt{a_2}\lambda }{q}\left(1\frac{\alpha _1^2\lambda ^2}{3a_2q^2}\right)\right]=|f|^2|U|.$$
$`(B18)`$
In the same approximation, Eq.(12) takes the form

$$N\frac{dK}{dt}k|f|^2|U|+\frac{2\alpha _1}{3}\frac{d}{dt}\left(\frac{\lambda ^3}{\sqrt{a_2}q^2}\right),$$
$`(B19)`$
At $`\alpha _10`$, Eqs. (B18) and (B19) turn into equations obtained in Ref.. From (B19), (B14) and Eq. (60) it follows an important qualitative result (64).

Now, consider Eq. (B15). Taking into account (B2),(B10), (B11) and (B19) we see that at

$$t\lambda k^2|f|^4U^2$$
$`(B20)`$
and small $`\alpha _1`$, Eq. (B15) is approximately equivalent to (B18). Performing calculations similar to those in Ref.  and using (B14), we conclude that the soliton width $`\lambda ^1`$ and velocity $`V(t)`$ are changing logarithmically slow. The same holds for the soliton amplitude $`u_(9t)`$, determined from (39) and (37), in agreement with numerical solution of Eq.(1) (see Figs.3). The chracteristic time of the variations of soliton parameters has the same order of magnitude as in Ref., which is much less than the right hand side of (B20).

If $`\alpha _1`$ is not small, the analysis is rather tedious and will not be considered here.

References

 A. Hasegawa and Y. Kodama, Solitons in Optical Communications (Oxford University Press, 1995).

 G.P. Agrawal, Nonlinear Fiber Optics (Academic Press, New York, 1995).

 M.J. Potasek and M.Tabor, Phys.Lett. A 154 , 449 (1991).

 R.Hirota, J. Math.Phys. 14 ,805 (1973).

 N.Sasa and J.Satsuma, J. Phys.Soc. Jpn. 60 ,409 (1991).

 E.M. Gromov, L.V. Piskunova and V.V.Tutin , Phys. Lett. A 256 , 153 (1999) .

 J.Yang, B.A.Malomed and D.J. Kaup, Phys. Rev. Lett. 83, 1958 (1999).

 P.K.A.Wai, C.R.Menyuk, Y.C.Lee and Chen, Opt.Lett.11, 464 (1986).

 P.K.A.Wai, H.H.Chen and Y.C.Lee, Phys.Rev. A 41, 426 (1990).

 H.H. Kuehl and C.Y. Zhang, Phys.Fluids B 2, 889 (1990).

 V.I. Karpman, Phys. Rev.E , 47, 2073 (1993) .

 J.P. Boyd,Weakly Nonlocal Solitary Waves and Beyond-All\_Order Asymptotics (Kluver, Dodrecht,1998).

 V.I. Karpman, In : Nonlinearity, Integrability and all that:Twenty years after NEEDS’79. Editors M.Boiti, L.Martina, F.Pempinelli, B. Prinari & G. Soliani, (”World Scientific”, 2000), p.396.

 V.I. Karpman, Phys.Rev. E 62 , 5678 (2000).

 P.Deeskow, H.Schamel, N.N.Rao, M.Y. Yu, R.K. Varma and P.K.Shukla, Phys. Fluids 30 , 2703 (1987).

 V.I. Karpman and H. Schamel, Phys. Plasmas 4,120 (1997); ibid.4 (1997) 2778 (E) .

 Y. Pomeau, A.Ramani and B. Grammaticos, Physica D 31 (1988) 127 .

 V.I.Karpman, Phys. Lett.A 181 (1993) 211.

 V.I. Karpman, Phys.Lett.A 186 (1994) 300.

 V.I.Karpman, Phys.Rev.Lett. 74 (1995) 2455.

 J. Weiss, M.Tabnd G. Carnevale, J. Math. Phys. 24 , 522 (1983) .

 M.A. Ablowitz and H. Segur Solitons and the Inverse Scattering Transform (SIAM, Philadelphia, 1981)

 M. Gedalin, T.C. Scott and Y.B. Band, Phys. Rev. Lett. 78, 448 (1997).

 V. E. Zakharov and E.A. Kuznetsov, JETP 86,1036 (1998).

 V.I. Karpman and A.G. Shagalov, Phys. Lett. A 254 , 319 (1999).

 V.I.Karpman, Phys. Lett. A 244 , 397 (1998) .

 Zhonghao Li, Lu Li, Huipping Tian and Guosheng Zhou, Phys.Rev.

Lett. 84 ,4096 (2000).

 K. Portezian and K. Nakkeeran, Phys.Rev.Lett. 76, 3955 (1996)