{'fraction_non_alphanumeric': 0.033501481303547125, 'fraction_numerical': 0.010441188245656178, 'mean_word_length': 4.636429280620995, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The stochastic Gross-Pitaevskii equation

We show how to adapt the ideas of local energy and momentum conservation in order to derive modifications to the Gross-Pitaevskii equation which can be used phenomenologically to describe irreversible effects in a Bose-Einstein condensate. Our approach involves the derivation of a simplified quantum kinetic theory, in which all processes are treated locally. It is shown that this kinetic theory can then be transformed into a number of phase-space representations, of which the Wigner function description, although approximate, is shown to be the most advantageous. In this description, the quantum kinetic master equation takes the form of a GPE with noise and damping added according to a well-defined prescription--an equation we call the stochastic GPE. From this, a very simplified description we call the phenomenological growth equation can be derived. We use this equation to study i) the nucleation and growth of vortex lattices, and ii) nonlinear losses in a hydrogen condensate, which it is shown can lead to a curious instability phenomenon.


Introduction
There is a growing consensus that in some sense the dynamics of a trapped Bose-Einstein condensate can be properly described by the Gross-Pitaevskii equation, provided appropriate statistical assumptions are introduced. This was first suggested by Svistunov [1,2,3], and numerical experiments were carried out by Damle and Sachdev [4], Marshall, New and Burnett [5], and most recently by Davis et al [6]. All of these considered the Gross-Pitaevskii equation for a Bose gas inside a box, with random initial conditions, such that the total energy and number of particles would correspond to those of a partially condensed system in thermal equilibrium. The numerical experiments then verified that this initial ensemble developed with time into a thermodynamic equilibrium ensemble as a result of the time evolution induced on each member of the ensemble by the Gross-Pitaevskii equation. The recent work of Davis et al [6] also checked the dynamics of the resulting equilibrium state against Morgan's [7] theoretical predictions from quantum field theory, and found that these were in agreement with each other.
The work of Steel et al. [8] took a different, and more rigorous view, that the correct physical relationship between the condensate wavefunction with random initial conditions and the physical field operator is made by means of a truncated Wigner function representation. It is known [9] that for the standard Bose-Einstein condensate Hamiltonian, the Wigner function obeys a generalized Fokker-Planck equation, in which derivatives up to third order occur. However, if the amplitudes of the Wigner function phase-space variables are large, the third order derivative terms are small, and can be neglected-this "truncation" gives the method its name. However, we also find that the expected second-order derivative terms, which correspond to noise, are for this Hamiltonian exactly zero, and thus what remains corresponds to a Liouville equation for an ensemble of wavefunctions which obey the Gross-Pitaevskii equation.
However the description can only be valid provided the amplitudes can be regarded as large, and this must be the case for all modes. Clearly, the higher energy (or smaller wavelength) modes will always have a very small amplitude, and cannot be described by a truncated Wigner function. Indeed, if no truncation is done, then the result of the quantum mechanical constraints-essentially Heisenberg's uncertainty principle-on the Wigner function leads to a delta function spatial correlation function, corresponding to infinite fluctuations at each point in space, which cannot be simulated. In all simulations, however, the implementation of the spatial differentiation requires a spatial grid, whose spacing corresponds to the inverse of the highest spatial frequency used, thus providing a cutoff in both space and energy. The results must then be cutoff dependent.
The work of Davis et al [6] is the first to put down explicitly that a cutoff in energy or spatial frequency must be imposed on this Gross-Pitaevskii equation equation. However, the basic idea, that high and low energies require different treatments is that which forms the basis of our quantum kinetic theory, which gives a formalism for coupling together a fully thermalized noncondensate band, consisting of particle with energies more than a certain value E R , to a condensate band, consisting of all lower energy particles. The non-condensate band can be treated by the quantum Boltzmann equation, while the condensate band should in principle be treated fully quantum mechanically.
In this paper we will combine the Wigner function ideas of Steel et al. [8] with quantum kinetic theory, and as well introduce the idea of local energy conservation as suggested by Zaremba et al. [10,11] to derive a method of coupling the thermalized noncondensate band to the Gross-Pitaevskii equation. This results in a Gross-Pitaevskii equation with added noise-a stochastic Gross-Pitaevskii equation.
The paper consists of three main parts: Sect.2 implements the idea of local energy and momentum conservation in quantum kinetic theory and results in a relatively simple master equation for the interaction of a condensate with a fixed bath of non-condensed atoms; Sect. 3 shows how to implement the Wigner function method to transform the master equation into a c-number equation with random initial conditions, and to compare with P-and Q-function methods; Sect.4 develops a very much simplified phenomenological equation, which we show can be applied to the two-body loss problem in the hydrogen condensate experiments, and to the stabilization of quantized vortex arrays. In Sect.6.1 we adapt the phenomenological growth equation to include losses from dipolar relaxation in a hydrogen condensate, and show that the nonlinear losses so introduced could lead to a "boom and bust" instability in which the condensate grows and collapses repeatedly, although the conditions under which hydrogen condensates are presently formed make it difficult to say whether this could be observed in practice.

Application of local energy conservation to quantum kinetic theory
In a formulation of the kinetic theory of non-condensed vapour in interaction with a condensate, Zaremba, Nikuni and Griffin [11] have used the concepts of local energy and momentum conservation to develop appropriate equations of motion for the system in a formulation based on Hartree-Fock-Popov methodology [12,13,14,15]. This is in contrast to our own [16,17,18,19,20,21,22,23] formulation of quantum kinetic theory, in which energy conservation is expressed in terms of transitions between eigenstates of the condensate. While it is clear that a description in terms of eigenfunctions must be more accurate, it yields equations which are not easy to handle exactly, and thus the greater accuracy can in practice be an illusion.
The full formulation of our quantum kinetic theory is to be found in QKV ‡, which depends strongly on QKIII. The major improvement in QKV is a full consideration of the mean field effects of condensate on vapour and conversely, but the general methodology is very similar. In these we divide the excitation spectrum at an energy E R , above which the excitations can be treated as particle like. The field operator is written (in the Schrodinger picture) as where the first operator represents the modes with excitation energies above E R (and is called the noncondensate band field operator), while φ(x) represents the excitations with energies below E R (and is called the condensate band field operator). This division is chosen and fixed for any given problem, and is to some extent arbitrary, since it is expected that the higher modes of the condensate band will be very close to particle-like, and in practice will also be thermalized.
Let us look at the derivation of energy conservation used in QKIII, using the notation used there. On p539 we examine a term in the master equation involving H (1) I,C as defined in (QKIII.8), which contains one φ(x) or φ † (x), where φ(x) is the field operator for the condensate band. Explicitly, the term in Hamiltonian can be written In this equation we have defined a notation Substituting into the master equation (QKIII.20), terms arise which are of the form Here we have used the notation In QKIII we expanded the operators φ(x) in eigenoperators of the condensate band Hamiltonian H 0 , so as to be able to perform the integral over the time τ , and thus arrive at the final master equation, which would not involve φ operators at different times. This method does not require that there be a condensate. If we assume that there is a condensate present, we use this knowledge to make an estimate of the time development in the time τ .

Local energy conservation
The basis for this concept is given by the works of Zaremba, Nikuni, Griffin, and co-workers [10,11,26], which themselves draw on the work of Kirkpatrick and Dorfman [27,28,29], and can be seen phenomenologically by use of the density-phase description of the Gross-Pitaevskii equation By writing now the Gross-Pitaevskii equation takes the form in which the condensate velocity v C (x, t), and the local chemical potential µ C (x, t) are given by This transcription of the Gross-Pitaevskii equation is the appropriate basis for a description of condensate behaviour when the condensate density n C (x, t) is slowly varying in space and time, so that we may consider the dependence on space and time to arise largely from the phase Θ(x, t). This leads to the hydrodynamic approximation in which the Laplacian term in (13) is dropped and thence to the energy conservation equation The static solution of the hydrodynamic description is the Thomas-Fermi description of the condensate wavefunction, in which v C (x, t) = 0, In the hydrodynamic regime, the local chemical potential µ C (x, t) is a strictly local quantity, and is a simple function of the density n C (x, t). The local energy density ǫ C (x, t), defined in (11) is explicitly equal to the derivative of the energy density (also evaluated using the hydrodynamic approximation) with respect to n C (x, t). Thus, when considering the possible addition of a particle to the condensate, it is intuitively appealing to say that this takes place at a position x, and that the energy added to the condensate by this process is ǫ(x, t). This is given some further backing by the knowledge that energy conservation in quantum mechanics depends on frequency matching, and that to the extent that n C (x, t) and ǫ(x, t) depend slowly on time, the relevant frequency will be ǫ(x, t). If n C (x, t) and v C (x, t) are both slowly dependent on space, the wavelength of the wavefunction will correspond to a momentum mv C (x, t), and thus momentum conservation will involve this momentum.

Application to quantum kinetic theory
The major task in developing a quantum stochastic description is to find a simple way of expressing the term φ(x ′ , −τ ) in terms of φ(x, 0), and we would like to follow Zaremba et al [11] in doing this. We assume that the evolution over this short time can be approximated by simply the phase evolution of the condensate according to (10,11). We also want to express the result in terms of operators at the same location, which we take as the midpoint u = (x + x ′ )/2. Thus, using the symbol y = x − x ′ , we write There are two main assumption here: i) There is an assumption that the state of the condensate band is dominated by a single condensate wavefunction. However, the wavefunction may have any form, and of course may also be time dependent. It is possible that a more accurate method could be developed using hydrodynamic quantization, such as that developed by Marques and Bagnato [30], which yields operator equations almost the same as the standard hydrodynamic formalism outlined above.
ii) The approximation (21) requires that only small y and τ contribute. This will be valid when the function Tr (4), is almost local in space and τ . This kind of behaviour is expected if the noncondensate band is fully thermalized at a sufficiently high temperature.
Both of these assumptions are used by Zaremba et al [11]

Final form of the master equation
By carrying out the procedures used in QKIII, we can finally get the master equation: For compactness, we have defined Here F (K, x) is the phase space density of the noncondensate band in terms of wavenumber K and position x. The transition rates E (±) , M , G (±) are defined by In the above equations we use the notations 2.3.1. Effect of the cutoff at E R The fact that φ(x) has an upper energy cutoff means that the commutator [φ(x), φ † (x ′ )] is only an approximate delta function. This has no influence on the form of the master equation, but equations of motion derived from it for averages of φ will automatically ensure that solutions remain in the correct subspace.

Conservation of energy and momentum
The choice made in Eqs (20,21) is very drastic, assigning as it does a single phase to all condensate band modes, and this has the consequence that the terms (25)(26)(27) are in fact zero. For the last two this is obvious, since the energy and momentum conservation delta functions in the definitions of E (±) would require the equality of the energy and momentum of a particle in the condensate band with those of a particle in the noncondensate band, which is by definition not possible. The vanishing of the term (25) is not so immediately obvious, but momentum and energy conservation here mean that in fact this is a forward scattering term, which the methodology of QKIII and QKV explicitly removes from the irreversible part. The origin of this is the requirement that, at a given position x, there is a unique energy and momentum which a particle in the condensate band can have. This is actually only true for particles contained in the condensate itself, so that we can expect that there is a spread of values of ǫ C (x, t) and v C (x, t). As well as this, it should be borne in mind that the hydrodynamic approximation is the foundation of these results, and this is indeed an approximation. Inclusion of corrections to this approximation would also yield a spread in local momenta and energies available for particles in the condensate band.
The effect of including a spread in values of ǫ C (x, t) and v C (x, t) would be principally to give a nonzero value to the terms (25)(26)(27) in the master equation. The last two terms would certainly be very small, but the scattering term (25) should become significant. The growth terms (22,23) would be modified, but not greatly. However the changes are only in the nonoperator coefficients-the basic structure of the master equation remains the same. The major point of simplification compared with the treatment of QKIII and QKV is the expression of the irreversible terms directly in terms of the field operators φ(x), which enables a much simpler analysis to be carried out.

Time development equations for the noncondensate band
The noncondensate band is described by the phase space density F (K, x), whose equation of motion is given in QKV, in the form of a quantum Boltzmann equations and with added terms to take account of transfer between the condensate and noncondensate bands. For simplicity, in this paper we assume that the noncondensate band is fully thermalized with temperature T and chemical potential µ.

Stochastic differential equations for the Wigner function
A master equation of the form we have produced can be treated using phase space representation methods [9]. These will give rise to equivalent stochastic differential equations for a c-number phase space variable α(x, t) whose averages are related to those of the operator φ(x) using symmetric ordering: thus Using the operator correspondences in [9] Eq. (4.5.12), we can deduce the stochastic differential equation Here the last two terms are Gaussian white noise terms. The quantity dW M (x, t) is a real Wiener noise term, to be interpreted in the Ito sense [31], whose correlation functions are while the other noise is complex, with correlation functions

Neglect of third order noise terms
The equation is also approximate in that third order noise terms do arise along with the trilinear term. These terms arise as third order partial derivatives in the Fokker-Planck equation generated by the Wigner function representation of the master equation (22)(23)(24)(25)(26)(27), and have no straightforward stochastic interpretation. When the variable α(x) is large, corresponding to large phase space occupation, these terms become small, as pointed out in [8]. Stoof [32] has developed a similar form of Wigner function stochastic Gross-Pitaevskii equation using path integral methods, and this has been used to analyse damping of condensate modes and the reversible formation of a condensate [33,34]. The actual form of the damping is not the same as ours, but it is easy to see that to lowest order in the damping it is equivalent, and this is all we would claim for our derivation. The noise terms are equivalent to ours.

Stoof's equation
In [33] a case is considered where one can take a "classical" limit of this Wigner approach-a low energy approximation and essentially always valid for the description of the dynamics of the condensate and its low-lying excitations-and in this approximation the P-Q-and Wigner function forms of the noise become the same. This leads to an approximate noise term which is of the the same form as that of the Q-function.

Use of the P-function
It is more conventional in quantum optics to use either the Glauber P-function or the positive P-function, instead of the Wigner function. The main advantage is that the approximation, necessary for the Wigner function, that third order noise terms must be neglected is not required-the resulting stochastic differential equations are exact.
For either P-function form, we would replace the symmetrized product averaging rule with one involving normal products, so that we would obtain (52) Using the operator correspondences in [31] Eq. (4.5.9), we can deduce the same stochastic differential equation (43), but with the correlation functions (The same as for the Wigner function), while the other noise has the altered correlation functions (58) The interpretation of (43) as a genuine stochastic differential equation requires that the matrix of noise coefficients should have only nonnegative eigenvalues. For higher temperature situations in which there is a substantial thermal component, this will certainly be true for all values of the variable α(x) which would turn up in a stochastic simulation. When this is not so, a Positive Prepresentation would be necessary. The experience of Drummond and co-workers [35,36] has shown this is in principle feasible, but application to experimentally realistic problems would be very difficult.

Use of the Q-function
In the case of the Q-function one should replace the symmetrized product averaging rule with one involving antinormal products, so that we would obtain Using the operator correspondences in [9] Eq. (4.5.10), we can deduce again the same stochastic differential equation (43), but with the correlation functions (Again the same as for the Wigner function), while the other noise has the altered correlation functions The condition for this to be a genuine stochastic differential equation is that the matrix Thus, although the Q-function always exists and is positive, this does not necessarily mean the stochastic differential equation is valid-see [31] Ch.6.

Approximations and simplifications
3.3.1. Effect of the cutoff These equations are expressed in a simplified form, since by definition the condensate band contains a finite range of energies, and hence of wavelengths. This means that the Laplacian, the delta functions, and the form of the trilinear term are modified by a projection into this band, as explained in QKIII and QKV. As well, the delta functions become delocalized, and the noise correlation functions are then well defined at x = x ′ . The actual implementation of this procedure is somewhat technical, but is essentially straightforward.
A reasonable way to do this is to express the δ(x − y) terms as a summation over a set of orthogonal condensate eigenfunctions. This amounts to truncating an exact expression for a delta function at the cutoff energy E R . But for sufficiently high energies the trap eigenfunctions approach the harmonic oscillator eigenfunctions, so we may instead use the expression in terms of harmonic oscillator eigenfunctions y * This representation necessarily gives an added noise which is only significant for x, y inside the region classically allowed for energy equal to E R , and thus permits the use of fast Fourier transform methods.

Effect of the cutoff on means and variances
In application of the Wigner function to a field theory, the fact that each mode contributes an extra half of a quantum to averages-as shown in (42)-means that there will be a contribution proportional to the number of modes retained below the cutoff. For the Q-function the effect is similar, but the contribution is a whole quantum. Suppose M modes contribute to the summation in (71). If we define the particle number operators for the field theory and the various phase space representationŝ then the means and variances are related by These results directly present the dilemma one is faced with in making a choice of representation. If we wish to consider only positive phase-space distributions, which therefore have a probabilistic interpretation, we can see immediately from these results that a P-function is only possible if there is at least a Poissonian number distribution-that is, the P-function does not always exist. (However the Positive P-representation can evade this restriction, at the cost of doubling the number of independent variables.) In contrast, there is no restriction on the Wigner and Q-functions, which are known to exist for all physical states. On the other hand, only the P-function has no "vacuum noise" contribution, proportional to M , the number of modes, unlike the Wigner and Q-functions. For these, the existence of the "vacuum noise" contribution means that a simple interpretation of the stochastic differential equation (43) as the equation of motion for the condensate wavefunction must require that the actual mean number of real particles N not be swamped by the vacuum contribution; that is we must require N ≫ M at the very least. This means that the mean occupation per mode is very much greater than 1, and that a Bose-Einstein condensate must already be present. For the P-function this is not a problem-M can be as large as we wish-even infinite if we desire. Finally, one must emphasize that the stochastic differential equation (43) is approximate for the Wigner function interpretation, even with the appropriate noise properties, because of the neglect of third-order noise terms, as previously noted. Fortunately, the condition that we van neglect these terms is the same as the condition that the vacuum modes not be dominant-that the occupation per mode is high for the modes of interest. On the other hand provided we have an initial distribution which is broader than the Poisson, and the temperature is sufficiently high for for the positivity condition of the noise matrix (59) to be satisfied, the P-function interpretation of the stochastic differential equations (43) is valid and exact.

The low temperature case
For sufficiently low temperatures it is clear that all the transition rates (31)(32)(33)(34)(35) are negligible, but the noise matrices (59,69) no longer have positive eigenvalues. This leaves us with only two feasible choices-either a Positive P-function interpretation, which is numerically very difficult, or a Wigner function interpretation, which is approximate, but has the attractive property that the noise contributions all vanish, as can be seen from (44)(45)(46)(47)(48)(49). This choice has been extensively investigated by Sinatra et al. [37].

The stochastic Gross-Pitaevskii equation for a pure condensate
To illustrate the procedure, we omit the noise and damping terms in (43), and consider the resulting time dependent Gross-Pitaevskii for a random function α(x, t). This takes the form i) α(x, t) satisfies the Gross-Pitaevskii equation: The prescription of statistics in terms of symmetrically ordered products as in (40)(41)(42) means that the initial conditions must be random even in the case of an initial pure state.
There are two possibilities, associated with the fact that the total number of particles is now exactly conserved, because of the vanishing of the terms (31)(32)(33)(34)(35). A straightforward interpretation would be to assume that the initial wavefunction of the condensate mode is ξ(x, 0) and use random initial conditions corresponding to all N particles being in the condensate, and none in the other M − 1 modes. Unfortunately, the Wigner function for an N particle state involves the Laguerre polynomial: and for N = 0 this is a highly oscillatory function, as shown in Fig. 1, which has no probability interpretation. However, the cumulative integral behaves very like a step-function, and this means that the mean values of any smooth functions of α, such as the first few powers, will be very well approximated by using a Gaussian approximation of the form α → √ n e iθ (85) where θ has a uniform distribution on (0, 2π) and n is a real Gaussian random variable with (87) It should be borne in mind that in practice we are quite unable to even contemplate measuring any moments of α(x) higher than the fourth, and this form does correctly give all moments up to and including the fourth order.
ii) The remaining modes would contain no particles, and thus for these the Wigner function (83) is to be evaluated with N = 0, which is then is positive and Gaussian. iii) This means that we set Here α k are complex Gaussian random variables, independent of each other, such that It follows that that ϕ(x, 0) is a Gaussian random function with statistics given by iv) The statistics of the initial random condensate wavefunction are Gaussian, with the nonzero means and correlations given by α * (x, 0)α(y, 0) = N ξ * (x, 0)ξ(y, 0) v) If there are M wavefunctions y n included in the summations (71), and hence M −1 terms in the summation over k in (88), then the mean number of particles N W corresponding to the ensemble in (97,98) is given by This can be seen correspond to N real particles, plus the effect of "half a quantum" per mode for the M modes included.
vi) If one runs a simulation for a time, and wishes to extract condensate and non-condensate number from the resulting ensemble, then the condensate wavefunction is and the condensate number is The number of particles not in the condensate is that is, the "vacuum" particles must be subtracted.
It is clear that the "vacuum particles" will contribute to the evolution by means of the nonlinear mixing arising from the particle interactions, and that their effect is quite possibly cutoff dependent. This cutoff dependence must disappear when one includes the full coupling to the noncondensate band by restoring the damping and noise terms in the second and third lines of (43).

Treatment of a non-pure condensate
It is not difficult in principle to include a condensate and its quasiparticles, and this was done in the original treatment [8]. One simply expresses the operators α k in terms of the quasiparticle operators β k for the system in the Bogoliubov theory, and then then one assigns to these the appropriate mean values instead of (90) Of course there will be considerable effort involved in solving the appropriate Bogoliubovde Gennes equations in order to even construct the operators and their corresponding wavefunctions. However, this can be made much easier by using the method of Sinatra et al [38], as explained in their later paper [37].

The choice of E R
It is tempting to choose E R to be so large that the non-condensateband phase space density F (K, x) is negligible, and all of the noise and damping terms in (43) can be neglected. However, this may not be possible, since this would mean that the higher energy parts of the non-condensate band would also need to have negligible occupation, and in this case the neglect of third order derivative terms would certainly not be permissible. Thus it turns out that the principal criterion for the choice of E R should be that there is significant occupation of all quantum states up to this level.

Approximate phenomenological equations
The methods outlined above are quite complicated, and do not entirely answer the need for a simple set of equations which can give a quick estimate of the effects of being studied. Therefore we will consider in here a very simplified equation obtained by making some drastic approximations to the already approximate methodology we have developed.

Mean value equations
Let us therefore neglect the terms involving E ± , for the reasons noted above in Sect.2.3.2 and consider the resulting equations of motion for the mean values: The resulting equation includes a modified Gross-Pitaevskii equation, and a local growth equation. We include the relationship between the backward and forward rates G ± , which arises from the definitions (31,32) Notice that i The terms involving G ± give local growth or decay ofφ andn C , and they wall also cause some dephasing, as was analysed in QKIII and QKIV. When the forward and backward rates balance. ii The spontaneous term-the +1 in (111)-can initiate the condensate growth, while the term proportional ton C gives the difference between stimulated growth and decay, which occur only with nonzeron C . iii The coefficient M [x] has no effect on the densityn C -it is a pure dephasing term of the kind well known in quantum optics and, as can bee seen from (109), it makes the amplitude of the coherent componentφ decay. Since it has no effect on the total number n C , this means that the the effect is to transfer particles from the coherent component into the incoherent or thermalized component of the condensate band. iv) This kind of equation cannot be expected to give a good description of the full condensate growth process, from no atoms in the condensate up to a fully developed condensate in equilibrium with a thermal vapour. As shown in [21,24,25,39] the full growth theory requires a much more detailed description of the kinetics of the thermal cloud, and one cannot simply assume the thermal cloud is always in equilibrium at a definite temperature and chemical potential. However, in [21] we noted that such a simple picture is reasonably valid once the condensate is reasonably large, say 50% of its final value, and in [25] it was noted that the growth process is indeed well described as being at definite temperature, but that the chemical potential one should use is an "effective chemical potential", evaluated from the lower lying energy levels of the thermal vapour. Thus our description is probably a reasonable description of the process of matter exchange between condensate and thermal cloud in a situation not unreasonably far from equilibrium.

The phenomenological growth equation
The hydrodynamic approximation in practice should be defined in terms of the mean wavefunctionφ(x, t), so that, to the extent that we can neglect derivatives of n C (x, t), we can write We further make three further approximations: i) The argument of the exponent in (109) is sufficiently small to use e x ≈ 1 + x. This is in practice almost always true.
ii) We factorize all averages of products of φ operators. iii) We neglect completely the term M [x]φ(x, t) in (110), since it is expected to be small, and its effect is merely to cause a small change in µ.
Doing these, we can replace (110) by where, for consistency with the notation of our previous papers, we have written The function G + is slowly varying in space and time, so that W + can probably be approximated by a constant, since the most essential time-and space-dependence has been included in the factor µφ(x, t) − ih∂φ(x, t)/∂t. The correct equilibrium results from the vanishing of the coefficent of W + , which ensures that the time dependence of the wavefunction is exp(−iµt/h)-the remainng terms in the equation then reduce to the time independent Gross-Pitaevskii equation with chemical potential µ. Thus we get a stationary condensate wavefunction with the same chemical potential as the noncondensate.

Application to vortex array stabilisation
It is numerically known [40,41,42,43] that merely stirring a condensate described by the Gross-Pitaevskii equation produces vortices, but these do not stabilise into an regular array of vortices-as is experimentally observed [44,45,46]-merely as the result of the ongoing progress of the solution of the time-dependent Gross-Pitaevskii equation. To create a vortex array numerically, one solves the time-dependent Gross-Pitaevskii equation with an imaginary time and an appropriately time-dependent added chemical potential, while constantly renormalizing the wavefunction so as to maintain a constant number of particles. This process finds a state with a energy minimum, but of course the method is purely an artefact, and does not represent the true underlying physics. However the form of the phenomenological growth equation (114) already includes what amounts to an imaginary time term, as well as a real time term. We first add an angular momentum term to transform to the rotating frame, and then use a simple form for W + as given in [18] where a is the scattering length for the nonlinear interaction term (so that u = 4πah 2 /m), k is Boltzmann's constant, and T the temperature of the noncondensate). We have included the correction factor g ≈ 3 to give an approximate match with the more detailed treatment, as suggested in [21]. The physics of this situation is the growth of a condensate in a frame rotating with angular velocity Ω about the z-axis from a vapour cloud which is itself stationary in the rotating frame; thus in the laboratory frame, this is vortex nucleation from a rotating vapour cloud, such as has been experimentally implemented by the JILA group [47] with no rotating trap potential. Thus we find that we get Here and using the values a ≈ 10 −8 m and g = 3, we find the factor γ ≈ 0.01. As noted in the previous section, the stationary solution of this eqaution will come from the equality of the coefficients of γ on both sides, ensuring that the time dependence of the wavefunction in the rotating frame is exp(−iµ NC t/h), leaving the wavefunction to satisfy the stationary Gross-Pitaevskii equation (modified by the angular momentum term) with chemical potential µ NC . This stationary solution will be a vortex lattice when it exists.
More generally, we can consider a trap rotating with angular velocity Ω nucleating from a cloud rotating with angular velocity α. In this case we get, in the frame rotating with the trap, The first equation (117) is similar to that presented in the paper of Tsubota et al. [48], but there are significant differences: i) There is no physical justification given in [48], and therefore their results, though qualitatively attractive, cannot be accepted as an explanation of the vortex lattice stabilization process. In contrast, our reasoning shows their choice of γ = 0.03 is of the same order of magnitude as that expected from quantum kinetic theory. ii) They have a real term −µψ on the RHS, whereas this equation has an imaginary term iγµ NC ψ.
iii) Their µ is adjusted with time to preserve the number of atoms in the condensate. Our µ NC is the physical chemical potential of the surrounding vapour, for which we have no particular model at the moment. Even in the experiments [45,46] where an almost pure condensate is stirred by a rotating trap potential, it is known that as the stirring starts, there is considerable heating, and then at the end possibly 20% of the atoms form a thermalized rotating cloud, which appears to rotate at a lower speed than the trap. Thus we can expect that (119) should give a semiquantitative description of the actual process of vortex formation and decay when appropriate (possibly time-dependent) values of µ NC and α are chosen. (However, we note that recent experimental work [49] suggests that the stabilization process is only weakly dependent on temperature, which is surprising, since the dissipation in our equation is represented by γ, which is proportional to temperature. A resolution of this anomaly must await our more detailed calculations which are presently in progress.) The basic mechanism arises from the spatially dependent local chemical potential of a random lattice of vortices. The equilibrium lattice is characterized by a uniform local chemical potential, which must balance that of any surrounding vapour. The phenomenological growth equation includes this irreversible process.
We can expect a model with fixed µ NC to explain the actual observed dissipation of the vibration of vortices which is observed, but the full process of formation and heating could well be a formidable task, involving both kinetic and GP methodologies. Preliminary calculations have shown that qualitatively, the results are very similar to those of Tsubota et al. [48]-that is, the formation and stabilization of a vortex lattice are observed in very much the same way. The numerical techniques necessary for both equations are almost identical 6. Application to hydrogen condensate system

Nonlinear losses
The experiments on hydrogen [50,51,52] give rise to a system with a very large proportion of non-condensed hydrogen, which feeds a relatively small condensate as the condensate atoms themselves leave the condensate because of two-body dipolar relaxation. The data indicate that the condensate appears to have a density profile which matches that of a solution of the Gross-Pitaevskii equation. In contrast, the non-condensed vapour has a profile which does not match that of a zero chemical potential Bose-Einstein distribution, but instead seems to have significantly more population at lower energies than was expected from the zero chemical potential Bose-Einstein distribution. From the point of view of the quantum kinetic theory of condensate growth processes the zero chemical potential model chosen is not very appropriate, since any condensate has a positive chemical potential, and that of the noncondensed atoms must be even higher in order to maintain a net inflow of atoms into the condensate to replace those which are continually leaving by two-body processes. Thus, to the extent that a condensate in this situation can indeed be described by a chemical potential µ C , one expects that the main body of the noncondensate would have a chemical potential µ NC > µ C , and that there would be a transition region for the lower energy particles in the noncondensate. More precisely, what one expects is a distribution function f (E) of the form where the energy-dependent chemical potential µ(E) takes the value µ NC for higher E, but approaches µ C as E → µ C -the situation is illustrated in figure 2. 10 5
Thermodynamic equilibrium, µ = 0 Interpolating Figure 2. Illustration of the interpolation between a phase-space distribution form at high chemical potential µ NC at high energies, and the phase-space distribution at low energies, whose chemical potential µ C is the same as that of the condensate

Phenomenological growth equation including two-body losses
We are most interested in the case of dynamic equilibrium, where the stationary state of the condensate is maintained by the flow of the atoms from the high chemical potential µ NC through the vapour at chemical potential µ C and thence through the process of two body dipolar relaxation to what amounts to chemical potential of −∞, since this process is irreversible. In this case we have to consider three principal processes, namely i) The linear growth process, by which an atom leaves the vapour and enters the condensate.
ii) The linear loss process, by which an atom returns to the vapour from the condensate.
iii) The nonlinear loss process, by which pairs of atoms leave the condensate by dipolar relaxation.
The first two of these terms are already included in the phenomenological growth equation (114), The nonlinear loss is obviously modelled by a term proportional to −|ψ| 2 ψ, leading to a phenomenological growth-loss equation for the hydrogen system: It is a major advantage of the phenomenological approach which we have developed that the nonlinear losses can be accommodated so straightforwardly, with an effectively complex scattering length. The simplicity of the modified Gross-Pitaevskii equation allows quantitative investigation of questions which might well be difficult even to formulate, if one had to proceed directly from first principles without the phenomenological theory as an intermediate step. Yet among such questions are those with potentially dramatic physical consequences, such as the one we now pursue.

Instability of the Thomas-Fermi stationary state under nonlinear loss
Under currently typical experimental conditions, the timescales of condensate growth and decay, including two-body loss, are much longer than those related to the harmonic trapping potential, the mean field energy µ C , or the non-condensate chemical potential µ NC . To take advantage of this small ratio of frequency scales, it is usual to introduce length and time scales related to µ C (the healing length at peak density, and the mean field time scale); but since in the hydrogen context we want to see how the condensate maintains dynamic equilibrium in the presence of a much larger thermal cloud, we will consider µ NC as the fixed parameter, which determines µ C . Since in dynamic equilibrium the two chemical potentials are comparable, this difference is not important. So we will convert to dimensionless variables This also motivates the dimensionless parameters so that we can re-write (121) as Since W + is on the order of the Boltzmann scattering rate σρ tc v T for the thermal cloud, for the hydrogen condensate at temperature approximately 50 µK this gives ε ∼ 10 −6 . And with the published values for the two-body loss rates in this system, the dimensionless parameter Λ is on the order of 10. The analysis below suggests that this is quite large enough to raise serious questions about the stability of the Thomas-Fermi stationary state; and it is possible that the large occupations in the lower quasiparticle levels may mean that Λ ∼ 10 is an underestimate.
In the limit where ε → 0, we obviously obtain the standard Gross-Pitaevskii equation. Since the trapping potential varies very slowly on the healing length scale, in this limit the Thomas-Fermi approximation is excellent. If ε is of order unity or larger, the Thomas-Fermi stationary solution breaks down badly, and it may not even be possible to find a stationary solution to (128). But in the experimental regime of small ε, there are indeed stationary solutions in which the density profile is very close to the standard Thomas-Fermi form. However, unless Λ is small enough, they are dynamically unstable. By itself, the two-body loss term drives the total number of particles towards an equilibrium value. And it even tends to correct against collective excitations of the steady state condensate. In bulk, it is easy to show that this corrective effect strongly stabilizes against perturbations; but as we show below, in a harmonic trap it overcorrects, and the collective excitations grow in amplitude. Numerical integration reveals that these excitations continue in 'boom-and-bust' cycles (where 'bust' means a collapse of central density to a fraction of its maximum value); but in this regime, the assumptions that justify the phenomenological mean field theory may well be breaking down. In an actual experiment, it is not clear whether we should expect coherent but non-stationary states like those assumed in our theory, or some kind of quasi-condensate with degraded phase coherence, or merely accelerated decay of the condensate.

Stationary states
We begin by extending the Thomas-Fermi stationary state to first order in ε. Since it is still true thatṼ is a slowly varying function, we are in the hydrodynamic limit, and we can approximate (128) by the hydrodynamic equations for ψ = √ ρe iθ and v = ∇θ, namelyρ where we drop a term −ερ/(4ρ) in the second equation. (We drop this term because, even when below we allow time-dependent density perturbations around the stationary state, this term will be of order ε 2 . It does become large as we approach the Thomas-Fermi surface, but this only leads to a small correction to the usual boundary layer theory that must be applied to match the hydrodynamic approximation to the near-surface region.) A stationary solution allowsθ = −µ C /µ NC as the only time dependence in ψ, giving for cylindrical co-ordinates and an axisymmetric harmonic trap. So to first order in ε, we have the familiar Thomas-Fermi density profile for an axially symmetric harmonic trap: where we introduce the dimensionless trap frequenciesω j ≡ ω j /µ NC and condensate chemical potentialμ≡ µ/µ NC . But there is now also a velocity field∇θ of order ε.
Since a very small velocity field is difficult to observe, we are mainly interested in the density profile. Actually finding v to first order is therefore only indirectly necessary (it fixes µ C ); but it is not so hard. Assume an ansatz of the form ∂θ ∂r = εr(a r + b rrr for unknown constants a j , b jk . Since∇θ must be irrotational, b rz = b zr . This leads to six equations for only five free parameters in the ansatz, and so forces a specific value on the heretofore undetermined term in ρ, namely the condensate chemical potentialμ (which makes sense, since we do expect the growth and loss terms to determine µ C ). The result for the axisymmetric trap is The uniqueness of this ansatz solution can be proved straightforwardly in the 1D limit, and also for spherical symmetry; we conjecture that it is more generally unique.
The results of (136-138) are somewhat complicated, but they make excellent physical sense. Our expectation that the nonlinear losses will lower µ C below µ NC is borne out. And the velocity field is inward from the edges of the condensate towards the core. Since the loss rate scales faster with the density than the growth rate, the periphery of the condensate receives from the thermal cloud more particles than are needed to maintain the local condensate density, and so it is able to donate a flux to the central regions, where loss rates exceed growth. The velocity field becomes radial and proportional to ρ in the limit of spherical symmetry, but otherwise it has no simple description; see Fig.3 for an illustrative example. In the limit ω z ≪ ω r , we find ∂ z θ becoming independent of r, and the same statement with r and z interchanged is also true, so we do obtain 1D and 2D limits where we expect them.
And it is always true that ∇θ · ∇ρ = 0 at the edge of the Thomas-Fermi cloud, which is required in the hydrodynamic approximation. (Density vanishing outside the TF surface means no flux can cross, since gain and loss terms both vanish at zero density).
So with the addition of the slow flux from periphery to core, the Thomas-Fermi solution is essentially maintained in the presence of slow two-body losses. But is the Thomas-Fermi solution stable?

Collective excitations
To answer this question we need to find the frequencies of linear collective excitations, by linearizing (129,130) about the stationary state (131,132). Using a method which can be dignified as multiple scale analysis, but which is essentially the same as time-independent perturbation theory in ordinary quantum mechanics, we compute the frequencies of the collective modes to first order in ε. The first order corrections are in general imaginary (as is not surprising given that our perturbations are imaginary). What may be surprising is that if Λ exceeds a threshold of order unity, the amplitudes of some collective excitations grow, indicating that the two-body losses cause an over-corrective instability.
The linearized hydrodynamic equations, to first order in ε, are We can combine these into a single second order equation for δθ (again dropping higher order terms in ε): The idea now is to look for δθ = e iΩt δθ(x), and expand Ω = Ω 0 + εΩ 1 + ..., and δθ = δθ 0 + εδθ 1 + .... At zeroth order we have Then using the orthonormality of the eigenfunctions of the zeroth order RHS, we obtain as in the time-independent perturbation theory of quantum mechanics where the integral is over the Thomas-Fermi cloud. (The terms on the RHS with∇θ in them add up to∇ · (δθ 2 0∇ θ), which by Stokes' Theorem is a surface integral that, as we mentioned above, vanishes, since the normal to the surface of the TF cloud is ∝ ∇ρ.) If Λ < 0.7, then our stationary solution is stable because the integrand in (142) is positive throughout the volume of integration; but if Λ > 0.7, then the integrand is negative in a volume around origin, and so the integral might perhaps be negative at least for some modes.
To investigate further we must learn the δθ 0 for the various collective modes.
Solving (141) for the hydrodynamic modes of a harmonic trap can be reduced to diagonalizing finite matrices using a Frobenius series approach. In an axially symmetric trap it is not hard to find modes up to third order polynomials analytically and exactly, and with spherical symmetry we can find all the modes. In an extremely prolate (cigar-shaped) trap, we can also find all the modes to zeroth order in (ω z /ω r ) 2 . We will proceed here with the spherical calculation. The extremely prolate case is analysed in the Appendix, and here we will only quote its results.

Spherical trap
In the spherical case it is convenient to solve (141) in spherical polar co-ordinates, with the angular dependence of δθ 0 obviously being a spherical harmonic Y lm .
With ω z = ω r = ω, define the rescaled spherical radius and write δθ 0 = Y lm (ϑ, φ)f nl (r) so that (141) becomes (The reason for the additional subscript n on f nl will appear shortly.) Taking a Frobenius ansatz f nl = a knl r k+2 for the radial dependence then yields the recursion relation Regularity at the origin then forces s = l, and further requires that a knl vanish for all k with parity opposite to l. Convergence as we approach the Thomas-Fermi surface r = 1 forces the sequence of non-zero a knl to terminate at some k = n − l ≥ 0. This introduces the quantum number n, which must have the same parity as l, and be greater than or equal to it. We recover the familiar result the point of repeating the calculation to this point has been to recall the recursion relation (145), from which we can evaluate the RHS of (142). In our rescaled spherical co-ordinates the equation for Ω 1 becomes (147) We can evaluate the ratio of moments in this expression by using the orthogonality under the weight r 2 of the eigenfunctions of (144) of different n, together with a fact we can extract from the recursion relation (145). Since the f nl are polynomials of rank n and definite parity, it follows that f n+2,l = n−l k=0 a k,n+2,l r k+l (148) = a n+2−l,n+2,l a n−l,nl r 2 f nl (149) + a n+2−l,n+2,l a n−l,nl a n−l,n+2,l a n+2−l,n+2,l − a n−2−l,nl a n−l,nl f nl (150) for some constants C j . This then means that r 2 f nl = a n−2−l,nl a n−l,nl − a n−l,n+2,l a n+2−l,n+2,l f nl + j =n so that by orthogonality we have 1 0 dr r 4 f 2 nl 1 0 dr r 2 f 2 nl = a n−2−l,nl a n−l,nl − a n−l,n+2,l a n+2−l,n+2,l (153) where the evaluation follows straightforwardly from (145). This immediately yields our desired result .
¿From this we can see that the so-called 'surface' modes, with l = n, are always stable. (The case n = l = 0 is not an excitation, but simply the stationary state, so this case does not count.) Since if Λ < 1/2 then ImΩ sphere 1 will be positive because 1 − 6Λ/7 > 1 − 2Λ, we need only look for instabilities in cases where 2Λ − 1 > 0. In these cases it is clear that instability can only occur for Λ > 7/6, in which case it occurs first for smaller l and larger n. Equation (155) indicates that for Λ > 7/6 there will always be instabilities at sufficiently high n, but the hydrodynamic approximation on which the equation is based will break down for n larger than some limit that depends on µ C and ω, and so in general the actual instability threshold will be somewhat higher. But for l = 0, all modes with n ≥ 2 will be unstable if Λ exceeds only 77/64, and so the instability threshold is not actually very sensitive to post-hydrodynamic corrections.
where |m|, n, p are whole numbers. The azimuthal quantum number m must have the same parity as n, and |m| ≤ n. In this case, all of the axial modes (n = m = 0) are stable even if Λ → ∞. And the dipole modes (n = m = 0, p = 1 and |m| = n = 1, p = 0) are not only stable, but have Ω 1 independent of Λ. (This is also true in the spherical case, and indeed for all harmonic traps, because the dipole modes merely translate the entire condensate, and so do not disturb the local balance between growth, loss, and flux.) Raising |m| and p tends to stabilize; the most unstable modes will have p = m = 0. For these modes we can recognize that a dynamical instability occurs for (6n + 8)Λ > (7n + 14), which will occur at high n for Λ > 7/6, just as in the spherical case. The same caution applies, that for n too large the hydrodynamic approximation breaks down. But for Λ > 7/5, all modes with p = m = 0, n ≥ 2 will be unstable, and so again post-hydrodynamic effects cannot shift the instability threshold very much. In the hydrogen condensate experiment at MIT, where the trap has an aspect ratio of 400, the extremely prolate limit definitely obtains, and with Λ ∼ 10 there are clearly a lot of unstable modes. We conclude that instability of the Thomas-Fermi stationary state is a general phenomenon if there are sufficiently strong two-body losses. The important point here is that 'sufficiently strong' means only that the two-body rate constant be slightly greater than the growth rate constant.

Implications
Since the MIT hydrogen condensate experiments are well above the instability threshold, our phenomenological growth equation seems to indicate that the quasi-steady state of these condensates cannot be the Thomas-Fermi density profile with perfect phase coherence, which all other condensates exhibit well. The problem is not so much that the Thomas-Fermi state should be oscillating unstably, but that the system cannot be expected to settle down to an unstable state in the first place. How seriously should we take this conclusion, which is after all derived from approximate solutions to a phenomenological theory? We would argue that, even if the theory might not be quantitatively precise, the overcorrective instability which it predicts under twobody losses is a simple physical phenomenon. When loss and gain processes scale differently with density, and density is inhomogeneous, particle flow is needed to maintain a steady state. In effect the condensate velocity field is a control system that tries to maintain a constant density profile. But moving atoms have momentum, and so the problem of overcorrection can obviously arise in this control system; whether this leads to instability is a detailed matter of length and time scales. It is possible that our approximate theory exaggerates the overcorrection problem, or that it may be counteracted by factors not included in the theory; but it is a genuinely physical possibility, and not a mere artifact of our phenomenological approach.
Among factors neglected in this Section, which might tend to suppress instabilities in real systems, is the thermal component of the velocity field. In the hydrogen experiments the thermal fluctuations in the condensate velocity field can be estimated to be much larger than the systematic flow which compensates for the two-body losses, and so these should really be taken into account, by solving the full, stochastic version of our modified Gross-Pitaevskii equation.
It is a great merit of the phenomenological growth equation that it allows explicit calculation of the instability question, which would be difficult even to formulate directly from first principles. But given the experimental data showing a high-density component of ultracold hydrogen, whose density profile seems consistent with those of condensates in the alkali vapours, the phenomenological growth equation cannot be correct when it concludes that the Thomas-Fermi ground state is impossible in hydrogen. The results of this Section do indicate, however, that the balance between thermal fluctuations and dissipation (loss) must be much less trivial in hydrogen than in the other ultracold Bose gases. Thus in this Section we have provided evidence that further development of non-equilibrium theory for cold, dilute, trapped bosons will provide insights into qualitatively new regimes, and not just small corrections to familiar results. At the same time the stochastic Gross-Pitaevskii approach has proven itself as a workable tool which can, without requiring unreasonable effort, yield definite answers to nontrivial physical questions.

Conclusion
In this paper we have given indications of how a stochastic Gross-Pitaevskii equation should be defined and used in practice. There are two major issues to be resolved in this formulation a) What is the appropriate definition of the "condensate wavefunction" ψ(x, t) for which the stochastic Gross-Pitaevskii equation should provide an equation of motion? Our conclusion is that a multimode Wigner function definition is the most appropriate, and although this cannot be done exactly, or even for all possible quantum states, the approximations and restrictions which are thus implicit in this definition appear to us to be unlikely to cause in serious problems when applied to practical situations. The most irksome complication is the necessity to include "vacuum noise" in the initial conditions, basically in order to make sure the Heisenberg uncertainty principle is not violated, as detailed in Sect.3.3. However, at all but the very lowest temperatures this is unlikely to be an issue-the influence of invalid initial conditions will very soon be eliminated by the noise induced by the thermalized atoms. b) The formulation a quantum kinetic master equation using the local formulation of energy and momentum conservation introduced by Zaremba et al [11] creates a very much simpler master equation than that of our earlier formulations, and its transformation into the stochastic Gross-Pitaevskii equation is then relatively painless. There are, however, many technical details involving the best way of quantitatively specifying the noise coefficients which are still in need of refinement. We do not see these as being too important quantitatively in practical situations, but they will eventually need to be attended to more precisely than we have done here.
The major result is that, subject to a number of caveats, including the correct choice of initial conditions, we can modify the Gross-Pitaevskii equation to include damping as a result of exchange with non-condensed atoms, in order to obtain a semiphenomenological description of the interacting systems which includes all of the quantum effects in at least an approximate form.
The development of the phenomenological growth equation, is similar to one presented by Williams and Griffin [26], is the most useful immediate application of this work. This equation can be seen as a very simplified description of the interaction of a condensate and a thermal cloud, which includes the major processes in a consistent but only semiquantitative way. The formulae for vortex lattice nucleation and stabilization given in Sect.5 are new, and in work which we shall publish elsewhere, we will show how a rich variety of behaviours arise out of this simple and elegant formulation.
We have also exploited this equation investigate the description of growth and loss in a hydrogen condensate, and have shown that there are hitherto unremarked hydrodynamic instabilities in the Thomas-Fermi stationary state of a cold condensate suffering two-body losses. We shall publish more detailed numerical studies of this system elsewhere.
We can use (162) to show that solutions F n,m with different n are orthogonal in the interval X ∈ [0, 1] with the weight X. Consequently, we can evaluate moments such as dX X 3 F 2 nm / dX XF 2 nm using the recursion relation (163), by observing that X 2 F n,m = k a n;k−2 X k = a n,n a n+2,n+2 F n+2,m + a n,n−2 a n,n − a n+2,n a n+2,n+2 F n,m + n−2 l c l F l,m for some c l . This implies that . (165) We can now go back to (160) to first order in η, and by integrating over X with XF nm , project out the equation for G(Z). (We have to integrate by parts several times, and apply (165) as well.) In fact we find a simpler equation for g(Z) ≡ G(Z)(1 − Z 2 ) −n/2 (which is reassuring since after all X n = R n (1 − Z 2 ) −n/2 , and so F nm (X)(1 − Z 2 ) n/2 g(Z) will be regular at R = 1, Z → ±1 if g is regular at Z = ±1). The result is where a fairly complicated bunch of terms involving m, n and ν 1 has been absorbed into the eigenvalue ξ. Similar analysis to that above shows ξ = p(p + 2n + 3) for non-negative integer p, determining ν 1 and hence Ω 0 as a function of the mode indices m, n, p. We have thus found the complete hydrodynamic spectrum of an extremely elongated harmonic trap. By a similar analysis as already performed for the F nm above, we can also compute the second moment of g(Z) (which are orthogonal over the interval Z ∈ [−1, 1] with the weight (1 − Z 2 ) n+1 ). So we also finally have the correction to the frequency to first order in ε: p, n|(1 − Z 2 )|p, n Z p, n|p, n Z n, m|(1 − X 2 )|n, m X n, m|n, m X (167) = iμ 2 4 7 Λ + (1 − 2Λ) p 2 + p(2n + 3) + (n + 2)(2n + 1) (2p + 2n + 1)(2p + 2n + 5) 1 − m 2 n(n + 2) . (168)