text:
# An Exactly Solvable Model For The Interaction Of Linear Waves With Korteweg-de Vries Solitons

## 1 Introduction

This paper is concerned with solving the coupled system of equations:

$`_tA+_x\left[{\displaystyle \frac{1}{2}}A^2+_x^2A\right]`$ $`=`$ $`0`$ (1)
$`_tB+_x\left[\kappa AB+_x^2B\right]`$ $`=`$ $`0,`$ (2)
where $`\kappa `$ is a real parameter. Of course, the nonlinear equation for $`A(x,t)`$ is simply the Korteweg-de Vries equation, and it can be solved independently by the inverse-scattering transform \[GGKM67\]. The coupled system (1) and (2) is a partially linearized version of the system proposed by Hirota and Satsuma \[HS81\] as a model for the dynamics of coupled long waves.

The coupled system (1) and (2) can be solved exactly when $`\kappa =1`$ and when $`\kappa =1/2`$. The case of $`\kappa =1`$ is well-known, for then the equation (2) is just the KdV equation itself linearized about the solution $`A(x,t)`$. An elementary exact solution of the linear equation (2) in this case is given by $`B(x,t)=_xA(x,t)`$. Further solutions can be expressed in terms of derivatives of the squared eigenfunctions of the related Schrödinger operator with potential $`A`$ \[S83\].

The case of $`\kappa =1/2`$ is essentially different. In this case, the linear equation (2) is no longer the linearization of KdV about any solution. An elementary exact solution of the linear equation in this case is given simply by $`B(x,t)=A(x,t)`$. The main goal of this paper is to construct the general solution of the initial value problem for this linear equation when $`A(x,t)`$ is a multisoliton solution of KdV.

One way to make clear the difference between the cases $`\kappa =1`$ and $`\kappa =1/2`$ is to consider $`A(x,t)`$ to be the simple soliton solution of KdV (1):

$$A(x,t)=12\eta ^2\mathrm{sech}^2(\eta (x4\eta ^2t\alpha ))=V(\chi ),$$
(3)
where $`\chi =xct\alpha `$ and the velocity is $`c=4\eta ^2`$. If we look for solutions of the linear equation (2) that are traveling waves with speed $`c`$, we find the equation

$$\left[\kappa V(\chi )B(\chi )+B^{\prime \prime }(\chi )\right]^{}=cB^{}(\chi ).$$
(4)
Integrating once, using vanishing boundary conditions at $`\chi =\pm \mathrm{}`$, yields a Schrödinger eigenvalue problem for $`B`$:

$$B^{\prime \prime }(\chi )+\kappa V(\chi )B(\chi )=EB(\chi ),$$
(5)
where $`E=c`$. For $`\kappa =1/2`$, it follows from the fact that $`B(x,t)=A(x,t)`$ is a solution of (2) that the function $`B(\chi )=V(\chi )`$ is an eigenfunction of the Schrödinger operator with eigenvalue $`E=c=4\eta ^2`$. Since it has no zeros, it is the ground state eigenfunction. We will see below that there is also one excited state for $`\kappa =1/2`$, although it is not relevant here since it corresponds to a different velocity. On the other hand, for $`\kappa =1`$, $`B(x,t)=_xA(x,t)`$ is a solution of (2), which implies that the function $`B(\chi )=_xV(\chi )`$ is an eigenfunction of the Schrödinger operator with the same eigenvalue $`E=c=4\eta ^2`$. In this case, the eigenfunction has a single zero and therefore is the first excited state. It follows that there are at least two eigenvalues for $`\kappa =1`$. In fact, there are exactly three states in this case. A final observation is that from the construction of the one-soliton solution of KdV (see below) it follows that for $`\kappa =1/6`$, the function $`B(\chi )=V(\chi )`$ is an eigenfunction of the Schrödinger operator with eigenvalue $`E=c/4=\eta ^2`$. It is the ground state, and the only eigenfunction. These relationships are summarized in Figure 1.

We will have more to say about this picture when we discuss the trapping of linear waves by solitons for general values of $`\kappa `$ in § 6.

The rest of this paper is primarily concerned with developing the general solution of the initial value problem for (2) with $`\kappa =1/2`$ when $`A(x,t)`$ is an $`N`$-soliton solution of KdV (1). In § 2 we show how for $`\kappa =1/2`$ a large family of exact solutions of (2) can be obtained from the simultaneous solutions of the Lax pair for KdV. When the solution of KdV contains only solitons and no radiation, the construction of Lax eigenfunctions is completely algorithmic and algebraic, and consequently the corresponding family of solutions of (2) for $`\kappa =1/2`$ can be obtained with great practicality. In § 3 we then establish that in the $`N`$-soliton case there are enough of these exact solutions of (2) for $`\kappa =1/2`$ to expand for fixed $`t`$ any absolutely continuous $`L^1()`$ function of $`x`$. This fact then leads to a general solution formula for (2) simply by expanding the initial data. We present and discuss this formula in § 4. There will turn out to be $`N`$ linearly independent solutions that are asymptotically confined to the union of soliton trajectories and can therefore be considered to be bound states. In § 5 we compute the scattering matrix that relates the asymptotic behavior of bound states for $`t\mathrm{}`$ to the corresponding behavior for $`t+\mathrm{}`$. In § 6 we consider general values of the coupling constant $`\kappa `$ and describe the behavior of some solutions of (2) when $`A(x,t)`$ is a one-soliton solution of KdV (1). These calculations indicate the exceptional nature of the two values $`\kappa =1/2`$ and $`\kappa =1`$. In the Appendix, we describe several physical applications of the coupled system (1) and (2) to topics in molecular dynamics, mechanics, soliton theory, and the fluid dynamics of internal waves.

## 2 Exact Solution Formulas for $`\kappa =1/2`$

As is well-known \[GGKM67\], the KdV equation (1) is the compatibility condition for a pair of linear equations involving a complex parameter $`\lambda `$ for an auxiliary function $`f(x,t,\lambda )`$. This pair of linear equations is

$$_x^2f=\frac{\lambda ^2}{4}f\frac{1}{6}Af,\text{ and }_tf=\frac{1}{6}_xAf+\left(\lambda ^2\frac{1}{3}A\right)_xf,$$
(6)
and is called a Lax pair. A simultaneous solution $`f(x,t,\lambda )`$ of these linear equations exists if and only if the function $`A(x,t)`$ satisfies KdV (1). Suppose that this is the case. Then, it is a direct matter to verify that for fixed but arbitrary $`\lambda `$, the two functions defined by

$$B(x,t):=_x\left[f(x,t,\lambda )\mathrm{exp}\left(\pm \frac{i}{2}(\lambda x+\lambda ^3t)\right)\right],$$
(7)
are solutions of the linear equation (2) when $`\kappa =1/2`$. Note here an important point of departure from the other solvable case, namely $`\kappa =1`$, where (2) is the linearized KdV equation. In the latter case, particular solutions are expressed in terms of the $`x`$-derivative of the square of the Lax eigenfunction $`f(x,t,\lambda )`$ \[GGKM74, S83\]. By contrast, the formula (7) for solutions of (2) for $`\kappa =1/2`$ is linear in $`f(x,t,\lambda )`$. This fact leads to some important simplifications.

The formula (7) is only really practical to use if one can explicitly compute the function $`f(x,t,\lambda )`$. This will be the case if the solution $`A(x,t)`$ of KdV (1) is a pure $`N`$-soliton solution. For each fixed $`t`$, $`A(x,t)`$ is then a reflectionless potential of the Schrödinger equation in the Lax pair (6). The multisoliton solutions of KdV and the associated solutions of the Lax pair are constructed as follows \[KM56\]. Let $`f_+(x,t,\lambda )`$ be given by

$$f_+(x,t,\lambda ):=\left(1+\underset{n=0}{\overset{N1}{}}\lambda ^{nN}f_n(x,t)\right)\mathrm{exp}\left(\frac{i}{2}(\lambda x+\lambda ^3t)\right),$$
(8)
where the $`f_n(x,t)`$ are unknown coefficients. Choose $`N`$ positive numbers $`\eta _1>\eta _2>\mathrm{}>\eta _N`$, and $`N`$ arbitrary real numbers $`\alpha _1,\mathrm{},\alpha _N`$, and insist that $`f_+(x,t,\lambda )`$ satisfy the relations

$$f_+(x,t,2i\eta _n)=(1)^{n+1}\mathrm{exp}(2\eta _n\alpha _n)f_+(x,t,2i\eta _n),$$
(9)
for all $`n=1,\mathrm{},N`$. It is easy to see that these relations imply a square linear algebraic system for the coefficients $`f_n(x,t)`$. The determinant of the system is always nonzero, and so the coefficients $`f_n(x,t)`$ are determined uniquely from the soliton eigenvalues $`\{\lambda _n=2i\eta _n\}`$ and the norming constants $`\{\alpha _n\}`$ in terms of exponential functions. From this construction, it can be shown that if one chooses

$$A(x,t):=6i_xf_{N1}(x,t),$$
(10)
then $`f_+(x,t,\lambda )`$ and $`f_{}(x,t,\lambda ):=f_+(x,t,\lambda )`$ are two simultaneous solutions of the Lax pair (6), and the function $`A(x,t)`$ defined by (10) satisfies KdV (1). The two functions $`f_\pm (x,t,\lambda )`$ are linearly independent for all nonzero $`\lambda \pm 2i\eta _n`$. According to the linear relations (9) that determine the coefficients, at the exceptional values of $`\lambda `$ the two functions are proportional.

The solution $`A(x,t)`$ of KdV so constructed represents the interaction of $`N`$ solitons. In particular, as $`t\pm \mathrm{}`$, the solution can be represented in the form

$$A(x,t)\underset{n=1}{\overset{N}{}}A_n^\pm (x,t),\text{ where }A_n^\pm (x,t):=12\eta _n^2\mathrm{sech}^2(\eta _n(x\alpha _n^\pm )4\eta _n^3t),$$
(11)
where the asymptotic phase constants $`\alpha _n^\pm `$ are functions of the $`\eta _n`$ and $`\alpha _n`$.

Below, we will need the asymptotic behavior of the functions $`f_\pm (x,t,\lambda )`$ as $`x\pm \mathrm{}`$ for $`\lambda `$ and $`t`$ fixed. It can be shown that the coefficient functions $`f_n(x,t)`$ remain bounded as $`x\pm \mathrm{}`$. Then, letting $`x`$ tend to $`\pm \mathrm{}`$ in the linear relations (9), one finds from dominant balance arguments that

$$\underset{x\pm \mathrm{}}{lim}\left(1+\underset{n=0}{\overset{N1}{}}\lambda ^{nN}f_n(x,t)\right)|_{\lambda =\pm 2i\eta _n}=0.$$
(12)
These relations imply that

$$\underset{x\pm \mathrm{}}{lim}\left(1+\underset{n=0}{\overset{N1}{}}(\pm \lambda )^{nN}f_n(x,t)\right)=\lambda ^N\underset{n=1}{\overset{N}{}}(\lambda 2i\eta _n).$$
(13)
Therefore, for all $`\lambda `$,

$$\underset{x\pm \mathrm{}}{lim}f_+(x,t,\lambda )\mathrm{exp}\left(\frac{i}{2}(\lambda x+\lambda ^3t)\right)=\lambda ^N\underset{n=1}{\overset{N}{}}(\lambda 2i\eta _n).$$
(14)
The large $`|x|`$ asymptotics for the other solution $`f_{}(x,t,\lambda )`$ follow from the definition $`f_{}(x,t,\lambda )=f_+(x,t,\lambda )`$. From these asymptotics, it is easy to construct the appropriate linear combinations of $`f_\pm (x,t,\lambda )`$ that correspond to the Jost functions of the Schrödinger equation, normalized at $`x=\pm \mathrm{}`$.

In the formula (7) we have a choice of sign in the exponent. In fact, it is easy to see that if one considers the totality of solutions obtained for all complex $`\lambda `$, the choice of sign is redundant. In what follows, we adopt a particular choice of the sign and maintain generality by using both Lax eigenfunctions $`f_+(x,t,\lambda )`$ and $`f_{}(x,t,\lambda )`$. Thus, the particular solutions of the linear equation (2) for $`\kappa =1/2`$ that we will consider below will be denoted by $`h_\pm (x,t,\lambda )`$, given by

$$h_\pm (x,t,\lambda ):=_xg_\pm (x,t,\lambda ),\text{ where }g_\pm (x,t,\lambda ):=f_\pm (x,t,\lambda )\mathrm{exp}\left(\frac{i}{2}(\lambda x+\lambda ^3t)\right).$$
(15)
From the Schrödinger equation (6) for $`f_\pm (x,t,\lambda )`$ it follows that $`g_\pm (x,t,\lambda )`$ satisfies the ODE

$$i_x^2g_\pm i\frac{A}{6}g_\pm =\lambda _xg_\pm .$$
(16)
This ODE plays an important role in suggesting the completeness relation for the solutions $`h_\pm (x,t,\lambda )`$.

## 3 The Completeness Relation for $`\kappa =1/2`$

Having in hand a large family of exact solutions of the linear equation (2) for $`\kappa =1/2`$ is certainly useful, but we may then ask whether there are enough of these solutions to construct the general solution of the initial value problem by superposition. A completeness relation is a formula that gives the expansion of arbitrary initial data in terms of such a collection of functions. In this section, we will establish the completeness relation for the exact solutions $`h_\pm (x,t,\lambda )`$ obtained in § 2.

The form of the completeness relation is suggested by a similar argument to that used by Sachs \[S83\] in his investigation of the completeness of squared eigenfunction solutions to the linearized KdV equation. The idea is that, ideally we would like to have a differential eigenvalue problem in standard form satisfied by the functions $`h_\pm (x,t,\lambda )`$:

$$L(t)h_\pm (x,t,\lambda )=\lambda h_\pm (x,t,\lambda ),$$
(17)
where $`\lambda `$ is the eigenvalue and $`L(t)`$ is some second-order linear differential operator in $`x`$. Then, using the two explicit solutions $`h_\pm (x,t,\lambda )`$ of this problem, we could solve the inhomogeneous problem

$$L(t)\psi \lambda \psi =\varphi $$
(18)
by variation of parameters, i.e. by writing $`\psi `$ as a linear combination of $`h_\pm (x,t,\lambda )`$ with nonconstant coefficients, and substituting into (18). For a fixed function $`\varphi (x)`$, this determines $`\psi (x,t,\lambda )`$, and we have thus constructed the resolvent of the operator $`L(t)`$,

$$\psi (x,t,\lambda )=(L(t)\lambda 𝕀)^1\varphi (x).$$
(19)
If the spectrum of $`L(t)`$ is contained in a bounded region of the complex plane (and also under some milder conditions), then the Dunford-Taylor integral of the resolvent on a positively oriented contour enclosing the spectrum yields the identity operator:

$$\frac{1}{2\pi i}\psi (x,t,\lambda )𝑑\lambda =\frac{1}{2\pi i}(L(t)\lambda 𝕀)^1\varphi (x)𝑑\lambda =\varphi (x).$$
(20)
But we do not have a second-order eigenvalue problem for $`h_\pm (x,t,\lambda )`$. Instead we have the second-order equation (16) for $`g_\pm (x,t,\lambda )`$. However, we make the guess that a similar procedure will apply here. Namely, for appropriate side conditions (see below) we solve the inhomogeneous equation

$$i_x^2\psi i\frac{A}{6}\psi \lambda _x\psi =\varphi ,$$
(21)
for $`\psi (x,t,\lambda )`$ using variation of parameters with the two functions $`g_\pm (x,t,\lambda )`$ solving the homogeneous equation (16), and then we differentiate the resulting formula with respect to $`x`$. Formally speaking only, we have thus constructed the resolvent of the “operator”

$$L(t)=i_xi\frac{A}{6}_x^1.$$
(22)
The obstruction to rigor here is that $`_x^1`$ is not well-defined. Nonetheless, we are guided to hypothesize that

$$\frac{1}{2\pi i}_x\psi (x,t,\lambda )d\lambda =\varphi (x),$$
(23)
for an appropriate contour of integration. This formula turns out to be correct, although a direct proof must be supplied. The proof we use follows Miller and Akhmediev \[MA98\].

### 3.1 Solving the inhomogeneous problem.

We express the solution of the inhomogeneous problem in the form

$$\psi (x,t,\lambda )=C_+(x,t,\lambda )g_+(x,t,\lambda )+C_{}(x,t,\lambda )g_{}(x,t,\lambda ),$$
(24)
subject to the usual “reduction of order” condition

$$_xC_+(x,t,\lambda )g_+(x,t,\lambda )+_xC_{}(x,t,\lambda )g_{}(x,t,\lambda )=0.$$
(25)
Substituting (24) into the equation for $`\psi `$, and using (25), one finds

$$_xC_+(x,t,\lambda )=i\frac{\varphi (x)g_{}(x,t,\lambda )}{W(g_+,g_{})},\text{ and }_xC_{}(x,t,\lambda )=i\frac{\varphi (x)g_+(x,t,\lambda )}{W(g_+,g_{})},$$
(26)
where $`W(g_+,g_{}):=g_+_xg_{}g_{}_xg_+`$ is the Wronskian.

From the differential equation (16) satisfied by $`g_\pm (x,t,\lambda )`$, it follows that

$$_xW(g_+,g_{})=i\lambda W(g_+,g_{}).$$
(27)
Using the large $`|x|`$ asymptotics of $`f_\pm (x,t,\lambda )`$ obtained in § 2 one then solves (27) uniquely and finds that

$$W(g_+,g_{})=i\lambda ^{12N}\mathrm{exp}(i(\lambda x+\lambda ^3t))\underset{n=1}{\overset{N}{}}(\lambda ^2+4\eta _n^2).$$
(28)
In solving the inhomogeneous equation (21) for $`\psi `$, we really should impose appropriate side conditions. Here, the side conditions we use are not related to boundary conditions in $`x`$ as much as analyticity conditions in $`\lambda `$. It is easy to check that for each $`x_{0,U}`$, the function $`\psi _U(x,t,\lambda )`$ defined by

$$\begin{array}{c}\psi _U(x,t,\lambda )=_{x_{0,U}}^x\frac{g_{}(z,t,\lambda )\mathrm{exp}(i(\lambda z+\lambda ^3t))\varphi (z)}{\lambda ^{12N}{\displaystyle \underset{n=1}{\overset{N}{}}}(\lambda ^2+4\eta _n^2)}𝑑zg_+(x,t,\lambda )\hfill \\ \\ +_{\mathrm{}}^x\frac{g_+(z,t,\lambda )\mathrm{exp}(i(\lambda z+\lambda ^3t))\varphi (z)}{\lambda ^{12N}{\displaystyle \underset{n=1}{\overset{N}{}}}(\lambda ^2+4\eta _n^2)}𝑑zg_{}(x,t,\lambda ),\hfill \end{array}$$
(29)
is a solution analytic in $`\lambda `$ for $`\mathrm{}(\lambda )>0`$ and $`|\lambda |`$ sufficiently large. Similarly, $`\psi _L(x,t,\lambda )`$ defined for each $`x_{0,L}`$ by

$$\begin{array}{c}\psi _L(x,t,\lambda )=_{x_{0,L}}^x\frac{g_{}(z,t,\lambda )\mathrm{exp}(i(\lambda z+\lambda ^3t))\varphi (z)}{\lambda ^{12N}{\displaystyle \underset{n=1}{\overset{N}{}}}(\lambda ^2+4\eta _n^2)}𝑑zg_+(x,t,\lambda )\hfill \\ \\ _x^{\mathrm{}}\frac{g_+(z,t,\lambda )\mathrm{exp}(i(\lambda z+\lambda ^3t))\varphi (z)}{\lambda ^{12N}{\displaystyle \underset{n=1}{\overset{N}{}}}(\lambda ^2+4\eta _n^2)}𝑑zg_{}(x,t,\lambda ),\hfill \end{array}$$
(30)
is a solution analytic for $`\mathrm{}(\lambda )<0`$ and $`|\lambda |`$ sufficiently large. The qualification of $`|\lambda |`$ being sufficiently large is necessary because the expressions have poles at the soliton eigenvalues in the respective half-planes where the two functions $`g_\pm (x,t,\lambda )`$ become proportional. However, these are the only finite singularities and both solutions $`\psi _U(x,t,\lambda )`$ and $`\psi _L(x,t,\lambda )`$ are meromorphic in the whole of their respective open half-planes.

The arbitrariness of the parameters $`x_{0,U}`$ and $`x_{0,L}`$ would seem to be a problem; but it will turn out that these terms contribute nothing to the Dunford-Taylor integral that we will prove gives the required completeness relation.

### 3.2 Integrating the resolvent.

Here we show that the guess we made is indeed correct.

###### Theorem 1

Let $`\varphi (x)`$ be an absolutely continuous function in $`L^1()`$. Let $`x_{0,U}`$ and $`x_{0,L}`$ be constants and let $`t`$ be fixed. Then,

$$\varphi (x)=\frac{1}{2\pi i}\underset{R\mathrm{}}{lim}\left[_{C_U}_x\psi _U(x,t,\lambda )d\lambda +_{C_L}_x\psi _L(x,t,\lambda )d\lambda \right],$$
(31)
where $`C_U`$ is the positively oriented half-circle from $`R`$ to $`R`$ in the upper half-plane, and $`C_L`$ is the positively oriented half-circle from $`R`$ to $`R`$ in the lower half-plane.

Proof: First, we show that the terms depending on the arbitrary parameters $`x_{0,U}`$ and $`x_{0,L}`$ converge to zero as $`R\mathrm{}`$. This will justify calling the function $`\psi _U`$ or $`\psi _L`$ a “resolvent” even though the inverse is not unique. Consider the integral

$$\begin{array}{ccc}\hfill J_U& :=& _{C_U}_x\left[_{x_{0,U}}^x\frac{g_{}(z,t,\lambda )\mathrm{exp}(i(\lambda z+\lambda ^3t))\varphi (z)}{\lambda ^{12N}{\displaystyle \underset{n=1}{\overset{N}{}}}(\lambda ^2+4\eta _n^2)}𝑑zg_+(x,t,\lambda )\right]d\lambda \hfill \\ & & \\ & =& _{C_U}\frac{\lambda ^{2N}h_+(x,t,\lambda )}{\lambda {\displaystyle \underset{n=1}{\overset{N}{}}}(\lambda ^2+4\eta _n^2)}_{x_{0,U}}^xg_{}(z,t,\lambda )\mathrm{exp}(i(\lambda z+\lambda ^3t))\varphi (z)𝑑z𝑑\lambda ,\hfill \end{array}$$
(32)
where we have used the relation (25). Recall that

$$\begin{array}{ccc}\hfill h_+(x,t,\lambda )& =& \underset{n=1}{\overset{N}{}}\lambda ^n_xf_{Nn}(x,t),\hfill \\ & & \\ \hfill g_{}(z,t,\lambda )\mathrm{exp}(i(\lambda z+\lambda ^3t))& =& 1+\underset{n=1}{\overset{N}{}}(\lambda )^nf_{Nn}(z,t).\hfill \end{array}$$
(33)
Therefore, for all $`\lambda `$ with $`|\lambda |=R>1`$,

$$|h_+(x,t,\lambda )|\frac{1}{R}\underset{n=1}{\overset{N}{}}|_xf_{Nn}(x,t)|,$$
(34)
and

$$\underset{z}{sup}|g_{}(z,t,\lambda )\mathrm{exp}(i(\lambda z+\lambda ^3t))|1+\underset{z}{sup}\underset{n=1}{\overset{N}{}}|f_{Nn}(z,t)|.$$
(35)
This latter relation assumes the uniform boundedness of the functions $`f_k(z,t)`$ in $`z`$. Finally, it is clear that for $`|\lambda |=R>sup_n2\eta _n`$,

$$\left|\lambda ^{12N}\underset{n=1}{\overset{N}{}}(\lambda ^2+4\eta _n^2)\right|R\underset{n=1}{\overset{N}{}}\left(1\frac{4\eta _n^2}{R^2}\right),.$$
(36)
It follows that for all $`\lambda `$ with $`|\lambda |=R`$ sufficiently large,

$$|J_U|\frac{K(x,t)}{R}\varphi _1,$$
(37)
where

$$K(x,t)=\pi \underset{n=1}{\overset{N}{}}\left(1\frac{4\eta _n^2}{R^2}\right)^1\left(\underset{n=1}{\overset{N}{}}|_xf_{Nn}(x,t)|\right)\left(1+\underset{z}{sup}\underset{n=1}{\overset{N}{}}|f_{Nn}(z,t)|\right).$$
(38)
The bound (37) clearly vanishes as $`R\mathrm{}`$. A nearly identical argument shows that the integral

$$\begin{array}{ccc}\hfill J_L& :=& _{C_L}_x\left[_{x_{0,L}}^x\frac{g_{}(z,t,\lambda )\mathrm{exp}(i(\lambda z+\lambda ^3t))\varphi (z)}{\lambda ^{12N}{\displaystyle \underset{n=1}{\overset{N}{}}}(\lambda ^2+4\eta _n^2)}𝑑zg_+(x,t,\lambda )\right]d\lambda \hfill \end{array}$$
(39)
satisfies the same bound (37) as $`J_U`$.

Now we consider integrating the second terms of $`_x\psi _U(x,t,\lambda )`$ and $`_x\psi _L(x,t,\lambda )`$ respectively. For brevity, define

$$Y(x,z,t,\lambda ):=\frac{g_+(z,t,\lambda )\mathrm{exp}(i(\lambda z+\lambda ^3t))h_{}(x,t,\lambda )}{\lambda ^{12N}{\displaystyle \underset{n=1}{\overset{N}{}}}(\lambda ^2+4\eta _n^2)}.$$
(40)
Note that this can be written as

$$\begin{array}{c}Y(x,z,t,\lambda )=\mathrm{exp}(i\lambda (xz))\frac{\left(1+{\displaystyle \underset{n=1}{\overset{N}{}}}{\displaystyle \frac{f_{Nn}(z,t)}{\lambda ^n}}\right)}{\lambda {\displaystyle \underset{n=1}{\overset{N}{}}}\left(1+{\displaystyle \frac{4\eta _n^2}{\lambda ^2}}\right)}\hfill \\ \\ \times \left(i\lambda \left(1+\underset{n=1}{\overset{N}{}}\frac{f_{Nn}(x,t)}{(\lambda )^n}\right)+\underset{n=1}{\overset{N}{}}\frac{_xf_{Nn}(x,t)}{(\lambda )^n}\right),\hfill \end{array}$$
(41)
and therefore,

$$Y(x,z,t,\lambda )=i\mathrm{exp}(i\lambda (xz))(1+\mathrm{\Delta }(x,z,t,\lambda ))),$$
(42)
where $`\mathrm{\Delta }(x,z,t,\lambda )=O(\lambda ^1)`$ uniformly in $`x`$ and $`z`$ for fixed $`t`$. It also follows from additional cancellation that for $`z=x`$, $`\mathrm{\Delta }(x,x,t,\lambda )=O(\lambda ^2)`$ uniformly in $`x`$. Finally, derivatives of $`\mathrm{\Delta }`$ are controlled as well: $`_z\mathrm{\Delta }(x,z,t,\lambda )=O(\lambda ^1)`$ uniformly. The integral we need to compute for the contribution of $`_x\psi _U(x,t,\lambda )`$ is

$$\begin{array}{c}_{C_R}_{\mathrm{}}^xY(x,z,t,\lambda )\varphi (z)𝑑z𝑑\lambda =i_{C_R}_{\mathrm{}}^x\mathrm{exp}(i\lambda (xz))\varphi (z)𝑑z𝑑\lambda \hfill \\ \\ +i_{C_R}_{\mathrm{}}^x\mathrm{\Delta }(x,z,t,\lambda )\mathrm{exp}(i\lambda (xz))\varphi (z)𝑑z𝑑\lambda .\hfill \end{array}$$
(43)
Note that since the integrand is analytic in the upper half-plane, the first term can be written as:

$$i_{C_R}_{\mathrm{}}^x\mathrm{exp}(i\lambda (xz))\varphi (z)𝑑z𝑑\lambda =i_R^R_{\mathrm{}}^x\mathrm{exp}(i\lambda (xz))\varphi (z)𝑑z𝑑\lambda .$$
(44)
In order to control the error term, it is necessary to integrate by parts once:

$$\begin{array}{c}i_{C_R}_{\mathrm{}}^x\mathrm{\Delta }(x,z,t,\lambda )\mathrm{exp}(i\lambda (xz))\varphi (z)𝑑z𝑑\lambda =\varphi (x)_{C_R}\frac{\mathrm{\Delta }(x,x,t,\lambda )}{\lambda }𝑑\lambda \hfill \\ \\ +_{C_R}_{\mathrm{}}^x\frac{\mathrm{exp}(i\lambda (xz))}{\lambda }_z(\mathrm{\Delta }(x,z,t,\lambda )\varphi (z))dzd\lambda .\hfill \end{array}$$
(45)
The boundary term at $`z=\mathrm{}`$ vanishes because $`\mathrm{\Delta }(x,z,t,\lambda )`$ is bounded there, $`\varphi `$ is continuous and integrable, and $`\mathrm{exp}(i\lambda (xz))`$ is exponentially small for $`\mathrm{}(\lambda )>0`$. Since the exponential is bounded in magnitude by unity for $`\mathrm{}(\lambda )>0`$, and the contour is of length $`\pi R`$, the above estimates of $`\mathrm{\Delta }`$ imply that there exist $`K_0(x,t)`$, $`K_1(x,t)`$ and $`K_2(x,t)`$ all positive, such that

$$\begin{array}{c}\left|_{C_R}_{\mathrm{}}^x\mathrm{\Delta }(x,z,t,\lambda )\mathrm{exp}(i\lambda (xz))\varphi (z)𝑑z𝑑\lambda \right|=\hfill \\ \\ \frac{K_0(x,t)}{R^2}|\varphi (x)|+\frac{K_1(x,t)}{R}\varphi _1+\frac{K_2(x,t)}{R}\varphi ^{}_1.\hfill \end{array}$$
(46)
This proves that

$$\underset{R\mathrm{}}{lim}_{C_R}_{\mathrm{}}^xY(x,z,t,\lambda )\varphi (z)𝑑z𝑑\lambda =i\underset{R\mathrm{}}{lim}_R^R_{\mathrm{}}^x\mathrm{exp}(i\lambda (xz))\varphi (z)𝑑z𝑑\lambda .$$
(47)
Similar arguments applied to the contribution of $`_x\psi _L(x,t,\lambda )`$ show that

$$\underset{R\mathrm{}}{lim}_{C_L}_x^{\mathrm{}}Y(x,z,t,\lambda )\varphi (z)𝑑z𝑑\lambda =i\underset{R\mathrm{}}{lim}_R^R_x^{\mathrm{}}\mathrm{exp}(i\lambda (xz))\varphi (z)𝑑z𝑑\lambda ,$$
(48)
and therefore

$$\begin{array}{c}\frac{1}{2\pi i}\underset{R\mathrm{}}{lim}\left[_{C_R}_x\psi _U(x,t,\lambda )d\lambda +_{C_L}_x\psi _L(x,t,\lambda )d\lambda \right]=\hfill \\ \\ \frac{1}{2\pi }\underset{R\mathrm{}}{lim}_R^R_{\mathrm{}}^{\mathrm{}}\mathrm{exp}(i\lambda (xz))\varphi (z)𝑑z𝑑\lambda =\varphi (x),\hfill \end{array}$$
(49)
with the last equality following from Fourier inversion. This establishes (31) and the theorem. $`\mathrm{}`$

As it stands, the completeness relation given in Theorem 1 is not really an expansion of $`\varphi (x)`$ in terms of the functions $`h_{}(x,t,\lambda )`$, because the expansion coefficients themselves depend on $`x`$. This is easily remedied by casting the right-hand side of the completeness relation into a more useful form. We do this now.

###### Theorem 2

Let $`\varphi (x)`$ be an absolutely continuous function in $`L^1()`$. Let $`t`$ be fixed, and choose any $`w\{\mathrm{},+\mathrm{}\}`$. Define the mode function

$$H(x,t,\lambda ):=\lambda ^Nh_{}(x,t,\lambda ),$$
(50)
which is an entire function of $`\lambda `$, and the amplitudes

$$\begin{array}{ccc}\hfill b^+(t,\lambda )& :=& _w^{\mathrm{}}\frac{\lambda ^Ng_+(z,t,\lambda )\mathrm{exp}(i(\lambda z+\lambda ^3t))}{\lambda {\displaystyle \underset{n=1}{\overset{N}{}}}(\lambda ^2+4\eta _n^2)}\varphi (z)𝑑z,\hfill \\ & & \\ \hfill b^{}(t,\lambda )& :=& _{\mathrm{}}^w\frac{\lambda ^Ng_+(z,t,\lambda )\mathrm{exp}(i(\lambda z+\lambda ^3t))}{\lambda {\displaystyle \underset{n=1}{\overset{N}{}}}(\lambda ^2+4\eta _n^2)}\varphi (z)𝑑z,\hfill \end{array}$$
(51)
and set $`b(t,\lambda ):=b^+(t,\lambda )+b^{}(t,\lambda )`$. The amplitudes have simple poles at $`\lambda =0`$ and $`\lambda =\pm 2i\eta _n`$ for $`n=1,\mathrm{},N`$. Finally, set

$$b_0(t):=\frac{1}{2}\begin{array}{c}\\ \mathrm{Res}\\ \lambda =0\end{array}(b^+(t,\lambda )b^{}(t,\lambda )),\text{ and }b_n^\pm (t):=\begin{array}{c}\\ \mathrm{Res}\\ \lambda =\pm 2i\eta _n\end{array}b^{}(t,\lambda ).$$
(52)
Then we have the expansion

$$\begin{array}{c}\varphi (x)=\underset{R\mathrm{}}{lim}\frac{1}{2\pi i}\mathrm{P}.\mathrm{V}._R^Rb(t,\lambda )H(x,t,\lambda )𝑑\lambda \hfill \\ \\ +b_0(t)H(x,t,0)+\underset{n=1}{\overset{N}{}}\left[b_n^{}(t)H(x,t,2i\eta _n)+b_n^+(t)H(x,t,2i\eta _n)\right].\hfill \end{array}$$
(53)
Remark: Since $`w`$ is now fixed and not a function of $`x`$, this expansion (53) is a true completeness relation, expressing an arbitrary given function $`\varphi (x)`$ as a sum of known functions $`H(x,t,\lambda )`$. From the exact formulas (50), (15), and (8), it is clear that the part of the expansion (53) represented by the singular integral is Fourier-like, with the corresponding components of the solution, $`H(x,t,\lambda )`$ for $`\lambda `$ being bounded oscillatory functions tending to complex exponentials for large $`x`$. On the other hand, the discrete contributions to the solution represent bound states. The $`2N+1`$ bound state terms in (53) are not linearly independent. From the fact that at the eigenvalues $`\pm 2i\eta _n`$ the functions $`g_{}(x,t,\lambda )`$ are all linear combinations of the same $`N`$ functions $`f_0(x,t),\mathrm{},f_{N1}(x,t)`$, it is clear that only $`N`$ of the bound states are linearly independent. These facts are easiest to see when one takes $`w`$ to $`\mathrm{}`$ or $`\mathrm{}`$. Then, half of the contributions from the eigenvalues disappear, and it only remains to express the bound state at zero, $`H(x,t,0)`$, in terms of $`H(x,t,\pm 2i\eta _n)`$. This can be done directly. From the exact formulas (50), (15) and (8), we see that

$$H(x,t,0)=(1)^N_xf_0(x,t),$$
(54)
and making use of the relations (9) satisfied by $`f_+`$ at the eigenvalues,

$$H(x,t,2i\eta _n)=(1)^{n+1}\mathrm{exp}(2\eta _n\alpha _n)\underset{p=0}{\overset{N1}{}}(2i\eta _n)^p_xf_p(x,t).$$
(55)
Expressing $`H(x,t,0)`$ in terms of $`H(x,t,2i\eta _n)`$ is therefore a polynomial interpolation problem. Introduce the polynomial

$$P(\lambda )=\underset{p=0}{\overset{N1}{}}_xf_p(x,t)\lambda ^p.$$
(56)
Given isolated values of this polynomial:

$$P(2i\eta _n)=(1)^{n+1}\mathrm{exp}(2\eta _n\alpha _n)H(x,t,2i\eta _n),\text{ for }n=1,\mathrm{},N,$$
(57)
we are to find $`P(0)`$ and thus $`H(x,t,0)=(1)^N_xf_0(x,t)=(1)^NP(0)`$. Expressing $`P(\lambda )`$ explicitly in terms of Lagrange polynomials gives

$$P(\lambda )=\underset{n=1}{\overset{N}{}}(1)^{n+1}\mathrm{exp}(2\eta _n\alpha _n)H(x,t,2i\eta _n)\underset{kn}{}\frac{\lambda 2i\eta _k}{2i\eta _n2i\eta _k},$$
(58)
and therefore

$$H(x,t,0)=\underset{n=1}{\overset{N}{}}\left[(1)^n\mathrm{exp}(2\eta _n\alpha _n)\underset{kn}{}\frac{\eta _k}{\eta _n\eta _k}\right]H(x,t,2i\eta _n).$$
(59)
There is of course a similar expression for $`H(x,t,0)`$ in terms of $`H(x,t,2i\eta _n)`$. $`\mathrm{}`$

Remark: An important distinction between the completeness relation stated in Theorem 53 and that found by Sachs \[S83\] for derivatives of squared Schrödinger eigenfunctions is in the nature of the singularity at $`\lambda =0`$. Sachs shows that in the expansion of $`\varphi `$ in terms of derivatives of squared eigenfunctions, there is an apparent singularity at $`\lambda =0`$ that it is in fact removable. On the other hand, the integral in Theorem 53 is essentially singular and the residue contribution of $`b_0(t)`$ is nonzero. $`\mathrm{}`$

Proof of Theorem 53: We first establish that in the formulas (29) for $`_x\psi _U(x,t,\lambda )`$ and (30) for $`_x\psi _L(x,t,\lambda )`$ we may replace $`x`$ in the limits of integration by any other value without changing the result of the theorem. That is, we will now show that the integral

$$_{C_U}_{\mathrm{}}^wY(x,z,t,\lambda )\varphi (z)𝑑z𝑑\lambda _{C_L}_w^{\mathrm{}}Y(x,z,t,\lambda )\varphi (z)𝑑z𝑑\lambda ,$$
(60)
which we have already seen converges as $`R`$ tends to infinity to $`2\pi i\varphi (x)`$ in the case that $`w=x`$, is in fact independent of $`w`$. Holding $`R`$ fixed and differentiating with respect to $`w`$, we must show that for sufficiently large $`R`$,

$$\varphi (w)_{|\lambda |=R}Y(x,w,t,\lambda )𝑑\lambda 0,$$
(61)
identically in $`x`$, $`w`$, and $`t`$. Being as the integrand is meromorphic in the finite $`\lambda `$ plane, we can evaluate the integral by residues. There are simple poles at $`\lambda =0`$ and $`\lambda =\pm 2i\eta _n`$ for $`n=1,\mathrm{},N`$. Using the linear relations (9) satisfied by $`f_\pm `$ at the eigenvalues $`\lambda =2i\eta _n`$, we find

$$\begin{array}{ccc}\hfill \begin{array}{c}\\ \mathrm{Res}\\ \lambda =2i\eta _k\end{array}Y(x,w,t,\lambda )& =& \underset{p=0}{\overset{N1}{}}\frac{(2i\eta _k)^{N+p1}}{D_k}_xf_p(x,t)\hfill \\ & & \\ & & +\underset{p,q=0}{\overset{N1}{}}(1)^{Nq}\frac{(2i\eta _k)^{p+q1}}{D_k}f_q(w,t)_xf_p(x,t),\hfill \\ & & \\ \hfill \begin{array}{c}\\ \mathrm{Res}\\ \lambda =2i\eta _k\end{array}Y(x,w,t,\lambda )& =& \underset{p=0}{\overset{N1}{}}(1)^{Np}\frac{(2i\eta _k)^{N+p1}}{D_k}_xf_p(x,t)\hfill \\ & & \\ & & +\underset{p,q=0}{\overset{N1}{}}(1)^{Np}\frac{(2i\eta _k)^{p+q1}}{D_k}f_q(w,t)_xf_p(x,t),\hfill \end{array}$$
(62)
where

$$D_k:=\underset{nk}{}(2i\eta _k2i\eta _n)\underset{n=1}{\overset{N}{}}(2i\eta _k+2i\eta _n).$$
(63)
Similarly, for the residue at zero,

$$\begin{array}{c}\\ \mathrm{Res}\\ \lambda =0\end{array}Y(x,w,t,\lambda )=(1)^N\frac{f_0(w,t)_xf_0(x,t)}{{\displaystyle \underset{n=1}{\overset{N}{}}}4\eta _n^2}.$$
(64)
Adding all the residues and collecting coefficients of the terms $`_xf_p(x,t)`$ and $`f_q(w,t)_xf_p(x,t)`$, we find that the sum of the residues will be zero if:

$$I_p:=\underset{k=1}{\overset{N}{}}\frac{(2i\eta _k)^p}{D_k}=0$$
(65)
for all odd $`p=1,3,5,\mathrm{},2N3`$, and if:

$$\frac{1}{{\displaystyle \underset{n=1}{\overset{N}{}}}4\eta _n^2}+2\underset{k=1}{\overset{N}{}}\frac{1}{2i\eta _kD_k}=0.$$
(66)
These expressions are themselves sums of residues of meromorphic differentials. Thus, by inspection, one finds that for $`p=1,3,5,\mathrm{},2N3`$,

$$I_p=\frac{1}{2\pi i}_C\frac{\lambda ^pd\lambda }{{\displaystyle \underset{n=1}{\overset{N}{}}}(\lambda ^2+4\eta _n^2)},$$
(67)
where $`C`$ is any simple counter-clockwise oriented contour that encircles the points $`\lambda =2i\eta _n`$ for $`n=1,\mathrm{},N`$ (but without enclosing the conjugate eigenvalues or $`\lambda =0`$). With $`p`$ bounded by $`2N3`$, the path of integration can be blown out to infinity in the upper half-plane, and then brought down to the real axis so that

$$I_p=\frac{1}{2\pi i}_{\mathrm{}}^{\mathrm{}}\frac{\lambda ^pd\lambda }{{\displaystyle \underset{n=1}{\overset{N}{}}}(\lambda ^2+4\eta _n^2)}=0,$$
(68)
with the last equality following from the oddness of the integrand for odd $`p`$. Finally, consider the integral $`I_1`$ defined by

$$I_1:=\frac{1}{2\pi i}_C\frac{d\lambda }{\lambda {\displaystyle \underset{n=1}{\overset{N}{}}}(\lambda ^2+4\eta _n^2)}.$$
(69)
Evaluating the residues inside $`C`$, we find

$$I_1=\underset{k=1}{\overset{N}{}}\frac{1}{2i\eta _k{\displaystyle \underset{nk}{}}(2i\eta _k2i\eta _n){\displaystyle \underset{n=1}{\overset{N}{}}}(2i\eta _k+2i\eta _n)}.$$
(70)
On the other hand, we can again blow the contour $`C`$ out to infinity in the upper half-plane and bring it down to the real axis. This time, there is a singularity at $`\lambda =0`$, so the Plemelj formula must be used. We find

$$I_1=\frac{1}{2}\frac{1}{{\displaystyle \underset{n=1}{\overset{N}{}}}4\eta _n^2}+\frac{1}{2\pi i}\mathrm{P}.\mathrm{V}._{\mathrm{}}^{\mathrm{}}\frac{d\lambda }{\lambda {\displaystyle \underset{n=1}{\overset{N}{}}}(\lambda ^2+4\eta _n^2)}.$$
(71)
Once again, by oddness, the principal value integral vanishes identically, and then combining this result with the previous expression we obtain the required vanishing.

This shows that for any $`w`$,

$$\varphi (x)=\frac{1}{2\pi i}\underset{R\mathrm{}}{lim}\left[_{C_U}_{\mathrm{}}^wY(x,z,t,\lambda )\varphi (z)𝑑z𝑑\lambda _{C_L}_w^{\mathrm{}}Y(x,z,t,\lambda )\varphi (z)𝑑z𝑑\lambda \right].$$
(72)
Establishing (53) and therefore the Theorem now amounts to using the residue theorem once again to deform the integration paths $`C_U`$ and $`C_L`$ in (72) to the real axis. One finds discrete contributions at the poles $`\lambda =\pm 2i\eta _n`$ and then applying the Plemelj formula to contract the contour to the real axis in the neighborhood of $`\lambda =0`$ gives a discrete contribution proportional to $`H(x,t,0)`$ and the principal value regularization of the singular integral over the continuous spectrum. $`\mathrm{}`$

## 4 Solution of the Initial Value Problem for $`\kappa =1/2`$

It is easy to see that when $`A(x,t)`$ is an $`N`$-soliton solution of KdV (1), one can use the completeness relation to solve the initial value problem:

$$_tB+_x\left[\frac{1}{2}AB+_x^2B\right]=0,B(x,0)=\varphi (x).$$
(73)
Setting $`t=0`$, and picking a convenient value of $`w`$, say, $`w=+\mathrm{}`$, one computes the amplitudes (51) and discrete coefficients (52). Then, because the function $`H(x,t,\lambda )`$ satisfies (2) for $`\kappa =1/2`$ and for each complex $`\lambda `$, the expression

$$\begin{array}{c}B(x,t):=\underset{R\mathrm{}}{lim}\frac{1}{2\pi i}\mathrm{P}.\mathrm{V}._R^Rb(0,\lambda )H(x,t,\lambda )𝑑\lambda \hfill \\ \\ +b_0(0)H(x,t,0)+\underset{n=1}{\overset{N}{}}\left[b_n^{}(0)H(x,t,2i\eta _n)+b_n^+(0)H(x,t,2i\eta _n)\right],\hfill \end{array}$$
(74)
provides the solution of the initial value problem (73), generally in the sense of distributions. That is, $`B(x,0)=\varphi (x)`$ by Theorem 53, and for each test function $`\phi (x,t)`$ that is differentiable in $`t`$ and three times differentiable in $`x`$ and has compact support in $`(x,t)\times _+`$, one shows by exchanging the order of integration that

$$_0^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\left[_t\phi (x,t)+\frac{1}{2}A(x,t)_x\phi (x,t)+_x^3\phi (x,t)\right]B(x,t)𝑑x𝑑t=0.$$
(75)
The solution will be classical in as much as it is possible to differentiate with respect to $`x`$ and $`t`$ under the integral sign in the solution formula (74). This requires additional smoothness and decay assumptions on the initial data $`\varphi (x)`$ that we do not consider here.

## 5 Scattering of Bound States for $`\kappa =1/2`$

Of particular interest in applications is the $`N`$-dimensional (recall that $`A(x,t)`$ is an $`N`$-soliton solution of KdV) subspace of solutions of (2) for $`\kappa =1/2`$ consisting of bound states. This subspace represents linear waves that are trapped by the solitons of the potential $`A(x,t)`$. For large $`|t|`$, these bound state solutions are all confined to the trajectories of the solitons. Therefore, it follows that each bound state $`B(x,t)`$ has two asymptotic representations:

$$B(x,t)\underset{n=1}{\overset{N}{}}\beta _n^\pm A_n^\pm (x,t),t\pm \mathrm{},$$
(76)
for some constants $`\beta _n^\pm `$ depending on $`B(x,t)`$, where $`A_n^\pm (x,t)`$ are defined by (11). Since there are exactly $`N`$ linearly independent bound states, it follows that the constants $`\beta _n^+`$ are completely determined from the constants $`\beta _n^{}`$. In particular, there exists an invertible $`N\times N`$ matrix $`𝐓`$ with entries depending only on the data specifying the $`N`$-soliton solution $`A(x,t)`$, such that

$$\beta _j^+=\underset{k=1}{\overset{N}{}}T_{jk}\beta _k^{}.$$
(77)
The matrix $`𝐓`$ is called the bound state scattering matrix. In this section, we compute the scattering matrix explicitly, and show that its elements only depend on the soliton eigenvalues $`\eta _1,\mathrm{},\eta _N`$.

If $`A(x,t)`$ is an $`N`$-soliton solution of KdV (1), then a family of solutions to (2) for $`\kappa =1/2`$, parametrized by complex $`\lambda `$, is given by

$$h_+(x,t,\lambda )=\frac{A(x,t)}{6i\lambda }+\underset{n=0}{\overset{N2}{}}\lambda ^{nN}_xf_n(x,t).$$
(78)
We want to analyze these solutions in the limit of large $`|t|`$, in a frame of reference traveling with constant velocity $`c`$.

The first step is to see how the coefficients $`f_n(x,t)`$ behave for large $`|t|`$. Let $`\chi =xct`$ be fixed as $`\tau =t`$ goes to either $`+\mathrm{}`$ or $`\mathrm{}`$. Begin by taking $`\eta _m^2<4c<\eta _{m1}^2`$ to see how the coefficients behave in between the solitons. In the limit $`\tau +\mathrm{}`$, the equations (9) imply that

$$\begin{array}{ccccccc}\hfill 1& +& \underset{k=0}{\overset{N1}{}}(2i\eta _n)^{kN}f_k& & 0,& n=1,\mathrm{},m1,\hfill & \\ & & & & & & \\ \hfill 1& +& \underset{k=0}{\overset{N1}{}}(2i\eta _n)^{kN}f_k& & 0,& n=m,\mathrm{},N.\hfill & \end{array}$$
(79)
This is an invertible Vandermonde system for the coefficients $`f_k`$, so that as $`\tau +\mathrm{}`$, the $`f_k`$ all become constants, independent of $`\chi `$ and $`\tau `$. Thus, $`_xf_k(x,t)`$ vanishes between the solitons for all $`k`$. The analogous result holds as $`\tau \mathrm{}`$. This shows that the solutions of (2) for $`\kappa =1/2`$ described by the formula (78) are asymptotically confined to the individual frames of reference of the moving solitons in the potential field $`A(x,t)`$.

Now set $`c=4\eta _m^2`$ to go into the moving frame of reference of one of the solitons. Taking the limit $`\tau +\mathrm{}`$ yields

$$\begin{array}{ccccccc}\hfill 1& +& \underset{k=0}{\overset{N1}{}}(2i\eta _n)^{kN}f_k& & 0,& n=1,\mathrm{},m1,\hfill & \\ & & & & & & \\ \hfill 1& +& \underset{k=0}{\overset{N1}{}}(2i\eta _n)^{kN}f_k& & 0,& n=m+1,\mathrm{},N.\hfill & \end{array}$$
(80)
This is a system of $`N1`$ equations in $`N`$ unknowns, so it can be used to asymptotically eliminate $`_\chi f_0`$ through $`_\chi f_{N2}`$ in favor of $`_\chi f_{N1}`$, which we know is proportional to the $`N`$-soliton solution of KdV, $`A(x,t)`$. Thus, as $`\tau +\mathrm{}`$, with $`c=4\eta _m^2`$,

$$_\chi f_k=Q_{mk}_\chi f_{N1}=\frac{1}{6i}Q_{mk}A,$$
(81)
for $`k=0,\mathrm{},N2`$ where the numbers $`Q_{mk}`$ are the unique solution of the inhomogeneous system of linear algebraic equations

$$\begin{array}{ccccccc}\hfill (2i\eta _n)^1& +& \underset{k=0}{\overset{N2}{}}(2i\eta _n)^{kN}Q_{mk}& =& 0,& n=1,\mathrm{},m1,\hfill & \\ & & & & & & \\ \hfill (2i\eta _n)^1& +& \underset{k=0}{\overset{N2}{}}(2i\eta _n)^{kN}Q_{mk}& =& 0,& n=m+1,\mathrm{},N.\hfill & \end{array}$$
(82)
One can similarly show that as $`\tau \mathrm{}`$, with $`c=4\eta _m^2`$,

$$_\chi f_k=Q_{mk}^{}_\chi f_{N1}=\frac{1}{6i}Q_{mk}^{}A,$$
(83)
for $`k=0,\mathrm{},N2`$, where the star denotes complex conjugation.

Now consider particular solutions $`B_j(x,t)`$ of (2) for $`\kappa =1/2`$ obtained as linear combinations of $`N`$ others expressed by the formula (78) evaluated on the $`N`$ soliton eigenvalues. The formula for $`B_j(x,t)`$ is

$$B_j(x,t)=\underset{k=1}{\overset{N}{}}F_{jk}h_+(x,t,2i\eta _k)=\underset{k=1}{\overset{N}{}}F_{jk}\left[\frac{A(x,t)}{12\eta _k}+\underset{n=0}{\overset{N2}{}}(2i\eta _k)^{nN}_xf_n(x,t)\right],$$
(84)
where $`𝐅=\{F_{jk}\}`$ is a matrix of arbitrary constants. From the asymptotics of $`f_n(x,t)`$, we have as $`\tau \mathrm{}`$ with $`c=4\eta _m^2`$,

$$B_jA\underset{k=1}{\overset{N}{}}F_{jk}G_{km}^{},\mathrm{where}G_{km}^{}:=\frac{1}{12\eta _k}+\frac{1}{6i}\underset{n=0}{\overset{N2}{}}(2i\eta _k)^{nN}Q_{mn}^{}.$$
(85)
So, with the choice that the matrix $`\{F_{jk}\}`$ is the inverse of the matrix $`𝐆^{}=\{G_{km}^{}\}`$, the particular solution $`B_j(x,t)`$ of (2) for $`\kappa =1/2`$ will be completely confined as $`t\mathrm{}`$ to the frame of reference moving with speed $`c=4\eta _j^2`$, where it will be locally indistinguishable from the solution $`A(x,t)`$ of KdV. Let us now determine how $`B_j(x,t)`$ will behave in the various soliton frames as $`t+\mathrm{}`$. Passing to the limit of $`\tau +\mathrm{}`$ in the frame with velocity $`c=4\eta _m^2`$ gives

$$B_jA\underset{k=1}{\overset{N}{}}F_{jk}G_{km}^+,\mathrm{where}G_{km}^+:=\frac{1}{12\eta _k}+\frac{1}{6i}\underset{n=0}{\overset{N2}{}}(2i\eta _k)^{nN}Q_{mn}.$$
(86)
These asymptotics give us a formula for the bound state scattering matrix:

$$𝐓:=\left[(𝐆^{})^1𝐆^+\right]^T.$$
(87)
It is clear that the elements of $`𝐓`$ depend only on the $`N`$ soliton eigenvalues $`\eta _1,\mathrm{},\eta _N`$. There is no dependence on the soliton phase variables $`\alpha _1,\mathrm{},\alpha _N`$. Therefore, the asymptotic scattering properties of linear waves in (2) with $`\kappa =1/2`$ are insensitive to phase shifts among the solitons in the potential $`A(x,t)`$. As a concrete example of the scattering matrix, we compute it explicitly for $`N=2`$ for arbitrary $`\eta _1>\eta _2>0`$:

$$𝐓=\frac{1}{\eta _1^2\eta _2^2}\left[\begin{array}{cc}(\eta _1\eta _2)^2& 2\eta _2(\eta _1\eta _2)\\ 2\eta _1(\eta _1\eta _2)& (\eta _1\eta _2)^2\end{array}\right].$$
(88)
The fact that $`T_{22}`$ is negative means that it is possible for the interactions of the solitons in $`A(x,t)`$ to convert trapped linear waves of elevation into waves of depression, and vice-versa.

## 6 General Values of $`\kappa `$

We expect that for most values of $`\kappa `$, the linear waves satisfying (2) will not be permanently trapped by solitons present in the potential $`A(x,t)`$. This is suggested by considering the simplest case, namely taking $`A(x,t)`$ to be the one-soliton solution of KdV (1). The soliton travels with velocity $`c=4\eta ^2`$, so that $`A=V(\chi )`$ with $`\chi =xct\alpha `$. Corresponding traveling wave solutions $`B(\chi )`$ of the linear problem that propagate with the same velocity and decay as $`\chi \pm \mathrm{}`$ satisfy

$$B^{\prime \prime }(\chi )+\kappa V(\chi )B(\chi )=cB(\chi ).$$
(89)
Since $`c`$ is fixed, we can view this as an eigenvalue equation with $`\kappa `$ as the eigenvalue. We therefore expect that only isolated values of $`\kappa `$ will admit nontrivial decaying solutions $`B(\chi )`$. We have already seen that $`\kappa =1/2`$ and $`\kappa =1`$ are indeed eigenvalues. For $`\kappa =1/2`$ the eigenfunction $`B(\chi )`$ is an even function of $`\chi `$, while for $`\kappa =1`$ the eigenfunction $`B(\chi )`$ is odd in $`\chi `$. Since eigenfunctions of (89) must be nondegenerate and therefore have either odd or even parity in $`\chi `$, there cannot exist a nontrivial bound state eigenfunction of (89) for all $`\kappa [1/2,1]`$ since the eigenfunction would have to change parity from one endpoint to the other. Therefore, at least one value of $`\kappa [1/2,1]`$ is not an eigenvalue. For such $`\kappa `$, there is no bound state traveling wave solution of (2) that is trapped in the soliton trajectory.

We can be more precise about this phenomenon. The left-hand side of (89) can also be viewed as a Schrödinger operator $`L(\kappa )`$ depending on a coupling constant $`\kappa `$, and the condition for wave trapping by solitons is simply that $`c\mathrm{\Sigma }_\mathrm{p}(L(\kappa ))`$, where $`\mathrm{\Sigma }_\mathrm{p}`$ denotes the point spectrum. The number of discrete eigenvalues is a nondecreasing function of $`\kappa >0`$, corresponding to the deepening of the potential well. There is an infinite unbounded sequence of cutoff values $`\kappa _n^{\mathrm{cut}}`$ of $`\kappa `$ at which the number of eigenvalues changes by one, and the new eigenvalue is born from the continuum. Each eigenvalue, once born, is distinct and is a decreasing function of $`\kappa `$. From these arguments, it follows that there exists an infinite unbounded sequence of bifurcation values $`\kappa _n^{\mathrm{bif}}`$ of $`\kappa `$ at which one eigenvalue crosses the level $`E=c`$, and a bound state traveling wave solution of (2) exists.

It is easy to find the bifurcation points because the hyperbolic secant squared potential is so well understood. The potential $`\kappa V(\chi )`$ is exactly reflectionless for $`12\kappa =n(n+1)`$ for $`n=1,2,3,\mathrm{}`$. The corresponding energy levels are $`E_{n,k}=k^2\eta ^2`$, for $`k=1,\mathrm{},n`$. Therefore, for $`n>1`$ in this sequence, there is always one eigenvalue that is exactly equal to $`c=4\eta ^2`$. The corresponding eigenstate is always the $`n1`$st state and therefore has $`n2`$ zeros. It follows that the bifurcation points are $`\kappa =\kappa _n^{\mathrm{bif}}=(n+1)(n+2)/12`$ for $`n=1,2,3,\mathrm{}`$.

The fact that some linear waves may be permanently trapped by isolated solitons at a bifurcation point $`\kappa =\kappa _n^{\mathrm{bif}}`$ does not necessarily imply that there will be no losses to radiation when solitons in the field $`A(x,t)`$ interact with one another. Such a lossless interaction might suggest the “integrability” of the linear equation (2). We have indeed seen that this is the case for the first two bifurcation points, $`\kappa =\kappa _1^{\mathrm{bif}}=1/2`$ and $`\kappa =\kappa _2^{\mathrm{bif}}=1`$, but it is by no means clear that the trend continues for higher-order bifurcation points. For the rest of this section, we therefore restrict attention to the case $`N=1`$, that is, we take the nonlinear field $`A(x,t)`$ to be a one-soliton solution of KdV (1).

Using (3) and the change of variables $`x^{}=\eta (x4\eta ^2t\alpha )`$ and $`t^{}=\eta ^3t`$, the equation (2) becomes, after dropping primes,

$$_tB+_x[4B+12\kappa \mathrm{sech}^2(x)B+_x^2B]=0.$$
(90)
This equation is of course solved by separation of variables. We seek separated solutions $`B(x,t)=b_\sigma (x)\mathrm{exp}(\sigma t)`$ and obtain the third-order eigenvalue problem

$$[4b_\sigma (x)12\kappa \mathrm{sech}^2(x)b_\sigma (x)b_\sigma ^{\prime \prime }(x)]^{}=\sigma b_\sigma (x),$$
(91)
where the prime denotes differentiation with respect to $`x`$. In this context, what we have been calling “trapped linear waves” correspond to bound-state eigenfunctions of (91) with $`\sigma =0`$. Such solutions have finite mass and energy and are stationary in the moving frame of reference of the soliton $`A(x,t)`$. As we know, such eigenfunctions with $`\sigma =0`$ exist only at the bifurcation values of $`\kappa =\kappa _n^{\mathrm{bif}}`$. However, it is clear that for general values of $`\kappa `$ there are other possibilities. There may be eigenvalues $`\sigma `$ that are purely imaginary, giving rise to oscillating modes that travel in the soliton frame. More generally, if an eigenvalue has a nonzero real part for some $`\kappa `$, then there will be a mode that is either amplified or exponentially damped as it propagates with the soliton.

The eigenvalue problem (91) has two simple symmetries. Whenever $`b_\sigma (x)`$ is an eigenfunction with eigenvalue $`\sigma `$, then $`b_\sigma (x)`$ is an eigenfunction with eigenvalue $`\sigma `$ and $`b_\sigma (x)^{}`$ is an eigenfunction with eigenvalue $`\sigma ^{}`$. Therefore, the eigenvalues either come in purely real pairs $`(|\sigma |,|\sigma |)`$, purely imaginary pairs $`(i|\sigma |,i|\sigma |)`$, or in complex quartets $`(\sigma ,\sigma ,\sigma ^{},\sigma ^{})`$. These symmetries indicate the distinguished role of $`\sigma =0`$ as a point that if it appears in the spectrum for some $`\kappa `$ can signal a bifurcation in the number of eigenvalues. This explains our terminology and notation for the values $`\kappa =\kappa _n^{\mathrm{bif}}`$.

Most points on the imaginary $`\sigma `$ axis correspond to continuous spectrum. This can be seen by the following argument. Let $`\kappa `$ be fixed. Suppose $`\sigma =i\omega `$ with $`\omega `$. For large $`|x|`$, the solutions of (91) have the form of linear combinations of $`\mathrm{exp}(ik_{\omega ,j}x)`$ where $`k=k_{\omega ,j}`$ are the three roots of $`k^3+4k\omega =0`$. Exactly one of these roots, say $`k=k_{\omega ,0}`$, is real, while the other two form a complex-conjugate pair. If we seek a generalized eigenfunction normalized to $`\mathrm{exp}(ik_{\omega ,0}x)`$ as $`x\mathrm{}`$ through a “shooting” method, we have three complex constants to exploit: the coefficient of the decaying mode for large negative $`x`$, and the coefficients of the decaying mode for large positive $`x`$ and the finite amplitude contribution for large positive $`x`$. Matching the values of $`b_{i\omega }(x)`$, $`b_{i\omega }^{}(x)`$, and $`b_{i\omega }^{\prime \prime }(x)`$ at $`x=0`$ gives three complex equations in three complex unknowns. If this system of equations is solvable at all, one expects it to be solvable for almost all real $`\omega `$, yielding a generalized eigenfunction. For exceptional values of $`\omega `$ where there is not a generalized eigenfunction, there will be a genuine bound-state eigenfunction, since the spectrum is a closed set.

We have used a numerical Fourier-based collocation (pseudospectral) method to find the discrete eigenvalues of (91) over a range of values of the coupling constant $`\kappa `$. Essentially this involves approximating the continuous function $`b_\sigma (x)`$ by a periodic discrete series. Then the derivative $`_x`$ can be approximated to exponential accuracy by a derivative matrix $`𝐃`$, for which an explicit formula is given in \[CHQZ88\]. This is then used to construct a discrete approximation to the operator on the left-hand side of (91). Standard techniques can then be used to obtain the eigenvalues and eigenvectors of this matrix. The corresponding eigenfunctions always decay exponentially, but sometimes they decay very slowly for large $`x`$ of one or the other sign — luckily not both. To obtain accurate results it was necessary in these cases to change the dependent variable by multiplying by an appropriate exponential function of $`x`$ to enhance the decay on the slowly-decaying side without changing decay into growth on the other side. Our results over the range $`0<\kappa <5`$ are shown in Figure 2.

The bifurcation values $`\kappa _n^{\mathrm{bif}}`$ appear to be of two different types. If $`n`$ is odd, then when $`\kappa `$ increases through the value $`\kappa =\kappa _n^{\mathrm{bif}}`$, a new pair of real eigenvalues is born from the origin $`\sigma =0`$. As $`\kappa `$ is further increased, the pair of eigenvalues moves at first away from the imaginary axis and then changes direction and contracts toward the origin. When $`\kappa `$ increases through the even bifurcation value $`\kappa =\kappa _{n+1}^{\mathrm{bif}}`$, the pair enters the origin and re-emerges as a complex eigenvalue quartet. Further increasing the value of $`\kappa `$ causes the quartet of eigenvalues to move through a maximum in the magnitude of the real part and then toward the imaginary axis with the magnitude of the real part decreasing to zero while the magnitude of the imaginary part increases without bound. It does not appear that the quartet of eigenvalues ever disappears into the continuous spectrum, although it comes arbitrarily close as $`\kappa `$ increases. This scenario is repeated again and again as $`\kappa `$ increases through each odd bifurcation value. Representative eigenfunctions are plotted in Figure 3.

Here, one can see that when the eigenvalue $`\sigma `$ has a nonzero imaginary part, the decay of the eigenfunction can be quite slow on the “downstream” side of the soliton. As remarked above, this effect is compensated for in our numerics by working with a modified eigenfunction.

There are no discrete eigenvalues at all for $`\kappa <1/2`$ (and in particular for $`\kappa <0`$), and for all $`\kappa `$ satisfying $`1/2<\kappa <1`$ and $`\kappa >1`$, there is always at least one eigenvalue with a nonzero real part, which corresponds to an exponentially growing eigenfunction and therefore instability. The values $`\kappa =1/2`$ and $`\kappa =1`$ are distinguished as the only values for which there exist discrete eigenvalues and at the same time all eigenvalues have zero real parts, so the system is neutrally stable. For all other values of $`\kappa `$, either there are no discrete eigenvalues at all in which case all initial conditions for (2) disperse away algebraically in time, or there are discrete eigenvalues with positive real parts in which case the linear waves are amplified by the soliton in the field $`A(x,t)`$.

## 7 Conclusion

The coupled system consisting of the KdV equation (1) and the linear equation (2) is integrable for two distinct values of the coupling parameter $`\kappa `$. The integrable case of $`\kappa =1`$ has been studied by other authors \[GGKM74, S83\]. In this paper, we have given new results for the other integrable case, namely $`\kappa =1/2`$. In particular, we have shown how to construct the general solution of (2) for $`\kappa =1/2`$ when the nonlinear field $`A(x,t)`$ is an $`N`$-soliton solution of the KdV equation. This general solution is represented in terms of a number of bound states (equal to the number $`N`$ of solitons in the field $`A(x,t)`$) and a continuum superposition of radiative states given by a singular integral. With the help of numerical computations, we have placed the integrable cases in context by examining the behavior of the linear equation (2) for general values of $`\kappa `$, when $`A(x,t)`$ is a one-soliton solution of KdV. These calculations show that the linear equation (2) behaves as an unstable dynamical system for most positive $`\kappa `$. The integrable value of $`\kappa =1`$, for which the equation (2) is the linearized KdV equation, is an isolated stable point, since a small change of either sign in the value of $`\kappa `$ will lead to the presence of exponentially growing modes. The other integrable value of $`\kappa =1/2`$ represents the boundary between a stable system without any bound states for $`\kappa <1/2`$ and an exponentially unstable system for $`\kappa >1/2`$.

In physical applications of the coupled system (1) and (2) as discussed in the Appendix, the presence of instabilities indicates that more terms need to be included in the model. However, in the stable cases the model is indeed expected to be physically meaningful. And in this regard, the two integrable cases can provide useful starting points for perturbation theory.

As a final remark, let us indicate the kind of calculations that are possible for the coupled system (1) and (2) for $`\kappa =1/2`$ with the aid of the completeness relation. For a family of relevant initial data for the linear equation, one can explicitly compute the projection onto the bound states, and consequently determine the long time behavior of the corresponding solution of (2). Also, the long time behavior of the dispersive part of the solution can be computed from the explicit representation of this component of the solution as a singular integral. We leave such applications of the completeness relation for further investigations.

## Appendix Appendix: Some Applications

It is useful to keep in mind some applications in which the coupled system (1) and (2) might arise. In fact, such equations appear in the modeling of coupling of acoustic phonons in long polymer molecules. Many organic polymers (e.g. DNA and $`\alpha `$-helix proteins like acetanelide) may be considered from the mechanical point of view as long chains of nearly identical masses. This “backbone” of the molecule supports a longitudinal vibrational mode in which the masses are all moving in tandem with zero frequency (i.e. simple translation) in the long-wave limit; the associated quanta are called acoustic phonons. In the presence of intrinsic weak nonlinearity, the Korteweg-de Vries equation describes these vibrations in the long-wave limit. Usually, the masses making up the chain contain internal degrees of freedom (e.g. the “breathing” modes of base-pairs in DNA, and the so-called amide I exciton modes of the C=O bond in each peptide group of an $`\alpha `$-helix protein). The coupling of these internal degrees of freedom to the motion of the backbone leads to a variety of interesting dynamical models (e.g. the discrete sine-Gordon equation for DNA, and the discrete nonlinear Schrödinger equation for $`\alpha `$-helix proteins).

We may consider a situation in which the internal degrees of freedom are themselves acoustic phonons associated with transverse vibrational modes. This can be visualized with the help of a concrete mechanical model, whose equilibrium configuration is shown in Figure 4.

The backbone is made of heavy masses $`M`$ connected by stiff springs. Mounted on each heavy mass is a transversely-oriented frictionless track in which rides a small mass $`m`$. The mass $`M`$ is assumed to include the mass of the track and small mass $`m`$. The small masses are themselves connected by weaker springs. Assigning longitudinal displacements $`u_n`$ to the large masses $`M`$, and transverse displacements $`v_n`$ to the small masses $`m`$ in the frictionless tracks, the Hamiltonian of the mechanical model is:

$$H=\underset{n}{}\left[\frac{1}{2}M\dot{u}_n^2+\frac{1}{2}m\dot{v}_n^2+W(L+u_{n+1}u_n)+V\left(\sqrt{(L+u_{n+1}^2u_n)^2+(v_{n+1}v_n)^2}\right)\right],$$
(92)
where $`W`$ is the potential energy of the stiff springs connecting the large masses and $`V`$ is the potential energy of the weaker springs.

The associated equations of motion are

$$\begin{array}{ccc}\hfill M\ddot{u}_n& =& W^{}(L+u_{n+1}u_n)W^{}(L+u_nu_{n1})\hfill \\ & & \\ & & +S(D_{n+1})(L+u_{n+1}u_n)S(D_n)(L+u_nu_{n1})\hfill \\ & & \\ \hfill m\ddot{v}_n& =& S(D_{n+1})(v_{n+1}v_n)S(D_n)(v_nv_{n1}),\hfill \end{array}$$
(93)
where we have set $`S(D):=V^{}(D)/D`$ and $`D_n:=\sqrt{(L+u_nu_{n1})^2+(v_nv_{n1})^2}`$. It is clear that one may take the undisturbed state of the internal modes $`v_{n+1}=v_n`$ for all $`n`$ to hold exactly, in which case only the backbone motion is relevant. We will be interested in small amplitude, linear motions of the $`v_n`$ and how they are affected by the motion of the backbone.

The disparity between the masses $`M`$ and $`m`$, and that of the strengths of the associated springs, is introduced by letting $`\mu `$ be a small parameter and assuming $`m=\mu M`$ and $`V=\mu U`$ (and correspondingly, $`S=\mu Z`$). We make the small-amplitude long-wave ansatz:

$$u_n(t)=hu(X=hn,T=ht),\text{ and }v_n(t)=hv(X=hn,T=ht),$$
(94)
where $`h`$ is a small lattice-spacing parameter. Expanding the functions $`W`$ and $`Z`$ in Taylor series about the equilibrium position, the equations of motion become

$$\begin{array}{ccc}\hfill Mh^3_T^2u& =& h^3W^{\prime \prime }(L)_X^2u+\frac{h^5}{12}W^{\prime \prime }(L)_X^4u+h^5W^{\prime \prime \prime }(L)_Xu_X^2u+O(h^6)+O(h^3\mu ),\hfill \\ & & \\ \hfill Mh^3_T^2v& =& h^3Z(L)_X^2v+\frac{h^5}{12}Z(L)_X^4v+h^5Z^{}(L)_X^2u_Xv+h^5Z^{}(L)_Xu_X^2v+O(h^6).\hfill \end{array}$$
(95)
Trapping of $`v`$-waves by $`u`$-waves becomes possible if the wave speeds are equal. Therefore, we assume that $`W^{\prime \prime }(L)=Z(L)=Mc^2`$. Changing variables to $`\chi =XcT`$ and $`\tau =h^2T`$ yields

$$\begin{array}{ccc}\hfill Mh^2_\tau ^2u2Mc_\chi _\tau u& =& \frac{Mc^2}{12}_\chi ^4u+W^{\prime \prime \prime }(L)_\chi u_\chi ^2u+O(h)+O(\mu /h^2),\hfill \\ & & \\ \hfill Mh^2_\tau ^2v2Mc_\chi _\tau v& =& \frac{Mc^2}{12}_\chi ^4v+Z^{}(L)_\chi ^2u_\chi v+Z^{}(L)_\chi u_\chi ^2v+O(h).\hfill \end{array}$$
(96)
As $`h0`$ with $`\mu h^2`$, we find the coupled system

$$\begin{array}{ccccc}\hfill _\tau A& +& _\chi \left[\frac{1}{2}A^2+\frac{c}{24}_\chi ^2A\right]& =& 0,\hfill \\ & & & & \\ \hfill _\tau B& +& _\chi \left[\frac{Z^{}(L)}{W^{\prime \prime \prime }(L)}AB+\frac{c}{24}_\chi ^2B\right]& =& 0,\hfill \end{array}$$
(97)
as a formal limit, where $`A=W^{\prime \prime \prime }(L)_\chi u/(2Mc)`$ and $`B=_\chi v`$. After a simple rescaling of $`\chi `$ and $`\tau `$, this becomes (1) and (2) with $`\kappa =Z^{}(L)/W^{\prime \prime \prime }(L)`$. As described in § 6, the influence of solitons on linear waves can be qualitatively different for different values of the coupling constant $`\kappa `$, with important bifurcations occurring at the values $`\kappa =\kappa _n^{\mathrm{bif}}=(n+1)(n+2)/12`$, for $`n=1,2,3,\mathrm{}`$. As we have seen, this coupled system can be solved exactly in (at least) two cases: $`\kappa =1`$, and $`\kappa =1/2`$. The former case is just the linearized KdV; see Sachs \[S83\]. The latter case is the one that is solved in the main text of this paper.

Consider this example with the potential of the strong and weak springs given respectively by

$$W(D):=\frac{1}{2}\kappa _wD^2+\frac{1}{24}\alpha D^4,\text{ and }V(D):=\mu (\frac{1}{2}\kappa _vD^2+\frac{1}{24}\beta D^4).$$
(98)
Thus $`\beta =3\alpha \kappa `$ and the condition that the wave speeds are equal is

$$\frac{1}{2}\alpha L^2(\kappa 1)=\kappa _w\kappa _v.$$
(99)
The effect of each bifurcation point in $`\kappa `$ is now clear. At $`\kappa =1/2`$ the first harmonic of the $`v`$-waves begins to resonate with the $`u`$-waves. As $`\kappa `$ increases through $`\kappa =1`$ we pass through a transition from supercritical resonance to subcritical resonance. Similarly, at the odd bifurcation points, $`\kappa =\kappa _{2m1}^{\mathrm{bif}}`$, for $`m=1,2,3,\mathrm{}`$, the $`m^{\mathrm{th}}`$ harmonic $`v`$-wave begins to resonate with the $`u`$-waves. Then at the even bifurcation points, $`\kappa =\kappa _{2m}^{\mathrm{bif}}`$, for $`m=1,2,3,\mathrm{}`$, the nature of the resonance for this mode changes from supercritical to subcritical.

Coupled systems of equations like the pair (1) and (2) often arise as formal asymptotic reductions of mechanical models for complicated one-dimensional waves. Often these asymptotic models are integrable. For example, in an elastic rod, the interaction between axial twist waves and helical deformation waves gives rise to an integrable Manakov system of coupled nonlinear Schrödinger equations \[LG99\].

The coupled system (1) and (2) for $`\kappa =1/2`$ is also intimately connected with an integrable multicomponent (an arbitrary number of components, all appearing symmetrically) coupled KdV equation \[MC99\]. Indeed, from one point of view it is this connection that yields the solvability of (1) and (2) for $`\kappa =1/2`$ described in detail in this paper. The solution method presented here also can be used to give the complete solution of the coupled KdV system. That system in turn can be interpreted as a phenomenological model for the transport of the mass integral through an $`N`$-soliton solution of KdV \[MC99\].

Finally, we would like to point out that there are also some applications in which linear equations of the form (2) occur with $`\kappa A(x,t)`$ being a given function. In this case, $`c(x,t)=\kappa A(x,t)`$ represents a given spatiotemporal modulation of the speed of linear dispersive waves, say due to propagation in an inhomogeneous medium. Such problems arise in the modeling of the propagation of weak internal waves in a channel of varying width \[CG99\]. For such applications, we may view the solvability of the coupled system (1) and (2) for $`\kappa =1`$ and $`\kappa =1/2`$ as a kind of (big) catalog of special cases of the function $`c(x,t)`$ for which the linear equation (2) is solvable in its own right. For other values of $`\kappa >1/2`$ the linear wave system is unstable, while for all values of $`\kappa <1/2`$ it is stable.


format: text

text:
# Vortices in Bose-Einstein condensates confined in a multiply connected Laguerre-Gaussian optical trap.

## I Introduction

Superfluidity, and its characteristic manifestation as a state with quantized circulation (a vortex), is intimately connected with the phenomenon of Bose-Einstein condensation, yet the precise relation between superfluid persistent currents and Bose-Einstein condensation (BEC) is only beginning to be elucidated . Before 1995, the link between superfluidity and BEC was almost exclusively studied in the context of liquid helium (4He and 3He) , where the study of the relation between superfluidity and BEC is complicated by the strong interatomic interactions between the atoms in the liquid. In 1995, Bose-Einstein condensation was realized in magnetically trapped clouds of alkali atoms . In these novel condensates, the bosonic atoms are weakly interacting (contrary to the case of liquid helium), and as such these systems have the potential to shed new light on superfluidity. Soon after the initial creation of alkali gas condensates, several experimental groups set out to create a vortex – a quantum of superfluid circulation – in this novel system .

The initial attempts to create a vortex in a condensate by stirring the trapped condensate with a blue detuned laser beam  were unsuccessful, even though early theoretical work by Dalfovo et al. indicated that persistent superfluid currents can indeed manifest themselves in Bose-Einstein condensates as vortices . Subsequent analysis  showed that vortices are unstable in the simply-connected, not-stirred, spin-polarized condensates formed in the original experiments .

This can be understood as follows. Along the vortex line, the order parameter of the Bose-Einstein condensate has to vanish. Phrased metaphorically, a “hole” has to be “drilled” in the condensate along the vortex line. In the magnetic trap, the modulus square of the order parameter of a condensate without a vortex is largest in the center of the trap. As a consequence, a vortex line through the center of the trap will perturb the order parameter more than a vortex line at the edge of the condensate (since, using our metaphor again, it will cost more energy to “drill” through the center of the condensate than through its edge). It will be energetically favorable for the vortex line to be at the edge of the condensate. In the presence of dissipation, this will cause the vortex line to migrate to the edge of the condensate so that the vortex condensate will decay into a non-vortex state. This argument, sketched here with some roughness, has been worked out with precision by Fetter and co-workers in .

Several schemes have been proposed, both theoretically and experimentally, to stabilize vortices: rotating the trapping potential (analogous to rotating a bucket containing 4He)  or stirring the condensate with an off-resonance laser , raising the temperature (to ‘pin’ the vortex in the potential created by the non-condensate fraction at the center of the vortex) , phase imprinting methods , and various other techniques . Recently, vortices were created experimentally, both with a ‘rotating bucket’ experiment  and with the use of a two-component condensate . In the latter experiment, one of the components of the spinor condensate ‘pins’ the vortex present in the other component and a Ramsey type interference between the two components is used to detect the vortex. The long lifetimes of the two-component condensates in ref.  is due to similar singlet and triplet scattering lengths of rubidium, resulting in an anomalously low inelastic loss rate .

Nevertheless, stable vortices have not yet been realized in a spin-polarized condensate in non-rotating traps, and new methods of stabilization and detection must be developed. A promising scheme, based on an analysis similar to that of Fetter , is the use of multiply-connected condensates. A candidate trap to create a multiply connected condensates consists of a red-detuned laser beam in a Laguerre-Gauss mode , which we discuss in section II. If the condensate order parameter has a toroidal geometry, a vortex line threaded through the cylindrical symmetry axis of the torus will not perturb the condensate order parameter strongly, and moreover create a metastability barrier for vortex decay. One of the goals of this paper is to verify this statement about vortex stability (in section III). In a condensate with a toroidal order parameter, the vortex line can no longer be detected as a line along which the density of Bose-Einstein condensed atoms vanishes. To detect vorticity, we propose (in section IV) a method based on interference, similar to that proposed in refs.  for simply connected condensates.

## II Bose-Einstein condensation in Laguerre Gauss traps

### A Trapping geometry

Toroidal confinement for ultracold atoms can be obtained by an optical dipole trap  which consists of a laser beam in a Laguerre-Gaussian mode. The Laguerre-Gaussian mode $`\{n,m\}`$ is characterized by an intensity profile given by :

$$I_{n,l}(r,z)\frac{(2r^2/W_0^2)^l}{1+(2z/kW_0^2)^2}L_{n,l}^2\left(2r^2/W_0^2\right)\mathrm{exp}\left\{\frac{2r^2}{W_0^2}\right\},$$
(1)
where $`r`$ is the radial distance from the center of the beam, $`z`$ is the position along the propagation direction of the beam, $`W_0`$ (referred to as the ‘waist parameter’) is a parameter controlling the minimal width of the beam, $`k`$ is the wave number of the laser, and $`L_{n,l}`$ is the Laguerre polynomial of order $`\{n,l\}`$. Such laser beams have modes which show a node in the center, and trap the atoms in a cylindrical shell around the axis of propagation of the beam. The intensity profile of the laser beam in the Laguerre-Gaussian mode is illustrated in Figure 1. Laguerre-Gaussian (LG) laser beams have already been used to successfully trap atoms , and subsequent theoretical work has shown that toroidal traps formed by a red-detuned LG beam can be loaded from initial conditions similar to those of conventional magnetic traps . Once the trap has been loaded, one possibility to create a vortex in the trapped gas would be by a phase imprinting method, which already successfully resulted in the creation of solitons in condensates . Alternatively, a rotating perturbing potential can be used to stir the condensate and set up persistent flow in the toroidal geometry .

In this paper we investigate the properties of the vortex Bose-Einstein condensate, optically trapped by a laser beam in a Laguerre-Gaussian propagation mode $`\{n,l\}=\{0,1\}`$. Along the $`z`$-axis (the direction of propagation of the laser beam) an additional magnetic trap  results in a parabolic $`z`$-axis confinement with frequency $`\mathrm{\Omega }`$, which prevents the atoms from escaping along the direction of propagation of the laser beam. This ‘plugging’ of the optical trap was achieved in  by using blue-detuned ‘plugging beams’. The condensate in the Laguerre-Gauss geometry  will be denoted by ‘Laguerre-Gaussian condensate’ (LG condensate), in contrast with the condensate in a parabolic confinement.

The intensity profile of such a Laguerre-Gauss beam in the $`zy`$-plane (where $`z`$ is the axis of propagation of the laser beam), given by (1), is shown in Figure 1. In the remainder of this paper, we use units so that $`\mathrm{}=m=\mathrm{\Omega }=1`$ (where $`m`$ is the mass of the atoms). In these units, the laser beam parameters for Figure 1 are chosen as follows: $`k=2,W_0=5`$. Figure 2 illustrates a surface of constant intensity of the laser beam; a hollow cylindrical shell. This will also be the shape of the cloud of trapped atoms . The extension of the cloud along the $`z`$-axis can be tuned by selecting the $`\mathrm{\Omega }`$ frequency of the magnetic trap along the $`z`$-axis.

### B Mean-field approach

Confined Bose-Einstein condensates are well described by a mean-field theory where the properties of the condensate are derived from a complex function $`\mathrm{\Psi }`$. This function (the order parameter) is interpreted as a macroscopic wave function and obeys the Gross-Pitaevskii equation :

$$\frac{\mathrm{}^2}{2m}\mathrm{\Delta }\mathrm{\Psi }+V_{conf}(\rho ,z)\mathrm{\Psi }+U_0|\mathrm{\Psi }|^2\mathrm{\Psi }=E\mathrm{\Psi }.$$
(2)
The term nonlinear in $`\mathrm{\Psi }`$ in (2) arises from the interparticle interaction potential, which is treated as a contact potential with scattering length $`a_{scat}`$ so that $`U_0=4\pi \mathrm{}^2a_{scat}/m`$ where $`m`$ is the mass of an atom. The optical confinement is generated by a red-detuned Laguerre-Gaussian laser beam in the $`\{0,1\}`$ mode, with waist parameter $`W_0`$ and wave number $`k`$, and propagating along the $`z`$-axis . In addition to this, a harmonic confinement (characterized by a frequency $`\mathrm{\Omega }`$) is present which confines the atoms along in the $`z`$-direction to a region $`z<kW_0^2`$ . The potential energy corresponding to this geometry is given, in cylindrical coordinates $`\{\rho ,\theta ,z\}`$, by:

$$V_{conf}(\rho ,\theta ,z)=A_0\frac{2\rho ^2}{W_0^2}\mathrm{exp}\left\{2\rho ^2/W_0^2\right\}+\frac{m\mathrm{\Omega }^2}{2}z^2.$$
(3)
In the present treatment, the Gross-Pitaevskii equation (2) will be solved variationally . Consider for this purpose the following two trial functions, expressed in cylindrical coordinates $`\{\rho ,\theta ,z\}`$:

$`\psi (\rho ,\theta ,z)`$ $`=`$ $`𝒩_0\text{ }\mathrm{exp}\{\sigma _0(\rho \rho _0)^2/2\zeta _0z^2/2\},`$ (4)
$`\psi _\text{v}(\rho ,\theta ,z)`$ $`=`$ $`𝒩_1\text{ }\rho \mathrm{exp}\{\sigma _1(\rho \rho _1)^2/2\zeta _1z^2/2\}e^{i\theta }.`$ (5)
The variational parameters in the trial functions are $`\sigma _0,\sigma _1`$, controlling the width of the functions in the radial direction; $`\zeta _0`$,$`\zeta _1`$, controlling the widths in the $`z`$-direction; and $`\rho _0,\rho _1,`$ controlling the radial displacement of the maximum of the function away from the $`z`$-axis. $`𝒩_0`$ and $`𝒩_1`$ are normalization constants, which are determined by demanding that the norm of the trial function equals the number of particles in the condensate.

The circulation around a closed loop in the trapped, interacting Bose gas is defined by $`\kappa =\frac{\mathrm{}}{m}S,`$ where $`S`$ represents the phase of the order parameter solving the Gross-Pitaevskii equation (2). The single-valuedness of the order parameter ensures that the circulation in the trapped Bose gas is quantized: $`\kappa =nh/m`$ with $`n=0,1,2,\mathrm{}`$. A vortex is present whenever $`n>0`$. The trial function $`\psi (\rho ,\theta ,z)`$ has a constant phase and hence there is no circulation, no vortex. On the other hand, the phase of $`\psi _\text{v}(\rho ,\theta ,z)`$ changes by $`2\pi `$ along any closed loop encircling the $`z`$-axis. Thus $`\psi _\text{v}`$ is a trial function for the order parameter of a condensate with one quantum of vorticity: its circulation equals $`\kappa =h/m`$. The trial function $`\psi _\text{v}`$ must have a node along the $`z`$-axis (the factor $`\rho `$ in expression (5)): if this node were not present, the variational kinetic energy would diverge along the $`z`$-axis.

We have used $`\psi `$ and $`\psi _\text{v}`$ as variational trial functions to solve the Gross-Pitaevskii equation for a condensate without and with a vortex, respectively. The resulting variational energy is shown in Figure 3 as a function of $`Na_{scat}/a_{HO}`$ where $`N`$ is the number of atoms and $`a_{HO}=\sqrt{\mathrm{}/(m\mathrm{\Omega })}`$. The quantity $`Na_{scat}/a_{HO}`$ is a dimensionless measure of the strength of the interaction: upon increasing $`Na_{scat}/a_{HO}>0`$ the interatomic interaction becomes more repulsive. The energy of the LG condensate with a vortex is larger than the energy of the LG condensate without a vortex for all interaction strengths investigated. The difference is smallest for the non-interacting gas and increases monotonically as $`Na_{scat}/a_{HO}`$ increases. This means that for all investigated interactions strengths the vortex LG condensate is either metastable or unstable. In inset (a) of Figure 3, the optimal value of the variational parameters $`\sigma ,\zeta `$ are shown for both trial functions as a function of the interaction strength. Both the variational trial function of the LG condensate with vortex and without vortex are broadened under the influence of the repulsive interactions.

## III Vortex metastability barriers

A Laguerre-Gaussian condensate with a vortex can be metastable if there exists an energy barrier separating this state from a state without a vortex. In this section we estimate the height of the energy barrier both by the method proposed by Benakli et al.  for 2D traps, and by the method proposed by Fetter and co-workers .

### A Hydrodynamic and microscopic instability

The hydrodynamical and the microscopic instabilities studied by Fetter and co-workers  involve the displacement of the vortex core relative to the center of the trap. As discussed in the introduction, for a parabolically trapped condensate it is energetically favorable for the distance between the vortex core and the center of the trap to increase. Thus, if dissipation is present, the vortex condensate in the parabolic trap can decay into a non-vortex condensate through a migration of the vortex to the edge of the cloud. This effect is difficult to observe experimentally, since the vortex core can tilt and the image contrast between the vortex core and the cloud of atoms is reduced .

To estimate the height of the energy barrier separating the vortex condensate from the condensate without a vortex in our present case, we introduce a third variational function to represent the order parameter of a vortex LG-condensate where the vortex core is at a given distance $`R`$ from the cylindrical symmetry axis of the Laguerre-Gauss trapping beam:

$$\psi _\text{v}(R;\rho ,\theta ,z)=𝒩_2\text{ }r(\rho ,\theta )\mathrm{exp}\{\sigma _R(\rho \rho _R)^2/2\zeta _Rz^2/2\}e^{i\mathrm{\Theta }(\rho ,\theta )}.$$
(6)
The function $`r(\rho ,\theta )`$ gives the distance between the point $`\{\rho ,\theta ,0\}`$ and the vortex core:

$$r(\rho ,\theta )=\sqrt{(\rho \mathrm{cos}\theta R)^2+\rho ^2\mathrm{sin}^2\theta },$$
(7)
and the function $`\mathrm{\Theta }(\rho ,\theta )`$ gives the angle between the line connecting the vortex core with the point $`\{\rho ,\theta ,0\}`$ and the $`x`$-axis:

$$\mathrm{\Theta }(\rho ,\theta )=\mathrm{arctan}\left(\frac{\rho \mathrm{sin}\theta }{\rho \mathrm{cos}\theta R}\right).$$
(8)
The parameters $`\{\sigma _R,\rho _R,\zeta _R\}`$ are determined variationally for every studied distance $`R`$ between the vortex core and the center of the trapping geometry (the axis of propagation of the LG beam), and $`𝒩_2`$ is a normalization constant. The limit $`R0`$ retrieves the results for the vortex condensate studied earlier. In the limit $`R\mathrm{}`$, the result tends to the result for the condensate without a vortex. Using the variational approach, we calculate the energy $`E(R)`$ of the displaced vortex state for any intermediate $`R`$. If the energy $`E(R)`$ decreases monotonously with increasing $`R`$, the vortex is unstable - in the presence of dissipation the vortex condensate will decay. If there exists a maximum energy $`E(R_{\text{max}})`$ for an $`R_{\text{max}}`$ different from zero, there is an energy barrier $`E_{\text{barrier}}=E(R_{\text{max}})E(R=0)`$ which will hinder the vortex decay.

In Figure 4, the energy barrier $`E_{\text{barrier}}`$ is shown as a function of $`Na_{\text{scat}}/a_{\text{HO}},`$ the dimensionless measure of the interaction strength used in the previous section. The calculations were performed for a configuration so that $`A_0=5`$ $`\mathrm{}\mathrm{\Omega },`$ $`W_0^2=20`$ $`a_{HO}^2`$. We found that a metastability barrier exists ($`E_{\text{barrier}}>0`$) for $`Na_{scat}/a_{HO}>1.6\pm 0.1.`$ This means that vortices are metastable only if the effective interaction strength (controlled by the scattering length, the number of atoms and the confinement strength) is large enough: the metastability is induced by interactions. Upon further increasing the interaction strength above the threshold value, the metastability barrier increases. The variational energy $`E(R)`$ is shown as a function of $`R`$ in the insets: once for an interaction strength such that the vortex is not stable, and once for a situation where the vortex is metastable. A reasonable estimate of the experimental parameters required to realize this trapping geometry gives $`a_{HO}=4`$ $`\mu `$m and $`\mathrm{\Omega }=50`$ Hz. This implies that for 87Rb, the critical number of trapped atoms necessary to make vortices metastable is of the order of $`10^3`$. Since one can reasonably expect to trap more than $`10^5`$ atoms in the trap, vortices created in a Bose gas in the confinement potential mentioned above, will be well into the metastable regime ($`Na_{scat}/a_{HO}625`$ for $`N=500000`$).

### B Uniform transition to a non-vortex state

The other method which we used to study the metastability of the vortices in Laguerre-Gauss condensates, was proposed by Benakli et al.  for 2D traps with a an axial hole, punctured by an off-resonance laser beam. In this method, a trial solution for the Gross-Pitaevskii equation is constructed as a superposition between the vortex state and the state without a vortex : $`\mathrm{\Psi }(\rho ,\theta ,z)=C_0\psi (\rho ,\theta ,z)+C_1\psi _\text{v}(\rho ,\theta ,z)`$ where $`|C_0|^2+|C_1|^2=1.`$ Adapting the coefficients $`C_0`$ and $`C_1`$ of this superposition, the function changes from that of a LG condensate with a vortex to that without a vortex. If the energy of the intermediate states is found to be higher than the energy of the vortex state, this constitutes an energy barrier against uniform transition from vortex state to the state without a vortex, making the former metastable.

Using the variationally optimized trial functions for the order parameter of a LG condensate with a vortex and without a vortex, we found with the method of Benakli et al. that such a metastability barrier exists for $`Na_{scat}/a_{HO}>1.9\pm 0.1.`$ The dependence of the energy barrier on $`Na_{scat}/a_{HO}`$ is shown in Figure 5. In the inset of Figure 5, the variational energy of the trial function $`\mathrm{\Psi }`$ is shown as a function of $`|C_0|^2`$ for an interaction strength below (a) and above (b) the threshold value for metastable vortices. The calculations were performed for a configuration so that $`A_0=5`$ $`\mathrm{}\mathrm{\Omega },`$ $`W_0^2=20`$ $`a_{HO}`$. Note that the method of ref.  followed here does not describe the decay of vortices mediated by the excitations, for example by nucleation of vortex loops. However, for a two-dimensional system punctured by an off resonance laser beam, Benakli et al.  could show that excitation-mediated decay of vortices only becomes appreciable for vortices in higher angular momentum states. Furthermore, for a purely one-dimensional toroidal system, Kagan et al.  showed that for a weakly interacting system at low temperature, the relaxation time for the decay the superfluid persistent current due to phonons is strongly suppressed.

Both the method of Fetter  and the method of Benakli et al.  give qualitatively the same result, namely that above a critical strength of interaction, expressed by $`Na_{scat}/a_{HO}`$, the vortex state of the LG condensate is stabilized by a metastability barrier. Quantitative differences remain between the two methods, which is to be expected since they describe different possible mechanisms of vortex decay. The general conclusion – the stabilization of the vortex state due to the toroidal geometry – is further supported by a calculation for a vortex condensate in a mexican hat potential .The metastability barrier calculated as by the method of ref.  is smaller than the barrier calculated used the method of ref.  for $`Na_{scat}/a_{HO}3.5`$, which indicates that the latter method  overestimates the real metastability barrier.

## IV Interference and detection of vorticity

Phase coherence, a property necessary for the existence of vortices, was demonstrated experimentally using interference experiments . Two parabolically trapped condensates, displaced by a given distance, are allowed to expand freely. In the region where the two expanding condensates overlap, an interference pattern is observed. This type of experiment has also been proposed to observe vorticity in parabolically trapped condensates : the presence of a vortex leads to an observable edge dislocation in the pattern of otherwise parallel interference fringes.

Consider a parabolically trapped condensate positioned on the symmetry axis of the LG condensate. The LG condensate forms a cylindrical shell surrounding the parabolically trapped condensate in the center. As the trapping potentials are switched off, both condensates will expand and exhibit an interference pattern in the region of overlap. This interference pattern is the subject of the present section.

To find the function representing the order parameter at a given time $`t`$ after the start of the free expansion, the original function $`\psi _\text{v}(\rho ,\theta ,z;t=0)`$ is expanded in free particle eigenfunctions. These eigenfunctions acquire a phase factor as time elapses, so that $`\psi _\text{v}`$ at a time $`t`$ is given by the resummed expansion with the ‘time evolved’ eigenfunctions. More explicitly, for the LG condensate with a vortex, first the coefficients of the plane wave expansion are evaluated:

$$\psi _\text{v}(\rho ,\theta ,z;t=0)=c(𝐤)\frac{e^{i𝐤.𝐫}}{(2\pi )^{3/2}}𝑑𝐤,$$
(9)
$`c(k_\rho ,\varphi ,k_z)`$ $`=`$ $`{\displaystyle \psi _\text{v}(\rho ,\theta ,z)\frac{\mathrm{exp}\{ik_\rho \rho \mathrm{cos}(\theta \varphi )ik_zz\}}{(2\pi )^{3/2}}𝑑𝐤}`$ (10)
$`=`$ $`{\displaystyle \frac{\mathrm{exp}\left\{k_z^2/2\zeta \right\}}{\sqrt{2\pi \zeta }}}\times e^{i(\varphi +\pi /2)}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑\rho \text{ }\rho ^2J_1(k_\rho \rho )\mathrm{exp}\{\sigma (\rho \rho _0)^2/2\}.`$ (11)
In this expression $`J_1(x)`$ is the Bessel function of first order of the first kind and the wave number $`𝐤`$ is expressed in cylindrical coordinates $`\{k_\rho ,\varphi ,k_z\}`$. As time elapses, the free particles eigenfunctions acquire a phase factor $`\mathrm{exp}\{i\mathrm{}k^2t/(2m)\}`$. The function at time $`t`$ after the start of the free expansion is found by resumming the eigenfunctions at time $`t`$ :

$`\psi _\text{v}(\rho ,\theta ,z;t)`$ $`=`$ $`{\displaystyle c(𝐤)\frac{e^{i𝐤.𝐫+i\mathrm{}k^2t/(2m)}}{(2\pi )^{3/2}}𝑑𝐤}`$ (12)
$`=`$ $`\sqrt{{\displaystyle \frac{1}{1+i\zeta t}}}\mathrm{exp}\left\{{\displaystyle \frac{\zeta z}{2(1+i\zeta t)}}\right\}\times e^{i\theta }`$ (14)
$`\times {\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dk_\rho \text{ }(\rho ^{})^2{\displaystyle \frac{J_1(\rho ^{}\rho /t)}{t}}e^{\sigma (\rho ^{}\rho _0)^2/2}\mathrm{exp}\left\{i{\displaystyle \frac{\rho ^2+(\rho ^{})^2}{2t}}\right\}.`$
The time evolution for the free expansion of a condensate prepared in a parabolic confinement is derived analogously . The total measured density generated by the two condensates is then given by $`|\psi _\text{v}(\rho ,\theta ,z;t)+\psi (\rho ,\theta ,z;t)|^2`$ \- remember that each condensate function is normalized to the number of particles in the given condensate. Several time frames of the resulting evolution of the density of the expanding condensates are shown in Figure 3. The frames in Figure 3 show a cross-section of the density along the xy-plane.

Immediately after the traps are switched off, the density is that of a cylindrical, Laguerre-Gaussian condensate with a parabolically trapped condensate in the middle. As time goes by, both condensates expand: the parabolically trapped condensate expands radially, and the cylindrical shell of the LG condensate broadens. As the expanding condensates start to overlap the fringe pattern appears. If the LG condensate does not contain a quantum of superfluid circulation (i.e. no vortex), the interference pattern consists of a series of concentric circles with linearly increasing radius. If however the LG condensate does contain a vortex, the interference pattern is an Archimedean spiral. At higher vorticity, the number of arms in the Archimedean spiral equals the number of vortex quanta in the LG condensate.

Figure 3 shows that, as a function of time, the spiral interference pattern rotates around the cylindrical symmetry axis of the trapped condensates, with a frequency of the order of the frequencies characterizing the parabolic approximation to the trapping potentials. Furthermore, as can be seen from Figure 3, the distance between two successive windings of the spiral increases as time increases. Hence, to detect the spiral interference pattern, the density has to be measured on a time scale shorter than that given by the inverse of the characteristic trapping frequencies.

## V Conclusions

In conclusion, we have studied the properties of vortices in a condensate in an optical trap generated by a laser beam in a Laguerre-Gaussian mode and proposed a method to detect these vortices. The energy of the vortex state was calculated variationally in a mean-field framework, and a threshold interaction strength was found beyond which there exists a metastability barrier stabilizing the LG vortex state against a transition to a non-vortex state. Both the hydrodynamic instability  and the uniform transition to a non-vortex state  were considered in the investigation of the vortex stability. Below the critical interaction strength the metastability barrier vanishes, irrespective of the fact that the present confined geometry allows for the vortex core to be positioned in a region where the condensate is nonzero.

The time evolution under free expansion of a LG condensate was derived, and subsequently used to calculate the interference pattern which arises from the overlap of an expanding LG condensate with an expanding parabolically trapped condensate in its center. When superfluid circulation is present, the interference pattern consists of an Archimedean spiral, clearly distinct from the series of concentric cylindrical fringes which arise in the interference pattern when no vorticity is present. The observation of such a spiral interference pattern would constitute clear evidence for vorticity in a spin-polarized Bose-Einstein condensates.

## Acknowledgments

We thank F. Brosens and L. F. Lemmens for intensive discussions. Also discussions with T. Kuga are gratefully acknowledged. This work is performed within the framework of the FWO projects No. 1.5.545.98, G.0287.95, 9.0193.97, G.0071.98 and WO.073.94N (Wetenschappelijke Onderzoeksgemeenschap, Scientific Research Community of the FWO on “Low Dimensional Systems”), the “Interuniversitaire Attractiepolen – Belgische Staat, Diensten van de Eerste Minister – Wetenschappelijke, Technische en Culturele Aangelegenheden”, and in the framework of the GOA BOF UA 2000 projects of the Universiteit Antwerpen. One of the authors (J.T., aspirant bij het Fonds voor Wetenschappelijk Onderzoek - Vlaanderen) acknowledges the FWO for financial support.

## Figure captions

Figure 1 : The Laguerre-Gauss (LG) trap is an optical dipole trap consisting of a red-detuned laser beam. A typical intensity profile of the laser beam in the LG-mode $`\{0,1\}`$ is shown in this figure, in a cross section through the axis of propagation of the beam (the $`z`$-axis). The atoms will feel the optical dipole force attracting them to the (toroidal) region of highest intensity of the laser beam.$`\sqrt{20}=\mathrm{4.\hspace{0.17em}472\hspace{0.17em}14}`$

Figure 2 :A typical surface of constant intensity of the laser beam in the Laguerre-Gauss (LG) mode $`\{0,1\}`$ is shown. The atoms in the red-detuned LG beam will be attracted to a toroidal or cylindrical region such as that within the surface of Figure 2. An additional parabolic magnetic confinement potential can be added in the $`z`$-direction .

Figure 3 : The variational result for the energy of the condensate in the Laguerre-Gaussian optical trap is shown as a function of the interaction strength. The full curve shows the energy of the condensate without a vortex, the dashed curve shows the energy with a vortex. The inset depicts the results for the variation parameters $`\sigma _0,\sigma _1`$ in the trial function for the condensate with vortex (dashed line) and without (full line), as a function of interaction strength. The trapping (beam) parameters were chosen as follows: $`A_0=5`$ $`\mathrm{}\mathrm{\Omega },`$ $`W_0^2=20`$ $`a_{HO}^2`$. In all graphs of this Figure, energies are expressed in units $`\mathrm{}\mathrm{\Omega }`$ and lengths in units $`a_{HO}=\sqrt{\mathrm{}/(m\mathrm{\Omega })}`$.

Figure 4 : The energy barrier for removing a vortex from a condensate in a Laguerre-Gaussian optical trap, is depicted as a function of interaction strength. This metastability barrier $`E_{\text{barrier}}`$ was calculated using the formalism of Fetter and co-workers . In this formalism, $`E_{\text{barrier}}`$ is found by deriving the energy $`E(R)`$ as a function of the distance $`R`$ between the vortex core and the center of the trapping potential, here the axis of propagation of the laser beam. This is illustrated in the insets. For interaction strengths lower than a critical value $`Na_{scat}/a_{HO}<1.6\pm 0.1`$ the barrier vanishes and vortices in Laguerre-Gaussian condensates are unstable with respect to the non-vortex state. For interaction strengths above this critical value, a metastability barrier exists.

Figure 5 : The energy barrier for a uniform transition from a vortex condensate to a condensate without a vortex, in a Laguerre-Gaussian optical trap, is depicted as a function of interaction strength. The energy barrier for this mechanism was calculated by the method of Benakli et al. . In the inset, the energy per particle of the variational function $`\mathrm{\Psi }=C_0\psi +C_1\psi _v`$ is given as a function of $`|C_0|^2`$, for two different interaction strengths: one below (a) and one above (b) the threshold for metastability of vortices, and the relation to the energy barrier for a uniform transition is shown.

Figure 6 : Several time frames in the evolution of freely expanding, overlapping condensates are shown. The gray scale represents the density in a cross-section orthogonal to the symmetry axis of the trap (the direction of propagation of the Laguerre-Gaussian laser beam), with black being the maximum density and white the minimum density. The initial situation is depicted in the top left panel: a parabolically trapped condensate in the center is surrounded by a cylindrical shaped Laguerre-Gaussian condensate containing a vortex. The time evolution of the density in the region within the rectangle in the top left panel is shown in more detail in the subsequent panels. When the confinement of both condensates is switched off, they expand and overlap - the density in the gap between both condensates increases - and a spiral interference pattern is formed. Units are chosen so that lengths are expressed in $`a_{HO}=\sqrt{\mathrm{}/(m\mathrm{\Omega })}`$ and the time is in units of $`1/\mathrm{\Omega }`$.