text:
# The GPE and higher order theories in one-dimensional Bose gases

## I Introduction

Recent advances in experimental techniques Görlitz et al. (2001); Schreck et al. (2001) have stimulated interest in the topic of ultra-cold dilute Bose gases in one dimension. The one-dimensional (1D) case leads to some interesting physics which does not occur in higher dimensions, and is also of theoretical importance due to the existence of exactly soluble models Girardeau (1960); Lieb and Liniger (1963); Lieb (1963). In this paper we investigate the 1D Bose gas in the limit of zero temperature. We extend our earlier work Lee et al. (2002) on the interactions in two dimensions to the 1D case, and obtain a 1D form of Gross-Pitaevskii equation (GPE). We also apply standard methods of mean-field theory developed to describe three-dimensional Bose-Einstein condensates (BECs) to the 1D trapped case. In both cases we compare the results obtained from these methods to the exact solutions.

It is well known that a true BEC cannot occur at any finite temperature in an interacting homogeneous Bose gas in 1D Hohenberg (1967), due to long wavelength fluctuations; nor can it occur in the limit of zero temperature Pitaevskii and Stringari (1991), due to quantum fluctuations. However, the presence of a trapping potential changes the density of states at low energies, and in the weakly-interacting limit a BEC may be formed Ketterle and van Druten (1996). In a recent paper, Petrov *et al.* discussed three different regimes of quantum degeneracy which can occur in a confined 1D interacting system Petrov et al. (2000): BEC, quasi-condensate, and Tonks-Girardeau gas. In the weakly interacting limit a BEC exists, but as the interactions become stronger the mean-field energy becomes important. When this energy becomes of the order of the energy level spacing of the trap, fluctuations again become significant. In this regime, the system forms a quasi-condensate, with local phase coherence, rather than the global coherence associated with a true BEC. In this regime the gas is similar to the homogeneous case, in which it was shown that at $`T=0`$ the phase coherence has only a power-law decay, rather than the exponential decay occurring at higher temperatures Schwartz (1977). Finally, in the limit of very strong interactions, the particles become impenetrable and the system is known as a Tonks-Girardeau gas Tonks (1936); Girardeau (1960). This limit is exactly soluble, due to a mapping relationship between the wave function for the bosons and that of a non-interacting fermion gas Girardeau (1960).

In the following section we briefly outline the relevant exact results for a 1D gas. Section III then presents a form of the GPE in which we use an approximate many-body T-matrix to describe the interactions in a 1D gas. The GPE gives the correct results in both the weakly-interacting and Tonks-Girardeau limits, and interpolates smoothly between these regimes. We confirm that the non-linear dependence in the GPE changes from cubic in the weakly-interacting regime to quintic in the strongly-interacting regime, and that this arises due to the many-body effects on the scattering of particles. In Sec. IV, we then discuss several standard mean-field theories and present the results of these for the trapped 1D gas. In the weakly interacting limit these theories should be valid since a true BEC exists in this regime. Our aim is to push several different theories into the strongly-interacting limit with the dual objectives of describing the 1D system and investigating the strengths and weaknesses of these theories.

## II Exact results in 1D

The 1D homogeneous system of $`N`$ bosons which interact via a delta-function potential of the form

$$V(x_1x_2)=V_0\delta (x_1x_2),$$
(1)
is exactly integrable. In the limit $`V_0\mathrm{}`$ the particles are impenetrable (the Tonks-Girardeau gas), and exact results can be obtained even in spatially inhomogeneous systems. The exact solutions are found by applying the appropriate boundary conditions to the wave function describing the system.

In the impenetrable limit the wave function must necessarily vanish wherever two particles meet (i.e. wherever $`x_i=x_j`$). Girardeau showed in Girardeau (1960) that the exact $`N`$-boson ground state wave function subject to such a boundary condition is equal to $`|\psi _F(x_1,\mathrm{},x_N)|`$, where $`\psi _F(x_1,\mathrm{},x_N)`$ is the wave function of $`N`$ non-interacting spinless fermions governed by the single-particle Hamiltonian $`\widehat{H}_{\mathrm{sp}}`$ (including terms from any external potential). Since the wave functions of such fermions must be antisymmetric with respect to co-ordinate exchanges, their wave functions automatically vanish whenever $`x_i=x_j`$, and so the Tonks-Girardeau gas boundary condition is implicitly obeyed by such a gas.

A number of bulk properties of a 1D Tonks-Girardeau gas may be calculated quite trivially as a result of this Fermi-Bose mapping theorem Girardeau (1960); Girardeau et al. (2001). For example, the ground state energy of a Tonks-Girardeau gas is simply found from a summation of the energies of the first $`N`$ non-interacting single-particle states (i.e. those that would be occupied by $`N`$ ideal fermions in the ground state configuration). Thus, for a homogeneous gas Girardeau (1960)

$$\frac{E_0}{N}=\frac{\mathrm{}^2}{m}\frac{\pi ^2}{6}n^2,$$
(2)
where $`n`$ is the density. The chemical potential is then easily calculated as

$$\mu =\frac{E_0}{N}=\frac{\mathrm{}^2}{m}\frac{\pi ^2}{2}n^2.$$
(3)
Similarly, the chemical potential for a Tonks-Girardeau gas in a harmonic potential of frequency $`\omega `$ is

$$\mu =\mathrm{}\omega (N+1/2).$$
(4)
The system of $`N`$ bosons interacting with an arbitrary value of $`V_0`$ was solved in 1963 by Lieb and Liniger Lieb and Liniger (1963); Lieb (1963) for a homogeneous gas. The Tonks-Girardeau gas results are recovered in the limit $`V_0\mathrm{}`$. The ground-state energy of the Lieb-Liniger gas is given (at $`T=0`$) by

$$\frac{E_0}{N}=\frac{\mathrm{}^2}{2m}n^2e(\gamma ),$$
(5)
where $`e(\gamma )`$ is a solution to the Lieb-Liniger system of equations, and is a function of the dimensionless quantity $`\gamma =V_0m/\mathrm{}^2n`$ which was found to be the important parameter in the system. The function $`e(\gamma )`$ has the following asymptotic behaviour

$$e(\gamma )\{\begin{array}{cc}\gamma \hfill & \text{ for }\gamma 0,\hfill \\ \frac{\pi ^2}{3}\left(\frac{\gamma }{\gamma +2}\right)\hfill & \text{ for }\gamma \mathrm{},\hfill \end{array}$$
(6)
and the energy per particle agrees with Eq. (2) in the impenetrable limit. The chemical potential is given by the equation

$$\mu =\frac{\mathrm{}^2}{2m}n^2\left[3e(\gamma )\gamma \frac{de(\gamma )}{d\gamma }\right].$$
(7)
In the $`\gamma 0`$ weakly-interacting limit, the exact results were found to be in agreement with the predictions of standard Bogoliubov theory. The Bogoliubov predictions begin to fail as $`\gamma `$ is increased and for $`\gamma 2`$ the results do not agree well with the exact theory Lieb and Liniger (1963). Interestingly, for fixed $`V_0`$ the $`\gamma 0`$ limit corresponds to the high-density limit, whilst a high value of $`\gamma `$ implies a small density. This is then the reverse of the 3D case where Bogoliubov theory is valid in the low density limit. In a subsequent section we will show to what extent higher-order modifications to Bogoliubov theory can improve the agreement for higher $`\gamma `$.

For a system where we keep the density constant, and increase the interaction strength $`V_0`$ from the non-interacting limit to the impenetrable limit (i.e. increasing $`\gamma `$), we expect the chemical potential to vary in the following manner. In the low-$`V_0`$ limit the chemical potential increases linearly with $`V_0`$, but as $`V_0`$ is increased further $`\mu `$ starts to saturate, and in the high-$`V_0`$ limit it becomes independent of $`V_0`$, in agreement with Eq. (3). In the following section we introduce a modified GPE which follows this behaviour.

## III The GPE in 1D

We consider the model system of a gas of one-dimensional point particles interacting with an arbitrary strength $`V_0`$, corresponding to the Lieb-Liniger gas. We emphasize here that we are interested in the limit in which the interparticle interactions in the system are genuinely one-dimensional, because it is only in this limit that we can compare to the exact results mentioned earlier.

We can write a one-dimensional Gross-Pitaevskii equation for a system of $`N`$ particles by a straightforward generalisation of the usual 3D equation, giving

$$\mu \psi (x)=\left[\frac{\mathrm{}^2}{2m}\frac{d^2}{dx^2}+V_{\mathrm{trap}}(x)+g_{1\mathrm{D}}N|\psi (x)|^2\right]\psi (x),$$
(8)
where $`\psi (x)`$ is the condensate wave function, and $`V_{\mathrm{trap}}(x)`$ represents the trapping potential. The coupling parameter $`g_{1\mathrm{D}}`$ describes the interactions between particles, and we discuss the nature of this quantity in the following section before presenting the results predicted by the GPE.

### III.1 Interactions and the GPE

The GPE can easily be derived by functional differentiation of the Hamiltonian for a gas with a general interparticle interaction $`V(x)`$. The coupling parameter is then given by the matrix element $`0|V(x)|0`$, where the particles in both the incoming and outgoing states have zero relative momentum (see, for example, Morgan (2000)). However, a representation in terms of matrix elements of the bare interparticle potential is undesirable, since for a realistic potential they are often very large and frequently the exact potential itself is unknown. In order to avoid this it is common for 3D problems to use the two-body T-matrix $`T_{2\mathrm{B}}`$ in place of the interparticle potential. This T-matrix describes the scattering of two particles in a vacuum, and it is possible to rewrite the interaction term in the Hamiltonian in terms of $`T_{2\mathrm{B}}`$ Morgan (2000). Carrying out the functional differentiation again leads to a GPE of the above form, where the coupling parameter is now $`g=0|T_{2\mathrm{B}}(E=0)|0`$. Here $`E`$ is the energy of the interacting particles, and for two condensate particles this is zero. In three dimensions this leads to the familiar result $`g_{3\mathrm{D}}=4\pi \mathrm{}^2a/m`$, where $`a`$ is the *s*-wave scattering length. However, in 1D (and 2D) the two-body T-matrix vanishes in the zero energy and momentum limit, and has significant imaginary terms at a finite real energy (see, for example, Morgan et al. (2002)).

In earlier papers Lee et al. (2002); Lee and Morgan (2002), we have argued that the many-body effects of the surrounding atoms on particle interactions, while being a small perturbation in 3D, are crucial in lower dimensions. The many-body effects can be incorporated into a many-body T-matrix $`T_{\mathrm{MB}}`$ Stoof et al. (1996); Bijlsma and Stoof (1997); Proukakis et al. (1998); Morgan (2000), and the coupling parameter in the GPE becomes

$$g_{1\mathrm{D}}=0|T_{\mathrm{MB}}(E=0)|0.$$
(9)
Unfortunately, the many-body T-matrix is very difficult to solve exactly. However, in Lee et al. (2002) we presented an approximate solution, valid in the zero temperature limit, in which the many-body T-matrix was approximated by the two-body T-matrix evaluated at a shifted effective interaction energy (off the energy shell). The GPE in 1D is therefore given by Eq. (8) with the coupling parameter

$$g_{1\mathrm{D}}=0|T_{\mathrm{MB}}(E=0)|00|T_{2\mathrm{B}}(E^{})|0,$$
(10)
where $`E^{}`$ is the shifted interaction energy, and was found in Lee et al. (2002) to be negative. As an illustration of the utility and accuracy of this approximation, we discuss its application to the 3D case in the appendix.

### III.2 The hard-core limit

We will look first at the GPE in the impenetrable, $`V_0\mathrm{}`$, limit in which the system is a Tonks-Girardeau gas. The off-shell two-body T-matrix for such a 1D gas, calculated in Morgan et al. (2002), is

$`0|T_{2\mathrm{B}}^{\mathrm{Tonks}}(E)|0=`$ (13)
$`\{\begin{array}{cc}\frac{2\mathrm{}^2}{ma}\left(i\sqrt{\frac{mE}{\mathrm{}^2}}a+\frac{mE}{\mathrm{}^2}a^2\right)\hfill & \text{for }E>0,\hfill \\ \frac{2\mathrm{}^2}{ma}\left(\sqrt{\frac{m|E|}{\mathrm{}^2}}a\frac{mE}{\mathrm{}^2}a^2\right)\hfill & \text{for }E<0,\hfill \end{array},`$
where $`a`$ is the size of the impenetrable particle, and for the delta-function gas we take the limit that $`a0`$. Combined with Eq. (10) this gives a coupling parameter of the form Lee et al. (2002)

$$g_{\mathrm{Tonks}}(\mu )=\sqrt{4C\mathrm{}^2\mu /m}$$
(14)
to leading order for a homogeneous system. Here we have taken $`E^{}=C\mu `$ and assumed that $`C>0`$ Lee et al. (2002), (note that if this is not the case then the coupling parameter is imaginary to leading order). Using the homogeneous solution to the GPE $`\mu =g_{\mathrm{Tonks}}n`$ we can rewrite the coupling parameter in terms of the density $`n`$ as

$$g_{\mathrm{Tonks}}(n)=\frac{4C\mathrm{}^2}{m}n.$$
(15)
Combining this with $`\mu =g_{\mathrm{Tonks}}n`$ gives $`\mu n^2`$, and a comparison with the exact result found by Girardeau and given in Eq. (3), shows that we have obtained the correct functional dependence of $`\mu `$ on $`n`$. This stems from the dependence of the coupling parameter on $`\sqrt{\mu }`$ in Eq. (14). This is further verification of our approximation for the many-body T-matrix in terms of the off-shell two-body T-matrix. A comparison with the Girardeau result also shows that $`C=\pi ^2/8`$.

To apply this result to a trapped system, we use a local density approximation, assuming that the scattering between two particles occurs on a length scale much shorter than the variation in the potential. First we solve the homogeneous expression $`\mu =g_{\mathrm{Tonks}}(\mu )n`$ to obtain an expression for the chemical potential in terms of the density $`n`$. Applying the local density approximation then requires $`\mu (n(x))=g_{\mathrm{Tonks}}(x)n(x)`$ to hold for all $`x`$, and we obtain a spatially-dependent coupling parameter

$$g_{\mathrm{Tonks}}(x)=\frac{\pi ^2\mathrm{}^2N}{2m}|\psi (x)|^2,$$
(16)
where we have used $`n(x)=N|\psi (x)|^2`$.

We can now write the GPE for a 1D trapped condensate as

$$\mu \psi (x)=\left[\frac{\mathrm{}^2}{2m}\frac{d^2}{dx^2}+V_{\mathrm{trap}}(x)+\frac{\pi ^2\mathrm{}^2N^2}{2m}|\psi (x)|^4\right]\psi (x).$$
(17)
This GPE differs significantly from the 3D version by virtue of the fact that the nonlinear term is now proportional to $`|\psi |^4\psi `$ instead of $`|\psi |^2\psi `$. In the Thomas-Fermi limit we can neglect the kinetic energy term and the solution for the wave function can easily be shown to be

$$|\psi (x)|=\frac{\left(\mu V_{\mathrm{trap}}(x)\right)^{1/4}}{(\mathrm{}^2\pi ^2N^2/2m)^{1/4}}\theta \left(\mu V_{\mathrm{trap}}(x)\right)$$
(18)
This differs substantially from the 3D and 2D cases, where the Thomas-Fermi wave function is instead proportional to $`(\mu V_{\mathrm{trap}}(x))^{1/2}`$, due to the different form of the non-linearity.

The Tonks-Girardeau GPE and the results presented in this section are in exact agreement with the earlier work of Kolomeisky *et al.* Kolomeisky et al. (2000), who obtained the same form of coupling parameter from a renormalisation group approach Kolomeisky and Straley (1992)111Note that in the work of Kolomeisky *et al.* , the authors predicted the correct functional form of the coupling parameter, but they also needed to correct their numerical factor $`C`$ in order to agree with the exact results Kolomeisky and Straley (1992), as we do..

In Fig. 1a we present a comparison of the density $`N|\psi |^2`$ obtained from the solution of the GPE of Eq. (17) with the exact result from the Bose-Fermi mapping theorem calculated by a summation of the first $`N`$ single-particle states of the trap. It can be seen that the GPE solution is in very good agreement with the exact result, and differs only in the absence of the oscillations seen in the exact result. These oscillations rapidly average out as the number of atoms is increased, and for a moderately large number the difference between the GPE results and the exact results will be negligible.

### III.3 The delta-function potential

We now turn to a discussion of the Lieb-Liniger gas. We are interested here in describing the gas from the non-interacting limit right through to the Tonks-Girardeau limit by means of an appropriate GPE.

The same argument which we used in deriving the Tonks-Girardeau GPE can again be used, namely that the coupling parameter should be the many-body T-matrix. Again we approximate this by a two-body T-matrix evaluated off the energy shell at an energy of $`C\mu `$, where we use the value of $`C=\pi ^2/8`$ obtained above.

The two-body T-matrix for a delta-function potential is given by Morgan et al. (2002)

$$0|T_{2\mathrm{B}}^{\mathrm{LL}}(E)|0=\frac{V_0}{1+iV_0\left({\displaystyle \frac{m}{2\mathrm{}^2}}\right)\sqrt{{\displaystyle \frac{\mathrm{}^2}{mE}}}}.$$
(19)
Note that for $`V_0\mathrm{}`$, this reduces to the Tonks limit of Eq. (13). It follows that our approximation for the many-body T-matrix in a homogeneous 1D system is

$$g_{\mathrm{LL}}(\mu )=\frac{V_0}{1+V_0\sqrt{{\displaystyle \frac{2m}{\mathrm{}^2}}{\displaystyle \frac{1}{\pi ^2\mu }}}}.$$
(20)
This energy-dependent form of the coupling parameter can be transformed to a density-dependent expression by means of the local density approximation, as outlined above. This leads to the expression

$`g_{\mathrm{LL}}(x)=`$
$`{\displaystyle \frac{1}{n(x)}}\left(\sqrt{{\displaystyle \frac{m}{2\pi ^2\mathrm{}^2}}V_0^2+V_0n(x)}\sqrt{{\displaystyle \frac{m}{2\pi ^2\mathrm{}^2}}}V_0\right)^2,`$
which describes an inhomogeneous system. Note that this expression gives the correct results for the homogeneous gas in the extreme limits $`V_00`$ and $`V_0\mathrm{}`$, giving $`\mu =nV_0`$ and Eq. (3) respectively.

Thus we can now use the above density-dependent expression for the coupling parameter in Eq. (8) to obtain a GPE for a trapped 1D gas for any arbitrary $`V_0`$.

We use the GPE here to describe the system over the full range of $`V_0`$, despite the fact that, as mentioned earlier, a true BEC only exists in a trap in the weakly-interacting limit Petrov et al. (2000). Once the interactions cause sufficient perturbation to the density of states that their discrete nature is smeared out then no condensate will form in 1D, and in the very high-$`V_0`$ (Tonks) limit it is known that there is no condensate. We can still use the GPE in the high-$`V_0`$ regime however, as our previous results have shown, provided that we now reinterpret the “wave function” as simply the square root of the density $`\psi (x)=\sqrt{n(x)}`$, rather than the coherent part of the field operator.

#### III.3.1 The Homogeneous Case

In the homogeneous case we can compare the results of the GPE with the interaction strength given by Eq. (III.3) to the exact results for a Lieb-Liniger gas. In Fig. 2 we compare the chemical potentials resulting from the GPE with the exact result of Eq. (3). The results from the GPE show the correct behaviour in the two extreme limits of $`\gamma `$, and the overall agreement is quite good. However, the chemical potential of the exact results can be seen to rise more steeply in the intermediate region. We note that the results have been obtained using a value of $`C`$ which strictly applies only in the large $`\gamma `$ limit. A larger value of $`C`$ gives a curve which agrees very well with the exact results in the weakly-interacting regime, but predicts the wrong value in the large $`\gamma `$ limit. It seems possible therefore that the value of $`C`$ should actually be dependent on $`\gamma `$, but a detailed investigation of this is outside the scope of this paper.

#### III.3.2 The Trapped Case

We now turn to the case of a trapped gas, where we will parameterise our results both by the interaction strength $`V_0`$ and by the dimensionless parameter $`\gamma _{\mathrm{trap}}=V_0m/\mathrm{}^2\stackrel{~}{n}_0`$ where $`\stackrel{~}{n}_0`$ is the peak density of the system (analogous to the parameter used by Lieb and Liniger).

We expect that in the low-$`V_0`$ (small $`\gamma _{\mathrm{trap}}`$) limit the chemical potential will vary as

$$\mu V_0^{2/3},$$
(22)
which is the Thomas-Fermi result obtained from the GPE with a cubic non-linearity (this is not quite the correct prediction in the extremely low-$`V_0`$ limit where the Thomas-Fermi approximation is invalid, but it should describe the majority of the low-$`V_0`$ regime adequately) . In the opposite high-$`V_0`$ limit, the Bose-Fermi mapping theorem for the trapped Tonks-Girardeau gas led to Eq. (4) which is completely independent of $`V_0`$. Note that our derivation of the GPE in Eq. (8) neglected terms of order $`1/N_0`$ since it was derived in the high-$`N_0`$ limit. For this reason we should expect only to reach an asymptotic chemical potential of $`\mu =N_0\mathrm{}\omega `$, which is approximately equal to Eq. (4) in this limit.

Figure 3 shows sample solutions to the GPE for a range of $`V_0`$ for $`N_0=15`$. In the weakly-interacting limit the wave function can clearly be seen to be a gaussian, as expected. However as the interactions are increased the wave function broadens until in the high-$`V_0`$ limit we have the shape discussed earlier in the hard-core case. Figure 4 shows the increase in the chemical potential over the same range. As expected from the exact results for the hard-core gas, in the high-$`V_0`$ limit the chemical potential becomes almost independent of the interaction strength. For the data shown the chemical potential does not appear to quite reach the limit $`\mu =N_0\mathrm{}\omega `$ expected from the exact results. However, in a calculation with the extreme value of $`V_0=10^4\mathrm{}\omega \mathrm{}_{\mathrm{trap}}`$ we obtained the result $`\mu =15.00\mathrm{}\omega `$ as expected. This agrees with the behaviour found by Lieb and Liniger in the homogeneous case, who found that the results reached the asymptotic limit very slowly (e.g. for $`\gamma 4`$ the chemical potential had only reached $`0.9`$ of its asymptotic value Lieb and Liniger (1963)).

We have therefore presented a form of the GPE in which the coupling parameter is an approximate many-body T-matrix given by the off-shell two-body T-matrix. We have shown that it agrees with the exact results for the Lieb-Liniger gas in the weakly-interacting and strongly-interacting limits, and it provides an interpolation between these limits in the intermediate regime. A drawback of using this GPE is that we lose any information about the coherences involved in the system, since we reinterpret the equation simply in terms of the density in the high-$`V_0`$ limit. The concept of using the many-body T-matrix in the GPE stemmed initially from higher order theories in which it appears naturally via the anomalous average (defined below). We therefore turn, in section IV, to an exploration of these theories in 1D and see how the changes in the interaction strength arise naturally.

#### III.3.3 Quasi-1D gases

Although we have been dealing with the case of a model delta-function potential gas, the results can to some extent be mapped onto those for a quasi-1D system. In practice, a 1D system could be created by squeezing a 3D trapped gas condensate in two dimensions. As the confinement is increased, the system progresses from being 3D to being 1D, via an intermediate regime we call quasi-1D for the purposes of this paper. In this intermediate regime, the dynamics of the system as a whole are frozen out in two dimensions, but the scattering between particles is not yet genuinely one-dimensional.

In the 3D regime the system is described by a GPE with a $`|\psi |^2\psi `$ nonlinearity, as for the low $`V_0`$ Lieb-Liniger gas. However, in the strict 1D limit the atoms will be unable to scatter past one another, and will form a Tonks-Girardeau gas. In this limit the appropriate GPE will be that of Eq. (17). The progression from a 3D gas to a 1D gas therefore mirrors in many respects the progression which occurs if the potential strength $`V_0`$ is increased for the Lieb-Liniger gas. The 3D to 1D progression has been investigated in more detail in references Olshanii (1998); Dunjko et al. (2001); Bergeman et al. (2002).

## IV Higher order theories in 1D

In the previous section we used the GPE in conjunction with an approximate many-body T-matrix in order to describe the Lieb-Liniger gas. In this section we examine the suitability of higher-order theories to describe this gas. Ideally we would like a theory which makes the transition between the weakly interacting and strongly-interacting regimes smoothly, and which accurately describes the coherence of the system as well as the densities. There are several higher order theories which are used to good effect to describe 3D condensates, and in this section we will apply some of these to the 1D trapped gas case.

Since we can always choose a value of $`V_0`$ to be sufficiently weak that a condensate does exist in a trapped zero-temperature 1D gas, there will always be a regime in which the 3D theories can be applied. As we increase $`V_0`$ the condensate is replaced by a quasicondensate with only short range phase coherence, and this will lead to problems in the higher order theories. However, we have described in the previous section the behaviour that we expect in the high-$`V_0`$ limit. It is therefore of interest to apply the higher order theories to see how and why they fail in this limit, and to gain insight into any improvements needed to describe the 1D gas adequately.

As we have seen, the many-body effects on interactions are fundamental in understanding the behaviour of low dimensional BECs. We must therefore use a theory in which these effects are accounted for accurately. Such a theory is provided by the gapless-Hartree-Fock-Bogoliubov (GHFB) theory Proukakis et al. (1998); Hutchinson et al. (2000). In this theory, the condensate wave function $`\psi _0(x)`$ is described by the generalised GPE in 1D

$`\mu \psi _0(x)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{d^2}{dx^2}}\psi _0(𝐫)+V_{\mathrm{trap}}(x)\psi _0(x)`$ (23)
$`+N_0U_{\mathrm{con}}(x)|\psi _0(x)|^2\psi _0(x)`$
$`+2U_{\mathrm{ex}}(x)\rho (x)\psi _0(x).`$
Here $`\rho (x)`$ is the non-condensate density and it is defined, along with the associated anomalous average $`\kappa `$ \[needed in Eq. (29)\], by

$`\rho (x)`$ $`=`$ $`{\displaystyle \underset{i0}{}}\left(|u_i(x)|^2+|v_i(x)|^2\right)n_i+|v_i(x)|^2,`$ (24)
$`\kappa (x)`$ $`=`$ $`{\displaystyle \underset{i0}{}}u_i(x)v_i^{}(x)\left(2n_i+1\right).`$ (25)
The Bogoliubov functions $`u_i(x)`$ and $`v_i(x)`$ are the solutions of

$$\begin{array}{ccc}(x)u_i(x)+(x)v_i(x)\hfill & =& \hfill \epsilon _iu_i(x),\\ (x)v_i(x)+^{}(x)u_i(x)\hfill & =& \hfill \epsilon _iv_i(x),\end{array}$$
(26)
solved in a basis orthogonal to $`\psi _0(x)`$, and with

$`(x)`$ $`=`$ $`\widehat{H}_{\mathrm{sp}}\mu +2N_0U_{\mathrm{con}}(x)|\psi _0(x)|^2`$ (27)
$`+2U_{\mathrm{ex}}(x)\rho (x),`$
$`(x)`$ $`=`$ $`N_0U_{\mathrm{con}}(x)\psi _0^2(x).`$ (28)
The quasiparticle population factors $`n_i`$ are given by the Bose-Einstein distribution function, and for the zero temperature case we are concerned with here these factors vanish.

This theory distinguishes between the interactions which occur between two condensate atoms (governed by $`U_{\mathrm{con}}(x)`$), and those which occur between a condensate atom and a non-condensate atom (governed by $`U_{\mathrm{ex}}`$). The various higher order theories which we will explore in this section can now be summarised thus:

By setting $`U_{\mathrm{ex}}=0`$ (i.e. neglecting interactions between the condensate and the non-condensate), and using $`U_{\mathrm{con}}(x)=V_0`$, we obtain the familiar GPE and Bogoliubov-de Gennes equations.

The HFB-Popov theory corresponds to the choice $`U_{\mathrm{con}}(𝐫)=V_0`$ and $`U_{\mathrm{ex}}(𝐫)=V_0`$. Thus both condensate-condensate interactions and condensate-non-condensate interactions are included, although many-body effects on the interaction strength are not.

The first version of the GHFB theory is given by setting $`U_{\mathrm{con}}(𝐫)=T_{\mathrm{MB}}(𝐫)`$ \[defined in Eq. (29)\] and $`U_{\mathrm{ex}}(𝐫)=V_0`$. This version of GHFB therefore improves upon the description of the BdG or HFB-Popov descriptions by including many-body effects on the scattering of two condensate atoms.

From the above arguments it is clear that another GHFB theory may be proposed. This consists of setting $`U_{\mathrm{con}}(𝐫)=T_{\mathrm{MB}}(𝐫)`$ and $`U_{\mathrm{ex}}(𝐫)=T_{\mathrm{MB}}(𝐫)`$, which upgrades the description of both condensate-condensate interactions and condensate-non-condensate interactions to the many-body T-matrix.

In the GHFB theories, the many-body T-matrix is given by Proukakis et al. (1998); Morgan (2000)

$$T_{\mathrm{MB}}(x)=V_0\left(1+\frac{\kappa (x)}{N_0\psi _0^2(x)}\right),$$
(29)
and so we can see that the inclusion of the anomalous average leads to the inclusion of the many-body T-matrix in the theories.

Note that we may use the bare potential $`V_0`$ throughout these equations, since we consider the *model* Lieb-Liniger system in which the interaction potential truly is a delta function. Furthermore, in 1D there is an absence of the ultra-violet divergences which occur in other dimensions when a contact potential is used. Since the above equations are always solved self-consistently, we still generate the many-body T-matrix (where used) to all orders from the bare potential Proukakis et al. (1998).

### IV.1 Condensate and quasi-condensate properties

We start our simulations in the weakly-interacting limit with a small value of $`V_0`$. In this limit Lieb and Liniger Lieb and Liniger (1963) found that the Bogoliubov theory gave a good description of the homogeneous system, and so we expect our model should also be valid. We solve the above equations self-consistently to obtain the wave function, chemical potential, many-body T-matrix, $`\rho (x)`$, and $`\kappa (x)`$. We then ramp up the value of $`V_0`$, using the results from the previous iteration as our starting point. By this method we have managed to solve the equations from the weakly-interacting limit into the strongly-interacting limit where the system should become a Tonks-Girardeau gas. Again, we expect the chemical potential to change initially as $`\mu V_0^{2/3}`$ at low $`V_0`$, before reaching a plateau at $`\mu =N_{\mathrm{tot}}\mathrm{}\omega `$ in the strongly-interacting limit.

We present the results from the different forms of the theory below. In all cases we have used $`N_0=15`$, and we start our simulations with $`\stackrel{~}{V_0}=V_0/\mathrm{}\omega \mathrm{}_{\mathrm{trap}}=0.2`$ which is sufficiently weak an interaction that a solution is obtained if $`\rho `$ and $`\kappa `$ are initially zero.

In Fig. 5a the chemical potential has been plotted as a function of the potential strength $`V_0`$, for the various different models, including the GPE results from the previous section. Figure 5b also shows $`\mu /N_{\mathrm{tot}}`$ the chemical potential per atom, in order to compare more easily to the expected high-$`V_0`$ limit $`\mu =N_{\mathrm{tot}}\mathrm{}\omega `$. Note that $`N_{\mathrm{tot}}`$ is the total population of the system, consisting of both the condensate and non-condensate populations. In the low-$`V_0`$ limit it can be seen that all of the models provide very similar results, as should be expected since the gas is weakly-interacting in this limit. As $`V_0`$ is increased, however, the different models show quite different behaviour.

In terms of the chemical potential, neither the BdG nor the HFB-Popov models give the required behaviour, instead $`\mu `$ increases rapidly as $`V_0`$ is increased. In neither model does the graph show a plateau in the high-$`V_0`$ limit, and the chemical potential greatly exceeds the Tonks-Girardeau limit $`\mu =N_{\mathrm{tot}}\mathrm{}\omega `$ in the HFB-Popov case. Instead the chemical potential roughly follows the $`V_0^{2/3}`$ behaviour for all $`V_0`$. The differences between these two models is seen in Fig. 5b, in which the results differ due to the different total populations. In the BdG case the non-condensate population increases very rapidly, such that in the very high-$`V_0`$ limit $`\mu /N_{\mathrm{tot}}`$ is actually decreasing. In the HFB-Popov model, the interactions between the condensate and non-condensate slow the growth of the non-condensate population, but the chemical potential still rises unchecked.

The disagreement of these models with the exact results is not at all surprising. In both models the anomalous average is completely neglected, and so the many-body T-matrix is not introduced at all. We have argued here, and in earlier papers Lee et al. (2002); Lee and Morgan (2002), about the importance of many-body effects in interactions in low-dimensions, and so we expect that neglecting these effects will lead to an inadequate theory. The BdG and HFB-Popov results presented here are similar to those obtained in the homogeneous gas by Lieb and Liniger Lieb and Liniger (1963), who found that the standard Bogoliubov method failed at high interaction strengths, when $`\gamma 2`$. In our results we also see quite good agreement with the GPE of the previous section for $`\gamma _{\mathrm{trap}}1`$ and the differences between the theories become much greater when $`\gamma _{\mathrm{trap}}`$ is increased above this point 222$`\gamma _{\mathrm{trap}}`$ has not been plotted on Fig. 5 since it varies slightly differently for each curve. However a good idea of its approximate value may be obtained by comparing with Fig. 4.. Note that $`\gamma _{\mathrm{trap}}2`$ corresponds to $`\stackrel{~}{V}_010`$ in the figure, which is roughly the point at which the curves separate.

The anomalous average, and therefore the many-body T-matrix, can first be included by solving the GHFB1 equations. Figure 5 shows that there is little difference between these results and those of the HFB-Popov theory. The chemical potential still easily exceeds the limit of the Tonks-Girardeau gas, and so the GHFB1 theory also fails rapidly as $`V_0`$ increases.

However, the GHFB1 results do allow us to examine the form of the many-body T-matrix. In Fig. 6 the spatial dependence of the many-body T-matrix is plotted as $`V_0`$ is increased. In fact, the value being plotted is the T-matrix divided by the condensate density $`|\psi (x)|^2`$. From the previous section we expect that the T-matrix becomes proportional to the condensate density as $`V_0\mathrm{}`$, as in Eq. (16). Figure 6 shows that as the interaction strength is increased this becomes an increasingly good approximation. The very good agreement with the form expected for large $`V_0`$ is further evidence that the coupling parameter $`g_{\mathrm{LL}}(x)`$ derived earlier is accurate in this limit. This verifies our approximation for the many-body T-matrix in terms of the off-shell two-body T-matrix numerically.

That the GHFB1 theory fails in the high-$`V_0`$ limit is not particularly surprising. In the low-$`V_0`$ limit almost all the atoms are in the condensate. As the interaction strength is increased, however, we reach the regime of the quasicondensate discussed by Petrov *et al.* Petrov et al. (2000). As $`V_0`$ is increased the number of atoms in excited states (i.e. $`\rho `$) increases, and can become a sizable proportion of the total number of atoms in the system. It therefore becomes important to describe the interactions between these atoms accurately.

The description of interactions between excited atoms is more complicated than those between condensate atoms. However, two limits may be investigated relatively straightforwardly. If the majority of the collisions occur at a high energy then the correct description would be given by the two-body T-matrix (since the many-body T-matrix tends towards the two-body T-matrix in this limit). On the other hand, if the majority of the collisions are at low energy, then it is more appropriate to describe them by the same many-body T-matrix as used to treat the condensate-condensate interactions. In 2D and 3D it might be expected that the large density of states at high energies should mean that the majority of collisions are in this regime. However, in 1D the density of states is independent of energy in a trap. In the zero-temperature limit in which we are working it therefore seems probable that the appropriate description of the interactions involving non-condensate atoms is the same many-body T-matrix used to describe the condensate interactions.

The many-body T-matrix is used to describe non-condensate interactions in the GHFB2 version of the theory. It can be seen from Fig. 5 that the results using this theory are quite different from those of the GHFB1 theory. In the low-$`V_0`$ limit the results are similar, but as the interactions are increased the chemical potential becomes independent of $`V_0`$, as we expect in the Tonks-Girardeau limit. The figure clearly shows that it reaches an asymptotic value of $`0.7N_{\mathrm{tot}}\mathrm{}\omega `$. This is lower than the expected value of $`N_{\mathrm{tot}}\mathrm{}\omega `$, but nonetheless the results are in qualitative agreement with the behaviour expected for the Tonks gas.

Again, we can compare the spatial shape of the many-body T-matrix with our earlier approximation, and the results look essentially the same as those calculated from the GHFB1 results (Fig. 6). The major change in the GHFB2 theory is therefore to provide a better description of the condensate-non-condensate interactions rather than any implicit change in the form of the many-body T-matrix.

We can also compare the predicted density distribution in the high-$`V_0`$ limit with the solutions of the quintic GPE which was shown to give good agreement with the Tonks-Girardeau results in section III. In Fig. 7 we compare the total densities $`[N_0|\psi (x)|^2+\rho (x)]`$ for the different models with the Thomas-Fermi prediction from Eq. (18). It can easily be seen that the best agreement is given by the GHFB2 theory, while the BdG, HFB-Popov and GHFB1 results all give entirely the wrong shape. It must be said however, that the GHFB2 results are still not in good agreement with the predicted density.

Thus the GHFB2 theory interpreted in the normal manner fails to quantitatively describe the 1D system in the impenetrable limit. This is not surprising since it assumes the existence of a condensate while we know that a condensate does not exist in this limit. Nevertheless, the GHFB2 theory does give qualitatively the correct behaviour for the chemical potential and density in the high-$`V_0`$ limit, which is a significant improvement on the other forms of mean-field theory used for 3D condensates.

The failure of the GHFB1 theory and the relative success of the GHFB2 theory indicate that the description of the non-condensate is the key to solving the 1D problem in the strongly-interacting limit. The GHFB1 model fails at high $`V_0`$ since an increase in $`V_0`$ leads to an increase in $`\rho `$. The $`\rho `$ term obviously becomes far too dominant in GHFB1 and therefore needs to be suppressed by some means. In GHFB2 the many-body T-matrix is introduced to describe interactions involving $`\rho `$, and since $`T_{\mathrm{MB}}/V_0`$ decreases with increasing $`V_0`$ this achieves the required suppression. As confirmation of this argument, the GHFB1 results were rerun with $`\rho `$ artificially set to zero. With this adjustment, the GHFB1 model gave results which looked qualitatively similar to the GHFB2 results, with the chemical potential reaching a plateau in the high-$`V_0`$ limit.

### IV.2 Quasiparticle spectrum

Having seen qualitative agreement between the GHFB2 theory and the expected behaviour of the chemical potential and the density, we turn now to investigate the quasiparticle spectrum. In the high-energy limit the quasiparticles are predicted to be single-particle like, and therefore for a 1D harmonic potential the energy will rise linearly with the quasiparticle quantum number. However, for the low-energy quasiparticle states (with energies less than $`\mu `$), we can use the hydrodynamic theory to obtain an expression for the quasiparticle energies $`\epsilon _n`$ in the Thomas-Fermi limit Stringari (1996). Using the same method as in Stringari (1996), we find that the low-energy quasiparticle spectrum for a BEC in 1D described by a cubic GPE is

$$\epsilon _n=\mathrm{}\omega \sqrt{\frac{n(n+1)}{2}},$$
(30)
where $`n=1,2,\mathrm{}`$ and $`n=1`$ is the Kohn mode. However, we have shown earlier that in the impenetrable limit the system is described instead by a quintic GPE. For such a GPE the same hydrodynamic approach instead predicts Minguzzi et al. (2001)

$$\epsilon _n=n\mathrm{}\omega ,$$
(31)
which is obviously what one would expect for non-interacting fermions, and so we would expect this relationship from a correct theory in the high-$`V_0`$ limit.

Figure 8 shows the GHFB2 results for the energies of the first four quasiparticle modes as a function of the interaction strength. In the weakly-interacting limit the energies are all close to the prediction of Eq. (30), as expected. However, as the interaction strength is increased, the quasiparticle energies do not go over to the limit of Eq. (31). Instead they all fall to around 60% of that value and, even in the high-$`V_0`$ limit, the relationship between the different energies scales more like Eq. (30) than Eq. (31). The conclusion which must be drawn is that the GHFB2 theory does not predict the quasiparticle energies at all accurately in the high-$`V_0`$ limit. This can be seen particularly strongly in the case of the Kohn mode ($`\epsilon _1`$), which should be equal to $`\mathrm{}\omega `$ at all interaction potentials since it simply corresponds to a centre-of-mass oscillation of the entire system in the harmonic potential.

The reason that GHFB2 does not predict the quasiparticle spectrum well is because it does not include any dynamics of the non-condensate. In the strongly-interacting limit a large proportion of the atoms in the system are in the non-condensate, which is assumed to be completely static in all of the theories discussed here. Since the Kohn mode corresponds to the motion of all of the atoms together, it is not surprising that a model which does not include the dynamics of a large proportion of those atoms predicts the wrong frequency.

To see this in more detail, all of the theories discussed here can be derived by starting from the generalised time-dependent GPE

$`i\mathrm{}{\displaystyle \frac{\mathrm{\Psi }(x,t)}{t}}=`$ (32)
$`{\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{^2}{x^2}}\mathrm{\Psi }(x,t)+N_0U_{\mathrm{con}}(x)|\mathrm{\Psi }(x,t)|^2\mathrm{\Psi }(x,t)`$
$`+V_{\mathrm{trap}}(x)\mathrm{\Psi }(x,t)+2U_{\mathrm{ex}}(x)\rho (x)\mathrm{\Psi }(x,t),`$
and expanding $`\mathrm{\Psi }(x,t)`$ in terms of the unperturbed ground-state $`\psi (x)`$ and time-dependent fluctuations $`\mathrm{\Delta }\psi (x,t)`$, and linearising (see, for example, Ruprecht et al. (1996)). The trouble occurs when it is recognised that in 1D the appropriate coupling parameter $`U_{\mathrm{con}}(x)`$ is the many-body T-matrix, as argued earlier. In the high-$`V_0`$ limit, the many-body T-matrix is proportional to $`N_0|\psi (x)|^2`$ in the static case as we have shown. Applying the local density approximation, the time-dependent form of the many-body T-matrix should then be $`N_0|\mathrm{\Psi }(x,t)|^2`$ (assuming that the time-scale for interparticle interactions is much quicker than that for condensate dynamics). In the linearisation procedure, one should therefore also expand the many-body T-matrix as

$$T_{\mathrm{MB}}(x)|\psi (x)+\mathrm{\Delta }\psi (x,t)|^2.$$
(33)
Making this replacement for $`U_{\mathrm{con}}(x)`$ in Eq. (32) while setting $`\rho =0`$ leads directly to the quintic GPE discussed earlier and hence to the result of Eq. (31). We see therefore that including the dynamics of the anomalous average in the GHFB2 theory will lead to the correct prediction for the excitation spectrum in the high-$`V_0`$ limit. For a general value of $`V_0`$, including the fluctuations of the many-body T-matrix requires a dynamic treatment of the anomalous average and hence of the non-condensate generally.

### IV.3 Coherences

It is of interest to calculate the coherence functions for the Lieb-Liniger gas from the GHFB2 results. We have argued that we are dealing with a gas in which a condensate may form in the low-$`V_0`$ limit, but where it disappears when the interactions are increased. Since coherence is a major feature of a condensate we should be able to see a loss of coherence as $`V_0`$ is increased, and again we can try to relate this to the expected results in the Tonks limit.

#### IV.3.1 First-order coherence

The first order coherence function $`G^{(1)}(𝐫,𝐫^{})`$ is defined in terms of the quantum field operators by Naraschewski and Glauber (1999)

$$G^{(1)}(𝐫,𝐫^{})=\widehat{\mathrm{\Psi }}^{}(𝐫)\widehat{\mathrm{\Psi }}(𝐫^{}),$$
(34)
and the normalised first order coherence function is therefore

$$g^{(1)}(𝐫,𝐫^{})=\frac{G^{(1)}(𝐫,𝐫^{})}{\sqrt{G^{(1)}(𝐫,𝐫)G^{(1)}(𝐫^{},𝐫^{})}}.$$
(35)
Decomposing the field operator in terms of the condensate and its fluctuations leads to an expression for the coherence function at $`T=0`$ in terms of the condensate wave function and Bogoliubov quasiparticle amplitudes Dodd (1997)

$$G^{(1)}(𝐫,𝐫^{})=N_0\psi _0^{}(𝐫)\psi _0(𝐫^{})+\underset{i0}{}v_i(𝐫)v_i^{}(𝐫^{}).$$
(36)
Figure 9 shows the normalised first order coherence function obtained from the GHFB2 results for various values of $`V_0`$. Note that the figure has been plotted in terms of the distance scale $`\stackrel{~}{x}=x/R_{\mathrm{TF}}`$, where $`R_{\mathrm{TF}}`$ is the Thomas-Fermi radius of the condensate, in order to take into account the greater spatial extent of a condensate with stronger interactions. Figure 9 clearly shows that the system is coherent over its entire length in the weakly-interacting limit. However, as the interaction strength is increased the coherence can be seen to diminish substantially, and the width of the central maximum decreases as well. A diminishing range of coherence is consistent with the absence of a true condensate in the strongly-interacting regime. However, the rate at which the coherence calculated here decreases is surprisingly slow. Gangardt and Shlyapnikov Gangardt and Shlyapnikov (2002) have presented results for the coherence function obtained from hydrodynamic theory applied to a 1D gas, and they see a much greater fall in the range of the coherence length - dropping by almost an order of magnitude as $`\gamma _{\mathrm{trap}}`$ decreases from $`0.1`$ to $`10`$. However, their results were obtained for a system of $`10^4`$ atoms, many more than we deal with here. Also, even for the largest value of $`V_0`$ plotted, the condensate depletion is still only about $`50`$%, in our case.

#### IV.3.2 Second-order coherence

We now look briefly at the second-order coherence function defined by Naraschewski and Glauber (1999)

$$g^{(2)}(𝐫,𝐫^{})=\frac{\widehat{\mathrm{\Psi }}^{}(𝐫)\widehat{\mathrm{\Psi }}^{}(𝐫^{})\widehat{\mathrm{\Psi }}(𝐫^{})\widehat{\mathrm{\Psi }}(𝐫)}{\widehat{\mathrm{\Psi }}^{}(𝐫)\widehat{\mathrm{\Psi }}(𝐫)\widehat{\mathrm{\Psi }}^{}(𝐫^{})\widehat{\mathrm{\Psi }}(𝐫^{})}.$$
(37)
We will consider this function only in the equal-position limit (i.e. when $`𝐫=𝐫^{}`$), in which case it may be written as

$`g^{(2)}(x,x)=`$
$`1+{\displaystyle \frac{2N_0|\psi (x)|^2[\rho (x)+\kappa (x)]+|\rho (x)|^2+|\kappa (x)|^2}{\left[N_0|\psi (x)|^2+\rho (x)\right]^2}},`$
following a similar argument as for the first-order coherence function. The second-order equal-position coherence function is of interest because of the various limits in which exact results are known. For a pure BEC with perfect second-order coherence we have $`g^{(2)}(x,x)=1`$, whilst for a thermal cloud of bosons $`g^{(2)}(x,x)=2`$ Dodd et al. (1997); Naraschewski and Glauber (1999). In the strong-coupling limit in 1D however, the Bose-Fermi mapping theorem predicts that $`g^{(2)}(x,x)=0`$ Girardeau et al. (2001) for impenetrable bosons. We would therefore expect to find a second-order coherence of close to $`1`$ in the weakly-interacting limit, which decreases to $`0`$ as the interaction strength is increased, in contrast to the higher dimensional case where the coherence function would increase to $`2`$ as the non-condensate became significant 333This is not entirely correct, in fact, since we have ignored the strong repulsion between atoms at very small distances, which means that $`g^{(2)}(𝐫,𝐫^{})`$ must actually vanish for $`𝐫=𝐫^{}`$ in all dimensions. The actual function $`g^{(2)}(𝐫,𝐫^{})`$ increases quite rapidly, however, as a function of $`|𝐫𝐫^{}|`$ (on the scale of the scattering length), increasing to values close to those given in the text. Equation (IV.3.2) is therefore the result for an effective Hamiltonian involving a pseudopotential interaction term, rather than the bare interparticle potential. This is discussed in more detail in Naraschewski and Glauber (1999)..

Figure 10 shows $`g^{(2)}(x,x)`$ as calculated from the results of the GHFB2 theory. The results can be seen to be close to unity, as expected, for the weakly-interacting case (and rising to $`2`$ outside the region where the condensate exists). As the interactions are increased, $`g^{(2)}(x,x)`$ initially decreases closer to the $`g^{(2)}(x,x)=0`$ limit expected. However, in the very strongly-interacting limit the coherence function rises again to become greater than $`1`$. Thus the GHFB2 theory, whilst in qualitative agreement with the behaviour of the chemical potential and the density, does not give the expected results for the coherence of the system in the impenetrable limit.

## V Reinterpreting the results

As seen in the previous section, although the GHFB2 theory seems to go further than other standard theories in describing the 1D Bose gas into the strongly-interacting limit, there are significant differences between the results in the impenetrable limit and the exact results of the Tonks-Girardeau gas. However, we will see in this section that the agreement can be improved significantly by reinterpreting the results in a different manner. The disagreement with the expected values for the asymptotic value of the chemical potential and with the shape of the density are both related to the non-condensate population in the form of $`\rho (x)`$. Looking again at Fig. 5 one can see that the chemical potential in the GHFB2 case rises to an asymptotic value of approximately $`15\mathrm{}\omega `$. Although this result was obtained for a condensate population of $`N_0=15`$, it is in excellent agreement with the expected result for a system with $`15`$ atoms in total. By itself this may appear to be a co-incidence; however the shape of the density resulting merely from the condensate “wave function” $`\psi (x)`$ is plotted against the Thomas-Fermi exact prediction in Fig. 11, and the two shapes can be seen to agree quite well. It seems that the best agreement would be midway between that in Fig. 11, and that in Fig. 7b (where the total density has been plotted).

This seems to imply that the GHFB2 theory predicts the correct results in the strongly-interacting limit if one partly ignores the contributions from the non-condensate (these contributions are sizeable - the non-condensate population for the results plotted in Fig. 11 is about $`N_0/2`$). This would have to be combined with a re-interpretation of the condensate density as the quasi-condensate density and, ultimately, as the total density. Some procedure of this kind is clearly required in 1D because the various theories discussed in the previous section are based on the existence of a condensate whereas in fact the condensate disappears as $`V_0`$ increases.

As another illustration of the good agreement with exact results which can be obtained with such a reinterpretation, in Fig. 12 the second-order equal-position coherence function is plotted with results obtained from the GHFB1 theory with $`\rho (x)`$ artificially set to zero (as mentioned earlier, this approach gave qualitative agreement with the expected behaviour of the chemical potential). The results in this case clearly show the behaviour predicted in the previous section, with $`g^{(2)}(x,x)`$ decreasing from unity to zero as the interaction strength is increased. Such a theory (with $`\rho (x)=0`$) would result from a reinterpretation of $`N_0|\psi (x)|^2`$ as the total density of the system, rather than merely that of the condensate. The agreement with the expected behaviour for $`g^{(2)}(x,x)`$ can easily be seen from the following argument. In the high-$`V_0`$ limit we have shown that $`T_{\mathrm{MB}}`$ becomes independent of $`V_0`$. From Eq. (29), this must mean that $`\kappa (x)N_0\psi _0^2(x)`$ in the same limit. Equation (IV.3.2) then clearly shows that the second order coherence function vanishes in the high-$`V_0`$ limit if $`\rho =0`$.

That such a reinterpretation of the condensate is necessary and will lead to the correct mean-field theory at large $`V_0`$ can be seen by comparison with the GPE results in section III.3. If we set $`\rho =0`$ and reinterpret $`n_0(x)`$ as $`n_{\mathrm{tot}}(x)`$ then the GHFB theory reduces to the earlier GPE theory, the only difference being that $`g_{1\mathrm{D}}(x)`$ is calculated from Eq. (29) instead of Eq. (III.3). Assuming the validity of the local density approximation, however, these two definitions will lead to the same results. We note that a theory has recently been proposed by Anderson *et al.* Andersen et al. (2002); Khawaja et al. (2002) which in effect leads to a similar reinterpretation of the quantities $`\rho (x)`$, $`\kappa (x)`$ and $`n_0(x)`$ by considering the effects of phase fluctuations to all orders, and this approach may provide a microscopic justification for this procedure.

## VI Conclusions

In this paper we have given a form for an approximate many-body T-matrix to describe the scattering which occurs in a zero temperature 1D Bose gas. This approximate T-matrix enables the 1D Bose system to be described by a Gross-Pitaevskii equation which is easily solved. We have shown that the solutions to this equation for a trapped system are in excellent agreement with the density of a 1D Lieb-Liniger model gas in the two limits in which exact results are known, and that the density distribution of such a gas changes markedly as the interactions become stronger.

We have also attempted to push various higher-order theories of BECs from the weakly-interacting regime into the strongly-interacting regime for the Lieb-Liniger gas. We have shown that the theories which neglect many-body effects on the interactions (BdG and HFB-Popov) fail to provide adequate results at even quite low interaction strengths. The GHFB1 theory includes the many-body T-matrix for condensate interactions, but it also fails to describe the strongly-interacting limit since it does not account for the condensate-non-condensate interactions properly and hence the terms involving $`\rho `$ become too large.

However, the GHFB2 theory, which includes many-body effects in the condensate-non-condensate interactions, gave results which were qualitatively correct even in the strongly-interacting limit, predicting the correct behaviour of the chemical potential and density profile. Quantitative agreement was not found, but the results show that GHFB2 describes the strongly-interacting limit to a much better degree than the other mean-field theories tested.

GHFB2 theory still fails to adequately describe some aspects of the system however. The first-order coherence function shows qualitatively the behaviour expected, but the second-order coherence, while initially showing the correct behaviour in the low- to mid-$`V_0`$ range, deviates from the expected limit at high-$`V_0`$. Furthermore, the quasiparticle excitation spectrum does not go over to the correct behaviour at high-$`V_0`$.

The quasiparticle spectrum can be corrected if the theory is modified to allow for fluctuations in the many-body T-matrix, as well as in the condensate density. This requires a dynamical treatment of the non-condensate. At the same time, it seems that the key to pushing the GHFB2 theory into the high-$`V_0`$ limit (obtaining the correct coherences and quantitative agreement in general) lies in a reinterpretation of the non-condensate density as the total density, and we have presented some heuristic evidence that this will lead to the expected results.

###### Acknowledgements.

This research was supported by the Engineering and Physical Sciences Research Council of the United Kingdom, and by the European Union via the “Cold Quantum Gases” network. S. M. thanks Trinity College, Oxford, and the Royal Society for financial support. K. B. thanks the Royal Society and the Wolfson Foundation for support.

## Appendix A The many-body T-matrix in 3D

As an illustration of our approximation of the many-body T-matrix in terms of the off-shell two-body T-matrix, we show in this appendix that in the three-dimensional case it reproduces the known results for the many-body T-matrix.

The 3D case differs from 1D and 2D because the many-body T-matrix is only a small correction to the two-body T-matrix. As mentioned earlier, the 3D two-body T-matrix is a constant $`U_0=4\pi \mathrm{}^2a/m`$ to first order in the low-energy limit. The first approximation to the 3D coupling parameter is therefore $`g_{3\mathrm{D}}=U_0`$. We can then ask if the method used in this paper for 1D (and used in reference Lee et al. (2002) for 2D) accurately predicts the effects of the medium on the scattering.

In Lee et al. (2002) we showed that our approximation for the many-body T-matrix in 3D is

$$g_{3\mathrm{D}}=0|T_{\mathrm{MB}}(E=0)|00|T_{2\mathrm{B}}(E=\frac{16}{\pi ^2}\mu )|0,$$
(39)
and the off-shell two-body T-matrix for a gas of hard-spheres of radius $`a`$ was shown in Morgan et al. (2002) to be

$`\mathrm{𝟎}|T_{2\mathrm{B}}^{3\mathrm{D}}(E)|\mathrm{𝟎}=`$ (42)
$`\{\begin{array}{cc}U_0\left[1i\sqrt{\frac{mE}{\mathrm{}^2}}a\frac{1}{3}\frac{mE}{\mathrm{}^2}a^2\right],\hfill & \text{for }E>0,\hfill \\ U_0\left[1+\sqrt{\frac{m|E|}{\mathrm{}^2}}a\frac{1}{3}\frac{mE}{\mathrm{}^2}a^2\right],\hfill & \text{for }E<0,\hfill \end{array}`$
Combining these equations gives (to first order in the small parameter $`k_\mu a=\sqrt{m\mu a^2/\mathrm{}^2}`$)

$$g_{3\mathrm{D}}=U_0\left[1+\left(\frac{16\mu ma^2}{\pi ^2\mathrm{}^2}\right)^{1/2}\right]+\text{O}[(k_\mu a)^2].$$
(43)
We can rewrite this in terms of the density $`n_0`$, since in a homogeneous system $`\mu =g_{3\mathrm{D}}n_0`$ to leading order. Making this substitution and solving for $`g_{3\mathrm{D}}`$ gives the result

$$g_{3\mathrm{D}}=U_0\left[1+\frac{8}{\sqrt{\pi }}(n_0a^3)^{1/2}\right]+\text{O}(n_0a^3),$$
(44)
where we have expanded in terms of the dilute gas parameter $`n_0a^3`$, which is small for a dilute gas BEC. This result agrees with other solutions Shi and Griffin (1998); Morgan (1999) for the many-body T-matrix in 3D homogeneous systems. At $`T=0`$ we can, therefore, account for many-body effects on particle scattering using an off-shell two-body T-matrix.


format: text

text:
# text

## 1 Introduction

The essential postulates of the deBroglie Bohm causal theory  can be summarized as:

* An individual physical system comprises a wave together with a point particle. The wave and the particle are distinct parts of the same system. The wave and the particle are never separated.
* The wave $`\psi `$ is a solution to the Schrödinger equation.
* Particle motion is obtained as the solution to a modified Hamilton - Jacobi equation. The modification to the classical Hamilton - Jacobi equation is the addition of a term Q, the quantum potential. The functional form of Q is determined mathematically by $`\psi `$.

The causal theory has many attractive features. For example:

* It is deterministic.
* It allows for the possibility of a complete description of a physical system, i.e., it addresses one of the fundamental complaints Einstein had about the Copenhagen interpretation \[Ref. 2, pp. 11-15\].
* Particle dynamics is written as the sum of classical and quantum terms \[Ref. 1, p. 29\].
* There is a natural relationship between the classical and quantum regimes \[Ref. 2, pp. 270-274\].

Together these allow for the possibility that the Causal Theory may be able to address some of the deepest problems associated with the Copenhagen interpretation. The theory may also lead to new physics.

There are, however, weaknesses in the original theory. One of the most obvious of these relates to the quantum potential Q: What is its source? Typically in physics a force, and its associated potential, have a source. However, nowhere in the literature is this fundamental question addressed in a physically reasonable way.

To compound the problem, the adherents of the theory have gone further and made conjectures leading to results which are in direct contradiction with experiment.

These conjectures, which we will discuss later, have now become part of the orthodox view of the deBroglie Bohm causal interpretation. Unfortunately, they have done much to undermine the appealing aspects of the original theory.

Our approach is as follows:

1. Return to the original theory as stated by deBroglie and Bohm.
2. Reject all unphysical conjectures.
3. Let experiment guide us as much as possible.
4. For those things about which Copenhagen makes a prediction, our theory must be in agreement.
5. Start with a many-body approach because it is closer to physical reality than a one-body approach.
6. Address two fundamental problems: What is the physical reason for the forces associated with the quantum potential? What is the source of these quantum forces?

Our goal is to then go on in future work to see if the causal theory can be used to address some of the compelling questions of the our day.

## 2 Basics Of The Causal Theory

We start by considering a physical system consisting of n interacting particles. The time dependent Schrödinger equation for this system is,

$`ı\mathrm{}\left({\displaystyle \frac{\psi }{t}}\right)=\left[{\displaystyle \underset{i=1}{\overset{n}{}}}\left({\displaystyle \frac{\mathrm{}^2}{2m_i}}\right)_i^2+V\right]\psi `$ (1)
where

$$\psi =\psi (\text{x}_1,\text{x}_2,\text{x}_3,\mathrm{}\text{x}_n,t),$$

$$V=V(\text{x}_1,\text{x}_2,\text{x}_3,\mathrm{}\text{x}_n,t),$$

$$_i=(\frac{}{x_i},\frac{}{y_i},\frac{}{z_i}),$$

and $`(\text{x}_1,\text{x}_2,\text{x}_3,\mathrm{}\text{x}_n)`$ provides a set of rectangular Cartesian coordinates of the n particles. Also V is the classical potential energy which includes the interparticle and external potentials from the traditional forces (gravity, E&M, strong and weak interactions.)

Assume that the energy spectrum is discrete, so that the wavefunction $`\psi `$ is localized in configuration space. Without loss of generality, the wavefunction can be written as

$`\psi =R\mathrm{exp}ı(S/\mathrm{})`$ (2)
where R and S are real and

$$R=R(\text{x}_1,\text{x}_2,\text{x}_3,\mathrm{}\text{x}_n,t)$$

$$S=S(\text{x}_1,\text{x}_2,\text{x}_3,\mathrm{}\text{x}_n,t)$$

Equation (1) is equivalent to the pair of equations

$`{\displaystyle \frac{S}{t}}+{\displaystyle \underset{i=1}{\overset{n}{}}}(_iS)^2(2m_i)^1+Q+V=0`$ (3)
$`{\displaystyle \frac{R^2}{t}}+{\displaystyle \underset{i=1}{\overset{n}{}}}_i\left({\displaystyle \frac{R^2_iS}{m_i}}\right)=0`$ (4)
where

$$Q=Q(\text{x}_1,\text{x}_2,\text{x}_3,\mathrm{}\text{x}_n,t)$$

$`Q={\displaystyle \underset{i=1}{\overset{n}{}}}\left({\displaystyle \frac{\mathrm{}^2}{2m_iR}}\right)_i^2R`$ (5)
is defined as the quantum potential \[Ref. 2, p. 279\]. If it were not for the presence of the quantum potential Q, Eq. (3) would be the classical Hamilton - Jacobi equation.

Just as the causal form of the Hamilton - Jacobi equation contains the additional term Q, the causal form of Newton’s second law contains an additional term involving Q \[Ref. 2, p. 279-280\],

$`{\displaystyle \frac{d_i}{dt}}=_iV_iQ`$ (6)
here $`\text{P}_i`$ is the momentum of the i-th particle, $`V_i`$ is the sum of all the conventional forces on the i-th particle, and $`_iQ`$ is interpreted as the quantum force on the i-th particle. The total momentum of our n-particle system is given by

$`={\displaystyle \underset{i=1}{\overset{n}{}}}_i`$ (7)
with dynamical equation \[Ref. 2, p. 285\]

$`{\displaystyle \frac{d}{dt}}={\displaystyle \underset{i=1}{\overset{n}{}}}_iV{\displaystyle \underset{i=i}{\overset{n}{}}}_iQ`$ (8)
## 3 Studies of the quantum potential

The purpose of this paper is to attempt to answer the following questions. What is the physical reason for the forces associated with the quantum potential? Furthermore, what is the source of these forces? We suggest that Q is the result of a nonholonomic nonlocal constraint on the system, the requirement that the wavefunction $`\psi `$ not exhibit deterministic chaos (extreme sensitivity to initial conditions.) The forces associated with Q are forces of constraint. And we furthermore suggest that the source of the quantum force on a given particle is the other particles of our n particle system.

Why are we interested in preventing deterministic chaos? In general, quantum mechanical systems localized in configuration space are described by Hermitian Hamiltonians. With such a Hamiltonian, the wavefunction is never subjected to extreme sensitivity to initial conditions. However, removing Q from the causal Hamilton-Jacobi equation (Eq. 3) opens the possibility that the corresponding wavefunction be subjected to extreme sensitivity to initial conditions. All this is shown in Appendix A. In Appendix B ways to detect the presence of deterministic chaos are discussed.

While deterministic chaos of the wavefunction should not occur in a localized quantum system, it is quite possible that such chaos can be associated with the particle trajectories of the system, as calculated by deBroglie Bohm theory . Furthermore, in a macroscopic quantum system the wavefunction may exhibit deterministic chaos, as is discussed in Appendix C.

We now turn to our second question: What is the source of the quantum forces? On this issue the leading proponents of the causal theory make an assumption (conjecture) which we find truly mysterious. Consider the following quotes:

1. On page 170 of Ref.3, Bohm states that the field $`\psi `$ “exerts a force on the particle in a way that is analogous to, but not identical with, the way in which an electromagnetic field exerts a force on a charge, and a meson exerts a force on a nucleon111It should be noted that Bohm and Hiley \[Ref. 1, pp. 29, 30, 37, 38, 40\] imply, in contrast, that in our Eq. (6) , $`_iQ`$ is a force being exerted on the i-th particle, but that it is not directly exerted by the wavefunction $`\psi `$. Rather, the wavefunction controls the action and reaction of the particle with something like the fluctuations of the vacuum, causing “acceleration of the particle in its self-movement.”.”
2. On page 91 of Ref. 2, Holland states “it should be emphasized that while the quantum field \[ $`\psi `$ \] does not push on the particle as we might expect a classical wave to, it does nevertheless guide the particle by exerting a direct force on it via Q.”
We infer from these statements that Bohm and Holland make the assumption that the wavefunction, not only is an independent entity distinct from the particle itself, but also is the ”source” of the quantum force! Consequently, it is no surprise when Holland goes on to say :
3. Ref. 2, p. 79: “the particle simply responds to local values of the field in the vicinity (via Q) - there is no reciprocal action of the particle on the wave.”
4. Ref. 2, p. 286: “A classically closed system may not be closed quantum mechanically.”

$$\frac{d\text{P}}{dt}=\underset{i=1}{\overset{n}{}}_iQ0$$

The assumption that the source of the quantum force is the wavefunction itself has led the orthodox proponents of the theory to unphysical results. For example, the third quote above suggests that Newton’s third law (to every action there is an equal and opposite reaction) is violated. Furthermore, the fourth quote allows for the possibility that an isolated system may self-accelerate in the absence of any known force. Not only are these results in disagreement with predictions made using the Copenhagen interpretation of quantum mechanics, but they are also in disagreement with all known experiments.

These unphysical results are a consequence of combining the assumption that “the wavefunction itself exerts a force on the particle” with the fact that the Schrödinger equation is homogeneous. Rather than questioning the assumption on the wavefunction, it is apparently the orthodox deBroglie Bohm approach to accept it and instead speculate on the possibility of some inhomogeneous version of the Schrödinger equation222 Bohm \[Ref. 3, p. 179\] suggests that at a distance less than $`10^{13}`$ cm, Schrödinger’s equation becomes inhomogeneous allowing a particle to react back on the field, in analogy with the inhomogeneous Maxwell’s equations allowing a charge to act on the electromagnetic field.. There is, however, no experimental evidence to support these speculations.

We take a different approach. To obtain a physically reasonable causal theory we add two assumptions to the original deBroglie Bohm postulates:

1. The wavefunction exerts no forces on the particles of the system. Although it is apparent from Eq. (1) through (5) that the wavefunction determines the functional form of the quantum potential, there is no indication that $`\psi `$ itself exerts forces on the particles of the system.
2. For an isolated system, the total momentum P is a constant of the motion; i.e., it does not change with time. This is consistent with experiment. It is also consistent with the Copenhagen interpretation; because for an isolated system, the wavefunction is an eigenfunction of the total momentum operator.

Using these ideas, and returning to our many-body equations, we are led to a very different result for the source of the quantum force than that in the orthodox deBroglie Bohm view.

Consider the situation in which our n-particle system is an isolated system. Since the traditional forces obey Newton’s third law we have

$$\underset{i=1}{\overset{n}{}}_iV=0$$

For an isolated system we have

$$\frac{d\text{P}}{dt}=0$$

so that

$$\underset{i=1}{\overset{n}{}}_iQ=\underset{i=1}{\overset{n}{}}\text{F}_i=0$$

That is, the sum of all the quantum forces on the n-particles is zero. From this we have

$`\text{F}_i={\displaystyle \underset{ji}{\overset{n}{}}}(\text{F}_j)`$ (9)
$`\text{F}_i`$ is the quantum force on the i-th particle. Eq. (9) states that the quantum force on the i-th particle of our isolated system is equal and opposite to the sum of the quantum forces on all the other particles. For n = 2, this gives $`\text{F}_1=\text{F}_2`$, i.e., the quantum force on particle 1 is equal and opposite to the quantum force on particle 2. This strongly suggests that the source of the quantum force on one particle is the other particle. (Which, conveniently, is also consistent with Newton’s Third Law.)

Eq. (9) implies that the net quantum force on the i-th particle is the result of all the other particles of the system exerting forces on the particle through the intermediary of the quantum potential. If the system is acted on by external forces, the simplest thing to do conceptually is to enlarge the system to include the bodies creating the external forces, so that one has an enlarged isolated system. We see that the source of the quantum force on any particle is nothing more that all the other particles of the system. To us, this seems to be a much more satisfying explanation than is that of Holland or of Bohm and Hiley \[Ref. 1, pp. 30, 37, 38, 40\]. The nonlocality of the quantum potential discussed by Holland \[Ref. 2, p. 282\] is the cause of the nonlocal interactions between the particles of the system.

Consider the situation of an isolated hydrogen atom in the ground state. The total force on the atom vanishes. At the same time there are equal and opposite Coulomb forces of attraction on the electron and on the proton. Since the internal wavefunction of the system is real, the causal theory says that the electron is motionless relative to the proton. This results from the fact that there are equal and opposite quantum forces of repulsion between the two particles, such that there is no net force on either particle; and their relative position stays fixed.

With this understanding of the physical reason for, and the source of, the quantum forces, the authors feel that the concept of the quantum potential is less strange and more reasonable than would appear at first sight.

## 4 Acknowledgments

The one of the authors (ALD) would like to thank Harold Brooks and Richard Shapiro for several helpful discussions.

## 5 Appendix A

When we remove Q from the Hamilton - Jacobi equation \[Eq. (3)\], we are adding to the Schrödinger equation the term \[Ref. 2, p. 56\]

$`\left({\displaystyle \frac{\mathrm{}^2}{2m}}\right)\mathrm{exp}ı(S/\mathrm{})^2R`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{}^2}{2m}}\right)\psi ^1\psi ^2\psi `$ (10)
$`=`$ $`Q\psi `$
so that the effective Hamiltonian becomes

$`H_{eff}=\left({\displaystyle \frac{\mathrm{}^2}{2m}}\right)^2+V+\left({\displaystyle \frac{\mathrm{}^2}{2m}}\right)\psi ^1^2\psi `$ (11)
Note that $`H_{eff}=H_{eff}(\psi )`$ is a functional of $`\psi `$, the state upon which it is operating. This means that the superposition principle no longer holds. When $`\phi \psi `$,

$`{\displaystyle (\phi ^{}H_{eff}\psi \psi H_{eff}^{}\phi ^{})𝑑\tau }`$
$`=`$ $`\left({\displaystyle \frac{\mathrm{}^2}{2m}}\right){\displaystyle \phi ^{}\psi \left[\psi ^1^2\psi \phi ^1^2\phi \right]𝑑\tau }`$
$``$ $`0`$
so that $`H_{eff}`$ is not Hermitian, $`H_{eff}H_{eff}^{}`$. This means that the time development operator $`\mathrm{exp}[(ı/\mathrm{})H_{eff}t]`$ is not unitary. The time dependent Schrödinger equation is a nonunitary flow. Since

$`ı\mathrm{}\left({\displaystyle \frac{}{t}}\right)\psi ^{}\psi =\psi ^{}H_{eff}\psi \psi H_{eff}\psi ^{}`$ (13)
we have

$`ı\left({\displaystyle \frac{}{t}}\right){\displaystyle \psi ^2𝑑\tau }=`$
$`\left({\displaystyle \frac{\mathrm{}}{2m}}\right){\displaystyle \psi ^2\left[\psi ^1^2\psi \psi ^1^2\psi \right]𝑑\tau }`$
and the normalization of the wavefunction is time independent. When $`\phi \psi `$

$`ı\left({\displaystyle \frac{}{t}}\right){\displaystyle \phi ^{}\psi 𝑑\tau }`$
$`=\left({\displaystyle \frac{\mathrm{}}{2m}}\right){\displaystyle \phi ^{}\psi \left[\psi ^1^2\psi \phi ^1^2\phi \right]𝑑\tau }`$
$`=0`$
so that

$`\left({\displaystyle \frac{}{t}}\right){\displaystyle \psi \phi ^2𝑑\tau }`$
$`=\left({\displaystyle \frac{\mathrm{}}{2m}}\right){\displaystyle ı\left[\psi ^{}\phi \phi ^{}\psi \right]\left[\psi ^1^2\psi \phi ^1^2\phi \right]𝑑\tau }`$
$`0`$
The normalization of $`(\psi \phi )`$ is real and time dependent.

Consider the case where two initial conditions for the time dependent Schrödinger equation differ only infinitesimally. As time progresses the two corresponding wavefunctions can become quite different. This indicates the possibility of deterministic chaos of the wavefunction (extreme sensitivity to initial conditions.) All this is a consequence of $`H_{eff}`$ being a functional of the state upon which it is acting.

If we remove the term

$$\left(\frac{\mathrm{}^2}{2m}\right)\psi ^1^2\psi $$

from Eq. (11) , we are left with a Hermitian Hamiltonian, as demonstrated by removing the corresponding terms in Eq. (12). Removing the corresponding terms from Eqs. (13) and (14) , we see that the normalization of $`(\psi \phi )`$ is time independent. This shows that there can be no deterministic chaos of the wavefunction associated with a Hermitian Hamiltonian.

## 6 Appendix B

In order to detect the presence of deterministic chaos, it is necessary to transform the time dependent Schrödinger equation so that it looks like a set of classical equations of motion . By expanding the wavefunction $`\psi (t)`$ in terms of the members of some orthonormal complete set $`\phi _i`$

$`\psi (t)={\displaystyle \underset{j}{}}a_j(t)\phi _j`$ (17)
the Schrödinger equation becomes

$`{\displaystyle \frac{da_j}{dt}}={\displaystyle \underset{k}{}}M_{jk}a_k`$ (18)
with

$`M_{jk}(ı\mathrm{})^1jHk.`$ (19)
This is a complex N dimensional flow, where N is the dimension of the Hilbert space appropriate to the problem. (The equivalent real flow is 2N dimensional.) Under ordinary circumstances the flow is linear and there is no reason to expect deterministic chaos. However, if the Hamiltonian is a functional of the state of the system, then $`M_{jk}`$ will be a function of the coefficients $`a_i`$ and the flow is nonlinear. Thus, if $`2N>3`$, the Poincare - Bendixson theorem  indicates that deterministic chaos may occur. While numerically integrating Eq. (18), one can simultaneously calculate the Lyapunov exponents appropriate to the solution . If the largest exponent is positive, there is deterministic chaos.

## 7 Appendix C

Consider the quantum mechanical treatment of a macroscopic system. By definition, such a system is sufficiently large that intensive properties are invariant to size increase. This suggests considering such systems in the limit of infinite size and constant density, know as the thermodynamic limit, where phase transitions are sharp (not broadened) and Poincare recurrences do not occur. In the last three decades there has been considerable progress in the algebraic treatment of infinite systems, both with regard to quantum field theory and with regard to quantum statistical mechanics .

There can be many representations of the quantum mechanical operators corresponding to the observables of a physical system. For a finite system, all representations are unitarily equivalent, so that the choice of representation is a matter of convenience. In contrast, for an infinite system there can be many unitarily inequivalent representations. Each representation can have associated with it a complete set of microscopic states and microscopic variables. The various inequivalent representations can be distinguished by the values of intensive macroscopic variables associated with them. For example, in a ferromagnet the components of the magnetization vector are macroscopic variables. There is a three-fold continuous infinity of inequivalent representations associated with the values of these components.

For our purposes, a most important property is that the Hamiltonian is a functional of the representation \[ Ref. 9, p. 32\]. A classic example is the BCS theory of superconductivity , where the so-called reduced Hamiltonian contains the complex superconducting order parameter $`\mathrm{\Delta }`$ , whose magnitude is a measure of the number of Cooper pairs in the system (and thus a function of temperature) and whose phase plays a crucial role in the Josephson effect .

Since the Hamiltonian is a functional of the representation, and thus the state of the system, the matrix elements $`M_{jk}`$ of Eq. (19) are functions of the coefficients $`a_i`$ . Thus the equations of motion, Eq. (18), are nonlinear in the coefficients, and deterministic chaos may occur. More specifically, the Hamiltonian H and the matrix elements $`M_{jk}`$ are functions of n macroscopic variables indexing the various inequivalent representations, the macroscopic variables in turn being functions of the $`a_i`$ . Thus the N dimensional flow of Eq. (18) can be used to generate an n dimensional nonlinear flow involving the n macroscopic variables. The Poincare-Bendixson theorem  requires $`n3`$ for the possibility of chaos.

Since the Hamiltonian is a functional of the state of the system, the superposition theorem breaks down. This may have something to do with the fact that no one has ever observed a coherent superposition of the two macroscopic states of Schrödinger’s cat .

## 8 References

 D. Bohm and B. J. Hiley, The Undivided Universe (Routledge, London 1993).

 P. R. Holland, The Quantum Theory of Motion (Cambridge University Press,Cambridge 1993).

 D. Bohm, Phys. Rev. 85, 166 (1952).

 R. H. Parmenter and R. W. Valentine, Phys. Lett.  A 227, 5 (1997).

 R. H. Parmenter and L. Y. Yu, Physica D 80, 289 (1995).

 M. W. Hirsch and S. Smale, Differential Equations, Dynamical Systems, and Linear Algebra (Academic Press, New York 1974).

 A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vistano, Physica D 16, 285 (1985).

 D. Ruelle, Statistical Mechanics (W. A. Benjamin, New York 1969), p. 2.

 G. L. Sewell, Quantum Theory of Collective Phenomena (Oxford University Press,Oxford 1986).

 J. Von Neumann, Math. Annalen 104, 570 (1931).

 J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 108, 1175 (1957).

 B. D. Josephson, Adv. Phys. 14, 410 (1965).

 E. Schrödinger in Quantum Theory and Measurement, edited by J. A. Wheeler and W. H. Zurek (Princeton University Press, Princeton 1983).


text:
# The hexatic phase of the two-dimensional hard disk system

## Acknowledgments

We thank Kurt Binder and Lothar Schülke for helpful comments on our draft and Innovista in Hannover for providing computer time. We also wish to thank Claudia Rivière for correcting grammatical errors and typos. Especially we benefited from discussions with Conny Tollp. Deiters.


text:
# Robust Steering of n-level Quantum Systems

## 1 Introduction

We consider an isolated n-dimensional quantum system with time evolution described by the following Schrödinger equation:

$`i\mathrm{}|\dot{\psi }(t)`$ $`=`$ $`H(t)|\psi .`$ (1)
Here $`|\psi (t)`$ is a vector of unit norm in $`^n`$ representing the state of the system at time $`t`$. The unitary time evolution of the system is governed by the system Hamiltonian:

$$H(t)=H_0+\underset{j=1}{\overset{m}{}}H_ju_j(\theta ,t).$$
(2)
The *internal Hamiltonian* $`H_0^{n\times n}`$ is an Hermitian matrix describing the free evolution of the system. The *control Hamiltonian*

$$H_c(t)=\underset{j=1}{\overset{m}{}}H_ju_j(\theta ,t),$$

where $`H_j^{n\times n}`$ are also Hermitian matrices, accounts for the effects of the control inputs $`u_1(\theta ,t),\mathrm{},u_m(\theta ,t)`$ on the dynamics of the system. We assume that these control functions depend on a finite number of real parameters $`\theta =(\theta _1,\mathrm{},\theta _p),\theta _k𝒯`$, with $`𝒯`$ being an open set in $`^p`$. This kind of assumption is reasonable if we think of the small set of parameters we can control in an experimental setting.

We consider the problem of steering the system from a given initial state $`|\psi _0=|\psi (0)`$ to a final state $`|\psi _1`$, where $`|\psi _0`$ and $`|\psi _1`$ are unit vectors in $`^n`$. We assume that the transition occurs (at t=T) when $`\theta =\theta ^{}`$ which we take as “nominal” value of the parameters. Clearly, if $`\theta \theta ^{}`$, the transition will, in general, not occur.

It is convenient to introduce the *error probability* for each control strategy. Consider the normalized final state for the time evolution, $`|\psi (T,\theta )`$. It can be written as $`|\psi (T,\theta )=\psi _1|\psi (T,\theta )|\psi _1+|\psi ^{}(\theta ,T)`$ with $`|\psi ^{}(\theta ,T)`$ orthogonal to $`|\psi _1`$. If we imagine to perform a discrete measure111Quantum measure fundamental postulates can be founded in standard quantum mechanics textbooks, see e.g. , or . on an observable that has $`|\psi _1`$ as eigenstate, the probability to obtain the eigenvalue associated $`|\psi _1`$(that corresponds to the probability of finding the system in $`|\psi _1`$ immediately after the measure) is: $`P_{|\psi _1}=|\psi _1|\psi (T,\theta )|^2.`$ Then the *error probability* corresponding to the value $`\theta `$ is:

$`P_{err}(T,\theta )`$ $`=`$ $`1|\psi _1|\psi (T,\theta )|^2`$ (3)
$`=`$ $`\psi ^{}(\theta ,T)|\psi ^{}(\theta ,T),`$
thanks to the fact that $`|\psi (T,\theta )`$ is normalized. By assumption, we have $`P_{err}(T,\theta ^{})=0`$.

## 2 Robustness of the control strategy

In the quantum control field, the expression “robustness of the control strategy” means that the control performance is insensible with respect to errors in the control implementation. In , a control strategy is considered robust “*if significant local changes in the amplitude and the form of the pulse and of the chirp do not change significantly the final transfer probability*.” The pulse and the chirp, in the setting described in Section 3, are the system inputs parameters. A quantitative definition of robustness is, however, missing.

The contribution of the present paper is to provide such a quantitative definition reformulating the problem as a robust control-theoretic problem, and then to analyze the robustness properties of some significant strategies considered in the relevant literature. In terms of our model, this concept of robustness can be qualitatively formulated as follows: A control strategy is robust when, for values of the parameters $`\theta `$ different from the nominal ones, the final state $`|\psi (T,\theta )`$ is close to the desired one $`|\psi _1`$. This robustness request is satisfied if $`P_{err}(T,\theta )`$ is small in the parameter set $`𝒯`$.

In classical control theory, plant uncertainty is described by a set $`𝒫`$ of possible plants . This uncertainty can be either structured (parametrized by a finite number of scalar parameters or a discrete set of plants) or unstructured (disc-like uncertainty). A controller is said to be robust with respect to some property if this property holds for every plant in $`𝒫`$. It is quite simple to reformulate our problem as a particular case of structured-like classical robustness problem. First, we notice that our quantum “plant” is determined by the matrices $`(H_0,H_1,\mathrm{},H_m)`$. These matrices determine the system Hamiltonian (2), given the control strategy. Let $`P_0=(H_0,H_1,\mathrm{},H_m)`$ be our nominal plant.

As in , we can transfer the uncertainty from the control parameters to the internal Hamiltonian. In fact, by defining $`\delta u_i(\theta )=u_i(\theta )u_i(\theta ^{})`$ we can write:

$`H(t)`$ $`=`$ $`H_0+{\displaystyle \underset{i=1}{\overset{m}{}}}H_i\left(u_i(\theta ^{})+\delta u_i(\theta )\right)`$ (4)
$`=`$ $`H_0+{\displaystyle \underset{i=1}{\overset{m}{}}}H_i\delta u_i(\theta )+{\displaystyle \underset{i=1}{\overset{m}{}}}H_iu_i(\theta ^{})`$
$`=`$ $`\left(H_0+\mathrm{\Delta }H_u(\theta )\right)+{\displaystyle \underset{i=1}{\overset{m}{}}}H_iu_i(\theta ^{}).`$
where $`\mathrm{\Delta }H_u(\theta )=_{i=1}^NH_i\delta u_i(\theta )`$. Such a cosmetic transformation shows that our control strategy uncertainty can be seen as a particular case of the plant uncertainty (with control inputs $`u_i(\theta ^{})`$). The plant set $`𝒫`$ is here given by:

$$𝒫=\{(H_0+\mathrm{\Delta }H_u(\theta ),H_1,\mathrm{},H_m)|\theta 𝒯\}.$$

The property we are interested in, as e.g. in , is the error probability. We require this probability not to exceed a fixed threshold $`ϵ[0,1)`$ at a given $`T`$. All the ingredients of a classical robustness problem have now been specified. Introduce the *$`ϵ`$-robustness set* $`_ϵ`$ as

$$_ϵ=\{\theta 𝒯|P_{err}(\theta ,T)ϵ\}.$$
(5)
We give the following definition.

###### Definition 1.

A control strategy $`\{u_1(\theta ^{},t),\mathrm{},u_m(\theta ^{},t)\},t[0,T]`$ is *$`ϵ`$-robust* with respect to parameters uncertainty if:

$$_ϵ=𝒯.$$
(6)
Notice that only the 0-robustness case ensures us an exact steering of the system state for all $`\theta 𝒯`$.

## 3 Some applications

In this section we analyze, in the light of the above definition, the robustness of some prototype examples. We will compare our results with qualitative observations and robustness claims in the relevant literature. To do so, we introduce a particular form of (1) frequently used in NMR (Nuclear Magnetic Resonance) quantum control problems. We consider a two level quantum system, and the associated a bi-dimensional Hilbert space. The time evolution is described by a scaled time Schrödinger equation in the form:

$$i\mathrm{}\frac{}{s}|\psi (s)=TH(s)|\psi (s),$$
(7)
where $`s=t/T`$ and

$$H(s)=\left(\begin{array}{cc}\mathrm{\Delta }(s)& \mathrm{\Omega }(s)\\ \mathrm{\Omega }(s)& \mathrm{\Delta }(s)\end{array}\right),$$

is represented in the canonical (*diabatic*) base. The control functions:

$$\{\begin{array}{c}\mathrm{\Delta }(s)=\mathrm{\Delta }_0\mathrm{\Phi }(s)\hfill \\ \mathrm{\Omega }(s)=\mathrm{\Omega }_0\mathrm{\Lambda }(s)\hfill \end{array},$$

are the inputs, with $`\mathrm{\Phi }(s)`$, $`\mathrm{\Lambda }(s)`$ fixed envelops and $`\mathrm{\Delta }_0`$, $`\mathrm{\Omega }_0^+`$ amplitude parameters. In this picture we have $`\mathrm{\Delta }(s)=u_1(s,\mathrm{\Delta }_0)`$ and $`\mathrm{\Omega }(s)=u_2(s,\mathrm{\Omega }_0)`$. Thus $`\theta =(\theta _1,\theta _2)=(\mathrm{\Delta }_0,\mathrm{\Omega }_0)`$ are the parameters we are interested in. In the contest of particle-laser field interaction and the RWA (Rotating Wave Approximation ), these functions depend on the chirp (*detuning*) and the amplitude (*time-dependent Rabi frequency*) of the active pulse222To find some detailed information about the physical meaning of these parameters and about the resonance phenomenon see ,.. This model can be seen as a particular case of model (1), and it is suitable to describe control techniques based both on magnetic resonance and adiabatic passage. At each time the unitary transformation

$$U(s)=\left(\begin{array}{cc}\mathrm{cos}\theta (s)& \mathrm{sin}\theta (s)\\ \mathrm{sin}\theta (s)& \mathrm{cos}\theta (s)\end{array}\right),$$
(8)
with $`\theta (s)=\frac{1}{2}\mathrm{arctan}\left(\mathrm{\Omega }(s)/\mathrm{\Delta }(s)\right)`$, diagonalizes the Hamiltonian

$$U^{}(s)H(s)U(s)=\left(\begin{array}{cc}\epsilon (s)& 0\\ 0& \epsilon (s)\end{array}\right)=D(s).$$
(9)
Here

$$\epsilon (s)=\sqrt{\mathrm{\Delta }^2(s)+\mathrm{\Omega }^2(s)}.$$
(10)
is the energy eigenvalue. Applying $`U(s)`$ as a time dependent change of basis and defining $`\varphi `$, we obtain:

$`i\mathrm{}{\displaystyle \frac{|\varphi (s)}{s}}`$ $`=`$ $`\left[TD(s)iU^{}(s){\displaystyle \frac{}{s}}U(s)\right]|\varphi (s)`$ (13)
$`=`$ $`\left(\begin{array}{cc}T\epsilon (s)& i\gamma (s)\\ i\gamma (s)& T\epsilon (s)\end{array}\right)|\varphi (s).`$
The new basis vectors are called *adiabatic* states. In the adiabatic limit, $`T\mathrm{}`$, the $`\gamma (s)`$ terms can be neglected, as shown by the adiabatic approximation theory.

In the following subsections, we will investigate the robustness properties of different control strategies in a typical steering problem, the state-flipping: Transfer the system state from one basis vector to the other. The standard resonance technique and two adiabatic models will be discussed and compared.

### 3.1 Magnetic Resonance

A simple way to obtain such a transfer is by using the magnetic resonance phenomena: Under properly tailored oscillating fields, the state vectors rotate between the two basis states ,. This kind of effect can be generated by the following fields-control functions:

$$\{\begin{array}{c}\mathrm{\Delta }(s)=0\hfill \\ \mathrm{\Omega }(s)=\mathrm{\Omega }_0\mathrm{\Lambda }(s),\hfill \end{array}$$
(14)
where $`\mathrm{\Lambda }(s)`$ is the $`\mathrm{\Omega }`$-pulse envelope. This parametrization, and some easy calculations , lead to the following expression for the error probability:

$$P_{err}(T,\mathrm{\Omega }_0,A_\mathrm{\Lambda })=\mathrm{cos}^2\left[T\mathrm{\Omega }_0_{s_i}^{s_f}\mathrm{\Lambda }(s)𝑑s\right]=\mathrm{cos}^2\mathrm{\Omega }_0TA_\mathrm{\Lambda },$$
(15)
with $`A_\mathrm{\Lambda }`$ the $`\mathrm{\Omega }`$-pulse area.

This probability is equal to zero for:

$$\mathrm{\Omega }_{0,k}^{}=\frac{(k+\frac{1}{2})\pi }{TA_\mathrm{\Lambda }},k=0,1,2,\mathrm{},$$

or, equivalently, for:

$$A_\mathrm{\Lambda }^{}=\frac{(k+\frac{1}{2})\pi }{T\mathrm{\Omega }_0},k=0,1,2,\mathrm{}.$$

Thus, $`P_{err}(T,\mathrm{\Omega }_0,A_\mathrm{\Lambda })=0`$ in a family of hyperbolas parameterized in $`k`$ (zero-measure set in the parameters space). Consider $`𝒯=[\mathrm{\Omega }_0^{}\beta ,\mathrm{\Omega }_0^{}+\beta ]\times [A_\mathrm{\Lambda }^{}\sigma ,A_\mathrm{\Lambda }^{}+\sigma ]`$, a common setting if we are working with nominal values and error intervals, with $`\beta ,\sigma `$ such that $`\mathrm{sin}2T\mathrm{\Omega }_0A_\mathrm{\Lambda }`$ is monotone in every direction. Then, the maximum absolute value for the error probability in $`𝒯`$ is:

$$P_{max}=\mathrm{cos}^2T\overline{\mathrm{\Omega }}\overline{A}_\mathrm{\Lambda },$$
(16)
where $`\overline{A}_\mathrm{\Lambda }=A_\mathrm{\Lambda }^{}+\sigma `$ and $`\overline{\mathrm{\Omega }}=\mathrm{\Omega }_0^{}+\beta `$. Then $`𝒯`$ is $`\overline{ϵ}`$-robust, with $`\overline{ϵ}=P_{max}`$.

According to qualitative evaluation found, the magnetic resonance strategy doesn’t seem to ensure enough insensibility towards errors in control implementation and can be sensitively improved by adiabatic passage techniques.

### 3.2 Landau-Zener Model

The Landau-Zener model is one of the simplest choice of control function leading to an adiabatic transition. We will consider:

$$\{\begin{array}{c}\mathrm{\Delta }(s)=\frac{\mathrm{\Delta }_0^2}{T}s\hfill \\ \mathrm{\Omega }(s)=\mathrm{\Omega }_0\hfill \end{array},$$
(17)
the *detuning* varies linearly with a zero crossing, while the Rabi frequency is maintained constant. For $`s=0`$ we have a minimum in the difference between energy levels that leads to a state inversion if the evolution satisfies the condition needed for the adiabatic approximation. The error probability is estimated with the Landau-Zener formula:

$$P_{err}(T,\mathrm{\Omega }_0,\mathrm{\Delta }_0)e^{\pi T\frac{\mathrm{\Omega }_0^2}{\mathrm{\Delta }_0^2}}.$$
(18)
This probability goes to zero in the adiabatic limit $`T\mathrm{}`$ for any choice of $`\mathrm{\Omega }0`$ and . Thus, the robustness set for this strategy is the whole open first quadrant without its boundary ($`\mathrm{\Omega }_0=0,\mathrm{\Delta }_0=0`$). The advantages given by this adiabatic technique are evident, as long (18) estimates correctly the error probability. Even if the $`T\mathrm{}`$ condition is not realizable, we can take a $`T`$ large enough to maintain $`P_{err}`$ arbitrary small for (almost-)every parametrization of the control strategy. We will call this behavior *intrinsically robust*.

### 3.3 Allen-Eberly Model

We now analyze an adiabatic control strategy more complex than the previous one. The Allen-Eberly  parametrization allows to obtain an exact expression for the error probability and, in the $`\mathrm{\Omega }_0=\mathrm{\Delta }_0`$ case, forces the state time evolution along the energy *level lines*, maintaining $`\epsilon (\mathrm{\Omega }(s),\mathrm{\Delta }(s))=c`$, $`c`$ constant for every $`s`$ . This kind of choice leads to good results in terms of error probability even quite far from the ideal $`T\mathrm{}`$ condition, as we are going to show. In terms of control functions, we consider:

$$\{\begin{array}{c}\mathrm{\Delta }(s)=\mathrm{\Delta }_0\sqrt{1\text{sech}^2(s)}=\mathrm{\Delta }_0\mathrm{tanh}(s)\hfill \\ \mathrm{\Omega }(s)=\mathrm{\Omega }_0\text{sech}(s).\hfill \end{array}$$
(19)
Then, the exact expression for the error probability is:

$$P_{err}(T,\mathrm{\Omega }_0,\mathrm{\Delta }_0)=\mathrm{cosh}^2\left(\pi T\sqrt{\mathrm{\Delta }_0^2\mathrm{\Omega }_0^2}\right)\text{sech}^2\left(\pi \mathrm{\Delta }_0T\right),$$
(20)
for every regime, adiabatic or not. We can notice that, for large $`T`$ and for $`\mathrm{\Delta }_0\mathrm{\Omega }_0`$, the error probability can be bounded by:

$$P_{err}(T,\mathrm{\Omega }_0,\mathrm{\Delta }_0)4e^{2\pi T(\mathrm{\Delta }_0\sqrt{\mathrm{\Delta }_0^2\mathrm{\Omega }_0^2})}.$$
(21)
Thus, for every $`\mathrm{\Delta }_0`$ and $`\mathrm{\Omega }_0`$, $`\mathrm{\Delta }_0\mathrm{\Omega }_0`$, the error probability decreases exponentially to zero in the adiabatic limit. The best choice for the parameters values is to take the largest $`\mathrm{\Delta }_0=\mathrm{\Omega }_0`$. In the case $`\mathrm{\Omega }_0>\mathrm{\Delta }_0`$, the error probability becomes:

$$P_{err}(T,\mathrm{\Omega }_0,\mathrm{\Delta }_0)=\mathrm{cos}^2\left(\pi T\sqrt{\mathrm{\Omega }_0^2\mathrm{\Delta }_0^2}\right)\text{sech}^2\left(\pi \mathrm{\Delta }_0T\right).$$
(22)
This expression tends to zero with dumped oscillations, due to the term $`\mathrm{cos}^2\left(\pi T\sqrt{\mathrm{\Omega }_0^2\mathrm{\Delta }_0^2}\right)`$. Again, larger $`\mathrm{\Delta }_0`$ make $`P_{err}`$ converge faster. Thus, for each fixed $`ϵ`$, we can compute a $`T_ϵ`$ such that the error probability $`P_{err}(T,\mathrm{\Omega }_0,\mathrm{\Delta }_0)<ϵ`$ for every $`T>T_ϵ`$. Indeed, it is easy to see that

$$T_ϵ=\mathrm{max}\{\frac{\mathrm{ln}\frac{ϵ}{4}}{2\pi (\mathrm{\Delta }_0\sqrt{\mathrm{\Delta }_0^2\mathrm{\Omega }_0^2})},\frac{\mathrm{ln}ϵ}{2\pi \mathrm{\Delta }_0}\}.$$

This control strategy is therefore *intrinsically robust* for $`T`$ sufficiently large. According to the Landau-Zener case, every choice of $`\mathrm{\Omega }_00`$ and $`\mathrm{\Delta }_00`$ drives the system to the target state. The level line condition $`(\mathrm{\Delta }_0=\mathrm{\Omega }_0)`$ and large $`\mathrm{\Omega }_0`$ give faster convergence to the desired state.

## 4 Discussion

Comparing the results, the advantages given by the adiabatic strategies are evident. They can be effectively used, however, when the transfer time is not critical: Their intrinsic robustness is exhibited only with a large time.

The examples analyzed are also treated in  to illustrate that control strategies based on the level lines are optimal for adiabatic population transfer (the level line strategies minimize the error probability and corresponds to the minimum pulse area). In , robustness of the control is also taken in account: Contour plots of error probability with respect to parameters variations are obtained thanks to numerical simulations for the system evolution. Different strategies are qualitatively compared. It is shown that the simple resonance case generates larger error probability than the adiabatic optimal techniques, once a parameter variation is fixed.

Here we have given a formal definition of the robustness property, reformulating the problem in the control theoretical setting. We have obtained quantitative results consistent to the qualitative ones just described, and we have provided an analysis tool useful to evaluate and compare robustness behavior of different strategies. From a control theoretic viewpoint, we have analyzed a specific robustness problem for open-loop control of a bilinear system.