text:
# Stochastic dynamics of electronic wave packets in fluctuating laser fields

## I Introduction

The advancement of sophisticated trapping techniques and the development of powerful laser sources has stimulated numerous theoretical and experimental investigations on the dynamics of wave packets in elementary quantum systems. An understanding of their dynamics is important for our conception of quantum mechanics and its relation to classical mechanics. So far most of the research in this context has concentrated on coherent aspects of wave packet dynamics which may be traced back theoretically to semiclassical aspects originating from the smallness of the de Broglie wave lengths involved. However, for an understanding of the emergence of classical behaviour also a detailed understanding of the destruction of quantum coherence is required. Such a destruction of coherence may arise from the coupling of a quantum system to a reservoir or from stochastic external influences. Though by now many aspects of the coherent dynamics of wave packets are well understood still many questions concerning the influence of stochastic perturbations on elementary quantum systems with a high level density are open.

A paradigm of a quantum system in which many of these latter aspects can be investigated in great detail are Rydberg systems interacting with fluctuating laser fields. Due to the inherent stochastic nature of laser light a detailed understanding of optical excitation processes with fluctuating laser fields is of vital interest for laser spectroscopy. Rydberg systems  are of particular interest in this context due to their high level density close to an ionization threshold. Any laser-induced excitation process which involves Rydberg and continuum states close to an ionization threshold typically leads to the preparation of a spatially localized electronic Rydberg wave packet . Under the influence of a fluctuating laser field the coherence of such an electronic wave packet is disturbed whenever it is close to the ionic core where the electron-laser interaction is localized . Eventually these random perturbations are expected to lead to a stochastically dominated Brownian motion of the electronic wave packet. Indeed it has been demonstrated recently  that such a transition to a diffusive behaviour takes place and that this diffusive dynamics are dominated by characteristic power laws which govern the time evolution of the Rydberg system. However, these previous investigations were restricted to a particular type of phase fluctuations of laser fields which can be described by the so called phase diffusion model (PDM) . This PDM implies a Lorentzian spectrum for the fluctuating laser field. It is known from the dynamics of three level systems that the somewhat unrealistic asymptotic frequency dependence of a Lorentzian laser spectrum may lead to unphysical predictions . However, from our previous investigations it remained open to which extent these characteristic, diffusive long time dynamics of an excited Rydberg electron depend on details of the fluctuations of the exciting laser field. Such a dependence might be expected on intuitive grounds as the diffusion of the excited Rydberg electron in energy space eventually also reaches the far-off resonant regions of the laser spectrum.

In this paper we tackle these open questions by generalizing our previous results to arbitrary types of laser fluctuations. For this purpose we derive rate equations for the relevant density matrix elements of the excited Rydberg electron which are averaged over the fluctuations of the laser field. These rate equations are based on a decorrelation of the relevant electron-field averages. This decorrelation approximation (DCA) is valid as long as the characteristic correlation time of the fluctuating laser field, i.e. its inverse bandwidth, is much smaller than all other relevant intrinsic dynamical time scales. Within this framework it will become apparent that it is the laser spectrum only which determines the time evolution of the excited Rydberg electron. On the basis of this approach it will be demonstrated which aspects of the diffusive long time dynamics of an excited electronic Rydberg wave packet depend on which details of the laser spectrum.

The range of validity of these Pauli-type rate equations is investigated in detail for a special class of phase fluctuations  of the exciting laser field. This special class of phase fluctuations implies non-Lorentzian spectra which for large laser frequencies decrease more rapidly than a Lorentzian. These phase fluctuations might be considered as a realistic model for a single mode laser field which is operated well above the laser threshold. In order to access the range of validity of the Pauli-type rate equations for this special class of laser fluctuations a more general master equation is derived for the averaged density operator of the Rydberg electron. This more general approach is also capable of dealing with all coherent aspects of the laser excitation process.

The paper is organized in the following way: In Sec. II the theoretical models for describing the laser fluctuations and the dynamics of the Rydberg electron are presented. For the sake of simplicity we restrict our subsequent discussion to Rydberg systems which can be described within the framework of a one-channel approximation . Typically this approximation is well satisfied for Alkali atoms. In Sec.III we derive rate equations which describe the dynamics of the excited Rydberg electron averaged over the laser fluctuations. Self consistent validity conditions for the applicability of the decorrelation approximation (DCA) are discussed on which these rate equations are based on. In Sec. IV characteristic aspects of the time evolution predicted by the rate equations of Sec. III are exemplified. The different dynamical long time regimes and their characteristic power laws are discussed in detail. From the resulting analytical expressions for these power laws it is apparent to which extent they depend on details of the laser spectrum. In Sec.V we derive a more sophisticated master equation for the averaged dynamics of the excited Rydberg electron. This master equation is capable of describing also coherent aspects of the dynamics of the excited Rydberg electron but its validity is restricted to a particular class of phase fluctuations only. In Sec. VI solutions of this master equation are compared with the corresponding results of the rate equations. Thus we are able to determine the range of validity of the DCA. In our subsequent discussions we use atomic units ($`m_e=e=\mathrm{}=1`$).

## II Theoretical framework

In this section the theoretical models are introduced with which the fluctuating laser field and the excited Rydberg system are described.

### A The fluctuating laser field

We consider an atomic or molecular Rydberg system which is driven by a laser field with electric field strength

$$𝑬(t)=𝒆\epsilon (t)e^{i\omega t}+\text{c.c..}$$
(1)
We assume that this laser field can be described by a classical stochastic process . The mean frequency of this laser field is denoted $`\omega `$ and $`𝒆`$ is its polarization vector. The fluctuations of this laser field are described by the envelope function $`\epsilon (t)`$ which is assumed to be slowly varying on time scales of the order of $`1/\omega `$. The associated spectrum of this laser field is defined by

$$S(\mathrm{\Omega })=\frac{1}{\pi }\text{Re}_0^{\mathrm{}}𝑑\tau K(\tau )e^{i\mathrm{\Omega }\tau }$$
(2)
with the two-time correlation function of the slowly varying amplitude

$$K(\tau )=\epsilon (t+\tau )\epsilon ^{}(t).$$
(3)
Thereby $`\mathrm{}`$ denotes statistical averaging over the fluctuations of the laser field.

For a single mode laser which is operated well above laser threshold to a good degree of approximation the amplitude of the laser field is stable. In these cases fluctuations of a realistic laser field can be described by a classical electromagnetic field whose phase is fluctuating, i.e.

$$\epsilon (t)=\epsilon _0e^{i\mathrm{\Phi }(t)}.$$
(4)
The fluctuating phase $`\mathrm{\Phi }(t)`$ obeys the (Ito-)stochastic differential equation

$$d\varphi (t)=\varphi (t)\beta dt+\sqrt{2b}\beta dW(t)$$
(5)
with $`\varphi (t)=\dot{\mathrm{\Phi }}(t)`$. In Eq.(5) $`1/\beta `$ determines the correlation time of the stochastic frequency $`\varphi (t)`$ and $`dW(t)`$ is the differential of a real-valued Wiener process with zero mean and unit variance, i.e. $`dW(t)=0`$, $`dW(t)^2=dt`$ . Eqs.(4) and (5) imply the relation

$$K(\tau )=|\epsilon _0|^2\mathrm{exp}\left[b\tau +b/\beta \left(1e^{\tau \beta }\right)\right].$$
(6)
In the limit of large values of $`\beta `$ $`\mathrm{\Phi }(t)`$ itself approaches a real-valued Wiener process, i.e.

$$d\mathrm{\Phi }(t)=\sqrt{2b}dW(t).$$
(7)
This limiting case constitutes the so called phase diffusion model (PDM). It implies a Lorentzian laser spectrum of the form

$$S(\mathrm{\Omega })=|\epsilon _0|^2\frac{1}{\pi }\frac{b}{b^2+\mathrm{\Omega }^2}.$$
(8)
For $`\beta b`$ the parameter $`\beta `$ may be interpreted as a cut-off parameter of the laser spectrum. This becomes apparent by noting that for frequencies $`\mathrm{\Omega }\beta `$ the spectrum is always approximately Lorentzian whereas for large frequencies, i.e. $`\mathrm{\Omega }\beta `$, it tends to zero more rapidly as can be seen in Fig.1. More precisely, for $`\beta b`$ we obtain from Eqs.(6) and (2) the approximate relation

$`S(\mathrm{\Omega })`$ $`=`$ $`{\displaystyle \frac{|\epsilon _0|^2}{\pi }}{\displaystyle \frac{b}{\mathrm{\Omega }^2+b^2}}{\displaystyle \frac{1}{1+\left(\frac{\mathrm{\Omega }}{\beta }\right)^2}}.`$ (9)
### B The interaction Hamiltonian

Let us assume that the considered atom or molecule is prepared initially in an energetically low lying bound state $`|g`$ with energy $`ϵ_g`$. Furthermore, it is excited resonantly by the fluctuating laser field to Rydberg and/or continuum states $`|n`$ with energies $`ϵ_n`$ close to an ionization threshold. In the dipole and rotating wave approximation this excitation process can be described by the Hamiltonian

$$H=\underset{j=g,n}{}|jj|ϵ_j\underset{n}{}\left[\epsilon (t)e^{i\omega t}d_{ng}|ng|+\text{h.c.}\right].$$
(10)
The dipole matrix elements between the initial state $`|g`$ and the excited Rydberg states $`|n`$ are denoted $`d_{ng}=n\left|𝒅\mathbf{}𝒆\right|g`$. It is understood that the sum over the excited states $`|n`$ appearing in Eq.(10) includes also an integration over the adjacent continuum states. The energies of the excited Rydberg states and the energy dependence of the dipole matrix elements entering the Hamiltonian of Eq.(10) can be determined with the help of quantum defect theory (QDT) . In the case of excited Rydberg states which can be described within the framework of the one-channel approximation we find, for example,

$`ϵ_n`$ $`=`$ $`{\displaystyle \frac{1}{2(n\alpha )^2}},`$ (11)
$`d_{ng}`$ $`=`$ $`d_{ϵg}(n\alpha )^{3/2}d_{ϵg}|ϵ_nϵ_{n+1}|^{1/2}.`$ (12)
Thereby $`\alpha `$ denotes the quantum defect of the excited Rydberg series and $`d_{ϵg}`$ is the dipole matrix element between the initial state $`|g`$ and an energy normalized continuum state $`|ϵ`$ with energy $`ϵ0`$. Within the framework of QDT $`\alpha `$ as well as $`d_{ϵg}`$ are approximately energy independent close to the ionization threshold, i.e. for energies $`|ϵ|1`$. A one channel approximation is appropriate for all cases in which excited states of the ionic core are located far away from the excited energy region. Typically this condition is fulfilled for Alkali atoms.

The main problem is to solve the stochastic Schrödinger equation associated with the Hamiltonian of Eq.(10). In general this is a complicated task due to the simultaneous presence of the laser fluctuations and of the threshold effects arising from the infinitely many bound Rydberg states converging to the ionization threshold. By the resulting intricate interplay between laser fluctuations and threshold phenomena it is difficult to apply stochastic simulation methods which typically become unreliable due to numerical inaccuracies in particular for long interaction times.

## III Decorrelation approximation (DCA)

In this section an approximation method is developed for determining the dynamics of the Rydberg system in the fluctuating laser field. This approximation method is based on a decorrelation of atom-field averages and leads to a Pauli-type master equation for the density operator of the Rydberg system which is averaged over the laser fluctuations. In this master equation all coherences between different energy levels have been eliminated adiabatically. This decorrelation approximation (DCA) is valid for arbitrary types of laser fluctuations provided these fluctuations are sufficiently fast (compare with conditions (23), (24) and (26) derived below).

Let us start by determining first of all an approximate equation of motion for the probabilities $`\rho _{nn}(t)=n|\psi (t)\psi (t)|n`$ of observing the excited Rydberg system in one of the Rydberg states $`|n`$. Neglecting coherences $`\rho _{nn^{}}(t)`$ with $`nn^{}`$ we find from the stochastic Schrödinger equation with Hamiltonian (10) the relation

$`\underset{nn}{\overset{\mathbf{}}{\rho }}(t)`$ $`=`$ $`2d_{ng}^2\mathrm{Re}\{{\displaystyle _{t_0}^t}𝑑t^{}\epsilon (t)\epsilon ^{}(t^{})e^{i(\overline{ϵ}ϵ_n)(tt^{})}\}`$ (15)
$`\times \left[\rho _{gg}(t^{})\rho _{nn}(t^{})\right]`$
$`\left[id_{ng}\epsilon (t)e^{i(\overline{ϵ}ϵ_n)(tt_0)}\rho _{gn}(t_0)+\text{h.c.}\right]`$
with $`0t_0t`$ and with $`\rho _{kj}(t)=k|\psi (t)\psi (t)|j\mathrm{exp}[it(ϵ_kϵ_j)]`$. The mean excited energy is denoted $`\overline{ϵ}=ϵ_g+\omega +\delta \omega `$ with the quadratic Stark shift contribution of all other (non-resonant) states. In the subsequent discussion we assume for the sake of simplicity that the intensity dependence of $`\delta \omega `$ does not affect the dynamics of the excited Rydberg system. This is valid either for pure phase fluctuations of the laser field or in the case of arbitrary laser fluctuations for laser bandwidths which are much larger than $`\delta \omega `$. The self-consistency condition for the omission of the coherences $`\rho _{nn^{}}(t)`$ with $`nn^{}`$ will be discussed later (compare with Eq.(24)). Taking the integration interval $`[t_0,t]`$ to be smaller than the characteristic time scale over which $`[\rho _{nn}(t^{})\rho _{gg}(t^{})]`$ changes significantly this latter term can be approximated by its value at time $`t`$. If on the other hand the interval $`[t_0,t]`$ is assumed to be larger than the correlation time of the fluctuating laser field we can replace the lower integration limit $`t_0`$ by $`\mathrm{}`$ in Eq.(15). Thus we obtain

$`\underset{nn}{\overset{\mathbf{}}{\rho }}(t)`$ $`=`$ $`2|d_{ng}|^2\left[\rho _{gg}(t)\rho _{nn}(t)\right]`$ (18)
$`\times \text{Re}{\displaystyle _0^{\mathrm{}}}𝑑\tau \epsilon (t)\epsilon ^{}(t\tau )e^{i\tau (\overline{ϵ}ϵ_n)}`$
$`\left\{id_{ng}\epsilon (t)e^{i(tt_0)(\overline{ϵ}ϵ_n)}\rho _{gn}(t_0)+\text{h.c.}\right\}.`$
Now we are able to carry out the statistic average $`\mathrm{}`$ over the laser fluctuations. Due to the above mentioned conditions on the integration interval the involved density matrix elements and the laser field $`\epsilon (t)`$ decorrelate. As $`\epsilon (t)=0`$ the contribution of the last term on the right hand side of Eq.(18) vanishes. The remaining terms yield the rate equation

$$\underset{nn}{\overset{\mathbf{}}{\rho }}(t)=_{ng}[\rho _{gg}(t)\rho _{nn}(t)]$$
(19)
with the time independent rates

$$_{ng}=2\pi |d_{ng}|^2S(\overline{ϵ}ϵ_n)$$
(20)
and with the laser spectrum $`S(\mathrm{\Omega })`$ as defined by Eq.(2).

In order to work out quantitative criteria for the validity of Eq.(19) let us define an effective bandwidth $``$ of the fluctuating laser field $`S(\mathrm{\Omega })`$ by the relation

$$S(0)\pi _{\mathrm{}}^{\mathrm{}}𝑑\mathrm{\Omega }S(\mathrm{\Omega })=|\epsilon |^2.$$
(21)
The quantity $`1/`$ measures the correlation time of the fluctuating laser field. In the case of the PDM, for example, this effective bandwidth equals the width of the Lorentzian spectrum $`b`$. In all other cases which are described by Eq.(5) it characterizes the effective frequency width of the laser spectrum of Eq.(2). According to Eq.(19) the inverse rates $`1/_{ng}`$ define the characteristic time scale over which $`[\rho _{gg}(t)\rho _{nn}(t)]`$ varies significantly. Thus, the above decorrelation approximation applies to cases only for which $`_{ng}`$. Typically one finds $`_{n_{res}g}_{ng}`$ with $`ϵ_{n_{res}}=\overline{ϵ}`$ so that one of the validity conditions for the decorrelation condition becomes

$$2\pi |d_{n_{res}g}|^2S(0).$$
(22)
Eliminating $`S(0)`$ by Eq.(21) we finally arrive at the equivalent condition

$$^2\frac{1}{2}\mathrm{\Omega }_R^2$$
(23)
with the average Rabi-frequency $`\mathrm{\Omega }_R=2|d_{n_{res}g}|\sqrt{|\epsilon |^2}`$. What remains to be found is a validity condition for neglecting the coherences $`\rho _{nn^{}}(t)`$ with $`nn^{}`$ in Eq.(15). As long as

$$|ϵ_{n_{res}}ϵ_{n_{res}+1}|_{n_{res}g}=\frac{\mathrm{\Omega }_R^2}{2}$$
(24)
the coherences $`\rho _{nn^{}}(t)`$ with $`nn^{}`$ are rapidly oscillating functions in comparison with the slowly varying probabilities $`\rho _{nn}(t)`$ and $`\rho _{gg}(t)`$ entering Eq.(15) so that their influence averages to zero approximately. Therefore the inequalities (23) and (24) are the required conditions for the validity of the DCA. According to these conditions we may distinguish two limiting cases. In the limit of small laser bandwidths for which $`|ϵ_{n_{res}}ϵ_{n_{res}+1}|`$ they reduce to the requirement $`^2\mathrm{\Omega }_R^2/2`$. In two-level systems which are excited resonantly by a fluctuating laser field this is the well known limit of large laser bandwidths in which the dynamics are dominated by rate equations. In the opposite limit where the bandwidth is large enough to affect many excited Rydberg states, i.e. $`|ϵ_{n_{res}}ϵ_{n_{res}+1}|`$, the conditions for the applicability of the DCA reduce to the relation $`\gamma /\pi `$. Thereby we have introduced the laser-induced rate

$$\gamma =2\pi |d_{ϵg}|^2|\epsilon |^2\frac{\pi }{2}\mathrm{\Omega }_R^2|ϵ_{n_{res}}ϵ_{n_{res}+1}|^1$$
(25)
which characterizes ionization of the initial state $`|g`$ into continuum states close to the ionization threshold according to Fermi’s Golden rule.

Up to now, our arguments for the derivation of the rate equation (19) and of conditions (23) and (24) apply for discrete excited states only. However, our previous arguments can be generalized easily also to continuum states by viewing these continuum states as infinitesimally spaced discrete energy levels. According to quantum defect theory close to an ionization threshold the energy dependence of the discrete dipole matrix elements $`d_{ng}`$ is described by Eq.(12). Thus condition (24) reduces to

$$2|d_{ϵg}|^2|\epsilon |^2\frac{\gamma }{\pi }.$$
(26)
In the limit of an infinitesimally small level spacing between the excited states Eqs. (19) and (12) imply that the probability of finding the excited Rydberg system in a continuum state becomes vanishingly small. Thus we find

$$\dot{\rho }_{ϵϵ}(t)=_{ϵg}\rho _{gg}(t)$$
(27)
with $`_{ϵg}=2\pi |d_{ϵg}|^2S(\overline{ϵ}ϵ)`$. Integration over the whole electron continuum finally yields

$$\underset{ion}{\overset{\mathbf{}}{P}}(t)=\mathrm{\Gamma }\rho _{gg}(t)$$
(28)
with the mean ionization probability $`P_{ion}(t)=_0^{\mathrm{}}𝑑ϵ\rho _{ϵϵ}(t)`$ and with the effective ionization rate

$$\mathrm{\Gamma }=2\pi _0^{\mathrm{}}𝑑ϵ|d_{ϵg}|^2S(\overline{ϵ}ϵ).$$
(29)
If the mean excited energy $`\overline{ϵ}`$ is located well above threshold and if $`d_{ϵg}`$ is still energy independent over the energy region over which $`S(\overline{ϵ}ϵ)`$ is significant, this effective ionization $`\mathrm{\Gamma }`$ reduces to the previously introduced ionization rate $`\gamma `$ of Eq.(25).

The rate Eqs. (19) and (28) together with the conservation of probability, i.e.

$$\rho _{gg}(t)=1P_{ion}(t)\underset{n}{}\rho _{nn}(t)$$
(30)
and together with the initial condition $`\rho (t=0)=|gg|`$ determine the time evolution of a laser excited Rydberg electron within the framework of the DCA.

## IV Stochastic dynamics of Rydberg systems within the DCA

In this section the dynamics of a Rydberg system is investigated with the help of the DCA on the basis of Eqs.(19), (28) and (30).

Within the framework of a one-channel approximation the excited energies and the energy dependence of the relevant dipole matrix elements of a Rydberg system can be described by Eqs.(12). They are characterized by a quantum defect $`\alpha `$ and by an energy-normalized dipole matrix element $`d_{ϵg}`$ which are both approximately energy independent for $`|ϵ|1`$. The laser-induced coupling between the initial state $`|g`$ and the excited Rydberg- and continuum states is characterized by the ionization rate $`\gamma `$ of Eq.(25). Typically this description is adequate for Rydberg states of Alkali atoms.

The rate equations for the averaged density operator of the Rydberg system (compare with Eqs.(19),(28) and (30)) can be analyzed in a convenient way with the help of Laplace transformations. Defining the Laplace transformed density operator by

$$\stackrel{~}{\rho }(z)=_0^{\mathrm{}}𝑑te^{izt}\rho (t)$$
(31)
the associated inverse transformation is given by

$$\rho (t)\frac{1}{2\pi }_{\mathrm{}+i0}^{\mathrm{}+i0}𝑑ze^{izt}\stackrel{~}{\rho }(z).$$
(32)
Thus the Laplace transformed rate equations (19), (28) and (30) imply the relations

$`\stackrel{~}{\rho }_{gg}(z)`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Gamma }iz\sigma (z)}}`$ (33)
$`\text{and }\stackrel{~}{P}_{ion}(z)`$ $`=`$ $`{\displaystyle \frac{i\mathrm{\Gamma }}{z\left[\mathrm{\Gamma }iz\sigma (z)\right]}}`$ (34)
with

$$\sigma (z)=\underset{n}{}\frac{_{ng}}{_{ng}iz}.$$
(35)
The rates $`_{ng}`$ entering Eq.(35) characterize the transitions between states $`|g`$ and $`|n`$ within the DCA and are defined by Eq.(20). In the derivation of Eqs.(33) and (34) $`\rho _{gg}(t)`$ has been neglected in comparison with $`_n\rho _{nn}(t)`$ and $`P_{ion}(t)`$ in Eq.(30).

We may distinguish various dynamics regimes which are treated subsequently.

### A Asymptotic long time behaviour

The time evolution of the averaged density operator of the Rydberg system can be obtained from Eqs.(33) and (34) and from the inversion formula (32). In general the time evolution will exhibit both exponential decays originating from poles of the Laplace transforms (33) and (34) in the complex $`z`$-plane and power law decays which originate from cut contributions starting from the branch point of $`\sigma (z)`$ at $`z=0`$. As the asymptotic long time behaviour will be dominated by power law decays we have to investigate the structure of the characteristic kernel $`\sigma (z)`$ around the branch point $`z=0`$ in more detail. From Eqs.(35) and (20) it follows that for $`z0`$ its main contributions arise from the infinitely many Rydberg states very close to the ionization threshold. Hence in the long time limit we may approximate $`S(\overline{ϵ}ϵ_n)`$ by $`S(\overline{ϵ})`$ in expression (20). Furthermore we may replace the sum over all Rydberg states in Eq.(35) by an integration. So finally in the limit $`z0`$ we obtain the relation

$$\sigma (z)\frac{2\pi }{3\sqrt{3}}\left(\frac{i2\pi |d_{ϵg}|^2S(\overline{ϵ})}{z}\right)^{1/3}(z0).$$
(36)
Inserting Eq.(36) into Eqs.(32), (33) and (34) one obtains the asymptotic long time behaviour

$`\rho _{gg}(t)`$ $`=`$ $`\left({\displaystyle \frac{S(\overline{ϵ})}{|\epsilon |^2}}\right)^{1/3}{\displaystyle \frac{\mathrm{\Gamma }(\frac{5}{3})}{3(\mathrm{\Gamma }/\gamma )^2}}(t\gamma )^{5/3}`$ (37)
$`P_{ion}(t)`$ $`=`$ $`1\left({\displaystyle \frac{S(\overline{ϵ})}{|\epsilon |^2}}\right)^{1/3}{\displaystyle \frac{\mathrm{\Gamma }(\frac{2}{3})}{3(\mathrm{\Gamma }/\gamma )}}(t\gamma )^{2/3}`$ (38)
with $`\mathrm{\Gamma }(x)=_0^{\mathrm{}}𝑑ue^uu^{x1}`$ denoting the gamma function . Eqs.(37) and (38) describe the time evolution of the mean initial state probability and of the mean ionization probability for sufficiently long interaction times. They are generalizations of our previous results of Ref.  which were only valid for phase fluctuations of the PDM. Within the framework of the DCA these asymptotic laws are valid for arbitrary fluctuations of the laser field provided $`\gamma /\pi `$ (compare with Eq.(26)). Obviously this asymptotic long time behaviour is independent of the quantum defect which characterizes the influence of the ionic core of the Rydberg system. Furthermore, the characteristic exponents of the long time behaviour are a peculiar property of the Coulomb problem and do not depend on details of the laser spectrum. However, the time independent pre-factors of these power laws depend on $`S(\mathrm{\Omega })`$ and on the effective ionization rate $`\mathrm{\Gamma }`$ of Eqs.(2) and (29).

After which interaction time do we expect the asymptotic power laws of Eqs.(37) and (38) to become valid? As apparent from Eq.(38) this diffusive long time dynamics finally leads to complete ionization of the Rydberg system. Thus it is reasonable to characterize the onset of this asymptotic long time dynamics by a stochastic ionization time $`t_c`$ which is defined by the condition $`P_{ion}(t_c)=1/2`$ which yields

$$t_c=\frac{1}{\gamma }\left[\left\{\frac{\gamma \mathrm{\Gamma }(\frac{2}{3})}{\mathrm{\Gamma }}\right\}^3\frac{8}{27}\frac{S(\overline{ϵ})}{|\epsilon |^2}\right]^{1/2}.$$
(39)
### B Intermediate interaction times

In this section we deal with characteristic aspects of the dynamics described by Eqs.(19) and (28) in cases in which the interaction times are large enough so that the initial state $`|g`$ is depleted significantly but which are still much smaller than the stochastic ionization time $`t_c`$ of Eq.(39). For these interaction times we may distinguish two characteristic regimes depending on whether the excited states are close to the ionization threshold or above or whether they are located well below this threshold.

#### 1 Excitation at or above threshold

In this dynamical regime the significantly excited states are located at the ionization threshold or above. This implies that all rates $`_{ng}`$ which describe the coupling between states $`|g`$ and $`|n`$ are small in comparison with the total ionization rate $`\gamma `$. Thus considering interaction times $`t`$ which are not very much larger than $`1/\gamma `$ implies that we may take $`t_{ng}1`$ for all quantum numbers $`n`$. Hence, considering the formal solution of Eq.(19)

$$\rho _{nn}(t)=_{ng}_0^t𝑑se^{(ts)_{ng}}\rho _{gg}(s)$$
(40)
we may replace the exponential by unity. Performing the summation over all Rydberg states we find with the help of $`_n_{ng}\gamma \mathrm{\Gamma }`$ (the sum over all Rydberg states has been replaced by an integral)

$$\underset{n}{}\underset{nn}{\overset{\mathbf{}}{\rho }}(t)=(\gamma \mathrm{\Gamma })\rho _{gg}(t)\gamma ^2\mathrm{\Lambda }_{Sp}_0^t𝑑s\rho _{gg}(s)$$
(41)
where we used the quantity

$`\mathrm{\Lambda }_{Sp}`$ $`=`$ $`{\displaystyle \frac{1}{|\epsilon |^2^2}}{\displaystyle _{\mathrm{}}^0}𝑑ϵ_n({\displaystyle \frac{dϵ_n}{dn}})S^2(\overline{ϵ}ϵ_n)`$ (43)
$`{\displaystyle \frac{1}{|\epsilon |^2^2}}{\displaystyle _0^{\mathrm{}}}𝑑ϵS^2(\overline{ϵ}+ϵ)(2ϵ)^{3/2}`$
which characterizes the spectral properties of the fluctuating laser field.

We take the time-derivative of Eq.(30) and eliminate the ionization probability with Eq.(28). Inserting Eq.(41) in this equation we find an integro-differential equation for $`\rho _{gg}(t)`$, namely

$$\gamma \rho _{gg}(t)\gamma ^2\mathrm{\Lambda }_{Sp}_0^t𝑑s\rho _{gg}(s)+\underset{gg}{\overset{\mathbf{}}{\rho }}(t)=0.$$
(44)
As $`\mathrm{\Lambda }_{Sp}1`$ one finally obtains the relations

$$\rho _{gg}(t)=\frac{1}{1+\mathrm{\Lambda }_{Sp}}\left[\mathrm{exp}(\gamma t)+\mathrm{\Lambda }_{Sp}\mathrm{exp}(\gamma \mathrm{\Lambda }_{Sp}t)\right]$$
(45)
and

$$P_{ion}(t)=\frac{\mathrm{\Gamma }}{\gamma (1+\mathrm{\Lambda }_{Sp})}\left[\mathrm{exp}(\gamma \mathrm{\Lambda }_{Sp}t)\mathrm{exp}(\gamma t)\right].$$
(46)
These equations even apply to interaction times $`t<1/\gamma `$. Note that consistent with our approximations the interaction times always fulfill the inequality $`\gamma \mathrm{\Lambda }_{Sp}t1`$. In the special case of the PDM $`\mathrm{\Lambda }_{Sp}`$ reduces to

$$\mathrm{\Lambda }_{Sp}=\sqrt{\frac{b}{2\pi ^2}}\text{Re}\left\{\left(12i\frac{\overline{ϵ}}{b}\right)\sqrt{\frac{\overline{ϵ}}{b}i}\right\}.$$
(47)
Very close to threshold, i.e. for $`\overline{ϵ}0`$, one obtains $`\mathrm{\Lambda }_{Sp}=\sqrt{b}/(2\pi )`$ so that in this special case we obtain again our previous results of Ref. .

#### 2 Excitation well below threshold

If the fluctuating laser field excites Rydberg states well below the ionization threshold, i.e. $`|\overline{ϵ}|`$ and $`\overline{ϵ}<0`$, the considerations of Sec.IV B 1 have to be modified. Since the interaction time is assumed to be smaller than the stochastic ionization time $`t_c`$ we can approximate the characteristics of the dominantly excited Rydberg states by

$`ϵ_n`$ $``$ $`\overline{ϵ}+2\pi (nn_{res})/T_{\overline{ϵ}},`$ (48)
$`d_{ng}`$ $``$ $`d_{ϵg}\sqrt{2\pi /T_{\overline{ϵ}}}.`$ (49)
Thereby $`T_{\overline{ϵ}}2\pi (2\overline{ϵ})^{3/2}`$ denotes the classical Kepler period of the mean excited Rydberg state of energy $`\overline{ϵ}<0`$. Whereas in the previous subsection the stochastic influence could be characterized by a single average spectral property for arbitrary types of fluctuations, namely by $`\mathrm{\Lambda }_{Sp}`$, the excitation dynamics well below threshold turns out to be much more sensitive to the details of the laser spectrum. This is easily demonstrated by considering phase fluctuations which can be described by the laser spectrum of Eq.(9) as a particular example. This spectrum describes fluctuations of a single mode laser field well above the laser threshold in the limit $`\beta b`$. In addition, if only Rydberg states well below threshold are excited significantly we may neglect the effective ionization rate $`\mathrm{\Gamma }`$ in the denominator of Eqs.(33) and (34). In the limit $`\beta b`$ we thus arrive at the relations

$`\rho _{gg}(t)`$ $`=`$ $`{\displaystyle \frac{2}{T_{\overline{ϵ}}\beta }}f^{}\left({\displaystyle \frac{2\gamma b}{T_{\overline{ϵ}}\beta ^2}}t\right)`$ (50)
$`P_{ion}(t)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}\left[{\displaystyle \frac{\beta }{|\overline{ϵ}|}}\mathrm{arctan}\left({\displaystyle \frac{\beta }{|\overline{ϵ}|}}\right)\right]f\left({\displaystyle \frac{2\gamma b}{T_{\overline{ϵ}}\beta ^2}}t\right).`$ (51)
Thus within this limit for arbitrary values of $`\beta `$ and $`b`$ the influence of the phase fluctuations of the laser field is described by the single scaling function $`f(\tau )`$ which is defined by the equation

$$\frac{df(\tau )}{d\tau }=\text{Im}_0^{\mathrm{}}\frac{d\zeta }{\zeta }e^{i\zeta \tau }\left\{_{\mathrm{}}^{\mathrm{}}\frac{dx}{1i\zeta (x^2+x^4)}\right\}^1.$$
(52)
To end up with Eq.(52) we had to apply the further approximation $`\mathrm{\Omega }b`$ in the laser spectrum of Eq.(9). Physically speaking this approximation means that we consider cases in which the essential dynamics are dominated by energy states which are located in the wings of the laser spectrum.

In the limits $`\tau 1`$ and $`\tau 1`$ asymptotic expressions are easily obtained from Eq.(52). The limit of small values of $`\tau `$ is realized in the PDM where $`\beta \mathrm{}`$ and where the spectrum of Eq.(9) reduces to a Lorentzian form. In this case one obtains the expression

$$f(\tau )2\sqrt{\frac{\tau }{\pi }}(\tau 1).$$
(53)
Consequently Eqs.(50) and (51) yield

$`\rho _{gg}(t)`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{\pi bT_{\overline{ϵ}}\gamma t}}},`$ (54)
$`P_{ion}(t)`$ $`=`$ $`2\left[{\displaystyle \frac{1}{|\overline{ϵ}|}}{\displaystyle \frac{1}{\beta }}\mathrm{arctan}\left({\displaystyle \frac{\beta }{|\overline{ϵ}|}}\right)\right]\sqrt{{\displaystyle \frac{2t\gamma b}{\pi ^3T_{\overline{ϵ}}}}}.`$ (55)
From the numerical data shown in Fig.2 it is apparent that Eqs.(54) and (55) are good estimates for interaction times $`t<t_{\text{PDM}}=T_{\overline{ϵ}}\beta ^2/(200b\gamma )`$.

In the extreme opposite limit of large values of $`\tau `$ we obtain the relations

$$f(\tau )\frac{4}{3\pi }\mathrm{\Gamma }(\frac{1}{4})\tau ^{3/4}(\tau 1)$$
(56)
and

$`\rho _{gg}(t)={\displaystyle \frac{\mathrm{\Gamma }(\frac{1}{4})}{\pi }}\left({\displaystyle \frac{8}{T_{\overline{ϵ}}^3\beta ^2t\gamma b}}\right)^{1/4},`$ (57)
$`P_{ion}(t)=`$ (58)
$`{\displaystyle \frac{4}{3\pi ^2}}\mathrm{\Gamma }(1/4)\left[{\displaystyle \frac{\beta }{|\overline{ϵ}|}}\mathrm{arctan}\left({\displaystyle \frac{\beta }{|\overline{ϵ}|}}\right)\right]\left({\displaystyle \frac{2\gamma bt}{T_{\overline{ϵ}}\beta ^2}}\right)^{3/4}.`$ (59)
In this case the power law decays which characterize the diffusive dynamics of the excited Rydberg electron differ from the corresponding results of the PDM significantly. Even the characteristic exponents are changed. According to Fig.2 this dynamical regime is realized for interaction times $`t`$ which fulfill the relation $`t_{1/4}<tt_c`$ with $`t_{1/4}=T_{\overline{ϵ}}\beta ^2/(2b\gamma )`$.

## V Full master-equation

Due to their simplicity and their applicability to all types of laser spectra the DCA rate equations are ideal for understanding the dynamics of Rydberg systems in the case of large laser bandwidths (compare with Eqs.(23),(24) and (26)). However the DCA approximation is not capable of describing coherent aspects of the laser-induced excitation process. In order to investigate the limits of applicability of the DCA rate equations in this section a more general approach is developed which is also capable of describing all coherent aspects of the excitation process. For the sake of simplicity we shall restrict our subsequent discussion to the case of phase fluctuations of the exciting laser field which deviate only slightly from a Lorentzian spectrum and which can be modelled by Eqs.(4), (5) and (9). For these type of laser fluctuations we shall derive an approximate master equation involving the density matrix elements of the excited Rydberg electron which are averaged over the fluctuations of the laser field. This procedure is a generalization of previous approaches which so far have been applied to atomic few level system only . In the special case of the PDM these subsequently derived density matrix equations reduce to our previous results of Ref..

We start from the Schrödinger equation with Hamiltonian (10) with a fluctuating laser field as given by Eqs. (4) and (5). For the effective density operator

$$\rho (t)=\underset{k,jg,n}{}|kj|k|\psi (t)\psi (t)|je^{i(\mathrm{\Phi }(t)+\omega t)(\delta _{jg}\delta _{kg})}$$
(60)
we obtain the equation of motion

$$\stackrel{\mathbf{}}{\rho }(t)=i[H_{dr.},\rho (t)]i\varphi (t)[|gg|,\rho (t)].$$
(61)
The self-adjoint Hamiltonian

$$H_{dr.}=\underset{n,g}{}|nn|ϵ_n+|gg|\overline{ϵ}\epsilon _0\underset{n}{}\left(d_{ng}|ng|+\text{h.c.}\right)$$
(62)
describes the dynamics of the Rydberg system in the absence of phase fluctuations and the stochastic process $`\varphi (t)`$ is defined by Eq.(5). In order to average Eq.(61) over all possible realizations of the stochastic process $`\varphi (t)`$ it is convenient to introduce the averaged operators

$$\rho ^{(n)}(t)=\left(\frac{\beta }{b}\right)^{n/2}\frac{(i)^n}{\sqrt{n!}}_{\mathrm{}}^{\mathrm{}}𝑑\varphi Q_n(\varphi )\rho (t)p(\varphi ,t)$$
(63)
with $`n=0,1,2,3,\mathrm{}`$ . The (conditional) probability distribution $`p(\varphi ,t)`$ obeys the Fokker-Planck equation

$$[\frac{}{t}+]p(\varphi ,t)=0$$
(64)
with the Fokker-Planck operator

$$=\beta \frac{}{\varphi }\varphi +b\beta \frac{^2}{\varphi ^2}.$$
(65)
This equation has to be solved with the initial condition

$$p(\varphi ,0)=\frac{1}{\sqrt{2\beta b\pi }}\mathrm{exp}\left(\frac{\varphi ^2}{2\beta b}\right)$$
(66)
which represents the stationary solution of the Fokker-Planck equation. The quantities

$$Q_n(\varphi )=H_n\left(\varphi \left[2^{n+1}n!\beta b\right]^{1/2}\right)$$
(67)
with the Hermite polynomials $`H_n`$ are eigenfunctions of the adjoined Fokker-Planck operator $`^{}`$ with eigenvalues $`\mathrm{\Lambda }_n=n\beta `$ . Starting from Eq.(61) we obtain a set of coupled differential equations for the operators $`\rho ^{(n)}(t)`$, namely

$`\stackrel{\mathbf{}}{\rho ^{(n)}}(t)`$ $`=`$ $`i[H_{dr.},\rho ^{(n)}(t)]n\beta \rho ^{(n)}(t)`$ (68)
$`+`$ $`b(n+1)[|gg|,\rho ^{(n+1)}(t)]\beta [|gg|,\rho ^{(n1)}(t)].`$ (69)
These equations have to solved subject to the initial condition

$$\rho ^{(n)}(0)=\delta _{n0}\rho (0).$$
(70)
According to Eqs.(63) and (67) $`\rho ^{(0)}(t)`$ is the required density operator which is averaged over the phase fluctuations of the laser field.

We may derive an approximate equation of motion for the averaged density operator $`\rho ^{(0)}(t)`$. As a first step we Laplace transform Eqs.(68), i.e.

$`iz\stackrel{~}{\rho }^{(n)}(z)`$ $`=`$ $`\rho ^{(n)}(0)n\beta \stackrel{~}{\rho }^{(n)}(z)`$ (71)
$``$ $`i[H_{dr.},\stackrel{~}{\rho }^{(n)}(z)]+`$ (72)
$`+`$ $`[|gg|,b(n+1)\stackrel{~}{\rho }^{(n+1)}(z)\beta \stackrel{~}{\rho }^{(n1)}(z)].`$ (73)
In view of the large values of $`\beta `$ we are interested in we may neglect terms containing $`\epsilon _0`$ in comparison with terms containing $`\beta `$. Thus with the definition

$$\stackrel{~}{\alpha }_l^n(z)=\frac{\stackrel{~}{\rho }_{lg}^{n+1}(z)}{\stackrel{~}{\rho }_{lg}^n(z)}$$
(74)
we arrive at the recursion relations

$$n+1+\frac{i}{\beta }(ϵ_l\overline{ϵ}z)=\frac{b}{\beta }(n+2)\stackrel{~}{\alpha }_l^{n+1}(z)+\frac{1}{\stackrel{~}{\alpha }_l^n(z)}.$$
(75)
From these recursion relations we find

$$\stackrel{~}{\alpha }_k^0(z)=\frac{1}{1+i{\displaystyle \frac{ϵ_k\overline{ϵ}z}{\beta }}+{\displaystyle \frac{2b/\beta }{2+i{\displaystyle \frac{ϵ_k\overline{ϵ}z}{\beta }}+{\displaystyle \frac{3b/\beta }{3+\mathrm{}}}}}}.$$
(76)
Using Eqs.(74) and (76) we may now eliminate $`\stackrel{~}{\rho }^{(1)}(z)`$ in Eq.(72) for $`n=0`$. Performing the Laplace back transformation (32) and using the definition $`\rho (t)=\rho ^{(0)}(t)`$ we finally obtain the master equation

$`\stackrel{\mathbf{}}{\rho }(t)`$ $`=`$ $`i[H_{dr.},\rho (t)]`$ (77)
$``$ $`b{\displaystyle _0^t}𝑑\tau {\displaystyle \underset{kg}{}}\left\{\alpha _k^0(\tau )|kk|\rho (t\tau )|gg|+\text{h.c.}\right\}`$ (78)
with the memory function

$$\alpha _k^0(\tau )=\frac{1}{2\pi }_{\mathrm{}+i0}^{\mathrm{}+i0}𝑑ze^{iz\tau }\stackrel{~}{\alpha }_k^0(z).$$
(79)
In the limit of the PDM, i.e. for $`\beta \mathrm{}`$, this master equation reduces to the well known form

$`\stackrel{\mathbf{}}{\rho }(t)`$ $`=`$ $`i[H_{dr.},\rho (t)]`$ (80)
$`+`$ $`{\displaystyle \frac{1}{2}}\left\{[L,\rho (t)L^{}]+[L\rho (t),L^{}]\right\}`$ (81)
with the Lindblad operator

$$L=\sqrt{2b}|gg|.$$
(82)
## VI Numerical results

In this section numerical solutions of the master equation (77) are compared with the corresponding solutions of the DCA rate equations (19),(28) and (30). Details of the numerical technique for solving Eq.(77) are summarized in the appendix. On the basis of these comparisons the validity conditions for the applicability of the DCA and its accuracy can be tested. For this purpose we consider the laser excitation of a Rydberg system which can be described by quantum defect theory in a one-channel approximation (compare with Eqs.(12)). Typically this is a good approximation for Alkali atoms.

The time evolution of the mean initial state probability $`\rho _{gg}(t)`$ and of the mean ionization probability $`P_{ion}(t)`$ are depicted in Fig.3 for excitation at and well below the ionization threshold for different values of $`\beta `$. In both cases it is assumed that the exciting laser field has a well defined amplitude and a fluctuating phase.

Let us first turn to the case depicted in Fig.3a: The spectrum of the laser field is close to Lorentzian ($`\beta b`$), so that the asymptotic form Eq.(9) applies well. Thus the effective bandwidth $``$ as defined by Eq.(21) is approximately equal to the parameter $`b`$ which characterizes the spectrum of Eq.(9) and the parameter $`\beta `$ might be interpreted as an effective cut-off frequency of the laser spectrum. Rydberg states are excited by the fluctuating laser field well below the ionization threshold. The mean excited energy corresponds to a quantum number $`n_{res}=(2\overline{ϵ})^{1/2}=200`$. The laser bandwidth $`b`$ and the laser-induced rate $`\gamma `$ are so small that the excited Rydberg states are located well below threshold, i.e. $`\overline{ϵ}b,\gamma `$. However, the values of $`b`$ and $`\gamma `$ are large enough so that more than one Rydberg state around energy $`\overline{ϵ}`$ is affected significantly by the laser field, i.e. $`T_{\overline{ϵ}}\gamma ,T_{\overline{ϵ}}b>1`$. The three curves of Fig.3a (solid, dashed and long dashed) correspond to different values of the effective cut-off frequency $`\beta `$ of the laser spectrum. As many excited states are involved in the depletion of state $`|g`$ the initial stage of the time evolution is governed by an approximate exponential decay of state $`|g`$ with rate $`\gamma `$ . This initial stage of the time evolution is independent of the fluctuations of the laser field. At larger interaction times with $`t>1/\gamma `$ a coherent oscillation starts to appear in $`\rho _{gg}(t)`$ with the classical Kepler period $`T_{\overline{ϵ}}`$. This oscillation reflects the time evolution of the electronic Rydberg wave packet which has been prepared by the fast depletion of the initial state $`|g`$. With each return to the core region this Rydberg wave packet might undergo a transition to state $`|g`$ thus increasing $`\rho _{gg}(t)`$. These coherent oscillations cannot be described by the DCA rate equations. However, due to laser fluctuations after a few Kepler periods these coherent oscillations are damped out and merge into diffusive dynamics which is characterized by power law decay of the initial state $`|g`$. From this time on the dynamics of the Rydberg system under the influence of the fluctuating laser field is well described by the DCA rate equations. This is apparent from Fig.3a by comparing the numerical solutions of the master equation (solid, dashed and long dashed curves) with the asymptotic solutions of the DCA rate equations (circles and thin dashed curves). According to the discussion presented in Sec. IV B 2 this diffusive dynamical behaviour appears when all coherent effects are damped out and disappears again at interaction times $`t>t_c`$ at which stochastic ionization starts to dominate. Physically speaking for these intermediate interaction times the excited electronic Rydberg wave packet starts to diffuse in energy space towards the ionization threshold. It reaches the ionization threshold roughly at time $`t_c`$ at which the ionization probability rises significantly from vanishingly small values to values close to unity. The early stages of this diffusion towards the ionization threshold are governed by the power law decay of Eq.(54) for $`\rho _{gg}(t)`$ which is characterized by the exponent $`(1/2)`$. In the case of laser fluctuations which can be described by the PDM to a good degree of approximation this power law decay governs the time evolution up to the stochastic ionization time $`t_c`$. However, for non-Lorentzian spectra this is no longer the case. In cases in which the non-Lorentzian effects can be described by the spectrum of Eq.(9) the characteristic exponent of this power law decay is changed to a value of $`(1/4)`$ as soon as the interaction times become larger than the characteristic time $`t_{1/4}=T_{\overline{ϵ}}\beta ^2/(2b\gamma )`$ provided $`t_{1/4}<t_c`$ (compare with Eq.(57)). This non-Lorentzian effect is clearly apparent in Fig.3a where the characteristic times $`t_{1/4}`$ are indicated for $`b/\beta =0.03`$ and $`b/\beta =0.2`$. With increasing values of $`\beta `$ this characteristic time increases and this cross over phenomenon disappears for sufficiently large values of $`\beta `$ as soon as $`t_{1/4}>t_c`$. At interaction times exceeding the stochastic ionization time $`t_c`$ the excited Rydberg wave packet has already reached the ionization threshold and the mean ionization probability rises to a value of unity. This asymptotic long time behaviour of the excitation dynamics is described to a good degree of approximation by Eqs.(37) and (38). That is apparent from a comparison of the thin dashed lines of Fig.3a with the corresponding numerical solutions of the master equation (77).

In Fig.3b the laser bandwidth is so large that the significantly excited energy region $`[\overline{ϵ},\overline{ϵ}+]`$ (before the onset of the electronic diffusion process) contains already the ionization threshold. Thus the excited Rydberg system is ionized significantly already in the early stages of the time evolution. As $`\gamma `$ this early stage of the ionization process is well described by the DCA rate equations which yield (compare to Eqs.(45) and (46))

$`\rho _{gg}(t)`$ $`=`$ $`e^{\gamma t},`$ (83)
$`P_{ion}(t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }}{\gamma }}(1e^{\gamma t}).`$ (84)
These approximate solutions are obtained from Eqs.(19) and (28) by neglecting $`\rho _{nn}(t)`$ in comparison with $`\rho _{gg}(t)`$. According to Eqs.(45) and (46) this initial ionization process saturates as soon as the mean initial state probability and the mean ionization probability have reached the values $`\mathrm{\Lambda }_{Sp}`$ and $`\mathrm{\Gamma }/\gamma `$. At these interaction times we still have $`\gamma \mathrm{\Lambda }_{Sp}t1`$. These characteristic aspects of the laser-induced excitation process are clearly apparent in Fig.3b. Physically speaking in this initial stage of the excitation process the Rydberg electron is ionized with probability $`\mathrm{\Gamma }/\gamma `$. With a probability of $`(1\mathrm{\Gamma }/\gamma )`$ an excited electronic Rydberg wave packet is prepared after a time of the order of $`1/\gamma `$. This wave packet is formed by a coherent superposition of all Rydberg states within the dominantly excited energy interval $`[\overline{ϵ},0]`$. Depending on the actual value of the laser bandwidth the coherent dynamics of this electronic wave packet is damped sooner or later. After the destruction of all coherences to a good degree of approximation the subsequent dynamics is governed by the DCA rate equations. In Fig. 3b small coherence oscillations are visible in the time evolution of $`\rho _{gg}(t)`$. As soon as the interaction time exceeds the stochastic ionization time the excited Rydberg electron starts to ionize significantly. The time evolution of this stochastic ionization process is well described by Eqs.(37) and (38) within the framework of the DCA rate equations.

## VII Summary and conclusion

The dynamics of an electronic Rydberg wave packet under the influence of a fluctuating cw-laser field has been discussed. It has been shown that for large laser bandwidths its dynamics can be described by Pauli-type rate equations for the relevant density matrix elements of the excited Rydberg electron averaged over the laser fluctuations. These rate equations are valid for arbitrary types of laser fluctuations and their dynamics is determined by the spectrum of the laser field only and not by any of the higher order correlation functions. The validity of these rate equations has been investigated in detail for a special class of phase fluctuations of the laser field.

With the help of these rate equations we have investigated the dynamics of a laser excited Rydberg electron for long interaction times. At these interaction times the dynamics of the Rydberg electron are dominated by stochastic diffusion in energy space towards the ionization threshold which leads finally to stochastic ionization. This diffusion process is accompanied by a characteristic scenario of power law decays. Analytical expressions have been derived for these power laws and their associated characteristic exponents. These analytical expressions exhibit in a clear way that the asymptotic power laws are independent of the quantum defect of the excited Rydberg states and to which extent they depend on details of the laser spectrum. In particular, it has been demonstrated that the characteristic exponents which describe the process of stochastic ionization are completely independent of the laser spectrum. However, the initial stages of the diffusion of the excited Rydberg electron depend on details of the laser spectrum.

Support by the Deutsche Forschungsgemeinschaft within the SPP ‘Zeitabhängige Phänomene und Methoden’ is acknowledged.

## A Numerical solution technique for the Master equation

In this appendix an efficient numerical method for solving the master equation (77) is outlined.

First, we make a rearrangement of Eq.(77) splitting it into a PDM-part and an additional part that disappears in the limit $`\beta \mathrm{}`$, namely

$$\stackrel{\mathbf{}}{\rho }(t)=i\left[H_{eff}\rho (t)\rho (t)H_{eff}^{}\right]+𝐋\rho (t).$$
(A1)
Thereby we have introduced the effective non-Hermitian Hamiltonian

$$H_{eff}=H_{dr.}ib|gg|,$$
(A2)
the damping operator

$`𝐋\rho (t)`$ $`=`$ $`2b|gg|\rho (t)|gg|`$ (A3)
$`+`$ $`b{\displaystyle _0^t}𝑑\tau {\displaystyle \underset{ng}{}}\left\{|nn|\rho (t)|gg|w_n(\tau )+\text{h.c.}\right\}`$ (A4)
and a memory function

$$w_n=\underset{\mathrm{\Omega }\mathrm{}}{lim}\mathrm{\Omega }e^{\mathrm{\Omega }\tau }\alpha _n^0(\tau ).$$
(A5)
In the PDM-limit, i.e. $`\beta \mathrm{}`$, the second term on the right side of Eq.(A3) disappears and the master equation reduces to Eq.(80). Integration of Eq.(A1) yields

$`\rho (t)=U(t)\rho (0)U^{}(t)+`$ (A6)
$`{\displaystyle _0^{\mathrm{}}}𝑑t^{}\mathrm{\Theta }(tt^{})U(tt^{})𝐋\rho (t^{})U^{}(tt^{})`$ (A7)
where $`U(t)=\mathrm{exp}[iH_{eff}t]`$ is a non unitary time evolution operator. If the initial condition is taken to be $`\rho (0)=|gg|`$, the Laplace transformed matrix elements of the density operator $`\stackrel{~}{\rho }_{ij}(z)_z\rho _{ij}(t)_0^{\mathrm{}}𝑑te^{izt}\rho _{ij}(t)`$ become

$`\stackrel{~}{\rho }_{lk}(z)=H_{lggk}(z)\left[1+2b\stackrel{~}{\rho }_{gg}(z)\right]+`$ (A8)
$`b{\displaystyle \underset{ng}{}}\left\{H_{lngk}(z)\stackrel{~}{w}_n(z)\stackrel{~}{\rho }_{ng}(z)+H_{lgnk}(z)\stackrel{~}{w}_n^{}(z)\stackrel{~}{\rho }_{gn}(z)\right\},`$ (A9)
$`H_{abcd}(z)=`$ (A10)
$`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑z_1_{z_1+i0}a\left|U(t)\right|b\left[_{z_1z+i0}d\left|U(t)\right|c\right]^{}`$ (A11)
with $`\stackrel{~}{w}_n(z)=1\stackrel{~}{\alpha }_n^0(z)`$. The Laplace transformed transition amplitudes $`_zi\left|U(t)\right|j`$ appearing in Eq.(A10) are easily calculated .

$`_zg\left|U(t)\right|g={\displaystyle \frac{i}{z+ib\overline{ϵ}\mathrm{\Sigma }(z)}},`$ (A12)
$`_zg\left|U(t)\right|n=_zn\left|U(t)\right|g`$ (A13)
$`={\displaystyle \frac{i\epsilon _0d_{ng}}{(zϵ_n)[z+ib\overline{ϵ}\mathrm{\Sigma }(z)]}},`$ (A14)
$`_zn\left|U(t)\right|l=`$ (A15)
$`{\displaystyle \frac{1}{zϵ_n}}\left\{i\delta _{nl}+{\displaystyle \frac{i\epsilon _0^2d_{ng}d_{lg}}{(zϵ_l)[z+ib\overline{ϵ}\mathrm{\Sigma }(z)]}}\right\}`$ (A16)
where $`\mathrm{\Sigma }(z)=_{ng}\frac{|\epsilon _0d_{ng}|^2}{zϵ_n}`$ is the self energy of state $`|g`$ and the Kronecker-symbol $`\delta _{nl}`$ turns into a Dirac-delta function $`\delta (ϵ_nϵ_l)`$ for energy normalized continuum states $`|ϵ_n`$ and $`|ϵ_l`$. For the subsequent treatment it is convenient to introduce the expressions

$`\alpha _n(z)`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{\rho }_{ng}(z)}{[1+2b\stackrel{~}{\rho }_{gg}(z)]\epsilon _0d_{ng}}},`$ (A17)
$`\beta _n(z)`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{\rho }_{gn}(z)}{[1+2b\stackrel{~}{\rho }_{gg}(z)]\epsilon _0d_{ng}}},`$ (A18)
$`E_n(z)`$ $`=`$ $`{\displaystyle \frac{H_{gngg}(z)}{\epsilon _0d_{ng}}},F_n(z)={\displaystyle \frac{H_{ggng}(z)}{\epsilon _0d_{ng}}}.`$ (A19)
Thus Eq.(A8) and Laplace transformation of Eq.(A1) yield

$`\stackrel{~}{\rho }_{gg}(z)={\displaystyle \frac{𝒦(z)}{12b𝒦(z)}},`$ (A20)
$`\stackrel{~}{P}_{ion}(z)=`$ (A21)
$`{\displaystyle \frac{1}{2\pi (z+i0)(12b𝒦(z))}}{\displaystyle _0^{\mathrm{}}}𝑑ϵ\left[\alpha _ϵ(z)\beta _ϵ(z)\right]`$ (A22)
with the characteristic kernel

$$𝒦(z)=H_{gggg}(z)+b\left(\underset{ng}{}+_0^{\mathrm{}}𝑑n(ϵ)\right)|\epsilon _0d_{ng}|^2\left\{E_n(z)\stackrel{~}{w}_n(z)\alpha _n(z)+F_n(z)\stackrel{~}{w}_n^{}(z)\beta _n(z)\right\}$$
(A23)
and with

$`\alpha _l(z)`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{w}_l^{}(z)C_l(z)|\epsilon _0d_{lg}|^2\left(F_l(z)+S_l(z)\right)+\left(1J_l(z)\right)\left(E_l(z)+T_l(z)\right)}{\left(1G_l(z)\right)\left(1J_l(z)\right)(\stackrel{~}{w}_l(z)\stackrel{~}{w}_l^{}(z))^2C_l(z)^2|d_{lg}|^4}},`$ (A24)
$`\beta _l(z)`$ $`=`$ $`{\displaystyle \frac{w_l(z)C_l(z)|\epsilon _0d_{lg}|^2\left(E_l(z)+T_l(z)\right)+\left(1G_l(z)\right)\left(F_l(z)+S_l(z)\right)}{\left(1G_l(z)\right)\left(1J_l(z)\right)(\stackrel{~}{w}_l(z)\stackrel{~}{w}_l^{}(z))^2C_l(z)^2|d_{lg}|^4}}.`$ (A25)
The non-diagonal couplings $`(\alpha _l,\beta _l\alpha _n,\beta _n,nl)`$ due to the sum in Eq.(A8) give rise to the expressions $`S_l(z)`$ and $`T_l(z)`$ appearing in Eqs.(A24) and (A25), namely

$`S_l(z)`$ $`=`$ $`b\left\{{\displaystyle \underset{n\{g,l\}}{}}+{\displaystyle _0^{\mathrm{}}}𝑑n(ϵ)\right\}|\epsilon _0d_{ng}|^2\left[\beta _n(z){\displaystyle \frac{F_n(z)F_l(z)}{ϵ_lϵ_n}}\alpha _n(z){\displaystyle \frac{F_l(z)E_n(z)}{ϵ_nϵ_lz2i0}}\right],`$ (A26)
$`T_l(z)`$ $`=`$ $`b\left\{{\displaystyle \underset{n\{g,l\}}{}}+{\displaystyle _0^{\mathrm{}}}𝑑n(ϵ)\right\}|\epsilon _0d_{ng}|^2\left[\alpha _n(z){\displaystyle \frac{E_n(z)E_l(z)}{ϵ_lϵ_n}}\beta _n(z){\displaystyle \frac{E_l(z)F_n(z)}{ϵ_nϵ_l+z+2i0}}\right].`$ (A27)
The diagonal couplings of the $`\alpha _n`$ and $`\beta _n`$ yield $`J_l(z),G_l(z)`$ and $`C_l(z)`$, i.e.

$`C_l(z)`$ $`=`$ $`b\mathrm{\Theta }(ϵ_l){\displaystyle \frac{E_l(z)F_l(z)}{z+2i0}},`$ (A28)
$`G_l(z)`$ $`=`$ $`{\displaystyle \frac{b\stackrel{~}{w}_l(z)}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dz_1\left\{1+\mathrm{\Theta }(ϵ_l)|\epsilon _0d_{lg}|^2\left[z_1+ib\overline{ϵ}\mathrm{\Sigma }(z_1+i0)\right]^1(z_1ϵ_l+i0)^1\right\}}{\left[z_1zib\overline{ϵ}\mathrm{\Sigma }(z_1zi0)\right](z_1ϵ_l+i0)}},`$ (A29)
$`J_l(z)`$ $`=`$ $`{\displaystyle \frac{b\stackrel{~}{w}_l^{}(z)}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dz_1\left\{1+\mathrm{\Theta }(ϵ_l)|\epsilon _0d_{lg}|^2\left[z_1zib\overline{ϵ}\mathrm{\Sigma }(z_1zi0)\right]^1(z_1zϵ_li0)^1\right\}}{\left[z_1+ib\overline{ϵ}\mathrm{\Sigma }(z_1+i0)\right](z_1zϵ_li0)}}.`$ (A30)
All information about the quantum system is contained in the self energy $`\mathrm{\Sigma }(z)`$. In order to solve Eqs.(A24-A27) for a given value of $`z`$, the coefficients $`H_{gggg},E_l,F_l,J_l`$ and $`G_l`$ have to be calculated. Using the energy levels and dipole matrix elements of Eqs.(12) the self energy becomes

$`\mathrm{\Sigma }(z)`$ $`=`$ $`\delta \omega i{\displaystyle \frac{\gamma }{2}}+i\gamma {\displaystyle \frac{1}{1\mathrm{exp}(2\pi i\nu (z))]}}`$ (A31)
$`\text{with}\nu (z)`$ $`=`$ $`(2z)^{1/2}+\alpha `$ (A32)
and with the (non-resonant) quadratic Stark-shift contribution $`\delta \omega `$. In we calculated the quantity $`H_{gggg}`$ ($`f(z)`$ in that work) by contour integration. In an analogous way the quantities $`E_l,F_l,J_l`$ and $`G_l`$ can be calculated but for sake of brevity we do not give them here explicitly. Starting with $`S_n=T_n=0`$, Eqs.(A20-A27) are solved by iteration. Actually we found the non-diagonal coupling terms $`S_l,T_l`$ to be very small in comparison with the diagonal couplings so that this iteration converges very rapidly.


format: text

text:
# Instabilities of dispersion-managed solitons in the normal dispersion regime

## 1 Introduction

New ways in optimization of existing telecommunication systems based on dispersion management technology attracted recently wide research interest from soliton-based groups (see reviews ). The main idea was to combine a high local group-velocity dispersion with a low path-average dispersion. The former feature results in the reduction of the four-wave mixing while the latter one reduces the Gordon–Haus timing jitter effects. When the path-average dispersion is small and normal, i.e. the defocussing segment in the fiber is dominant over the focussing one, a new phenomenon of branching of soliton solutions was discovered . The soliton propagation in this regime is not supported by an uniform-dispersion optical fiber and seems to be one of the remarkable achievement of the dispersion management with sufficiently high local dispersion.

The stability of branching soliton solutions in the normal dispersion regime was a subject of intense studies which lead to contradictory conclusions. Grigoryan and Menyuk announced the linear and nonlinear stability of both the branches , while Berntson et al. conjectured instability of one of the branches .

In this paper, we intend to shed light on the complicated issue of existence and stability of soliton signals in the normal regime of the dispersion map. We find, in the small-amplitude approximation, that there exist two branches of soliton solutions for different levels of energy and different pulse durations at a fixed propagation constant. The short pulses with larger energy are proved to be linearly unstable, while the other (long) pulses with smaller energy are neutrally stable. We show that the transition from large-energy unstable solitons to the stable soliton signals occurs via long-term transient oscillations. The two branches of soliton solutions correspond to a single (small-energy) branch B in Fig. 1 of . Depending on a normalization condition (see Section 4), this branch may be either stable or unstable.

Our strategy to develop the small-amplitude approximation is based on the combination of two analytical approaches: the Gaussian variational approximation and the integral evolution model.

The Gaussian variational approximation, being inaccurate in details, is still useful for a quick and rough analysis (see  for review and references). Also it was shown that the method can be extended to a rigorous Gauss-Hermite expansion of the basic model . We improve the previous results summarized in  by deriving a new dynamical system from the variational equations of a Gaussian pulse. The system clearly displays the linear and nonlinear instability of the short Gaussian pulse with larger energy.

More rigorous analysis of the problem is based on the integral evolution model obtained by Gabitov and Turitsyn  and by Ablowitz and Biondini . Although this model is more complicated from computational point of view (see recent papers ), we managed to study numerically the construction of the linear spectrum of dispersion-managed solitons. Our results confirm the instability and transition scenarios predicted within the variational model. We also deduce from this model that the soliton signals in the normal regime of the dispersion map are in resonance with the wave continuum of linear excitations of the map. The resonance implies isually the generation of wave packets from stable pulses oscillating in time. The latter effects are beyond the accuracy of the analytical model and are left for further studies.

## 2 Gaussian approximation: New dynamical model

We study the NLS model in the dimensionless form ,

$$iu_z+\frac{1}{2}D(z)u_{tt}+ϵ\left(\frac{1}{2}D_0u_{tt}+|u|^2u\right)=0,$$
(1)
where $`u(z,t)`$ is the envelope of an optical pulse in the retarded reference frame of the fiber. The small parameter $`ϵ`$ measures the length of the dispersion’s map and the inverse variance of the local dispersion. After normalization, $`D_0`$ and $`D(z)`$ are assumed to be of order of $`\mathrm{O}(1)`$, and

$$D=_0^1D(z)𝑑z=0,D(z+1)=D(z).$$
(2)
Further physical motivations for derivation and verification of the model (1) can be found in . Soliton-like optical pulses are solutions of the model in the form,

$$u(z,t)=\psi (z,t)e^{iϵ\mu z},$$
(3)
where $`\mu `$ is the propagation constant and $`\psi (z,t)`$ is a soliton profile satisfying the boundary conditions,

$$\psi (z+1,t)=\psi (z,t)$$
(4)
and

$$\underset{t\mathrm{}}{lim}\psi (z,t)=0.$$
(5)
One of the conventional approximation for soliton solutions of NLS-type equations is based on averaging the Gaussian anzatz in the Lagrangian density and further varying the Lagrangian density with respect to parameters of the Gaussian pulse (see  for review). The Gaussian approximation is the first term of the Gauss-Hermite expansions when solving the NLS equation (1) in the limit $`ϵ0`$ . We apply the Gaussian anzatz in the form,

$$u(z,t)=c(z)\mathrm{exp}\left(\frac{(\alpha (z)2i\beta (z))}{\alpha (z)^2+4\beta (z)^2}t^2+i\varphi (z)\right).$$
(6)
Here the four parameters of the Gaussian pulse are: $`c(z)`$ \- the amplitude, $`\varphi (z)`$ \- the gauge parameter, $`\alpha (z)`$ \- the pulse duration, and $`\beta (z)`$ \- the chirp. It was found that the four equations for variations of the Lagrangian density can be decoupled into a system for $`\alpha (z)`$ and $`\beta (z)`$ of the form,

$`{\displaystyle \frac{d\alpha }{dz}}`$ $`=`$ $`{\displaystyle \frac{4ϵE\alpha ^{5/2}\beta }{(\alpha ^2+4\beta ^2)^{3/2}}},`$ (7)
$`{\displaystyle \frac{d\beta }{dz}}`$ $`=`$ $`D(z)+ϵ\left(D_0{\displaystyle \frac{E\alpha ^{3/2}(\alpha ^24\beta ^2)}{2(\alpha ^2+4\beta ^2)^{3/2}}}\right).`$ (8)
The phase factor $`\varphi (z)`$ is expressed in terms of $`\alpha (z)`$ and $`\beta (z)`$,

$$\frac{d}{dz}\left(\varphi +\frac{1}{2}\mathrm{arctan}\frac{2\beta }{\alpha }\right)=\frac{ϵE\alpha ^{1/2}(3\alpha ^2+20\beta ^2)}{4(\alpha ^2+4\beta ^2)^{3/2}},$$
(9)
while the amplitude $`c(z)`$ is given in terms of the input energy constant $`E`$ as

$$E=\frac{\sqrt{\alpha ^2+4\beta ^2}}{\sqrt{2\alpha }}c^2=\frac{1}{\sqrt{\pi }}_0^1𝑑z_{\mathrm{}}^{\mathrm{}}𝑑t|u|^2(z,t)>0.$$
(10)
The stationary pulse (3) - (5) corresponds to the periodic solutions of the system (7) and (8) in the form,

$$\alpha (z+1)=\alpha (z),\beta (z+1)=\beta (z),\varphi (z+1)=\varphi (z)+ϵ\mu .$$
(11)
For simplicity, we study the periodic solutions in the limit $`ϵ0`$ by using a two-step dispersion map with zero average,

$$D(z)=\{\begin{array}{c}D_1,\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}<z<L\\ D_2,L<z<1\end{array},$$
(12)
where

$$m=D_1L=D_2(1L)>0.$$

The asymptotic solution in the limit $`ϵ0`$ can be sought in the regular form,

$$\alpha (z)=\alpha _s+\mathrm{O}(ϵ),\beta (z)=_0^zD(z^{})𝑑z^{}+\beta _s+\mathrm{O}(ϵ),$$

where $`\alpha _s`$, $`\beta _s`$ are constant. The periodic solutions appear when $`\beta _s=m/2`$ and $`\alpha _s`$ is a root of the equation,

$$D_0=E\alpha _s^{3/2}\left[\frac{1}{(m^2+\alpha _s^2)^{1/2}}\frac{1}{2m}\mathrm{log}\left(\frac{m+(m^2+\alpha _s^2)^{1/2}}{\alpha _s}\right)\right].$$
(13)
In addition, the propagation constant $`\mu `$ can be obtained as a function of $`E`$ and $`\alpha _s`$ according to the equation,

$$\mu =\frac{1}{4}E\alpha _s^{1/2}\left[\frac{2}{(m^2+\alpha _s^2)^{1/2}}+\frac{5}{m}\mathrm{log}\left(\frac{m+(m^2+\alpha _s^2)^{1/2}}{\alpha _s}\right)\right].$$
(14)
These equations have been already derived in the literature, see  for (13) and  for (14). However, the relations (13) and (14) were viewed typically under the normalization condition,

$$m=1,E=\frac{1}{\sqrt{2S}},\alpha _s=\frac{1}{S},$$
(15)
where $`S`$ is called the map strength. In this normalization, the expression (13) gives a small-amplitude limit of the results of , i.e. the slope $`E/D_0`$ is a function of $`S`$. The existence of solitons was identified both for $`D_0>0`$ (when $`S<S_{thr}`$) and for $`D_0<0`$ (when $`S<S_{thr}`$), where $`S_{thr}3.32`$.

In this paper, we develop a different frame to view the soliton solutions (13) and (14). Guided by the stability analysis of solitons in generalized NLS equations , we fix the parameters of the model ($`D_0,m`$) and construct periodic solutions as a one-parameter family in terms of the propagation constant $`\mu `$. As a result, the parameters $`\alpha `$ and $`E`$ can be found from (13)-(14) as $`\alpha _s=\alpha _s(\mu )`$ and $`E=E(\mu )`$. These functions are shown in Fig. 1(a,b) for $`D_0=0.02`$ and $`m=2`$ and in Fig. 2(a,b) for $`D_0=0.02`$ and $`m=2`$. Obviously, the branching occurs at $`D_0<0`$ (i.e. at the normal regime of the dispersion map), when the dispersion map is defocussing on average. The two solutions coexist for $`\mu >\mu _{thr}(D_0,m)`$ and $`E>E_{thr}(D_0,m)`$. Both the branches I and II correspond to a single branch B in the small-amplitude approximation under the normalization condition (15) .

In order to describe non-stationary dynamics of the Gaussian pulse near the periodic solutions, we derive a dynamical model from Eqs. (7) and (8) by setting,

$$\alpha (z)=\alpha _0(\zeta )+ϵ\alpha _1(z,\zeta )+\mathrm{O}(ϵ^2)$$

and

$$\beta (z)=_0^zD(z^{})𝑑z^{}+\beta _0(\zeta )+ϵ\beta _1(z,\zeta )+\mathrm{O}(ϵ^2).$$

Here $`\zeta =ϵz`$ is the distance to measure the evolution of a Gaussian pulse over many map’s periods. The coupled system (7) and (8) can be averaged over the map’s period subject to the periodic conditions: $`\alpha _1(z+1,\zeta )=\alpha _1(z,\zeta )`$ and $`\beta _1(z+1,\zeta )=\beta _1(z,\zeta )`$. Then, the non-stationary system reduces to the dynamical model for $`\alpha _0(\zeta )`$ and $`\beta (\zeta )`$,

$`{\displaystyle \frac{d\alpha _0}{d\zeta }}=F_\alpha (\alpha _0,\beta _0)`$ $``$ $`{\displaystyle \frac{E\alpha _0^{5/2}}{m}}\left[{\displaystyle \frac{1}{(\alpha _0^2+4\beta _0^2)^{1/2}}}{\displaystyle \frac{1}{(\alpha _0^2+4(\beta _0+m)^2)^{1/2}}}\right],`$ (16)
$`{\displaystyle \frac{d\beta _0}{d\zeta }}=F_\beta (\alpha _0,\beta _0)`$ $``$ $`D_0{\displaystyle \frac{E\alpha _0^{3/2}}{4m}}[{\displaystyle \frac{4(m+\beta _0)}{(\alpha _0^2+4(m+\beta _0)^2)^{1/2}}}{\displaystyle \frac{4\beta _0}{(\alpha _0^2+4\beta _0^2)^{1/2}}}`$ (17)
$`+`$ $`\mathrm{log}\left({\displaystyle \frac{2\beta _0+(\alpha _0+4\beta _0^2)^{1/2}}{2(m+\beta _0)+(\alpha _0^2+4(m+\beta _0)^2)^{1/2}}}\right)].`$
This system has of course the same stationary solutions $`\alpha _0=\alpha _s`$ and $`\beta _0=\beta _s=m/2`$ as those given in (13). The stationary solutions appear as equilibrium states in the dynamical system, whose stability can be found by linearizing,

$$\alpha _0(\zeta )=\alpha _s+\mathrm{\Delta }\alpha e^{i\lambda \zeta },$$

$$\beta _0(\zeta )=\beta _s+\mathrm{\Delta }\beta e^{i\lambda \zeta },$$

where the eigenvalue $`\lambda `$ is

$$\lambda ^2(\mu )=\frac{F_\alpha }{\beta }(\alpha _s)\frac{F_\beta }{\alpha }(\alpha _s)=\frac{2E\alpha _s^{3/2}}{(m^2+\alpha _s^2)^{3/2}}\left[3D_0+\frac{E\alpha _s^{3/2}(m^2\alpha _s^2)}{(m^2+\alpha _s^2)^{3/2}}\right].$$
(18)
We plot $`\lambda ^2(\mu )`$ in Fig.3 to confirm that $`\lambda ^2>0`$ for branch I of the periodic solutions and $`\lambda ^2<0`$ for branch II (cf. Fig. 2). Thus, the linear analysis predicts the instability of the short Gaussian pulses with larger energy at a fixed propagation constant $`\mu `$ (branch II). In the limit $`\mu \mu _{thr}(D_0,m)`$, the instability disappears, i.e.

$$\underset{\mu \mu _{thr}}{lim}\lambda ^2(\mu )=0.$$

This property follows from Eq. (17) in the limit $`\mu \mu _{thr}(D_0,m)`$, when

$$\frac{F_\beta }{E}(\alpha _s,\beta _s)\frac{dE}{d\mu }+\frac{F_\beta }{\alpha }(\alpha _s,\beta _s)\frac{d\alpha }{d\mu }=0.$$
(19)
Connecting Eqs. (18) and (19), we find the following asymptotic approximation,

$$\lambda ^2\left(\frac{F_\beta /EF_\alpha /\beta }{d\alpha /d\mu }\right)\frac{dE}{d\mu }.$$
(20)
Taking into account that $`dE/d\mu `$, $`F_\beta /E`$, and $`F_\alpha /\beta `$ are all positive for $`D_0<0`$, and $`\alpha (\mu \mu _{thr})^{1/2}`$ (see Fig. 2(a,b)), the asymptotic approximation (20) produces the result, $`\lambda (\mu \mu _{thr})^{1/4}`$.

The nonlinear dynamics of the system (16) and (17) is shown in Fig. 4 for $`D_0=0.02`$ and $`m=2`$. At a fixed value of the energy $`E`$, there are two stationary Gaussian pulses of different durations: a short pulse is a saddle point, while a long one is a center. Inside the separatrix loop, there are oscillations of the pulse trapped by the center point. Outside the separatrix, the Gaussian pulse transfers to chirped linear waves.

We notice that the transition scenario resembles the nonlinear dynamics of unstable solitons in generalized NLS equations . The only difference is that the unstable branch in generalized NLS equations is located for those values of soliton propagation constant $`\mu `$, where $`dE/d\mu <0`$. Although this conventional stability criterion failed for the dispersion-managed solitons (see Fig. 2(a)), the instability development shows up to be alike (cf. Fig. 4 here and Fig. 2(b) in ).

## 3 Integral evolution model: Numerical analysis

The Gaussian approximation of the optical pulse in the NLS model (1) can be improved by summating all higher-order Gauss-Hermite solutions of the linear equation, $`iu_z+0.5D(z)u_{tt}=0`$ as shown in . However, this analysis results in a complicated infinite-dimensional system of algebraic equations for parameters and coefficients of the Gauss-Hermite expansion. Instead, we adopt a direct asymptotic method to deduce an integral evolution model valid in the limit $`ϵ0`$. This method is based on Fourier expansion of solutions of the linear equation above  as well as on the asymptotic expansion,

$$u(z,t)=u_0(z,t)+ϵu_1(z,t)+\mathrm{O}(ϵ^2),$$

where $`u_0(z,t)`$ is given in the Fourier form as

$$u_0=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑\omega W(\omega ,\zeta )\mathrm{exp}\left(\frac{i}{2}\omega ^2\left(_0^zD(z^{})𝑑z^{}\right)+i\omega t\right).$$
(21)
Here $`W(\omega ,\zeta )`$ is a complex Fourier coefficient and $`\zeta =ϵz`$ is the distance to measure the pulse evolution over many map’s periods. By supplying the periodic condition $`u_1(z+1,t)=u_1(z,t)`$ in the Fourier form, the NLS equation (1) can be reduced to the integral evolution model,

$$iW_\zeta \frac{1}{2}D_0\omega ^2W+_{\mathrm{}}^{\mathrm{}}𝑑\omega _1𝑑\omega _2r(\omega _1\omega _2)W(\omega +\omega _1)W(\omega +\omega _2)\overline{W}(\omega +\omega _1+\omega _2)=0,$$
(22)
where

$$r(\omega _1\omega _2)=\frac{1}{4\pi ^2}_0^1𝑑z\mathrm{exp}\left(i\omega _1\omega _2_0^zD(z^{})𝑑z^{}\right).$$

For the two-step dispersion map (12), the integral kernel $`r(x)`$ can be computed explicitly as

$$r(x)=\frac{1}{4\pi ^2}\frac{\mathrm{sin}\left(\frac{mx}{2}\right)}{\frac{mx}{2}}.$$
(23)
It is obvious that the dynamical system (16) and (17) studied in the previous section can be deduced from (22) within the same Gaussian approximation. This correspondence implies that the qualitative results on instability of short Gaussian pulses for $`D_0<0`$ can be reconfirmed within a more systematic theory.

In this section we present numerical results consisting of three subsections. In the first subsection, we construct a numerical solution of the stationary problem identifying optical solitons in the normal regime, when $`D_0<0`$. In the second part, we analyze the linearized problem and locate numerically the linear spectrum in the problem, indicating possible instability of optical solitons. Then, we simulate numerically the non-stationary problem described by (22) and discuss the transformation routes of the unstable dispersion-managed solitons.

### 3.1 Stationary solutions

The periodic-type localized solutions of the NLS equation in the form (3) are equivalent to stationary solutions of (22) in the form,

$$W(\omega ,\zeta )=\mathrm{\Phi }(\omega )e^{i\mu \zeta },$$
(24)
where $`\mathrm{\Phi }(\omega )`$ is the real Fourier coefficient which defines $`\psi (z,t)`$ according to (21). This function satisfies a nonlinear integral boundary-value problem,

$$\left(\mu +\frac{1}{2}D_0\omega ^2\right)\mathrm{\Phi }(\omega )=R(\omega )_{\mathrm{}}^{\mathrm{}}𝑑\omega _1𝑑\omega _2r(\omega _1\omega _2)\mathrm{\Phi }(\omega +\omega _1)\mathrm{\Phi }(\omega +\omega _2)\mathrm{\Phi }(\omega +\omega _1+\omega _2),$$
(25)
where $`\mathrm{\Phi }(\omega )=\mathrm{\Phi }(\omega )`$ (the symmetry condition) and $`lim_\omega \mathrm{}\mathrm{\Phi }(\omega )=0`$ (the boundary condition).

For numerical analysis, we intend to use the Petviashvili’s iteration scheme :

$$\mathrm{\Phi }^{(n)}(\omega )\mathrm{\Phi }^{(n+1)}(\omega )$$

for $`n=0,1,2,\mathrm{}`$. Within this scheme, the right-hand-side $`R(\omega )`$ can be approximated at the $`n`$-th approximation by $`\mathrm{\Phi }^{(n)}(\omega )`$ provided a certain stabilizing factor is introduced for convergence (see (27) and (28) below). However, the numerical scheme breaks down for $`D_0<0`$ due to resonances at $`\omega =\pm \omega _{res}`$, where

$$\omega _{res}=\sqrt{\frac{2\mu }{|D_0|}}.$$
(26)
Indeed for $`\omega =\pm \omega _{res}`$, the left-hand-side of (25) vanishes. \[Here we notice that $`\mu >0`$ for the Gaussian pulse solutions (6) of the NLS model (1).\] In order to avoid resonances in the numerical scheme, we add and subtract a dummy positive dispersion term $`0.5|D_0|\omega ^2\mathrm{\Phi }(\omega )`$ to the left-hand-side of (25). As a result, the scheme converts to the following map,

$$\mathrm{\Phi }^{(n)}(\omega )\mathrm{\Phi }^{(n+1)}(\omega )=S_n^{3/2}\left(\frac{R^{(n)}(\omega )+\frac{1}{2}(|D_0|D_0)\omega ^2\mathrm{\Phi }^{(n)}(\omega )}{\mu +\frac{1}{2}|D_0|\omega ^2}\right),D_0<0,$$
(27)
where $`S_n`$ is the Petviashvili’s stabilizing factor given by

$$S_n=\frac{_{\mathrm{}}^{\mathrm{}}𝑑\omega (\mu +\frac{1}{2}|D_0|\omega ^2)\mathrm{\Phi }^{(n)}(\omega )}{_{\mathrm{}}^{\mathrm{}}𝑑\omega \mathrm{\Phi }^{(n)}(\omega )\left(R^{(n)}(\omega )+\frac{1}{2}(|D_0|D_0)\omega ^2\mathrm{\Phi }^{(n)}(\omega )\right)}.$$
(28)
The factor $`S_n`$ is unity at the stationary solution and serves therefore as an indicator for termination of the iterating procedure. We stop iterations when $`|S_n1|<10^5`$.

To use the map (27), we apply the Simpson’s integration method, reducing complexity due to the symmetry: $`\mathrm{\Phi }(\omega )=\mathrm{\Phi }(\omega )`$. As a starting solution, the profile $`\mathrm{\Phi }(\omega )`$ can be approximated by the Gaussian pulse with parameters corresponding to the periodic solution (13) and (14),

$$\mathrm{\Phi }^{(0)}(\omega )=\sqrt{\pi E\sqrt{2\alpha _s}}\mathrm{exp}\left(\frac{1}{4}\alpha _s\omega ^2\right).$$
(29)
| Number of iterations | $`S_n`$: $`D_0=0.02`$ | $`S_n`$: $`D_0=0.02`$ (I) | $`S_n`$: $`D_0=0.02`$ (II) |
| --- | --- | --- | --- |
| 1 | 0.9897 | 0.9957 | 0.9920 |
| 2 | 0.9971 | 0.9981 | 0.9917 |
| 3 | 0.9994 | 0.9988 | 0.9921 |
| 4 | 0.9998 | 0.9991 | 0.9942 |
| 5 | 0.9999 | 0.9993 | 0.9989 |
| 6 |  | 0.9994 | 1.0069 |
| 7 |  | 0.9995 | 1.0167 |
| 8 |  | 0.9996 | 1.0173 |
| 9 |  | 0.9997 | 0.9642 |
| 10 |  | 0.9997 | 0.8133 |
| 11 |  | 0.9998 | 0.7025 |
| 12 |  | 0.9998 | 0.6833 |
| 13 |  | 0.9999 | 0.6773 |
| 14 |  | 0.9999 | 0.6695 |
| 15 |  | 0.9999 | 0.6601 |
| 16 |  |  | 0.6499 |
| 17 |  |  | 0.6394 |
| 18 |  |  | 0.6292 |
| 19 |  |  | 0.6192 |

Table I. Iterations of the stabilizing factor $`S_n`$ for $`D_0=0.02`$ and $`D_0=0.02`$ (branches I and II).

Table 1 shows iterations for the stabilizing factor $`S_n`$ in the three different cases: (i) $`D_0=0.02`$, (ii) $`D_0=0.02`$ (branch I), and (iii) $`D_0=0.02`$ (branch II). For all the cases, the other parameters are $`\mu =1`$ and $`m=2`$. In the first case, the convergence is monotonic and the profile for stationary soliton $`\mathrm{\Phi }(\omega )`$ is shown in Fig. 5(a). The numerical value for energy of the stationary soliton is shown in Fig. 1(a) by a bullet. In the second case, the iterations converge slowly to the stationary soliton shown in Fig. 5(b). Sometimes, the convergence is accompanied by a single oscillation of $`S_n`$ near unity. The numerical value for the energy is shown in Fig. 2(a) by a bullet at branch I. In the last case, however, the iterations oscillate and finally diverge. Inspection of the profile $`\mathrm{\Phi }^{(n)}(\omega )`$ at a final iteration shows that the iterations change the initial pulse drastically leading to its disappearance. Two conjectures follow from this fact. Either the shorter soliton with larger energy at branch II does not exist as a stationary solution of (25) or it is unstable within the iterational scheme (27). Since the short Gaussian pulse does exist (see Fig. 2), we incline to work along the second conjecture. The iterational scheme (27) is not relevant for the time-evolution problem and rigorous analysis of linearized problem is needed to confirm predictions of the instability of the short stationary pulse.

### 3.2 Linear spectrum

There are several forms of the linear problem associated to the NLS-type equations. We will use the matrix form which was studied in our previous paper  subject to certain simplifications. The matrix form appears upon perturbations of the stationary solutions of (22) as

$$W(\omega ,\zeta )=e^{i\mu \zeta }\left[\mathrm{\Phi }(\omega )+w_1(\omega )e^{i\lambda \zeta }+\overline{w}_2(\omega )e^{i\lambda \zeta }\right].$$
(30)
The vector $`𝐰(\omega )=(w_1,w_2)^T`$ can be shown to satisfy the matrix linear problem,

$$\lambda 𝐰(\omega )=\left(\mu +\frac{1}{2}D_0\omega ^2\right)\sigma _3𝐰(\omega )+_{\mathrm{}}^{\mathrm{}}𝑑\omega _1\left(2K_1(\omega ,\omega _1)\sigma _3+K_2(\omega ,\omega _1)\sigma _2\right)𝐰(\omega _1),$$
(31)
where

$$\sigma _3=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],\sigma _2=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right],$$

and the integral kernels are

$`K_1(\omega ,\omega _1)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega _2r[(\omega \omega _1)(\omega \omega _2)]\mathrm{\Phi }(\omega _2)\mathrm{\Phi }(\omega _1+\omega _2\omega ),`$
$`K_2(\omega ,\omega _1)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega _2r[(\omega _2\omega _1)(\omega _2\omega )]\mathrm{\Phi }(\omega _2)\mathrm{\Phi }(\omega _1\omega _2+\omega ).`$
Using the stationary solution $`\mathrm{\Phi }(\omega )`$ from the previous subsection and implementing the Simpson’s integration method again, we solve the linear problem by using the linear algebra packages built in Matlab 5.2. We identify two types of modes of the linear spectrum: symmetric eigenfunctions, when $`𝐰(\omega )=𝐰(\omega )`$, and anti-symmetric eigenfunctions, when $`𝐰(\omega )=𝐰(\omega )`$.

The linear spectrum for the soliton of Fig. 5(a) is shown in Fig. 6(a,b). It consists of three main parts: (i) continuous spectrum, (ii) neutral (zero) modes, and (iii) internal (oscillatory) modes.

When $`D_0>0`$ (the anomalous regime of the dispersion map), the continuous spectrum is located at the real axis for $`|\lambda |>\mu `$ (see Fig. 6(a,b) where $`\mu =1`$). Indeed, the continuous modes are singular in the Fourier representation, i.e. $`𝐰(\omega )\delta (\omega \mathrm{\Omega })`$. Then, the linear problem (31) has the continuous spectrum at $`\lambda =\pm \lambda _\mathrm{\Omega }`$, where

$$\lambda _\mathrm{\Omega }=\mu +\frac{1}{2}D_0\mathrm{\Omega }^2,$$
(32)
provided the following integral kernels are not singular,

$`\underset{\omega \mathrm{\Omega }}{lim}K_1(\omega ,\mathrm{\Omega }\omega )`$ $`=`$ $`{\displaystyle \frac{1}{4\pi ^2}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega \mathrm{\Phi }^2(\omega ),`$
$`\underset{\omega \mathrm{\Omega }}{lim}K_2(\omega ,\mathrm{\Omega }\omega )`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega r(\omega ^2)\mathrm{\Phi }(\omega +\mathrm{\Omega })\mathrm{\Phi }(\omega \mathrm{\Omega }).`$
The neutral (zero) modes always appear at $`\lambda =0`$ as double degenerate states for both symmetric and asymmetric eigenfunctions. However, the inaccuracy of the numerical method destroys the degeneracy of the zero modes. As a result, the two zero modes may appear either for small real or for small imaginary values of $`\lambda `$. Fig. 7(b) displays two imaginary eigenvalues of order of $`\mathrm{O}(10^2)`$ which appear to be shifted from the origin of $`\lambda `$ due to this numerical effect. Since the zero modes are not of interest from stability point of view, we neglect this effect and leave the scheme without any additional modification.

The internal (oscillatory) modes are located in the gap of the continuous spectrum as $`0<|\lambda |<\mu `$. The set of internal modes may contain different number of eigenvalues. We have shown in previous paper  that the set is empty in the NLS limit (which corresponds to the limit $`\mu 0`$ at fixed $`D_0>0`$ and $`m`$). \[Note that in  we used the Gaussian pulse (29) for approximating $`\mathrm{\Phi }(\omega )`$ while here we substitute the numerical result from Eq. (25).\] Then, we showed that the number of internal modes increases with the map’s strength (if the parameter $`\alpha _s`$ is set to unity, the map’s strength is proportional to $`m`$ and vice versa). In Fig. 6(a,b) for $`\mu =1`$, we identify 14 internal modes for symmetric eigenfunctions and $`12`$ internal modes for anti-symmetric eigenfunctions. Still complex eigenvalues are absent for $`D_0>0`$ which confirms stability of dispersion-managed solitons in the anomalous regime.

The linear spectrum for the soliton of Fig. 5(b) is shown in Fig. 7(a,b). The continuous spectrum is seen to have changed drastically. When $`D_0<0`$ (normal regime of the dispersion map), the continuous spectrum covers the segment $`|\lambda |<\mu `$ twice according to (32). As a result, neutral and internal modes, if any, become embedded in the wave continuum as seen in Fig. 7(a,b). This indicates a resonance of stationary soliton with the linear spectrum in the normal regime of the dispersion map. However, this resonance does not result to any instability of solitons of branch I within the linear theory. We discuss the resonance issue in Section 4. The two imaginary eigenvalues on Fig. 7(b) appear from the origin as artifacts of the numerical scheme as it was explained above.

At last, we would like to construct the linear spectrum for the solitons of branch II. However, the stationary solutions were not identified within the Petviashvili’s numerical method. Therefore, assuming that the solutions still exist, the profile $`\mathrm{\Phi }(k)`$ can only be approximated by the Gaussian pulse (29) as we did in . The linear spectrum in this approximation is shown in Fig. 8(a,b) for $`D_0=0.02`$, $`m=2`$, and $`\mu =1`$. The same type of the continuous spectrum is clearly seen not to possess any gap in the origin. In addition to this, we identify new complex eigenvalues both for symmetric and anti-symmetric eigenfunctions. These complex eigenvalues has relatively large, order of $`\mathrm{O}(10^1)`$, imaginary part and they are associated with the instability of the solitons of branch II. The numerical result for the instability eigenvalues is in a reasonable comparison with the asymptotic predictions which follow from the dynamical model (16) and (17) (see Fig. 3, branch II). The eigenvectors $`𝐰(\omega )`$ for the unstable (imaginary) eigenvalues are shown in Fig. 8(c,d) for the symmetric and anti-symmetric eigenfunctions, respectively. The numerical approximations of the eigenvectors, being inaccurate in details, display clearly that the unstable modes are localized at the intermediate wave frequencies $`\omega `$ of the pulse spectrum (at $`\omega 4.5`$). Thus, the development of the unstable eigenvectors would affect the duration of the soliton pulse in the nonlinear stage as described in the next subsection.

### 3.3 Non-stationary evolution

To confirm the transition scenario, we simulate the non-stationary dynamics of unstable solitons in the integral model (22) by using the central-difference scheme,

$$\frac{V^{(n+1)}(\omega )V^{(n1)}(\omega )}{2\mathrm{\Delta }\zeta }=i_{\mathrm{}}^{\mathrm{}}𝑑\omega _1𝑑\omega _2r(\omega _1\omega _2)e^{iD_0\omega _1\omega _2\zeta _n}V^{(n)}(\omega +\omega _1)V^{(n)}(\omega +\omega _2)\overline{V}^{(n)}(\omega +\omega _1+\omega _2)$$
(33)
where

$$V^{(n)}(\omega )=W(\omega ,\zeta _n)\mathrm{exp}\left(\frac{i}{2}D_0\omega ^2\zeta _n\right)$$

and $`\zeta _n=n\mathrm{\Delta }\zeta `$. An initial iteration can be done within a forward scheme starting with the initial Gaussian pulse (29) with the parameter $`\alpha _s`$ and $`E`$ corresponding to branches I and II at $`\mu =1`$ on Fig. 2(a,b). Evolution of a stable long pulse (branch I) is shown in Fig. 9(a,b), while that of an unstable short pulse (branch II) is shown in Fig. 9(c,d). The stable long Gaussian pulse quickly transits to the stationary pulse given by Eq. (24) (Fig. 5(b)) which propagates later without visible distortions. The second and third peeks in the signal spectrum and profile (Fig. 9(a,b)) appear in complete agreement with the profile of the stationary DM soliton (cf. Fig. 1 from ). Some oscillations along the distance $`\zeta `$ are excited due to the difference between the Gaussian pulse (29) and the exact stationary solution of Eq. (25). These oscillations are small compared to the soliton profile and they do not change the duration of the soliton pulse (see Fig. 9(a)). The nonlinear resonance at $`\omega _{res}`$ ($`\omega _{res}=10`$ for Fig. 9(a)) is not seen to be excited during the signal propagation.

Evolution of the short (unstable) pulse differs drastically from the previous picture. The pulse is being broaden during the evolution, it generates the strong radiation tail and tends to the long (stable) signal which has the first node at $`\omega 4.5`$ (cf. Fig. 9(a) and Fig. 9(c)). This transformation is accompanied by the intermediate oscillations around the soliton’s shape. Thus, we confirm the analytical predictions that the short pulses are linearly unstable and switch into long stable solitons via long-term oscillatory dynamics. The unstable eigenvectors at Fig. 8(c,d) clearly match at Fig. 9(c) with the growing deformations of the localization of the pulse spectrum.

## 4 Discussion: Resonance of dispersion-managed solitons

We have shown that both the Gaussian variational approximation and the integral model prescribe the instability of short nonlinear signals in the normal regime of the dispersion map. This instability broadens the signal’s profile through some intermediate oscillations. The long signals propagate stably then.

We point out that the small-amplitude approximation considered here corresponds to the asymptotic limit $`E/D_00`$ (see branch B at Fig. 1 in ). The stability of the small-amplitude branch for $`D_0<0`$ was previously under discussion in the literature . This discussion can be resolved if one takes into account the normalization condition (15), which was used in the previous works. We have checked that the normalization condition (15) selects the pulse solution along the stable branch I for $`S_{thr}<S<S_{stab}`$, where $`S_{stab}6.76`$. For $`S>S_{stab}`$, it selects solutions along the unstable branch II (see Figs. 2(a,b)). Thus, for the intermediate map strengths (when $`S<S_{stab}`$), the soliton signal propagates stably in the limit of small energies, as reported in . However, if the map strength $`S`$ exceeds the value $`S_{stab}`$, the soliton signal breaks down and switches to a longer signal, as conjectured in .

For both the branches, we observe the resonance appearing between the linear wave spectrum and the stationary signals. This resonance is related to the fact that the origin $`\lambda =0`$ is absorbed in the continuous spectrum of the linear problem. Another way to find the resonance is to construct the linear spectrum for the integral model (22), $`W(\omega ,\zeta )W_0e^{i\mathrm{\Omega }(\omega )\zeta }`$, where

$$\mathrm{\Omega }=\frac{1}{2}D_0\omega ^20$$

and $`D_0<0`$. Therefore, for any $`\mu >0`$ there exists $`\omega _{res}`$ such that $`\mathrm{\Omega }(\omega _{res})=\mu `$. This resonance does not show up in the linear theory since the discrete and continuous spectra are separated. However, in the nonlinear stage, the resonance generally leads to emission of wave packets and soliton’s decay.

Since the transformation or decay of long stable solitons for $`D_0<0`$ have never been observed numerically (nor in our simulations reported in Figs. 9), it is likely that the effective gap in the spectrum still appears in the nonlinear theory. Also the truncated approximation given by the integral model (22) may not be valid for a correct description of the resonance. The latter issues remain open for further analytical consideration.

## Aknowledgements.

The author thanks S. Turitsyn and G. Biondini for bringing the present problem to his attention as well as T. Lakoba and C. Sulem for collaboration and valuable remarks. Fig. 9 was prepared by using computational capacities of Institute of Optics at Rochester, NY.