text:
# Dynamics of a trapped 2D Bose-Einstein condensate with periodically and randomly varying atomic scattering length

## I Introduction

The problem of the oscillations of two and three dimensional trapped Bose-Einstein condensates (BEC) have attracted a great deal of recent attention, particularly the problem where the trapping potential is temporally modulated. The main physical motivation for considering such a problem is the experiments of Jin, et. al., which observed resonant response of the BEC oscillations to periodically perturbed traps. Some theoretical investigations of the problem have included the work of Castin and Dum, in which the resonances in the width oscillations have been analyzed using scaling theory; the work of Garcia, in which it was argued that the oscillations of 2D condensate can be described by the parametric resonances of linear oscillator; and the papers by Pitaevski and Rosch and Kagan et. al., where it was shown that a 2D condensate exhibits a harmonic mode with frequency $`2\omega _0`$ when the harmonic trap is driven with frequency $`\omega _0`$. Note that the related problem of a soliton interacting with an impurity in 2D molecular crystals had been studied earlier by Gadidei et al..

The excitation of resonant oscillations is also interesting as a possible mechanism for the stabilization of a BEC with attractive interactions. Such a possibility has been demonstrated for the BEC with negative atomic scattering length under temporally modulated trap potential if the surface mode is exited . In this work a variational approach is used, and the trap was anisotropic.

It is also interesting to consider the response of a condensate to temporally periodic variations in the atomic scattering length. Such a variation of the atomic scattering length can be achieved experimentally by varying the magnetic field  or using optically induced Feshbach resonances . The case where the scattering length varies monotonically in time, in particular the case where the sign changes from repulsive to attractive, has been studied recently by Dalfovo et. al, and Fedichev, et. al.. Similarly the effect of a periodic variation of atomic scattering length on the tunneling between two condensates in the double-well trap, a resonant tunneling, has been studied by Abdullaev and Krankel.

These observations suggest the possibility of new phenomenon in Bose-Einstein condensate width oscillations when the atomic scattering length $`a_s`$ is allowed to fluctuate about some mean value. In this paper we study the influence of both periodic and random fluctuations of the atomic scattering length on the condensate dynamics. As noted in , at temperatures $`Tng,g=4\pi \mathrm{}^2a_s/m`$, (where $`n`$ is the gas density, $`a_s`$ is the atomic scattering length, $`m`$ is the atom mass), only the condensate evolution is pronounced. The perturbation of the thermal cloud is small at these temperatures, and in this case the eigenfrequencies of small oscillations of the thermal cloud are close to those of the condensate.

The specific problem that we consider is the 2D Gross-Pitaevskii equation with periodic or random variations in the coefficient of the nonlinear term. This problem is interesting not only for BEC but also for nonlinear optics too - for example the dynamics of optical beams in nonlinear layered waveguides .

The structure of the remainder of the article is follows: In Sect.2 we describe the model and derive the reduced ODE system to describe the dynamics of the condensate width. In Sect.3 we analyze the effect of time-period perturbations on the width dynamics using the action-angle variables for the reduced ODEs. We explore the resonances in condensate oscillations in Sect.4, and the dynamics of width oscillations under fluctuating trap potential and atomic scattering length in last section.

## II The Reduced ODE Model

It is well-known that wavefunction for a 2D Bose-Einstein condensate in a trap potential $`V(r)`$ is described by the Gross-Pitaevskii equation:

$$i\mathrm{}\psi _t=\frac{\mathrm{}^2}{2m}\mathrm{\Delta }\psi +V_{tr}(r)\psi +g(t)|\psi |^2\psi .$$
(1)
Here $`V_{tr}(r)=m\omega ^2r^2/2`$ is the trap potential and $`g(t)=4\pi \mathrm{}^2a_s/m`$, with $`a_s`$ being the atomic scattering length. We will assume that the time-dependent scattering length is constant to leading order, with

$$a_s=a_0(1+ϵ_0(t)).$$

In this paper we analyze the case where $`ϵ_0(t)`$, representing the small fluctuations, is a periodic function of time as well as the case where $`ϵ_0(t)`$ is a mean zero white noise random process: $`<ϵ_0(t)ϵ_0(t^{})>=2\sigma ^2\delta (tt^{})`$.

Under the scaling $`\tau =t\omega `$,$`x^{}=x\sqrt{\mathrm{}/(2m\omega )}`$, $`u=\psi \sqrt{4\pi a_0/\omega }`$ the Gross-Pitaevskii equation can be written in the dimensionless form

$$iu_\tau +\mathrm{\Delta }u\frac{r^2}{2}u+s|u|^2u=0,$$
(2)
where $`s>0,s<0`$ correspond to attractive or repulsive interactions between atoms in the condensate respectively.

A number of methods can be employed to investigate the dynamics of the Bose-Einstein condensate under temporal variations of the atomic scattering length. One of the simplest to apply is the averaged Lagrangian approach. According to this method we take the Gaussian anzatz for the field .

$$u(r,\tau )=A(\tau )\mathrm{exp}(\frac{r^2}{2a(\tau )^2}+\frac{ib(\tau )r^2}{2}+i\varphi (\tau )),$$
(3)
To derive the equations for the wavepacket parameters $`A(\tau ),a(\tau ),b(\tau ),\varphi (\tau )`$ one should calculate the averaged Lagrangian

$$\overline{L}(\tau )=r𝑑rL(r,\tau ).$$

For the Gaussian ansatz in Eq.(3) the averaged Lagrangian is given by

$$\overline{L}=\frac{\pi }{2}A^2a^2(a^2b_\tau +2\varphi _\tau +\frac{2}{a^2}+2a^2b^2+a^2\frac{s(\tau )A^2}{2}).$$
(4)
The Euler-Lagrange equations for the functional $`\overline{L}(\tau )`$ lead to the following equations for the the width $`a(\tau )`$ and chirp $`b(\tau )`$.

$`a_\tau =2ab`$ (5)
$`b_\tau ={\displaystyle \frac{2s(\tau )N}{a^4}}2b^21,`$ (6)
with $`N=\frac{1}{2\pi }|u|^2d^2r=\frac{a^2A^2}{2}.`$ It is easy to see that the above is a Hamiltonian system, with the width $`a`$ playing the role of a position variable, and the chirp $`b`$ the conjugate momentum. One can eliminate $`b`$ from the above system which leads to the following evolution equation for $`a`$:

$$a_{tt}+a=\frac{2s(t)N}{a^3},$$
(7)
where $`t=\sqrt{2}\tau `$. For the attractive problem with constant atomic scattering length $`s(t)=1`$ the variational approach predicts the critical threshold for collapse $`N_c=2`$, compared with the the exact value $`N_c=1.862`$.

We comment that, if the initial condition is close to that of the ground state (Townes) soliton, $`|NN_c|ϵ1`$, then the modulation theory of Fibich and Papanicolaou leads naturally to a modulation equation of the same form as the above, where again the constant $`2`$ is replaced by the critical value $`1.862`$. The qualitative types of behavior exhibited by the above equation do not, of course, depend on the exact values of the constants, since all such constants can be scaled out.

It is useful to rewrite Eq.(7) in the form

$$a_{tt}+a=\frac{Q(t)}{a^3},$$
(8)
where $`Q(t)=Q(1+ϵ(t)),Q=2N,ϵ=Nϵ_0/(2N)`$. In the next section we construct the action-angle variables for the above problem.

## III Perturbation theory in the action-angle variables

We begin our analysis of the Eq.(7) by constructing the action-angle variables. Of course any two dimensional Hamiltonian system can always be reduced to quadrature, but the above system is particularly nice because the action-angle variables can be expressed in terms of elementary functions. First we note that Eq.(7) is Hamiltonian:

$$H(a_t,a)=\frac{a_t^2}{2}+U(a),U(a)=\frac{a^2}{2}+\frac{Q}{2a^2}.$$
(9)
For $`N<2,Q>0`$ and the motion of the effective particle is bounded, and bounded away from $`0`$ due to the “angular momentum” barrier $`\frac{Q}{2a^2}`$. If the energy is $`E`$ then the width oscillates between $`a_{min}=\sqrt{E\sqrt{E^2Q}}`$ and $`a_{max}=\sqrt{E+\sqrt{E^2Q}}`$. The minimum of the potential $`U`$ occurs at $`a_c=Q^{1/4}`$, with a minimum energy of $`E=U_c=\sqrt{Q}`$. When $`N>2,Q<0`$ there is no local minimum, and there exist solutions for which $`a0`$, corresponding to collapse of a condensate (see Fig.1).

Since the Hamiltonian is conserved, $`H=E`$, we have

$`{\displaystyle \frac{da}{\sqrt{2Ea^2Qa^2}}}={\displaystyle 𝑑t}`$

which can be integrated up to the following solution for $`Q>0:`$

$$a(t)=\sqrt{E+\sqrt{E^2Q}\mathrm{sin}(2t+\psi _0)}.$$
(10)
The action variable for a Hamiltonian system is defined to be

$$J=\frac{1}{2\pi }p𝑑q,$$
(11)
where $`q`$ is the position variable and $`p`$ the conjugate momentum. For the above oscillator the position variable is $`a`$, with momentum $`p=\sqrt{2Ea^2a^4Q}/a.`$ The integration is easily done via contour integration in the complex $`a`$ plane leading to the following expression for the action vaiable:

$$J=\frac{1}{2}(E\sqrt{Q}).$$
(12)
The bottom of the potential well, $`E=\sqrt{Q}`$, corresponds to $`J=0`$. Since the energy is linear in the action, $`E=2J+\sqrt{Q},`$ the frequency of the unperturbed oscillations is constant

$$\omega =\frac{dE}{dJ}=2.$$
(13)
This can, of course, also be seen directly from the solution given in Eq. (10).

The Hamiltonian for the perturbed problem is given by

$$H=p^2/2+U(a)+ϵ(t)V=p^2/2+U(a)+\frac{Qϵ(t)}{2a^2}$$

, with $`V(a)`$ the perturbation Hamiltonian. In the action-angle coordinates $`J,\theta `$ the perturbation becomes

$$V_\theta =\frac{Q\sqrt{E^2Q}\mathrm{cos}(\theta )}{2(E+\sqrt{E^2Q}\mathrm{sin}(\theta ))^2}.$$

and the perturbed equations of motion are given by

$`{\displaystyle \frac{dJ}{dt}}=ϵ(t){\displaystyle \frac{V}{\theta }},`$ (14)
$`{\displaystyle \frac{d\theta }{dt}}=2+ϵ(t){\displaystyle \frac{V}{J}}.`$ (15)
## IV Resonances of BEC oscillations under periodically varying atomic scattering length

In this section we consider perturbations $`ϵ(t)`$ which are periodic in time. For simplicity we discuss the case $`ϵ(t)=\mathrm{sin}(\mathrm{\Omega }t)`$, though this can easily be generalized. To analyze the resonances of a BEC under such perturbations we use the multiscale expansion method. We introduce the slow time $`T=ϵt`$ and assume a multiple-scales ansatz of the form

$`\theta =\theta ^{(0)}(t,T)+ϵ\theta ^{(1)}(t,T)+\mathrm{},`$ (16)
$`J=J^{(0)}(t,T)+ϵJ^{(1)}+ϵJ^{(1)}(t,T)+\mathrm{}`$ (17)
The solution to the evolution on the fast scale is obviously $`\theta ^{(0)}=2t+\mathrm{\Phi },J=\mathrm{constant}`$. At the next order we find the following evolutions for the action $`J`$ and the slow angle $`\mathrm{\Phi }`$ on the slow scale:

$`{\displaystyle \frac{dJ}{dT}}={\displaystyle \frac{\overline{H}}{\mathrm{\Phi }}}`$ (18)
$`{\displaystyle \frac{d\mathrm{\Phi }}{dT}}={\displaystyle \frac{\overline{H}}{J}},`$ (19)
where $`\overline{H}`$, the slow Hamiltonian, is given by

$$\overline{H}(J,\mathrm{\Phi })=\frac{1}{2\pi }_0^\pi \frac{\mathrm{sin}(\mathrm{\Omega }t)dt}{(2J+Q^{1/2}+2J^{1/2}(J+Q^{1/2})^{1/2}\mathrm{sin}(2t+\mathrm{\Phi }))}.$$

It is clear that we have a resonant response whenever $`\omega =2n`$, and the period of the perturbation is commensurate with the natural period of the condensate. Note that, even though this is a nonlinear oscillator the period is independent of the amplitude. The effect of this is that there is no ”detuning” from the resonant frequency as the amplitude increases. Of course when $`J`$ becomes large and the width oscillations become significant we no longer expect the averaged Lagrangian ODE’s to provide an accurate description of the dynamics of the condensate width.

The first (nontrivial) resonance occurs for $`\mathrm{\Omega }=2.`$ In this case the slow Hamiltonian is given by

$$\overline{H}(J,\mathrm{\Phi })=\frac{\mathrm{cos}(\mathrm{\Phi })J^{1/2}}{2(J+Q^{1/2})^{1/2}Q^{1/2}}.$$
(20)
and the equations of motion are

$`{\displaystyle \frac{dJ}{dt}}`$ $`=`$ $`{\displaystyle \frac{\overline{H}}{\mathrm{\Phi }}}={\displaystyle \frac{\mathrm{sin}(\mathrm{\Phi })J^{1/2}}{2(J+Q^{1/2})^{1/2}Q^{1/2}}},`$ (21)
$`{\displaystyle \frac{d\mathrm{\Phi }}{dt}}`$ $`=`$ $`{\displaystyle \frac{\overline{H}}{J}}={\displaystyle \frac{\mathrm{cos}(\mathrm{\Phi })}{4J^{1/2}(J+Q^{1/2})^{3/2}}}.`$ (22)
The phase plane for the oscillator is depicted in Fig. (2).

Note that the line $`\mathrm{\Phi }=\pi /2`$ is invariant under the dynamics, and along this line $`J`$ evolves according to

$$\frac{dJ}{dt}=\frac{J^{1/2}}{2(J+Q^{1/2})^{1/2}Q^{1/2}}.$$
(23)
This solution corresponds to a resonant driving, where the variations in the scattering length reinforce the width oscillations. It is clear that for large $`t`$ the action grows linearly, $`Jt/2Q^{1/2}`$ and thus amplitude of oscillations grows like $`a\sqrt{t}`$. It is also clear from the phase portrait that all orbits for which $`\mathrm{\Phi }\pi /2`$ are asymptotic to the invariant manifold $`\mathrm{\Phi }=\pi /2`$, so for generic initial conditions one expects that the width will grow like $`a\sqrt{t}.`$ There is, of course, also the solution $`\mathrm{\Phi }=\pi /2`$ in which the variations in the scattering length are anti-resonant with the variations in the width of the condensate, and act to damp the width oscillations. It is easy to see that when $`\mathrm{\Phi }=\pi /2`$ the action $`J`$ goes to zero in finite time. Since generic initial conditions are asymptotic to $`\mathrm{\Phi }=\pi /2`$ the $`\mathrm{\Phi }=\pi /2`$ solutions are unlikely to be observed experimentally, though it is possible that they could be realized with some kind of control.

Next resonance occurs at $`\mathrm{\Omega }=4`$. In this case the Hamiltonian is given by

$$\overline{H}=\frac{J\mathrm{sin}(2\mathrm{\Phi })}{(J+Q^{1/2})Q^{1/2}}.$$
(24)
Again it is easy to see that the line $`\mathrm{\Phi }=\pi /2`$ is invariant under the dynamics, corresponding to width oscillations which are in phase with the variation of the scattering length. Since the action evolves according to

$$\frac{dJ}{dT}=\frac{2J\mathrm{cos}(2\mathrm{\Phi })}{(J+Q^{1/2})Q^{1/2}}$$
(25)
when $`\mathrm{\Phi }=\pi /2`$ we have that the growth of the action $`J`$ is asymptotically linear - $`J2t/Q^{1/2}`$.

### A Numerical simulations

We have conducted some numerical simulations to test the validity of the calculations of the last section. We discretize the problem in the standard way, with time step $`\mathrm{\Delta }t`$ and spatial step $`h`$, so the $`u_j^k`$ approximates $`u(jh,k\mathrm{\Delta }t)`$. More specifically we approximate Eq.(2) with the following second order accuracte semi-implicit Crank-Nicholson scheme,

$`{\displaystyle \frac{i(u_j^{k+1}u_j^k)}{\mathrm{\Delta }t}}={\displaystyle \frac{1}{2h^2}}[(u_{j1}^{k+1}2u_j^{k+1}+u_{j+1}^{k+1})+(u_{j1}^k2u_j^k+u_{j+1}^k)]`$ (26)
$`{\displaystyle \frac{1}{4r_jh}}[(u_{j+1}^{k+1}u_{j1}^{k+1})+(u_{j+1}^ku_{j1}^k)]+{\displaystyle \frac{1}{2}}[{\displaystyle \frac{r_{j}^{}{}_{}{}^{2}}{2}}(1+ϵ(t))s|u_j^k|^2](u_j^k+u_j^{k+1}),`$ (27)
where $`ϵ(t)`$ is the perturbation term. In the numerical simulations, as in the analysis, the perturbation was chosen to be $`ϵ(t)=ϵ\mathrm{sin}(\mathrm{\Omega }t).`$ Eq. (26) represents a tridiagonal set of equations for unknowns $`u_{j1}^{k+1}`$, $`u_j^{k+1}`$ and $`u_{j+1}^{k+1}`$ $`[j=1,2\mathrm{}(N1)]`$ in a lattice of $`N`$ points, with the values of $`u_0^{k+1}`$ and $`u_{N+1}^{k+1}`$ being determined from the boundary conditions $`\frac{u}{r}_{r=0}=0`$ and $`u(r)_{r=\mathrm{}}0`$.

The set of algebraic equations (26) is solved by the vectorial sweep method. In actual calculations the typical space step $`h`$ ranged from 0.01 to 0.005 and time step $`\mathrm{\Delta }t`$ from 0.005 to 0.002 depending on the closeness of $`\mathrm{\Omega }`$ to the points of resonance.

The first experiment, shown in Fig. 3, depicts the solution of the Gross-Pitaevski equation when the driving frequency is slightly off resonance: the natural frequency of the width oscillations is $`\mathrm{\Omega }_0=2`$, while the perturbation has frequency $`\mathrm{\Omega }=1.9`$ and amplitude $`ϵ=0.1`$. In this plot, as well as all subsequent plots, the solid line represents the solution to the full PDE while the dotted line represents the solution to the reduced ODEs. We observe oscillations with beats, representing the superposition of low and high frequency oscillations. From simulations we seen the good agreement between the full PDE and the reduced ODE model.

Fig.4 represents numerical simulations of the ODE and PDE models at the resonance point $`\mathrm{\Omega }=\mathrm{\Omega }_0=2`$. The agreement is quite good for the period of the width oscillations, although there is clearly some discrepancy in the actual value of oscillation amplitude. Note that the same phenomenon has been observed in numerical simulations of the resonances in condensate oscillations under periodically varying trap potential . The graph in Fig.5 depicts the energy versus time for the ODE (dotted line) and PDE (solid line) simulations for the same values of parameters as in Fig.4. The agreement between simulations of the full PDE and ODE is very good for time less than 30 or so, though for times between $`30`$ and $`40`$ the oscillations in the energy of the PDE are smaller than the analogous oscillations of the ODE, probably due to radiative damping. As was argued earlier the variational approach is unlikely to be valid when the oscillations have large amplitude and other effects, such as radiative damping, become important.

Fig.6 depicts the oscillations of the square width of the condensate at the $`2:1`$ resonance $`\mathrm{\Omega }=2\mathrm{\Omega }_0`$, where the frequency of the perturbation is twice the natural period of the width oscillation. As in previous case the frequencies agree very well, but the amplitude of oscillations is larger for PDE in comparison with ODE. This discrepancy grows with time. In this case it seems clear that the Gaussian anzatz does not correctly capture the behavior of the underlying PDE. The rate of growth of the energy, which is not depicted, is still linear, in agreement with the theoretical estimates.

## V Evolution of BEC under random fluctuations

The action-angle formulation also provides a suitable framework for the analysis of condensate oscillations under random (white-noise) perturbations. One interesting quantity which can be calculated is the mean time to achieve a given distortion under a stochastic perturbation. From the ODE point of view this problem is equivalent to the problem of the mean time to achieve the given level of the oscillations amplitude in the effective potential.

The starting point for this calculation is Eqn(14,15), where $`ϵ(t)`$ is taken to be a white noise process $`ϵ(t)dt=dB`$, with $`dB`$ the increment of a Brownian motion $`B`$. This pair of stochastic ordinary differential equations has an associated Fokker-Planck equation for $`P(J,\theta ,t)`$, the probability of having values $`J,\theta `$ at time $`t`$. In the weak noise limit one can do a straightforward multiple scale expansion on this Fokker-Planck equation. Upon doing so and averaging over the fast angle variable one is lead to the following equation for $`P(J,t)`$:

$$\frac{P}{t}=\sigma ^2\frac{}{J}(A(J)\frac{P}{J})P(J,0)=\delta (JJ_0),$$
(28)
where $`A(J)`$ the average diffusivity for diffusion across energy levels is given by

$$A(J)=\frac{1}{2\pi }_0^{2\pi }V_\theta ^2𝑑\theta .$$

For details of this calculation see the paper of Abdullaev, et. al This one-dimensional diffusion can be analyzed in some detail. For instance the mean time to reach action $`J`$ staring from $`J_0`$ is given by

$$<t>_{J_0,J}=_{J_0}^J\frac{JdJ}{A(J)}.$$
(29)
Note that if $`A(J)`$ grows sufficiently rapidly, so that the above integral converges, then the mean time to random walk to $`J=\mathrm{}`$ is actually finite.

### A Fluctuating trap potential

One physically interesting random perturbation is when the strength of the trap is allowed to fluctuate randomly. Since the trapping potential is imposed optically, by a laser,it is important to take into account the fluctuations of the effective trap potential due to the fluctuations of the laser field intensity. The fluctuations of the laser field intensity lead to random variations of the frequency of the effective harmonic trap $`\omega ^2=\omega _0^2(1+ϵ(t))`$, where $`ϵ(t)`$ is the white noise process . The function $`ϵ(t)`$ should be thought of as the fluctuations of the laser intensity around its mean value $`E_0`$:

$$ϵ(t)=\frac{E(t)E_0}{E_0}$$

This problem was studied in  using a moment expansion, though without numerical simulations of the full stochastic GP equation. Here we consider this problem both as an illustration of effectiveness of our technique, and to provide numerics for this problem. In this case the perturbation Hamiltonian is

$$V(a)=\frac{1}{2}(E+\sqrt{E^2Q}\mathrm{sin}\theta ).$$

Then from Eqn.(V) the diffusivity is given by

$$A(J)=\frac{J(J+\sqrt{Q})}{2}.$$

Calculating the mean exit time $`<t>`$ we find the energy as a function of the mean exit time

$$E=2\sqrt{Q}(e^{\sigma ^2<t>/2}\frac{1}{2}).$$
(30)
While the mean exit time is $`\mathrm{n}\mathrm{o}\mathrm{t}`$ the same as the physical time in the limit of weak noise we expect that, to leading order in the noise parameter, the mean energy should have the same dependence on $`t`$. In Fig.7 we present the results of comparison of the theory and the numerical simulations of the full stochastic NLS equation (2) for $`\sigma =0.04,N=1.0`$. For the PDE we have averaged over 50 realizations, and for ODE we have averaged over 1000 realizations. The figure shows good agreement between predictions of the theory and numerical simulations. For large times the solution to the ODE overestimates the actual energy. This discrepancy is probably due to effects such as dispersive radiation, excitation of higher modes, etc. which are difficult to incorporate into the variational ansatz.

From the above we can easily estimate the mean exit time from the potential. In optical traps the depth of the effective potential is $`U_0=2E_R=2\mathrm{}^2k^2/2m`$, with $`k=2\pi /\lambda `$, where $`\lambda `$ is the wavelength, The condensate lifetime can be estimated as the time when become $`EU_0`$. Then

$$t_{BEC}\frac{2}{\sigma ^2}\mathrm{ln}(\frac{U_0}{2\sqrt{Q}}+\frac{1}{2}).$$

Substituting the typical values for the laser intensity fluctuations we find that the mean exit time is of the order of seconds, $`t_{BEC}sec`$.

### B Fluctuating atomic scattering length

In the case where the atomic scattering length is allowed to fluctuate we find that the perturbation Hamiltonian is given by

$$V_\theta =\frac{Q\sqrt{E^2Q}\mathrm{cos}(\theta )}{2(E+\sqrt{E^2Q}\mathrm{sin}(\theta ))^2}.$$

From this is follows that the the effective diffusivity is

$$A(J)=\frac{Q^2(E^2Q)}{8\pi }_0^{2\pi }\frac{\mathrm{cos}^2(\theta )}{(E+\sqrt{E^2Q}\mathrm{sin}(\theta ))^4}𝑑\theta =\frac{E(E^2Q)}{8\sqrt{Q}},$$
(31)
and expressing the energy $`E`$ in terms of the action $`J`$ we find

$$A(J)=\frac{J(2J+\sqrt{Q})(J+\sqrt{Q})}{2\sqrt{Q}}.$$

Note that the expected time to random walk to infinity is actually finite, since $`A(J)J^3`$. Substituting this expression for into Eq.(29) we obtain for the time to pass from the bottom of the potential well where the action $`J_0=0`$ to the state with the action $`J`$

$$J=\frac{\sqrt{Q}(e^y1)}{2e^y}y=\frac{\sigma ^2t}{2},$$
(32)
or in the terms of the total energy

$$E=\frac{\sqrt{Q}e^y}{2e^y}.$$
(33)
From the above it is easy to see that the expected time for the width of the condensate to grow to infinity - the mean time for the condensate to break up - is given by

$$t^{}=\frac{2\mathrm{ln}2}{\sigma ^2}.$$
(34)
In Fig. 8 we compare the theoretical expression Eq.(33) with numerical simulations of the full stochastic GP equation with fluctuating scattering length in the case where $`\sigma =0.04`$ and $`N=1.0`$. The energy has been numerically calculated for the localized part of the condensate wavefunction. We observe good agreement between the theory and numerical simulations for $`t80`$.

## VI Conclusions

We have considered the oscillations of a 2D BEC with radial symmetry under periodic and random modulations of the atomic scattering length, as well as random fluctuations of trap potential. We have calculated the position of resonances and the energy growth using a reduced ODE for the condensate width. We have also confirmed the analytical predictions with numerical simulations of the 2D Gross-Pitaevskii equation. In the resonant case the frequency of oscillations agrees very well with the predictions of ODE, though the amplitude shows some discrepancy at large times that is probably due to the approximate nature of the variational approach. For the random modulations of the trap potential and the atomic scattering length we have calculated the mean exit time corresponding to the time for the amplitude of oscillations to exceed a given value and estimated the magnitude of the escape time for real experiments with 2D Bose-Einstein condensates.

## VII Acknowledgements

We acknowledge partial support by the US Civilian Research and Development Foundation (Award ZM2-2095) and by the Uzbekistan Fundamental Science Support Foundation (Award 15-02).


format: text

text:
# Field-Induced Magnetic Ordering in the Quantum Spin System KCuCl3

## I Introduction

Recently, the physics of coupled antiferromagnetic spin-dimer systems has been attracting considerable attention from the viewpoints of magnetic excitations , the magnetization processes including magnetization plateaus  and field-induced three-dimensional (3D) magnetic ordering. This study is concerned with the magnetization process of and the field-induced 3D magnetic ordering in KCuCl3, which is an $`S=\frac{1}{2}`$ coupled antiferromagnetic spin-dimer system and has a singlet ground state with an excitation gap (spin gap) $`\mathrm{\Delta }/k_\mathrm{B}=31`$ K .

KCuCl3 has a monoclinic structure (space group $`P2_1/c`$) . The crystal structure is composed of planar dimers of Cu2Cl6. The dimers are stacked on top of one another to form infinite double chains parallel to the crystallographic $`a`$-axis. These double chains are located at the corners and center of the unit cell in the $`bc`$ plane as shown in Fig. 1. From the structural point of view, KCuCl3 was first assumed to be a double chain spin system . The magnetic excitations in KCuCl3 have been extensively investigated through neutron inelastic scattering  and ESR measurements . The effective dimer approximation  and the cluster series expansion  were successfully applied to analyze the dispersion relations obtained. From these analyses, it was found that the origin of the spin gap is the strong antiferromagnetic interaction in the chemical dimer Cu2Cl6, and that the neighboring dimers are weakly coupled not only along the double chain, but also in the $`(1,0,2)`$ plane, in which the hole orbitals of Cu2+ spread. Consequently, KCuCl3 was characterized as a weakly and three-dimensionally coupled spin-dimer system.

KCuCl3 differs from isostructural TlCuCl3 and NH4CuCl3 in magnetic character. TlCuCl3 also has a gapped ground state, however, the gap is considerably suppressed due to strong interdimer interactions . On the other hand, NH4CuCl3 has a gapless ground state at zero magnetic field . The remarkable feature of NH4CuCl3 is the existence of magnetization plateaus at one-quarter and three-quarters of the saturation magnetization .

When a magnetic field is applied in the spin gap system, the gap $`\mathrm{\Delta }`$ is suppressed and closes at $`H_\mathrm{g}=\mathrm{\Delta }/g\mu _\mathrm{B}`$. For $`H>H_\mathrm{g}`$ the system can undergo magnetic ordering due to 3D interactions with decreasing temperature. Such field-induced magnetic ordering was observed in the isostructural TlCuCl3 . For the field-induced magnetic ordering, two characteristic features which cannot be explained by the mean-field approach from the real space  have been observed, irrespective of the applied field direction, when the applied field $`H`$ is slightly higher than $`H_\mathrm{g}`$. One is that the magnetization has a cusplike minimum at the transition temperature $`T_\mathrm{N}`$. The other is that the phase boundary between the paramagnetic phase and the ordered phase can be expressed by the power law

$$[H_\mathrm{N}(T)H_\mathrm{g}]T^\varphi ,$$
(1)
with $`\varphi =2.1(1)`$. Here, $`H_\mathrm{N}(T)`$ denotes the transition field at temperature $`T`$. These features can be understood in terms of the Bose-Einstein condensation (BEC) of excited triplets (magnons) .

When $`H`$ is slightly higher than $`H_\mathrm{g}`$, the density of magnon $`n`$ which corresponds to the uniform magnetization is small. In this case, the system can be mapped onto the dilute magnon model. Using the Hartree-Fock approximation, Nikuni et al.  demonstrated that the magnetization has the cusplike minimum at the transition temperature where the BEC of magnons occurs and that the phase boundary can be described by eq. (1) with $`\varphi =\frac{3}{2}`$. Although the theoretical exponent is somewhat smaller than the experimental value $`\varphi =2.1`$, the magnon BEC theory gives a good description of the experimental results for TlCuCl3.

The magnetic anisotropy is negligible in KCuCl3  as it is in TlCuCl3 . Therefore, the previously mentioned characteristic features for the magnon BEC should also be observed in KCuCl3, because the magnon BEC is universal in the isotropic spin gap system. With this motivation, we carried out magnetization measurements in static high magnetic fields above the gap field $`H_\mathrm{g}20`$ T using a hybrid magnet. Since the entire magnetization curve in KCuCl3 has not been observed, we also performed magnetization measurements in pulsed high magnetic fields of up to 60 T. As shown below, the saturation of the magnetization was observed at $`H_\mathrm{s}50`$ T. From the saturation field $`H_\mathrm{s}`$, it was found that the exchange parameters determined from the analyses of the dispersion should be reduced. This suggests the importance of the renormalization effect in the magnetic excitations in KCuCl3.

## II Experimental Details

The single crystals of KCuCl3 were prepared by the Bridgman method. The details of sample preparation were almost the same as those for TlCuCl3 . For KCuCl3, the temperature of the center of the furnace was set at 500C. Single crystals of size $`0.55`$ cm3 were obtained. The crystals are easily cleaved along the $`(1,0,\overline{2})`$ plane. The second cleavage plane is $`(0,1,0)`$. These cleavage planes are perpendicular to each other.

The high-field magnetization process in KCuCl3 was measured at $`T=1.3`$ K using an induction method with a multilayer pulse magnet at the Research Center for Materials Science at Extreme Conditions (KYOKUGEN), Osaka University. The magnetic fields were applied perpendicularly to the $`(1,0,\overline{2})`$ plane and along the $`[2,0,1]`$ direction which is parallel to both cleavage $`(1,0,\overline{2})`$ and $`(0,1,0)`$ planes. We used a single crystal in the measurement for $`H[2,0,1]`$. For $`H(1,0,\overline{2})`$ measurement, we stacked several single crystals along the cylindrical axis of the sample holder.

The temperature dependence of the magnetization of KCuCl3 was measured in static magnetic fields of up to 30 T by a sample-extraction method with a 40 T-class hybrid magnet at the National Institute for Materials Science, Tsukuba. The magnetic fields were applied perpendicularly to the cleavage $`(0,1,0)`$ and $`(1,0,\overline{2})`$ planes, and along the $`[2,0,1]`$ direction. These field directions were perpendicular to one another.

## III Results and Discussions

### A Magnetization Process

Figure 2 shows the magnetization curves of KCuCl3 measured at $`T=1.3`$ K for $`H(1,0,\overline{2})`$ and $`H[2,0,1]`$. Here, the magnetization curves are normalized with the $`g`$-factors, i.e., $`g=2.26`$ for $`H(1,0,\overline{2})`$ and $`g=2.04`$ for $`H[2,0,1]`$, which were obtained by ESR measurements. The magnetization curves for two different field directions almost coincide with each other when normalized by the $`g`$-factor, although the data for $`H[2,0,1]`$ are somewhat scattered. This implies that the difference between the transition fields $`H_\mathrm{g}`$ and the saturation fields $`H_\mathrm{s}`$ for $`H(1,0,\overline{2})`$ and $`H[2,0,1]`$ is attributed to the anisotropy of the $`g`$-factor, and that the magnetic anisotropy is negligible, as previously concluded . Due to the spin gap, the magnetization is almost zero up to the transition field $`(g/2)H_\mathrm{g}=23`$ T, and increases rapidly and monotonically, and then saturates at $`(g/2)H_\mathrm{s}=54`$ T. As shown in Fig. 2, the slope of the magnetization in the vicinity of $`H_\mathrm{g}`$ and $`H_\mathrm{s}`$ is steeper than that in the intermediate field region. This behavior arises due to the quantum fluctuation.

In the magnetization slope region between $`H_\mathrm{g}`$ and $`H_\mathrm{s}`$, transverse spin components have a long-range order. Since the lowest magnetic excitations occur at $`=(0,0,1)`$ and its equivalent reciprocal lattice points, we assume that the spin ordering equivalent to that observed in TlCuCl3 is realized in the ordered state (see Fig. 3). For the exchange interactions shown in Fig. 3, we use the notation given in ref. : the main intradimer exchange is denoted as $`J`$. $`J_{(lmn)}`$ and $`J_{(lmn)}^{}`$ denote the exchange interactions between dimers separated by a lattice vector $`l𝕒+m𝕓+n𝕔`$. Applying the mean-field approximation to the interdimer interactions, we calculate the magnetization curve for the ground state. The details of the calculation are given in the Appendix.

In Table I, we list the exchange parameters in KCuCl3 as determined by Müller and Mikeska . They applied a cluster series expansion to analyze the dispersion relations observed in KCuCl3 . Their theory describes the experimental results, and the individual interdimer exchange parameters which cannot be obtained by the effective dimer approximation were determined. However, they used the same intradimer exchange interaction, $`J=4.25`$ meV, as that obtained from the effective dimer approximation .

The dotted line in Fig. 4 is the magnetization curve calculated with the exchange parameters obtained by Müller and Mikeska . We also plotted the experimental result obtained at $`T=1.3`$ K for $`H(1,0,\overline{2})`$. The calculated gap field $`H_\mathrm{g}`$ coincides with the experimental value, while the calculated saturation field $`H_\mathrm{s}`$ is considerably larger than the experimental value. We infer that the saturation field of eq. (A9) in the Appendix is close to the rigorous solution, because even for the $`S=\frac{1}{2}`$ 1D antiferromagnet, the classical calculation gives the same saturation field as the rigorous solution . Therefore, we have to reduce the exchange parameters. However, it is difficult to adjust all of the exchange parameters, because we have only two equations for $`H_\mathrm{g}`$ and $`H_\mathrm{s}`$. Thus, we make a uniform reduction for all of the interdimer exchange parameters. The thin solid line in Fig. 4 is the result calculated with the parameters listed in Table I. The intradimer interaction $`J`$ and the interdimer interaction $`J_{(lmn)}`$ or $`J_{(lmn)}^{}`$ are reduced by factors 0.9 and 0.8, respectively. Both the calculated gap field $`H_\mathrm{g}`$ and saturation field $`H_\mathrm{s}`$ coincide with experimental values. Hence, we suggest that the renormalization effect in the magnetic excitations cannot be neglected, although the present system is 3D. The renormalization factor $`\pi /2`$ in the $`S=\frac{1}{2}`$ 1D antiferromagnet is well known . The difference between the calculated and experimental results for $`H_\mathrm{g}<H<H_\mathrm{s}`$ is due to the quantum fluctuation, which is not taken into account in the present calculation.

### B Field-induced magnetic ordering

As mentioned in the Introduction, the present system can undergo magnetic ordering at a low temperature, when the magnetic field $`H`$ is higher than the gap field $`H_\mathrm{g}`$. Figure 5 shows the low-temperature magnetization measured at various external fields above $`H_\mathrm{g}`$ for $`Hb`$, $`H(1,0,\overline{2})`$ and $`H[2,0,1]`$. With decreasing temperature, the magnetization exhibits a cusplike minimum, irrespective of the field direction, as observed in isostructural TlCuCl3. We assign the temperature with a cusplike minimum in the magnetization to the transition temperature $`T_\mathrm{N}`$. However, there is a certain degree of error in determining the transition points, because the experimental data are somewhat scattered.

When the magnetic field is slightly higher than the gap field $`H_\mathrm{g}`$, the number of created triplets (magnons) is small. Through the transverse components of the interdimer exchange interactions, the magnons can hop to the neighboring dimers in the same way as particles, which have bosonic natures. The longitudinal component of the the interdimer exchange interaction gives rise to the interaction between magnons. Hence, the system can be mapped onto the interacting dilute boson system, and the phase transition can be described by the formation of the coherent state of dilute magnons, i.e., Bose-Einstein condensation (BEC) . The cusplike minimum of the temperature variation of the magnetization is characteristic of the BEC of the dilute magnons , and cannot be described by the mean-field approach from the real space . It is known that the mean-field approximation gives a reasonable description of the physical quantities at $`T=0`$, while it does not at finite temperatures. The mean-field approach gives the temperature-independent magnetization below $`T_\mathrm{N}`$. The density of magnons $`n`$ corresponds to the magnetization per Cu2+ ion $`m`$, i.e., $`m=(g/2)\mu _Bn`$. The reason that the magnetization increases again below $`T_\mathrm{N}`$ is because the increase of the number of condensed magnons is greater than the decrease of the number of thermally excited noncondensed magnons, so that the total number of magnons increases below $`T_\mathrm{N}`$. The present experimental result is consistent with the magnon BEC theory of the phase transition.

The increase of the magnetization below $`T_\mathrm{N}`$ becomes smaller with increasing magnetic field. At $`H=30`$ T, the density of magnons $`n`$ is $`n0.2`$. For this large value of $`n`$, the condition of dilute magnons is no longer satisfied. In the dense magnon region, the hopping of magnon is significantly suppressed, and thus, the mean-field approach from the real space may give a better description of the phase transition.

The transition temperature $`T_\mathrm{N}`$ and transition field obtained for $`Hb`$, $`H(1,0,\overline{2})`$ and $`H[2,0,1]`$ are plotted in Fig. 6. Since the phase boundary depends on the $`g`$-factor, we normalize the phase diagram by the $`g`$-factor. The $`g`$-factors used are $`g=2.05`$ for $`Hb`$, $`g=2.26`$ for $`H(1,0,\overline{2})`$ and $`g=2.04`$ for $`H[2,0,1]`$, which were determined by ESR measurements. Figure 7 shows the phase diagram normalized by the $`g`$-factor. The phase boundaries for three different field directions almost coincide when normalized by the $`g`$-factor, although the boundary for $`H(1,0,\overline{2})`$ tends to deviate from the others at around $`(g/2)H30`$ T. This result may reconfirm that the magnetic anisotropy is negligible in KCuCl3, as in TlCuCl3.

The phase boundary near the gap field can be described by the power law of eq. (1), as predicted by the magnon BEC theory . We fit eq. (1) to the data for $`T<3`$ K, for which $`n<0.03`$. The solid line in Fig. 7 is the fitting with $`\varphi =2.3(1)`$ and $`(g/2)H_\mathrm{g}=22.6(1)`$ T. This exponent $`\varphi `$ is close to $`\varphi =2.1(1)`$ observed in TlCuCl3 . The exponent $`\varphi =2.3`$ obtained by the present measurements is somewhat larger than the value $`\varphi =\frac{3}{2}`$ predicted by the magnon BEC theory . We note that the exponent depends on the temperature region used for fitting. If we use the data up to higher temperatures, we have a larger exponent. This is because the ordered phase exists in a closed area in the magnetic field vs temperature diagram, and thus the entire phase boundary cannot be described by the power law with a single exponent. Since the low-temperature data contribute greatly to the determination of the exponent $`\varphi `$ for the phase boundary near the gap field, the measurements at lower temperatures are needed.

The increase of the transition field obeying the power law with increasing temperature cannot be described in terms of the mean-field approach from the real space  or the temperature dependence of the excitation gap . The mean-field result gives the transition field almost independent of temperature at low temperatures. The temperature dependence of the excitation gap was investigated by means of neutron scattering . The gap is almost independent of temperature up to approximately 10 K, and then increases significantly.

## IV Conclusion

We have presented the results of magnetization measurements on the coupled spin-dimer system KCuCl3. The entire magnetization curve including the saturation region was obtained at $`T=1.3`$ K, using pulsed magnetic fields up to 60 T. Applying mean-field approximation on the interdimer interactions, we analyzed the magnetization curve. It was found that the exchange parameters determined from the dispersion relations of the magnetic excitations should be reduced to fit both the gap field and the saturation field. This suggests the importance of the renormalization effect in the magnetic excitations.

The field-induced magnetic ordering was investigated in static magnetic fields up to 30 T. As observed in isostructural TlCuCl3, the magnetization exhibits a cusplike minimum at the transition temperature. The phase boundary is almost independent of field direction, and is described by the power law. These features are compatible with the magnon BEC theory. Therefore, we conclude that the field-induced 3D magnetic ordering near the gap field in isotropic spin gap systems such as KCuCl3 and TlCuCl3 is universally represented by the BEC of dilute magnons.

###### Acknowledgements.

This work was supported by the Toray Science Foundation and a Grant-in-Aid for Scientific Research on Priority Areas (B) from the Ministry of Education, Culture, Sports, Science and Technology of Japan. A. O. was supported by the Research Fellowships of the Japan Society for the Promotion of Science for Young Scientists.

## A Derivation of the magnetization curve

We calculate the magnetization curve at $`T=0`$, assuming that the triplet excitations are created only on the dimer site. The interdimer interactions are treated by the mean-field approximation. In the magnetization slope region, the transverse spin components form the long-range order. Since the lowest magnetic excitations in KCuCl3 and TlCuCl3 occur at the same reciprocal lattice points, i.e., $`=(h,0,l)`$ with integer $`h`$ and odd $`l`$ in the $`a^{}c^{}`$ plane , we can assume that the spin structure in the ordered phase in KCuCl3 is equivalent to that observed in TlCuCl3 (see Fig. 3) . We express the spin state of a dimer with the total spin $`S=S_1+S_2`$ and the $`z`$ component $`S^z=S_1^z+S_2^z`$ as $`|S,S^z`$, where subscript numbers 1 and 2 distinguish spins in a dimer. In the ordered phase, the transverse components of two spins on a dimer should be antiparallel and their magnitudes should be the same, i.e., $`S_1^{x,y}=S_2^{x,y}`$. This indicates that the spin state $`|1,0`$ does not contribute to the basis state. Thus, we write the basis state of the $`j`$-th dimer as

$$\psi _j=|0,0\mathrm{cos}\theta +\left(|1,1\mathrm{cos}\phi \mathrm{e}^{\mathrm{i}\varphi _j}|1,1\mathrm{sin}\phi \mathrm{e}^{\mathrm{i}\varphi _j}\right)\mathrm{sin}\theta .$$
(A1)
Angles $`\theta `$ and $`\phi `$ were introduced to satisfy the normalization condition. The phase $`\varphi _j`$ corresponds to the angle between the $`x`$-direction and the transverse component of the spin on the $`j`$-th dimer.

The average values of spin operators are given by

$`S_{j,1}^z`$ $`=`$ $`S_{j,2}^z={\displaystyle \frac{1}{2}}\mathrm{sin}^2\theta \mathrm{cos}2\phi ,`$ (A2)
$`S_{j,1}^+`$ $`=`$ $`S_{j,2}^+={\displaystyle \frac{1}{\sqrt{2}}}\mathrm{cos}\theta \mathrm{sin}\theta \left(\mathrm{cos}\phi +\mathrm{sin}\phi \right)\mathrm{e}^{\mathrm{i}\varphi _j},`$ (A3)
$`S_{j,1}^{}`$ $`=`$ $`S_{j,2}^{}={\displaystyle \frac{1}{\sqrt{2}}}\mathrm{cos}\theta \mathrm{sin}\theta \left(\mathrm{cos}\phi +\mathrm{sin}\phi \right)\mathrm{e}^{\mathrm{i}\varphi _j}.`$ (A4)
From the spin structure shown in Fig. 3, we see that all of the phases $`\varphi _j`$ of dimers are the same in a chemical double chain, and that the phases of dimers in the double chains located at the corner and the center of the unit cell in the $`bc`$ plane differ by $`\pi `$. With eq. (A2) and the exchange interactions shown in Fig. 3, the energy per dimer is expressed as

$`E`$ $`=`$ $`{\displaystyle \frac{3}{4}}J\mathrm{cos}^2\theta +{\displaystyle \frac{1}{4}}J\mathrm{sin}^2\theta g\mu _\mathrm{B}H\mathrm{sin}^2\theta \mathrm{cos}2\phi `$ (A5)
$`+`$ $`\stackrel{~}{J}\mathrm{sin}^2\theta \mathrm{cos}^2\theta \left(1+\mathrm{sin}2\phi \right)+{\displaystyle \frac{1}{2}}\overline{J}\mathrm{sin}^4\theta \mathrm{cos}^22\phi ,`$ (A6)
where

$$\stackrel{~}{J}=J_{(100)}^{\mathrm{eff}}2J_{(1,\frac{1}{2},\frac{1}{2})}^{\mathrm{eff}}+J_{(2,0,1)}^{\mathrm{eff}},$$
(A7)
with

$`J_{(100)}^{\mathrm{eff}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(2J_{(100)}J_{(100)}^{}\right),`$ (A8)
$`J_{(1,\frac{1}{2},\frac{1}{2})}^{\mathrm{eff}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(J_{(1,\frac{1}{2},\frac{1}{2})}J_{(1,\frac{1}{2},\frac{1}{2})}^{}\right),`$ (A9)
$`J_{(2,0,1)}^{\mathrm{eff}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}J_{(2,0,1)}^{}`$ (A10)
and

$$\overline{J}=J_{(100)}+\frac{1}{2}J_{(100)}^{}+J_{(1,\frac{1}{2},\frac{1}{2})}+J_{(1,\frac{1}{2},\frac{1}{2})}^{}+\frac{1}{2}J_{(2,0,1)}^{}.$$
(A11)
Angles $`\theta `$ and $`\phi `$ can be determined by $`E/\theta =0`$ and $`E/\phi =0`$, which lead to

$$Jg\mu _\mathrm{B}H\mathrm{cos}2\phi +\stackrel{~}{J}\mathrm{cos}2\theta \left(1+\mathrm{sin}2\phi \right)+\overline{J}\mathrm{sin}^2\theta \mathrm{cos}^22\phi =0,$$
(A12)
and

$$g\mu _\mathrm{B}H\mathrm{sin}2\phi +\stackrel{~}{J}\mathrm{cos}^2\theta \mathrm{cos}2\phi \overline{J}\mathrm{sin}^2\theta \mathrm{cos}2\phi \mathrm{sin}2\phi =0,$$
(A13)
for the ordered state between the gap field $`H_\mathrm{g}`$ and the saturation field $`H_\mathrm{s}`$. The magnetization curve is obtained by solving eqs. (A7) and (A8) self-consistently.

At the saturation field $`H_\mathrm{s}`$, $`\mathrm{sin}\theta =1`$ and $`\mathrm{sin}\phi =0`$. Thus, $`H_\mathrm{s}`$ is given by

$$g\mu _\mathrm{B}H_\mathrm{s}=J\stackrel{~}{J}+\overline{J}.$$
(A14)
At the gap field $`H_\mathrm{g}`$, $`\mathrm{sin}\theta =0`$. Substituting this condition into eqs. (A7) and (A8), we obtain

$$g\mu _\mathrm{B}H_\mathrm{g}=\sqrt{J^2+2J\stackrel{~}{J}}.$$
(A15)
It is noted that the right-hand side of eq. (A10) is equivalent to the lowest excitation energy given by the effective dimer approximation, in which the individual interdimer interactions $`J_{(lmn)}`$ and $`J_{(lmn)}^{}`$ are reduced to an effective interdimer interaction $`J_{(lmn)}^{\mathrm{eff}}`$ as shown in eq. (A5) . If we neglect the highest $`|1,1`$ state, i.e., $`\mathrm{sin}\phi 0`$, we obtain

$$g\mu _\mathrm{B}H_\mathrm{g}=J+\stackrel{~}{J}.$$
(A16)
This result can also be derived from the Tachiki-Yamada theory .


text:
# Shape Alignments of Satellite Galaxies

## 1 Introduction

A number of the satellite galaxies of the Milky Way are distorted, in some cases dramatically (e.g. Dohm-Palmer et al. (2001)), by the tidal interactions with the Galaxy. The obvious manifestations of tidal distortion are streams of material that are liberated by tidal forces to lead or trail the satellite in its orbit. One might also expect the tide to stretch the satellites’ gravitational potential wells along the radial direction, and hence induce a tendency for the bound parts of satellite galaxies to be extended toward their primaries. Since our vantage point precludes a measurement of the satellites’ extent along the Galactic radius vector, it is not possible to test this hypothesis for the Milky Way satellites. Given a large sample of satellites around other galaxies, however, we can test whether the (projected) satellite major axes have a tendency to be circumferential or radial rather than randomly oriented.

While a detection of satellite shape alignments would serve as a crude gauge of the effect of tidal forces on satellite morphologies, it would be a significant nuisance for measurements of galaxy masses using weak gravitational lensing. In such studies, the halo mass distribution of the primary (foreground) galaxy is inferred by measuring the slight tendency toward circumferential alignment that the gravitational lensing distortion of the primary galaxy induces upon the images of distant background galaxies. This “galaxy-galaxy lensing” signal was first proposed and sought by Tyson et al. (1984) and tentatively detected by Brainerd, Blandford, & Smail (1996). More recently it has been detected at high significance by Fischer et al. (2000), Smith et al. (2001)\[Sm01\], and McKay et al. (2001)\[McK01\]. The latter two studies use foreground galaxies with redshifts measured by the Las Campanas (Shectman et al., 1996) and Sloan (York et al., 2000) redshift surveys, respectively, and are able to determine the halo masses as a function of luminosity and morphological type. Such information is extremely valuable—but subject to error if the measured mean circumferential ellipticities are due to intrinsic alignments of the primaries’ satellites rather than by the lensing effect. In galaxy-galaxy lensing measurements to date, there has been no attempt to cull the satellites from the sample of lensed background galaxies.

In this work we test for shape alignment among a sample of satellites of isolated galaxies found in the 2dF Galaxy Redshift Survey (2dFGRS; Colless et al. (2001)), in an effort to see if such alignments are a significant difficulty for galaxy-galaxy lensing studies. The following sections describe the satellite sample; the determination of their mean elongations; and the implications for weak lensing studies and for tidal effects. Throughout the paper, we adopt a flat $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Lambda }_0=0.7`$ cosmology to convert redshift into comoving distance and we assume the value of the Hubble constant to be $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$.

## 2 Satellite Sample

The sample of satellite galaxies around isolated primaries is obtained from the 2dFGRS dataset as of December 2001. The details of how we identify the satellites and primaries within the survey and the tests on the robustness on the method used are given in Norberg et al. (in prep.). However in order to be complete, we briefly repeat the steps of this satellite selection scheme below, which follows the one developed by Zaritsky et al. (1993). The criteria for isolation of the primary galaxy are more relaxed in this work than for studies of satellite galaxy dynamics. Using the full 2dFGRS redshift catalogue containing over 200,000 galaxies, we first select regions of sufficiently high redshift completeness. We then loop over all galaxies which are at least 2.2 magnitudes brighter than the magnitude limit of the 2dFGRS, and localize the galaxies which are isolated. A bright galaxy is isolated and considered as a primary if all its neighbouring galaxies with $`\mathrm{\Delta }V=|V^{\mathrm{prim}}V^{\mathrm{gal}}|\mathrm{\hspace{0.17em}1000}`$ km/s and within a projected distance $`500`$ $`\mathrm{kpc}`$ are such that $`b_\mathrm{J}^{\mathrm{gal}}b_\mathrm{J}^{\mathrm{prim}}+2.2`$

Finally all galaxies around a primary which are within a projected distance $`500`$ $`\mathrm{kpc}`$ and with relative velocity $`\mathrm{\Delta }V=|V^{\mathrm{prim}}V^{\mathrm{gal}}|500`$ km/s are considered as detected satellites around bright isolated galaxies. Only galaxies with a good redshift measurement will be taken up in our catalogue of primaries and satellite galaxies. We also make sure, by using the full photometric 2dFGRS catalogue, that no bright galaxy without redshift measurement lies within the projected radius of an isolated primary. In other words, we eliminate any source of contamination due to spectroscopic determination failure.

We note here that the minimum required magnitude difference between primaries and satellites is about the same as the typical foreground/background magnitude difference in the Sm01 and McK01 lensing studies. This selection makes this satellite sample particularly appropriate for determination of the lensing contamination. However, we have to point out that for most primaries we have only detected one or two satellite galaxies, this mainly being due to the limiting depth of the 2dFGRS, of typically bJ=19.45. As a matter of fact, only 12 primaries in our catalogue of 1227 primaries have 6 or more galaxy satellites around them, and we have also detected 5886 primaries without any spectroscopically confirmed satellite galaxy.

## 3 Ellipticity Determinations

### 3.1 Shape Measurements

Shapes for all satellites are taken from the APM scans of photographic survey plates Maddox et al. (1990). We rotate the APM-determined shapes into a coordinate system with the $`y`$ axis along the radius vector from primary to secondary. There result two ellipticity components which are defined by

$`e_+`$ $``$ $`{\displaystyle \frac{I_{xx}I_{yy}}{I_{xx}+I_{yy}}}`$ (1)
$`e_\times `$ $``$ $`{\displaystyle \frac{2I_{xy}}{I_{xx}+I_{yy}}}`$ (2)
where $`I_{ij}`$ are the quadratic central moments of the galaxy intensity integrated within a bounding isophote. Figure 1 illustrates the meaning of $`e_+`$ and $`e_\times `$. We expect a tidal effect to be manifested as non-zero $`e_+`$. The mean $`e_\times `$ must be zero if the Universe has inversion symmetry, hence we can use $`e_\times `$ as a test for systematic errors.

The unweighted isophotal moments produced by the APM software are noisier and more difficult to correct for point-spread-function (PSF) distortions than more recently devised ellipticity estimators. The satellites are, however, all at least 1 mag above the $`b_J=20.5`$ limit of the APM catalog, and should hence have fairly high signal-to-noise and be well resolved in the APM images. In fact we make no correction at all for PSF effects. Crudely, the measured ellipticity $`𝐞`$ will be related to the true image ellipticity $`𝐞^i`$ and the PSF ellipticity $`𝐞^{}`$ by

$$𝐞=(1R)𝐞^i+R𝐞^{},$$
(3)
where $`R`$ is the ratio of the intensity-weighted moment $`r^2`$ of the PSF to that of the measured image. In our application, the second term has no effect, since $`𝐞^{}=0`$ if the PSF orientation has no correlation with the primary-satellite vector. We can estimate the effect of the circularization factor $`(1R)`$ in the first term by noting that the measured satellite shapes have $`e_+^2=e_\times ^2=(0.31\pm 0.01)^2`$. This measured value will be a factor $`(1R)^2`$ smaller than the intrinsic $`e_+^2`$ of the population. Bernstein & Jarvis (2002)\[BJ02\] find the RMS intrinsic ellipticity of high-surface-brightness galaxies in the nearby Universe to be 0.30. If the intrinsic shapes of satellite are similar to those of high-surface-brightness galaxies in an apparent-magnitude-limited sample, then the circularization factor $`(1R)`$ is unity to within 10% or so. BJ02 find an RMS ellipticity of 0.48 for the lowest-surface-brightness subset of nearby galaxies. If the satellites have this intrinsic dispersion, then $`1R0.6`$ But the higher shape dispersion of low-surface-brightness galaxies is due to a preponderance of disk galaxies, which should be largely absent from the satellite sample. This prejudice is reinforced by the observation that the distribution function for $`|e|`$ values for the satellites is very similar to BJ02’s distribution function for high-surface-brightness galaxies (their Figure 4). It therefore seems more likely that there is negligible PSF circularization in the APM ellipticities, so we will not make any correction for the $`1R`$ term.

As a further check on the possibility of PSF circularization, we split the satellite galaxies into a bright half ($`b_J<18.8`$) and faint half ($`b_J>18.8`$). The brighter half is presumably nearer and better resolved than the fainter half, so a substantial circularization effect would be manifested as a lower RMS ellipticity for the fainter half. The RMS ellipticities $`\sqrt{e_+^2}`$ of the brighter and fainter subsamples are $`0.355\pm 0.008`$ and $`0.298\pm 0.007`$, respectively. The circularization of the fainter half is hence a $`15\%`$ effect; the fainter galaxies could also just be intrinsically rounder.

### 3.2 Mean Ellipticities

An unweighted average of the 1819 satellite ellipticities yields

$$\begin{array}{ccc}\hfill e_+& =& +0.004\pm 0.008\hfill \\ \hfill e_\times & =& +0.005\pm 0.008.\hfill \end{array}$$
(4)
In the weak lensing measurements, one ideally weights low-$`e`$ galaxies more heavily as they offer better sensitivity to lensing shear. Such a weighted average is not justified in considering the effect of tidal forces on intrinsic satellite shapes. It is, however, relevant for assessing the contamination of weak lensing data by satellite alignments since the lensing data are weighted. Following BJ02, we apply a weight $`w(e)=(e^2+0.01)^{1/2}`$ to each satellite. In these scheme we obtain

$$\begin{array}{ccc}\hfill e_+& =& +0.002\pm 0.007\hfill \\ \hfill e_\times & =& +0.007\pm 0.007.\hfill \end{array}$$
(5)
There is clearly no significant detection of net satellite alignment, in either component, and the $`e_\times `$ values show that our uncertainties are sensible.

As illustrated in Figure 2, no alignment is detected even when we restrict the sample to the smaller projected radii where tidal effects should be stronger. We have also tested for alignment in the half of the galaxies that have primaries brighter than $`M_{\mathrm{b}_\mathrm{J}}=21.5`$. The unweighted mean $`e_+=+0.025\pm 0.012`$, a $`2\sigma `$ signal, but the weighted mean is only $`1.4\sigma `$ non-zero and there is no detection of real significance. Further restricting the sample to small projected radii does not increase the measured mean.

## 4 Implications

### 4.1 Tidal Distortions of Satellites

Equations (4) indicate that, for satellites within $`500`$ $`\mathrm{kpc}`$ projected radius, the mean radial or circumferential elongation is less than 2% in projection at 95% confidence. We wish to place some bound on the three-dimensional satellite shapes from the limit on the projected shape. We now define a parameter $`ϵ`$ that characterizes the mean departure of satellite shapes from spherical symmetry. We place a Cartesian coordinate system on each satellite such that the $`y`$ axis is along the primary-satellite radius vector, and the $`z`$ axis is perpendicular to both the line of sight and the $`y`$ axis. The luminosity distribution of each satellite has second central moments $`I_{xx}`$, $`I_{yy}`$, and $`I_{zz}`$ in this system. Because the ensemble of primary-satellite systems must be isotropic with respect to the line of sight, the mean transverse moments must satisfy $`I_{xx}=I_{zz}`$. We can then quantify the mean tidal distortion by defining $`\sigma ^2`$ and $`ϵ`$ such that

$`I_{xx}=I_{zz}`$ $`=`$ $`\sigma ^2(1ϵ/2)`$ (6)
$`I_{yy}`$ $`=`$ $`\sigma ^2(1+ϵ)`$ (7)
A positive value of $`ϵ`$ would indicate the presence of a radial elongation from tidal effects. It is also possible that a nonzero $`ϵ`$ could arise if satellites tend to be found along and align with the dark matter filaments that are believed to underly the large-scale distribution of galaxies. We are not aware of any theoretical estimates of the size (or sign!) expected of such an effect.

If there is no absorption in the galaxy, then the mean projected second moments can be simply related to the intrinsic elongation $`ϵ`$ via

$$e_+=ϵ/2.$$
(8)
We may thus bound the mean radial asymmetry of satellites to be $`ϵ<4\%`$ at 95% CL.

### 4.2 Weak Lensing Measurements

Equations (5) imply that the population of satellites will produce a false lensing distortion of $`|\delta |<0.016`$ at 95% CL. F00 measure the angular correlation between foreground and “background” samples and infer that the pollution of the “background” sample by neighbors of the foreground galaxy—i.e. satellites—is about 15% within impact parameters of 60 $`h^1\mathrm{kpc}`$. The contamination in Sm01 is significantly lower, 10% at most. Hence the deleterious effect of satellite alignment will be significantly diluted at these radii, and will be even weaker at larger radii. For the Sm01 data, we may conclude that the false distortion signal produced by aligned satellites satisfies

$$|\delta _{\mathrm{false}}|<0.0016$$
(9)
at $`95\%`$ confidence within 60 $`h^1\mathrm{kpc}`$ impact parameter. Since Sm01 detect a mean distortion of $`\delta 0.008`$ at this impact parameter around an $`L_{}`$ galaxy at $`z=0.1`$, the error due to satellite alignments is at most a 20% effect. We cannot yet conclude that the effect is totally negligible, but this is likely the case, especially at larger impact parameters, where the satellite fraction is smaller, and the tidal effect drops more rapidly than the lensing shear.

Larger and deeper satellite samples containing more primaries with many satellite galaxies each may reveal a significant alignment, and future galaxy-galaxy lensing measurements will be aiming for higher precision. If the satellite alignments do ultimately pose a barrier to high precision, then their deleterious effect could be eliminated by using photometric redshift estimates to exclude the bulk of the satellite population from the sample of lensed background galaxies. Satellite alignments do not, however, seem to be a problem for current galaxy-galaxy lensing data.

We are grateful to the 2dFGRS team for allowing us to use their data prior to publication and to Carlos Frenk for useful discussions. The 2dFGRS is being carried out using the 2 degree field facility on the 3.9m AngloAustralian Telescope (AAT). We thank all those involved in the smooth running and continued success of the 2dF and the AAT. This work was supported in part by grant AST-9624592 from the National Science Foundation and in part by a PPARC rolling grant at Durham. PN is supported by the Swiss National Science Foundation and an ORS award.


text:
# Ground state of finite nuclei evaluated from realistic interactions

## I Introduction

One of the central aims of theoretical nuclear physics is the attempt to determine the bulk properties of nuclear systems like its binding energy and its density or radius from a realistic model of the nucleon-nucleon (NN) interaction, i.e. a NN interaction which yields an accurate fit of the NN scattering data below the threshold of $`\pi `$-production. Quite some progress has been made during the last ten years in the determination of such realistic NN interaction models. A family of charge-dependent NN potentials has been generated, which all fit the empirical scattering data with high precision. Nevertheless, using such NN interactions in nuclear structure calculations, one obtains results, which exhibit significant differences.

The NN interactions containing non-local terms like the charge-dependent Bonn potential, CDBonn, or one of the Nijmegen interaction model tend to be “softer” than the purely local interactions like the Argonne V18 potential. Here the expression softer interaction is used to identify those interactions which induce weaker two-nucleon correlations in the nuclear many-body wave function. This implies that the total energy of nuclear matter calculated at a given density without correlations, i.e. using the Hartree-Fock approximation, yields less repulsive results for a soft, non-local potential as compared to the energy calculated for a stiffer, local NN interaction. As an example we mention that Hartree-Fock calculations for nuclear matter at the empirical saturation density yield 4.6 MeV per nucleon using the CDBonn interaction while 30.3 MeV are obtained if the Argonne V18 potential is used. Including effects of NN correlation as it is done e.g. in a Brueckner-Hartree-Fock (BHF) calculation leads to results which are attractive and much closer to each other.

All these modern interactions, however, tend to predict too much binding energy for nuclear matter and saturation densities which are too large. This is the case for the BHF approach but is also observed if variational calculations are performed Softer non-local interaction yield larger binding energies and saturation densities than local interactions. This over-binding is often compensated by a three-nucleon force, which can be adjusted in such a way that the empirical saturation point of nuclear matter is reproduced. These three-nucleon forces can be understood to simulate the relativistic effects as incorporated e.g. in the Dirac-BHF approach. These three-nucleon forces, however, can also be interpreted to describe the effects of virtual excitations of nucleons to the $`\mathrm{\Delta }`$(3,3) or $`N^{}`$ Roper resonance.

Phenomenological three-nucleon forces are also employed to describe the properties of light nuclei with particle number up to A=8 using the Green function Monte-Carlo method. However, at first sight, the situation in calculating bulk properties of finite nuclei seems to be different from the corresponding situation in nuclear matter. While microscopic calculations for nuclear matter, using the modern NN interactions, yield too much binding energy, BHF calculations but also calculations using the coupled cluster approach predict binding energies for finite nuclei, which are too small as compared to the experimental data.

It is the aim of the present work to investigate this situation a bit more in detail. For that purpose we consider the BHF approximation paying special attention to the particle-particle excitations at low energies and including the effects of hole-hole scattering terms. We are investigating the predictions originating from various models of realistic NN interactions. These include local, like the Argonne V18, and non-local models, like the CDBonn interaction, but also the realistic meson-exchange models based on chiral perturbation theory, which have recently been developed by the Idaho group.

After this introduction we will describe in section 2 a technique for the calculation of bulk properties of nuclei which account in a consistent way for short-range and long-range correlations. Results of such calculation for the nuclei 16O and 40Ca will be presented and discussed in section 3, which is followed by a short summary in section 4.

## II Consistent treatment of long- and short-range correlations

The BHF approximation is one of the most popular approaches to account for effects of correlations beyond the mean-field approximation in calculating bulk properties of infinite nuclear matter and finite nuclei from realistic NN interactions. It is characterized by a solution of the Bethe-Goldstone equation, leading to the $`G`$-matrix $`𝒢_\mathcal{0}`$

$$𝒢_\mathcal{0}(\omega )=V+V\frac{Q_0}{\omega Q_0H_0Q_0}𝒢_\mathcal{0},$$
(1)
and the self-consistent evaluation of the BHF single-particle energies

$$\epsilon _i=i\left|\frac{p^2}{2m}\right|i+\underset{j<F}{}ij\left|𝒢_\mathcal{0}(\omega =\epsilon _i+\epsilon _j)\right|ij.$$
(2)
In these equations $`V`$ stands for the bare NN interaction, $`Q_0`$ represents the Pauli operator, which restricts the propagator in the Bethe-Goldstone Eq.(1) to particle states with energies above the Fermi energy, and $`i`$ and $`j`$ refer to hole-states, i.e. the eigenstates of the BHF single-particle hamiltonian with energies $`\epsilon _i`$ and $`\epsilon _j`$ below the Fermi energy. The self-consistent definition of the starting energy $`\omega `$ in terms of the single-particle energies is determined by the Bethe-Brandow-Petschek (BBP) theorem. Since, however, the BBP theorem can only be used to justify the definition of the single-particle potential for the hole states, the optimal definition of the single-particle energies for the particle states, which enter the Bethe-Goldstone Eq.(1) because they define the eigenvalues of $`H_0`$, has been discussed over many years in the literature. The conventional choice has been to ignore a single-particle potential or self-energy contributions for the particle states completely and approximate $`H_0`$ by the kinetic energy only. This conventional choice is supported by the coupled cluster or exponential S method, which using the $`S_2`$ approximation essentially leads to the same approach, This conventional choice for the single-particle spectrum of nuclear matter is not very appealing as it leads to a large gap at the Fermi surface.

Mahaux and collaborators argued that it would be more natural to choose the propagator according to the Green function method, i.e. define the single-particle propagator with a single-particle energy which includes the real part of the self-energy as a single-particle potential for particle and hole states. This leads to a spectrum which is continuous at the Fermi momentum, which provided the name “continuous choice” for this approach. This continuous choice leads to an enhancement of correlation effects in the medium and tends to predict larger binding energies for nuclear matter than the conventional choice. Inclusion of three-hole-line contribution indicate that the continuous choice seems to lead to a better convergence of the hole-line expansion and is therefore preferable. In fact, recent studies in nuclear matter show that the result is rather sensitive to details of the single-particle spectrum around the Fermi energy.

In order to investigate these effects of the low-energy particle-particle excitations, which should correspond to long-range correlations, also in finite nuclei, we follow the concept of a double partitioned Hilbert space, as it has been used before for the study of infinite nuclear matter (see e.g. ) as well as finite nuclei (see e.g. ). The long range correlations are taken into account by means of the Green function approach within a finite model space. This model space shall be defined in terms of shell-model configurations including oscillator single-particle states up certain shell. For our studies of 16O and 40Ca we have chosen to include configurations up to the pf shell and sdg shell, respectively. It turned out that the final results are not very sensitive to this choice. This is due to the fact that the effect of short-range correlations, i.e. originating from configurations outside this model space are not ignored but taken into account by means of an effective interaction, i.e. a $`G`$-matrix appropriate for this model space. This effective interaction $`𝒢`$ is determined as the solution of the Bethe-Goldstone equation

$$𝒢(\omega )=V+V\frac{Q_{\text{mod}}}{\omega Q_{\text{mod}}H_0Q_{\text{mod}}}𝒢.$$
(3)
The Pauli operator $`Q_{\text{mod}}`$ is defined to exclude two-particle states with one of the particles below the Fermi level of the nucleus considered or with both nucleons in states inside the model space. With this $`G`$-matrix we solve the BHF Eq.(2) using the self-consistent definition of the starting energy $`\omega `$ for hole states but also for the oscillator particle states, which are inside the model space. We consider two different choices for the spectrum of the high-lying particle states outside the model space, defined by $`H_0`$ in Eq.(3): The conventional choice, i.e. $`H_0`$ is just the operator of the kinetic energy for the interacting particles, and the continuous choice, for which we add an attractive constant to the kinetic energy. This constant is determined in such a way, that it corresponds to the mean value of the potential energies of the low-lying particle states inside the model space.

Note that this approximate BHF scheme only accounts for particle-particle correlations outside the model space. We will call these correlations short-range correlations in the discussion below. A measure for these short-range correlations is given by the depletion coefficient

$$\kappa _i=\underset{j<F}{}ij\left|\frac{𝒢}{\omega }\right|ij$$
(4)
or the corresponding occupation probability

$$\rho _i^{\text{sr}}=1\kappa _i.$$
(5)
Using the single-particle energies we define an effective interaction for the model space in terms of the oscillator matrix elements

$$ij\left|𝒱\right|kl=\frac{1}{2}\left[ij\left|𝒢(\omega =\epsilon _i+\epsilon _j)\right|kl+ij\left|𝒢(\omega =\epsilon _k+\epsilon _l)\right|kl\right].$$
(6)
The effects of long-range correlations or correlations inside the model space shall be evaluated by means of the Green function method. To determine the correlated single-particle Green function for a nucleon with isospin $`\tau `$, orbital angular momentum $`l`$ and total angular momentum $`j`$ one has to solve a Dyson equation, which in the case of a finite system with spherical symmetry and within a discretize model space takes the form

$$g_{\tau lj}(n,n^{};\omega )=g_{\tau lj}^{(BHF)}(n,n^{};\omega )+\underset{n^{\prime \prime },n^{\prime \prime \prime }}{}g_{\tau lj}^{(BHF)}(n,n^{\prime \prime };\omega )\mathrm{\Delta }\mathrm{\Sigma }_{\tau lj}(n^{\prime \prime },n^{\prime \prime \prime };\omega )g_{\tau lj}(n^{\prime \prime \prime },n^{};\omega ),$$
(7)
where $`g_{\tau lj}^{(BHF)}(n,n^{\prime \prime };\omega )`$ refers to the BHF propagator, which we assume to be diagonal in the radial quantum numbers $`n`$, $`n^{}`$

$$g_{\tau lj}^{(BHF)}(n,n^{};\omega )=\delta _{n,n^{}}\frac{1}{\omega \epsilon _{n\tau lj}\pm i\eta }.$$
(8)
The correction to the BHF self-energy in terms of two particle one hole ($`2p1h`$) and two hole one particle ($`2h1p`$) configurations inside the model space

$$\mathrm{\Delta }\mathrm{\Sigma }_{\tau lj}(n,n^{},\omega )=\mathrm{\Sigma }_{\tau lj}^{(2p1h)}(n,n^{},\omega )+\mathrm{\Sigma }_{lj}^{(2h1p)}(n,n^{},\omega ),$$
(9)
is defined in terms of the effective interaction $`𝒱`$ of Eq.(6). As an example we consider the $`2p1h`$ contribution of second order in $`𝒱`$

$$\mathrm{\Sigma }_{\tau lj}^{(2p1h)}(n,n^{},\omega )=\frac{1}{2}\underset{h<F}{}\underset{p_1,p_2>F}{}\frac{<nh|𝒱|p_1p_2><p_1p_2|𝒱|n^{}h>}{\omega (\epsilon _{p1}+\epsilon _{p2}\epsilon _h)+i\eta }.$$
(10)
In order to evaluate the total energy of the system and expectation values of one-body operators one has to rewrite the single-particle Green function in the Lehmann representation

$$g_{\tau lj}(n,n^{};\omega )=\underset{\alpha }{}\frac{S_{\tau lj}^p(n,n^{},\omega _{\alpha \tau lj})}{\omega \omega _{\alpha \tau lj}+i\eta }+\underset{\beta }{}\frac{S_{\tau lj}^h(n,n^{},\omega _{\beta \tau lj})}{\omega \omega _{\beta \tau lj}i\eta }.$$
(11)
The single-particle density matrix is then defined in terms of the hole-spectral function

$$\rho _{\tau lj}(n,n^{})=\underset{\beta }{}S_{\tau lj}^h(n,n^{},\omega _{\beta \tau lj})$$
(12)
and the total energy can be evaluated from

$$E=\underset{\tau lj\beta n,n^{}}{}\frac{(2j+1)}{2}S_{\tau lj}^h(n,n^{},\omega _{\beta \tau lj})n\left|\frac{p^2}{2m}+\omega _{\beta \tau lj}\right|n^{}.$$
(13)
In order to obtain the spectral function and the positions of the poles in the single-particle Green function, $`\omega _{\beta \tau lj}`$, we reformulate the Dyson equation into an eigenvalue problem

$$\left(\begin{array}{ccccccccc}\epsilon _1& & 0& a_{11}& \mathrm{}& a_{1P}& A_{11}& \mathrm{}& A_{1Q}\\ & \mathrm{}& & \mathrm{}& & \mathrm{}& \mathrm{}& & \mathrm{}\\ 0& & \epsilon _N& a_{N1}& \mathrm{}& a_{NP}& A_{N1}& \mathrm{}& A_{NQ}\\ a_{11}& \mathrm{}& a_{N1}& e_1& & & 0& & \\ \mathrm{}& & \mathrm{}& & \mathrm{}& & & & \\ a_{1P}& \mathrm{}& a_{NP}& 0& & e_P& & & 0\\ A_{11}& \mathrm{}& A_{N1}& & & & E_1& & \\ \mathrm{}& & \mathrm{}& & & & & \mathrm{}& \\ A_{1Q}& \mathrm{}& A_{NQ}& 0& \mathrm{}& 0& & \mathrm{}& E_Q\end{array}\right)\left(\begin{array}{c}X_{0,1}^\alpha \\ \mathrm{}\\ X_{0,N}^\alpha \\ X_1^\alpha \\ \mathrm{}\\ X_P^\alpha \\ Y_1^\alpha \\ \mathrm{}\\ Y_Q^\alpha \end{array}\right)=\omega _\alpha \left(\begin{array}{c}X_{0,1}^\alpha \\ \mathrm{}\\ X_{0,N}^\alpha \\ X_1^\alpha \\ \mathrm{}\\ X_P^\alpha \\ Y_1^\alpha \\ \mathrm{}\\ Y_Q^\alpha \end{array}\right),$$
(14)
where for simplicity we have dropped the corresponding conserved quantum numbers for isospin, parity and angular momentum ($`\tau lj`$). The matrix to be diagonalized contains the BHF hamiltonian defined in terms of the $`N`$ BHF single-particle energies of the symmetry assumed within the model space, the coupling to the $`P`$ different $`2p1h`$ configurations and $`Q`$ $`2h1p`$ states. These couplings are expressed in terms of the matrix elements

$$a_{mi}=<mh|𝒱|p_1p_2>\text{and}A_{mj}=<mp|𝒱|h_1h_2>.$$
(15)
As long as we are still ignoring any residual interaction between the various $`2p1h`$ and $`2h1p`$ configurations the corresponding parts of the matrix in (14) are diagonal with elements in terms of single-particle energies and denoted by $`e_i`$ ($`E_j`$) for $`2p1h`$ ($`2h1p`$) configurations, respectively. One can easily improve this approach and incorporate the effects of residual interactions between the $`2p1h`$ configurations or $`2h1p`$ configurations. One simply has to modify the corresponding parts of the matrix in Eq.(14) and replace where $`_{2p1h}`$ and $`_{2h1p}`$ contain the residual interactions in the $`2p1h`$ and $`2h1p`$ subspaces

$$\left(\begin{array}{ccc}e_1& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \\ 0& \mathrm{}& e_P\end{array}\right)_{2p1h},\text{and}\left(\begin{array}{ccc}E_1& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \\ 0& \mathrm{}& E_Q\end{array}\right)_{2h1p}.$$
(16)
For the calculations which we are going to discuss here, we only consider the matrix elements for the particle-particle interaction in $`_{2p1h}`$ and the hole-hole interaction for $`_{2h1p}`$. This implies that the corresponding particle-particle and hole-hole ladder diagrams are taken into account. A more complete treatment of the residual interaction requires the treatment of three-body terms and has recently been discussed e.g. by Barbieri et al.

The eigenvalues $`\omega _\alpha `$ of Eq.(14) correspond to the poles of the single-particle Green function in Eq.(11) and the spectral function is given in terms of the components of the eigenvectors by

$$S^h(n,n^{};\omega _\alpha )=X_{0,n}^\alpha X_{0,n^{}}^\alpha $$
(17)
for eigenvalues $`\omega _\alpha `$ below the Fermi energy $`E_F`$, while a corresponding equation holds for the spectral function $`S^p`$ for energies above $`E_F`$.

## III Results and discussion

As a first example we would like to consider some results displayed in table I evaluated for the nucleus 16O using the CDBonn potential supplemented by the Coulomb interaction between protons. For the results displayed in this table the single-particle wave functions have been constrained to wave functions of the harmonic oscillator defined in terms of an harmonic oscillator constant

$$b=\frac{1}{\sqrt{2}\alpha },$$
(18)
with $`\alpha `$ = 0.4 fm-1. This corresponds to an oscillator frequency of

$$\mathrm{}\omega =\frac{(\mathrm{}c)^2}{mc^2b^2}$$
(19)
of 13.27 MeV and leads to a radius for 16O, assuming simple shell model occupancies, of 2.65 fm, which is close to the empirical value. The first two columns of table I refer to BHF calculation assuming the conventional and a continuous choice for the spectrum of the particle states in the Bethe-Goldstone equation (1). It should be recalled that we define the continuous spectrum in terms of the kinetic energy shifted by a constant such that the single-particle energies for low-lying particle states is identical to the mean value of the corresponding BHF single-particle energies calculated according to (2). Beside the single-particle energies also the occupation probabilities $`\rho _i^{\text{sr}}`$ calculated according to Eq. (5) are listed. Comparing these two columns one can see that the use of a continuous single-particle spectrum leads to more attractive single-particle energies and a binding energy, which is enhanced by about 2 MeV per nucleon. This enhancement of the attraction is accompanied by lower occupation probabilities $`\rho _i^{\text{sr}}`$, which implies larger values for the corresponding depletion coefficients. This demonstrates that the lowering of the particle-state spectrum, going from the conventional to the continuous choice, leads to a substantial enhancement of correlations.

The third and fourth column of table I, labeled BHF0 and BHF(mod), refer to the model space approach as introduced in the preceding section. The BHF0 approach identifies the BHF calculation using the G-matrix defined in (3) with a continuous choice for the particle state spectrum. This means that due to the Pauli operator $`Q_{\text{mod}}`$ the contribution of low-lying particle states (for the present example those in the 1s0d and 1p0f shell) to the G-matrix are suppressed. This leads to larger occupation probabilities (compare BHF0 with BHF (cont)) and less attractive single-particle energies. Also the binding energy per nucleon is reduced by about 1 MeV.

The contribution of the low-lying particle-particle configurations are again included in the BHF(mod) approach by considering the corresponding (2p1h) component in the definition of $`\mathrm{\Delta }\mathrm{\Sigma }`$ in Eq.(9) and solving Eq.(14) ignoring the coupling to the (2h1p) configurations but allowing for the residual interaction between the particle states according to Eq.(16). In this way the effects of the low-lying particle-particle configurations are taken into account allowing for individual single-particle energies $`\epsilon _{nlj}`$ for all sub-shells with quantum numbers $`nlj`$ and not just a replacement of these single-particle energies with the kinetic energy shifted by a global constant as it was done in the BHF (cont) approach discussed above. This more sophisticated treatment of the low-lying particle spectrum leads to some additional attraction in the single-particle energies ranging from 0.9 MeV in the case of the $`s_{1/2}`$ shell over 1.6 MeV ($`p_{3/2}`$) to 1.8 MeV in the case of the $`p_{1/2}`$ shell. The effect is obviously larger for states closer to the Fermi energy as these states are more sensitive to the details of the long-range correlations. The same feature can also be observed in the occupation probabilities $`\rho _{lj}`$. It should be pointed out that the results listed for the BHF (mod) approach accounts for depletion due to the excitation of particle state configurations inside the model space by means of Eq.(12) while the depletion due to short-range correlations leading to excitations outside the model space are accounted for by means of Eq.(5). The more specific treatment of the low-lying particle-particle configuration reduces the spin-orbit splitting for the protons in the $`p`$-shell from 3.4 MeV in the case of BHF(cont) to 3.1 MeV in the BHF(mod) approach, both values being much smaller than the experimental one (6.3 MeV).

For all approaches discussed so far, the single-particle strength for the hole states is just concentrated in one quasiparticle state. The corresponding Green functions exhibit only one pole below the Fermi energy. A distribution of this hole-strength is obtained if we also account for the 2h1p contribution in the definition of the self-energy in Eq.(9). This distribution is defined in terms of the eigenvalues $`\omega _\alpha `$ in Eq.(14) and the strength defined in Eq.(17). Results for the spectral distribution are displayed in Fig. 1. In the upper panel of this figure, referring to the removal of protons with $`s_{1/2}`$ quantum numbers, the strength is widely distributed. In the column labeled “Green” of table I we give the mean value of this spectral distribution defined by

$$\epsilon _{\tau lj}=\frac{_{\beta n}\omega _{\beta \tau lj}S_{\tau lj}^h(n,n,\omega _{\beta \tau lj})}{_n\rho _{\tau lj}(n,n)}$$
(20)
but also the energy of the quasiparticle state (in brackets), which is defined by that eigenvalue $`\omega _{\beta \tau lj}`$, which carries the largest strength. Inspecting Fig. 1 and the corresponding numbers in table I, one finds that the quasiparticle state for the $`s_{1/2}`$ state carries only rather little strength.

The quasiparticle contribution is much more important for the $`p_{3/2}`$ and $`p_{1/2}`$ states. For these more weakly bound states the quasiparticle states are also those which are closest to the Fermi energy, which means that they correspond to the removal of a nucleon leading to the ground state or lowest excited state of the daughter nucleus of given parity ($`l`$) and angular momentum ($`j`$). Therefore one should consider those quasiparticle energies in comparing with experimental removal energies for these states (values presented in the last column of table I). The spin-orbit splitting resulting from these quasiparticle energies is slightly larger than the one derived from BHF(mod) approximation but still too small as compared with the experimental value. Here it should be recalled that a substantial enhancement of the spin-orbit splitting is obtained if the relativistic features of the Dirac Brueckner Hartree Fock approximation are taken into account.

Comparing the results of the Green function approach with those obtained in the BHF(mod) approximation, one finds that the hole-hole terms which are included in the Green function approach tend to reduce the binding energy of the quasiparticle state but lead to more attractive mean values of the spectral distribution for the hole states. Since this mean value enters the evaluation of the total binding energy (see Eq.(13)) we also get a slightly larger binding energy in the Green function approach as compared to the BHF(mod) approximation. It is worth mentioning that the Green function approach with inclusion of 2h1p terms in the self-energy also provides non-vanishing occupation probabilities for states which are unoccupied in the simple shell model. Also these occupancies contribute to the total binding energy.

At this stage it is useful to make a first comparison with the situation in corresponding calculations of nuclear matter. It has been observed also for nuclear matter that the results of BHF calculations are rather sensitive to the choice of the spectrum for the particle states, in particular for those with momenta close to the Fermi momentum. Using the precise single-particle spectrum rather than a quadratic parametrisation leads to an increase of the binding energy in nuclear matter around saturation density of about 1.5 MeV per nucleon. This is even more than the gain in binding energy from the BHF(cont) to the BHF(mod) approach displayed in table I. Also in nuclear matter one observes that the inclusion of the 2h1p terms in the self-energy lead to less attractive quasiparticle energies. The spread of the single-particle strength to lower energies, however, leads to additional binding energy. This feature, less attractive quasiparticle energies but more binding energy per nucleon, helps to fulfill the Hugenholtz - Van Hove theorem for nuclear matter. The gain in binding energy due to the 2p1h components in the self-energy is about 0.5 MeV per nucleon in nuclear matter at saturation density. Also this is a result rather similar to the difference between the BHF(mod) approach and the result of the complete Green function listed in table I.

For a first comparison of results originating from different NN interactions, we list in table II results of the Green function calculation for 16O using the charge-dependent Bonn interaction (CD Bonn), the version A of the Idaho interaction (Idaho A), which is based on features of chiral perturbation theory, and the Argonne V18 interaction. All interaction models yield quasiparticle energies which are more attractive than the experimental removal energies but a total binding energy which is slightly smaller than the empirical energy of 7.98 MeV per nucleon.

The prediction for the spectral strength of the $`p_{3/2}`$ and $`p_{1/2}`$ quasiparticle states range between 0.71 and 0.75. These values are larger than spectroscopic factors deduced from nucleon knock-out experiments, $`(e,e^{})p`$, which are around 0.6. The local NN interaction Argonne V18 is stiffer than the non-local meson-exchange interaction models CD Bonn and Idaho A (see also the comparison in nuclear matter ). This means that it predicts stronger correlations, as indicated by the smaller spectroscopic factors, and smaller binding energy.

Altogether one may argue that the bulk properties of 16O are very well reproduced. All interactions reproduce the total energy within 1 MeV per nucleon and predict a radius for the proton distribution of 2.72 fm, which is in very good agreement with the experiment. Here, however, one must keep in mind that the radius is to a large extent determined by the choice of the oscillator parameter for the model space. The choice of $`\alpha =0.4`$ has been made to obtain a radius of 2.65 fm for the uncorrelated shell-model wave function. Correlations within the model space lead to small enhancements, only. In order to test the sensitivity of the calculation on the oscillator parameter, we have performed calculations for various values of $`\alpha `$. Results for the binding energy of 16O obtained in the Green function approach are displayed in Fig. 2 as a function of the oscillator parameter $`\alpha `$. These figures show that an energy minimum is obtained for $`\alpha `$ around 0.475 fm-1 (0.45 fm-1, 0.44 fm-1) using the CD Bonn, the Idaho A and the Argonne V18 interaction, respectively. The corresponding radii of the proton distribution are 2.29 fm, 2.41 fm and 2.48 fm, all of them to small as compared to the experimental value. Therefore the situation for calculating bulk properties of 16O is similar to the attempts of evaluating the saturation point of nuclear matter. The calculation tend to predict a radius which is to small or a density which is too large. It is worth noting, however, that the disagreement for nuclear matter is larger.

As a second example we also consider the nucleus 40Ca. In this case we consider a model space, which is defined in terms of harmonic oscillator wave functions including all shells up to the 2s1d0g shell. For a first comparison we fix the oscillator constant to $`\alpha `$ = 0.35 fm-1. This corresponds to simple shell model prediction for the radius of 3.5 fm, in reasonable agreement with the empirical value of 3.48 fm.

Examples for the spectral distribution are displayed in Fig. 3 (assuming the CD Bonn interaction) and numerical results obtained in the Green function approach are listed in table III. Since we have two quasiparticle states for the $`s_{1/2}`$ channel in the BHF approximation (see upper panel of Fig. 3), we also list two energies and occupation probabilities for this partial wave in table III. Also for this nucleus one obtains a broad distribution of the spectral strength. A well defined quasiparticle peak only shows up for the states in the 1s0d shell.

Also in this case, the predictions for the removal energy (absolute value of the quasiparticle energy for $`d_{3/2}`$, are larger than the experimental value of 8.3 MeV. Note, however, that the inclusion of 2h1p contributions in the Green function approach reduces the discrepancy by more than 3 MeV as compared to the BHF(mod) approximation. On the other hand, the calculated binding energies are in good agreement with the experimental value of 8.55 MeV per nucleon. So again, the calculation of bulk properties of 40Ca is rather successful, if one fixes the radius with the appropriate choice of the oscillator constant $`\alpha `$, which defines the basis of the model space.

If one releases this constraint and considers various values for $`\alpha `$ one obtains results as displayed in Figs. 4 and 5. Also in this case we observe that the minima occur for oscillator parameters which are larger than $`\alpha `$ = 0.35 fm-1, which means that the corresponding radii are smaller than the empirical one. This is true for the various approximations (see Fig. 4) but also for the various interactions (see Fig. 5).

The comparison of the various approximation schemes in Fig. 4 confirms the conclusions, which have been given above for the case of 16O. In the case of 40Ca, however, all interactions predict too much binding energy and a radius, which in the case of CD Bonn interaction is about 25 percent smaller than the experimental value. This corresponds to an average density, which is too large by about a factor of two, a situation which is approaching the situation of nuclear matter. Also in this case, the stiffer interaction Argonne V18 yields less energy and smaller densities than the softer interactions (Idaho A and CD Bonn).

## IV Conclusions

Starting from the Brueckner-Hartree-Fock (BHF) approach various approximation schemes have been investigated to derive bulk properties of finite nuclei from realistic NN interactions. It is observed that the results of BHF calculations, for finite nuclei as well as infinite nuclear matter, are rather sensitive to the spectrum of particle states, in particular those with energies close to the Fermi energy. Therefore a technique has been applied, to separate a detailed description of the low-energy excitations corresponding to long range correlations from the treatment of short-range correlations. While the effects of short-range correlations are taken into account by means of the G-matrix approximation of the BHF scheme, the long-range correlations are considered within the framework of the self-consistent evaluation of single-particle Green function. This approach includes the effects of particle-particle correlations but also corresponding hole-hole scattering terms.

If the basis of the model space is constrained to obtain the empirical value for the radius, one is able to reproduce the binding energy of nuclei (we consider 16O and 40Ca) within 0.5 MeV per nucleon, which is very satisfactory for a calculation based on realistic NN interaction. This success is due to the attraction which is obtained from the careful treatment of the low-lying particle-particle excitations but also due to the inclusion of two-hole one-particle configurations in the definition of the self-energy for the single-particle propagator. These terms yield a distribution of the single-particle strength and a shift of the quasiparticle energy. The repulsive shift of the quasiparticle energy for states close to the Fermi level improves the agreement with the experimental removal energies. The distribution of the strength, on the other hand, allows a gain in binding energy although the quasiparticle energies are less attractive. This helps to improve the fulfillment of the Hugenholtz - Van Hove theorem. The calculated spectroscopic factors for the quasiparticle states (around 0.7) are slightly above the empirical values (around 0.6) derived from $`(e,e^{})p`$ experiments.

If the constraining condition for the radius is released, the minima for the calculated energies lead to radii which are too small and densities which are too large. The calculated binding energies tend to be larger than the empirical values. This is more pronounced for the nucleus 40Ca than for 16O. This deficiency could be cured by considering the repulsive three nucleon forces, which are adjusted to describe the saturation point of infinite nuclear matter also for the calculation of finite nuclei.