text:
# Chaos and the continuum limit in the gravitational 𝑁-body problem. I. Integrable potentials

## I INTRODUCTION AND MOTIVATION

Many astronomical objects, including, e.g., globular clusters, are typically modeled by bulk gravitational potentials which manifest a high degree of symmetry and which, being integrable, lead to completely regular characteristics with no possibility of chaotic behaviour. One knows, however, that such bulk potentials constitute idealisations, the true system corresponding (at least approximately) to a realisation of the gravitational $`N`$-body problem. The important point, then, is that motion in the $`N`$-body problem, even for an $`N`$-body system which samples a smooth, time-independent phase space distribution corresponding to an integrable potential, is typically chaotic in the sense that orbits exhibit exponential sensitivity towards small changes in initial conditions. This perhaps is not surprising. The true potential associated with a collection of point masses no longer possesses the symmetries of the original integrable potential, so that there is no reason why the orbits should not be chaotic.

However, what is, perhaps, surprising is the expectation, derived both from theoretical arguments  and from numerical simulations, that the $`N`$-body problem remains chaotic even for very large $`N`$. Suppose, e.g., that a system of total mass $`M=1`$ is represented by a collection of $`N`$ objects of mass $`m=1/N`$, so distributed as to sample the density distribution corresponding to an integrable potential. The claim then is that, when expressed in units of a natural dynamical, or crossing, time $`t_D\mathrm{\hspace{0.33em}1}/\sqrt{G\rho }`$, with $`\rho `$ a typical density, the characteristic time scale $`\tau `$ on which an initial perturbation in any given orbit tends to grow will not diverge for $`N\mathrm{}`$. In this sense, the degree of chaos manifested by individual orbits is not expected to ‘turn off’ for very large $`N`$. There is an apparent consensus, motivated both from theory and numerical experiments, that $`\tau `$ should not increase without bound for $`N\mathrm{}`$, although there is some disagreement in the literature as to whether $`\tau (N)`$ should converge towards an $`N`$-independent value  or whether $`\tau `$ should instead slowly decrease with increasing $`N`$.

If, however, this be true, one is confronted with subtle questions of principle regarding the nature of the continuum limit. It is generally assumed  that, for sufficiently large $`N`$, a self-gravitating system of discrete point masses can be characterised adequately by a smooth phase space density that solves the collisionless Boltzmann equation (CBE), i.e., the gravitational analogue of the Vlasov equation from plasma physics. The obvious point, then, is that time-independent solutions to this equation which manifest a high degree of symmetry correspond typically to bulk potentials which are integrable or, even if they be nonintegrable, admit large measures of regular orbits. But how is one to reconcile integrable or near-integrable behaviour in such bulk potentials with the presumed fact that, even for very large $`N`$, individual orbits in the true $`N`$-body problem typically manifest chaotic behaviour on a time scale $`t_D`$?

Strictly speaking, there is no logical contradiction: It is completely possible for collective properties of an $`N`$-body system to be described correctly by the CBE, even if the characteristics associated with the self-consistent potential do not coincide, even approximately, with real $`N`$-body trajectories. However, it would seem important to pin down carefully what is actually going on:

$``$ Is it really true that individual trajectories in the $`N`$-body problem are chaotic for very large $`N`$, even if the bulk potential associated with the system is integrable? The indications are that the answer to this is: yes. However, most of the work done to date on chaos in the $`N`$-body problem has focused on systems with comparatively small $`N`$ and/or a hierarchy of masses, or, for larger systems, on comparatively short time behaviour. Little if any work has been done to provide estimates of honest Lyapunov exponents over intervals $`t_D`$ for large $`N`$ systems comprised of bodies of comparable mass.

$``$ Does this chaos reflect simply the fact that the system is grainy, or does it reflect the details of the full $`N`$-body dynamics? If, e.g., one were to replace a smooth time-independent, integrable potential by the irregular potential associated with an $`N`$-particle sampling of the smooth mass density that is frozen in time, to what extent will motion in that grainy time-independent potential be chaotic?

$``$ Even presuming that the $`N`$-body problem is chaotic on a time scale $`t_D`$, is there some well defined sense in which the ‘average’ properties of individual $`N`$-body trajectories track characteristics given by the CBE?

These conceptual issues are also related directly to the problem of ‘softening.’ It is generally recognised that, for small $`N`$, close encounters between individual masses are more important dynamically than for larger $`N`$. For this reason, $`N`$-body simulators interested in exploring the physics of the $`N`$-body problem for larger $`N`$ often suppress the effects of close encounters artificially by replacing the true $`1/r`$ potential by a softened potential $`V(r)(r^2+ϵ^2)^{1/2}`$ for some ‘softening parameter’ $`ϵ`$. This certainly suppresses encounters with impact parameters $`<ϵ`$ which, presumably, is a good thing. However, there are strong indications that orbits integrated with such a softened potential tend to be ‘less chaotic’ in their behaviour, so that the introduction of softening also has the potentially undesirable effect of removing $`N`$-body chaos which really ought to be present, even for very large $`N`$. In any event, earlier investigations of chaos in the $`N`$-body problem based on simulations which incorporate a large amount of softening must be viewed with suspicion, since such simulations could suppress precisely the effects which one might wish to explore!

This paper summarises a detailed exploration of chaos in time-independent potentials generated by sampling the smooth density $`\rho (𝐫)`$ associated with a time-independent solution to the CBE to create a frozen $`N`$-body realisation of that equilibrium. Most of the work focuses on the particularly simple case of an integrable Plummer potential, which derives from a spherically symmetric mass distribution. However, it was also confirmed that, modulo one point discussed in the concluding section, the same qualitative results obtain for the potential associated with a constant density spherical configuration.

Section II begins by describing the numerical experiments that were performed. Section III summarises a computation of honest Lyapunov exponents in frozen $`N`$-body realisations of the Plummer potential, exploring how the largest exponent $`\chi `$ associated with representative initial conditions varies as a function of $`ϵ`$ and $`N`$. The principal conclusion here is that, at least for small values of $`ϵ`$, orbits in such potentials are invariably chaotic; and that, even for particle number as large as $`N=10^{5.5}`$, there is no sense in which increasing $`N`$ ‘turns the chaos off.’ Section IV demonstrates that, even though the Lyapunov exponents do not decrease with increasing $`N`$, there is a well defined sense in which, as $`N`$ increases, orbits in frozen-$`N`$ potentials remain ‘close’ to smooth potential characteristics with the same initial condition for progressively longer times. Section V concludes by summarising the principal conclusions, providing a simple physical interpretation, and then commenting on potential implications.

The principal conclusion of this paper is that, for integrable smooth potentials which admit no chaos, the continuum limit makes sense even at the level of pointwise properties of individual trajectories. The possibility of chaotic characteristics leads necessarily to very different behaviour and, for this reason, the case of nonintegrable potentials that admit both regular and chaotic characteristics will be considered in a separate paper.

## II DESCRIPTION OF THE NUMERICAL EXPERIMENTS

The numerical computations reported here were performed for a so-called Plummer potential,

$$\mathrm{\Phi }(r)=\frac{GM}{\sqrt{r^2+b^2}}.$$
(1)
This potential is generated via Poisson’s equation from a density profile

$$\rho (r)=\left(\frac{3M}{4\pi b^3}\right)\left(1+\frac{r^2}{b^2}\right)^{5/2},$$
(2)
and corresponds to an equilibrium solution to the CBE satisfying

$$f(E)=\{\begin{array}{cc}A(E)^{7/2}\hfill & \text{if }\mathrm{\Phi }(r=0)<E=\frac{1}{2}v^2+\mathrm{\Phi }<0\text{;}\hfill \\ \mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}\hfill & \text{if }E=\frac{1}{2}v^2+\mathrm{\Phi }>0\text{.}\hfill \end{array}$$
(3)
Units were so chosen that $`G=M=b=1`$.

The principal aim was to compare orbits generated in the smooth potential with orbits evolved in time-independent $`N`$-body realisations of the potential. For a variety of fixed values of $`N`$ and $`ϵ`$, $`20`$ different time-independent $`N`$-body potentials were constructed. Each of these was associated with a random sampling of the smooth density distribution generated using a von Neumann rejection algorithm (cf. ). This entailed constructing singular density distributions

$$\rho _N(𝐫)=\frac{M}{N}\underset{i=1}{\overset{N}{}}\delta _D(𝐫𝐫_i),$$
(4)
which, allowing for a softening parameter $`ϵ`$, yielded potentials of the form

$$V_N(𝐫)=\frac{GM}{N}\underset{i=1}{\overset{N}{}}\frac{1}{\sqrt{(𝐫𝐫_i)^2+ϵ^2}}.$$
(5)
The objective then was to select individual initial conditions $`(𝐫_0,𝐯_0)`$ and to evolve these initial conditions in both the smooth potential and the $`20`$ ‘frozen’ $`N`$-body potentials, while simultaneously tracking the evolution of a small initial perturbation, periodically renormalised at fixed intervals $`\delta t`$, so as to extract an estimate of the largest (short time) Lyapunov exponent .

The integrations were performed for a time corresponding physically to $`\mathrm{\hspace{0.33em}100}t_D`$ using a Runge-Kutta integrator that typically conserved energy to at least one part in $`10^4`$. The value $`100t_D`$ was selected (i) because it corresponded to an interval sufficiently long that one began to see convergence towards a well-defined Lyapunov exponent $`\chi `$ and, perhaps more importantly, (ii) because, for physical systems like real galaxies, $`100t_D`$ corresponds to an interval comparable to the age of the Universe.

The total particle number $`N`$ in the ‘frozen’ $`N`$-body potentials was allowed to vary between $`N=10^{2.5}`$ and $`N=10^{5.5}`$. Physical interest focuses primarily on the limit $`ϵ0`$, this corresponding to an ‘honest’ $`N`$-body calculation. However, the effects of a nonzero $`ϵ`$ were also considered in some detail, with the aim of ascertaining possible undesirable consequences for conventional $`N`$-body simulations, which typically involve a substantial softening. The experiments with variable $`ϵ`$ indicated that, for $`ϵ<10^4`$ or so, the precise value of $`ϵ`$ was largely immaterial, at least statistically.

## III SHORT TIME LYAPUNOV EXPONENTS

The principal diagnostic here was the mean (short time) Lyapunov exponent $`\chi `$, generated, for a given choice of initial condition and for specified values of $`ϵ`$ and $`N`$, as the average value of $`\chi `$ at $`t=100t_D`$ for $`20`$ different frozen-$`N`$ potentials. The fundamental question was how, for fixed initial condition, this $`\chi `$ depends on $`ϵ`$ and $`N`$. FIGURE 1 exhibits $`\chi `$ as a function of $`\mathrm{log}_{10}ϵ`$ for multiple integrations of one representative initial condition, which corresponded initially to a roughly isotropic distribution of velocities, allowing for several different values of $`N`$. FIGURE 2 gives $`\chi `$ as a function of $`\mathrm{log}_{10}N`$ for the same initial condition, now allowing for several different values of $`ϵ`$.

It is evident from FIGURE 1 that, at least for comparatively large values of softening parameter, decreasing $`ϵ`$ tends to make the orbit more chaotic. This is hardly surprising: Since the smooth potential is integrable, one anticipates that the chaos is associated completely with close encounters between the test mass and individual frozen masses. The introduction of a nonzero smoothing corresponds de facto to the introduction of a minimum impact parameter (since the potential is bounded in magnitude by $`V_{max}=GM/Nϵ`$) but the existence of such a minimum impact parameter limits the maximum effect that can arise from a close encounter.

However, for sufficiently small values of $`ϵ`$, the precise value of $`ϵ`$ appears to be largely immaterial. This again is hardly surprising: As long as $`ϵ`$ is small compared with the value of the closest separation between the test particle and any of the frozen particles during the course of the integration, the test particle feels an essentially unsoftened potential and should behave (at least statistically) as if $`ϵ=0`$. The point then is that, for $`N\mathrm{\hspace{0.33em}10}^6`$ and an integration time as short as $`100t_D`$, the minimum separation associated with the closest encounter between the test mass and any of the frozen masses should be greater than or comparable to $`ϵ\mathrm{\hspace{0.33em}10}^4`$. Indeed, a simple geometric argument indicates that the time scale $`t_ϵ`$ for a close encounter with minimum separation as small as $`ϵ`$ scales as

$$\frac{t_ϵ}{t_D}\frac{R_{sys}^2}{Nϵ^2},$$
(6)
where $`R_{sys}`$ is the size of the system in question.

One obvious implication of these results is that the introduction of a large amount of softening into a numerical simulation can have the unnatural result of significantly decreasing the amount of chaos manifested by individual orbits in a real astronomical system.

For comparatively large values of $`ϵ`$, $`\chi `$ decreases rapidly with increasing $`N`$ but, for sufficiently small values of $`N`$, it appears that $`\chi `$ is nearly independent of $`N`$ (although there are hints that $`\chi `$ may continue to increase very slowly). The fact that, for large $`ϵ`$, $`\chi `$ should decrease with increasing $`N`$ can again be explained by comparing the magnitude of $`ϵ`$ with the typical distance between masses in the system, which is of order $`n^{1/3}R_{sys}/N^{1/3}`$, with $`n`$ a characteristic number density. If $`ϵ`$ is larger than, or comparable to, $`n^{1/3}`$, even weak close encounters are essentially ‘turned off,’ so that the source of chaos has been largely reduced, if not completely removed. The fact that $`\chi `$ should be essentially independent of $`N`$ in the limit $`ϵ0`$ has been argued by various authors in a number of different ways. A simple physical explanation is provided in the concluding section.

If a single orbit be integrated for progressively longer times, how quickly will the short time Lyapunov exponent $`\chi (t)`$ converge towards the true time-independent $`\chi `$? Studies of orbits in smooth nonintegrable potentials reveal that, when the phase space is highly complex and, because of the Arnold web, orbits can be ‘stuck’ temporarily in regions where the short time Lyapunov exponents are especially small or especially large, the time required for a reasonable level of convergence can be extremely long, $`\mathrm{\hspace{0.33em}10}^5t_D10^6t_D`$ or even more. If, however, the phase space is simpler in the sense that the Arnold web forms less of an impediment and such trapping is comparatively rare, the time required is typically much shorter. One way in which to quantify the overall rate of convergence is by performing a simple time series analysis: An orbit segment of length $`T`$ can of course be divided into $`k`$ segments of length $`\mathrm{\Delta }t=T/k`$ and a short time Lyapunov exponent $`\chi (\mathrm{\Delta }t)`$ computed for each segment. The dispersion $`\sigma _\chi (\mathrm{\Delta }t)`$ then provides a useful probe of the degree to which, on time scales $`\mathrm{\Delta }t`$, the degree of chaos exhibited by different orbit segments is more or less the same. Determining $`\sigma _\chi `$ as a function of $`\mathrm{\Delta }t`$ provides a quantitative characterisation of the rate of convergence towards a unique $`\chi _{\mathrm{}}`$. A simple argument based on the Central Limits Theorem suggests that, if the accessible phase space regions are simple and trapping is rare, so that the amounts of chaos manifested at times $`t`$ and $`t+\mathrm{\Delta }t`$ are essentially uncorrelated,

$$\sigma _\chi (\mathrm{\Delta }t)^p,$$
(7)
with $`p=1/2`$. If, alternatively, the phase space is complex and trapping is important, one would expect that $`\sigma _\chi `$ decreases much more slowly with increasing $`\mathrm{\Delta }t`$.

Such a time series analysis was performed for the data sets used to generate the mean exponents $`\chi `$. For each set of $`20`$ integrations, each orbit segment of length $`T=100t_D`$ was separated into $`k`$ segments of length $`\mathrm{\Delta }t=T/k`$. A short time Lyapunov exponent $`\chi (\mathrm{\Delta }t)`$ was then computed for each of the resulting $`20k`$ segments, and these were used to compute the dispersion $`\sigma _\chi (\mathrm{\Delta }t)`$. Allowing for $`k=2^q`$, for $`q=0,1,2,3,4,5,`$ and $`6`$ was equivalent to varying $`\mathrm{\Delta }t`$ between $`\mathrm{\Delta }t=(100/64)t_D`$ and $`\mathrm{\Delta }t=100t_D`$. This time series analysis led to the conclusion that the dispersion $`\sigma _\chi `$ is typically well fit by a power law dependence of the form given by eq. (3.2), although the exponent $`p`$ tends to be somewhat smaller than $`p=1/2`$, the best fit value typically satisfying $`p\mathrm{\hspace{0.33em}0.4}`$. Several examples are exhibited in FIGURE 3. The fact that $`p`$ is comparatively close to $`1/2`$, rather than the much smaller values that are often observed in very ‘sticky’ nonintegrable potentials, corroborates the intuition that, because the chaos in this problem is associated exclusively with close encounters, trapping is rare and the degree of chaos exhibited at different times tends to be statistically uncorrelated.

## IV COMPARISON OF SMOOTH AND N-BODY ORBITS

It is clear that, for sufficiently short times, a frozen-$`N`$ orbit will coincide almost exactly with the smooth potential characteristic associated with the same initial condition. And similarly, it is clear that, at sufficiently late times, the irregularities in the frozen-$`N`$ potential will cause the frozen-$`N`$ orbit to deviate significantly from the smooth characteristic. Probing the validity of the continuum limit at the level of individual orbits thus devolves into determining the rate at which the frozen-$`N`$ orbits and smooth characteristics diverge. In this connection, two obvious questions arise. Do frozen-$`N`$ orbits diverge from the smooth characteristics exponentially or as a power law in time? And how does the overall rate of divergence depend on $`N`$?

Such probes of the validity of the continuum limit differ from the ordinary point of view, where convergence is typically defined in terms of quantities like bulk moments of the system, ignoring completely the behaviour of individual trajectories. A possible intermediate characterisation is to focus not on the pointwise behaviour of the chaotic orbits but, instead, on quantities which might be less sensitive to the $`N`$-body chaos. In particular, one can also ask: How do frozen-$`N`$ orbits deviate from smooth characteristics in terms of quantities which, in the smooth potential, correspond to time-independent constants of the motion, like angular momentum in a spherically symmetric system?

These questions were addressed here both visually and quantitatively through a computation of the statistical properties of frozen-$`N`$ orbits. Given $`n=20`$ different trajectories $`\{(𝐫_i(t),𝐯_i(t))\}`$, $`i=1,\mathrm{},n`$, and the smooth characteristic $`(𝐫_s(t),𝐯_s(t))`$ associated with the same initial condition, there are two types of moments which one might choose to consider. Quantities like

$$𝐫=\frac{1}{n}\underset{i=1}{\overset{m}{}}𝐫_i$$
(8)
and

$$Dr^2=|𝐫_i𝐫|^2\frac{1}{n}\underset{i}{}|𝐫_i𝐫|^2$$
(9)
and the corresponding quantities generated from $`𝐯`$ and $`𝐉=𝐫\times 𝐯`$ focus on the frozen-$`N`$ orbits in and of themselves. Alternatively, such moments as

$$\mathrm{\Delta }r^2|𝐫𝐫_s|^2=\frac{1}{n}\underset{i}{}|𝐫_i𝐫_s|^2$$
(10)
and

$$\delta r^2=|𝐫𝐫_s|^2$$
(11)
compare the frozen-$`N`$ orbits with the smooth potential characteristic and, as such, their behaviour as a function of $`N`$ is particularly relevant in understanding the continuum limit. Overall, the quantities $`Dr`$, $`\delta r`$, and $`\mathrm{\Delta }r`$ were found to exhibit comparatively similar evolutions, so that attention below focuses on the moments $`𝐫`$ and $`\mathrm{\Delta }r`$, which seem especially natural physically.

The most striking conclusion is that individual frozen-$`N`$ orbits typically diverge from the smooth characteristic as a power law in time, not exponentially. This is true both for comparatively large values of $`ϵ`$, where the frozen-$`N`$ orbits are nearly regular, and for smaller values of $`ϵ`$, where the orbits are much more chaotic. This result is perhaps surprising. One might naively have supposed that, since the frozen-$`N`$ orbits are strongly chaotic, at least for small $`ϵ`$, they would tend to diverge exponentially from the smooth characteristics on a time scale $`\tau \chi ^1`$. However, such an exponential divergence is most definitely not observed.

For how long does this power law divergence persist? Does it cease when the distance between the frozen-$`N`$ orbit and the smooth characteristic is still small, or does the divergence continue until the frozen-$`N`$ orbit and the smooth characteristic tend to be widely separated in configuration space? If, e.g., this divergence terminated at comparatively small separations, much smaller than the size of the system, one could argue that, even though the frozen-$`N`$ orbits are chaotic, they still remain ‘close’ to the smooth characteristics. The answer here is that this divergence continues until the typical separation has become comparable to the size of the configuration space region to which the orbits are confined and the frozen-$`N`$ orbit has become completely ‘decorrelated’ in appearance from the smooth potential characteristic.

The same conclusion also obtains if one focuses on an ensemble of $`20`$ frozen-$`N`$ orbits and probes their statistical properties. The six panels of FIGURE 4, generated for the initial condition used in FIGURES 1 - 3, compare $`x`$ for frozen-$`N`$ ensembles with the smooth $`x_s`$ for orbits evolved with $`ϵ=10^5`$, allowing for six values of $`N`$ extending from $`N=316`$ to $`N=100000`$. In each case, one finds that, for sufficiently large $`t`$, $`x0`$, as would be expected if the frozen-$`N`$ orbits have become completely different from one another and move through configuration space with random orientations. FIGURE 5 compares the radial coordinates $`r`$ and $`r_s`$ for the extreme case of an initial condition corresponding in the smooth potential to a purely radial orbit.

The time scale $`t_G`$ on which the frozen-$`N`$ orbits diverge from the smooth characteristic, and hence the time scale on which $`\mathrm{\Delta }r`$ grows, increases with increasing $`N`$. Even though the frozen-$`N`$ orbits remain ‘equally chaotic’ in the sense that their Lyapunov exponents $`\chi `$ remain nearly constant, they remain close to the smooth characteristic for progressively longer times.

The left hand panels of FIGURE 6 exhibit $`\mathrm{\Delta }r/2^{1/2}R_s`$, with $`R_s^2`$ the mean value of $`r^2`$ associated with the smooth characteristic, as computed for the same initial condition evolved with $`ϵ=10^4`$ for $`N=1000`$ and $`N=100000`$. The right hand side exhibits the same data, recorded at intervals of $`0.025t_D`$, once they have been subjected to a boxcar averaging over an interval $`\delta t=1t_D`$. The large envelops associated with the curves in the left hand panels reflect, e.g., the fact that, at late times, individual orbits in the $`n=20`$ orbit ensembles are oscillating about a value of unity.

That $`\mathrm{\Delta }r/2^{1/2}R_s`$ converges towards unity is a reflection of the fact that the orbits have indeed become completely different one from another: Given that the frozen-$`N`$ orbits conserve energy, one might expect that their average distance from the origin should, on the average, be the same as for the smooth characteristic, so that

$$r^2(t)R_s^2\mathrm{for}t\mathrm{}.$$
(12)
Assuming, however, that this be true and that

$$𝐫(t)𝐫_s(t)0\mathrm{for}t\mathrm{},$$
(13)
one infers that $`\mathrm{\Delta }𝐫^22R_s^2.`$ Analogous behaviour is observed for the quantity $`\mathrm{\Delta }v/2^{1/2}V_s`$, with $`V_s^2`$ defined correspondingly for the smooth characteristic.

As is manifested by FIGURE 6, the growth of $`\mathrm{\Delta }r`$ and $`\mathrm{\Delta }v`$ is roughly linear in time. Indeed, when comparing an ensemble of frozen-$`N`$ orbits with a smooth orbit characteristic generated from the same initial condition, one finds that, for small $`ϵ`$, $`\mathrm{\Delta }r`$ and $`\mathrm{\Delta }v`$ are both reasonably well fit by a linear growth law of the form

$$\frac{\mathrm{\Delta }r}{r_s}=\frac{\mathrm{\Delta }v}{v_s}=\frac{t}{t_G}.$$
(14)
The growth time $`t_G`$, which is the same for both $`\mathrm{\Delta }r`$ and $`\mathrm{\Delta }v`$, satisfies

$$t_GA_GN^pt_D,$$
(15)
with $`A_G`$ of order unity and $`p\mathrm{\hspace{0.33em}1}/2`$. FIGURE 7 exhibits $`\mathrm{log}_{10}(t_G/t_D)`$ as a function of $`\mathrm{log}_{10}N`$ for two different initial conditions evolved with $`ϵ=10^5`$.

The fact that $`t_G`$ scales as $`N^{1/2}`$ would suggest that the divergence of the frozen-$`N`$ orbits from smooth characteristics reflects a diffusion process, associated with a collection of random close encounters. However, this might in turn suggest that $`\mathrm{\Delta }r`$ and $`\mathrm{\Delta }v`$ should grow as $`t^{1/2}`$, rather than the approximately linear growth that was observed in the numerical simulations. Interesting, though, such a $`t^{1/2}`$ behaviour does obtain for quantities like angular momentum, which are conserved absolutely in the smooth potential. Indeed, one finds that, for small $`ϵ`$, $`\mathrm{\Delta }J`$ satisfies

$$\frac{\mathrm{\Delta }J^2}{J_s^2}=\frac{t}{t_J},$$
(16)
where $`J_s`$ is the typical magnitude of the angular momentum associated with a characteristic with the specified energy. The growth time $`t_J`$ again scales as $`N^{1/2}`$, but tends to be somewhat larger than $`t_G`$, so that

$$t_JA_JN^pt_D,$$
(17)
with $`A_J\mathrm{\hspace{0.33em}3}A_G`$ and $`p\mathrm{\hspace{0.33em}1}/2`$. FIGURE 8 exhibits $`\mathrm{log}_{10}(t_J/t_D)`$ as a function of $`\mathrm{log}_{10}N`$ for the same integrations used to generate FIGURE 7.

There is also a clear visual sense in which, as $`N`$ increases, the frozen-$`N`$ orbits become progressively more regular in appearance. This is, e.g., evident in FIGURE 9, which exhibits the $`x`$-$`y`$ projections of representative frozen-$`N`$ orbits with $`N`$ varying between $`N=316`$ and $`N=316228`$, all generated from the same initial condition and integrated for a time $`t=25t_D`$ with $`ϵ=10^5`$. The final panel exhibits the smooth characteristic associated with the same initial condition. The most obvious point is that, as $`N`$ increases, the configuration space region to which the orbit is restricted more closely coincides with the region occupied by the characteristic. For example, only for the three largest values of $`N`$ is the orbit ‘centrophobic’ in the same sense as the characteristic.

Also evident is the fact that the orbit ‘looks smoother’ for larger values of $`N`$. This visual impression reflects the fact that, as $`N`$ increases, the power associated with the Fourier spectrum of an orbit tends to become more concentrated near a few special frequencies. (Since the smooth orbit associated with the same initial condition is regular, all its power is concentrated at a countable set of discrete frequencies.) This trend is illustrated in the eight panels of FIGURE 10, each of which exhibits $`|x(\omega )|`$ for a single frozen-$`N`$ orbit generated from the same initial condition. In each case, the data are so normalised that the peak frequency has $`|x(\omega )|=1`$. The spectra were generated from a time series of $`4001`$ points, recorded at intervals of $`0.025t_D`$.

The degree to which the orbits become more nearly regular with increasing $`N`$ can be quantified by determining the ‘complexity’ of the orbits, i.e., the number of frequencies in the discrete Fourier spectrum which contain an appreciable amount of power. Two such measures of complexity are illustrated in FIGURE 11, which was computed for ensembles of frozen-$`N`$ orbits with varying $`N`$, all evolved with $`ϵ=10^5`$ and generated from the same initial condition. The solid curve exhibits $`f_{0.1}`$, defined as the sum of the numbers of frequencies $`f_{0.1,x}`$, $`f_{0.1,y}`$, and $`f_{0.1,z}`$ which have more than 10% as much power as the peak frequencies for $`\omega _x`$, $`\omega _y`$, and $`\omega _z`$, i.e.,

$$f_{0.1}=f_{0.1,x}+f_{0.1,y}+f_{0.1,z}.$$
(18)
The dashed curve curve exhibits $`k_{0.95}`$, defined correspondingly as the sum of the numbers of frequencies required to capture 95% of the power in the $`x`$-, $`y`$-, and $`z`$-directions. In each case, the curve represents an average over different orbits in the ensemble, and the error bars represent the associated dispersions. The obvious point is that both these quantities decrease with increasing $`N`$.

The fact that, as $`N`$ increases, power becomes more concentrated near a few special frequencies has important implications for various physical processes which rely on resonances. For example, a variety of recent arguments in both galactic and solar system dynamics invoke a process of so-called of ‘resonant relaxation,’ which relies on the assumption that, in the presence of a large central object (a supermassive black hole in the center of a galaxy or the Sun at the center of the solar system), $`N`$-body orbits behave very nearly as if they were Keplerian trajectories in the fixed $`1/r`$ potential associated with the central object. If the chaos exhibited by individual orbits implied that these orbits were highly irregular, so that their power was not concentrated near the special Keplerian frequencies, resonant relaxation might seem quite implausible. Given, however, that the orbits become progressively more regular for increasing $`N`$, resonant relaxation would seem eminently reasonable, at least for systems in which $`N`$ is sufficiently large.

## V CONCLUSIONS AND DISCUSSION

Even though trajectories remain chaotic in the sense that the largest Lyapunov exponent does not decrease towards zero, there is a clear sense in which, for increasing $`N`$, orbits in frozen-$`N`$ potentials exhibit a pointwise convergence towards characteristics in the smooth potential which the frozen-$`N`$ potentials were chosen to sample. Viewed in configuration or velocity space, frozen-$`N`$ orbits tend to diverge linearly from the smooth characteristic with the same initial condition on a time scale $`t_G`$ that is proportional to $`N^{1/2}`$. In this sense, the continuum limit appears justified even at the level of individual trajectories. The fact that frozen-$`N`$ orbits remain chaotic for very large $`N`$ is completely consistent with the existence of a well-defined continuum limit.

It is easy to understand qualitatively why the frozen-$`N`$ orbits should remain chaotic even for very large $`N`$. Given that the chaos disappears completely in the continuum limit, where the orbits reduce to integrable characteristics, it would seem clear that the chaos must be associated with a sequence of ‘random’ interactions between a ‘test’ particle and a collection of ‘field’ particles. However, this would suggest that the time scale associated with the growth of a small initial perturbation can be estimated by considering the tidal effects associated with a pair of particles separated by a distance comparable to the typical interparticle separation. This tidal acceleration will of course scale as

$$\delta \ddot{𝐫}=(\delta 𝐫)𝐚\frac{Gm}{r^3}\delta 𝐫,$$
(19)
with $`r`$ the separation and $`m`$ the particle mass. Given, however, that $`rn^{1/3}N^{1/3}R_{sys}`$, with $`n`$ a characteristic number density and $`R_{sys}`$ the size of the system, it follows that the time scale $`t_{}`$ associated with the interaction should satisfy

$$t_{}\mathrm{\hspace{0.33em}1}/\sqrt{G\rho }.$$
(20)
In other words, the time scale associated with any orbital instability induced by the graininess of the system should be comparable to the dynamical time $`t_D`$, independent of particle number $`N`$. As $`N`$ increases, the size of the individual particle mass, $`m`$, and the cube of the typical separation between particles, $`n^1`$, both decrease as $`N^1`$ so that their ratio is independent of particle number.

The fact that the chaotic frozen-$`N`$ orbits appear to become “more nearly regular” as $`N`$ increases is consistent with recent claims that the “scale” associated with $`N`$-body chaos decreases with increasing $`N`$. Specifically, by comparing trajectories associated with two nearby initial conditions evolved in the same frozen-$`N`$ potential, Valluri and Merritt found that, when scaled in terms of $`R_{sys}`$, the size of the system, the typical separation $`R_{sat}`$ on which the initial exponential divergence saturates decreases with increasing particle number, so that $`R_{sat}/R_{sys}`$ is a decreasing function of $`N`$.

That the rate of divergence in the nonlinear regime slows more and more for larger $`N`$ can be quantified by tracking the actual evolution of two orbits generated from nearby initial conditions and determining the time required before their separation becomes ‘macroscopic.’ The result of such an investigation is illustrated in FIGURE 12, which was generated once again from ensembles of $`20`$ frozen-$`N`$ orbits all evolved with $`ϵ=10^5`$. In each case, the unperturbed orbits were identical to those used to generate FIGURE 1; the perturbed orbits involved changing the initial value of $`x`$ by an amount $`\delta x=10^6`$. FIGURE 12 exhibits as a function of $`N`$ the mean time $`\tau `$ required before the separation

$$\delta r=(\delta x^2+\delta y^2+\delta z^2)^{1/2}$$
(21)
had achieved the value $`\delta r=1`$. (For this initial condition, the average value of $`r`$ associated with the smooth characteristic was $`R_s\mathrm{\hspace{0.33em}1.83}`$.) The error bars were derived by considering the first and second ten orbits in the ensemble separately. Because individual orbits diverge at vastly different rates, the dispersion associated with a 20 orbit ensemble is much larger than reflected by these error bars. It is clear that $`\tau `$ increases systematically with increasing $`N`$, although considerably more slowly than with the $`N^{1/2}`$ dependence observed for the divergence time scales $`t_G`$ and $`t_J`$.

When viewed in terms of collisionless invariants like angular momentum, the divergence of frozen-$`N`$ orbits from smooth characteristics with the same initial condition is well approximated as a diffusion process, in which $`\mathrm{\Delta }J`$ grows as $`(t/t_J)^{1/2}`$ and where, for fixed $`t_D`$, the divergence time scale $`t_J`$ varies at least approximately as $`N^{1/2}`$. This reinforces the conventional wisdom that discreteness effects may be modeled as white, or nearly white, Gaussian noise in the context of a Langevin or Fokker-Planck description. It might, therefore, seem somewhat surprising that, although the divergence time scale $`t_G`$ in configuration or velocity space again scales as $`N^{1/2}`$, the quantities $`\mathrm{\Delta }r`$ and $`\mathrm{\Delta }v`$ grow linearly in time, rather than as $`t^{1/2}`$.

In this regard, it is significant that if the smooth Plummer potential be replaced by the smooth potential associated with a constant density configuration, the linear growth exhibited $`\mathrm{\Delta }r`$ and $`\mathrm{\Delta }v`$ is in fact replaced by the ‘expected’ diffusive behaviour. In this case $`\mathrm{\Delta }r`$ and $`\mathrm{\Delta }v`$ both grow as $`t^{1/2}`$, and, when expressed in units of the dynamical time $`t_D`$, the growth time $`t_G`$ is somewhat longer, corresponding more nearly to the time scale $`t_J`$ associated with $`\mathrm{\Delta }J`$.

This suggests strongly that the behaviour of $`\mathrm{\Delta }r`$ and $`\mathrm{\Delta }v`$ observed for the Plummer potential is associated with linear phase mixing. Because of finite number statistics, the same initial condition $`(𝐫_0,𝐯_0)`$ in different frozen-$`N`$ realisations of a Plummer potential will correspond to somewhat different energies, the values of which are conserved in the subsequent evolution. However, even neglecting discreteness effects, initially proximate orbits in a generic integrable potential will, if their energies be unequal, tend to diverge linearly. For example, two orbits evolved in a smooth Plummer potential with the same initial $`𝐫`$ but slightly different values of $`𝐯`$ and, hence, slightly different energies, will oscillate with somewhat different frequencies and, as a result, exhibit an overall linear divergence. If, however, the orbits are evolved instead in the potential associated with a constant density distribution, this is no longer true. A constant density sphere corresponds to a harmonic potential, where all orbits have the same unperturbed frequencies; and, for this reason, orbits in the smooth potential with slightly different energies will not exhibit such a systematic divergence.

The fact that frozen-$`N`$ orbits look “more nearly regular” for large $`N`$ suggests that the chaos associated with discreteness effects in the $`N`$-body problem should be viewed very differently from the chaos associated with a bulk nonintegrable potential. When evolved into the future, two nearby chaotic initial conditions in such a potential tend to diverge exponentially until they are separated by a distance comparable to the size of the easily accessible (i.e., not significantly impeded by the Arnold web) connected phase space region to which the orbits are confined, a region which tends, typically, to be macroscopic. By contrast, the scale associated with chaos induced by discreteness effects in the $`N`$-body problem is distinctly microscopic, at least for comparatively large $`N`$. It would appear that any single orbit with fixed energy can access a phase space region which is in fact very large; but the chaos which it experiences is a superposition of short range effects with characteristic scale $`R_{sys}`$.

###### Acknowledgements.

The authors acknowledge useful discussions with Alexei Fridman, Salman Habib, and Ilya Pogorelov. This research was supported in part by NSF AST-0070809 and by the Institute for Geophysics and Planetary Physics at Los Alamos National Laboratory.


format: text

text:
# Modeling two-state cooperativity in protein folding

## Abstract

A protein model with the pairwise interaction energies varying as local environment changes, i.e., including some kinds of collective effect between the contacts, is proposed. Lattice Monte Carlo simulations on the thermodynamical characteristics and free energy profile show a well-defined two-state behavior and cooperativity of folding for such a model. As a comparison, related simulations for the usual Gō model, where the interaction energies are independent of the local conformations, are also made. Our results indicate that the evolution of interactions during the folding process plays an important role in the two-state cooperativity in protein folding.

Due to the developments of experimental methods and theoretical models, many achievements of protein folding have been made recently. A protein can fold itself to its uniquely well-defined native structure in a biologically short time, regardless of the huge number of possible conformations, showing a highly cooperatively kinetic behavior. It is now clear that the cooperativity of folding may result from the backbone hydrogen bonding, sidechain packing and hydrophobic interactions, among them the hydrophobic interactions are believed to be the dominant driving force for folding. For many small single-domain proteins or lattice proteinlike models, there is a two-state behavior between the unfolded states and the folded native one. Recently, Chan and Kaya indicated that according to the calorimetric criterion, which is widely used in experiments as a condition for two-state folding, popular lattice models, e.g., two-letter HP and 20-letter MJ models, are far from two-state models. This may be due to some flawed assumptions in the potential functions used in these models. Lattice models usually use statistical potential functions extracted from the pairing frequencies of 20 kinds of amino acids in databases of protein structures. Although these knowledge-based potentials may be a good approximation to the relative strength of interactions between the residues in the native state, they provide no information about how the interactions evolve during the folding. For computational convenience, a common assumption in lattice models is that the interactions are additive, and they are the same during the folding as in the native state. This means that the interaction energies are conformation-independent. Clearly this is not relevant to the experimental situation. In fact, as Dill pointed out, the thermodynamic additivity principle which is widely used in chemistry may be unsuitable in biochemistry. Some recent experiments also indicated that the transition state is an expanded version of the native state, where the majority of interactions are partially formed, and their strengths are different from those in the native state (with $`\mathrm{\Phi }<1`$). That is, these interactions depend on the conformations, especially the local structures around the contacts as emphasized recently in Ref. . Previously, the non-additivity was built in a lattice model for packing effects; and the hydrophobic force depending on the local density of peptide atoms was also taken into account in an off-lattice model . Studies on these models show that the introduction of the non-additivity is significiant, but the two-state cooperativity of these models is not checked, and the effects of the non-additivity on thermodynamics and kinetics of folding need to be further studied.

In this paper, we develop a refined Gō model where the pairwise interaction energies vary as the local environment changes, i.e., some kinds of collective effects between contacts are introduced. Our purpose here is to study the two-state cooperativity of protein folding and its physical origin with such a model. Our results give a general picture about how the conformation-dependent interactions affect the folding kinetics, which is consistent with the phenomenological explanation based on experimental results.

We model a polypeptide chain as a self-avoiding chain on a cubic lattice. A contact is formed if two residues are space adjacent but not sequence adjacent. If two residues form a contact as the same as in the native state, we call this contact a native contact, otherwise a non-native contact. Following the Gō model, only native contacts are considered to contribute to the total energy. Different from the Gō model, we assume that the interaction energies between residues are conformation-dependent, and vary with changes of the local environment. To achieve this, we introduce a parameter S to describe the degree for a residue being ordered relative to the native state. For the $`i`$-th residue in a certain conformation, its degree of order $`S_i`$ is defined as

$$S_i=z_i/z_i^{nat},$$
(1)
where $`z_i`$ is the number of native contacts in this conformation, $`z_i^{nat}`$ is the number of contacts formed in the native state. Obviously, $`S_i`$ varies between $`0`$ (the $`i`$-th residue being fully disordered) and $`1`$ (being fully ordered). Thus, the interaction energy between residues $`i`$ and $`j`$, $`B_{ij}`$=-($`S_i`$+$`S_j`$)$`\epsilon `$/2 is defined, where $`S_i`$ and $`S_j`$ are the degrees of order for residues $`i`$ and $`j`$, respectively. $`\epsilon `$ is the unit of energy and is set to be $`1`$ in this work. The total energy of the conformation then is $`E=_{i<j}\mathrm{\Delta }_{ij}B_{ij}`$, where $`\mathrm{\Delta }`$ is unity when residues $`i`$ and $`j`$ form a native contact, and zero otherwise. Here, a contact formed between residues $`i`$ and $`j`$ may have different energies in different conformations, i.e., $`B_{ij}`$ may change from one conformation to another (for the Gō model, one always has $`B_{ij}=\epsilon `$). In general, a contact formed between residues $`i`$ and $`j`$ will stabilize, to some extent, other contacts that residue $`i`$ or $`j`$ formed with other residues. On the contrary, its breakage may destabilize those contacts as well. Therefore, the introduction of the degree of order for a residue into the potential function reflects the cooperativity between the residues. Although the correlation distance is small, only one lattice unit, the many-body effects are obviously included in our model. Figure 1 shows such a collective effect. The interaction energies of contact A-B (or B-C) are different when the other contact is present or not present. Clearly, the energy of state $`I_3`$ is lower than the sum of that of states $`I_1`$ and $`I_2`$, indicating the interaction non-additivity. Each contact is stabilized by the other contact due to the collective effect. Note that in this paper our model is called Gō+ model to distinguish from the Gō model.

Now let us present the Monte Carlo simulations on the thermodynamic and kinetic features for both models. The mean first passage time (MFPT), as a common measure of folding rate, is calculated by an average of the first passage time (FPT) over 1000 runs. Each run begins with a random conformation, and ends when the native state is reached for the first time. The FPT is the Monte Carlo steps (MCS) consumed in a run.

Generally, as the temperature $`T`$ decreases, the population of the native state, $`P_N`$, increases from zero to about unity. The degree of sharpness of changes in $`P_N`$, similar to the “rapidity” in Ref., is a measure of the cooperativity of the folding reaction. Figure 2 shows the population $`P_N`$ and the specific heat $`C_v`$ versus temperature $`T`$ for a $`36`$-mer chain for both models. $`P_N`$ is defined as $`P_N=e^{E_N/T}/_E\mathrm{\Omega }(E)e^{E/T}`$, where $`\mathrm{\Omega }(E)`$ is the density of states for energy $`E`$, $`E_N`$ is the energy of the native state. $`\mathrm{\Omega }(E)`$ is calculated with the Monte Carlo histogram method. From Fig.2 we can see that the folding transition for our Gō+ model is much sharper than that of the Gō model, i.e., a sharper change in $`P_N`$. There is also a single peak in $`C_v`$ curve, but it is narrower than that of the Gō model. For our Gō+ model, the maximum of $`C_v`$ occurs at a temperature nearly the midpoint temperature of transition with $`P_N`$=1/2, i.e., the difference between these two temperatures is quite small. This is consistent with recent studies on naturally occuring proteins, implying a good cooperativity of folding in this model. Differently, such a temperature difference is large for the Gō model (see Fig.2), indicating that the folding of the Gō model is much less cooperative than that of the Gō+ model. Since the sharpness is only a qualitative description for the transition, we further calculate the equilibrium energy distribution at the folding transition temperature, $`T_f`$. Figure 3 shows such distributions for both models. Clearly our Gō+ model shows a good bimodal behavior, and the denatured-state energy is distributed in a narrow region \[see Fig. 3(a)\]. This means clearly a two-state folding and there is basically no intermediate states at equilibrium. Differently, for the Gō model as shown in Fig. 3(b), there are many intermediate states and the bimodal behavior is not so significant as that in Fig.3(a). Thus for the Gō model the folding is not of a two-state. This is in agreement with Chan and Kaya’s argument.

In experiments, a well-established criterion for two-state folding is that the van’t Hoff enthalpy $`\mathrm{\Delta }H_{vH}`$ around the transition midpoint is equal, or very close, to the calorimetric enthalpy $`\mathrm{\Delta }H_{cal}`$ of the entire transition. In this work, we calculate the ratio $`\mathrm{\Delta }H_{vH}`$/$`\mathrm{\Delta }H_{cal}`$ as suggested in Ref. (here, the definition of $`\mathrm{\Delta }H_{vH}`$/$`\mathrm{\Delta }H_{cal}`$ is equal to $`(k_2)^2`$ in Ref.), and list the results in Table I. From Table I, we can clearly see the difference between the Gō model and our Gō+ model. The Gō model, which is considered as a model with minimal energetic frustrations, does not meet the calorimetric two-state criterion and gives out the value of $`\mathrm{\Delta }H_{vH}`$/$`\mathrm{\Delta }H_{cal}`$ far from 1. Nevertheless, our model satisfies the criterion quite well (for real proteins, the value of $`\mathrm{\Delta }H_{vH}`$/$`\mathrm{\Delta }H_{cal}`$ is 0.96$`\pm 0.03`$ ). This, again, implies the two-state folding and the good cooperativity of our Gō+ model.

Physically, the high cooperativity of our model may result from the narrow distribution of the denatured states and the high population of the native state at the folding temperature (see also Figs.2 and 3). In our model, the energy spectrum relating to various conformations is redistributed, comparing with that of the Gō model, due to the collective effect between interactions. As a result, the energies of non-native conformations are moved to higher energy levels and a larger energy gap is left between the non-native conformations and the native one (for the two models, the energies of the native state are the same). The large energy gap makes the native state paticularly stable, which is believed to be a necessary condition for cooperative folding. This may be the physical origin of the two-state cooperativity. It can be further explained from the viewpoint of the free energy profile. For our Gō+ model, as shown in Fig.4, the free energy profiles have broad activation barriers. The broad activation barriers can account for the large movement of transition state caused by mutation or temperature changes, and are considered as a common feature of the two-state folding. Our numerical results are surprisingly consistent with a phenomenological speculation for the existence of such a free energy profile in Ref. . It should be noted that the broad activation barriers are consistent with the narrow distribution of the denatured states.

Now let us make a comparison of the foldability based on the plots of the MFPT versus $`P_N`$ for both models. Note that we use $`P_N`$ instead of the commonly used temperature $`T`$ in the horizontal axis in Fig.5. This is because that an identical condition should be taken for the comparison. In lattice simulations, the temperature has an arbitrary unit and also has no direct relationship with the real temperature. The comparison between two different models at the same temperature may make no sense. Nevertheless, at an identical condition with the same $`P_N`$, the differences in the foldability can be well-defined. This is similar to other conditions used previously. From Fig. 5, we can see that the MFPT for our Gō+ model shows a slow decrease as $`P_N`$ increases, it reaches a minimum at $`P_N0.93`$, and then it increases. For the Gō model, there is also a minimum but at $`P_N0.71`$. It is clearly that when the native state is stable (say, $`P_N`$ 0.9), our Gō+ model folds significantly fast, i.e., the MFPT is smaller with one or two orders of magnitude than that of the Gō model. Physically, this can be explained as follows. From Eq.(1) we can easily see that the energy gain of forming a contact is usually smaller for our Gō+ model than that for the Gō model. At high temperatures, entropic contribution is dominant to the free energy barrier, and the loss of entropy is always undercompensated by the energy gain, thus the Gō+ model folds slower for its smaller energy gain. Whereas at low temperatures, folding is nearly a downhill process, and the loss of entropy is always overcompensated by the energy gain. Therefore, for the Gō+ model, it is easier to escape from kinetic traps, and the folding is faster. Finally, we note that for the two models the pathways of reaching the transition state from the denatured state are different. Due to the high cooperativity in our Gō+ model, a good core, the assembly of non-polar residues, is formed much earlier at low temperatures than that in the Gō model. Detailed kinetic results will be reported elsewhere. We also note that similar results are obtained for different chain sizes.

In conclusion, our Gō+ model, with many-body interactions depending on the local structures included, exhibits a good two-state folding behavior. Our results suggest that the evolution of interactions during the folding plays an important role in the two-state cooperativity in protein folding. We give a possible way how the interactions evolve in the folding, which may capture some essential features of the two-state folding. We expect further study could provide new insights into the mechanism of protein folding.

We thank H.S. Chan, A. Maritan and D. Thirumalai for useful suggestions. This work was supported by the Foundation of NNSF (No.19625409, and No.10074030).

$``$Email address: wangwei@nju.edu.cn

Table I: The ratios of $`\mathrm{\Delta }H_{vH}`$ / $`\mathrm{\Delta }H_{cal}`$ for the Gō model and our Gō+ model, respectively. Ten sequences are calculated for each chain size.

FIG. 1: Schematic illustration of collective effect between two interactions. From a state $`I_0`$ with three unstructured residues, the chain can be settled in a state $`I_1`$ (or $`I_2`$) with a contact A-B (or B-C) and an equilibrium constant $`K_1`$ (or $`K_2`$). A state $`I_3`$ with two contacts A-B and B-C can be reached from state $`I_1`$ or $`I_2`$, but with different equilibrium constants $`K_2\gamma `$ or $`K_1\gamma `$. In state $`I_3`$, each interaction is stronger by a factor $`\gamma `$ due to the existence of the other contact.

FIG. 2: Population $`P_N`$ and specific heat $`C_v`$ varying with the temperature $`T`$ for a $`36`$-mer chain.

FIG. 3: The energy distribution for the same 36-mer used in Fig.2, using (a) Gō+ potential and (b) Gō potential at respective folding transition temperature, $`T_f`$.

FIG. 4: The free energy profile $`F(E)=ETS(E)`$ of our Gō+ model at different temperatures, where entropy $`S(E)`$ is calculated by using entropy sampling Monte Carlo method. Here U, N and TS denote the unfolded state, native state and transition state, respectively. Note that the free energy profile at high temperature is overall shifted so that the unfolded states are overlapped.

FIG. 5: MFPT versus $`P_N`$ for a 36-mer chain.


text:
# Simple technique for superconducting joints quality estimation in bulk melt-processed high temperature superconductors

## Abstract

We propose an empirical approach to estimate the quality of superconducting joints (welds) between blocks of bulk high temperature superconductors (HTS). As a measuring value, we introduce a joint’s quality factor and show its natural correlation with joint’s critical current density. Being simple and non-destructive, this approach is considered to be quite important to solve the problem of utilization of HTS in large scale applications. The approach has been applied to characterize the joint’s quality of melt-processed Y-123 joined by Tm-123 solder.

Bulk melt-processed high temperature superconductors (MP HTS) have been appeared to be quite suitable for such large scale HTS applications as superconducting motors, contactless bearings, flywheels for energy storage, levitation transport. Remarkable values of critical current densities have been achieved for these quasi-single crystals and the main problem now, which is an obstacle in the way of their practical utilization, is superconducting joining of them into appropriate for this utilization blocks. While the problem being under active investigation  has not been completely resolved yet, it is vitally important to find a simple method for quantitative characterization of superconducting joints.

The most physical way to describe a superconducting joint quantitatively is to determine density of critical current that flows through it. While the transport measurements are quite difficult to exploit here, the contactless techniques appear to be the most attractive. It has been shown that the levitation force measurements, magneto-optical imaging techniques, or Hall-probe magnetometry can be successfully used to characterize superconducting joint’s quality. At last, in a very recent paper, an approach to estimate critical current density through superconducting joint for ring samples was proposed.

In this Letter we propose a simple non-destructive contactless method to evaluate the quality of superconducting joint and estimate critical current density through it. The method is based on the technique of critical current density determination from local levitation force measurements which, as has been shown in Ref. , allows to determine critical current density in a thin undersurface layer.

We test the method on Y-123 MP HTS samples with superconducting joints made by solidification technique (”soldered joints”) with Tm-123 as a solder. The details of the technique are described in Ref. . Two types of experimental procedure were used: ’three times procedure’ and ’four points procedure’. Within ’three times procedure’ we measure the levitation force that acts on a small permanent magnet (PM) three times at one point: above the uniform sample before cutting, $`F_{un}`$, above the ’crack’ (mechanical joint) after cutting, $`F_{cr}`$, and above the superconducting joint after soldering procedure, $`F_{jn}`$. ’Four points procedure’ (see Fig.1), being less precise due to samples’ inhomogeneity, was used for the samples which were not measured before cutting and joining, and $`F_{un}`$ was estimated as $`(F_1+F_2)/2`$. We used two PMs: spherical one, of 1.5 mm in diameter with magnetic moment $`\mu `$ = 1.9 G cm3, and cylindrical one, of 6.3 mm in diameter and 2.3 mm in width with $`\mu `$ = 38 G cm3. The readout distances between centers of the magnets and HTS surface were 1.25 mm and 6 mm respectively. The basic experimental setup and the technique for critical current density determination have been described in details in Ref. .

Based on such local levitation force measurements, one can introduce a joint’s quality factor

$$q=\frac{F_{jn}F_{cr}}{F_{un}F_{cr}}.$$
(1)
This formula was constructed to satisfy the natural asymptotic conditions: $`q0`$ when $`F_{jn}F_{cr}`$ (non-superconducting joint), and $`q1`$ when $`F_{jn}F_{un}`$ (an ideal superconducting joint). In the following we will find the physical meaning of thus introduced quality factor through its relation with joint’s critical current density $`J_{jn}`$.

Let us consider the next naive picture of circular currents which are sketched at Fig.2. In a uniform sample a circular current $`I_0`$ flows mainly over the top surface of the sample  and one can write $`F_{un}=\alpha I_0`$. If we crack the sample, the current (having the same value, until the size of the crack is much less than sample dimensions) cannot flow though the crack and has to close over the crack walls but one can consider this as an additional opposite ($`I_1=I_0`$) circular current flowing through the crack and over the crack walls: $`F_{cr}=\alpha I_0+\beta I_1=(\alpha \beta )I_0`$. Next, it is quite natural to assume that for a sample with a superconducting joint one can estimate the real current through the joint as

$$I_{jn}=I_0I_1=\frac{J_{jn}}{J_c}I_0,$$
(2)
where $`J_c`$ is the bulk critical current density. Then $`I_1=(1(J_{jn}/J_c))I_0`$ and $`F_{jn}=(\alpha \beta (1(J_{jn}/J_c)))I_0`$. Substituting this into (1) we obtain

$$q=\frac{J_{jn}}{J_c}.$$
(3)
Thus, the superconducting joint’s quality factor introduced by (1) represents the ratio of density of the critical current that flows through the joint to the bulk critical current density of the sample. The Eq.1 by itself shows a simple way to estimate this quality factor but in combination with the method of $`J_c`$ determination from levitation force measurements  gives a simple non-destructive method for determination of density of the critical current through the superconducting joint.

Though the given derivation of Eq.3 from Eq.1 is quite simple, it has some restrictions which should be considered. The first of them is a validity of consideration of levitation forces to be proportional to corresponding currents - we have assumed above that the introduced coefficients $`\alpha `$ and $`\beta `$ do not depend on $`I_0`$ and $`I_1`$, correspondingly. The first assumption, $`\alpha (I_0)=const`$, works when the depth $`\delta `$ of penetrated magnetic field at sample surface ($`\delta =cB_r/4\pi J_c`$ in CGS units, here $`c`$ is the velocity of light and $`B_r`$ is the tangential magnetic field at HTS surface which is twice bigger than PM’s field in ”zero approximation”) is much less than system dimensions $`L`$ (the distance from center of the magnet to sample surface in our geometry), which is the case in the experiment being considered and can always be achieved by choosing appropriate values of $`L`$ or $`\mu `$. Moreover, it was shown in Ref. that the condition $`\delta L`$ is too strong and the assumption works well even if $`\delta L`$. To fulfill the assumption that $`\beta (I_1)=const`$ is not so easy. In this case, the penetration depth $`\delta _{jn}`$ of intrajoint field must be less then joint’s width that impose strong restrictions on magnetic field value. For the samples which we investigated the joints were $`d`$ 100 $`\mu `$m in width that has motivated the readout distances we used to ensure $`\delta _{jn}<d`$. So, here we measure the superconducting joint’s quality for magnetic field 102–103 G but the method can be extended to higher fields by replacing levitation force measurements to resonance oscillations technique.

Another point, which should be mentioned, relates to sample inhomogeneity. Although every method of critical current density determination encounters the same problem a reasonable approximation here is to assume the validity of Eq.3 as integral relation over a region where the main amount of the induced current flows ($`2L`$ in our case). This sets a restriction on spatial resolution of the method but also gives a possibility to vary the integrated area using different magnets. For two magnets we used the integrated areas were about 2 mm and 10 mm in diameter respectively and 0.1 mm in depth.

The results of our measurements are summarized in the Table 1. The samples were treated at temperature $`T_{tr}`$ for 0.17 h in air (1,2) or for 0.5 h in oxygen (3-6). As a solder, the powder of Tm-123 was used for even samples and Tm-123 with addition of 10 wt% of Y-211 (’green phase’) for odd samples. The values $`q_I`$, $`q_{II}`$ and $`q_M`$ were determined with first and second magnets and from magnetic flux mapping technique respectively. The joint’s quality value from mapping technique was estimated as $`q_M=2B_{ju}/(B_1+B_2)`$, where $`B_{ju}`$, $`B_1`$, and $`B_2`$ are the local maximum trapped magnetic fields in junction and in two pieces of cut samples. The values of $`J_c`$ (determined from $`F_{un}`$, see Ref.) and $`J_{jn}`$ are also represented in the Table.

Although the mapping technique gives rather rough values it is possible to make some conclusions. First, the $`q_{II}`$ values are closer to $`q_M`$ than $`q_I`$ because they are integrated over lager volume. Second, we have compared the average values of quality factors for different series of measurements (here $`qq_I`$): over all samples, $`q=0.57`$, $`q_M=0.48`$; over all in-plane surfaces (parallel to $`ab`$-plane), $`q_{}=0.43`$, $`q_M_{}=0.50`$; and over all out-plane surfaces (perpendicular to $`ab`$-plane), $`q_{}=0.70`$, $`q_M_{}=0.47`$. Assuming that the quality factors determined from levitation force measurements give the information from a thin undersurface layer (undersurface quality) and the factors determined from mapping represent an integrated values over whole volume (bulk quality), it is possible to say, that the undersurface joint’s quality is slightly better than the bulk joint’s quality and that the undersurface quality is much better for out-plane surfaces in comparison to in-plane ones. In other words this means that in the tested samples: (i) the average quality of joints is slightly better at the edges than in the middle; (ii) the critical current density of the material which fills the joints is less anisotropic than the bulk critical current density of the samples itself.

In summary, we have proposed a simple empirical approach to estimate the quality of superconducting joints between blocks of bulk high temperature superconductors. As a measuring value, we introduced a joint’s quality factor and found its natural correlation with joint’s critical current density. Simultaneous use of the proposed method and flux mapping technique has allowed us to estimate space inhomogeneity of superconducting joints and even obtain information about anisotropy of the joint’s critical current density.


text:
# Superconducting properties of MgB2 probed by radiation-induced disordering

## Abstract

The electrical resistivity $`\rho `$ and upper critical field $`H_{c2}`$ were measured in MgB2 disordered by nuclear reactor neutrons. We found that superconducting temperature $`T_c`$ decreases under irradiation from 40 to 5 K. Despite strong disordering (more than 10 dpa (displacements per atom)), the initial crystal structure is preserved. The residual resistivity $`\rho _0`$ increases from 0.35 to 2 m$`\mathrm{\Omega }`$ cm while $`dH_{c2}/dT`$ remains approximately unchanged ($`dH_{c2}/dT0.5`$ T/K), and such behavior may be interpreted as decrease in the density of electronic states.

With the discovery of superconductivity (SC) in binary compound MgB2 with relatively high $`T_c40`$ K Akimitsu , a question arises about the possible origin of SC in such a simple compound. It is well known that high density of electronic states at the Fermi level $`N(E_\text{F})`$, similar to soft phonon frequencies, can lead to enhancement of $`T_c`$ in BCS phonon-mediate superconductivity McMillan , while some form of exotic coupling can take place in this case. One of the ways to determine the nature of coupling is to study the response of SC on disordering induced by irradiation with high-energy particles.

Disordering of intermetallic compounds, such as A15 (Nb3Sn, V3Si), with high $`N(E_\text{F})`$ leads to decrease of $`T_c`$ from 15 - 20 K to 1 - 3 K due to drop in $`N(E_\text{F})`$ at disordering Karkin . In low-$`T_c`$ compounds (Mo3Si and Mo3Ge), irradiation leads to increase of $`T_c`$ (from 1.5 to 7 K) due to rise in $`N(E_\text{F})`$ as well as due to softening of phonon frequencies. It means that individual features of the electronic structure disappear under disordering, and the superconducting properties of the disordered compounds become similar to the properties of amorphous materials obtained by various methods.

The influence of disordering on SC of high-$`T_c`$ compounds is not similar to that of A15: irradiation leads to very fast degradation of SC ($`T_c=0`$) Aleksashin . Although the origin of $`T_c`$-degradation is not clear (because the mechanism of SC in high-$`T_c`$ is unknown), there is a significant difference between intermetallic compounds and more exotic materials in the response to disordering.

In this paper we studied the transport properties of polycrystalline samples ($`0.0515`$ mm3) irradiated in nuclear reactor IVV-2M at $`T350`$ K up to fluences of thermal neutrons $`\mathrm{\Phi }=110^{19}`$ cm-2, and of fast neutrons, to $`510^{18}`$ cm-2 (to the total dose - of more than 10 dpa). The MgB2 powder was prepared by the technology described in Ref. Gerashenko . The samples were compacted under 9 Gpa at room temperature without any subsequent heat treatment. Such sample preparation produces a stressed structure, which results in broadening of X-ray peaks (Fig. 1) and resistivity $`\rho `$ (Fig. 2) and ac-susceptibility $`\chi `$ (Fig. 3) superconducting transitions. However the initial crystal structure (C32) is preserved under irradiation inducing such high damage to the material.

Fig. 1 shows X-ray diffraction patterns (Cu-Kα) for the initial and the irradiated samples of MgB2. The samples contain traces ($``$ 3 %) of MgO. The structure parameters were refined with Rietveld analysis. Some discrepancy between the measured and the calculated curves (see the insert in Fig. 3 ) was connected with the sample texture. Irradiation leads to anisotropic expansion of the crystal lattice (the cell volume increase is about 1.4 %) with faster increase of c, as is well seen in the insert in Fig. 3 and Table 1. The broadened diffraction peaks in the initial samples become narrower after irradiation, which is an evidence of partial relief of stress. The Rietveld analysis refinement yields some decrease in the Mg site occupation numbers (if B site occupation is fixed equal to unit), although we do not expect appearance of anti-structural defects because of the difference in the atomic sizes of Mg and B.

Figs. 2, 3 show transformations of $`\rho `$\- and $`\chi `$-transitions in irradiated and annealed samples. There are no essential changes of transition widths under irradiation, which testifies to homogeneous distribution of radiation defects. Annealing at 500 K recovers 85 % of $`T_c`$. The values of the upper critical field $`H_{c2}`$ determined by half-width of the normal state resistivity show almost the same slope for the initial, irradiated and post-annealed samples (Fig. 4). $`H^{}=dH_{c2}/dT`$ shows low dependence on normal state resistivity $`\rho _0`$ (see the insert in Fig. 4). The resistivity $`\rho (T)`$ increases after irradiation, decreases after annealing up to 300 C and grows again at subsequent annealing reaching 4 m$`\mathrm{\Omega }`$ cm (Fig. 5). We consider that such a behavior of $`\rho (T)`$ may be partly due to intergrain transport, therefore $`\rho (T)`$ is not completely an intrinsic property. As is expected, $`H^{}`$ within the dirty limit should increase with the decrease in the mean free path of electron under disordering, therefore the small change in $`H^{}`$ may be interpreted as the decrease of $`N(E_\text{F})`$ with decreasing $`T_c`$. It may be supposed that, similar to A15 compounds, there is a fine structure in $`N(E)`$ in the vicinity of the Fermi level, which is smeared under radiation-induced disordering.

Work supported by the Russian State contracts Nos 107-1(00)-P, 107-19(00)-P, 108-31(00)-P, the Russian State Program of Support to Leading Scientific Schools (Project No. 00-15-96581) and RFBR Grant No. 00-02-16877.


text:
# Atom laser dynamics

## Abstract

An ideal atom laser would produce an atomic beam with highly stable flux and energy. In practice the stability is likely to be limited by technical noise and nonlinear dynamical effects. We investigate the dynamics of an atom laser using a comprehensive one dimensional, mean-field numerical model. We fully model the output beam and experimentally important physics such as three-body recombination. We find that at high pump rates the latter plays a role in suppressing the high frequency dynamics, which would otherwise limit the stability of the output beam.

Optical lasers have had an enormous impact on science and technology, due to the intensity and coherence of the light they produce. We present a theoretical study of the output beam of the analogous matter wave device, known as the “atom laser”. These have been demonstrated in a number of laboratories, although so far only in an unpumped mode .

Theoretical studies have highlighted both the similarities and the differences between optical and atom lasers . The differences arise because atoms are more complex than photons: they have mass, giving a different free space dispersion relation, and they interact with each other, producing strong nonlinearities. The latter generates complex dynamics, potentially complicating certain experimental measurements, of quantum noise for example. Nevertheless, the output beam dynamics might be useful, in itself, as a probe of the excitations of the trapped condensate .

Fundamentally, it is the coherence properties, such as linewidth and intensity correlations, that best capture the unique physics of lasers . Coherence is also important for many practical applications, such as interferometry. A classical nonlinear model of the laser as a noise driven van der Pol oscillator shows that the noise power, in both the phase and the amplitude, decreases in inverse proportion to the laser power .

Schawlow and Townes showed that optical lasers’ first order coherence, or linewidth, is ultimately limited by spontaneous emission, and also decreases in inverse proportion to the laser power . Analogous results have been derived for atom lasers . In practice, however, the Schawlow-Townes limit is not achieved, because the linewidth is limited by technical noise and by dynamical effects rather than by quantum noise. For example relaxation oscillations, due to the nonlinear interaction between the inversion and the light, are the primary determinant of the low frequency noise spectrum in many optical lasers .

In this paper we focus on one aspect of the nonlinear dynamics of atom lasers: the frequency spectrum of the atom laser output beam density. In particular, we investigate how the spectrum depends on the pumping rate. This spectrum describes the atom laser nonlinear dynamics. An experimentally measured spectrum would also include genuine “noise” due to physical processes that we do not model, such as fluctuations in the trapping potential. This investigation is motivated by the inverse relation between noise power and pumping rate in the optical laser. We choose to focus on the density because it can be measured experimentally using sensitive optical techniques such as phase modulation spectroscopy .

A conclusion of our work is that three-body recombination plays a major role in the high frequency dynamics. As the pumping rate is increased the spectral peaks weaken and move to higher frequencies. This knowledge will facilitate atom laser applications, and measurements of their quantum noise.

Our atom laser model is schematically illustrated in Fig. 1. It consists of a dilute gas Bose-Einstein condensate trapped by a highly anisotropic, effectively one-dimensional, potential. This is pumped from a reservoir of non-condensed trapped atoms. The output atom laser beam is generated by a Raman transition to an untrapped atomic state. Our model builds on that of Kneer et al. by adding three-body recombination and a propagating output beam. In the following we will refer to the trapped Bose-Einstein condensate as the “condensate”, and to the output beam as the “beam”. We make the mean-field approximation in which the condensate and beam are described by the classical fields $`\mathrm{\Psi }_a(x,t)`$ and $`\mathrm{\Psi }_b(x,t)`$ respectively, which obey Gross-Pitaevskii type equations . Our model is more complete than previous mean-field treatments  since we include pumping, output coupling, three-body recombination, and explicitly model the dynamics of the pump reservoir, the condensate, and the beam. Although the extension to three dimensions is straightforward, we do not report it here. To describe the pumping we use a phenomenological model, which mimics the pumping of an optical laser . It only depends on the total number of atoms in the uncondensed component, $`N_u(t)`$, not on its spatial structure. The condensate atoms are coupled, by a reversible Raman transition , to an untrapped electronic state forming the atom laser beam. The Raman transition imparts a momentum kick $`\mathrm{}k`$ to the out-coupled atoms. The beam evolves under the influences of gravity and atom-atom interactions, which are dominated by those due to the condensate. We have only considered the particular case for which the Raman coupling is tuned to transfer atoms starting from the center of the harmonic trap . Hence the output beam overlaps with the lower half of the condensate. Experimentally, three-body recombination is well described as a loss process occurring at a rate proportional to the cube of the local density . We follow Kagan et al.  in incorporating three-body recombination into the atom laser Gross-Pitaevskii equations.

Our model is defined by the following dimensionless equations:

$$\begin{array}{c}\frac{dN_u}{dt}=r\gamma _uN_u\gamma _pN_uN_a,\hfill \\ i\frac{\mathrm{\Psi }_a}{t}=\frac{1}{2}\frac{^2\mathrm{\Psi }_a}{x^2}+\frac{1}{2}x^2\mathrm{\Psi }_a+U_a\mathrm{\Psi }_a|\mathrm{\Psi }_a|^2\hfill \\ +U_{ab}\mathrm{\Psi }_a|\mathrm{\Psi }_b|^2i\gamma _r\mathrm{\Psi }_a(|\mathrm{\Psi }_a|^4+|\mathrm{\Psi }_b|^4)\hfill \\ +\gamma _Re^{ikx}\mathrm{\Psi }_b+\frac{i}{2}\gamma _pN_u\mathrm{\Psi }_a,\hfill \\ i\frac{\mathrm{\Psi }_b}{t}=\frac{1}{2}\frac{^2\mathrm{\Psi }_b}{x^2}+Gx\mathrm{\Psi }_b+U_b\mathrm{\Psi }_b|\mathrm{\Psi }_b|^2\hfill \\ +U_{ab}\mathrm{\Psi }_b|\mathrm{\Psi }_a|^2i\gamma _r\mathrm{\Psi }_b(|\mathrm{\Psi }_a|^4+|\mathrm{\Psi }_b|^4)\hfill \\ +\gamma _Re^{ikx}\mathrm{\Psi }_a.\hfill \end{array}$$
(1)
The model is made dimensionless using the characteristic trap length $`l=(\mathrm{}/\omega m)^{1/2}`$ and angular frequency $`\omega `$, with $`m`$ the atomic mass. Hence time $`t`$, position $`x`$, and the fields are measured in units of $`\omega ^1`$, $`l`$ and $`l^{1/2}`$, respectively. Experimentally reasonable values of the parameters are: a trap frequency $`\omega 125`$ Hz , and the atomic mass of sodium $`m=3.8\times 10^{26}`$ kg. These give a time scale of $`\omega ^1=8`$ ms, and a length scale of $`l=4.7`$ $`\mu `$m.

Atoms are injected into the uncondensed fraction at the rate $`r`$, are lost at the rate $`\gamma _uN_u`$, and pump the condensate at the rate $`\gamma _pN_uN_a`$, where $`N_a=_{\mathrm{}}^{\mathrm{}}|\mathrm{\Psi }_a|^2𝑑x`$ is the total condensate population. $`U_a`$ and $`U_b`$ are the intra- and $`U_{ab}`$ the inter- species two-body interaction coefficients. The dimensionally correct coefficients can be written as $`U_{a,b,ab}=4\pi a_s/(\mathrm{}/\omega m)^{1/2}`$ with $`a_s`$ the appropriate s-wave scattering lengths . We have previously found that three-body recombination, with coefficient $`\gamma _r`$, is necessary for the system to reach a quasi-stationary-state on a timescale comparable to experimental condensate lifetimes . The three-body recombination rate $`\gamma _r=10^7`$ is chosen to produce $`N_a35009500`$ condensed atoms at steady state, as the pump rate varies from $`r=200`$ to $`1600`$. For simplicity we assume that the inter- and intra-species three-body recombination rates are the same. $`G`$ is the scaled gravitational acceleration . Gravity has not been explicitly included for the trapped atoms, since it is equivalent to a spatial shift of the trap potential. The Raman coefficient $`\gamma _R`$ is a function of the Raman laser amplitude and detuning .

Our numerical method is a split-step Fourier method  with a fourth order Runge-Kutta in place of the usual nonlinear step. We have absorbing boundaries at the ends of the spatial grids to prevent reflections.

Typical results for quasi-stationary-state spatial profiles of the condensate and beam densities are shown in Fig. 2. For our parameters, the condensate shape is well described in the Thomas-Fermi approximation. Note that the plotted output beam density has been scaled up by a factor of 100. The flux in the beam is about 90 atoms per unit dimensionless time, or about $`1.1\times 10^4`$ atoms per second with $`\omega =125`$ Hz. The prominent spatial oscillations of the output beam, occurring within the extent of the condensate, have been observed in previous work  and are related to the oscillatory Airy type eigenfunctions of particles in the linear gravitational potential.

In Fig. 3 we show the output beam density $`|\mathrm{\Psi }_b(x_0,t)|^2`$ as a function of time at a fixed position far from the trapped condensate. The time span has been chosen to be longer than that of the transient dynamics associated with the growth of the condensate from a seed, and short enough to be experimentally accessible. Quasi-stationary, nearly periodic dynamics develop, becoming approximately harmonic for large pump rates $`r`$. We have verified that the quasi-stationary dynamics is independent of the particular form of the seed used to initiate the condensate growth. The fundamental oscillation of the beam, at higher pump rates, arises from the well known Kohn mode , corresponding to a rigid “sloshing” back and forth of the condensate at the trap frequency.

Fig. 4 shows the frequency spectra of the density time series in Fig. 3, and is the central result of this paper. The spectra are those of the final, quasi-steady, parts of the time series. The frequencies are well approximated by the eigen-frequencies of an uncoupled 1D condensate. Kneer et al.  found these to be $`\mathrm{\Omega }=\omega \sqrt{n(n+1)/2}`$ where $`n=1,2,3,\mathrm{}`$, with the corresponding spatial modes given by the Legendre polynomials $`P_n(x)`$. The Kohn sloshing mode at $`\omega `$ is $`n=1`$, and $`n=2`$ is the breathing mode at $`\sqrt{3}\omega `$. In this mode the width of the condensate oscillates, with a corresponding density oscillation.

Fig. 4 shows that as the pumping rate $`r`$ increases, the slow beam oscillations change from the frequency of the condensate breathing mode to that of the sloshing mode. This is because the condensate itself changes from breathing to sloshing at higher pumping rates. As Kneer et al.  noted, the sloshing mode is spatially asymmetric and hence must be excited by a spatially asymmetric perturbation, such as the output beam. This is consistent with the observed increase of the sloshing mode power with the density of the output beam, that is with pumping rate. On the other hand, mode damping due to three-body recombination increases rapidly with density. We have verified that in a simple model of the Kneer et al.  form, the breathing mode decays with pumping rate.

Although the frequencies in the beam are primarily determined by the condensate, there is a complex relationship between the power spectra of the condensate and of the beam. In particular, the beam spectrum reflects the spatial dependence of the condensate spectrum. Consider, for example, the condensate Kohn sloshing mode. The magnitude of the local density change varies with the spatial derivative of the density profile. It is maximum where the condensate has greatest slope, and zero at the center where the slope is zero. That the Kohn mode frequency is prominent in the power spectrum of the output beam indicates that the coupling between the condensate and the beam is not confined to the center of the condensate. Similar observations apply for the high frequency components of the power spectrum, the envelopes of which differ in detail between the condensate and the beam.

As the pumping rate increases, Fig. 4 shows that the frequency of the high frequency group increases and its power decreases rather dramatically. A decrease of spectral power with pumping rate, due to gain saturation, is characteristic of optical lasers. Since three-body recombination contributes to gain saturation in atom lasers, we have investigated its role in a simplified pumped atom laser model similar to that of Kneer et al. . That is, with a spatially uniform phenomenological loss rather than an explicitly modelled output beam. Without three-body recombination the condensate density power spectrum simply rolls off at high frequencies. Adding three-body recombination produces a high frequency spectral group, such as seen in Fig. 4. Hence it is an important factor in the dynamics of the system.

In conclusion, we have explored some of the rich and complex behavior of the atom laser with a model incorporating important experimentally relevant physics, such as three-body recombination. The latter was found to strongly influence the spectrum of the atom laser output beam. This work is a step towards the future goal of understanding and measuring the quantum noise properties of the atom laser.