text:
# Travelling solitons in the parametrically driven nonlinear Schrödinger equation

## I Introduction

The parametrically driven damped nonlinear Schrödinger (NLS) equation,

$$i\psi _t+\psi _{XX}+2|\psi |^2\psi \psi =h\psi ^{}i\gamma \psi ,$$
(1)
was used to model the nonlinear Faraday resonance in a vertically oscillating water tank  and the effect of phase-sensitive parametric amplifiers on solitons in optical fibers . The same equation describes an easy-plane ferromagnet with a combination of a static and hf field in the easy plane , and the planar weakly anisotropic $`XY`$ model . It also serves as a continuum limit for the parametrically driven Frenkel-Kontorova chain (an array of diffusively coupled pendula). The Frenkel-Kontorova system is regarded as a fairly realistic model of a number of physical and biological systems and phenomena, including ladder networks of discrete Josephson junctions, charge-density-wave conductors, crystal dislocations in metals, DNA dynamics and proton conductivity in hydrogen-bonded chains .

The undamped undriven nonlinear Schrödinger equation exhibits soliton solutions which can travel with arbitrary velocity and transport physical characteristics such as mass, momentum and energy. (For the sake of brevity, we are making use of the hydrodynamical interpretation of the equation here.) In equation (1), the second term in the right-hand side accounts for dissipative losses which occur in all physical systems. The dissipation has two visible effects on the soliton: it attenuates its speed and damps its amplitude. The main purpose of the introduction of the pumping (represented by the first term in the right-hand side) is to compensate for these losses. The parametric forcing is well known to be capable of counterbalancing the damping of the quiescent soliton’s amplitude; a natural question now is whether it can sustain its motion with a nonzero velocity.

In fact the existence of travelling solitons is a nontrivial matter even in the absence of damping. The driving term $`h\psi ^{}`$ in Eq.(1) breaks the galilean invariance of the unperturbed nonlinear Schrödinger equation and hence one cannot obtain a moving soliton simply by boosting a static one. However the galilean or Lorentz symmetry is not a prerequisite for the existence of moving nonlinear waves. For example, dissipative systems do not possess any symmetries of this kind but are well known to support stably propagating fronts and pulses whose velocities are determined by parameters of the model. In particular, travelling domain walls arise in the parametrically-driven Ginsburg-Landau equations (where the motion is due to nongradient terms) . As far as solitons in Hamiltonian systems are concerned, the example of dark solitons in the nonlinear Schrödinger equations suggests that they have even greater mobility than dissipative fronts and pulses. Although in this case the galilean invariance is broken by the presence of the nonvanishing background, the dark solitons can propagate with arbitrary speeds bounded only by the velocity of sound waves .

A number of nonstationary regimes were reported in the water tank experiments, including the formation of oscillating soliton pairs , but no steadily moving solitons were detected so far. On the other hand, numerical simulations of the undamped parametrically driven NLS equation \[Eq.(1) with $`\gamma =0`$\] did exhibit travelling localised objects . It has remained an open question whether these moving objects preserve their speed and amplitude, or attenuate and decay slowly due to the emission of the second-harmonic radiation. The aim of the present paper is to study the existence of steadily propagating solitons, and examine their stability. Here we are confining ourselves to the undamped situation relegating the analysis of the effect of damping to future publications.

In addition to their role in transport phenomena, stably moving solitons are also of interest as alternative attractors which may compete with (static or oscillating) nonpropagating solutions. We will demonstrate that stable travelling solitons do exist in the parametrically driven nonlinear Schrödinger equation. Moreover, there are parameter ranges where moving solitons are stable whereas their quiescent counterparts are not. Unstable solitons are not meaningless either; they arise as long-lived transients and intermediate states in spatio-temporal chaotic regimes. In this paper we will identify oscillatory and translational instabilities of travelling solitons and simulate their nonlinear evolution near the transition curves.

The structure of the paper is as follows. In section II we derive an apriori bound for the existence domain of travelling solitons and introduce the linearised eigenvalue problem for their stability analysis. We also discuss some general properties of eigenvalues and eigenfunctions and formulate a simple criterion for the onset of the nonoscillatory instability: $`P/V=0`$, where $`V`$ is the velocity of the steadily moving soliton, and $`P`$ the associated momentum.

In section III we present several explicit quiescent ($`V=0`$) solutions, including a stationary complex of two solitons, and then derive the necessary condition for a static solution to be continuable to nonzero velocities. This condition is that the motionless solution should either not have any “free” parameters apart from the translational shift, or, if there is an additional parameter $`z`$, the equation $`P/z=0`$ should be satisfied. Here $`P`$ is the momentum of the motionless localised solution (which, contrary to one’s mechanical intuition, is not necessarily equal to zero). There are three static solutions satisfying the above condition, two of which being the well-known constant-phase $`\psi _+`$ and $`\psi _{}`$ solitons, respectively, while the third solution looks like a pulse with a bell-shaped modulus and twisted phase.

The most important results of this work are contained in section IV where we report on the numerical continuation of various branches of solutions and their stability analysis. In agreement with the analytical predictions of the preceding section, we find that each of the above static solutions admits the continuation to nonzero $`V`$. The stability properties of travelling solitons result from an intricate interplay of two types of instability, the oscillatory and translational instability. In accordance with the conclusions of section II, the numerical analysis of the linearised eigenvalues shows that the transition curves of the translational instability satisfy $`P/V=0`$. One interesting conclusion of the stability analysis is that although quiescent solitons are unstable for driving strengths larger than $`h=0.064`$, there are stable moving solitons for any $`0h1`$. We discuss in detail the soliton’s transformation as it is continued in $`V`$, paying special attention to the dynamics of the associated linearised eigenvalues on the complex plane. Two different scenarios of the transformation are identified, one occurring for small $`h`$ and the other one for larger driving strength, and we also describe an interesting cross-over from one to another.

Section V is devoted to the direct numerical simulations of the full time-dependent nonlinear Schrödinger equation. We show that the evolution of both types of the soliton instability leads, as $`t\mathrm{}`$, to the same asymptotic attractors. Finally, section VI contains concluding remarks and outlines some open problems.

## II Steadily travelling waves: existence and stability

### A Existence domain and integrals of motion

We will confine ourselves to localised travelling waves of the simplest form, $`\psi (X,t)=\psi (XVt)`$. Transforming to the moving frame, these correspond to time-independent soliton solutions of the equation

$$i\psi _tiV\psi _x+\psi _{xx}+2|\psi |^2\psi \psi =h\psi ^{},$$
(2)
where $`x=XVt`$. We will search for these static solutions by solving an ordinary differential equation

$$iV\psi _x+\psi _{xx}+2|\psi |^2\psi \psi =h\psi ^{}$$
(3)
under the vanishing boundary conditions $`|\psi (x)|0`$ as $`|x|\mathrm{}`$. Here $`h`$ is always taken positive; negative $`h`$’s can be recovered by the phase transformation $`\psi i\psi `$.

If $`\psi (x)`$ is a solution, then so is $`\psi (x)`$. Next, it is straightforward to notice that if the function $`\psi (x)`$ describes a soliton travelling with the velocity $`V`$, the function $`\psi ^{}`$(x) describes a soliton moving with the velocity $`V`$. We will try to restrict ourselves to positive $`V`$’s wherever possible; we will only present negative velocities where this may help visualising how different branches of solutions are connected. Since the soliton moving with the velocity $`V`$ is given by $`\psi (x)`$, the above observation tells us that either $`\psi ^{}(x)=\pm \psi (x)`$ (that is, one of the real and imaginary part of the solution is even and the other one odd), or there are two solutions associated with the same $`V`$. (Here we are not making any difference between solutions which are different just in the overal sign.) In the latter case the solutions will not exhibit the $`\psi ^{}(x)=\pm \psi (x)`$ symmetry.

Next, it is easy to show that solitons cannot travel faster than a certain limit speed. Indeed, as $`|x|\mathrm{}`$, the soliton’s asymptotic tail decays as $`\psi (x)e^{\kappa x}`$, where

$$2\kappa ^2=2V^2\pm \sqrt{\left(2V^2\right)^2+4(h^21)}.$$
(4)
Large driving strengths, $`h>1`$, are of little interest to us as in this case the zero background, $`\psi (x)=0`$, is unstable with respect to continuous spectrum waves . Therefore we are not going to discuss this case here. In the complementary region $`h<1`$, the complex structure of $`\kappa `$ depends on the value of the velocity. When $`V^2<22\sqrt{1h^2}`$, there are four real exponents; for $`22\sqrt{1h^2}<V^2<c^2,`$ where

$$c=\sqrt{2+2\sqrt{1h^2}},$$
(5)
we have a quadruplet of complex $`\kappa `$’s, and finally, for $`V^2>c^2`$ all four exponents are imaginary. Consequently, there can be no exponentially localised solitons travelling faster than $`c`$. Physically, $`c`$ represents the minimum phase velocity of linear waves governed by Eq.(1), and our condition $`V<c`$ is essentially an exclusion principle ruling out a resonance between solitons and linear waves.

In the undamped case the parametrically driven NLS equation (2) conserves the momentum,

$$P=\frac{i}{2}\left(\psi _x^{}\psi \psi _x\psi ^{}\right)𝑑x,$$
(6)
and energy:

$$E=\mathrm{Re}\left(|\psi _x|^2+|\psi |^2|\psi |^4+h\psi ^2\right)𝑑x.$$
(7)
Noting that stationary solutions satisfy

$`|\psi _x|^2=|\psi |^2|\psi |^4+h\mathrm{Re}\psi ^2,`$

Eq.(7) can be rewritten as

$$E=2\mathrm{Re}\left(|\psi |^2|\psi |^4+h\psi ^2\right)𝑑x.$$
(8)
Containing no derivatives, this formula for energy of stationary solutions has obvious advantages for the numerical implementation.

### B Linearised eigenvalue problem

In this paper we solve the equation (3) numerically and examine the stability of the resulting solutions by studying the associated eigenvalue problem. This eigenvalue problem arises by assuming a small perturbation of the form

$`\delta \psi (x,t)=y(x)e^{\lambda t},y(x)=\delta u(x)+i\delta v(x).`$

Substituting into (2) gives

$$Y=\lambda JY,$$
(9)
where the hermitean operator $``$ has the form

$`=I(_x^2+1)+VJ_x+`$ (10)
$`+\left(\begin{array}{cc}h6u^22v^2\hfill & \hfill 4uv\\ 4uv\hfill & \hfill h6v^22u^2\end{array}\right);`$ (13)
the matrix $`J`$ is given by

$$J=\left(\begin{array}{cc}0\hfill & \hfill 1\\ 1\hfill & \hfill 0\end{array}\right),$$
(14)
and the column-vector $`Y(x)=(\mathrm{Re}y,\mathrm{Im}y)^T=(\delta u,\delta v)^T`$. In Eq.(13) $`I`$ is the identity matrix, and we have decomposed the stationary solution as $`\psi (x)=u(x)+iv(x)`$.

For symmetric solutions satisfying $`\psi ^{}(x)=\pm \psi (x)`$ eigenvalues will always come in ($`\lambda ,\lambda `$)-pairs. This follows from the fact that for these solutions changing $`xx`$ in the operator (13) amounts to changing the sign of its off-diagonal elements, and hence if $`(\delta u(x),\delta v(x))^T`$ is an eigenfunction associated with an eigenvalue $`\lambda `$, the column $`(\delta u(x),\delta v(x))^T`$ will serve as an eigenfunction associated with an eigenvalue $`\lambda `$. As far as the zero eigenvalue is concerned, it will have a twin with the eigenfunction $`(\delta u(x),\delta v(x))^T`$ unless its eigenfunction $`y=\delta u+i\delta v`$ satisfies the symmetry $`y^{}(x)=e^{i\phi }y(x)`$, where $`\phi =const`$.

To complete the discussion of the spectrum structure, we mention that there are two branches of the continuous spectrum lying on the imaginary axis of $`\lambda `$: $`\lambda =i\omega _{1,2}(k)`$, where

$`\omega _{1,2}(k)=Vk\pm \sqrt{(k^2+1)^2h^2},`$

and $`\mathrm{}<k<\mathrm{}`$. (We are still assuming $`h<1`$). In the region $`V^2<c^2`$ which is of interest to us, the continuous spectrum has a gap: $`\omega _1(k)>\omega _0`$, $`\omega _2(k)<\omega _0`$, where $`\omega _0>0`$. This gap can harbour discrete eigenvalues representing stable oscillation modes.

### C Translational (nonoscillatory) instabilities

The aim of this subsection is to demonstrate that a pair of pure imaginary eigenvalues can collide at $`\lambda =0`$ and move onto the real axis only at the point where $`P/V=0`$. This criterion is known in the context of dark solitons of the undriven nonlinear Schrödinger equations; see ,. Here we simply adapt the proof given in  to the case of the equation with the parametric forcing. An important assumption that we make here, is that the solution whose stability is being examined, does not have any free parameters apart from the trivial translation parameter, $`x_0`$.

First of all we need to make a remark on the integrable case, $`h=0`$. In this case solutions of the ODE (3) can be obtained from a quiescent soliton of Eq.(1) by a Galilei transformation:

$$\psi (x)=e^{i(V/2)x}A\mathrm{sech}Ax,$$
(15)
where $`A=\sqrt{1V^2/4}`$. For $`h=0`$ and any $`|V|<2`$, the linearised operator $``$ has four zero eigenvalues associated with two eigenvectors. One of these eigenvectors originates from the translation symmetry and the other one results from the phase invariance of Eq.(3). The term $`h\psi ^{}`$ breaks the phase invariance and hence as $`h`$ is increased from zero, one pair of eigenvalues $`(\lambda ,\lambda )`$ moves away from the origin on the complex plane. As $`h`$ and $`V`$ are further varied, a pair of eigenvalues may return to the origin. If the solution of Eq.(3) at the point of their return is a member of a family parametrised by two free parameters, we will have, again, four zero eigenvalues with two eigenfunctions. (The eigenfunctions are simply derivatives of the solution with respect to the free parameters.) Our analysis will not be applicable in this case, and the equality $`P/V=0`$ does not have to be valid at the return point. (We will come across this type of a situation in section IV D below.) However, a more common situation is when the solution at the return point is a member of a one-parameter family. We will show that in this case the relation $`P/V=0`$ does have to be in place.

Let us denote $`V_c`$ the velocity for which the eigenvalue of the operator (9)-(13) vanishes. We can develop the solution $`\psi (V;x)`$ in powers of $`ϵ=VV_c`$:

$`\psi (V;x)=\psi _0(x)+ϵ\psi _1(x)+ϵ^2\psi _2(x)+\mathrm{},`$

where $`\psi _0=\psi (V_c;x)`$. Accordingly, the operator $``$ expands as $`=_0+ϵ_1+ϵ^2_2+\mathrm{}`$. If the eigenvalue $`\lambda `$ moves from imaginary to the real axis, it is natural to assume that it admits an expansion of the form

$$\lambda =ϵ^{1/2}\lambda _1+ϵ^{3/2}\lambda _3+ϵ^{5/2}\lambda _5+\mathrm{}.$$
(16)
The associated eigenfunction is then developed as

$$Y(x)=Y_0(x)+ϵ^{1/2}Y_1(x)+ϵY_2(x)+\mathrm{}.$$
(17)
When $`ϵ=0`$, we have $`_0Y_0=0`$, i.e. $`Y_0`$ is a null eigenvector at the bifurcation point $`V=V_c`$. Since we have assumed that $`\psi (V_c,x)`$ is a member of a one-parameter family of solutions, the operator $`_0`$ has only one null eigenvector, and we have to identify $`Y_0=\mathrm{\Psi }_0^{}(x)`$. Here $`\mathrm{\Psi }_0`$ is a column-vector formed by the real and imaginary part of the soliton $`\psi _0`$: $`\mathrm{\Psi }_0=(u_0,v_0)^T`$. The prime indicates differentiation with respect to $`x`$.

Next, setting the coefficient of $`ϵ^{1/2}`$ to zero yields

$`_0Y_1=\lambda _1JY_0.`$

Comparing this to the equation

$`_0{\displaystyle \frac{\mathrm{\Psi }}{V}}|_{V=V_c}=J\mathrm{\Psi }_0^{},`$

which arises from the differentiation of Eq.(3) with respect to $`V`$, we get

$`Y_1(x)=\lambda _1{\displaystyle \frac{\mathrm{\Psi }}{V}}|_{V=V_c}.`$

(In the above equations $`\mathrm{\Psi }=(u,v)^T`$.) The coefficient of $`ϵ^1`$ produces

$`_0Y_2=\lambda _1JY_1_1Y_0,`$

which has bounded solutions if the right-hand side is orthogonal to the null eigenvector of $`_0`$:

$$\lambda _1Y_0JY_1𝑑xY_0_1Y_0𝑑x=0.$$
(18)
The second term in equation (18) is readily shown to vanish — one only needs to expand the identity $`\mathrm{\Psi }^{}=0`$ in $`ϵ`$. (The coefficient of $`ϵ^1`$ gives $`_1\mathrm{\Psi }_0^{}=_0\mathrm{\Psi }_1^{}`$. Taking the scalar product with $`\mathrm{\Psi }_0^{}`$ yields the required $`\mathrm{\Psi }_0^{}_1\mathrm{\Psi }_0^{}𝑑x=0`$.) On the other hand, the first term in Eq.(18) is equal to $`(\lambda _1^2/2)P/V`$. Consequently, Eq.(18) gives either $`P/V=0`$ or $`\lambda _1=0`$. If we assume that $`\lambda _1=0`$, we will not be able to conclude that $`P/V=0`$ at this order of the expansion. However, the order $`ϵ^2`$ will then give us $`\lambda _2^2P/V=0`$, which implies either $`P/V=0`$ or $`\lambda _2=0`$. Proceeding by a similar token we will eventually arrive at the equation $`P/V=0`$ at some order $`ϵ^n`$ where $`n`$ is such that $`\lambda _n0`$. (Alternatively, we will have to conclude that all $`\lambda _n=0`$ and hence we are dealing with a symmetry eigenvalue which is equal to zero for all $`V`$.)

Thus a pair of real or pure imaginary eigenvalues of the same magnitude and opposite sign, can only collide for the value of $`V`$ which satisfies $`P/V=0`$. Here we wish to re-emphasise that we have obtained this conclusion under the assumption that the geometric multiplicity of the zero eigenvalue is not increased at the point of collision. A simple example when this assumption is not valid, is furnished by the case $`h=0`$. In this case the momentum corresponding to the soliton (15) is given by $`P=V\sqrt{1V^2/4}`$. Although $`P`$ has a maximum for $`V=\sqrt{2}`$, the stability properties of the undriven soliton do not change at this point. The reason is that for each $`V`$ the operator $``$ has two null eigenvectors in this case, and hence we cannot make the identification $`Y_0=\mathrm{\Psi }_0^{}`$. (Instead, $`Y_0`$ will be a linear combination of two zero modes.) Consequently, the above proof becomes invalid in this case.

Finally, one can easily check that the above result does not really depend on how the eigenvalue $`\lambda `$ expands in powers of $`ϵ`$. We assumed that the expansion (16) starts with terms of order $`ϵ^{1/2}`$. This assumption is natural and supported by the numerical evidence; however, even if we had postulated the expansion starting with terms of order $`ϵ^{1/4}`$, $`ϵ^{1/3}`$ or say, $`ϵ`$, we would have still arrived at the same necessary condition for the zero crossing: $`P/V=0`$.

## III Quiescent solutions ($`V=0`$)

### A The “twist” soliton

In order to continue in $`V`$ it is useful to have some starting solutions for $`V=0`$. Two such stationary solutions are given by

$`\psi _+(x)=A_+\mathrm{sech}(A_+x),`$ (20)
$`\psi _{}(x)=iA_{}\mathrm{sech}(A_{}x),`$ (21)
where $`A_\pm ^2=1\pm h`$. The soliton $`\psi _{}`$ is unstable with respect to a nonoscillatory mode for all $`h`$ . The $`\psi _+`$ is stable for $`h<h_0=0.063596`$ but developes an oscillatory instability as $`h`$ is increased beyond $`h_0`$ .

In this section we will produce several more explicit solutions of the undamped, parametrically driven NLS equation (2). Writing $`\psi =u+iv`$, the stationary equation (3) transforms into the system

$`u_{xx}uhu+2u(u^2+v^2)=0,`$ (22)
$`v_{xx}v+hv+2v(u^2+v^2)=0.`$ (23)
We can try to find explicit solutions of this system by imposing some plausible reductions, for instance by identifying $`u^2+v^2`$ with a function of $`u`$ and its derivatives: $`u_x`$, $`u_{xx}`$ and so on. In this case Eq.(22) is an equation for $`u`$ only, while Eq.(23) becomes a linear equation with variable coefficients. The simplest choice $`u^2+v^2=Cu^2`$ leads to the $`\psi _+`$ and $`\psi _{}`$ solitons (III A). Another simple possibility is to require that

$$u^2+v^2=Cu,C=\mathrm{const};$$
(24)
this converts the first equation into the stationary KdV with the well-known localised solution

$`u={\displaystyle \frac{3}{4}}{\displaystyle \frac{1+h}{C}}\mathrm{sech}^2\xi ,\xi ={\displaystyle \frac{\sqrt{1+h}}{2}}x.`$

The second equation has now the form of an eigenvalue problem for the Pöschl-Teller potential:

$$(_\xi ^2+16\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{h}^2\xi )v=Ev,$$
(25)
where $`E=1+4(h1)/(h+1)`$. The operator (25) has two localized eigenfunctions, $`v_0=\alpha \mathrm{sech}^2\xi `$ associated with an eigenvalue $`E_0=3`$, and $`v_1=\alpha \mathrm{sech}\xi \mathrm{tanh}\xi `$ with an eigenvalue $`E_1=0`$. The first eigenfunction can satisfy the constraint (24) for no $`\alpha `$ while the second one will satisfy it if we set $`C^2=\alpha ^2=\frac{3}{4}(1+h)`$.

Noticing that $`E=0`$ corresponds to $`h=\frac{3}{5}`$, we conclude that for $`h=\frac{3}{5}`$ there is an explicit solution of the form

$`u_T=\sqrt{{\displaystyle \frac{6}{5}}}\mathrm{sech}^2\xi ;v_T=\pm \sqrt{{\displaystyle \frac{6}{5}}}\mathrm{sech}\xi \mathrm{tanh}\xi ,`$ (26)
where $`\xi =\sqrt{\frac{2}{5}}x`$. Similarly to the solitons $`\psi _+`$ and $`\psi _{}`$, the modulus of the above solution is bell-shaped, but, unlike the constant phase of the $`\psi _+`$ and $`\psi _{}`$ solitons, the phase of this solution varies. (In the case of the positive sign in (26) the phase grows from $`\pi /2`$ at $`x=\mathrm{}`$ to $`\pi /2`$ at $`x=\mathrm{}`$.) The solution looks like a pulse twisted by $`180^{}`$ in the $`(u,v)`$-plane. For this reason we will be referring to Eq.(26) as the “twist” soliton.

### B The twist soliton as a bound state

The system (22)-(23) appeared previously as a stationary system governing light pulses in a birefringent optical fiber. Using Hirota’s approach, Tratnik and Sipe  have obtained the following exact solution to eqs.(3),(22)-(23):

$`\psi =\psi (z;x)=u+iv,`$ (28)
$`u=2A_+e^{\theta _2}D^1\left(1+e^{2(\theta _1\beta )}\right),`$ (29)
$`v=2A_{}e^{\theta _1}D^1\left(1e^{2(\theta _2\beta )}\right),`$ (30)
where

Fig.1 Solution (III B) for various $`z`$. (a): $`z=10`$; (b): $`z=\zeta `$ with $`\zeta `$ as in (33) (the twist soliton); (c): $`z=10`$. Solid curve: real part; dashed line: imaginary part.

$`D=1+e^{2\theta _1}+e^{2\theta _2}+e^{2(\theta _1+\theta _22\beta )},`$
$`\theta _1=A_{}(xz),\theta _2=A_+(x+z);`$
the constant

$$\beta =\frac{1}{2}\mathrm{ln}\left(\frac{A_++A_{}}{A_+A_{}}\right)>0,$$
(31)
the amplitudes $`A_\pm `$ are as in (III A): $`A_\pm =\sqrt{1\pm h}`$, and $`z`$ is a real parameter which can take arbitrary values. The solution (III B) with $`z=10`$ and $`10`$ is plotted in Fig.1, (a) and (c). As is clear from the figure, for sufficiently large $`|z|`$ the solution represents a complex of two solitons, $`\psi _+`$ and $`\psi _{}`$, with the separation approximately equal to $`2|z|`$. It is perhaps worth emphasizing here that the parameter $`z`$ has nothing to do with shifting the solution as a whole, $`xxx_0`$, which is possible due to the translation invariance of Eqs.(1)-(3). (This overall shift parameter, $`x_0`$, is disregarded in equations (26)-(III B) and in the remainder of this text). The parameter $`z`$ is nontrivial in the sense that the shape of the solution depends on $`z`$.

It is straightforward to verify that our “twist” solution is a particular case of Eq.(III B) with $`h=\frac{3}{5}`$. Indeed, choosing $`z=\frac{1}{8}\sqrt{\frac{5}{2}}\mathrm{ln}3`$ in Eqs.(III B), we obtain the soliton (26) (the one with the negative sign), centered at $`x=3z`$. Since the twist is a symmetric solution with sech-shaped modulus, it would be difficult to interpret it as a bound state of the $`\psi _+`$ and $`\psi _{}`$ without embedding it into a broader family of solitonic complexes. In fact, a similar symmetric solution exists for any $`h`$, not only for $`h=\frac{3}{5}`$. To see this, we notice a simple relation between two solutions of the form (III B) — one with the parameter value $`z=\zeta +\xi `$ and the other one with $`z=\zeta \xi `$:

$$\psi (\zeta +\xi ;\eta y)=\psi ^{}(\zeta \xi ;\eta +y).$$
(32)
Here $`\zeta `$ and $`\eta `$ are defined by the driving strength $`h`$:

$$\zeta =\frac{\beta }{2}\left(\frac{1}{A_{}}\frac{1}{A_+}\right)<0$$
(33)
and

$$\eta =\frac{\beta }{2}\left(\frac{1}{A_{}}+\frac{1}{A_+}\right),$$
(34)
while $`\xi `$ and $`y`$ can take arbitrary values. The relation (32) implies that the solution (III B) with $`z=\zeta `$ is symmetric about the point $`x=\eta `$:

$$\psi (\zeta ;\eta y)=\psi ^{}(\zeta ;\eta +y).$$
(35)
That is, the real part of this solution is even and imaginary part odd with respect to $`x=\eta `$:

$`u(\eta y)=u(\eta +y),v(\eta y)=v(\eta +y).`$

(See Fig.1 (b).) This particular representative of the family (III B) will play a special role in what follows. For the ease of reference we are retaining the name “twist” for this symmetric solution — for all $`h`$.

### C The moving soliton bifurcation

Suppose the equation (3) has a one-parameter family of quiescent solutions $`\psi (z;x)`$. Here $`z`$ can be any nontrivial parameter; the only requirement is that $`z`$ should not be just an overall shift in $`x`$. One such family is given by Eq.(III B) and there can also be other families for which $`\psi `$ is not available explicitly. We will show in this section that in order for a solution with some $`z=z_0`$ to be continuable to nonzero $`V`$, the corresponding momentum integral should satisfy

$$\frac{P}{z}|_{z=z_0}=0.$$
(36)
Let us assume that Eq.(3) with $`V0`$ has a solution $`\psi (x)`$, and that this solution is an analytic function of $`V`$ in some neighbourhood of $`V=0`$. Then we can expand it in the Taylor series

$$\psi (x)=\psi _0(x)+V\psi _1(x)+V^2\psi _2(x)+\mathrm{},$$
(37)
where $`\psi _0(x)=\psi (z_0;x)=u_0+iv_0`$ is some representative of the family of “motionless” solutions $`\psi (z;x)`$ with the parameter value $`z_0`$. Substituting (37) into (3) and equating coefficients of like powers of $`V`$, we get, at the order $`V^1`$:

$$\left(\begin{array}{c}u_1\\ v_1\end{array}\right)=J_x\left(\begin{array}{c}u_0\\ v_0\end{array}\right).$$
(38)
Here $`u_1+iv_1=\psi _1`$ and the operator $``$ is given by Eq.(13). The equation (38) is solvable in the class of square integrable functions if the vector in the right-hand side is orthogonal to all homogeneous solutions, i.e. to all null eigenvectors of the operator $``$. Since there is a family of “motionless” solutions parametrized by $`z`$ and by an arbitrary spatial shift $`x_0`$ (which we have disregarded so far), the operator $``$ has two zero modes. One is the translation mode $`_x\psi _0=_x(u_0+iv_0)`$; the corresponding solvability condition is trivially satisfied:

$`{\displaystyle _x(u_0,v_0)J_x\left(\begin{array}{c}u_0\\ v_0\end{array}\right)dx}=0.`$

The other zero mode is given by the derivative $`_z\psi _0=_zu_0+i_zv_0`$. The associated solvability condition reads

$`0={\displaystyle (_zu_0,_zv_0)J_x\left(\begin{array}{c}u_0\\ v_0\end{array}\right)dx}={\displaystyle \frac{1}{2}}{\displaystyle \frac{P}{z}},`$

where $`P`$ is the momentum integral (6). Consequently, a solution with nonzero $`V`$ can only detach from the $`V=0`$ branch at the point where $`P/z=0`$.

Coming back to our explicit solutions, the $`\psi _+`$ and $`\psi _{}`$ solitons do not have any free parameters apart from the trivial position shift. Consequently, both solutions are continuable to nonzero $`V`$. Next, we have a family of solitonic complexes (III B) with a nontrivial parameter $`z`$. As one can easily check, the momentum of the complex (III B) as a function of $`z`$ has a single minimum for some finite $`z=z_0`$ and tends to zero as $`z\pm \mathrm{}`$. To find $`z_0`$, we notice that the relation (32) implies

$`P(\zeta +\xi )=P(\zeta \xi ).`$

This means that the function $`P(z)`$ is even with respect to the point $`z=\zeta `$ and therefore, $`\zeta `$ is the point of the minimum: $`z_0=\zeta `$. Thus, the only representative of the family of the two-soliton complexes (III B) that can be continued to nonzero $`V`$, is our twist soliton, $`\psi (\zeta ;x)`$. (To be more precise, there are two twist solutions, one with positive and the other one with negative momentum. This is related to the fact that when $`V=0`$, we can generate new solutions to the system (22)-(23) by changing the sign of just one component $`u`$ or $`v`$.)

## IV Bifurcation diagram

We used a predictor-corrector continuation algorithm with a fourth-order Newtonian solver to continue solutions of equation (3) in $`V`$. Since derivatives of the momentum integral (6) determine stability and branching properties of solutions, the momentum was our natural choice for the bifurcation measure. Eq.(3) was solved under the vanishing boundary conditions $`\psi \left(\pm L/2\right)=0`$. We used $`L=200`$ (except the cases where we had to extend the interval to account for slow decay of solutions) and the discretisation stepsize $`\mathrm{\Delta }x=0.005`$. The eigenvalue problem (9) was solved on the interval $`(50,50)`$. Here we utilised the Fourier method, typically with 600 harmonics.

### A The travelling $`\psi _{}`$ soliton

We start our description with the branch departing from the quiescent soliton $`\psi _{}`$. For every $`h`$ this branch continues all way to $`V=c`$, where $`c`$ is the minimum phase velocity of linear waves given by Eq.(5). (See Fig.2). As $`Vc`$, the decay rate of $`\psi (x)`$ decreases and the soliton merges with the zero solution, with the momentum $`P`$ tending to zero. We plotted $`P(V)`$ for $`h=0.1,0.3,0.5,0.7`$ and $`0.9`$ in Fig.2. For technical reasons we could not connect the curve $`P(V)`$ to zero although we were able to approach the value $`V=c`$ as close as the fourth digit after the decimal point. (The problem is that since the decay rate of the solution decreases, one has to increase the length of the integration interval — and this cannot be done indefinitely.) The only curve which is connected to zero in Fig.2, is the one for the undriven case $`h=0`$. In this case we enjoy an explicit solution (15) with the momentum $`P=V\sqrt{1V^2/4}`$.

Fig.2 The momentum of the $`\psi _+`$ and $`\psi _{}`$ solitons as a function of their velocities. Thick respectively thin lines depict stable respectively unstable solutions. Note that as $`h0`$, the stability domain of the $`\psi _{}`$ tends to $`\sqrt{2}V2`$ and that of the $`\psi _+`$ to $`0V\sqrt{2}`$. The whole of the $`h=0`$ branch is stable.

For each $`h`$ the momentum of the soliton has a single maximum on this branch, at $`V=V_c`$ (Fig.2). To the left of $`V_c`$ the linearised operator (9)-(13) has a pair of real eigenvalues $`\pm \lambda `$ and consequently, the soliton $`\psi _{}`$, which is well known to be unstable for $`V=0`$ , remains unstable for small nonzero velocities. As $`V`$ approaches $`V_c`$, the two eigenvalues converge at the origin on the complex plane, with the associated eigenfunctions tending to the translation mode $`\mathrm{\Psi }_0^{}(x)`$. Increasing $`V`$ past $`V_c`$, the eigenvalues move onto the imaginary axis and hence the $`\psi _{}`$ soliton becomes stable for sufficiently large velocities (where $`P/V<0`$). This change of stability properties is in exact agreement with the scenario described in section II C.

### B The travelling $`\psi _+`$ soliton; $`h0.25`$

In agreement with conclusions of section III C, we have found that the soliton $`\psi _+`$ is also continuable to nonzero velocities. Unlike the $`\psi _{}`$-branch where the final product of the continuation process does not depend on the value of $`h`$, the transformation scenario of the soliton $`\psi _+`$ is different for different $`h`$. For $`h<0.28`$ the fate of the soliton $`\psi _+`$ is similar to that of the $`\psi _{}`$. As $`Vc`$, the soliton develops oscillations on its tails; the width of the resulting oscillatory “wavepacket” grows and the amplitude decreases, until the solution becomes equal to zero everywhere. The momentum $`P(V)`$ tends to zero as $`Vc`$ and has a single maximum at some $`V=V_c`$. Stability properties of the $`\psi _+`$ soliton depend on whether $`h`$ is smaller than 0.064, lies between 0.064 and 0.25, or is greater than 0.25.

Let, first, $`h0.25`$. As $`V`$ is increased to $`V_c`$, a pair of imaginary eigenvalues $`\pm \lambda `$ collides at the origin on the complex plane and moves onto the real axis. This is in contrast to the $`\psi _{}`$ soliton, where two real eigenvalues converged at the origin as $`V`$ was increased. For $`V>V_c`$ (i.e. where the slope $`P/V`$ is negative), the $`\psi _+`$ soliton becomes unstable. (Note that the $`\psi _{}`$ was unstable for positive slopes $`P/V`$.) As in the $`\psi _{}`$ case, this is the translational, or nonoscillatory, instability. Since, as we checked, the eigenfunctions associated with the colliding imaginary eigenvalues tend to $`\mathrm{\Psi }_0^{}(x)`$ as $`VV_c`$, this instability is of the type analysed in section II C. The numerically detected instability boundary, defined by the condition $`P/V=0`$, is in exact correspondence with our analytical predictions.

In addition, for some $`V`$ the $`\psi _+`$ soliton exhibits the oscillatory instability. (The oscillatory instability sets in when four eigenvalues collide, pairwise, on the imaginary axis and emerge into the complex plane.) For $`h0.064`$ the oscillatory instability does not arise for any $`V`$ and hence the whole range $`0<V<V_c`$ is stable. (See the $`h=0.05`$ curve in Fig.2). For $`0.064<h0.25`$ the oscillatory instability occurs for $`V`$ zero and small, but as $`V`$ is increased, both imaginary and real parts of the “unstable” eigenvalues decay, with the real parts decaying faster. Eventually the eigenvalues $`\pm \lambda ,\pm \lambda ^{}`$ converge, pairwise, on the imaginary axis and the soliton stabilizes. Increasing $`V`$ still further, two of the resulting imaginary eigenvalues, $`\lambda `$ and $`\lambda `$, start approaching each other. At $`V=V_c`$ where $`P/V=0`$, they collide and move onto the real axis. The soliton looses its stability once again — this time to a nonoscillatory mode. This scenario is exemplified by the curves $`h=0.1`$ and $`h=0.15`$ in Fig.2.

An interesting phenomenon occurs in the undriven case, $`h=0`$. As we mentioned, in this case the stationary solution of Eq.(2) is available in explicit form, Eq.(15). We also noted that the associated momentum $`P=V\sqrt{1V^2/4}`$ has a maximum for finite $`V=\sqrt{2}`$. On the other hand, in any small-$`h`$ neighbourhood of this solution there are solitons $`\psi _+`$ and $`\psi _{}`$, of which $`\psi _+`$ is stable only for $`P/V>0`$ (i.e. for $`V<\sqrt{2}`$) and $`\psi _{}`$ is stable only for $`P/V<0`$ (i.e. for $`V>\sqrt{2}`$). How can this be reconciled with the fact that the solution (15) is stable for all $`V`$ for which it exists ($`V2`$)?

The answer is related to the behaviour of the zero eigenvalue associated with the phase invariance of Eq.(2) with $`h=0`$. As $`h`$ increases from zero, this eigenvalue moves away from the origin. On the $`\psi _+`$ branch, it passes onto the imaginary axis, while on the $`\psi _{}`$ branch, the eigenvalue moves onto the real axis instead.

### C The travelling $`\psi _+`$ soliton; $`h>0.25`$

Let now $`0.25<h<0.28`$, and assume we are moving along the $`\psi _+`$ branch in the direction of larger $`V`$. For small $`V`$ we have a quadruplet of complex eigenvalues $`\pm \lambda `$, $`\pm \lambda ^{}`$ implying the oscillatory instability. As $`V`$ is increased, both imaginary and real parts decay — as in the $`h0.25`$ case. However, this time imaginary parts decay faster than the real parts, and the two pairs of eigenvalues converge on the real axis. For velocities above this point the oscillatory instability is replaced by the translational instability. As $`V`$ is increased further, one pair of the newly born real eigenvalues grows in absolute value whereas the other pair decreases in magnitude. At the point $`V=V_c`$ where $`P(V)`$ reaches its maximum, the latter pair converges at the origin and moves onto the imaginary axis. (This does not render the soliton stable though, as the other pair remains on the real axis.) This scenario is exemplified by the curve $`h=0.27`$ in Fig.2.

Next, let $`h`$ be greater than $`0.28`$. For these $`h`$ the branch $`P(V)`$ emanating from the origin, turns back at some $`V=V_{\mathrm{max}}`$ (Fig.2), with the derivative $`P/V`$ remaining strictly positive for all $`VV_{\mathrm{max}}`$. Below we will describe the transformation this solution undergoes when continued beyond the “turning point”, while here we only wish to emphasize that no new zero eigenvalues can appear at this point. The reason is that $`V=V_{\mathrm{max}}`$ is a bifurcation point of solutions of the ordinary differential equation (3) but not of the partial differential equation (1)-(2). (In other words, $`V`$ is an “internal” parameter characterising the solution and not an “external” control parameter.) Indeed, the soliton is a member of a two-parameter ($`x_0`$ and $`V`$) family of solutions of Eq.(1)-(2) and hence for any $`V`$ there are two zero eigenvalues in the spectrum of the linearised operator $``$. Consequently, although being a turning point for the ODE (3), the value $`V=V_{\mathrm{max}}`$ is in no way special as far as the PDE (1)-(2) and its linearisation are concerned. No changes of the soliton’s stability properties occur at this velocity.

How does one type of behaviour of the curve $`P(V)`$, occurring for $`h>0.28`$, replace the other one, arising for $`h<0.28`$? We scanned the interval $`0.28000<h<0.28020`$ and discovered a tiny region of transitional behavior, around $`h=0.28005`$. For this $`h`$, $`P(V)`$ grows until it reaches a maximum at $`V_c=1.051`$ and then starts decreasing, as in the case of $`h<0.28`$. However, the curve does not decay all way to $`P=0`$ as would be the case for $`h<0.28`$, but reaches a minimum at $`V_{cc}=1.0563`$. After that, the momentum starts growing, and, at $`V_{\mathrm{max}}=1.0565`$, the curve $`P(V)`$ turns back — as for $`h>0.28`$! To give an idea of how small this window of transitional behaviour is, it suffices to say that for $`h=0.28000`$ the momentum $`P(V)`$ tends to 0 as $`Vc`$, whereas for $`h`$ as close as $`0.28010`$, the curve $`P(V)`$ has a “turning point”.

Fig.3 Stability diagram for the $`\psi _+`$ and $`\psi _{}`$ solitons on the $`(h,V)`$-plane. In the region marked “stable” one of the two one-soliton solutions is stable whereas the other one is not. Across the solid line, the corresponding soliton looses its stability to an oscillatory or monotonically growing mode. No solitons exist in the dashed region.

What happens to the $`\psi _+`$ soliton with $`h>0.28`$ (more precisely, with $`h0.28010`$) as we continue it beyond the turning point? Fig.4 (a) shows the momentum as a function of $`V`$. The point of intersection with the vertical axis $`V=0`$ corresponds to the twist solution (Eq.(III B) with $`z=\zeta `$ given by Eq.(33).) On the diagram Fig.4, it is marked as $`\psi _T`$. As we mentioned before, the $`(V=0)`$-twist is a representative of a two-parameter family of stationary solutions of Eq.(2). Consequently, there should be four zero eigenvalues in the spectrum of the operator $``$ in this case, with two linearly independent eigenfunctions given by $`_x\mathrm{\Psi }(z;x)|_{z=\zeta }`$ and $`_z\mathrm{\Psi }(z;x)|_{z=\zeta }`$. (Here $`\mathrm{\Psi }(z;x)`$ is a two-component vector formed by the real and imaginary parts of Eq.(III B).) Numerically, we observed that as we approach the ($`V=0`$)-twist from the direction of positive $`V`$, a pair of opposite eigenvalues converges at the origin on the complex plane. The curve $`P(V)`$ does not have an extremum at this point and this may seem to be in contradiction with predictions of section II C. The paradox is resolved as soon as one recalls that the extremality condition $`P/V=0`$ was derived under the assumption that there is only one eigenvector associated with the zero eigenvalue whereas we have two linearly independent null eigenvectors in the case at hand.

As we continue further into the region $`V<0`$, the twist gives rise to a variety of multisoliton complexes; we shall describe them in the next subsection. Here we will restrict ourselves to the region $`V>0`$ where this branch can still be regarded as a branch of one-soliton solutions. Although these solutions undergo similar transformations for all $`h`$ in the interval $`(0.28,1)`$, there are a few differences with regard to the trajectories of eigenvalues on the complex plane. One difference worth to be mentioned is that for the driving strength $`h=0.3`$ and larger $`h`$, the quadruplet of complex $`\lambda `$ persists on the entire upper branch of $`P(V)`$ (i.e. for all $`V>0`$). This is in contrast to the case of $`0.25<h<0.28`$, where the complex quadruplet converges on the real axis. Near the left end of the interval $`0.28<h<1`$ (e.g. for $`h=0.2802`$), we have an intermediate pattern. Similarly to the case $`h<0.28`$, here the complex quadruplet converges on the real axis somewhere on the lower branch of $`P(V)`$ (i.e. before the turning point), but as we move onto the upper branch, the two emerging real pairs reunite quickly and the complex quadruplet reappears.

Next, as we know, there are only two points where a pair of eigenvalues can pass from the real onto the imaginary axis, or vice versa. One point is $`V=V_c`$ where $`P/V=0`$, and the other one is $`V=0`$. Therefore the dynamics of eigenvalues depends on which of the two points comes first, or, equivalently, whether the upper branch of $`P(V)`$ has the maximum for positive or negative $`V`$. For smaller values of $`h`$ in the interval $`(0.28,1)`$ (e.g. $`h=0.3`$), where $`V_c>0`$, two imaginary eigenvalues move to the real axis at $`V=V_c`$. These imaginary eigenvalues have detached from the continuous spectrum somewhere before the turning point (i.e. on the lower branch of $`P(V)`$.) The two newly born real eigenvalues first diverge from the origin but then reverse and, at $`V=0`$, move back onto the imaginary axis. For larger $`h`$ (e.g. $`h=0.7`$), where $`V_c<0`$, the pattern is different. Here, the two imaginary eigenvalues become real not at the point $`V_c`$ but at $`V=0`$. Subsequently, as we continue the branch to negative velocities, another pair of imaginary eigenvalues detaches from the continuum and at the point $`V_c<0`$ two (imaginary or real) eigenvalues pass through the origin.

Fig.3 shows the stability diagram of the $`\psi _+`$ and $`\psi _{}`$ solitons on the $`(h,V)`$-plane. For the $`\psi _+`$ soliton, the range of stable velocities approaches $`0V<\sqrt{2}`$ as $`h0`$, while the stability range of $`\psi _{}`$ tends to $`\sqrt{2}<V2`$. Finally, the domain of stability in the $`h=0`$ case is the union of the above two ranges: $`0V2`$.

### D Other branches; $`h>0.28`$

As we continue it to negative velocities, the twist (we are using this name here for $`V0`$-deformations of the quiescent twist solution) gradually transforms into a complex of two twists (plotted in Fig.5(a)). A further continuation of this branch takes us, via several “turning points”, to a solution that can be interpreted as an association of the twist and two $`\psi _{}`$ solitons of opposite polarities (denoted $`\psi _{(T)}`$). This solution is depicted in Fig.5 (b).

Fig.4 (a): The full bifurcation diagram for $`h=0.7`$. All branches shown in this figure are unstable. (For $`h=0.7`$ the soliton $`\psi _{}`$ is stable at large velocities but the stable region lies beyond the frame of this figure — see Fig.2.) More solution branches can be obtained by the reflection $`VV`$, $`PP`$. (b): The corresponding $`E(P)`$ dependence. The stable branch of $`\psi _{}`$ is depicted by a thick line at the bottom of the figure.

Fig.4 (b) shows the energy of different branches as calculated by Eq.(8). We have eliminated the dependence on the soliton’s velocity between $`P(V)`$ and $`E(V)`$ to obtain $`E`$ as a function of $`P`$. The purpose of this “Legendre transformation” is the following. The stationary equation (3) can be regarded as a condition that the energy (7) be stationary under the fixed momentum (6): $`(\delta E)_P=0`$ or, equivalently, $`\delta (EVP)=0`$, where $`V`$ is the Lagrange multiplier. Consequently, a steadily travelling soliton satisfies the relation of the Hamilton mechanics, $`V=E/P`$, and its velocity is given by the slope of the corresponding $`E(P)`$ branch in Fig.4 (b). The relation $`\delta E=V\delta P`$ implies that the curves $`E(V)`$ and $`P(V)`$ have extrema at the same set of points $`V=V_c`$ which will appear as cusps on the $`E(P)`$ plot. Since the function $`P(V)`$ has extrema precisely at points where linearised eigenvalues move from the imaginary to real axis, branches of the $`E(P)`$ curve separated by the cusps may have different stability properties. One such change of stabilities does indeed occur in Fig.4 (b) where the $`\psi _{}`$-curve is seen to be partitioned into stable and unstable branch. Solitons on the stable branch have lower energies than unstable solitons with the same values of the momentum.

Fig.5 (a): The $`\psi _{(TT)}`$ complex. (b): The $`\psi _{(T)}`$ solution. (c): A complex of two twist solitons arising from the continuation of the $`\psi _{(++)}`$ bound state. (Note the difference from the other two-twist complex shown in (a).) Solid line: real part; dashed line: imaginary part

Another branch emanating from the origin in Fig.4 (a), is a bound state of two solitons $`\psi _+`$. This solution was not obtained by the continuation from $`V=0`$ as Fig.4 may seem to be suggesting. Instead, we fixed a nonzero $`V`$ and continued in $`h`$ from the value $`h=0.05`$ where the complex $`\psi _{(++)}`$ arises from the $`V`$-continuation of the twist soliton (see section IV E). Omitting details of this procedure, we start the description of the resulting branch at some point $`(V,P)`$ away from the origin. As we approach the origin from this point, the separation between the solitons $`\psi _+`$ in the complex $`\psi _{(++)}`$ rapidly increases so that the field values between the two solitons become exponentially small. For example, for $`h=0.7`$, the (numerically calculated) separation at the point $`V=0`$ was equal to $`z21`$. The value of $`|\psi |`$ at the point on the $`x`$-axis, equally distanced from the left and right soliton, was of order $`10^6`$. Consequently, the nonlinear term in the equation (3) becomes negligible away from the solitons’ core and, in spite of an extremely small value of the residual that we used in our numerical algorithm ($`10^{10}`$), we were unable to distinguish between a genuine bound state and a linear superposition of two distant solitons. We conjecture that the complex $`\psi _{(++)}`$ exists all way to $`V=0`$ but as $`V0`$, the intersoliton separation $`z\mathrm{}`$. Another indication to this effect is that as $`V0`$, the imaginary part of the solution tends to zero, rapidly and uniformly. For example, the imaginary part of the above-mentioned numerical solution with the real part between the solitons $`10^6`$, was smaller than $`10^{13}`$ for all $`x`$. Since the only pure real solution that exists for $`V=0`$ is the (single) soliton $`\psi _+`$, the $`V0`$ limit of the $`\psi _{(++)}`$ complex should be an infinitely separated pair of the $`\psi _+`$’s.

If we, conversely, continue our solution away from the origin, the curve $`P(V)`$ turns left at some $`V`$ and the complex $`\psi _{(++)}`$ transforms into what can be interpreted as a bound state of two twists (denoted $`\psi _{(TT)}`$ in Fig.4.) This solution is depicted in Fig.5(c). As $`V0`$, the momentum of this bound state tends to zero (Fig.4 (a)). Unfortunately, we were only able to obtain this solution away from some small neighbourhood of $`V=0`$. (For $`h=0.7`$, the smallest value of the velocity for which we were still able to find the solution in question, was $`V=0.000283`$.) Whether this branch can be continued to $`V=0`$, remains an open question.

### E Other branches; $`h<0.28`$

As we have mentioned, for $`h<0.28`$ the branch $`\psi _+`$ extends all way to $`V=c`$ where it merges with the zero solution. No other solutions can be obtained from the $`\psi _+`$ soliton. However, in this case we can obtain new branches by continuing the (quiescent) twist soliton, Eq.(III B) with $`z=\zeta `$.

It is convenient to start our description with the motionless twist solution with the negative momentum. As we move in the direction of positive $`V`$, the twist gradually transforms into a bound state of two $`\psi _+`$-solitons (see Fig.6 (a)). At some $`V=V_{max}`$ the branch turns back, shortly after which, at the point $`V=V_c`$, the momentum reaches its maximum and starts decreasing. Adjacent to the turning point is a small range of velocities $`V_cVV_{max}`$ where we have two stable solutions corresponding to each $`V`$. If we continue the branch with positive-momentum twist solution, also in the direction of positive $`V`$, the solution gradually transforms into a complex of two $`\psi _{}`$ solitons. The momentum reaches its maximum, starts decreasing, then the branch turns back in $`V`$ and we find ourselves approaching the origin on the $`(V,P)`$-plane (Fig.6 (a)). As we move towards the origin along the $`\psi _{(++)}`$ or along the $`\psi _{()}`$ branch, the separation between two solitons constituing the corresponding complex grows while the imaginary part of the solution tends to zero. Similarly to what we had for larger $`h`$ (section IV D), we conjecture that the separation becomes infinite at $`V=0`$ in both cases.

Fig.6 (a): The full bifurcation diagram for $`h=0.05`$. Thick respectively thin lines depict stable respectively unstable branches. Note a region of stability of the complex $`\psi _{(++)}`$. Additional branches can be generated by employing the reflection symmetry $`VV`$, $`PP`$. (b): The corresponding $`E(P)`$ diagram.

Fig.6 (b) shows the corresponding $`E(P)`$ dependence. As in Fig.4 (b), the cusps mark points of the zero crossing by the stability eigenvalues. Note that as in the case of large $`h`$ (Fig.4 (b)), there are solitons with the same value of the momentum but different energies. Similarly to Fig.4 (b), the stable branch of $`\psi _{}`$ has the lowest energy; bound states on the stable $`\psi _T\psi _{(++)}`$ branch also have lower energies than their counterparts with the same $`P`$ and smaller $`|V|`$. However, in the case of the $`\psi _+`$ solitons we have an interesting reverse of fortunes: out of the two branches with the same $`P`$, the stable branch is the one with higher energy!

## V Nonlinear Stage of Instability

In this section we present results of our numerical simulations of the full time-dependent nonlinear Schrödinger equation (1) (with $`\gamma =0`$.) The objective was to study the nonlinear stage of the development of instabilities reported in the previous section and to identify the attractors emerging as $`t\mathrm{}`$. We utilised a split-step pseudospectral method, with $`2^{11}=2048`$ modes on the intervals $`40X40`$ and $`80X80`$, and with $`2^{12}=4096`$ modes on the interval $`(60,60)`$. The method imposes periodic boundary conditions $`\psi (L/2,t)=\psi (L/2,t)`$, $`\psi _X(L/2,t)=\psi _X(L/2,t)`$.

We have simulated the evolution of moving solitons unstable against an oscillatory mode and those with a positive, nonoscillatory, eigenvalue in their linearized spectrum. One of our conclusions here is that both types of instability give rise to the same asymptotic attractors. (This is in agreement with earlier simulations of motionless solitons .)

### A The decaying breather

Depending on the value of the driving strength, the initial conditions and the choice of the parameters of the numerical scheme, we observed one of the two scenarios. In the first scenario the soliton transforms into a bell-shaped structure, with a small amplitude and large spatial width, oscillating approximately as $`\psi e^{i\omega t}`$, with negative $`\omega `$. This localised solution was previously encountered in numerical simulations of Ref. where it was termed breather. The amplitude of the breather slowly decays with time and the width slowly grows.

We have detected this scenario for the driving strength $`h=0.1`$, with the initial condition in the form of the $`\psi _+`$ soliton travelling with the velocity $`V=0.05`$ and with $`V=0.8`$. (For both values of the velocity the $`\psi _+`$ soliton is unstable against an oscillatory mode.) Unlike earlier simulations  which started with the initial condition in the form of a quiescent unstable soliton and gave rise to a quiescent breather, the breather emerging from a travelling soliton has a nonzero speed.

One may naturally wonder whether the speed of the breather will decay to zero or approach a nonzero constant value as $`t`$ increases. Our simulations seem to support the latter hypothesis. In one run, the speed of the breather evolving out of the soliton travelling with the initial velocity of $`V=0.05`$, was seen to slowly grow and gradually approach the constant value of $`0.1`$. This simulation was repeated, with the same parameters of the numerical scheme and an initial condition which was only different from the previous one due to interpolation errors of order $`10^6`$. In this run the breather was first seen to slow down, stop but then start moving in the opposite direction with the velocity close to -0.2; see Fig.7 (a). (This remarkable sensitivity to the initial data deserves a separate comment; see below.) The velocity of the breather evolving out of the $`V=0.8`$ soliton, was tending to approximately $`2.1`$. However, for large $`t`$ the unambigous interpretation of the numerical data is hindered by the growth of the amplitude of the radiation background. The radiation waves emitted by the oscillating breather re-enter the interval via the periodic boundary conditions and at a certain stage their amplitudes become comparable with the amplitude of the breather. Consequently, the constant-velocity motion of the breather may have been induced by the interaction with the backround radiations.

Fig.7 The two types of asymptotic attractors resulting from the decay of the unstable steadily travelling solitons: (a) the decaying and (b) the growing breather. (a) corresponds to $`h=0.1`$ and the initial condition in the form of the $`\psi _+`$ soliton with $`V=0.05`$. In (b), $`h=0.05`$ and the initial condition was chosen as the $`\psi _{}`$ soliton with $`V=0.05`$. In both plots the emerging breather changes, spontaneously, its direction of motion. (Note that this happens not as a result of the reflection from the boundary, as the periodic boundary conditions are imposed.)

### B The growing breather

The decaying breather was detected in simulations on the interval $`(40,40)`$ with $`N=2^{11}`$ modes. Changing the numerical parameters produced an entirely different scenario, however, — for the same value of the control parameter in Eq.(1) ($`h=0.1`$), and for the same initial conditions ($`V=0.05`$ and $`V=0.8`$).

Namely, we increased the number of the Fourier modes to $`N=2^{12}`$ and the length of the interval first to $`L=120`$ and then to $`160`$. As in the case of $`L=80`$ and $`N=2^{11}`$, in simulations with the new values of $`N`$ and $`L`$ the unstable travelling soliton $`\psi _+`$ was seen to transform into a bell-shaped structure, oscillating roughly as $`\psi e^{i\omega t}`$. However, this time the emerging breather has a positive  frequency $`\omega `$; its amplitude is large and continues to slowly grow, while the width is narrow and keeps on decreasing (Fig.7 (b)).

This attractor was also observed previously in . It was found there that the decaying and growing breather coexist. Whether the evolution of the same unstable soliton settles to one or the other asymptotic attractor, was found to depend on the choice of the phase of a small perturbation applied to the initial condition. In our present simulations, the perturbation is introduced simply by changing the parameters of the numerical scheme.

We also examined initial conditions in the form of translationally unstable solitons, including the $`\psi _+`$-soliton with $`V=1.4`$ for the driving strength $`h=0.1`$ and the $`\psi _{}`$-soliton with initial velocities $`V=0.05`$ and $`V=1.4`$, for the driving strength $`h=0.05`$. For each of the above three situations the simulations were repeated with $`2^{11}`$ modes on the interval $`40X40`$, and with $`2^{12}`$ modes on the intervals $`(60,60)`$ and $`(80,80)`$. In all nine runs the unstable soliton was seen to evolve into the growing breather. (It is quite likely that some other choices of the numerical parameters may give rise to the decaying breather instead.)

The velocity of the growing breather may vary during its evolution. It can even wander erratically, changing the direction of its motion several times, but eventually, for $`t10^4`$ or even earlier, the speed of the breather locks on to some constant value. Since the amplitudes of radiation waves are comparable with the amplitude of the breather at that stage, this effect can be induced by the breather-radiation interactions.

## VI Concluding remarks and open problems

The main result of this paper is the demonstration of the existence of wide classes of travelling soliton solutions of the (undamped) parametrically driven nonlinear Schrödinger equation. We established the necessary conditions under which motionless solitons can be continued to nonzero velocities, and, in cases where these conditions were met, were indeed able to carry out the numerical continuation. The stability of all resulting branches of solutions was examined; oscillatory and translational instabilities identified, and the single-soliton stability chart compiled on the $`(h,V)`$-plane. The onset values of the translational instabilities, obtained numerically, were shown to verify the relation $`P/V=0`$ predicted by our theoretical analysis. As opposed to the case of the soliton $`\psi _{}`$, which undergoes similar transformations for any $`h`$, the result of the continuation of the $`\psi _+`$ has turned out to be sensitive to the value of the driving strength. We have identified two different transformation scenarios, one occurring for small and the other one for larger $`h`$, and described an interesting cross-over from one to the other.

In our analysis we paid a special attention to the trajectories of linearised eigenvalues on the complex plane. Apart from the information on the stability of different branches of solutions, the behaviour of the eigenvalues can give insight into the supercritical dynamics of solitons, i.e. dynamics beyond the instability threshold . The motion of the eigenvalues in the undamped case that we are currently concerned with, allows one to even predict the asymptotic attractors arising when there is a small but nonzero damping in the system . Relegating the corresponding bifurcation analysis to future publications, here we have restricted ourselves to a series of numerical simulations of the time-dependent NLS equation (1) (with $`\gamma =0`$.)

It is worth listing, separately, stable solutions obtained in this study. First, the quiescent soliton $`\psi _{}`$, which is always unstable with respect to a nonoscillatory mode, stabilizes when travels faster than a certain critical velocity. The stability boundary satisfies $`P/V=0`$. Second, the quiscent soliton $`\psi _+`$ is stable for $`h<0.064`$ and loses its stability to an oscillating soliton for $`h>0.064`$. For nonzero $`V`$ the stability region is shown in Fig.3. The lower boundary of this region corresponds to the onset of the oscillatory instability while along the upper boundary the soliton becomes unstable with respect to a nonoscillatory mode. The corresponding critical velocity satisfies $`P/V=0`$. Third, the bound state $`\psi _{(++)}`$ also displays a region of stability for small $`h`$ — see Fig.6.

We conclude this section by pointing out to several open questions. In the first place, it would be interesting to continue the $`\psi _{(T)}`$-branch in Fig.4. Will this multisoliton complex keep on attaching more solitons on its flanks? Another “open-ended” branch ($`\psi _{(TT)}`$) approaches the origin vertically down in Fig.4 (a); it would be interesting to continue this branch as well. The striking difference between the bifurcation diagrams for $`h>0.28`$ (Fig.4) and $`h<0.28`$ (Fig.6) is also worthy of a deeper analysis. Is the bifurcation diagram in Fig.6 complete, or there are other multisoliton branches similar to those arising for large $`h`$ (Fig.4)? The next open question concerns the decaying and growing breather solutions arising as a result of the growth of the instability of steadily travelling solitons. If the radiations were prevented from re-entering the interval of simulation, would the velocity of the breather decay to zero or approach some nonzero value determined by the initial velocity of the unstable soliton? Finally, a challenging problem is to find out whether travelling solitons can exist in the presence of damping. The effect of dissipation will be, of course, to attenuate the soliton. It is not obvious whether the spatially uniform parametric pumping is capable of compensating this type of losses and sustaining the soliton’s steady motion. It is fitting to note here that the parametrically driven, damped nonlinear Schrödinger equation has wide classes of stationary solitons , some of which may be continuable to nonzero $`V`$.

###### Acknowledgements.

We thank Dmitry Pelinovsky for drawing our attention to Ref. and solution (III B) presented therein. Useful conversations with Yuri Gaididei are also gratefully acknowledged. Special thanks go to Nora Alexeeva for providing us with a code for the time-dependent NLS equation and for her help with numerics. Finally, we are grateful to Professor Igor Puzynin for his continual encouragement and strong administrative support of this project. I.B. was supported by the NRF of South Africa and URC of the University of Cape Town. E.Z. was supported by an RFBR grant No.0001-00617.


format: text

text:
# Effects of non-perturbatively improved dynamical fermions in UKQCD simulations

## 1 SIMULATION STRATEGY

The first generation of full QCD simulations has confirmed that, for many purposes, the quenched approximation is indeed very good and is hard to improve upon with finite computational resources. For recent reviews see . Following an exploratory series of simulations using non-perturbatively improved Wilson fermions , the UKQCD collaboration has now completed the first phase of a sequence of matched simulations in which some attempt has been made to adjust the bare simulation parameters so that physically relevant quantities remain approximately constant . In particular, we have generated ensembles in which the Sommer scale parameter $`r_0`$, expressed in lattice units, is matched as closely as is practical.

Table 1 contains a summary of the simulation parameters and ensemble statistics.

This strategy is motivated partly by pragmatism and partly by theoretical considerations. With access to limited computational resources, it was clear that the Collaboration would be unable to address immediately all the limits required to recover a fully physical system: chiral, continuum and infinite volume limits. We have therefore adopted a partially unquenched approach in which, at finite lattice spacing, the dynamical quark mass (two degenerate flavours) is treated as an independently tunable parameter. Rather than hold the bare coupling fixed, we have chosen to hold $`r_0/a`$ fixed with the expectation that the discretisation and finite size effects which are inevitably present should be better under control. The choice of $`r_0`$ to set the scale removes the need to perform the chiral extrapolation in valence quark mass which use of a non-gluonic observable would entail. Of course an overall uncertainty in the phenomenological value of $`r_0`$ remains ($`2\%`$), but the use of $`r_0`$ provides a standard which is well-defined both theoretically and computationally. Here we use $`r_0=0.49`$ fm.

In the second phase, we are attempting to vary the lattice spacing while holding the ratio of pseudoscalar to vector meson mass fixed. Using fully non-perturbatively improved fermions , we expect that the discretisation errors are $`𝒪(a^2)`$. The fact that the improvement coefficient function $`c_{\mathrm{SW}}(\beta )`$ is only well-determined for $`\beta 5.2`$ limits the accessible lattice spacing from above. The spacing is limited from below by the need to have an acceptably large physical volume.

We have begun a cooperative programme of simulations with the QCDSF collaboration (see talks by Pleiter and Stüben ). It is planned that simulations on larger volumes and at complementary values of the bare parameters will be obtained. Some initial simulations from QCDSF are also summarised in Table 1. Simulation points are indicated in Figure 1 which has, superimposed on it, curves of constant $`r_0/a`$ and $`\xi m_\pi /m_\rho `$.

These curves show a simple polynomial interpolation which gives a very approximate description of the data obtained so far and provides guidance on the choice of future simulations points.

The lightest quark mass simulation which we have so far been able to achieve with useful statistics was at $`(\beta ,\kappa )=(5.20,0.1355)`$ (see Table 1). We are currently generating configurations which are expected to match the lattice spacing of this ensemble and the value of $`\xi `$ for the ensemble at $`(5.20,0.1350)`$. An attempt to simulate at $`(5.20,0.13565)`$ where $`\xi 0.56`$ (Table 1) was halted since the performance of the HMC algorithm was found to be unacceptably slow. For all ensembles studied, the measured integrated autocorrelation time for the plaquette (typically less $`\stackrel{<}{}20`$ HMC trajectories) showed no increase with decreasing quark mass. The same was true for measurements of the topological charge and glueball correlators where the integrated autocorrelation time was of the order of $`2530`$.

## 2 STATIC POTENTIAL

The static potential was extracted using an optimised variational method . Within the range $`0.35\stackrel{<}{}|𝐫|\stackrel{<}{}\mathrm{\hspace{0.17em}1.25}`$ fm, all matched data including quenched are well described by a simple bosonic string model ($`e=\pi /12`$).

$$[V(r)V(r_0)]r_0=(1.65e)\left(\frac{r}{r_0}1\right)e\left(\frac{r_0}{r}1\right)$$

when expressed in units of $`r_0`$. Figure 2 shows the deviations from this model.

For $`r\stackrel{<}{}\mathrm{\hspace{0.17em}3}a`$ there is clear evidence for (a) strong discretisation effects and (b) a systematic depression of the potential with decreasing dynamical quark mass. Parametric fits show that the effect corresponds to an increase in the effective strong coupling of $`18\pm 10\%`$ with respect to the quenched value. The simplest perturbation theory estimate of the effect of two light flavours is also $`18\%`$.

The values of $`r_0/a`$ used to obtain the effective lattice spacing shown in Table 1 were extracted from parametric fits to the potential (including off-axis measurements) but over a restricted range of $`|𝐫|`$. In fact the values of $`r_0/a`$ are rather stable with respect to the range used.

## 3 MESONS and GLUEBALLS

Figure 3 shows the valence quark mass dependence ($`m_{\mathrm{PS}}^2`$) of the vector meson mass in comparison with physical mesons having strange quark content.

New matched quenched data underline the conclusion that the lattice data are not far from the physical data even at finite lattice spacing and that there is little further improvement to be seen at the lightest sea quark masses reached (corresponding to $`m_\pi /m_\rho 0.6`$).

Results for baryons and for quark masses will be presented elsewhere .

In Figure 4 we show preliminary data for the mass of the scalar glueball ($`A_1^{++}`$ lattice state) obtained using optimised variational techniques.

The basis made use of 7 operators and 6 levels of fuzzing. Another important ingredient was the use of non-zero momentum operators $`|𝐩|^22a^2`$) to obtain effective energies from which effective masses could also be extracted. For these momenta, good plateaux were again obtained.

The figure shows the zero-momentum effective masses. Data for the lightest quark ensemble ($`\kappa ^{\mathrm{sea}}=0.1355`$) are compared with a ‘world average’ of quenched scalar masses interpolated to the relevant lattice scale. This provides reasonably convincing evidence of a decrease in the mass of order $`10\%`$ for the scalar glueball at fixed lattice spacing.

One possible source of (finite size) error could be mixing with a scalar ‘torelon’ which has an additional decay mode in the presence of fundamental charges. Direct measurements of the energy and vacuum expectation values of such states  show that the mixing is expected to be small on these lattices.

Measurements of the tensor glueball have also been made but larger statistical errors make the corresponding comparison with quenched results less fruitful.

## 4 TOPOLOGICAL SUSCEPTIBILITY

The matched ensembles have also been used to look for evidence of suppression of instanton effects by light quark modes . In Figure 5 we show the suitably scaled topological susceptibility obtained by standard cooling techniques applied to the gauge configurations after every tenth trajectory.

Interestingly, the measured decorrelation rate of the topological charge was little slower than for the plaquette. The data show a clear suppression with dynamical quark mass. Since the leading chiral behaviour is expected to be

$$\frac{r_0^2\chi }{m_\pi ^2}=(r_0f_\pi )^2/8+c_1(r_0m_\pi )^2,$$

one may use models to interpolate this data and hence extract an estimate of the pion decay constant $`f_\pi `$. Note that, in contrast to methods based on matrix elements, no non-perturbative renormalisation factors are required to obtain the value of $`f_\pi `$. Using $`r_0=.49`$ fm to set the scale and ignoring finite-$`a`$ corrections, one estimates $`f_\pi =148\pm 7_{14}^{+25}`$ MeV in comparison with $`132`$ MeV from experiment.

The use of matched data in this analysis is critical to its success in that, otherwise, the underlying suppression can easily be masked by the strong ($`r_0^4`$) change in scale expected when changing $`\kappa ^{\mathrm{sea}}`$ at fixed $`\beta `$.

## 5 SUMMARY and CONCLUSIONS

The use of matched ensembles of dynamical fermion configurations has allowed a preliminary study of possible unquenching effects prior to a full analysis of the continuum and chiral limits. Evidence is presented for a lowering of the scalar glueball mass, charge screening at short distances and of the suppression of instanton effects at light quark masses. For the quark masses studied ($`m_\pi /m_\rho \stackrel{>}{}\mathrm{\hspace{0.17em}0.6}`$) there is little evidence of dynamical quark effects in the light hadron spectrum. Further simulations exploring the lattice spacing dependence and quark mass dependence are underway.


text:
# Spatially homogeneous ground state of the two-dimensional Hubbard model

## I Introduction

Since the discovery of high-temperature superconductivity, the two-dimensional Hubbard model  has been the subject of a huge amount of work. Indeed, it is widely accepted that, despite the complicate structure of these materials, a key role is played by the electron correlation in the $`\mathrm{CuO}`$ planes, well represented by the Hubbard model.

This model and the related $`tJ`$ model reproduce, at least qualitatively, some of the physical properties of cuprates. For example, at half-filling (i.e., when the number of electrons is equal to the number of sites) the ground state (GS) is an antiferromagnetic insulator, and upon doping the antiferromagnetism is strongly suppressed. The possible insurgence of superconductivity upon doping and the symmetry of the order parameter in the model are still open questions. Another important point to study is the stability with respect to phase separation (PS) or charge-density waves (CDW) near half-filling. Indeed, both PS and CDW have been experimentally observed in different systems  and different authors  have pointed out a possible relation between charge instability and superconductivity.

Because of the strongly interacting nature of these systems, important insights came from nonperturbative numerical methods.

The exact diagonalization methods, e.g., Lanczos, are strongly limited by the exponential growth of the Hilbert space, and in practice it is possible to diagonalize only up to about 20 sites. A remarkable development of the exact diagonalization strategy is the density matrix renormalization group , which allows us to compute the GS using an iteratively improved basis. Despite the accuracy of this technique for one-dimensional and quasi-one-dimensional systems, there is no straightforward generalization for higher dimensions.

Attractive alternatives are stochastic methods such as Quantum Monte Carlo (QMC) that do not suffer from severe limitations in terms of lattice size and hence allow us to study fairly large systems. The simplest QMC simulation is the variational one, in which it is possible to compute the expectation values of the energy and all the correlation functions for a given wave function. The limitation of this method is that one has to guess a priori the form of the GS.

QMC methods allow us to sample directly the GS by using projection techniques, the simplest one being the so-called power method, namely the simple idea that an iterative application of the Hamiltonian $`H`$ filters out the ground state starting from an arbitrary state nonorthogonal to it. These approaches suffer, in fermionic systems, from the so-called sign problem. Indeed, due to the antisymmetry of the wave function under permutations of two particles, one gets after a few power iterations opposite contributions leading to the cancellation of large and fluctuating weights.

An alternative approach, which is quite efficient for the Hubbard model, is the auxiliary field quantum Monte Carlo (AFQMC) . The Coulomb interaction is linearized by using a Hubbard-Stratonovic transformation. In this way it is possible to apply exactly the operator $`e^{\beta H}`$ to a Slater determinant. In practice, AFQMC is very accurate only for small values of the Coulomb repulsion, whereas, when the interaction strength becomes comparable with the bare bandwidth, large fluctuations make the simulations rather unstable and inefficient, because the sign problem becomes particularly severe for strong Coulomb repulsion.

A sophisticated implementation of the power method is the Green-function Monte Carlo (GFMC) applied to lattice models . In principle, this method gives unbiased results but in practice, for fermion systems, the sign problem makes any simulation prohibitive, because of statistical errors. Recently, many attempts have been made to deal with fermions and overcome the sign problem . The fixed-node approximation  is introduced for the GFMC and replaces the original Hamiltonian with an effective one free of the sign problem, in which the nodes are fixed to be those of the so-called guiding wave function. In practice, the most common choice is the best variational wave function available, i.e., the one with the lowest energy. The fixed-node energy is an upper bound to the true GS energy and the method is variational. Moreover, it becomes exact if the guiding function is chosen to be the exact GS function. Although this approximation can be uncontrolled and needs an external input, it is quite accurate with respect to both energy and correlation functions for small sizes.

An alternative approach, recently proposed by one of us, is to use the GFMC with a stochastic reconfiguration (SR) . The main idea of this method is to use a suitable reference dynamics free of the sign problem to constrain the true dynamics. The fixed-node and a variational dynamics are natural candidates as references. At each reconfiguration the amplitudes of the GS are chosen as small perturbations of the reference wave function amplitudes such that the mixed averages of a few physically relevant quantities are conserved. The method is exact if all the operators in the Hilbert space are reconfigured, but it has been shown that a major role is played by a very small number of operators .

In this paper, we present a systematic study on the Hubbard model comparing different QMC techniques by studying in particular the stability of PS and CDW near half-filling for various values of the Coulomb repulsion. In Sec. II we introduce the Hubbard model and we discuss the implementation of our approximations, in Sec. III we report our results, and finally in Sec. IV we give a brief summary.

## II The Model

We consider the Hubbard model on a square lattice of $`L`$ sites with $`N=N_{}+N_{}`$ particles and $`N_{}=N_{}`$, where $`N_{}`$ ($`N_{}`$) is the number of spin-up (-down) particles. In order to study PS as close as possible to half-filling using only closed-shell configurations, we consider square lattices tilted by $`45^{}`$ with $`L^2=2l^2`$ and $`l`$ odd. In this way half-filling is a closed shell and the first doped closed shell has eight holes independent of $`L`$. The Hamiltonian reads

$$H=t\underset{i,j,\sigma }{}c_{i,\sigma }^{}c_{j,\sigma }+U\underset{i}{}n_{i,}n_{i,},$$
(1)
where $``$ stands for nearest neighbors, $`c_{i,\sigma }`$ ($`c_{i,\sigma }^{}`$) destroys (creates) an electron with spin $`\sigma `$ at site $`i`$, and $`n_{i,\sigma }=c_{i,\sigma }^{}c_{i,\sigma }`$. In the following, all energies are measured in units of $`t`$.

As already pointed out above, in the presence of the sign problem the choice of the guiding wave function is crucial. Our wave function reads

$$|\mathrm{\Psi }_G=𝒫_{Sz=0}𝒫_g𝒥|𝒟,$$
(2)
where $`|𝒟`$ is a Slater determinant in which the orbitals are suitably chosen (see below), $`𝒫_{Sz=0}`$ is the projector onto the subspace with $`N_{}=N_{}`$, i.e., with zero total spin component along the $`z`$ axis, $`𝒫_g`$ is a Gutzwiller operator that inhibits the double occupancies $`𝒫_g=exp(g_in_in_i)`$, where $`g`$ is a variational parameter, and $`𝒥`$ is a Jastrow factor $`𝒥=exp(\frac{\gamma }{2}_{i,j}v_{i,j}S_i^zS_j^z)`$, where $`\gamma `$ is another variational parameter and $`v_{i,j}`$ is taken from spin-waves theory . Care must be taken in the choice of the orbitals appearing in the Slater determinant. The most common choice is to take the orbitals from a Hartree-Fock (HF) approximation of the Hamiltonian breaking the SU(2) spin rotation symmetry along the $`z`$ axis.

$$H=t\underset{i,j,\sigma }{}c_{i,\sigma }^{}c_{j,\sigma }+\frac{U}{2}\underset{i,\sigma }{}\left[n_i\sigma (1)^{R_i}m_i\right]n_{i,\sigma },$$
(3)
where

$`n_i`$ $`=`$ $`n_{i,}+n_{i,},`$ (4)
$`m_i`$ $`=`$ $`(1)^{R_i}\left[n_{i,}n_{i,}\right].`$ (5)
We consider only fillings which are closed shells for $`U=0`$ and where a solution with constant density $`n_i=\frac{N}{L}=n`$ and staggered magnetization $`m_i=m`$ are found. In this case, the HF many-body wave function can be written as

$$|𝒟=\underset{k}{\overset{1,\mathrm{},N_{}}{}}\beta _{k,}^{}\underset{q}{\overset{1,\mathrm{},N_{}}{}}\beta _{q,}^{}|0,$$
(6)
where the quasiparticles have definite momentum modulo $`Q=(\pi ,\pi )`$ and definite spin, since the antiferromagnetic order parameter is along the $`z`$ axis

$$\beta _{k,\sigma }^{}=u_kc_{k,\sigma }^{}+\sigma v_kc_{k+Q,\sigma }^{},$$
(7)
$`k`$ is in the reduced magnetic Brillouin zone, and $`u_k`$ and $`v_k`$ are defined in Ref. . It is worth noting that for $`U/t\mathrm{}`$, $`u_k=v_k=\frac{1}{\sqrt{2}}`$, namely the spin up and the spin down are in different sub-lattices (classical Néel state).

In a previous work  the wave function (2) with $`|𝒟`$ given by Eq. (6) has been found to be a rather poor approximation for large $`U/t`$ at half-filling. In particular, in this representation the Jastrow factor $`𝒥`$ does not play any important role.

We propose a wave function which is a straightforward generalization of the one successfully used for the Heisenberg model . The fundamental ingredient is to allow spin fluctuations perpendicular to the staggered magnetization. An easy implementation of this idea is to put the magnetization in the $`xy`$ plane allowing transverse fluctuations along the $`z`$ axis through a Jastrow-like factor . This is achieved by a $`\frac{\pi }{2}`$ rotation $`U_y(\frac{\pi }{2})`$ around the $`y`$ axis of the canonical operators:

$`U_y^{}({\displaystyle \frac{\pi }{2}})c_{i,}^{}U_y({\displaystyle \frac{\pi }{2}})`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(c_{i,}^{}+c_{i,}^{}\right),`$ (8)
$`U_y^{}({\displaystyle \frac{\pi }{2}})c_{i,}^{}U_y({\displaystyle \frac{\pi }{2}})`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(c_{i,}^{}c_{i,}^{}\right).`$ (9)
The fermionic part of our guiding wave function is therefore defined as a Slater determinant of the transformed orbitals,

$`\beta _{k,+}^{}`$ $`=`$ $`U_y^{}({\displaystyle \frac{\pi }{2}})\beta _{k,}^{}U_y({\displaystyle \frac{\pi }{2}}),`$ (10)
$`\beta _{k,}^{}`$ $`=`$ $`U_y^{}({\displaystyle \frac{\pi }{2}})\beta _{k,}^{}U_y({\displaystyle \frac{\pi }{2}}),`$ (11)
namely it is given by

$$|𝒟=\underset{k}{\overset{1,\mathrm{},N_{}}{}}\beta _{k,+}^{}\underset{q}{\overset{1,\mathrm{},N_{}}{}}\beta _{q,}^{}|0.$$
(12)
Remarkably for $`U/t\mathrm{}`$ and half-filling, by construction Eq. (12) becomes the Néel state with spin quantization parallel to the $`x`$ axis, i.e. it has the correct Marshall sign on each of the $`2^L`$ configurations sampled by GFMC. In this limit, it is also clear why the Jastrow factor may be much more effective: being defined along the $`z`$ axis, it allows us to sample the quantum fluctuation perpendicular to the staggered magnetization. In the previous case instead both the Jastrow quantization axis and the order parameter were parallel, and for $`U/t\mathrm{}`$ there is no way to sample any fluctuation, the only possible configuration being the classical one.

It is worth noting that, although for the $`tJ`$ model the d-wave BCS wave function with a spin-rotationally invariant density-density Jastrow factor represents a very accurate variational state, for the Hubbard model at small and intermediate coupling ($`U10t`$) the best choice for the variational and guiding wave function is given by a Jastrow-Slater determinant with rotated orbitals Eqs. (10) and (11). Indeed, although the d-wave BCS wave function is a singlet and does not break the SU(2) symmetry, it has a very poor variational energy for the Hubbard model. The quality of the variational energy obtained with our Jastrow-Slater determinant remains considerably better than the BCS one at half-filling and $`U10t`$ even when the accuracy of the approximation is improved by the GFMC. Instead, in the doped case, the BCS wave function with GFMC is only slightly worse than the corresponding Jastrow-Slater determinant proposed in this work. This may suggest that antiferromagnetism is already suppressed at small finite doping, and d-wave superconductivity is a possible stable phase especially at large $`U/t`$.

An important systematic improvement of the wave function can be achieved by performing exactly one Lanczos step starting from $`|\mathrm{\Psi }_G`$,

$$|\mathrm{\Psi }_L=\left(1+\alpha H\right)|\mathrm{\Psi }_G,$$
(13)
with $`\alpha `$ free parameter chosen to minimize the energy. This technique has been successfully used for the $`tJ`$ model both to improve the variational calculation  and as a starting point for power methods . Henceforth, we will denote by VMC and LS the results obtained with the wave function Eqs. (2) and (13), respectively, using variational Monte Carlo. Analogously, the symbols FN and FNLS will indicate the fixed-node approximation applied to the wave function Eqs. (2) and (13), respectively. Finally, the symbols SR will denote the stochastic reconfiguration approximation applied to the wave function Eq. (13).

In Tables I and II we report the energies of 18 and 10 electrons on 18 sites, respectively. At half-filling we compare the results using Eq. (2) with the Slater determinant $`|𝒟`$ given by Eqs. (6) and (12) for different approximations and values of $`U/t`$. Using Eq. (12) we obtain a sizeable improvement for large $`U`$’s ($`U10t`$). Notice that for $`U=20t`$ the best variational result (the FNLS) with Eq. (6) is worse than the simple VMC with Eq. (12). For 10 electrons the two Slater determinants give the same results. Indeed, for this doping the antiferromagnetic order is strongly suppressed and the Jastrow factor $`𝒥`$ does not play any important role.

## III Results

One of the most debated issues in strongly correlated electron models is the nature of the charge distribution in their GS. Recently, many authors  have addressed the question of PS in the $`tJ`$ model. It is well accepted that for $`Jt`$, holes tend to cluster together leaving the rest of the system in an antiferromagnetic state without holes. Although most of the calculations lead to the conclusion that a critical value of $`J`$ below which the GS is homogeneous exists, it is not clear what this value is at low doping, ranging from $`0.5t`$ and $`1.2t`$. Moreover, some authors  have suggested that just before PS, the GS has charge modulations. Recently, an accurate numerical study of a few chains on the Hubbard model  has shown evidence of stripes, i.e., CDW oscillations. However, this result appears limited to quasi-one-dimensional geometry, as also suggested by the authors.

It is well known that PS is characterized by an infinite compressibility in the thermodynamic limit. The compressibility can be related to the curvature of the energy with respect to the electron density,

$$\chi =\left(\frac{^2E}{n^2}\right)^1.$$
(14)
From the above definition we have that a divergent $`\chi `$ corresponds to a vanishing curvature of the energy as a function of density. Therefore, PS can in principle be detected by means of energy measurements for various densities. This is an appealing property, since the accuracy on the energy is usually better than that of any other observable for most numerical methods. Many previous numerical studies of PS have therefore concentrated on the energy curve, or equivalently, on the energy per hole $`e_h(\delta )=[e(\delta )e_H]/\delta `$, where $`e(\delta )`$ is the energy per site at a hole density $`\delta =1n`$ and $`e_H=e(0)`$ is the energy at half-filling . If $`\chi `$ diverges, $`e_h(\delta )`$ is flat in the thermodynamic limit and develops a minimum for $`\delta =\delta _c`$ in finite systems, due to the finite positive surface energy at the phase boundary.

This approach has been pursued by Cosentini et al. , using FN calculations. They found that there is a large region of PS in the phase diagram, at least for $`U10t`$, clearly in contrast with what is found in the $`tJ`$ model. In this paper, we consider the Hubbard model and we show that a study of PS instability is very difficult using only energy calculation. Instead, a careful calculation of charge-correlation functions strongly indicates that the GS is homogeneous.

In order to show that the energy calculations may overestimate the tendency to a PS instability, it is important to compare the GFMC results with some exact reference results. Previous studies have shown that it is important to consider relatively large lattice sizes since finite-size effects favor PS . We need, therefore, a reference result for large lattices, where exact diagonalization is not available. In the case of the Hubbard model for small $`U`$, the AFQMC is almost exact and represents the reference we need. As stated in the preceding section, we consider only closed shell doping because, at least for small $`U`$, huge finite size effects affect the physical properties in a rather drastic way. For instance, the large bare density of states near half-filling determines an unphysical and spurious PS up to the first closed shell .

In Fig. 1, we show the accuracy of the GFMC results obtained with different approximations compared with the AFQMC ones for a 162-site lattice and $`U=4t`$. For this coupling value, AFQMC does not provide evidence for PS. We plot $`[E(\delta )E_{hs}(\delta )]/E_{hs}(\delta )`$, where $`E(\delta )`$ and $`E_{hs}(\delta )`$ are the energies of GFMC and AFQMC, respectively, for a doping $`\delta `$. Besides the improvement in the absolute accuracy, the curves get flatter and flatter improving the approximation, but only the SR accuracy is almost doping-independent. In other words, we need a very accurate calculation to eliminate the spurious dependence of the variational energy upon doping . Even for the best variational method, the FNLS, although the accuracy on the energy is for all doping less than $`1\%`$, the difference in accuracy between, for example, the half-filled case and the first closed shell is still sizable. This difference is very important, because it represents just the energy scale determining or ruling out PS.

In Fig. 2, the function $`e_h(\delta )`$ is shown for the FN, FNLS, SR, and AFQMC methods. We need to use SR to exclude the occurrence of PS, where even the FNLS data would imply PS. The reason for this disappointing situation is that all the known variational methods are still too dependent on the guiding wave function. With the previous analysis, the resolution in energy necessary to detect or rule out PS is very hard to reach with statistical methods, especially for large $`U/t`$.

On the other hand, GFMC methods have proven to be reliable not only for energy calculations but also for correlation functions such as $`N(q)=n_qn_q`$, where $`n_q`$ is the Fourier transform of the electron density . For a phase-separated system, there are strong fluctuations in the density for small momenta and $`N(q0)`$ is expected to be strongly enhanced for small momenta, that is, for $`|q|\frac{2\pi }{\xi }`$, where $`\xi `$ is the characteristic length of the phase-separated region. Moreover, if $`\chi `$ diverges, also $`N(q0)`$ diverges, yielding an alternative tool to probe PS.

This method turns out to be more reliable, since it is based on a single calculation for a given doping value, whereas the evaluation using $`e_h(\delta )`$ involves a comparison between energies obtained by different simulations for different fillings, with the corresponding guiding wave function having different accuracies. Indeed with GFMC, $`N(q)`$ has been proved to be a very sensitive tool to look at for detecting PS. In the $`tJ`$ model, $`N(q)`$ has a very different shape for stable and unstable systems. A clear peak at the smallest $`q`$ indicates PS even when, as shown in the inset of Fig. 3, $`J/t`$ is very close to the PS boundary . Moreover, from $`N(q)`$ it is also possible to extract information about charge fluctuations at finite $`q`$’s, related to CDW. In Ref.  we have shown that in the $`tJ`$ model for $`J=0.4t`$, $`N(q)`$ has some peaks at finite $`q`$’s, surprisingly near to what was found in recent experiments . The knowledge of $`N(q)`$ allows us to extract more general results with respect to the simple study of $`e_h(\delta )`$. Furthermore, $`N(q)`$ is found to be much less size dependent than $`e_h(\delta )`$. For the $`tJ`$ model there is no appreciable difference between 98-site and 162-site lattices and for $`J=0.4t`$. Both calculations suggest that there is no PS, whereas an analysis using the FN approximation of $`e_h(\delta )`$ should lead to PS for 98 sites and to a homogeneous state for 162 sites. Indeed, in this case PS is only a size effect and a homogeneous state is found by increasing the accuracy of the method or by increasing the lattice size.

We computed $`N(q)`$ by means of the forward-walking technique , within the FNLS approximation, at half-filling and for the first few closed-shell configurations on a 162- and a 98-site lattice. The evaluation of the density-density correlation function is in principle possible even within SR by numerical differentiation of the energy with respect to an external field coupled to $`N(q)`$. However, this approach is very demanding and does not give a significant improvement on the FNLS results, which are very accurate. Indeed for the smallest $`q`$ vector for 90 electrons on the 98-site lattice, i.e., $`q=(2\pi /7,2\pi /7)`$, we found the the AFQMC gives $`N(q)=0.0932(2)`$ and the FNLS gives $`N(q)=0.096(1)`$.

In Fig. 3, $`N(q)`$ is shown for $`U=4t`$ at half-filling for a 162-site and a 98-site lattice and for 154 electrons on a 162-site lattice and 90 electrons on a 98-site lattice. No sign of divergence, and consequently of PS or CDW, is seen in the data. We also notice that the two sets of points for the half-filled systems lie on the same curve, showing that we have substantially reached the thermodynamic limit. For this value of $`U`$, $`N(q)`$ is essentially featureless for all doping we considered, suggesting that there are no charge instabilities at any finite length. The smallest doping we considered is $`\delta 0.049`$ and we cannot exclude that for smaller doping PS or CDW are present.

In order to investigate smaller doping, we should consider larger lattices. Unfortunately, the accuracy of the approximations considered decreases when increasing the size of the system and the 162-site lattice represents the largest lattice where the accuracy is acceptable. In Fig. 4, we report $`[E(0)E_{hs}(0)]/E_{hs}(0)`$ for various sizes and for different approximations: from the 18 sites to the 162 sites, the accuracy of FNLS changes from less than $`0.1\%`$ to about $`0.5\%`$. These indications prevent us from considering sizes larger than the ones presented in this paper.

Now we turn to larger Coulomb interactions and consider $`U=10t`$, where the AFQMC results are not reliable due to large fluctuations. In principle, GFMC techniques do not suffer from intrinsic limitations in the large-coupling regime and it is possible to consider any value of $`U`$. In practice we need an accurate knowledge of the nodes, i.e., an accurate guiding wave function. Our choice, Eq. (2), with orbitals given by Eqs. (10) and (11) is a very good approximation for the half-filling case. In Table I, we report the energies for various methods for 18 electrons on 18 sites at different $`U`$’s. Although all the approximations are quite size-dependent, the wave function becomes more and more accurate by increasing the Coulomb potential. Therefore, we expect that it also gives a good starting point at least close to half-filling.

We present results for $`U=10t`$, for which previous FN calculations based on $`e_h(\delta )`$ and a less accurate wave function  have shown PS up to $`\delta 0.15`$. Indeed, if we use $`e_h(\delta )`$ as a probe for PS, we find that the phase diagram shows a large instability region, confirming the results of Ref. . As for the $`U=4t`$ case, this instability is very likely to be a spurious effect, a consequence of the different energy accuracy for different doping. This possibility, which cannot be proved without knowing the exact energies at strong coupling (at present impossible), is instead very clearly supported by the calculation of the charge-correlation functions.

Figure 5 displays $`N(q)`$ for the same fillings of Fig. 3 and for 138 electrons on the 162-site lattice, which corresponds to $`\delta 0.148`$. All the correlation functions are definitely nondivergent for $`q0`$ and are qualitatively similar to the $`U=4t`$ case, indicating that the system is far away from a PS instability. Furthermore $`N(q)`$ does not show peaks at any finite momenta for this Coulomb interaction. This finding shows that the Hubbard and the $`tJ`$ model may have different behaviors as far as charge correlations are concerned. Indeed, by diagonalizing exactly the 18-site lattice, we find that charge fluctuations for the Hubbard and for the $`tJ`$ model are quite different in the small doping region for $`U=10t`$ and $`J=0.4t`$, respectively .

## IV Conclusions

An extensive GFMC analysis of the Hubbard model at low hole doping has been carried out. In particular, we have focused on the possible instability of the model with respect to PS and CDW. Comparing GFMC results with AFQMC ones in the weak-coupling region ($`U=4t`$), we show that detecting PS by means of energy results requires a very accurate calculation at all electronic densities. Indeed, the accuracy of the energy is strongly dependent on the electron density, and the signature of PS based only on energy calculations is clearly affected by this bias, leading to a spurious region of PS instability. In the case of the Hubbard model, this is particularly relevant because, while we are able to give a very good description of the half-filled case, in which the GS is an antiferromagnetic insulator, we are not aware of equally accurate descriptions of the doped state. Even for $`U=4t`$, it is necessary to use the really accurate SR and AFQMC technique to eliminate the doping dependence of the accuracy and to rule out PS.

On the other hand, PS (and CDW) instability can be probed more easily using charge correlation functions. This approach has various advantages. First, it is found that $`N(q)`$ has very small size effects and the thermodynamic limit is reached with about 100 sites, both for the Hubbard and the $`tJ`$ model. Second, the information contained in $`N(q)`$ does not depend on different densities, implying that a different accuracy as a function of doping does not introduce any external bias.

Instead of using energy calculations, which are very expensive at moderate and large $`U`$’s, we calculate the charge correlation functions and we are able to find clear evidence for the absence of PS up to $`U=10t`$ in the low doping regime.

###### Acknowledgements.

It is a pleasure to acknowledge useful discussions with L. Capriotti, M. Calandra, G. Bachelet, E. Koch, and A. Parola, to whom we are also grateful for critical and careful reading of the manuscript. This work has been partially supported by MURST-COFIN99 and by Istituto Nazionale per la Fisica della Materia.