text:
# Bose-Einstein condensation in trapped bosons: A Variational Monte Carlo analysis

## I Introduction

The spectacular demonstration of Bose-Einstein condensation (BEC) in gases of alkali atoms $`{}_{}{}^{87}Rb`$, $`{}_{}{}^{23}Na`$, $`{}_{}{}^{7}Li`$ confined in magnetic trapsanderson95 ; davis95 ; bradley95  has led to an explosion of interest in confined Bose systems. Of interest is the fraction of condensed atoms, the nature of the condensate, the excitations above the condensate, the atomic density in the trap as a function of Temperature and the critical temperature of BEC, $`T_c`$. The extensive progress made up to early 1999 is reviewed by Dalfovo et al.dalfovo99 .

A key feature of the trapped alkali and atomic hydrogen systems is that they are dilute. The characteristic dimensions of a typical trap for $`{}_{}{}^{87}Rb`$ is $`a_{h0}=\left(\mathrm{}/m\omega _{}\right)^{\frac{1}{2}}=12\times 10^4`$ Å (Ref. 1). The interaction between $`{}_{}{}^{87}Rb`$ atoms can be well represented by its s-wave scattering length, $`a_{Rb}`$. This scattering length lies in the range $`85<a_{Rb}<140a_0`$ where $`a_0=0.5292`$ Å is the Bohr radius gardner95 . The definite value $`a_{Rb}=100a_0`$ is usually selected and for calculations the definite ratio of atom size to trap size $`a_{Rb}/a_{h0}=4.33\times 10^3`$ is usually chosen dalfovo99 . A typical $`{}_{}{}^{87}Rb`$ atom density in the trap is $`n10^{12}10^{14}`$ atoms/cm3 giving an inter-atom spacing $`\mathrm{}10^4`$ Å. Thus the effective atom size is small compared to both the trap size and the inter-atom spacing, the condition for diluteness (i.e., $`na_{Rb}^310^6`$ where $`n=N/V`$ is the number density). In this limit, although the interaction is important, dilute gas approximations such as the Bogoliubov theorybogoliubov47 , valid for small $`na^3`$ and large condensate fraction $`n_0=N_0/N`$, describe the system well. Also, since most of the atoms are in the condensate (except near $`T_c`$), the Gross-Pitaevskii equationgross61 ; pitaevskii61  for the condensate describes the whole gas well. Effects of atoms excited above the condensate have been incorporated within the Popov approximationhutchinson97 . One of the chief purposes of this paper is to go beyond the dilute limit, to test the limits of these approximations and to explore the properties of the trapped Bose gas as $`na^3`$ increases between the dilute limit and the dense limit. We use Variational Monte Carlo (VMC) methods. We increase the density by increasing both $`N`$ and the s- wave scattering length up to the value $`na^30.21`$ which describes liquid 4He at SVP when the 4He atoms are represented by hard spheres of diameter $`a=2.203`$ Åkalos74 .

In addition to the mean-field theories noted above, the trapped Bose gas at finite temperatures has been investigated using Path Integral Monte Carlo (PIMC) methods. Krauthkrauth96  simulated 10,000 atoms in a spherical trap with ratio of scattering length $`a`$ to trap length given by $`a/a_{h0}=4.3\times 10^3`$ noted above. He showed that the critical temperature $`T_c`$ is lowered compared to the ideal Bose gas as a result of interaction. $`T_c`$ is lower because the repulsion between the atoms spreads the atoms in the trap and lowers the density compared to the non- interacting case. The same result has been obtained in mean-field approximations dalfovo99 . Krauth also showed that, in the dilute limit, the condensed atoms are highly concentrated at the center of the trap while the uncondensed or thermal atoms are spread out over a wide range, are dilute and well approximated by a classical, ideal gas. There is little interaction between the condensed and uncondensed components (see also dalfovo99 ).

Grüter et al.gruter97  evaluated $`T_c`$ for a uniform, bulk, hard sphere Bose gas over a wide density range, from dilute to liquid 4He densities. They find that $`T_c`$ is increased above the ideal Bose gas value by interaction in the dilute range. In the uniform gas case the density is not changed by interactions. At liquid 4He densities, $`T_c`$ is decreased by interactionpollock92 ; wilks67 .

Holzmann et al.holzmann99  made a direct comparison between Hartree-Fock (HF) and PIMC calculations of the number density $`N(r)`$ of atoms in a trap. The atoms were again represented by hard spheres with $`a/a_{h0}=0.0043`$. From $`N(r)`$, they find for temperatures near $`T_c`$ that the condensate fraction $`N_0`$ is larger in the exact PIMC evaluation than in the HF approximation. The energy beyond the Hartree-Fock approximation is often denoted the “correlation” energy. This correlation apparently allows $`N_0`$ to increase at a given $`T`$ and allows condensation to begin at a higher temperature. At lower temperature $`T\stackrel{<}{}\text{ }0.75T_c`$, there is excellent agreement between the PIMC and HF $`N(r)`$. The increase in $`N_0`$ with exact representation of the interaction effects is consistent with the corresponding increase in $`T_c`$ with interaction in the uniform Bose gas.

Giorgini et al.giorgini  have evaluated the ground state energy $`E/N`$ and the condensate fraction $`N_0/N`$ at $`T=0`$ K of a uniform Bose gas over a wide density range $`(10^6na^310^1)`$ using Green Function Monte Carlo (GFMC) methods. They find that the mean field results of $`E/N`$ and the Bogoliubov result for $`N_0/N`$ agree well with the GFMC values in the density range $`10^6na^310^3`$. However, there are clear differences at higher densities $`na^310^3`$ (helium density is $`na^30.21`$). The results are not sensitive to reasonable variation of the inter-boson potentialgiorgini .

In this context, we have evaluated the ground state properties of a trapped, hard sphere Bose gas over a wide range of densities using Variational Monte Carlo (VMC) methods. We begin in the dilute limit, small $`N`$ and $`a/a_{ho}=0.0043`$ corresponding to $`{}_{}{}^{87}Rb`$ in a trap, and increase both $`N`$ and $`a`$ separately to increase the density up to liquid 4He densities. At the lower densities we compare the energy $`E/N`$ and root mean square amplitudes $`<x^2+y^2>`$ and $`<z^2>`$ of atoms in an anisotropic trap with Gross-Pitaevskii (GP) resultsdalfovo96 . The two methods agree well at low densities but even at densities $`na^310^5`$ small differences in $`E/N`$ are readily apparent. Also, at higher densities we find the effects of interaction depend separately on $`N`$ and $`a/a_{ho}`$, not simply in the product $`Na/a_{ho}`$ as appears in the GP theory. As density is increased still further, we find that the condensate is no longer concentrated at the center of the trap. Rather increased interaction (increased $`a/a_{ho}`$) depletes the condensate at the center and the condensate appears at the edges of the trap – as found in liquid 4He droplets (Lewart et al.lewart88 ). Also, as density increases, correlations in the single particle density appear, reflecting the interaction in a confined space, in the same way that interaction at higher density introduces correlation in the pair correlation function in the uniform case. We evaluate the momentum distribution and condensate fraction over the dilute to dense range and compare with mean-field resultsjavanainen96 .

In section 2 we introduce the Hamiltonian, the wave function and the MC method and the definition of the natural orbitals. The results are presented in section 3 and discussed in section 4.

## II Bosons in a Harmonic Trap

### II.1 The System

We consider N bosons of mass $`m`$ confined in an external trapping potential, $`V_{ext}(𝐫)`$, and interacting via a two-body potential $`V_{int}(𝐫_1,𝐫_2)`$. The Hamiltonian for this system is:

$$H=\underset{i}{\overset{N}{}}(\frac{\mathrm{}^2}{2m}_i^2+V_{ext}(𝐫_i))+\underset{i<j}{\overset{N}{}}V_{int}(𝐫_i,𝐫_j).$$
(1)
We consider both a spherically symmetric (S) harmonic trap and an elliptical (E) harmonic trap,

$$V_{ext}(𝐫)=\{\begin{array}{cc}\frac{1}{2}m\omega _{ho}^2r^2\hfill & (S)\hfill \\ \text{}\frac{1}{2}m[\omega _{ho}^2(x^2+y^2)+\omega _z^2z^2]\hfill & (E)\hfill \end{array}$$
(2)
Here $`\omega _{ho}^2`$ defines the trap potential strength. In the case of the elliptical trap, $`V_{ext}(x,y,z)`$, $`\omega _{ho}=\omega _{}`$ is the trap frequency in the perpendicular or $`xy`$ plane and $`\omega _z`$ the frequency in the $`z`$ direction. The mean square vibrational amplitude of a single boson at $`T=0K`$ in the trap (2) is $`<x^2>=(\mathrm{}/2m\omega _{ho})`$ so that $`a_{ho}(\mathrm{}/m\omega _{ho})^{\frac{1}{2}}`$ defines the characteristic length of the trap. The ratio of the frequencies is denoted $`\lambda =\omega _z/\omega _{}`$ leading to a ratio of the trap lengths $`(a_{}/a_z)=(\omega _z/\omega _{})^{\frac{1}{2}}=\sqrt{\lambda }`$.

We represent the inter boson interaction by a pairwise, hard core potential

$$V_{int}(r)=\{\begin{array}{cc}\mathrm{}\hfill & ra\hfill \\ 0\hfill & r>a\hfill \end{array}$$
(3)
where $`a`$ is the hard core diameter of the bosons. Clearly, $`V_{int}(r)`$ is zero if the bosons are separated by a distance $`r`$ greater than $`a`$ but infinite if they attempt to come within a distance $`ra`$.

The weak interaction limit is $`a<<a_{ho}`$ and $`a<<n^{\frac{1}{3}}`$ (where $`n=N/V`$ is the local number density), a hard core diameter small compared to the dimensions of the trap and compared to the inter-particle spacing $`l=(V/N)^{\frac{1}{3}}`$. For trapped alkali atoms we have typically $`na^3\stackrel{<}{}\text{ }10^5`$. Introducing lengths in units of $`a_{ho}`$, $`rr/a_{ho}`$, and $`\mathrm{}\omega _{ho}`$ as units of energy as in dalfovo99 , the Hamiltonian is:

$$H=\underset{i}{\overset{N}{}}\frac{1}{2}(_i^2+x_i^2+y_i^2+\lambda ^2z_i^2)+\underset{i<j}{}V_{int}(|𝐫_i𝐫_j|).$$
(4)
Since there is Bose condensation, we have $`n\lambda _T^3\stackrel{>}{}\text{ }2.616`$ where $`\lambda _T`$ is the atomic thermal wavelength. Thus we are in the regime where the atomic wavelength is long compared to the hard core diameter, $`\lambda _T>>a`$ or $`ka<<1`$ where $`k2\pi /\lambda _T`$. The scattering of two particles interacting via a hard core potential in the limit $`ka<<1`$ is purely $`s`$ wave with scattering length $`a`$. If we approximate the full potential between the two particles by a contact potential,

$$v(r)=g\delta (𝐫)=\frac{4\pi \mathrm{}^2a}{m}\delta (𝐫)$$
(5)
the scattering length in this limit between the two is again purely $`s`$ wave with scattering length $`a`$. Thus we may compare directly results calculated using a hard core potential (3) and with a contact potential approximation (5) in the regime $`a<<\lambda _T`$. Specifically, we may compare the present MC results with results calculated using (5) and the mean field, Gross- Pitaevskii (GP) equation. This comparison is especially interesting in the dilute limit $`na^3<<1`$. At high densities $`(na^3\stackrel{>}{}\text{ }0.1)`$, we expect short range, pair correlations induced by the hard core to be important and short range correlations are not well described by a mean field theory.

### II.2 Wave Function

To describe the ground state of the $`N`$ bosons, we introduce a variational trial wave function which is a product of a single particle function $`g(𝐫)`$ and a pair Jastrow function jastrow  $`f(|𝐫_1𝐫_2|)`$,

$$\mathrm{\Psi }_\nu (𝐫__1..𝐫__N,\alpha ,\beta )=\underset{i}{}g(\alpha ,\beta ,𝐫_i)\underset{i<j}{}f(a,|𝐫_i𝐫_j|)$$
(6)
where $`\alpha `$ and $`\beta `$ are the variational parameters. We select a single particle function,

$$g(\alpha ,\beta ,𝐫_i)=\mathrm{exp}\left[\alpha (x_i^2+y_i^2+\beta z_i^2)\right]$$
(7)
which is a HO ground state function having two variational parameters, $`\alpha `$ and $`\beta `$. For spherical traps, $`\beta =1`$, and for noninteracting bosons $`(a=0)`$, $`\alpha =1/2a_{ho}^2`$. For the pair function we select the exact solution of the Schroedinger equation for two particles interacting via the hard core potential (3) in the limit $`k0`$, i.e.

$$f(a,r)=\{\begin{array}{cc}(1\frac{a}{r})\hfill & r>a\hfill \\ 0\hfill & ra\hfill \end{array}$$
(8)
(see for example Huang huang63 ). The $`\mathrm{\Psi }_\nu (𝐫__1..𝐫__N)`$ therefore has the correct form for small $`|𝐫_i𝐫_j|`$ and has two variational parameters $`\alpha `$ and $`\beta `$ that describe the spread of the bosons in the trap as the hard core diameter is increased. By constructing the wave function in this way, we limit the number of variational parameters while preserving the correct functional form in the $`a0`$ limit. However, the lack of any variational parameters in the Jastrow term is a potential source of inaccuracy.

We then minimize the expectation value of the Hamiltonian as obtained from

$$<H>=\frac{d𝐫_1..d𝐫_N\mathrm{\Psi }_\nu ^{}\mathrm{\Psi }_\nu (\frac{H\mathrm{\Psi }_\nu }{\mathrm{\Psi }_\nu })}{d𝐫_1..d𝐫_N\mathrm{\Psi }_\nu ^{}\mathrm{\Psi }_\nu }$$

with respect to $`\alpha `$ and $`\beta `$ using the Metropolis Monte Carlo method of integration. This is accomplished by using a Metropolis random walk metropolis  to generate a set of $`M`$ particle configurations, $`\mathrm{\Omega }_1..\mathrm{\Omega }_M`$ which conform to the probability distribution $`|\mathrm{\Psi }_\nu |^2`$. We then approximate $`<H>`$ by summing over the ‘local energy’ as follows

$$<H>\stackrel{>}{}\text{ }\frac{1}{M}\underset{n=1}{\overset{M}{}}\left(\frac{H(\mathrm{\Omega }_n)\mathrm{\Psi }_\nu (\mathrm{\Omega }_n)}{\mathrm{\Psi }_\nu (\mathrm{\Omega }_n)}\right)$$

For a review of the Variational Monte Carlo method, see vmc\_Review .

### II.3 Condensate and Natural Orbitals

A goal here is to calculate the condensate fraction and condensate density in the ground state. To do this we require a definition of the condensate single particle state. Following Penrose and Onsager, Löwdin and others onsager56 ; lowdin55 , we take the one-body density matrix (OBDM) as the fundamental quantity for an interacting system and define the natural single particle orbitals (NO) in terms of the OBDM.

The OBDM is baym76

$$\rho (𝐫^{},𝐫)=<\widehat{\mathrm{\Psi }}^{}(𝐫^{}),\widehat{\mathrm{\Psi }}(𝐫)>,$$
(9)
where $`\widehat{\mathrm{\Psi }}(𝐫)`$ is the field operator that annihilates a single particle at the point $`𝐫`$ in the system. At $`T=0K`$, the expectation value is evaluated using the wave function $`\mathrm{\Psi }_\nu (𝐫)`$ in (6). To define the NO, we introduce a set of single particle states having wave functions $`\varphi _i(𝐫)`$ and expand $`\widehat{\mathrm{\Psi }}(𝐫)`$ in terms of the operators $`\widehat{a}_i`$ which annihilate a particle from state $`|i>`$,

$$\widehat{\mathrm{\Psi }}_0=\underset{i}{}\varphi _i(𝐫)\widehat{a}_i.$$
(10)
Requiring that the $`\widehat{a}_i`$ satisfy the usual commutation $`([\widehat{a}_i^{},\widehat{a}_j]=\delta _{ij})`$ and number relations $`(<\widehat{a}_i^{}\widehat{a}_j>=N_i\delta _{ij})`$, we have

$$\begin{array}{c}\hfill \rho (𝐫,𝐫^{})=_{ij}\varphi _j^{}(𝐫^{})\varphi _i(𝐫)N_i\delta _{ij}\\ \hfill =_{ij}\varphi _j^{}(𝐫)\varphi _i(𝐫^{})N_i\delta _{ij}\end{array}$$
(11)
This may be taken as the defining relation of the NO, $`\varphi _i(𝐫)`$. Specifically, we have from (11),

$$d𝐫..d𝐫^{}\varphi _i^{}(𝐫)\rho (𝐫,𝐫^{})\varphi _j(𝐫^{})=N_i\delta _{ij},$$
(12)
so that the NO may be obtained by diagonalizing the OBDM. The eigenvectors are the NO and the eigenvalues are the occupation, $`N_i`$, of the orbitals. The condensate is the orbital having the highest occupation, denoted $`\varphi _0(r)`$ and the condensate fraction is $`n_0=N_0/N`$.

The relations (11) and (12) involve the vector $`𝐫`$ and $`𝐫^{}`$ and cannot be solved directly as matrix equations. To obtain matrix equations, we restrict ourselves to spherical traps and seek equations for the radial component of the NO. Assuming the potential seen by a single particle, including inter-particle interaction, is spherically symmetric, the NO will have the form

$$\varphi _i(𝐫)=R_{nl}(r)Y_{lm}(\theta ,\varphi )$$
(13)
Where $`Y_{lm}(\mathrm{\Omega })`$ are the spherical harmonics and $`R_{nl}(r)`$ is the radial wave function and $`i=n,l,m`$ are the state indices. We expand the OBDM in its angular momentum $`(l)`$ components $`\rho _l(r,r^{})`$ as

$$\rho (𝐫,𝐫^{})=\underset{l}{}\frac{(2l+1)}{4\pi }P_l(\widehat{𝐫}\widehat{𝐫}^{})\rho _l(r,r^{})$$
(14)
where $`P_l(cos\theta )`$ are Legendre Polynomials in the angle between $`\widehat{𝐫}`$ and $`\widehat{𝐫}^{}`$ and

$$\begin{array}{c}\rho __l(r__1,r__1^{})=\hfill \\ \\ d\mathrm{\Omega }_1d𝐫__2..d𝐫__N\mathrm{\Psi }^{}(𝐫__1..𝐫__N)P_l(\widehat{r}__1\widehat{r}__1^{})\mathrm{\Psi }(𝐫__1^{}..𝐫__N).\hfill \end{array}$$
(15)
Substituting (13) and (14) into (11) and using the properties of the spherical harmonics, $`Y_{lm}(\mathrm{\Omega })`$, we obtain a relation equivalent to (11) for each $`l`$ component $`\varphi _{nl}(r)`$ of the NO,

$$\rho __l(r,r^{})=\underset{n}{}\varphi _{nl}(r)\varphi _{nl}(r^{})N_{nl}$$
(16)
where $`i=nl`$. To solve this equation readily as a matrix equation, we introduce the radial function

$$u_{nl}(r)=r\varphi _{nl}(r).$$
(17)
The $`u_{nl}(r)`$ is well behaved at $`r0`$ and has dimensions $`L^{\frac{1}{2}}`$ like a one-dimensional wave function. In terms of the $`u_{nl}(r)`$, the defining relations are:

$$[r\rho __l(r__1,r__1^{})r^{}]=\underset{n}{}u_{nl}(r)u_{nl}^{}(r^{})N_{nl}$$
(18)
and

$$𝑑r[r\rho __l(r__1,r__1^{})r^{}]u_{nl}(r^{})=u_{nl}(r)N_{nl}$$
(19)
The $`[r\rho __l(r__1,r__1^{})r^{}]`$ serves as a one-dimensional OBDM along the radial coordinate. This $`1D`$ matrix relation may be solved numerically on a grid in $`r`$ to obtain the $`u_{nl}(r)`$ as eigenvectors and the $`N_{nl}`$ as eigenvalues. The condensate orbital is

$$\varphi _0(r)=\frac{1}{\sqrt{4\pi }}u_{00}(r)/r.$$
(20)
The momentum distribution may be obtained from the OBDM as well given that the orbitals in momentum space are

$$\stackrel{~}{\varphi }_i(𝐤)=\frac{1}{(2\pi )^{\frac{3}{2}}}𝑑𝐫e^{i𝐤𝐫}\varphi _i(𝐫)$$
(21)
and the momentum distribution is

$$\stackrel{~}{\rho }(𝐤)=\underset{i}{}n_i|\stackrel{~}{\varphi }_i(𝐤)|^2.$$
(22)
Substituting (21) into (22) obtains an expression for $`\stackrel{~}{\rho }(𝐤)`$ in terms of the OBDM:

$$\stackrel{~}{\rho }(𝐤)=\frac{1}{(2\pi )^3}𝑑𝐫𝑑𝐫^{}\rho (𝐫,𝐫^{})e^{i𝐤(𝐫𝐫^{})}.$$
(23)
The momentum distribution for a spherically symmetric system is therefore

$$\begin{array}{c}\stackrel{~}{\rho }(𝐤)=\stackrel{~}{\rho }(|𝐤|)=\hfill \\ \text{}\text{}\frac{1}{2\pi ^2}d𝐫_\mathrm{𝟏}..d𝐫_𝐍|\mathrm{\Psi }_\nu (𝐫_1..𝐫_N)|^2\hfill \\ \times _0^{\mathrm{}}drr^2\frac{sin(kr)}{kr}\frac{\mathrm{\Psi }_\nu (𝐫_1+𝐫,𝐫_2..𝐫_N)}{\mathrm{\Psi }_\nu (𝐫_1..𝐫_N)}\hfill \end{array}$$
(24)
where the orientation of $`𝐫`$ may be chosen arbitrarily.

## III Results

In this section we present the results of the present Variational Monte Carlo (VMC) calculation of properties of bosons at $`T=0`$ K confined in a harmonic trap. To begin, Fig. 1 shows the condensate orbital $`\varphi _0(r)`$ for independent, non-interacting bosons in a spherical harmonic trap compared to the harmonic oscillator (HO) ground state function for the same trap. The $`\varphi _0(r)`$ is calculated by evaluating the one body density matrix (OBDM) and diagonalizing it numerically to obtain the single particle orbitals and their occupation as discussed in section 2. For no interaction, only $`\varphi _0(r)`$ is occupied $`(n_0=1)`$. The excellent agreement between $`\varphi _0(r)`$ and the HO ground state function for all $`r`$ provides a good check of the method for non-interacting bosons. The statistical sample is proportional to $`r^2`$ and the statistics become poor at $`r0`$.

Fig. 2 shows the energy per particle, $`E/N`$, of a gas of $`N`$ weakly interacting hard sphere bosons in an ellipsoidal trap. The ratio of the characteristic length of the short axis ($`z`$ axis) to the longer perpendicular ($`x`$ and $`y`$) axis of the trap is $`(a_z/z_{})=1/\sqrt{\lambda }`$;, $`\lambda =\sqrt{8}`$. The hard sphere diameter, $`a`$, (scattering length) of the bosons corresponds to $`{}_{}{}^{87}Rb`$ atoms in an ellipsoidal trap with $`a/a_{ho}=0.00433`$ with $`a_{}=a_{ho}`$. In Fig. 2 the dots are the present VMC $`E/N`$ of the whole Bose gas while the dashed line is the $`E/N`$ of the condensate calculated by Dalfovo and Stringaridalfovo96  using the Gross-Pitaevskii (GP) equation. Our purpose is to make comparisons between the present VMC calculation and GP equation results across the dilute regime for which GP is expected to be valid. The region $`100N20000`$ corresponds to atom densities at the center of the trap of $`2\times 10^6\stackrel{<}{}\text{ }na^3\stackrel{<}{}\text{ }2\times 10^5`$. In the earliest experimentsanderson95  with $`{}_{}{}^{87}Rb`$ there were typically $`N=10,000`$ atoms in the trap. In more recent experiments $`N`$ is larger, $`N10^510^6`$.

At small $`N`$ values in the dilute limit, there is excellent agreement between the GP and VMC energies. In this regime, the mean-field, GP equation is expected to be accurate. As N increases a clear difference between the VMC and GP $`E/N`$ values emerges. We find that for a fixed scattering length, the difference $`\delta (E/N)=(E_{MC}E_{GP})/N`$ is proportional to $`N^{\frac{3}{5}}`$ and can be well represented as

$$\delta (E/N)=[(3.0\times 10^5\pm 10^6)\times (N)^{\frac{3}{5}}0.052]\mathrm{}\omega _{}.$$
(25)
The difference is $`1.8\%`$ at $`N=10,000`$ and $`2.5\%`$ at $`N=20,000`$. The difference, we believe, arises largely because there is excitation of bosons above the condensate. $`E_{MC}/N`$ includes the excited atoms while $`E_{GP}/N`$ is the energy of the condensate alone. We return to this point in the Discussion.

Fig. 3 compares the root mean square displacement of hard sphere bosons from the center of the same anisotropic trap discussed in Fig. 2 calculated using VMC and the GP equation. The upper (lower) line is the radial displacement along the $`(z)`$ direction. The agreement between the GP and the VMC displacements is excellent, right up to $`N=20,000`$ for $`a=a_{Rb}`$.

It is interesting that the VMC and GP displacements agree well while the energies in Fig. 2 differ for $`N\stackrel{>}{}\text{ }10^4`$. Essentially, the $`E/N`$ is very sensitive to the few high energy bosons at the edges of the trap. That is, the small number of atoms having large displacements increase the energy significantly but change $`<r^2>`$ little. An example of this effect is found in the Thomas-Fermi approximation where the density is cut off to zero at a specific radius dalfovo99 . In the TF model there is no tail in the density reaching up to large $`r`$ values. The cut off changes $`<r^2>`$ little but $`E/N`$ is significantly affected dalfovo99 .

Fig. 4 again shows $`E/N`$ for hard sphere bosons in an anisotropic trap as a function of $`Na/a_{ho}`$. However, in this case the product $`Na/a_{ho}`$ is adjusted by varying both $`N`$ and $`a/a_{ho}`$. The star symbols show $`E/N`$ for $`a=a_{Rb}`$ and $`2000N20,000`$, the crosses for $`a=10a_{Rb}`$ and $`200N2,000`$ and the square $`a=20a_{Rb}`$ and $`100N1,000`$. If the impact of interaction and $`E/N`$ depended solely on the product $`Na/a_{ho}`$, as is the case in the mean-field, GP equation, all three lines in Fig. 4 would coincide. $`E/N`$ clearly depends separately on $`N`$ and $`a/a_{ho}`$, even in the region of $`N=10,00020,000`$. The separate dependence is not large at $`Na/a_{ho}20`$ but becomes increasingly large as $`Na/a_{ho}`$ increases. Also, at these and larger densities, the parameter which determines the magnitude of corrections arising from interactions is $`N^{\frac{3}{5}}(a/a_{ho})^{\frac{8}{5}}`$. Apparently the interaction effects depending on $`Na/a_{ho}`$ are valid only in the limit of small $`a`$ ($`a/a_{ho}<<1`$). We examine the functional form of the separate dependence of $`E/N`$ on $`N`$ and $`a/a_{ho}`$ in the discussion section.

Having investigated lower densities and made comparisons with results obtained using the GP equation, we now turn to higher densities and bosons represented by hard spheres having larger hard core diameters. We evaluate the OBDM and the density for these cases going up to densities comparable to liquid 4He droplets. Fig. 5 shows the condensate orbital (wave function) for 128 bosons in a spherical harmonic trap as the hard core radius, $`a`$, is increased from zero up to $`a=64a_{Rb}(a_{Rb}/a_{ho}=0.00433)`$. The case $`a/a_{Rb}=1`$ corresponds to $`N=128`$ Rb atoms in a spherical trap. Clearly the condensate orbital spreads out in the trap as $`a`$ increases. At $`a/a_{Rb}=64`$, the condensate density is effectively constant in the trap out to nearly three times the trap length parameter $`a_{ho}`$. For these larger core radii, the appropriate measure of the interaction is $`na^3`$ where $`n=N/V`$. For $`a=64a_{Rb}`$, $`na^32\times 10^3`$ at the center of the trap.

In Fig. 6 we show the total density $`\rho (r)`$ and the density of atoms in the condensate orbital as $`N|\varphi _0(r)|^2`$ for 64 bosons in a spherical trap with $`a`$ increased to $`128a_{Rb}`$, $`256a_{Rb}`$ and $`512a_{Rb}`$. Since we plot $`N|\varphi _0(r)|^2`$ rather than $`N_0|\varphi _0(r)|^2`$, “condensate density” can exceed the total density. At $`a=256a_{Rb}`$ the condensate is $`n_0=N_0/N=20\%`$. We note that as $`a`$ increases, the condensate moves away from the center of the trap. At $`a=512a_{Rb}`$, the condensate is at the edges of the trap. When $`na^3`$ is large at the center of the trap, the maximum condensate density is in the region of lower particle density found at the edge of the trap, as calculated for liquid 4He droplets lewart88 . Thus the location of the condensate is entirely different at small and large scattering length. In a slave boson approach, depletion of the condensate at the center of the trap has also been demonstrated ziegler97 .

Also, at large $`a/a_{ho}`$ the total density develops correlations. These correlations reflect the inter-boson correlations induced by the hard core interaction. For $`a=512a_{Rb}`$ we have $`a/a_{ho}2.2`$. With a peak in the density at $`r=0`$, we expect the first minimum in the density at $`ra/a_{ho}2.2`$ as in the upper frame of Fig. 6.

Fig. 7 gives the fraction of bosons in the condensate orbital, $`n_0`$, calculated by VMC and diagonalization of the OBDM corresponding to the condensate orbitals shown in Fig. 5. The $`n_0`$ values, as in Fig. 5, are for 128 hard sphere bosons in a spherical trap with hard sphere radius $`a/a_{Rb}=1,2,4,16,32`$ and 64 e.g. $`a=64\times a_{Rb}=64\times 4.33\times 10^3a_{ho}=0.277a_{ho}`$. The Bogoliubov bogoliubov47  result for $`n_0`$ adapted to a spherical trapjavanainen96  is shown as a dashed line. Visible departures of the VMC $`n_0`$ from the Bogoliubov result begin at $`n_00.96`$ (a depletion of 4%) corresponding to a density at the center of the trap $`na^31\times 10^3`$. In the uniform gas case, the Bogoliubov result remains accurate up to a condensate fraction $`n_00.89`$ (11% depletion) which occurs at a density $`na^35\times 10^3`$ giorgini . The Bogoliubov approximation has a more limited range of application for bosons in a trap because the interaction changes the density profile (and the shape of the condensate wave function) as well as simply depleting the condensate and $`n_0`$ depends on the density and shape of the condensate wave function. In the uniform case, the density cannot change.

In Fig 8., we again show the condensate fraction, $`n_0`$, as a function of the ratio of the scattering length, $`a`$, to the characteristic trap length, $`a_{ho}`$. Here, $`n_0`$ is given for three different numbers of particles: $`N=64,128`$ and $`256`$. The corresponding value of the particle density, $`na^3`$, for the $`64`$ particle case is shown on the top axis for reference. At liquid helium densities, the condensate fraction is roughly twice that of bulk liquid $`{}_{}{}^{4}He`$. This difference can be understood by noting the shape of the radial condensate density shown for the $`a=256a_{Rb}`$ case in Fig. 6. Here, we see that while the maximum particle density occurs at the center the trap, the condensate density is peaked in the low density region at the edge of the cloud. This dilute region allows for a larger fraction of particles to occupy the condensate orbital than in an uniform system at $`{}_{}{}^{4}He`$ densities.

Fig. 9 shows the effect of increased scattering length, $`a`$, on the momentum distribution of particles in a isotropic harmonic trap at $`T=0`$. The values for the scattering length and the trap configuration under consideration correspond to those shown for the spatial distribution in Fig. 5.

## IV Discussion

In this section we compare the present MC values for the density, condensate fraction and energy of bosons in a trap with mean field (MF) and Thomas-Fermi (TF) approximation results. The aim is to assess the limits of applicability of the MF and TF expressions and to investigate the origin of any differences between MF and MC values. We start with the density at the center of the trap, $`n(0)`$, since this density is needed in MF expressions for the depletion of the condensate and for the energy. Also, comparisons  of the density calculated in the TF approximation and in mean field, Gross-Pitaevskii (GP), approximations show that $`n_{TF}(0)`$ is accurate for large $`Na/a_{ho}`$ ($`Na/a_{ho}100`$ ).

The density of N independent bosons in an asymmetric trap is $`n_{ho}(𝐫)=\varphi _{ho}^2(𝐫)=N/(\pi ^{\frac{3}{2}}a_xa_ya_z)\mathrm{exp}[(x^2/a_x^2+y^2/a_y^2+z^2/a_z^2)]`$. For the elliptical trap considered above, $`a_x=a_y=a_{ho}`$ and $`a_z=a_{ho}/\sqrt{\lambda }`$, the density at the center of the trap is dalfovo99

$$n_{ho}(0)=Na_{ho}^2a_z/\pi ^{\frac{3}{2}}=\sqrt{\lambda }N/\pi ^{\frac{3}{2}}a_{ho}^3$$
(26)
The $`\varphi _0(r)`$ for a spherically symmetric trap $`(\lambda =1)`$ is shown in Fig. 1. As interaction is increased (e.g. $`a/a_{ho}`$ is increased) the $`\varphi _0(r)`$ spreads out and the density $`n(0)`$ at the center of the trap decreases, as depicted in Fig. 5. For $`\lambda =1`$, the ratio of the density at the center with interaction to the $`n_{ho}(0)`$ for no interaction in the TF approximation is

$$\frac{n_{TF}(0)}{n_{ho}(0)}=(\frac{15^{\frac{2}{5}}\sqrt{\pi }}{8})(Na/a_{ho})^{\frac{3}{5}}.$$
(27)
This ratio calculated using MC, $`n_{MC}(0)/n_{ho}(0)`$, for increasing scattering length, $`a/a_{ho}=Fa_{Rb}=F\times 4.33\times 10^3`$ with $`F=1,4,8,16,32`$ and $`64`$ can be obtained from Fig. 5 and is listed in Table 1. The TF ratio agrees well with the MC ratio in the range $`4<F<16`$ or $`2<Na/a_{ho}<10`$. We expect the TF result to be inaccurate at small $`Na/a_{ho}`$ because the TF limit is a large $`Na/a_{ho}`$ approximation.

In addition, at large $`Na/a_{ho}`$, the density exceeds the dilute limit so that all mean field theories become inaccurate. The dilute limit is exceeded at $`n(0)a^3\stackrel{>}{}\text{ }10^3`$ which corresponds to $`F16`$ for the gas considered in Table 1. Between these limits however, for $`Na/a_{ho}\stackrel{>}{}\text{ }5`$ and $`n(0)a^3\stackrel{<}{}\text{ }10^3`$, the $`n_{TF}(0)`$ gives a good estimate of the boson number density at the center of the trap.

For a uniform Bose gas, the MF, Bogoliubov theory predicts a condensate fraction $`n_0=N_0/N=1(8/3)(na^3/\pi )^{\frac{1}{2}}`$ where $`n=N/V`$ is the uniform density bogoliubov47 . That is, the fraction of bosons excited out of the condensate is

$$\frac{\delta N}{N}=\frac{NN_0}{N}=\frac{8}{3}(\frac{na^3}{\pi })^{\frac{1}{2}}.$$
(28)
The corresponding result for bosons in a spherical harmonic trap is dalfovo99

$$\frac{\delta N}{N}=\frac{5\pi }{8}(\frac{n_{TF}(0)a^3}{\pi })^{\frac{1}{2}}.$$
(29)
Fig. 7 shows that the mean field expressions (28) and (29) agree well with the present MC values of $`n_0`$ for $`N=128`$ bosons up to $`a/a_{ho}16a_{Rb}`$ or up to $`n(0)a^310^3`$. Specifically, for $`F=8`$ at $`n_{MC}(0)a^3=2.4\times 10^4`$, the MC value for the depletion in Fig. 7 is $`\delta N/N=1.9\%`$ while the Bogoliubov result from (29) is $`\delta N/N=1.7\%`$. For larger values of $`n(0)a^3`$ the depletion is significantly underestimated by the Bogoliubov expressions. This agrees with the findings of MC determinations for $`n_0`$ in a uniform Bose gasgiorgini . In the uniform case, (28) agrees with MC values up to $`na^35\times 10^3`$ and underestimates depletion at larger $`na^3`$. The MF results are accurate up to somewhat higher densities in the bulk probably because the density itself does not depend on the interaction in a uniform Bose gas. Recent experiments cornish00  have obtained stable condensates at densities corresponding to $`na^310^2`$. At this density, the condensate fraction is $`.85`$ and effects resulting from depletion are expected to be significant.

For a uniform Bose gas, the MF, Bogoliubov expression for the energy including the leading correction arising from depletion is bogoliubov47

$$\frac{E}{N}=4\pi na^3[1+\frac{128}{15}(\frac{na^3}{\pi })^{\frac{1}{2}}].$$
(30)
The corresponding TF expression for bosons in a spherical trap ($`\lambda `$ = 1) is dalfovo99

$$\frac{E}{N}=\frac{5}{7}\mu _{TF}[1+\frac{7\pi }{8}(\frac{n_{TF}(0)a^3}{\pi })^{\frac{1}{2}}]$$
(31)
where $`\mu _{TF}=\frac{1}{2}\mathrm{}\omega _{}(15Na/a_{ho})^{\frac{2}{5}}`$ and from (26) and (27) $`n_{TF}(0)a^3=(15^{\frac{2}{5}}/8\pi )N^{\frac{2}{5}}(a/a_{ho})^{\frac{12}{5}}`$

The above expressions may be used to estimate the energy of bosons in an anisotropic trap by replacing $`a_{ho}`$ with the geometric mean $`a_g=(a_xa_ya_z)^{\frac{1}{3}}`$ which is $`a_g=a_{ho}\lambda ^{\frac{1}{6}}`$ for the trap discussed in Fig. 2. The energy of independent, non-interacting bosons in this anisotropic trap is

$$(E/N)E_{ho}=\mathrm{}\omega _{ho}(1+\lambda /2)=\mathrm{}\omega _{ho}(2.414)$$
(32)
$`(\omega _{ho}=\omega _{})`$.

In Fig. 2, ($`E/N`$) has the non-interacting boson limit at $`N0`$. As N increases, ($`E/N`$) increases. As seen from (31), the mean field energy in the TF approximation $`(E/N)_{TF}=5\mu _{TF}/7`$ increases with N as $`(E/N)_{TF}(Na/a_{ho})^{\frac{2}{5}}`$. The $`(E/N)_{MC}`$ in Fig. 2 follows the dependence approximately. However, direct evaluation of (31) shows that $`(E/N)_{TF}`$ underestimates the energy significantly, by 25% at N=20,000 (N$`a/a_{ho}`$ = 86.6). Thus while the Thomas Fermi density $`n_{TF}(0)`$ at the center of the trap is accurate, the TF energy is a poor approximation in this density regime. This is because the TF approximation underestimates the density at large $`r`$ and the high energy particles at large $`r`$ contribute significantly to the energy.

In Fig. 2 we see that the present Monte Carlo $`E/N`$ lies above the Gross-Pitaevskii $`E/N`$, by $`2.5\%`$ at $`N=20,000`$. A difference could arise for two reasons.

Firstly, the MC energy is the E/N of the whole Bose gas while the GP equation calculates the $`E/N`$ of the condensate only. To investigate the effects of depletion on $`E/N`$ we note, comparing (28) and (30), that the fractional change in energy arising from depletion, $`\delta E/E`$, can be related to the fraction of atoms excited out of the condensate as,

$$\delta E/E=(128/15)(na^3/\pi )^{\frac{1}{2}}=(16/5)\delta N/N$$
(33)
in the bulk. The corresponding equation for the trap is, from (29) and (31),

$$\delta E/E=(7/5)\delta N/N.$$
(34)
These are lowest order expressions. In Table 2 we list the depletion of the condensate $`\delta N/N`$ for the bosons in the anisotropic trap considered in Fig. 2 for $`N=5000,10,000`$ and $`20,000`$ predicted by (29). We expect these predictions to be accurate since $`n_{TF}(0)`$ is reliable and $`n_{TF}(0)a^3`$ is small. The interaction depletes the condensate $`0.5\%`$ at $`N=20,000`$. From the above connection between $`\delta E`$ and $`\delta N`$, the energy $`(E/N)`$ including depletion is expected to lie 0.7-1.6% above the energy of the condensate. On this basis, VMC and GP energies are consistent.

To explain the connection between the difference in VMC and GP energies and depletion of the condensate more fully, we have plotted $`\delta (E/N)=(E_{MC}E_{GP})/N`$ vs $`Na/a_{ho}`$ in Fig. 10. From (31), the difference from the mean field energy arising from depletion is

$$\begin{array}{cc}\hfill \delta (E/N)_{TF}=& (5\pi \mu _{TF}/8)(n_{TF}(0)a^3/\pi )^{\frac{1}{2}}\hfill \\ \hfill & N^{\frac{3}{5}}(a/a_{ho})^{\frac{8}{5}}\hfill \end{array}$$
(35)
in the TF limit. In Fig. 10, we show the difference in the VMC and GP energies, $`\delta (E/N)`$, verses the product $`Na/a_{ho}`$. This difference is obtained from Fig. 4 for four values of $`a/a_{ho}=1,5,10,20`$. The lines in Fig. 10 are fits of $`\delta (E/N)=mN^{\frac{3}{5}}(a/a_{ho})^{\frac{8}{5}}+b`$ to the data where the fitting parameters $`m`$ and $`b`$ where allowed to change for each $`a/a_{ho}`$ value. The good fit of these lines shows that $`\delta (E/N)`$ reflects the dependence on N expected for a difference in energy arising from exciting bosons out of the condensate. Thus the difference in $`(E/N)_{MC}`$ from $`(E/N)_{GP}`$ is consistent in magnitude and dependence on N and $`a/a_{ho}`$ with that expected for an $`(E/N)_{MC}`$ that includes bosons both in and above the condensate while $`(E/N)_{GP}`$ is the energy of bosons in the condensate only. Thus we believe the MC and GP energies differ chiefly because the MC includes “excited” particles while the GP energy does not.

Secondly, the present VMC $`E/N`$ is a genuine upper bound for the whole gas and could lie above the whole gas energy. We have not tested the sensitivity of $`E/N`$ to different choices of the trial wave function. In addition, while the present pair Jastrow function (8) is exact in the dilute limit it does not contain any variational parameters and is therefor not optimized for trapped hard spheres at higher densities. As a result, at least some of the difference in $`(E/N)_{MC}`$ from $`(E/N)_{GP}`$ could arise from the present choice of trial wave function. For example, Fabrocini and Polls (FP) have evaluated the whole $`E/N`$ for bosons in a spherical trap $`a_x=a_y=a_z=4.33\times 10^3a_{ho}`$ using correlated basis function and hypernetted-chain (HNC) methodsfabrocini99 . Both methods provide estimates of $`E/N`$ which lie below the present VMC $`E/N`$. At $`N=10^5`$ the density ($`n(0)a^3=2.5\times 10^5`$) in this trap is similar to that for $`N=2\times 10^4`$ in the present elliptical trap. At this density the HNC $`E/N`$ lies 0.8% above the GP energy of the condensate while the present VMC $`E/N`$ is 2.5% above GP. FP also consider a mean field model incorporating quantum depletion which predicts an increase in $`E/N`$ of 1.2% above the GP result at this density. Generally, the HNC energy in the trap lies below the energy expected for depletion and may be too low. For example, the HNC energy for a uniform gas of bosons lies above the energy expected for depletion. Definite resolution of these differences awaits a model independent evaluation of $`E/N`$ by Diffusion Monte Carlo methods.

Finally, an important result of the present MC evaluation is that as the boson density $`na^3`$ increases, the condensate gradually moves from the center of the trap to the edges of the trap as shown in Fig. 6. At large $`na^3`$, the condensate is at the edges of the trap. In this limit, the depletion of the condensate is large and the condensate seeks the regions of lowest total density which are at the surface of the trap. Both the condensed and uncondensed atoms must be included in the calculation to obtain this effect. This result is consistent with the calculations in liquid 4He droplets lewart88  which find the condensate concentrated at the surface of the droplet.

###### Acknowledgements.

Stimulating discussions with Charles W. Clark and Timothy Ziman are gratefully acknowledged. J.L.D. acknowledges support from DE-SGGFP NASA Grant No. NGT5-40024.


format: text

text:
# RXTE observations of XTE J2012+381 during its 1998 outburst.

## 1 Introduction

The XTE J2012+381 was discovered as a new transient source on May 27, 1998 by All Sky Monitor (ASM) aboard the Rossi X-ray Timing Explorer, RXTE (Bradt, Swank $`\&`$ Rothschild 1993) observatory (Remillard et al. 1998). The X-ray source was localized with an accuracy of 1 arcmin in multiple scans of the region by PCA/RXTE experiment (Marshall $`\&`$ Strohmayer 1998): $`R.A.=20^h12^m43^s,Dec=38^{}11.0\mathrm{}`$ (equinox 2000.0). Observations in the optical band by 1.0-m Jacobus Kapteyn Telescope showed the presence of possible optical counterpart USNO 1275.13846761 : $`R.A.=20^h12^m37.8^s,Dec=+38^{}11\mathrm{}00.6\mathrm{}`$ (equinox 2000.0)(Hynes & Roche 1998). Observations in the radio band by VLA revealed the presence of the likely counterpart of XTE J2012+381: $`R.A.=20^h12^m37.67,Dec=+38^{}11\mathrm{}01.2\mathrm{}`$ (equinox 2000.0, uncertainty $`0.5\mathrm{}`$) as well. Significant variability of the new radio source supports its association with X-ray transient (Hjellming et al. 1998).

Pointed instruments of RXTE – PCA and HEXTE observed the source quasi evenly from May to July 1998, providing a good coverage of the whole outburst. The flux from the source was above ASM/RXTE detection limit until December 2, 1998. In this Letter we present the results of a spectral and timing analysis of RXTE pointed observations.

## 2 Observations and data analysis

In our analysis we used all the publicly available data obtained from RXTE archive including 24 pointed observations. The 24 observations quasi evenly cover the 1998 outburst of the source with a total exposure of $`46`$ ks. The list of observations is presented in Table 1.

For data reduction we used the standard FTOOLS package version 5.0 For the spectral analysis we used PCA data collected in the 3-20 keV energy range. The response matrix was constructed for every observation using pcarmf v3.5. For the PCA background subtraction we applied a Very Large Events (VLE)-based model. The standard dead time correction procedure was applied to the PCA data. In order to account for the uncertainties of the response matrix, a 1$`\%`$ systematic error was added to the statistical error for each PCA energy channel.

HEXTE data reduction was done with FTOOLS 5.0 standard tasks and according to RXTE GOF recommendations.

We generated the energy spectra of XTE J2012+381 averaging the data over the whole single observations. The spectral data were approximated with the simplest two-component model: the sum of the XSPEC “multicolor disk black body” (Mitsuda et al. 1984) and a simple power law model with inclusion of a broad gaussian emission line at 6.4 keV. The results of the spectral approximation of the PCA data with analytical models described above are presented in Table 2.

For the timing analysis of XTE J2012+381 the PCA Generic Event timing mode data were used. We generated power density spectra (PDS) in the $`0.0014096`$ Hz frequency range using combined data of several consequitive observations in order to improve statistical significance of the results. The resulting spectra were logarithmically rebinned to reduce scatter at high frequencies and normalized to square root of fractional variability rms. The white noise due to the Poissonian statistics corrected for the dead-time effects, was subtracted (Vikhlinin, Churazov & Gilfanov (1994); Zhang et al. (1995); Revnivtsev, Gifanov & Churazov (2000)).

## 3 Results

The X-ray flux history of the XTE J2012+381 outburst based on RXTE/PCA and ASM data in 1.3-12 keV energy range is shown in Fig.1. The evolution of the source flux in the soft X-ray band ($`1.312`$ keV as well as in $`320`$ keV energy band) was characterized by the fast initial rise to a level of $`220`$ mCrab 111Assuming a distance of 10 kpc the luminosity in the peak of the outburst is $`7\times 10^{37}`$ erg s-1 in the $`320`$ keV energy band. on a time scale of $``$ a week followed by a $``$3 day long maximum and relatively slow decay, interrupted by secondary maximum $``$30 days after the beginning of the outburst. The subsequent outburst evolution was unusual because of long and powerful tertiary peak (close to 150 days from the beginning of the outburst).

The energy spectra of the X-ray transient were quite typical – a dominant soft thermal component with a rather weak power law tail. Unfortunately, the transient was not bright enough for high energy instrument aboard RXTE (HEXTE), so we could not trace the energy spectrum of XTE J2012+381 to energies higher than $``$80 keV.

XTE J2012+381 one more time demonstrated that the soft spectral component in the spectra of X-ray Novae rise not instantly, but slightly after the power law component. The similar behavior was observed in Nova Muscae 1991 (Ebisawa et al. (1994)), KS 1730-312 (Trudolyubov et al. (1996)), GRS 1739-278 (Borozdin et al. (1998)) and XTE J1748-288 (Revnivtsev, Trudolyubov & Borozdin (2000)). The sequence of typical spectra is presented in Fig. 2. One can see that during the first observation the soft component is not very strong (May 27, 1998), the second spectrum already has much more dominant soft component (June 1, 1998) and the third spectrum shows the weakening of the power law component (July 7, 1998). Broad band spectrum, averaged over observations $`\mathrm{\#}2\mathrm{\#}17`$ is presented in Fig. 3.

As it is clearly seen from Fig. 4, the average power density spectrum of the source is dominated by Very Low Frequency noise (VLFN) component, reasonably approximated by a simple power law model ($`Pf^\alpha `$) with a slope $`\alpha =1.08\pm 0.06`$ ($`0.0010.5`$ Hz frequency range). We have not detected any statistically significant source flux variability at the frequencies higher than $``$ 1 Hz. The 2 $`\sigma `$ upper limit on the possible Lorentzian component at $`f5001000`$ Hz is approximately 2% for $`Q=10`$ and 1% for $`Q=1`$. In Fig. 5 the fractional rms amplitude of variability is shown as a function of the photon energy. In spite of a poor statistics (very low amplitude of X-ray variability), one can see that there is an indication on the rise of the rms with energy. Such dependence is quite typical for the emission of black hole candidates in the Very High State (VHS)(see discussion of such behavior in Churazov, Gilfanov & Revnivtsev (2000)).

## 4 Discussion

The the properties of the outburst of XTE J2012+381 are very similar to what of the X-ray Novae. It has a FRED (fast-rise-exponential-decay) light curve, with secondary maxima, as typical for X-ray transients (Chen, Shrader & Livio (1997)).

The spectrum of the source can be well described by standard two-component model which is typical for black hole candidate X-ray binaries (Tanaka & Lewin (1995)). The evolution of the energy spectrum of XTE J2012+381 can be characterized as a gradual increase of the contribution of the soft thermal component as the overall X-ray flux decreases.

Our calculations show that the value of model parameter $`R_{in}\mathrm{cos}i`$, inferred from the soft spectral component approximation is nearly constant during the decay phase of the outburst (Table 2) in general agreement with earlier results of Zhang et al. (1998). However we believe that it is premature to assume a constant radius of inner accretion disk during the outburst. 222Note that no corrections to the electron scattering and general relativity were made in the soft spectral component model (Shakura & Sunyaev (1973); Shimura & Takahara (1995)). This model assumes incorrect radial dependence of the disk effective temperature within $`5`$ gravitational radii. Thus the inferred value of the effective radius $`\mathrm{R}_{\mathrm{in}}`$ should not be treated as the actual size of the optically thick emitting region.

The PDS of XTE J2012+381 is also very typical for the VHS of black hole candidates with dominant soft spectral component. The energy dependence of the fractional rms provided another hint at the universal correlation between the amplitude of the variability and relative strength of the hard spectral component.

###### Acknowledgements.

This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center.


text:
# Abstract

## Abstract

The preliminary results are presented of the experimental determination of threshold reaction rates in experimental samples made of Al, Co, Bi, In, Au, to name but a few, placed both inside and on the surface of extended thick W-Na target irradiated with 0.8 GeV protons. The target consists of 26 alternating discs each 150 mm in diameter: 6 tungsten discs are 20 mm thick, 7 tungsten discs are 40 mm thick, and 13 sodium discs are 40 mm thick. The relative position of discs is matched with the aim of flattening the neutron field along the target surface. The comparison is made of the measured rates with results of their simulation using the LAHET and KASKAD-S codes, and the ENDF/B6, MENDL2, MENDL2P, SADKO-2, and ABBN-93 databases. The results are of interest both in terms of integral data collection and to test the up-to-the-date predictive power of the codes applied in designing of hybrid Accelerator Driven System (ADS) systems that use tungsten targets cooled with sodium.

## Foreword

Quantitative information on interaction of accelerated protons with different targets is necessary to design ADS. Application of hadron-nucleus process simulations should be tested by special experiments in which irradiation conditions, target material composition and location approximate the design type to the limit. There are several projects, , for example, where tungsten cooled by sodium is considered to be used as target material. This served as the basis for conducting experiments with micromodels of such targets. Comparison of experimental data obtained on such micromodels with corresponding calculated values will give us valuable information both for modifications of codes and databases and for assessment the calculation accuracy of the target part of the relevant ADS facility designs.

## Experiment plan

The target with alternating adherent tungsten and sodium discs was designed in our micromodel. The location of these discs, which was specially chosen, facilitates the maximum flattening of the neutron field along the target. Experimental samples made of Al, Co, In, Au, Bi, 63Cu, 65Cu, 93Nb, 64Zn, 19F (CF2), 12C, Ta, and Tm manufactured by punching the corresponding foils or by molding fine powders, 10.5 mm in diameter and 0.1–0.3 mm thick, were placed inside the target and on its surface.

The layout of tungsten and sodium discs is shown in Fig. 1. All discs are 150 mm in diameter, disc thickness and sequence is listed in Table 1. The tungsten discs have special design providing insertion of special bars with round recesses for experimental samples to be placed inside the target. The discs are located on a special adjustment table which provides alignment of target and proton beam axes with an accuracy of the order of 1 mm. The proton beam size examination has been performed by using in tentative exposures to radiation aluminum cut foils and Polaroid film.

W (97.5%), Ni (1.75%), Fe (0.75%), and less than 0.2 % of impurities are incorporated in tungsten discs. The average density of the tungsten discs is of 18.6 g/cm3. Sodium discs represent metallic sodium placed into a steel container with 0.4 mm thick walls. Impurities content in Na is less than 0.02%.

Target irradiation was performed with 0.8 GeV protons over a period of 10 hrs at an average intensity equal to 4.8$``$1010 p/cm$`{}_{}{}^{2}`$pulse using the ITEP synchrotron. Pulse repetition rate is of 15 pulses per minute. Changes in the proton beam intensity over the irradiation period is presented in Fig. 2. After a short decay lag, experimental samples were extracted from the target’s surface and inside volume and packaged into labeled polyethylene packages. Subsequent gamma spectra measurements were made using several spectrometers. The absolute value for different threshold reaction rates were determined using the PCNUDAT decay database after gamma-spectra processing with the GENIE2000 code (see Table 2).

## Simulation of reaction rates

To simulate the measured reaction rates, the LAHET Code System (LCS)  was used which involves the LAHET code for simulating hadron-nuclei interactions, the HMCNP code for simulating neutron transport at energies below 20 MeV, and the PHT code for simulating hard photon transport. Both neutron and proton spectra were calculated in all points where experimental samples are located.

Another code system applied to simulate the experiment is the KASKAD-S code  which uses a discrete ordinate algorithm for coupled charges/neutral particle transport calculations in 2D-pencil beam problems. The multigroup cross-section library SADKO-2 for nucleon-meson cascade calculations coupled with the CONSYST/ABBN-93 neutron and gamma-ray cross-section libraries below 20 MeV is used.

Reaction cross sections for neutrons with energies up to 100 MeV and protons with energies up to 200 MeV were taken from MENDL2  and MENDL2P  libraries, respectively. Reaction rates are obtained via integral product of spectra and cross sections ($`R=\varphi (E)\sigma (E)𝑑E`$). The discrepancy between calculated and experimental results was estimated using the root-mean-square discrepancy factor, $`<`$F$`>`$, defined in .

## Results

The preliminary processing of gamma spectra from Al samples placed on the target outside surface and in the target center makes it possible to estimate the generation rates of 24Na and 27Mg. Experimental and calculated values obtained using the LAHET code are presented in Fig. 3. Measurements of gamma-spectra for all the samples irradiated are still being in progress and all the reaction rates listed in Table 2 are expected to be determined after the completion of gamma-spectra processing.

## Convergence of calculated and experimental values

The rates for 24Na in Al foils predicted by LAHET are in satisfactory agreement with experimental values (except the first and last points on the surface), as shown in Fig. 3. Whereas, the predicted rates for 27Mg generation on the target surface are systematically overestimated compared to the measured values (by a factor of 1.6, on average), except the first point. At the same time, the 27Mg production rates predicted by KASKAD-S are underestimated on the average by a factor of 1.6. The observed discrepancies would be better understood only after the complete processing of gamma spectra for the rest of samples and determining both (n,p) – purely “neutron” and (p,n) – purely “proton” reactions.

## Acknowledgement

The work was carried out under the ISTC Project#1145 and was supported by the JAERI (Japan) and, partially, by the U. S. Department of Energy.