text:
# Metal-insulator crossover in superconducting cuprates in strong magnetic fields

\[

## Abstract

The metal-insulator crossover of the in-plane resistivity upon temperature decrease, recently observed in several classes of cuprate superconductors, when a strong magnetic field suppresses the superconductivity, is explained using the $`U(1)\times SU(2)`$ Chern-Simons gauge field theory. The origin of this crossover is the same as that for a similar phenomenon observed in heavily underdoped cuprates without magnetic field. It is due to the interplay between the diffusive motion of the charge carriers and the “peculiar” localization effect due to short-range antiferromagnetic order. We also calculate the in-plane transverse magnetoresistance which is in a fairly good agreement with available experimental data.

\]

The in-plane resistivity in heavily underdoped samples of cuprates (in particular LSCO) exhibits a minimum and a crossover from metallic to insulating behavior upon the temperature decrease. Recently a similar crossover was observed in several classes of superconducting cuprates when a strong magnetic field (up to 60 Tesla) suppresses the superconductivity. The “obvious” interpretation in terms of two-dimensional (2D) localization or 2D insulator-superconductor transition is ruled out, as the sheet resistance, defined as $`\rho _{sh}\rho _{ab}/a`$, where $`a`$ is the interlayer distance, at the crossover point is between 1/25 to 1/12 in units of $`h/e^2`$, or, using a free electron model, the estimated product $`k_Fl`$, where $`k_F`$ is the Fermi momentum, $`l`$ the mean free path, is between 12 and 25! The insulating ground state persists up to optimal doping in LSCO and in electron-doped superconductors Pr2-xCexCuO4 (PCCO), while in newly studied Bi2Sr2-xLaxCuO6+δ (BSLCO, or La-doped Bi-2201) it persists only up to 1/8 hole-doping without showing any signature of stripe formation. Thus ascribing this metal-insulator (MI) crossover to a quantum critical point related to charge density instability is open to objections. There were several attempts to interpret the insulating ground state in doped cuprates using various arguments on non-Fermi liquid behavior. However, the crossover phenomenon as temperature varies has not been addressed up to now, to the best of our knowledge. It was thought earlier that the MI crossover in the absence of magnetic field and that induced by magnetic field is of different origin (the sheet resistance in the first case was substantially higher). This is doubtful, because the recent experiments on YBCO show a MI crossover in nonsuperconducting compounds in the absence of magnetic field , and the same type of MI crossover (with comparable sheet resistance at the crossover) in the superconducting compounds of the same series doped with Zn in the presence of magnetic field. We believe the two crossover phenomena are of the same origin. We have used the $`U(1)\times SU(2)`$ Chern–Simons (C.S.) approach to the $`tJ`$ model, proposed by us earlier  to explain the MI crossover in heavily underdoped cuprates in the absence of magnetic field. In this Letter we generalize our formalism to include the effect of magnetic field and show that such a MI crossover is a universal feature of doped cuprates, and it is due to a “peculiar charge localization” effect (using the wording of Ref. ), resulting from the interplay of the spin-excitation gap (corresponding to short-range antiferromagnetic order (SRAFO)) and the holon induced anomalous dissipation. Moreover, we will show that the observed large positive in-plane transverse magnetoresistance (MR) at low temperatures  can be semi-quantitatively explained within this formalism.

The $`U(1)\times SU(2)`$ C.S. gauge field approach is a particular scheme of slave-particle formalism to treat the $`tJ`$ model based on an exact identity , introducing a $`U(1)`$ field gauging the global charge symmetry and a $`SU(2)`$ field gauging the global spin symmetry, both with C.S. actions. Using an optimization procedure of free energy, a careful mean field (MF) approximation gives the following results: The $`U(1)`$ gauge field for low doping $`\delta `$ develops a $`\pi `$-flux per plaquette converting holons into Dirac fermions with a Fermi energy $`ϵ_Ft\delta `$. The holons induce a vortex structure in the MF configurations of the $`SU(2)`$ gauge field, with vortices centered at the holon positions. These dressed holons in turn are seen as slowly moving impurities by spin waves giving rise to a spinon mass $`m_s\sqrt{\delta |\mathrm{ln}\delta |}`$. Notice that this feedback is self-consistent because for low $`\delta `$ we have $`ϵ_Ft\delta <<ϵ_sJ\sqrt{\delta |\mathrm{ln}\delta |}`$. We use $`J=0.1eV,t/J=3`$ in our numerical computations. Due to a special choice of “gauge fixing” (using the Néel gauge) the $`SU(2)`$ gauge field becomes physical, describing the spin fluctuations. The spinon action is given by a nonlinear $`\sigma `$-model with a mass term (spinon gap) which in the CP1 representation yields a new self-generated $`U(1)`$ gauge field $`A`$ coupling holons and spinons ( This field is analogous to the $`U(1)`$ gauge field in the standard slave-particle approaches). Due to coupling to holons (fermions in our approach), this gauge field acquires an anomalous dissipation term, the “Reizer singularity” , which dominates the low-energy action for the transverse component of the gauge field $`A^T`$. For $`\omega ,|\stackrel{}{q}|,\omega /|\stackrel{}{q}|0`$ we have $`A^TA^T(\omega ,\stackrel{}{q})(\chi |\stackrel{}{q}|^2+i\kappa \frac{\omega }{|\stackrel{}{q}|})^1`$, where $`\chi t/\delta `$ is the diamagnetic susceptibility and $`\kappa \delta `$ the Landau damping. The interplay of the two different energy scales, the spinon gap and the holon induced anomalous dissipation, is the key factor in our interpretation of the MI crossover. For the temporal component in the same limit, we have $`A^0A^0(\omega ,\stackrel{}{q})(\gamma +\omega _p)^1`$, where $`\gamma `$ is the fermion density of states and $`\omega _p`$ the plasmon gap.

Now we consider the introduction of a magnetic field $`H`$ perpendicular to the plane. The Ioffe–Larkin rule , $`R=R_s+R_h`$, i.e., the observed resistivity is the sum of the holon and spinon contributions, can be generalized to this case. The external electromagnetic potential $`A_{e.m.}`$, corresponding to the constant magnetic field $`H`$ can couple with coefficient $`\epsilon `$ to spinons and $`1\epsilon `$ to holons. In principle, $`0\epsilon 1`$ is arbitrary. However, to be consistent with the requirement that the physical inverse magnetitic susceptibility should be the sum of the inverse of that of holons and spinons, i.e., $`\chi ^1=(\chi _s^{})^1+(\chi _h^{})^1`$, where $`\chi _s^{}`$ and $`\chi _h^{}`$ are the renormalized spinon and holon susceptibility, respectively (this relation can be derived in the same way, as the Ioffe-Larkin rule), we find $`\epsilon =\chi _h^{}/(\chi _h^{}+\chi _s^{})`$. This value was argued earlier using variational considerations. Replacing these quantities by unrenormalized values, we find $`1\epsilon \chi _s/(\chi _h+\chi _s)\frac{J}{t}\sqrt{\delta /|\mathrm{ln}\delta |}<<1`$. In the Coulomb gauge $`A_{e.m.}^0=0`$, the effective action for the gauge field $`A`$ can be written as:

$`𝒮_{eff}(A)={\displaystyle }dx^0d^2x[{\displaystyle \frac{i}{2}}[A^0(\mathrm{\Pi }_h^0+\mathrm{\Pi }_s^0)A^0`$ (1)
$`+(A^T+(1\epsilon )A_{e.m.})\mathrm{\Pi }_h^{}(A^T+(1\epsilon )A_{e.m.})`$ (2)
$`+(A^T\epsilon A_{e.m.})\mathrm{\Pi }_s^{}(A^T\epsilon A_{e.m.})]+{\displaystyle \frac{i\sigma _h(H)}{2\pi }}A^0ϵ_{ij}^iA^j],`$ (3)
where $`\mathrm{\Pi }_h^{0,},\mathrm{\Pi }_s^{0,}`$ are unrenormalized polarization bubbles due to holons and spinons, respectively, $`\sigma _h(H)`$ is the Hall conductivity due to holons. Note $`A_{e.m.}`$ appears in two places in this low energy effective action: one is simply a shift of the transverse component of the gauge field $`A^T`$ by $`(1\epsilon )A_{e.m.}`$ and $`\epsilon A_{e.m.}`$ corresponding to the minimal coupling to holons and spinons, respectively, while the other is a C.S. term due to parity breaking induced by $`H`$.

As remarked in , the leading effect of the integration over $`A_0`$ is the renormalisation of the diamagnetic susceptibility: $`\chi \chi (H)=\chi +\frac{\sigma _h^2(H)}{4\pi ^2\gamma }`$ in the $`A^T`$ effective action. The holon contribution $`R_h`$ can be evaluated using the Boltzmann equation, taking into account the classical cyclotron effect, as in, obtaining

$`R_h=R_h^0[1+({\displaystyle \frac{(1\epsilon )H\tau }{m_h}})^2],`$ (5)
$`R_h^0{\displaystyle \frac{m_h}{8}}{\displaystyle \frac{1}{\tau }}\delta \left[{\displaystyle \frac{1}{ϵ_F\tau _{imp}}}+({\displaystyle \frac{T}{ϵ_F}})^{4/3}\right],`$ (6)
where $`\tau `$ is the transport relaxation time, $`\tau _{imp}`$ the impurity scattering time and $`m_h\delta /t`$ the holon mass. The spinon contribution $`R_s`$ is evaluated here using the Kubo formula for the spinon current: $`R_s=lim_{\omega 0}\omega (\mathrm{Im}\mathrm{\Pi }_s^{}(\omega ))^1`$, where $`\mathrm{\Pi }_s^{}`$ denotes the transverse polarization bubble at $`\stackrel{}{q}=0`$, renormalized by gauge fluctuations. At large scales, for $`x^0|\stackrel{}{x}|`$, $`\mathrm{\Pi }_s^{}(x)`$ is approximately given by $`_\mu G(x|A)^\mu G(x|A)`$, where $``$ denotes the $`A`$-expectation value and $`G(x|A)`$ is the spinon propagator. Using the Fradkin representation it can be transformed into a gauge invariant form $`_\mu G(x|F)^\mu G(x|F)`$, where $`G(x|F)`$ in terms of a path-integral over 3-velocities, $`\varphi ^\mu (t)\dot{q}_\mu (t),\mu =0,1,2`$, is given by

$`G(x|F)i{\displaystyle _0^{\mathrm{}}}𝑑se^{im_s^2s}{\displaystyle 𝒟\varphi ^\mu e^{\frac{i}{4}_0^s\varphi _\mu ^2(t)𝑑t}}`$ (7)
$`{\displaystyle }d^3pe^{ipxip^2s}e^{i𝒬^{ij}(p,s,\varphi |s^{},\lambda )[F_{ij}(p,s|s^{},\lambda )+ϵ_{ij}\epsilon H]}`$ (8)
with

$`𝒬^{ij}(p,s,\varphi |s^{},\lambda )[\mathrm{\#}]={\displaystyle _0^1}𝑑\lambda \lambda {\displaystyle _0^s}𝑑s^{}{\displaystyle _0^s^{}}𝑑s^{\prime \prime }`$ (9)
$`(\varphi ^i(s^{})2p^i)(\varphi ^j(s^{\prime \prime })2p^j)[\mathrm{\#}]`$ (10)
and

$`F_{ij}(p,s|s^{},\lambda )=F_{ij}(x+\lambda {\displaystyle _0^s^{}}𝑑s^{\prime \prime }\varphi (s^{\prime \prime })2ps),`$ (11)
$`F_{ij}_iA_j_jA_i.`$ (12)
After integration over $`A^T`$, using the effective action, and integration over the 3-velocities in the gaussian approximation, and over $`p`$ and $`s`$ by saddle point for low T, one obtains at large scales

$`\mathrm{\Pi }_s^{}(x)`$ $``$ $`[{\displaystyle \frac{}{x^\mu }}[e^{i\sqrt{m^2\frac{T}{\chi (H)}f(\alpha )+\frac{\alpha ^2\epsilon ^2H^2}{3q_0^2}}\sqrt{x_0^2|\stackrel{}{x}|^2}}`$ (14)
$`e^{\frac{T}{4\chi (H)}q_0^2g(\alpha )\frac{(x_0^2|\stackrel{}{x}|^2)}{m^2}}{\displaystyle \frac{1}{\sqrt{(x_0^2|\stackrel{}{x}|^2)}}}]]^2,`$
where $`\alpha =\frac{q_0|\stackrel{}{x}|}{2}`$, $`q_0`$ is the momentum scale associated with the anomalous skin effect due to Reizer singularity: $`q_0\left(\frac{\delta ^2T}{t}\right)^{1/3}`$. $`f`$ and $`g`$ are functions describing the effect of gauge fluctuations and for a real argument, $`f`$ is monotonically increasing, vanishing quadratically near the origin and $`g`$ is monotonically decreasing vanishing at large argument. (See  for explicit expressions.) The integration over $`|\stackrel{}{x}|`$ and $`x^0`$ appearing in the Fourier transformation are evaluated by saddle point for $`|\stackrel{}{x}|`$, at large $`x^0`$, and with scale renormalization by principal part evaluation for $`x^0`$. The integrals are dominated by a complex saddle point at $`|\stackrel{}{x}|=2q_0^1e^{i\frac{\pi }{4}}`$ for $`\chi (H)q_0|m_s(T,H)|\stackrel{<}{}T`$, Im$`(m_s^2(T,H))\stackrel{<}{}m_s^2`$, where

$$m_s^2(T,H)=m_s^2i\left(\frac{cT}{\chi (H)}\frac{\epsilon ^2H^2}{3q_0^2}\right),$$
(15)
with $`c=if(e^{i\pi /4}`$), a constant with real part $`3`$ and a small imaginary part. For the range of physical parameters considered here ($`H\stackrel{<}{}`$ 100 Tesla), these bounds gave a temperature range validity lying between a few tens and a few hundreds degrees.

The saddle point produces the following effects: it induces “renormalization” of the spinon mass yielding a $`T`$ and $`H`$ dependent damping: $`m_s^2m_s^2(T,H)`$; it gives rise to an attraction between spinon and antispinon leading to the formation of a damped spin wave; it introduces a multiplicative renormalization of the correlation functions which for $`R_s`$ is given by $`Z(T,H)[m_s^2(T,H)]^{\frac{1}{8}}`$, where $`Z(T,H)=(c^{}\frac{T}{\chi (H)}q_0^3\frac{2}{3}\epsilon ^2H^2q_0^5)^{1/2}`$ with $`c^{}`$ a new constant $`f^{\prime \prime }(e^{i\pi /4})`$.

The final result in the range of $`T`$ described above is given by

$$R_sZ(T,H)\frac{|m_s(T,H)|^{1/4}}{\mathrm{sin}\frac{\mathrm{\Theta }(T,H)}{4}},$$
(16)
where $`m_s(T,H)|m_s(T,H)|e^{i\mathrm{\Theta }(T,H)}`$. Basic features of our formulas can be summarized as follows: for low $`T`$, the effect of the spinon gap is dominating ($`\mathrm{\Theta }0`$), leading to an insulating behaviour; at higher temperatures one finds a metallic behavior due to the dissipation induced by gauge fluctuations, contained in $`|m_s(T,H)|`$, that becomes the dominant effect. Therefore a MI crossover is recovered, decreasing the temperature. The minimum of $`R`$ as a function of $`T`$, $`T_{MI}(\delta ,H)`$, is decreasing with $`\delta `$ and increasing with $`H`$ (see the MR curve $`(R(H)R(0))/R(0)`$ in Fig. 1), in agreement with experiment. In the absence of magnetic field the crossover is determined by the interplay between $`m_s^2=\xi ^2`$ and $`T/\chi Tm_h\lambda ^2`$, with $`\xi `$ the magnetic correlation length, $`\lambda `$ the thermal de Broglie wave length. When $`\lambda \xi `$, the “peculiar” localization effect due to SRAFO is not “felt”, and a metallic behavior is observed. In the opposite limit $`\lambda \xi `$ the cuprate is insulating. The external magnetic field effectively reduces the thermal energy, or increases the thermal wavelength, so the crossover temperature goes up. The resistivity is diverging at $`T=\frac{\epsilon ^2H^2\chi (H)}{3cq_0^2}`$, which approaches $`T=0`$ as $`H`$ vanishes. This divergence is lying outside the region of validity of our formulas and should be considered as an artifact.

However, the shift in MI crossover temperature leads to a large positive (in–plane transverse) MR at low $`T`$ which is our main new result, and it was absent in the earlier treatments. The derived MR scales quadratically with $`H`$ (See Fig. 2) in agreement, in particular, with data on LSCO, away from the doping $`\delta =1/8`$ where the stripe effects dominate. As remarked in , the shift of $`\chi `$ induced by the C.S. term reduces the damping and the $`H^2`$-term due to minimal coupling acts in the same direction. In the region of $`T`$ where dissipation dominates this induces a reduction of resistivity, a tendency contrasted by the classical cyclotron effect. One then has two possible types of MR curves: one is always positive but it exhibits a knee below the crossover temperature between the mass gap and the dissipation dominated regions (See Fig. 3). This behavior can be compared with the one observed in LSCO reported in  and we find a reasonably good agreement. If, on the contrast, the quantum effects related to $`\sigma _h(H)`$ are sufficiently strong,

a minimum develops, eventually leading to a negative MR in some region around it. This is illustrated in Fig.1. We should point out that the large positive MR at low temperatures is foreseen in this theory for both cases.

A comment on the limit $`H0`$ is in order, where we recover the results of , in particular $`m_s^2(T,0)=m_s^2icT/\chi ,Z(T,0)\frac{1}{\sqrt{\delta }}`$. In this limit the resistivity exhibits an inflection point at temperature $`T^{}(\delta )200300`$ K, (found also in the experimental curves), above which the theoretical curves start to deviate strongly from the experimental data. We propose to interpret this inflection point as a mid–point of a crossover to a new “phase” where our MF treatment is not valid anymore. If we identify our $`T^{}(\delta )`$ with the crossover temperature $`T^{}`$ found in experiments, both MI crossover temperature $`T_{MI}(\delta )T_{MI}(\delta ,0)`$ and $`T^{}(\delta )`$ are found in a reasonable agreement with experimental data (in the range $`0.02<\delta <0.08`$), due to a delicate cancellation of doping dependences: $`\frac{T}{\chi m_s^2}\frac{T\delta }{t\delta |\mathrm{ln}\delta |}=\frac{T}{t\mathrm{ln}\delta }`$. $`R_s`$ in this limit can be written in terms of a dimensionless variable $`xcT/\chi m_s^2`$ apart from an overall factor $`\sqrt{|\mathrm{ln}\delta |}`$. Hence, if we neglect the contribution $`T^{4/3}`$ due to holons and define a “normalised resistivity” by $`\stackrel{~}{R}(RR(T_{MI}))/(R(T^{})R(T_{MI})),`$ this is a function only of $`x`$, thus exhibiting a “universal” behaviour, as observed in YBCO.

As a final remark it might be worthwhile to notice that the same $`U(1)\times SU(2)`$ approach is able to reproduce qualitatively  the behavior of the spin lattice relaxation rate $`(1/T_1T)^{63}`$ found in underdoped YBCO and a structure of Fermi arcs around $`(\frac{\pi }{2},\frac{\pi }{2}`$) in the spectral density detected by ARPES.

We thank J.H. Dai for collaboration in the early stage of this work, and Y. Ando and T. Xiang for very helpful discussions. The work of P. M. is partially supported by RTN Programme HPRN-CT2000-00131, while L.Y. acknowledges the partial support by INTAS Georgia 97-1340.


format: text

text:
# Efficient N-body Realisations of axisymmetric Galaxies and Haloes

## 1 INTRODUCTION

The visible Milky Way and other disk galaxies may be contained in extended dark-matter haloes of largely unknown morphology. Spherical haloes are typically assumed, but there are indications that this matter component may be significantly flattened, or even triaxial (see Zaritsky 1998, Sackett 1999, Rusin & Tegmark 2000). An important indicator is the non-isotropic distribution of satellites around disc galaxies, or the Holmberg effect (see Holmberg 1969). This may be a result of dynamical friction of non-spherically symmetric dark haloes acting on the satellites (Peñarrubia, Boily & Kroupa 2000).

Equilibria of self-gravitating galaxies are achieved for a given distribution function ($`\mathrm{𝖽𝖿}`$) based on three integrals of the motion. This approach faces a difficulty in that only the first two classical integrals (total energy $`E`$ and one component of angular momentum $`𝑱`$) are known, the final integral taking an assumed form (see e.g. Lupton, Gunn & Griffin 1987; Einsel & Spurzem 1999). A greater difficulty still is that galaxies are heterogeneous systems of several components, and hence are hardly represented adequately with a unique $`\mathrm{𝖽𝖿}`$. A self-consistent equilibrium may be constructed by resorting to numerical integration of the equations of motion. Barnes (1988) introduced an iterative scheme by which particles relax in a fixed potential representing a galaxy component. This component is then given a particle representation in near-equilibrium.

Hernquist (1990, 1993) has introduced a by-now standard scheme which solves moments of the collisionless Boltzmann equation using tailored distribution functions. Three-component galaxies so constructed have the advantage over the Barnes approach of closer agreement with the sought equilibrium. In practice his algorithm is highly efficient when embedding a galactic disc in spherical components. (Benchmarking examples will follow in subsequent sections.) However we have found that for the purpose of constructing axisymmetric (non-spherical) bulges and haloes, the scheme requires significantly more cpu-time, mainly because one must solve numerically for the components of the velocity dispersion locally in three directions. This introduces an $`N^2`$ scaling of the computational time. By contrast, in spherical symmetry the velocity dispersion is known immediately from a one-dimensional integral (cf.  below). The algorithm of Hernquist (1993) remains efficient for up to tens of thousand of particles, however the computational time becomes prohibitive for larger numbers. Given that particle numbers adequate to galactic systems by far out-stretch available computer power, algorithmic problems such as this require solutions.

In this contribution we by-pass the Boltzmann equations approach, by noting that galactic components follow a well-defined hierarchy of sizes and masses. This suggests a perturbative treatment of particle potentials from the largest scales, down. Equilibrium solutions to Boltzmann equations in spherical symmetry are first transformed to the desired morphology directly, by mapping the velocity ellipsoid to the density isocontours. We construct equilibrium multi-component axisymmetric galaxy models by perturbing the velocity field of individual components to account for the added background potential (or, embedding). The computational time required for the whole procedure scales with particle number as for spherically symmetric models. Thus multi-component axisymmetric equilibrium configurations with millions of halo particles are easily constructed. In the following Section we review the problem and motivate our approach. This is then implemented and tested using a grid code. Benchmarking and comparisons with Hernquist’s standard buildgal follow. We conclude with suggestions for applications and further improvements.

## 2 METHOD

### 2.1 Scaling of Computational Time

Approximate solutions based on moments of the Jeans equations have proved useful when constructing equilibrium galaxies. Hernquist (1993) has introduced a scheme whereby equilibrium galaxy components are first constructed in isolation, then zipped together by adding the gravity of the other components in turn. The modified binding energy of individual particles is matched by a similar increase in kinetic energy to preserve the structure of each component.

The case when all components are spherically symmetric allows a simplification. Then the added binding energy is measured from the increase in mass inside a particle radius $`r`$. To add together two or more spherical components only requires to sort the particles radially. The computational time required by a scheme such as quicksort scales no faster than the $`3/2`$-power of particle number (see e.g. Press et al. 1992, §8).

For non-spherical potentials however, the modified velocity field (and square velocity dispersions) requires the computation of the local potential (and its spatial derivatives) for each particle. The exact, direct-summation algorithm to do this requires a number of operations increasing in proportion to $`N^2`$, the square particle number. The particle realisation of the density itself requires $`N`$ operations. Thus in general we may write the total computational time to realise a compound equilibrium of $`c`$ components as

$$T=\underset{i=1}{\overset{c}{}}\left(\stackrel{\mathrm{Density}}{\stackrel{}{N_it_{den,i}}}+\underset{\mathrm{Spherical}\mathrm{or}\mathrm{aspherical}}{\underset{}{N_it_{vel,i}\mathrm{or}N_i^2t_{vel,i}}}+\stackrel{\mathrm{Background}\mathrm{potentials}}{\stackrel{}{N_i\underset{j=1}{\overset{c}{}}N_jt_{cor,i}}}\right)$$
(1)
where $`N_i`$ is the number of particles for component $`i`$; $`t_{den},t_{vel}`$ the computational times per particle to realise the density and velocity fields for respective components; and $`t_{cor,i}`$ the characteristic time for adjusting the velocity field of $`i^{th}`$-component particles due to non-spherical potentials. For three-component galaxies made of a spherical bulge (b), disc (d) and aspherical halo (h), we would write

$$T=N_ht_{den,h}+N_h^2t_{vel,h}+N_h(N_d+N_b)t_{cor,h}+N_dt_{den,d}+N_d^2t_{vel,d}+$$

$$N_d(N_h+N_b)t_{cor,d}+N_bt_{den,b}+N_bt_{vel,b}+N_b(N_h+N_d)t_{cor,b}.$$
(2)
Clearly massive haloes call for large particle numbers and so the leading terms in (2) are $`N_h^2`$ and $`N_hN_d`$. We wish to remove these dependencies on $`N_h`$.

In any multi-component galaxy, the largest-scale structure would be closest to equilibrium self-gravity. This suggests that we look for a general scheme whereby the relative masses and scales of components are taken into account. Figure 1 provides a cartoon representation of our approach. Since details of the small scale components are lost on the much larger ones, a correction to the velocity field of larger structures would only involve truncated expansions of the small-component potentials. By contrast, a light structure requires more accurate treatment of the surrounding mass distribution: this is shown as light and heavy arrows on Fig. 1. We adopted this strategy by solving for the self-gravitating aspherical halo through a map of the velocity field of the spherical solution. As we shift down the scale ladder to smaller structures, higher-order expansions provide accurate corrections to the velocity field of individual components. Thus feedbacks from disc and bulge on the halo velocity field are given a low-oder series-expansion treatment, while the velocity field of disc and bulge particles are adjusted with a high-order expansion of the halo potential. Below we take a stepwise angle to implementing this algorithm.

### 2.2 Self-gravitating axisymmetric equilibria

The basic step consists in building axisymmetric systems from a transformation of spherically symmetric equilibria. If we consider as given the geometric aspect ratio and the mass profile of the equilibrium sought, it is then sufficient to construct a velocity field to match the gravity. In cylindrical coordinates $`(R=\sqrt{x^2+y^2},z,\theta )`$, the orbits of stars in an axisymmetric harmonic potential are recovered from integrating two decoupled equations of motion. In the case of a uniform oblate spheroid of density $`\rho `$, axes $`a,c;a>c`$, eccentricity $`e=\sqrt{1c^2/a^2}`$, the potential inside the bounding volume is (Table 2-1 in Binney & Tremaine 1987, hereafter BT+87; Chandrasekhar 1969, §3)

$$\mathrm{\Phi }_{sp}(R,z)=\pi G\rho \left(\mathrm{I}(e)a^2A_1(e)R^2A_3(e)z^2\right)$$
(3)
with known algebraic expressions for $`\mathrm{I},A_i`$. The motion of a test particle is found from solving for

$$\ddot{x}_i=_i\mathrm{\Phi }=2\pi G\rho A_ix_i,i=1\mathrm{}3$$
(4)
which admits a periodic solution in each component. Averaging the square velocity and position over one cycle, and substituting for $`A_i`$ from Table 2.2 of BT+87, we obtain

$$\frac{<v_3^2>}{<v_1^2>}=\frac{A_3}{A_1}\frac{<x_3^2>}{<x_1^2>}=\frac{A_3}{A_1}\left(1e^2\right)=2\frac{a/c\mathrm{sin}^1(e)/e}{\mathrm{sin}^1(e)/ec/a}\frac{c^2}{a^2}.$$
(5)
From (3) we find lines of constant potential intersecting the $`z=0`$ plane and the z-axis such that for these coordinates $`A_1x_1^2=A_3x_3^2`$. Defining the eccentricity of isopotential contours

$$e_\mathrm{\Phi }\sqrt{\left(1\frac{x_3^2}{x_1^2}\right)_\mathrm{\Phi }}=\sqrt{1\frac{A_1}{A_3}}$$
(6)
we may relate $`e_\mathrm{\Phi }`$ to the geometric eccentricity $`(e)`$ and that of the velocity ellipsoid, $`e_v`$, using (5) and (6)

$$e_v\sqrt{1\frac{<v_3^2>}{<v_1^2>}}=\left[\frac{e^2e_\mathrm{\Phi }^2}{1e_\mathrm{\Phi }^2}\right]^{{\scriptscriptstyle \frac{1}{2}}}.$$
(7)
The three quantities satisfy $`e_\mathrm{\Phi }<e_v<e`$, and hence the velocity ellipsoid is never as flat as the mass distribution that gives rise to it. This result may be derived directly from the virial theorem applied to the system as a whole (see BT+87, §4.3). Therefore, for the harmonic potential, (7) applies equally to averaged quantities and individual orbits.

In practice galaxies and haloes show peaked density profiles and hence the exact result (7) has limited application. An alternative to (7) is obtained by considering a star on a loop orbit in the mean potential of all the others. Since the matter is concentrated around the centre of gravity, a series expansion will be adequate when the orbit avoids the central region. Consider therefore the potential exterior to a homogeneous spheroid, of axes $`(a,c)`$ equal to the mass-weighted means of the peaked distribution. This is written in a form analogous to (3) (BT+87, Chandrasekhar 1969),

$$\mathrm{\Phi }_{sp}(R,z)=\pi G\rho \frac{a^2c}{(a^{})^2c^{}}\left(\mathrm{I}(a^{},c^{})(a^{})^2\underset{i}{}A_i(a^{},c^{})x_i^2\right).$$
(8)
In terms of $`(a,c)`$ we have

$$\frac{R^2}{(a^{})^2}+\frac{z^2}{(c^{})^2}=1;a^{}\sqrt{a^2+\lambda },c^{}\sqrt{c^2+\lambda }.$$
(9)
For oblate spheroids the solution $`\lambda (a,c,r)`$ is expressed in algebraic form:

$$2\lambda =r^2a^2(2e^2)+\sqrt{(a^2c^2)^2+2a^2e^2(z^2R^2)+r^4},$$

where $`r^2=R^2+z^2`$. (A relation similar to this one also exists for prolate spheroids but is omitted for reasons of clarity.) The solution $`\lambda >0`$ applies outside the bounding surface of the spheroid. For a loop orbit with $`r>a`$ a series expansion of (8) leads to (cf. Goldstein 1980, eq. 5-88)

$$\mathrm{\Phi }_{sp}(R,z)=\frac{GM_o}{r}+\frac{GM_o}{2r^3}(I_zI_R)P_2(\mathrm{cos}\theta )+\mathrm{O}([a/r]^4)$$
(10)
where $`\mathrm{tan}\theta =z/R`$, $`P_n(x)`$ is a Legendre polynomial, and $`I_x`$ an eigen-component of the inertia tensor per unit mass and we consider only the half-space $`z>0`$ in what follows. Integrating the equations of motion for the potential (10) leads to solutions of the form

$$x_i(t)=x_{o,i}\mathrm{cos}[w_1(tt_o)]\left(1+\frac{w_2^2}{w_1^2}\frac{x_{o,i}^2}{8r^2}\mathrm{cos}^2[w_1(tt_o)]+\mathrm{}\right)$$
(11)
where

$$w_1^2=G\overline{\rho }+\frac{9}{2}G\overline{\rho }\frac{I_zI_R}{r^2};w_2^2=\frac{9}{2}G\overline{\rho }\frac{I_zI_R}{r^2};\overline{\rho }\frac{M}{4\pi r^3/3}.$$

Thus when $`ra>c`$, $`w_20`$ with $`1/r^2`$ since the $`I_i`$’s are constants. The orbits describe ellipses of the same eccentricity in both $`\dot{𝒙}`$ and $`𝒙`$ . Note when $`I_z=I_R`$, (11) reduces to solutions for the right-hand side in (4). For such orbits we would therefore set

$$e_v=e.$$
(12)
Our solutions (7) and (12) bracket the range of possible aspect ratios the velocity ellipsoid may assume in equilibrium. Realistic galaxy equilibria will fall somewhere between strict homogeneity and point-mass distributions. We therefore sought a sensible interpolation for $`e_v`$ in the range $`e_\mathrm{\Phi }`$ to $`e`$ to account for this. The number of possibilities is a priori endless. We chose for each particle $`i`$

$$e_{v,i}^2=e_\mathrm{\Phi }^2+(e^2e_\mathrm{\Phi }^2)\sqrt{1\frac{<r^2>}{r_{g,i}^2}}$$
(13)
which was found to yield adequate equilibria. In computing $`e_v`$, we substitute in the above the mass-weighted average $`r^2`$ inside the spherical volume $`4\pi r_{g,i}^3/3`$. To obtain an eccentricity for the velocity ellipsoid that reflects the character of the gravitational field sampled by the star over one orbit, we compute the gravitational radius

$$r_{g,i}=\frac{GM}{E_i},$$

where $`E_i`$ is the binding energy of the $`i^{th}`$ star, and $`M`$ the total mass of the system. The velocity ellipsoid computed from (13) is the same for all stars of the same binding energy. We give some motivation for the interpolation (13) through an example.

To fix ideas, consider a Dehnen model for spherical galaxies. Dehnen (1993) introduced a family of potential-density pairs which provide a suitable fit to bulges and cuspy ellipticals :

$$M(r)=M_o\left(\frac{r}{r+r_c}\right)^{3\gamma },$$

$$\mathrm{\Phi }(r)=\frac{GM_o}{r_c}\times \{\begin{array}{cc}\frac{1}{2\gamma }\left(1\left[\frac{r}{r+r_c}\right]^{2\gamma }\right)\hfill & \hfill (\gamma 2)\\ \mathrm{ln}\frac{r}{r+r_c}\hfill & \hfill (\gamma =2)\end{array}.$$
(14)
Here $`(r_c,\gamma )`$ are two parameters and the total system mass $`=M_o`$. The case $`\gamma =1`$ corresponds to a spherical Hernquist (1990) density profile. Let $`r_c=1/10`$, and we truncate the distribution at $`r=1`$ for convenience. With $`M_o=1.21..`$ the mean square radius $`<r^2>=0.349\mathrm{}`$, enclosing 63% of the mass. We invoke the virial theorem and set $`E_i=\mathrm{\Phi }_i/2`$ for a given stellar orbit. Using this in (14) we compute $`r_{g,i}=GM_o/(\mathrm{\Phi }_i/2)=2(r_i+r_c)`$; the ratio of averaged values $`<r_g>/<r>=1/4`$. If we now squash the sphere down the z-axis to achieve an aspect ratio $`c/a=1/2`$, we compute from (13) an eccentricity $`e_v0.860..`$, close to $`e=\sqrt{1c^2/a^2}=0.866..`$. Setting $`\gamma =0`$ so that the central region is of uniform density, we find from (7) for this model $`e_v=0.832..`$, nearly unchanged from the previous value. This is because only $`12\%`$ of the mass now resides inside $`r=r_c`$. Thus unless the galaxy shows a broad uniform density region, the solution (13) on the mean lies close to $`e`$.

## 3 Stand-alone Axisymmetric Haloes and Bulges

We now implement the scheme described above. We begin by summarising the steps used in constructing spherical equilibria. Applications will centre on the Dehnen models (14), as they cover a range of central density gradients (through the index $`\gamma `$) from which to test the method. We concentrate on two aspects, namely the stability of equilibria so constructed and the range of cases to which this algorithm may be applied. We consider isolated self-gravitating haloes first; multi-component models are discussed in the next section.

### 3.1 Spherically Symmetric Equilibria

For the case of spherically symmetric mass profiles the potential at radius $`r`$ is

$$\mathrm{\Phi }(r)=_{\mathrm{}}^r\frac{GM(x)}{x^2}\text{d}x=_{\mathrm{}}^rG\frac{\text{d}x}{x^2}_0^x4\pi u^2\rho (u)\text{d}u$$
(15)
provided the run of mass $`M(r)`$ increases more slowly than $`r`$ at large distances. The mass distribution is a given quantity of the problem. A particle representation of any mass profile $`M(r)`$ can be obtained straightforwardly, first by filling a spherical volume uniformly by Monte Carlo method, then mapping the mass shells to the desired profile, i.e. $`M_h(\widehat{r})=M(r)`$, with $`M_h`$ the uniform-density spherical profile. The solution $`r(M)`$ is not known in analytic form in general and must be tuned by interpolation when only a discrete representation $`M_i(r_i)`$ is given on input. One advantage over solving directly by Monte Carlo method for $`r(M)`$ is that a unique realisation $`M_h(\widehat{r})`$ may be transformed to a variety of profiles. This essentially keeps any bias due to $`\sqrt{N}`$ fluctuations the same for different $`M(r)`$. This is especially important in systems where the effect of low-N statistics may be felt, such as in central galactic cusps. Interpolation errors would wipe out this advantage when the desired mass profile is known only on a coarse mesh $`\{r_i\}`$.

The problem is complete once the velocity field satisfying the Jeans equations under $`\mathrm{\Phi }(r)`$ is solved for. From the first of these equations (cf. Binney & Tremaine 1987, §4.2),

$$\frac{\text{d}\rho \overline{v_r^2}}{\text{d}r}+\frac{\rho }{r}\left[2\overline{v_r^2}\overline{v_\theta ^2}\overline{v_\varphi ^2}\right]=\rho \frac{\text{d}\mathrm{\Phi }}{\text{d}r},$$
(16)
for an isotropic velocity field the solution takes the form

$$\overline{v_r^2}=\frac{1}{\rho (r)}_r^{\mathrm{}}\rho (x)\frac{\text{d}\mathrm{\Phi }}{\text{d}x}\text{d}x=\frac{1}{\rho (r)}_r^{\mathrm{}}\rho (x)\frac{GM(x)}{x^2}\text{d}x.$$
(17)
With a known velocity dispersion at each radius it is customary to seek closure by picking a form for the velocity $`\mathrm{𝖽𝖿}`$. We will assume a velocity dispersion taking locally a Maxwellian form (Hernquist 1993),

$$F_v(𝒗,r)\text{d}d^3𝒗=4\pi \left(\frac{1}{2\pi \sigma ^2}\right)^{\frac{3}{2}}v^2\mathrm{exp}(v^2/2\overline{v_r^2})\mathrm{d}v.$$
(18)
N-body realisations of a Maxwellian distribution would require integrating to infinity in velocity space. In practice bound particles all have velocities below the local escape velocity $`(=\sqrt{2\mathrm{\Phi }(r)})`$, effectively setting an upper-limit for $`v`$ in (18). A proper implementation of this renormalises the velocities to preserve the Maxwellian character of the field, so

$$\overline{v^2}(r)=3\sigma ^2(r),$$
(19)
where $`\sigma `$ is the one-dimensional velocity dispersion at $`r`$ and $`\sigma ^2=\overline{v_r^2}`$ for non-rotating distributions.

### 3.2 From spheres to axisymmetry

To construct axisymmetric halo and bulge components, we proceeded as follows. Our starting point is a spherical equilibrium constructed as in the above. With $`M(r)`$ known, the potential $`\mathrm{\Phi }(r)`$ and velocity dispersion are evaluated for member stars from (15) and (17). Each particle is then given a velocity drawn from the distribution (18). Since only one-dimensional integrals are involved, the computational requirements for this part of the calculation scales with $`N`$, where $`N`$ is the total number of particles.

We then perform a homologous transformation of the position of each particle to achieve the desired morphology. The case where isodensity contours are concentric oblate spheroids is a good starting point. The particles are then mapped according to $`(x,y)(x,y)`$, $`zz\sqrt{1e^2}`$. The quantities $`<r^2>`$ and $`r_g`$ are computed directly from $`M(r)`$, prior to the mapping.

To modify the velocity field requires evaluating the potential energy of each particle in the new spheroidal potential, $`\mathrm{\Phi }_{sp}(R,z)`$. The quantity $`\sqrt{<r^2>}`$ is computed from mass lying inside any particle radius, $`r`$, for which we define a uniform density spheroid of axes $`(\sqrt{<r^2>},\sqrt{<r^2>}\sqrt{1e^2})`$, total mass $`M(r)`$. Equations (8) and (9) then yield the correct potential at $`(R,z)`$ up to quadrupolar order. By virtue of the virial theorem half of the potential energy accrued by the stars through the coordinate transformation should be invested in kinetic energy. Thus

$$T_i=E_i\mathrm{\Phi }(r_i)T_i^{}=T_i+\frac{\mathrm{\Phi }(r_i)\mathrm{\Phi }_{sp}(R_i,z_i)}{2}$$
(20)
for each particle $`i`$. The velocity vector $`𝒗_i^{}`$ must satisfy (13), which is the case when

$$𝒗_i^{}=\sqrt{\frac{2T_i^{}}{2T_iv_z^2e_v^2}}\times (v_x,v_y,v_z\sqrt{1e_v^2}).$$

In the above the primed quantities refer to the transformed, final state. The multiplicative factor ensures the new velocity gives the kinetic energy (20). The new potential could have been evaluated from (10), by computing the inertia tensor inside the spheroid $`(r,r\sqrt{1e^2})`$, for each particle. However (9) and (8) require only evaluating trigonometric functions and the algorithm is therefore optimised for the transformation we have performed here. When the eccentricity of the isodensity contours vary with position, one requires averaged eccentricities in (13) for every mass shell, which is provided through (10) by the tensor of inertia.

### 3.3 Numerical Setup

The algorithm we have described was tested using the nested-grid integrator superbox (Fellhauer et al. 2000). The tests we conducted used in total $`N=50000`$ particles. Since spherical equilibria provide the work horse of the scheme, we first constructed a spherical equilibrium to tune up the code.

The spherical profile used is the Dehnen family (14) with $`\gamma =1`$ and $`a=1/10`$. The models were all truncated at $`R_{cut}=1`$. We then evolved the distributions for 3 to 4 circular orbits at the half-mass radius, each revolution corresponding to time intervals $`4t_{\mathrm{cr}}=4\sqrt{3\pi /16(G\overline{\rho })^1}`$, where $`t_{\mathrm{cr}}`$ is the dynamical time and $`\overline{\rho }`$ the averaged density at the half-mass radius. Setting $`G=2`$ we have $`4t_{\mathrm{cr}}=0.727..`$ for a single revolution about the half-mass volume. At the edge of the system an orbit takes $`2.05..`$ units of time for a revolution.

superbox is a time-centred leap-frog integrator. The timestep we selected $`=0.006\times t_{\mathrm{cr}}`$; the radius where this timestep matches the local dynamical time for an Dehnen $`\gamma =1`$ model is at (roughly) $`1\times 10^3\times a`$, enclosing $`10^4`$ percent of the mass. The dynamics is therefore well resolved in time for all particles in the simulations. In general the central cusp is poorly resolved by grid codes due to finite linear resolution. However the three levels of resolution allowed by superbox means that we can focus on the central region and boost resolution as needed. Setting the total grid size $`=1/2`$ we then find for an 32-mesh grid a linear resolution $`l1/2\times 1/32=0.016`$ such that 2 percent of the total mass falls inside the innermost grid points. This would minimise structural evolution, however the coarse central density due to finite resolution means that velocity field and potential no longer match there, resulting in a small degree of relaxation at the centre.

The left-hand panels on Fig.2 show the results of evolving a spherical Dehnen $`\gamma =1`$ model. The figure displays the 10-percentile Lagrange radii of the mass sorted in spherical shells. We note that systematic fluctuations are low for the innermost mass shells. However near the edge appreciable evolution can be seen, at a level of 15% for the 90% mass radius. Bearing in mind this caveat, we then applied the same numerical setup to flattened Dehnen models. To preserve the good resolution of the dynamics in time, we enforced identical central density for all models by an homologous transformation of the particle positions, so the dynamical time $`t_{\mathrm{cr}}`$ remains the same.

### 3.4 Results for axisymmetric models with $`\gamma =1`$

Non-rotating self-gravitating equilibria of aspect ratios $`c/a>1/3`$ are susceptible to bending modes of instability (Combes et al. 1990, Sellwood & Wilkinson 1993). No elliptical galaxy has a projected aspect ratio flatter than E7 (Binney & Merrifield 1998). Furthermore the narrow stream of debris from the Sagittarius dwarf suggests that the halo of the Milky Way may be as round as or rounder than $`E1`$ (Ibata & Lewis 1998). Therefore we consider a flattened model with $`c:a=1:3`$ as a limiting case to test the algorithm. This is done first, before investigating the effects of varying aspect ratio and power index $`\gamma `$.

Figure2(b) graphs the Lagrange radii of an 1:3 flattened Dehnen model. The mass profile has been binned in spherical shells to facilitate comparisons with the spherical case. Assessing the overall setup, we find as expected a degree of relaxation in the innermost mass bins. Close inspection shows fluctuations of amplitude $`0.01`$ forming at all Lagrange radii. These fluctuations propagate inside out, and result from stars leaving the central region. For the model at hand we estimate that $`2.8\%`$ of the total mass gives rise to them. At the end of this simulation only 0.52% of stars had left the volume of the simulation, therefore we conclude from this exercise that $`2\%`$ of kinetic energy is not assigned correctly by our application of (20) and leads to large errors in the binding energy of relatively few particles.

We monitor in time the morphology of the model defining the mass-weighted principal axes

$$a_i^2\frac{I_j+I_kI_i}{2}$$
(21)
where $`\{I_i\}`$ are the eigenvalues of the inertia tensor and $`\{i,j,k\}=\{x,y,z\}`$ and permutations thereof (Goldstein 1980, §5-3; see also Kroupa 1997). Particles are first selected in spherical Lagrangian shells, from which the $`I_i`$’s are computed. A first estimate of the principal axes of the mass distribution is computed from (21). A second selection of particles is then made, with the bounding volume now assuming the ellipsoidal shape defined by the $`a_i`$’s. The inertia tensor and eigenvalues are recomputed from this new selection to obtain the final estimate of the system axes. Figure 3 (left-hand panel) graphs the evolution of the quantity $`2a_3/(a_1+a_2)`$ for two mass bins, the 30% and 60% Lagrange radii. For clarity the 30% results were divided by a factor two. We find the relaxed configuration giving nearly the exact same quantities after a half-unit of evolution. It is noteworthy that relaxation was more severe in the central region, where the computed aspect ratio drops from 0.18 initially to 0.16 at later times. The right-hand panel shows the velocity field for the 60% innermost particles. Here there are virtually no signs of evolution or trends beyond the level of root-$`N`$ fluctuations.

### 3.5 Models with $`\gamma 1`$

We investigated the range of applicability of our approach by evolving Dehnen models with different power indices $`\gamma `$. The solution (7) for the velocity ellipsoid applies to harmonic potentials of uniform density. Since we opted instead to implement (13), we ask whether equilibria can be setup this way with a harmonic central core, i.e. for models with $`\gamma =0`$.

Figure 4 (left-hand panel) illustrates the impact of the central harmonic core. We have graphed the same quantities as on Fig. 3, using the same numerical setup to validate the comparison. The oscillations of the 60% axis ratio seen early on in the simulation compare in magnitude to those observed for the $`\gamma =1`$ case. However in the inner region the ratio varies from 0.18 initially, to 0.15 after one time unit of evolution, or down 17%, somewhat more than what was observed for the same volume in the case of $`\gamma =1`$, where variations were of $`11\%`$ relative magnitude.

The multipolar expansion of the gravitational field suffers less from truncation if the matter is more concentrated about the centre of gravity. Therefore Jaffe’s (1983) profiles (Dehnen $`\gamma =2`$) provide another test since they are more concentrated. The right-hand panels on Fig. 4 show indeed that fluctuations due to dynamical relaxation are much reduced for the period $`t<1/2`$. The principal axes ratio fluctuates between $`0.17\pm 0.01`$ or 5.9% for the inner 30% mass radius. At larger volume it fluctuates initially by a similar relative fraction. We would expect this property to prove helpful when modelling dark haloes as isothermal gas, for which asymptotically $`\rho r^2`$ at large distances.

## 4 Multi-component Models

The map (13) provides equilibria for one-component self-gravitating spheroidal galaxies. We now generalise to multi-component galaxies. Our goal is to verify that the embedding of new components does not introduce further evolution in the structure of the galaxy. In an effort to bring out algorithmic, as opposed to dynamical, biases, we have put together a bulge-disc-halo model galaxy where the bulge and halo have the same spherical Dehnen $`\gamma =1`$ mass profiles. We then opted to combine a spherical bulge, exponential disc and axisymmetric halo, requiring respectively zero-, high- and low-order corrections to the monopole for a precise determination of the potentials. We modified a version of buildgal kindly provided to us by L. Hernquist and based on his 1993 algorithm. We implemented the map of Section 2.2 for axisymmetric galactic potentials, taking into account the hierarchy of individual components.

### 4.1 Disc, bulge & halo galaxy

While not attempting to model the Galaxy, we set up model galaxies with relative masses and lengths in rough agreement with those of Freudenreich (1998) for the Milky Way (see Sackett 1997, Binney & Merrifield 1999, §10). In model units, the bulge has a mass and radius $`=1`$; the core length $`r_c=0.10`$. The exponential disc has a radial scale length $`h=1.0`$, scale height $`z_o=0.3`$, total radius = 10 and mass = 3.0. We chose scale lengths $`h`$ and $`z_o`$ large compared with the linear resolution $`l=0.016`$ in the inner region of our grid code.

The heavier spheroidal halo was given an 1:2 aspect ratio as described in Section 2.2. We took $`z`$ as the symmetry axis so the spheroid equator lies in the x-y plane. The halo has a total radius = 20.00, mass = 30.00, and we chose a halo core length $`r_c=2.0`$. The halo and disc masses are in the ratio 4:1 at the edge of the disc, so the halo drives the dynamics there but not overwhelmingly so. This is done on purpose, with a view to perturb the halo equilibrium, as it must adjust to the disc gravity.

### 4.2 Computing the potentials

Experience and pragmatic considerations led to the following simplifying tricks. When computing the disc feedback on embedded bulge particles, we have kept the direct-summation algorithm from the original code, since the dimensions of the bulge are comparable to the disc scale length. Hence a multipole expansion would require high-order series expansion which we deemed not essential, as the number of bulge particles is less or comparable to the number of disc particles, which already gives $`tN_d^2`$ scaling of computational time (cf. equation 2). We found the direct-summation computation, contributing a term $`N_dN_b`$ in the total computational time, gives better result than a expanding the disc potential to high-order, although a crude evaluation using (8) also provides a sensible equilibrium - see below. The feed-back of the bulge on disc particles however is accounted for already when the monopole term of its potential is added to that of the disc.

The halo and disc interactions require some care. We found it sufficient for an axisymmetric halo to expand its potential to quadrupolar term (cf. ) and adding this to the disc’s (and the bulge’s) self-gravitating potentials. Equation (10) is easily differentiated in three dimensions. The kinetic energy is then added to the disc and bulge particles for each degree of freedom in the same proportion as the components of the halo gravitational forces (to account for the work done in reshaping the system). In this way, relatively more kinetic energy is invested in the x-y plane, when compared with the z-direction, to fight off the enhanced gravity due to the flattened halo in the equatorial plane of the spheroid. As the rms disc aspect ratio is close to 1:12, the boost in disc kinetic energy required to maintain equilibrium was added to the circular motion at cylindrical radius $`R`$ for each disc particle. Since the streaming motion at $`R`$ is known, energy can be directly added to each particle once (20) has been computed.

As halo particles spend little time near or in the plane of the disc, we treated the disc as a highly flattened oblate spheroid, for which the potential is computed from (8). To do so we averaged over the vertical structure of the disc, however the original scales of height and radial decay are preserved if we use as spheroid axes the mass-weighted disc parameters $`<R>`$ and $`<z>`$. This effectively removes the term $`N_dN_h`$ of the total computational time (2).

### 4.3 Results

For comparison purposes we display two models on Fig. 5. The first one has both spherical bulge and halo, while the second has an oblate spheroidal halo of aspect ratio 1:2. The motivation for doing this is clear: for spherically symmetric components, the multipole approach becomes exact, since only the monopole term contributes to the potential of each component. Models of a disc galaxy with spherical bulge and halo thus minimise relaxation effects.

The structure of the bulge is well accounted for by sorting the mass on spherical Lagrangian shells. For the case of the disc, of mean aspect ratio $`1:12`$, the Lagrangian shells effectively follow the cylindrical radii. Contrasting the two models, with and without spherically symmetric halo, we find little effect on the radial structure of disc and bulge introduced by our approximate treatment. In both the cases, radial oscillations are seen to propagate from the inside, outwards, while taking several dynamical times to dampen. On the figure the arrows indicate the dynamical time $`t_{\mathrm{cr}}`$ at the half-mass radius (bulge and halo) or at $`R=h`$ (disc) of individual components in isolation. These radial oscillations give rise to transient patterns on visual inspection. We note that after a few revolutions of the disc at $`R=h`$, much of the fluctuations have disappeared or are much attenuated.

The moments of inertia give insight into the structure of individual galaxy components. On Fig. 6 we display the moments of inertia, principal axes and velocity dispersion components of the bulge and halo for the case where the halo is an oblate spheroid. We find as on Fig. 2 slight changes in the ratio of principal axes of the halo. However, after settling, the structure remains close to the one specified on input. The bulge, on the other hand, is clearly stretched out of spherical symmetry, due to the combined pull of the disc and halo: at the edge of the bulge, $`R=1`$, the halo mass $`10/3`$ is already several times that of the bulge, which is not fully self-gravitating anymore. The aspect ratio of the bulge in equilibrium is close to $`0.65/0.75=0.87`$, when we would want unity. The external forces acting on the bulge leads to expansion, as seen from the increased moments of inertia (top left-hand side panel). A quick solution to counter this effect would be to stretch the bulge initially to mildly prolate structure of axis ratio $`a:c1:1.1`$, so as to relax to near sphericity. We have not resorted to this therapy and chose instead to illustrate the limitations of the algorithm.

Turning to the vertical structure of the disc, we averaged the height $`|z|`$ of stars inside $`R=2h`$, since this region is well resolved by our grid code and errors due to numerical integration are reduced. We computed both average vertical height and the one-dimensional velocity dispersion $`\sigma _z`$ for ten snapshots covering ten time units of evolution, or four revolutions. Table 1 gives the results. Early in the evolution of the disc, the vertical structure shrunk by $`15\%`$, as the mean height drops from $`0.2`$ to $`0.17`$. We observed a decrease of the standard deviation as well, of the same magnitude, which leads us to conclude in a global contraction of the vertical disc structure. This is matched by a rise in velocity dispersion of the same relative magnitude (see Table 1). Since $`<v_z>0`$ and is much less than the vertical velocity dispersion, the bulk motion remains negligible thoughout.

### 4.4 Benchmarks

Approximate methods such as the one we have presented would not be of interest without a handsome pay-out at the computer. We find the scaling properties of the method with a series of numerical tests. Specifically, we first bench-marked the original version of the code buildgal for spherical bulge + disc galaxy models. We then constructed an axisymmetric bulge using the full algorithm, which computes the response of the disc to bulge potential exactly. Since bulge and halo are constructed in the same way in buildgal , and the disc embedding done in the same exact fashion, we may relabel the bulge as halo. This is so when referring to buildgal below. In a second set of tests, we repeated these experiments, using the new version of the code. The same default disc parameters were used in all our tests, and haloes were given Hernquist (1990) mass profiles with $`r_c=0.1`$, and truncated at the disc edge.

All benchmarks refer to an 333 MHz G3 Apple processor operating under S.u.S.E. Linux 6.4 OS; we compiled buildgal with the widely available GNU Fortran compiler g77 -O4. The performance were done with disc and bulge particle numbers, $`N_d,N_b`$, in the range $`5,000`$ to $`5\times 10^4`$ (disc); and $`5000`$ to $`N_b=5\times 10^5`$ (bulge or halo).

We solved for the operation times $`t_{x,i}`$ in (2), first as a set of linear equations by constructing stand-alone disc and bulge (no feedback terms present), then re-doing the calculations with disc and bulge together. Our results are listed in Tables 2 and 3. The most important point to notice is the large increase in computer time for models with an axisymmetric bulge, which scales effectively with the square of particle number. The turn-over to quadratic dependence is approximately $`\mathrm{2\hspace{0.17em}10}^4`$ particles. The validity of (2) is inferred from the convergence of the operation times $`t_{x,i}`$’s with increasing $`N`$. We have extrapolated to infinite particle numbers using the polynomial fit routine pzextr of Press et al. (1992). The numbers are truncated to the last significant digit according to error estimates obtained from the results of the two highest-order polynomial fits.

Repeating this exercise with the new algorithm, we concentrated on the case of an axisymmetric halo + disc, since the case of a spherical halo would give essentially the same $`t_{x,i}`$’s. We therefore have to compute only the extra computational time required for the transformation of the halo to axisymmetry and its coupling with the disc. We find as expected that the time required to adjust the disc and halo velocity fields scale in proportion to their respective particle numbers. A comparison of total computational time shows that axisymmetric compound models require only a modest increase in computer time when compared with spherically symmetric components. The new code, MaGalie (for Mache Galaxie), can thus construct galaxies with bulge, disc and halo, all in axisymmetry, at virtually no extra costs.

## 5 CONCLUSIONS

We have presented a method of constructing approximate equilibria in good agreement with those sought. In particular the algorithm scales linearly with bulge and halo particle numbers, when these are axisymmetric or spherical. The approach takes into account implicitly the relative scales of components involved, and treats the interactions between components by expanding the potential to required order (Fig. 1).

The errors introduced from treating non-spherically symmetric components was found to be no worse than the case where only spherical components surround a disc. The required computational time for this algorithm scales with particle numbers as

$$t_{\mathrm{cpu}}=N_ht_{den,h}+N_ht_{vel,h}+N_ht_{cor,h}+N_dt_{den,d}+N_d^2t_{vel,d}+N_dt_{cor,d}$$
(22)
where the parameters $`t_{x,i}`$ are given in Table 3 for a G3 processor (see Section 4.4). With this algorithm non-spherical haloes surrounding galaxies can be constructed as a matter of routine. The results of a parameter survey of orbital decay in aspherical haloes will be presented in a separate contribution (Peñarrubia et al. 2000).

We have limited our discussion to prolate spheroidal potentials. The approach is readily applicable to prolate or triaxial structures however. To construct triaxial equilibria, one would map the axes of the velocity ellipsoid as in (13) to the geometry of the triaxial body. Equations (8) and (10) are easily generalised to triaxial systems (BT+87). Likewise, net rotation of individual components is implemented by alignment of particle angular momenta (see e.g. Lynden-Bell 1962, Hernquist 1993). Models of barred galaxies may possibly be constructed by converting random into bulk motion of a triaxial component, but we have not implemented this.

The disc self-gravity now represents the bottleneck through the $`N_d^2`$ term in equation (22). For situations where the disc substructure can be neglected, this dependence may be removed by constructing a disc from flattening a spheroid to small aspect ratio, then converting random motion to provide rotational support. When the disc structure matters to the dynamics however, this approach is unlikely to yield sensible disc equilibria, because the vertical and radial motions in thin structures require fine-tuning the velocities locally, hardly in line with such a brute-force approach. Note however that one may construct a disc similar to the Milky Way’s by considering the thin and thick discs as separate components. These may then be combined in the fashion described in Section 2. Sellwood & Wilkinson (1993) discuss methods of constructing stable thin galactic discs comprehensively.

Acknowledgements CMB and JPG acknowledge support from the German government through an SFB 439 grant at the Karls-Ruprecht-Universität Heidelberg.

Note added in proof - After this article was submitted, we found a natural, quick-fix solution to the ‘$`N_d^2`$’ scaling of the thin disc construction algorithm (cf. Tables 2 and 3): it consists in adding together $`m`$ realisations of an $`n`$-particle disc to reach the desired total particle number $`N_d=m\times n`$. The total cpu then scales with $`m\times n^2n\times N_d`$, or linearly with $`N_d`$. Though $`n`$ and hence the coefficient of proportionality should be reasonably large to dampen numerical noise, the approach is well-suited to parallel algorithms.


text:
# Splitting between quadrupole modes of dilute quantum gas in a two dimensional anisotropic trap

## I introduction

In a series of experiments, Bose-Einstein condensates (BEC) have been produced by cooling a vapor of alkali atoms to a temperature of a few nanokelvin . This system opens up interesting perspectives in the field of many body physics. There has been much progress in the theoretical understanding of this system . In particular, the low-energy collective excitation spectrum of a Bose condensed dilute gas in a trap has been discussed analytically by Stringari  using a hydrodynamic approximation and also by sum-rule approach. A few low-lying excitations have also been calculated analytically by using time-dependent variational approach . The low-energy excitation spectrum obtained by using time-dependent gaussian variational ansatz exactly coincides with the hydrodynamic results in the limit of large particle number N. Also a similar type of scaling ansatz has been used to describe the time evolution of the condensate in the large N limit . Experimentally the low-lying collective excitation frequencies of a condensate have been measured both at zero temperature  and at finite temperature . These observed values of the collective oscillation frequencies are in agreement with theoretical results at zero temperature.

After the discovery of BEC in alkali atomic gas, the behaviour of trapped Fermi gas is also in focus. It is also possible to trap the Fermionic atoms at very low temperature, where the quantum effects can be observed. There has been experimental progress towards cooling a Fermi gas into the degenerate regime $`(T<T_F)`$ . Several authors have studied the thermodynamic properties , collective excitation frequencies in the normal phase  as well as in the superfluid phase  of a three dimensional trapped Fermi gas.

The reduction in dimension of a quantum system is the subject of extensive studies in trapped Bose systems  as well as trapped Fermi systems . With present technology one can freeze the motion of the trapped particles in one direction to create a quasi two dimensional quantum gas by tuning the frequency in the $`z`$ direction. Recently, the lower dimensional BEC has been realized .

In a deformed quantum system, angular momentum is not a good quantum number and the angular momentum states mix with each other. Also in the presence of a small deformation of trap, the states with angular momentum quantum number $`+l`$ and $`l`$ split. The partition function of a deformed two-dimensional harmonic oscillator is exactly the same as the partition function of a rotating harmonic oscillator, where the rotation frequency $`\mathrm{\Omega }`$ is related to the trap anisotropy, $`\mathrm{\Omega }=(\omega _y\omega _x)/2`$. Here $`\omega _x`$ and $`\omega _y`$ are the oscillator frequencies of the deformed oscillator. From this analogy we find that the splitting occurs between the $`+l`$ and $`l`$ angular momentum states and it is proportional to the trap anisotropy. Similarly for a quantum gas in a deformed trap, the degenerate multipole modes with angular momentum quantum number $`+l`$ and $`l`$ also split. The dipole excitation frequencies of a quantum gas in deformed harmonic trap are the trap frequencies along x and y directions, $`\omega _x`$ and $`\omega _y`$. The frequencies and the splitting of the dipole mode neither depends on the statistics of the trapped gas nor on the interaction between the particles. In circular symmetric trap two quadrupole frequencies are degenerate and the degeneracy is lifted up by elliptic trap deformation. In the present work we consider quadrupole modes of quantum gas in a quasi two dimensional anisotropic trap. We calculate the frequencies of quadrupole like modes for both bosons and fermions in a deformed harmonic trap analytically. For both statistics, we study the effect of interaction on the splitting between the quadrupole modes for arbitrary deformation of the trap. First we consider the case of a two dimensional deformed trapped Bose gas at zero temperature. We analyze the nature of small oscillations of the confined gas within the time-dependent variational method. In case of deformed trap the quadrupole mode couples with the monopole mode whereas the scissors mode remains decoupled. This is due to fact that the Hamiltonian is invariant under reflection. We calculate the frequencies of the mixed type of modes using variational technique. Next we construct the most general type of excitation operator which couples the quadrupole and monopole mode. Then we consider only the minimum energy excitation in sum rule approach. The excitation energy obtained from the sum-rule approach exactly matches with that obtained from the time-dependent variational technique even for an arbitrary number of particles. In the case of a symmetric trap, these modes can be identified as a monopole mode and a quadrupole mode. There is a possibility of experimental verification of our results since the quasi two dimensional Bose condensed has been realized in MIT .

This sum-rule method allows us to calculate the monopole mode and the quadrupole modes of a Fermi system in a deformed trap where the equation of motion technique does not hold. These modes can be observed in nano structures like quantum dots. When the two-body interaction is long ranged, the sum rule method can be generalised.

There has been no systematic theoretical study on the collective excitations of a two dimensional deformed trapped quantum gas at zero temperature. The purpose of this paper is to give an analytic description of the quadrupole excitation frequencies of a two dimensional deformed trapped quantum gas at zero temperature and to calculate the splitting of the quadrupole modes for an arbitrary deformation of the trap.

The paper is organised as follows. In Sec.II, we model the quasi-two dimensional trapped Bose system. Using the time-dependent variational method we calculate the monopole and quadrupole excitation frequencies of a two dimensional deformed trapped interacting Bose gas. In Sec.III, using the sum-rule approach we derive a general dispersion relation for the quadrupole excitation frequencies of a two dimensional trapped quantum system interacting through the two-body potential. This relation is valid for both Fermi and Bose statistics. We apply this general dispersion relation to calculate the same excitation frequencies of a trapped Bose system. We show that the quadrupole excitation frequencies obtained from both methods are exactly the same. In Sec.IV, we consider trapped unpolarized fermions and apply the dispersion relation obtained from sum-rule approach to calculate the frequencies of quadrupole modes. In Sec. IV, we present the summary and conclusions of our work.

## II  Collective low-energy excitation frequencies of a two dimensional deformed trapped Bose gas

In BEC experiments, the trap potential can be approximated by an effective three dimensional harmonic oscillator potential, with tunable trap frequencies $`\omega _z`$ in the axial (z) direction and $`\omega _x`$, $`\omega _y`$ in the transverse (x-y) plane. The alkali-metal vapors used in experiment are very dilute and the interparticale interaction is well described by the short range pseudopotential and the interaction strength is determined by s-wave scattering length $`a`$. Here we consider the case when the interparticle interaction is strongly repulsive. The Gross -Pitaevskii(GP)  energy functional of the trapped boson of mass $`m`$ is given by,

$`E[\psi ]`$ $`=`$ $`{\displaystyle }d^2rdz[{\displaystyle \frac{\mathrm{}^2}{2m}}|\psi (\stackrel{}{r},z)|^2`$ (1)
$`+`$ $`{\displaystyle \frac{m}{2}}(\omega _x^2x^2+\omega _y^2y^2+\omega _z^2z^2)|\psi (\stackrel{}{r},z)|^2`$ (2)
$`+`$ $`{\displaystyle \frac{g}{2}}{\displaystyle }d^2r^{}dz^{}\delta ^2(\stackrel{}{r}\stackrel{}{r}^{})\delta (zz^{})|\psi (\stackrel{}{r}^{},z^{})|^2|\psi (\stackrel{}{r},z)|^2],`$ (3)
where, $`g=\frac{4\pi a\mathrm{}^2}{m}`$, $`\stackrel{}{r}`$ is the position vector in $`xy`$ plane and $`\psi (\stackrel{}{r},z)`$ is the condensate wave function.

It has been shown by Baganato et al.  that for an ideal two dimensional Bose gas under harmonic trap, a macroscopic occupation of the ground state can exist at temperature $`T<T_c=\sqrt{N/\zeta (2)}\mathrm{}\omega /k_B`$. With present technology it is possible to freeze the motion of the trapped particles in one direction to create a quasi-two dimensional Bose gas. In the frozen direction the particles execute zero point motion. To achieve this quasi-two dimensional system, the frequency in the frozen direction should be much larger than the frequency in the x-y plane and the mean interactions between the particles. Alternatively, the trap frequencies are such that $`\mathrm{}\omega _z\mu \mathrm{}\omega _0`$ and $`k_BT\mathrm{}\omega _z`$, where $`\mu `$ is the chemical potential of the two dimensional Bose gas.

For a quasi two dimensional system we may assume that the wave function in the z-direction is seperable and is given by,

$$\psi (z)=\frac{1}{(\sqrt{\pi }a_z)^{1/2}}e^{\frac{z^2}{2a_z^2}},$$
(4)
where $`a_z=\sqrt{\mathrm{}/m\omega _z}`$ is the oscillator length in the z-direction. Now we integrate out the z-component in the three dimensional GP energy functional, then we get the effective energy functional in two dimension:

$`E{\displaystyle \frac{\mathrm{}\omega _zN}{2}}`$ $`=`$ $`{\displaystyle }d^2r[{\displaystyle \frac{\mathrm{}^2}{2m}}|\psi (\stackrel{}{r})|^2`$ (5)
$`+`$ $`{\displaystyle \frac{m}{2}}(\omega _x^2x^2+\omega _y^2y^2)|\psi (\stackrel{}{r})|^2`$ (6)
$`+`$ $`{\displaystyle \frac{g_2}{2}}{\displaystyle }d^2r^{}\delta ^2(\stackrel{}{r}\stackrel{}{r}^{})|\psi (\stackrel{}{r}^{})|^2|\psi (\stackrel{}{r})|^2],`$ (7)
where $`g_2=2\sqrt{2\pi }\mathrm{}\omega _za_za`$ is the effective coupling strength in two dimension, $`a`$ is the $`s`$-wave scattering length in three dimension and N is the total number of particles in the condensate. The same effective coupling constant is obtained in Ref. . The effective interaction in two dimension is given by,

$$V_I=g_2\delta ^2(\stackrel{}{r}\stackrel{}{r^{}}).$$
(8)
The chemical potential of quasi two dimensional Bose condensed state is $`\mu =\mathrm{}\omega _0\sqrt{\sqrt{2\pi }Na/a_z}`$. Recently, the two dimensional Bose condensed state has been realized in MIT . In this experimental set up, they have loaded $`N10^4`$ number of sodium atoms in a trap with trap frequencies $`\omega _z/2\pi =790`$ Hz, $`\omega _0/2\pi =\sqrt{\omega _x\omega _y}20`$Hz and $`a=2.75`$ nm. One can easily calculate the chemical potential $`\mu 0.19\mathrm{}\omega _0\sqrt{N}`$ which satisfies the above mentioned inequality condition to be a quasi two dimensional Bose system.

In two dimensions, the equation of motion of the condensate wave function is described by the Gross-Pitaevskii equation,

$$i\mathrm{}\frac{\psi (\stackrel{}{r})}{t}=[\frac{\mathrm{}^2}{2m}^2+V(\stackrel{}{r})+g_2|\psi (\stackrel{}{r})|^2]\psi (\stackrel{}{r}),$$
(9)
where $`V(\stackrel{}{r})=\frac{1}{2}m(\omega _x^2x^2+\omega _y^2y^2)`$ is the deformed trap potential in two dimension. The normalization condition for $`\psi `$ is $`d^2r|\psi |^2=N`$. N is the number of particles in the condensate . One can write down the Lagrangian density corresponding to this system as follows:

$``$ $`=`$ $`{\displaystyle \frac{i\mathrm{}}{2}}(\psi {\displaystyle \frac{\psi ^{}}{t}}\psi ^{}{\displaystyle \frac{\psi }{t}})`$ (10)
$`+`$ $`({\displaystyle \frac{\mathrm{}^2}{2m}}|\psi |^2+V(\stackrel{}{r})|\psi |^2+{\displaystyle \frac{g_2}{2}}|\psi |^4),`$ (11)
where * denotes complex conjugation. One can get the non-linear Schroedinger equation (9) by minimizing the action related to the above Lagrangian density (10). In order to obtain the evolution of the condensate we assume the most general Gaussian wave function,

$$\psi (X,Y,t)=C(t)e^{\frac{1}{2}[\alpha (t)X^2+\beta (t)Y^2+\gamma (t)XY]},$$
(12)
where C(t) is the normalization constant. X and Y are the dimensionless variables, $`X=\frac{x}{a_0}`$, $`Y=\frac{y}{a_0}`$ where $`a_0=\sqrt{\frac{\mathrm{}}{m\omega _0}}`$ is the oscillator length and $`\omega _0=\sqrt{\omega _x\omega _y}`$ is the mean frequency. Further, $`\alpha =\alpha _1+i\alpha _2`$, $`\beta =\beta _1+i\beta _2`$ and $`\gamma =\gamma _1+i\gamma _2`$ are the time dependent dimensionless complex variational parameters. The $`\alpha _1`$ and $`\beta _1`$ are inverse square of the condensate widths in $`x`$ and $`y`$ direction respectively. The square of the normalization constant is $`|C(t)|^2=\frac{N\sqrt{D}}{\pi a_0^2}`$, where $`D=\alpha _1\beta _1\gamma _1^2`$. The Gaussian ansatz (eq. (12)) for the order parameter can also be generalized to three dimensional anisotropic trapped Bose system to study various scissors modes.

The gaussian variational ansatz becomes an exact ground state in the non interacting limit and in the presence of repulsive interaction it gives rise to spreading of the condensate wave function. To describe the quadrupoles and monopole oscillation, we consider the most general time dependent quadratic exponent of the variational ansatz.

We obtain the effective Lagrangian $`L`$ by substituting Eq. (12) into Eq. (10) and integrating the Lagrangian density over the space co-ordinates,

$`{\displaystyle \frac{L}{N\mathrm{}\omega _0}}`$ $`=`$ $`{\displaystyle \frac{1}{4D}}[(\beta _1\dot{\alpha _2}+\alpha _1\dot{\beta _2}2\gamma _1\dot{\gamma _2})`$ (13)
$`+`$ $`(\alpha _1+\beta _1)D+(\alpha _2^2+\gamma _2^2)\beta _1+(\beta _2^2+\gamma _2^2)\alpha _1`$ (14)
$``$ $`2(\alpha _2+\beta _2)\gamma _1\gamma _2+\lambda \beta _1+{\displaystyle \frac{\alpha _1}{\lambda }}+PD^{3/2}],`$ (15)
where $`\lambda `$ is the asymmetric ratio, $`\lambda =\omega _x/\omega _y`$ and $`P=\sqrt{\frac{2}{\pi }}2N\frac{a}{a_z}`$.

The variational energy of the static condensate at equilibrium is given in terms of the equilibrium values of the inverse square width of the condensate along x and y directions,

$$\frac{E}{N\mathrm{}\omega _0}=\frac{1}{4}[(\alpha _{10}+\beta _{10})+(\frac{\lambda }{\alpha _{10}}+\frac{1}{\lambda \beta _{10}})+P\sqrt{\alpha _{10}\beta _{10}}].$$
(16)
One can get the equilibrium value of the variational parameters, $`\alpha _{10}`$ and $`\beta _{10}`$ by minimizing the energy with respect to the variational parameters,

$$\alpha _{10}^2=\lambda \frac{P}{2}\alpha _{10}\sqrt{\alpha _{10}\beta _{10}},$$
(17)
$$\beta _{10}^2=\frac{1}{\lambda }\frac{P}{2}\beta _{10}\sqrt{\beta _{10}\alpha _{10}}.$$
(18)
From the above two relations, we obtain

$$\eta ^4+\frac{P\eta ^3}{2}\frac{P\eta }{2\lambda ^2}\frac{1}{\lambda ^2}=0,$$
(19)
where $`\eta `$ is the ratio of the condensate widths in the $`x`$ and $`y`$ direction, $`\eta =\sqrt{\frac{\beta _{10}}{\alpha _{10}}}`$. From Eq. (19) one can say how $`\eta `$ changes with the number of atoms N and the coupling constant $`g_2`$. The variation of $`\eta `$ with the the dimensionless effective interaction strength $`P`$ is shown in Fig. 1. The ratio between the widths of the condensate $`\eta `$ varies from $`1/\sqrt{\lambda }`$ to $`1/\lambda `$, as the interaction strength increases from zero to large value (Thomas-Fermi limit). In the Thomas-Fermi limit, the equilibrium values of the parameters $`\alpha _1`$ and $`\beta _1`$ are,

$$\alpha _{10}=\lambda \sqrt{\frac{2}{P}},\beta _{10}=\frac{1}{\lambda }\sqrt{\frac{2}{P}}$$
(20)
In this limit, the energy per particle is $`\frac{E}{N}=\mathrm{}\omega _0\sqrt{\frac{P}{2}}`$. In the non-interacting limit, $`\alpha _{10}^2=\lambda `$ and $`\beta _{10}^2=1/\lambda `$. The energy per particle is $`\frac{E}{N}=\frac{\mathrm{}\omega _0}{2}(\sqrt{\lambda }+\sqrt{1/\lambda })`$.

We are interested in the low-energy excitations of a Bose system. The low-energy excitations of the condensate correspond to the small oscillations of the cloud around the equilibrium configuration. Therefore we expand the time dependent variational parameters around the equilibrium points in the following way, $`\alpha _1=\alpha _{10}+\delta \alpha _1`$, $`\beta _1=\beta _{10}+\delta \beta _1`$ and $`\alpha _2=\delta \alpha _2,\beta _2=\delta \beta _2`$ $`\gamma _1=\delta \gamma _1`$, $`\gamma _{10}=0`$ and $`\gamma _2=\delta \gamma _2`$.

Using the Euler-Lagrange equation, the time evolution of the inverse square of the width around the equilibrium points are given by,

$$\delta \ddot{\alpha _1}+\lambda \frac{(8+3P\eta )}{(2+P\eta )}\delta \alpha _1+\frac{P\lambda \eta }{(2+P\eta )}\delta \beta _1=0,$$
(21)
$$\delta \ddot{\beta _1}+\frac{P}{\lambda (2\eta +P)}\delta \alpha _1+\frac{(8\eta +3P)}{\lambda (P+2\eta )}\delta \beta _1=0,$$
(22)
$$\delta \ddot{\gamma _1}+[\frac{4\lambda \eta ^2}{(2+P\eta )}+(\lambda +\frac{1}{\lambda })]\delta \gamma _1=0.$$
(23)
From the above equations we can see that the modes corresponding to the fluctuations of the average of $`x^2`$ and $`y^2`$ are coupled, but the mode associated with the fluctuation of the average value of $`xy`$ is decoupled. The Eqs. (21) and ( 22) are coupled equations of the modes $`\alpha _1`$ and $`\beta _1`$, where the mode $`xy`$ is decoupled. This is due to the fact that the Hamiltonian is invariant under reflection, $`xx`$ and $`yy`$ or $`xx`$ and $`yy`$, and the modes which are odd or even under this operation separate out.

Now we look for time dependent solutions of $`e^{i\omega t}`$ type, we obtain from Eqs. (21) and (22),

$`{\displaystyle \frac{\omega _\pm ^2}{\omega _0^2}}`$ $`=`$ $`{\displaystyle \frac{\lambda }{2}}{\displaystyle \frac{(8+3P\eta )}{(2+P\eta )}}+{\displaystyle \frac{1}{2\lambda }}{\displaystyle \frac{(8\eta +3P)}{(2\eta +P)}}`$ (24)
$`\pm `$ $`\sqrt{[{\displaystyle \frac{\lambda (8+3P\eta )}{2(2+P\eta )}}{\displaystyle \frac{(8\eta +3P)}{2\lambda (2\eta +P)}}]^2+({\displaystyle \frac{P\lambda \eta ^2}{2+P\eta }})^2}.`$ (25)
For an isotropic trap, $`\omega _+=2\omega _0`$ and $`\omega _{}^2=[\omega _0^2(8+2P)]/(2+P)`$. For large $`N`$ limit, $`\omega _{}=\sqrt{2}\omega _0`$. So $`\omega _+`$ and $`\omega _{}`$ may be identified as the monopole mode frequency and quadrupole mode frequency respectively. The monopole mode is coupled with the quadrupole mode in an anisotropic trap. However, the monopole mode frequency in an isotropic trap is independent of the interaction strength of the two-body potential and the number of particles in the condensate state. This is due to the underlying $`SO(2,1)`$ symmetry in the Hamiltonian , .

From Eq.(23) we obtain,

$$\frac{\omega _s^2}{\omega _0^2}=\frac{4\lambda \eta ^2}{(2+P\eta )}+(\lambda +\frac{1}{\lambda }).$$
(26)
In an isotropic trap, $`\omega _s`$ becomes $`\omega _{}`$. In the non-interacting limit, $`\omega _s=\omega _x+\omega _y`$. In the Thomas-Fermi limit, Eq. (26) reduces to $`\omega _s=\sqrt{\omega _x^2+\omega _y^2}`$. In an isotropic trap this mode corresponds to the quadrupole excitation. This excitation is also known as scissors mode , and this oscillation has been observed experimentally .

## III  Sum-Rules and Collective excitations

In this section, we study the quadrupole excitations of a two dimensional deformed trapped quantum gas at zero temperature within the sum rule approach. In the collisionless regime the collective excitation frequencies of a confined gas are well described by the sum rule method. The collective excitation of any system is usually probed by applying external fields. Given an excitation operator $`F`$, many useful quantities of the excited system can be calculated from the so-called strength function ,

$$S_\pm (E)=\underset{n}{}|<n|F_\pm |0>|^2\delta (EE_n),$$
(27)
where $`E_n`$ and $`|n>`$ are the excitation energy and the excited state respectively, and $`F_+=F`$, $`F_{}=F^{}`$. Various energy weighted sum rules are derived from the moments of the strength distribution function,

$$m_k^\pm =\frac{1}{2}E^k(S_+(E)\pm S_{}(E))𝑑E.$$
(28)
It is easy to see that, for a given $`k`$, the moments may be expressed in terms of the commutators of the excitation operator F with the many body Hamiltonian $`H`$. We give below some of the useful energy weighted sum rules,

$$m_0^{}=\frac{1}{2}<0|[F^{},F]|0>,$$
(29)
$$m_1^+=\frac{1}{2}<0|[F^{},[H,F]]|0>,$$
(30)
$$m_2^{}=\frac{1}{2}<0|[J^{},J]|0>,$$
(31)
$$m_3^+=\frac{1}{2}<0|[[F^{},H],[H,[H,F]]]|0>,J=[H,F],$$
(32)
where $`[,]`$ denotes the commutator between corresponding operators. Near the collective excitation frequency the strength distribution becomes sharply peaked, and the collective excitation energy is described by the moments of the strength distribution,

$$\mathrm{}\omega =\sqrt{\frac{m_3^+}{m_1^+}}.$$
(33)
Following Ref., one can derive the above form of collective excitation energy by using the variational principle. Given the many body ground state it is possible to find out the collective excitation energy and the excited state, if one is able to find an operator $`O^{}`$, which satisfies the following equation of motion :

$$[\widehat{H},O^{}]=\mathrm{}\omega _cO^{}.$$
(34)
The excitation energy is then given by the following expression,

$$\mathrm{}\omega _c=\frac{<0|[O,[\widehat{H},O^{}]]|0>}{<0|O,O^{}|0>}.$$
(35)
We may now take the variational ansatz for $`O^{}`$ as, $`O^{}=F+bJ`$ with the variational parameter $`b`$. By minimizing the energy with respect to the variational parameter, we obtain the collective excitation energy as $`E_c=\sqrt{m_3^+/m_1^+}`$, which is same as eqn.(33).

Similarly, we construct the most general excitation operator $`F=x^2+by^2`$ when monopole and quadrupole modes are coupled. $`b`$ is a variational parameter. In symmetric trap potential, if $`b=1`$, F is monopole mode and if $`b=1`$, F is the quadrupole mode. In the same way we can calculate the lowest energy excitation in this particular sector of excitations variationally. The lowest energy mode turns out to be the quadrupole mode.

Calculating the moments $`m_1`$ and $`m_3`$ by taking the excitation operator with the Hamiltonian, we obtain,

$$E_{coll}^2=\frac{4\mathrm{}^2}{m}\frac{E_x}{<x^2>}\frac{(1+Ab^2+Bb)}{(1+Cb^2)},$$
(36)
where $`A=\frac{E_y}{E_x},B=\frac{E_{int}}{2E_x},C=\frac{<y^2>}{<x^2>}`$. and $`E_x=<T_x>+<V_x>+\frac{1}{4}<E_{int}>`$, $`E_y=<T_y>+<V_y>+\frac{1}{4}<E_{int}>`$. Here $`<>`$ denotes the expectation value of the corresponding operators in the ground state and $`T_x`$, $`V_x`$ and $`T_y`$, $`V_y`$ represents the kinetic energy and potential energy along $`x`$ and $`y`$ coordinates respectively. The interaction energy is given by $`E_{int}=g_2/2|\psi |^4d^2r`$. Now we minimize this collective energy with respect to the variational parameter $`\beta `$. The value of $`\beta `$ for which the collective energy is minimum, is given by

$$b_0=\frac{2(CA)\pm \sqrt{4(CA)^2+4B^2C}}{2BC}.$$
(37)
It can be easily shown that for an isotropic trap, $`b_0=\pm 1`$. So we have identified that the monopole mode can be excited by the operator $`F=x^2+b_0y^2`$ and quadrupole mode can be generated by $`F=x^2b_0y^2`$. Inserting $`b_0`$ into Eq. (36), we obtain the following collective oscillation frequencies:

$`\omega _\pm ^2`$ $`=`$ $`{\displaystyle \frac{2}{m}}[({\displaystyle \frac{E_y}{<y^2>}}+{\displaystyle \frac{E_x}{<x^2>}})`$ (38)
$`\pm `$ $`\sqrt{({\displaystyle \frac{E_y}{<y^2>}}{\displaystyle \frac{E_x}{<x^2>}})^2+{\displaystyle \frac{E_{int}^2}{4<x^2><y^2>}}}].`$ (39)
Using the variational wave function of the ground state in deformed trap, one can easily get the excitation frequencies. The lowest energy excitation frequency in this sector is

$`{\displaystyle \frac{\omega _{}^2}{\omega _0^2}}`$ $`=`$ $`{\displaystyle \frac{\lambda }{2}}{\displaystyle \frac{(8+3P\eta )}{(2+P\eta )}}+{\displaystyle \frac{1}{2\lambda }}{\displaystyle \frac{(8\eta +3P)}{(2\eta +P)}}`$ (40)
$``$ $`\sqrt{[{\displaystyle \frac{\lambda (8+3P\eta )}{2(2+P\eta )}}{\displaystyle \frac{(8\eta +3P)}{2\lambda (2\eta +P)}}]^2+({\displaystyle \frac{P\lambda \eta ^2}{2+P\eta }})^2}.`$ (41)
It can be identified as quadrupole mode since in an isotropic trap, its excitation frequency exactly matches with the quadrupole mode frequency. The above expression for the excitation frequency Eq. (40) is exactly same as the mode frequency $`\omega _{}`$ in Eq. (24).

The higher energy excitation exactly matches within the monopole mode, although it is not the local minimum of the energy, Eq. (36). The dispersion relation for this monopole mode frequency is:

$`{\displaystyle \frac{\omega _+^2}{\omega _0^2}}`$ $`=`$ $`{\displaystyle \frac{\lambda }{2}}{\displaystyle \frac{(8+3P\eta )}{(2+P\eta )}}+{\displaystyle \frac{1}{2\lambda }}{\displaystyle \frac{(8\eta +3P)}{(2\eta +P)}}`$ (42)
$`+`$ $`\sqrt{[{\displaystyle \frac{\lambda (8+3P\eta )}{2(2+P\eta )}}{\displaystyle \frac{(8\eta +3P)}{2\lambda (2\eta +P)}}]^2+({\displaystyle \frac{P\lambda \eta ^2}{2+P\eta }})^2}.`$ (43)
For another quadrupole mode, the excitation operator is $`F=xy`$. Using the commutation relation, we obtain,

$$m_1=\frac{\mathrm{}^2}{2m}<(x^2+y^2)>,$$
(44)
$$m_3=m_3(T)+m_3(v)+m_3(ee),$$
(45)
where,

$$m_3(T)=\frac{\mathrm{}^4}{m^3}<(p_x^2+p_y^2)>,$$
(46)
$$m_3(V)=\frac{\mathrm{}^4}{2m}(\omega _x^2+\omega _y^2)<(x^2+y^2)>,$$
(47)
$`m_3(ee)`$ $`=`$ $`{\displaystyle \frac{g_2\mathrm{}^4}{2m^2}}[{\displaystyle }d^2r\rho (r)(y{\displaystyle \frac{}{x}}+x{\displaystyle \frac{}{y}})^2\rho (r)`$ (48)
$`+`$ $`{\displaystyle }d^2r(x{\displaystyle \frac{\rho }{y}}+y{\displaystyle \frac{\rho }{x}})^2].`$ (49)
Using the variational wave function, we can get all these moments. In this case $`m_3(ee)`$ exactly vanishes. So the frequency for the quadrupole mode is,

$$\frac{\omega _s^2}{\omega _0^2}=\frac{4\lambda \eta ^2}{(2+P\eta )}+(\lambda +\frac{1}{\lambda }).$$
(50)
This expression for the quadrupole frequency is also the same as Eq.(26). So $`\omega _{}`$ in Eq. (40) and $`\omega _s`$ in Eq.(50) shows the splitting occurs between two quadrupole modes in a two dimensional deformed trapped Bose gas and the dependence of the splitting on the interaction strength and trap anisotropy can be analyzed from the analytical expressions. For an isotropic trap, the two quadrupole modes are degenerate. The variation of the splitting between two quadrupole modes $`\mathrm{\Delta }_b=\omega _s\omega _{}`$, of a trapped interacting Bose gas, with the dimensionless interaction parameter $`P`$ is shown in Fig. 2.

We have checked that the sum rule method gives correct results for the excitation frequencies of the two quadrupole modes for a system of interacting bosons in an anisotropic trap. Now we apply this method to calculate the excitation energies of the quadrupole modes of a system of interacting Fermions in a deformed trap.

## IV two dimensional trapped anisotropic Fermi gas at zero temperature

In this section we discuss the collective oscillation of a two dimensional deformed trapped unpolarized Fermi gas at zero temperature within the sum-rule approach. Using this approach, the collective excitations have been studied in other finite Fermionic systems like atomic nuclei, metal clusters  and quantum dots .

We consider a two dimensional deformed trapped unpolarized Fermionic atoms at very low temperature. The two-body interaction of the dilute gas can be described by the short range pseudopotential $`V(\stackrel{}{r}\stackrel{}{r}^{})=g_2\delta ^2(\stackrel{}{r}\stackrel{}{r}^{})`$, where $`g_2`$ is the coupling constant and its form is given in Sec.II. The Hamiltonian of the trapped Fermionic atoms is given by,

$$H=\underset{i}{}\frac{p_i^2}{2m}+V_{ext}+g_2\underset{i<j}{}\delta ^2(\stackrel{}{r}_i\stackrel{}{r}_j),$$
(51)
where the confining potential is

$$V_{ext}=\frac{1}{2}m\omega _0^2(\lambda x^2+\frac{y^2}{\lambda }).$$
(52)
The Thomas-Fermi energy functional of this trapped interacting Fermi system is given by,

$$E[\rho (r)]=d^2r[\frac{\mathrm{}^2\pi }{2m}\rho ^2(r)+V_{ext}\rho (r)+\frac{\stackrel{~}{g_2}}{2}\rho ^2(r)],$$
(53)
where $`\stackrel{~}{g_2}=g_2/2`$. Here we assume the density of two spin components are same, $`\rho _{}=\rho _{}`$. The interaction energy density $`g_2\rho _{}\rho _{}`$, can be written as $`\frac{g_2}{4}\rho ^2`$, where $`\rho `$ is the total density. By minimizing the energy functional with respect to density, we obtain,

$$\rho (\stackrel{}{r})=\frac{R_F^2}{2K_0\pi a_0^4}[1\frac{r^2}{R_F^2}],rR_F,$$
(54)
where $`R_F=(4NK_0)^{1/4}a_0`$ is the radius of the atomic gas which is determined from the condition $`d^2r\rho (r)=N`$. $`K_0=1+\frac{\stackrel{~}{g_2}m}{\pi \mathrm{}^2}`$ is a dimensionless constant. At very low temperatures, collisions are suppressed due to Fermi statistics and system is in the collisionless regime. We study the collective excitation frequencies in this regime by sum rule approach.

In Sec.III we have derived the expressions for quadrupole excitation frequencies within sum-rule approach. We can use expression Eq.(38), to calculate one of the quadrupole mode frequencies for Fermi gas also.

We evaluate all the expectation values of the corresponding operators by using the Thomas-Fermi density (54),

$$\frac{E_x}{<x^2>}\pm \frac{E_y}{<y^2>}=\frac{m\omega _0^2}{2}\frac{(3K_0+1)}{2K_0}(\lambda \pm \frac{1}{\lambda }).$$
(55)
Using Eqs. (38) and (55), we obtain,

$`{\displaystyle \frac{\omega _\pm ^2}{\omega _0^2}}`$ $`=`$ $`{\displaystyle \frac{(3K_0+1)}{2K_0}}(\lambda +{\displaystyle \frac{1}{\lambda }})`$ (56)
$`\pm `$ $`\sqrt{[{\displaystyle \frac{(3K_0+1)}{2K_0}}(\lambda {\displaystyle \frac{1}{\lambda }})]^2+(1{\displaystyle \frac{1}{K_0}})^2}.`$ (57)
For an isotropic trap, the monopole mode frequency becomes $`\omega _+=2\omega _0`$.

There is another quadrupole mode for which the excitation operator is $`F=xy`$. Using the density for the trapped interacting Fermi gas at $`T=0`$, we get the following moments:

$$m_1=\frac{\mathrm{}^3R_F^6}{48m^2a_0^6\omega _0K_0}(\lambda +\frac{1}{\lambda }),$$
(58)
$$m_3=\frac{\mathrm{}^5\omega _0R_F^6}{12K_0m^2a_0^6}[\frac{1}{K_0}+\frac{(\lambda +\frac{1}{\lambda })^2}{4}].$$
(59)
In this case also, $`m_3(ee)`$ exactly vanishes. The quadrupole oscillation frequency is given by,

$$\frac{\omega _s}{\omega _0}=\sqrt{[(\lambda +\frac{1}{\lambda })+\frac{4}{K_0(\lambda +\frac{1}{\lambda })}]}.$$
(60)
For an isotropic trap, $`\omega _{}`$ in Eq. (56) and Eq. (60) becomes

$$\omega _s=\sqrt{2}\omega _0\sqrt{\frac{(1+K_0)}{K_0}}.$$
(61)
In Eq.(56) and Eq.(60), $`\omega _{}`$ exhibits the splitting of the quadrupole modes of a two dimensional deformed trapped Fermi gas. The scissors mode is also discussed in Ref.  for superfluid Fermi gas.

The monopole mode frequency of an isotropic trapped interacting Fermi system is $`2\omega _0`$ which is independent of the interaction strength of the two-body potential and the number of particles. This is because of the presence of SO(2,1) symmetry in the Hamiltonian (Eq.(51)) as discussed in the previous section.

The splitting between two quadrupole modes, $`\mathrm{\Delta }_f=\omega _s\omega _{}`$ , of a deformed trapped interacting Fermi gas is shown in Fig. 3. The frequencies of these two modes and the splitting between them are independent of the particle number for two dimensional Fermions within Thomas-Fermi approximation. This splitting decreases almost linearly with increasing interaction strength.

## V summary and conclusions

In this paper, we have mainly considered two non-degenerate quadrupole modes of a quantum gas in an anisotropic harmonic oscillator potential. We investigated the effect of interaction on the splitting between these quadrupole modes for arbitrary trap deformation. We have calculated a few low-lying collective excitation frequencies of a two dimensional trapped Bose gas in an anisotropic trap, by using time dependent variational method. We found that one quadrupole mode is coupled with the monopole mode in presence of trap deformation. Another quadrupole mode associated with the fluctuation of the average value of $`xy`$ (which is also known as scissors mode), is decoupled.

Using the energy weighted sum-rule approach we derived the general dispersion relation of the two quadrupole excitations. Using the same variational wave function for Bosons, we checked that the collective frequencies obtained from the sum-rule approach are exactly the same as those obtained from the variational method. The main advantage of the sum-rule method is that it can be applied to both trapped Bosons and Fermions to calculate the excitation frequencies in the collisionless regime. This method can be applied for any number of confined particles and also it can be generalised for long range interactions. This energy weighted sum-rule method can be extended for Coulomb interaction to study the quadrupole excitations in a deformed electronic nanostructure like an elliptic quantum dot. The splitting between the quadrupole modes obtained from this method is non perturbative in the trap anisotropy parameter.

We considered a system of two dimensional spin unpolarised interacting Fermions in an anisotropic harmonic oscillator potential within Thomas-Fermi approximation. Applying the sum-rule technique to this deformed Fermi gas, we obtain two quadrupole excitation frequencies and the splitting between them analytically. For both statistics, the amount of splitting between the quadrupole modes decreases with increasing interaction strength. For a two dimensional Fermi system the frequencies and the splitting are independent of the particle number. For an isotropic trap, the monopole mode frequency of a Bose gas as well as Fermi gas is the universal frequency $`2\omega _0`$. This monopole mode frequency is independent of the strength of the two-body interaction potential and the number of particles. This is due to the underlying $`SO(2,1)`$ symmetry in the Hamiltonian. Strictly speaking, our all the results are valid when the conditions $`\mathrm{}\omega _z\mu \mathrm{}\omega _0`$ and $`k_BT\mathrm{}\omega _z`$ are satisfied.

Recent experimental progress in MIT  on quasi two dimensional Bose condensed shows the possibilities of verification of our results. Above mentioned quadrupole modes are excited in the two dimensional plane and for simplicity we consider only the two dimensional trapped gas. This method and the most general Gaussian anstaz for the order parameter can also be extended to three dimensional anisotropic systems to study the various quadrupole modes. The splitting in these two quadrupole modes may be used to find trap anisotropy. It will be an interesting to study the splitting between the quadrupole modes of an anisotropic quantum system in presence of terms having definite chirality, like magnetic field or rotation.

We would like to thank G. Baskaran and M. V. N. Murthy for helpful discussions. LKB is a unité de recherche de l’Ecole normale supérieure et de l’Universit ’e Pierre et Marie Curie, associée au CNRS.