text:
# Local adiabatic quantum search with different paths

## 1 Introduction

One of the promising approaches to quantum computing is provided by the adiabatic evolution method introduced by Farhi et. al. and elaborated and developed by several groups . While the standard paradigm of quantum computation uses quantum gates (i.e., unitary operators) applied to quantum registers, the main idea of the adiabatic computing algorithms is to encode the initial state and solution state of the problem at hand in properly chosen Hamiltonians. Generally, the meaning of the ”proper” is the requirement that the initial and the solution states are the ground states of an initial and final Hamiltonians. Then these two Hamiltonians are combined in a time-dependent way and system is let to evolve under the influence of this combined Hamiltonian. The speed of mixing is determined by the adiabatic theorem of quantum mechanics . According to this theorem, if the evolution of a quantum system is governed by a Hamiltonian that varies slowly enough, the system will stay near its instantaneous ground state. The general time-complexity of adiabatic quantum computation is still an open problem . The complexity of the algorithm is determined by the time needed for switching from initial to solution state ”slowly enough”. The switching speed depends inversely on the size of the spectral gap between the instantaneous ground and first excited states of the Hamiltonian.

One of the successful quantum algorithms is the database search algorithm of Grover , which has been examined in great details in many studies . This algorithm shows a quadratic speedup in searching a database over classical algorithms. Adiabatic search algorithm in its local adiabatic version was shown to provide similar speedup in searching, namely $`O(\sqrt{N})`$ number of queries , where $`N`$ is the number of items in the database. Also, some numerical evidence suggests that adiabatic algorithm might efficiently solve some instances of hard combinatorial problems, outperforming classical methods .

Farhi et. al.  showed that one does not need to take a straight path from initial to final Hamiltonian. They showed that by applying a control Hamiltonain at the evolution stage an algorithmic failure can be turned to a success for 3-SAT problem. Recently, Boulatov and Smelyanskiy reported  a detailed theoretical study of using different paths for a version of 3-SAT problem.

In the present article, we show that one can get a better than square-root speedup by using different path approach of Ref. . We use part of ’sure success’ Hamiltonian of Bae and Kwon  as the control Hamiltonian. Added control Hamiltonian keeps the minimum spectral gap a constant plus a system size dependent term which goes to zero as the system size tends to get close to infinity. It is shown that complexity of the proposed algorithm approaches a constant almost independent of the size of the system.

## 2 Quantum Adiabatic Theorem

Consider a quantum system in state $`|\mathrm{\Psi }(t)`$ which evolves according to the Schrodinger equation

$`i{\displaystyle \frac{d}{dt}}|\mathrm{\Psi }(t)=H(t)|\mathrm{\Psi }(t)`$
where $`H(t)`$ is the Hamiltonian governing the dynamics of the system (we take $`\mathrm{}=1`$). If the Hamiltonian is independent of time and the system is at its ground state, it will remain there forever. According to adiabatic theorem  if the Hamiltonian varies slowly enough, the state of the system will remain at the instantaneous ground state of the Hamiltonian for all $`t`$.

Let $`|\mathrm{\Psi }_k(t)`$ be the eigenstate of $`H(t)`$, which satisfy

$`H(t)|\mathrm{\Psi }_k(t)=E_k(t)|\mathrm{\Psi }_k(t)`$
where $`E_k(t)`$ is the corresponding eigenvalue and k labels the eigenvalues ($`k=0`$ is the ground state).

The adiabacity conditions are related to the minimum energy gap of the spectrum of $`H(t)`$ and the transition to excited states of the system. These conditions can be made as follows: The minimum energy gap is defined as

$`g_{\mathrm{min}}=\underset{0tT}{\mathrm{min}}\left[E_1(t)E_0(t)\right]`$
and the maximum of the matrix element of $`dH/dt`$ between the ground and the first excited state;

$`D_{\mathrm{max}}=\underset{0tT}{\mathrm{max}}\left|\mathrm{\Psi }_1(t)\left|{\displaystyle \frac{dH}{dt}}\right|\mathrm{\Psi }_0(t)\right|`$ (1)
Adiabatic theorem states that, if we prepare the system at time $`t=0`$ in its ground state $`|\mathrm{\Psi }_0(0)`$ and let it evolve under the Hamiltonian $`H(t)`$ for a time $`T`$, then

$`\left|\mathrm{\Psi }_0(T)|\mathrm{\Psi }(T)\right|^21ϵ^2`$
provided that

$`{\displaystyle \frac{D_{max}}{g_{min}^2}}ϵ`$ (2)
where $`ϵ1`$.

## 3 Local Adiabatic Evolution Algorithm

The problem we consider here is to find a marked item in a set of $`N`$ items. Items in the system are labeled by $`n`$ qubits, so the dimension of the Hilbert space of the system is $`N=2^n`$. The initial state of the system is the uniform superposition of all basis states, which also includes the unknown marked state $`|m`$

$`|\psi _0={\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \underset{i=0}{\overset{N1}{}}}|i`$ (3)
where $`|i`$ are the basis states with $`i=0,\mathrm{},N1`$. $`|\psi _0`$ is the zero energy ground state of

$`H_0=I|\psi _0\psi _0|`$
whose ground state is $`|\psi _0`$ with energy $`0`$. Let the marked state Hamiltonian be

$`H_m=I|mm|`$
whose ground state, the marked state $`|m`$, is unknown. The time-dependent Hamiltonian that interpolates between these two Hamiltonians is

$`H(s)=(1s)H_0+sH_m`$ (4)
where the mixing parameter $`s=s(t)`$ is a monotonic function with $`s(0)=0`$ and $`s(T)=1`$, $`T`$ being the required running time of the algorithm. The algorithm is to prepare the system in the state $`|\psi (0)=|\psi _0`$ and then apply $`H(s)`$ for a time interval $`T`$ to end up with the state $`|m`$. Roland and Cerf  shown that complexity of this algorithm when applied in local adiabatic scheme scales with $`\sqrt{N}`$, similar to Grover’s algorithm and this result is optimal for any chose of the time evolution of $`s(t)`$.

One can generalize the time-dependent Hamiltonian of Eq. 4 by adding to it a Hamiltonian of the form

$`H_D=a(s)|\psi _0m|+b(s)|m\psi _0|`$
where we have chosen $`a(s)=b(s)=s(1s)`$ so that the Hamiltonian at $`s=0`$ is $`H_0`$ and at $`s=1`$ it is $`H_m`$. $`H_D`$ acts only during the adiabatic evolution of the system. So its effect is to change the path taken by the system for going from the initial state $`|\psi _0`$ to final state $`|m`$. A similar extra Hamiltonian was shown to convert an algorithmic failure to a success by Farhi et.al. . Now the time-dependent Hamiltonian of the system is

$`H(s)=(1s)H_0+sH_mH_D`$ (5)
which in matrix representation can be written as

$`H=E\left(\begin{array}{cc}(s1)(x(s+x)1)& (s1)(s+x)\sqrt{1x^2}\\ (s1)(s+x)\sqrt{1x^2}& s(1s)x^2\end{array}\right)`$
where $`x=m|\psi _0=1/\sqrt{N}`$. Two lowest eigenvalues of this Hamiltonian are plotted as function of $`s`$ in Fig. 1. These eigenvalues are separated by a gap

$`g(s)=\sqrt{14s(1s)\left(1x(x+1)s(1s)\right)}`$ (7)
note that the minimum of $`g(s)`$ is at $`s=1/2`$ with value $`g_{\mathrm{min}}=1/2+1/\sqrt{N}`$.

The eigenvectors corresponding to two lowest eigenvalues are

$`|E_0;s`$ $`=`$ $`{\displaystyle \frac{(1g(s)2s)|m>+\left(s1\right)\left(s+x\right)\left(1+\sqrt{1x^2}\right)|\psi _0>}{\sqrt{4\left(s1\right)^2\left(s+x\right)^2\left(1x^2\right)+\left(1g(s)2s+2\left(s1\right)x\left(s+x\right)\right)^2}}}`$
$`|E_1;s`$ $`=`$ $`{\displaystyle \frac{(1+g(s)2s)|m>+\left(s1\right)\left(s+x\right)\left(1+\sqrt{1x^2}\right)|\psi _0>}{\sqrt{4\left(s1\right)^2\left(s+x\right)^2\left(1x^2\right)+\left(1+g(s)2s+2\left(s1\right)x\left(s+x\right)\right)^2}}}`$
The matrix element in Eq. 1 can be written as

$`{\displaystyle \frac{dH}{dt}}_{1,0}`$ $`=`$ $`{\displaystyle \frac{ds}{dt}}{\displaystyle \frac{d\stackrel{~}{H}}{ds}}_{1,0}`$ (8)
$`=`$ $`{\displaystyle \frac{ds}{dt}}(H_mH_0+(12s)H_D_{1,0})`$
$`H_mH_0+(12s)H_D_{1,0}1`$ .

In the the local adiabatic approximation one tries to adjust the evolution rate of the $`s`$ for infinitesimal time intervals $`dt`$ and apply the adiabacity condition locally to each of these intervals. The adiabacity condition for the local adiabatic evolution would be

$`\left|{\displaystyle \frac{ds}{dt}}\right|ϵ{\displaystyle \frac{g^2(s)}{\left|\frac{d\stackrel{~}{H}}{ds}_{1,0}\right|}}`$
for all times $`t`$. Using Eq. 7, the evolution of the Hamiltonian is chosen to evolve at a rate that is solution of

$`{\displaystyle \frac{ds}{dt}}=ϵg^2(s)=ϵ\left(14s(1s)\left(1{\displaystyle \frac{1}{N}}{\displaystyle \frac{1}{\sqrt{N}}}+s(s1)\right)\right)`$ (9)
Integral of this equation involves the logs of the roots of $`g^2(s)`$ and is displayed in Eq. 10. $`s(t)`$ obtained from this integral is shown in Fig. 2 along with the $`s(t)`$ of Roland and Cerf’s algorithm.

$`t={\displaystyle \frac{1}{\sqrt{2}h(N)}}\left({\displaystyle \frac{1}{k_2(N)}}\mathrm{tanh}^1\left[{\displaystyle \frac{2\sqrt{2}k_2(N)s}{2k_2(N)^24s}}\right]{\displaystyle \frac{1}{k_1(N)}}\mathrm{tanh}^1\left[{\displaystyle \frac{2\sqrt{2}k_1(N)s}{2k_1(N)^24s}}\right]\right)`$ (10)
where

$`h(N)`$ $`=`$ $`{\displaystyle \frac{1}{N}}\sqrt{(1N)(2\sqrt{N}+1)}`$
$`k_1(N)`$ $`=`$ $`\sqrt{24h(N)+{\displaystyle \frac{4}{N}}+{\displaystyle \frac{4}{\sqrt{N}}}}`$
$`k_2(N)`$ $`=`$ $`\sqrt{2+4h(N)+{\displaystyle \frac{4}{N}}+{\displaystyle \frac{4}{\sqrt{N}}}}`$
As can be seen from the figure, added term $`H_D`$ changes $`H(s)`$ much faster compared to the original Hamiltonian of Ref. . The computation time of the algorithm is evaluated by taking $`s=1`$ in the solution of Eq. 9, and is plotted against $`N`$ in Fig. 3 along with the same quantity from Ref. . The asymptotic form of Eq. 9 for $`s=1`$ is given by

$`T=\left(1+{\displaystyle \frac{\pi }{4}}\right)+O\left(1/\sqrt{N}\right)`$ (11)
which show that as the system size $`N`$ approaches infinity, the required running time approaches the constant $`\left(1+\frac{\pi }{4}\right)`$.

The main factor that determines the required evolution time $`T`$ is the behavior of minimum gap $`g_{\mathrm{min}}`$. This energy difference usually reaches its minimum at an avoided level-crossing. If this approach is exponentially fast with a growing particle number then the algorithm is not efficient. The effect of extra term $`H_D`$ is to keep these two levels far apart. Our result seem to be similar to one obtained by Das et. al.  who use Eq. 4 with different coefficients for $`H_0`$ and $`H_m`$ terms. One can, actually have a minimum gap that is larger than Eq. 7 by choosing $`a(s)=b(s)=\sqrt{s(s+1)}`$ in the definition of the driving Hamiltonian $`H_D`$. In that case, the running time approaches $`1`$ for the large $`N`$ limit.

## 4 Conclusion

In summary, we have presented the results of an analysis of a generalized adiabatic quantum search algorithm. Our analysis shows that the optimal speed associated with Grover’s search algorithm can be improved by a careful choice of the time-dependant Hamiltonian. This analysis does not deal with the complexity of implementing the driving Hamiltonian. Hamiltonian of adiabatic algorithms includes searched state $`|m`$ terms. The construction of these terms are thought to be analogous to having access to an oracle in original Grover algorithm. It seems that the time-complexity of searching in present approach is shifted from searching process to construction of the Hamiltonian.


format: text

text:
# Staggered flux vortices and the superconducting transition in the layered cuprates

(September 19, 2003)

## Abstract

We propose an effective model for the superconducting transition in the high-$`T_c`$ cuprates motivated by the SU(2) gauge theory approach. In addition to variations of the superconducting phase we allow for local admixture of staggered flux order. This leads to an unbinding transition of vortices with staggered flux core that are energetically preferable to conventional vortices. Based on parameter estimates for the two-dimensional $`t`$-$`J`$-model we argue that the staggered flux vortices provide a way to understand a phase with a moderate density of mobile vortices over a large temperature range above $`T_c`$ that yet exhibits otherwise normal transport properties. This picture is consistent with the large Nernst signal observed in this region.

The nature of the pseudogap phase of the (hole-)underdoped high-$`T_c`$ cuprates is one of the central questions of correlated electron physics. A number of scenarios and descriptions are successful in capturing certain aspects of the problem. Yet many theories face substantial difficulties when it comes to combining the large number of experimentally established anomalies of the underdoped state. A new challenge in this context has been set out recently by the Nernst effect measurements in underdoped sampleswang . In these experiments a thermal gradient is applied in the copper-oxide planes. In the presence of a small out-of-plane magnetic field, a voltage drop in the direction perpendicular to magnetic field and thermal gradient is observed. This voltage is interpreted as the phase slip signal arising from vortices moving from hot to cold. Thus the Nernst effect experiment reveals the existence of substantial superconducting (SC) short range correlations over a sizable region, starting significantly below the onset temperature $`T^{}`$ of the pseudogap, but extending up to temperatures high above the low $`T_c`$s of underdoped samples.

At first sight the observation of vortices above $`T_c`$ fits well into a scenarioemkiv  where SC phase fluctuations destroy the long range coherence, and short range pairing correlations survive up to much higher temperatures. In view of the small superfluid weight of underdoped systems it is conceivable that the SC transition is driven by a vortex unbinding similar to the XY-transition, the critical temperature disappearing like the doping $`x`$ for $`x0`$. However a simple phase fluctuation scenario faces the following problem. The creation of a vortex comes at a price, as the SC order parameter goes to zero in the vortex core and condensation energy is lost. In disordered films of conventional superconductors, where a Berezinskii-Kosterlitz-Thouless(BKT) transition can be observed, the mean-field critical temperature and the vortex unbinding temperature are similar, and due to the mean free path $`\mathrm{}`$ entering the effective coherence length $`\xi =\sqrt{\xi _0\mathrm{}}`$ (with $`\xi _0v_F/\mathrm{\Delta }`$) the vortex core energy $`\xi ^2\mathrm{\Delta }^2/\epsilon _F`$ becomes small near $`T_c`$, as $`\mathrm{\Delta }`$ decreases faster for $`TT_c`$ than $`\xi `$ increases. Thus the loss of condensation energy is limited and the vortices can be cheap. In contrast with that the high-$`T_c`$ cuprates are in the clean limit and for underdoped samples the gap magnitude remains large up to temperatures far above the SC transition. Hence the vortex core energy for bringing the gap magnitude down to zero inside the vortex would normally be expected to be huge ($`\epsilon _F`$ in BCS theory and of order $`J`$ in our case), and only exponentially few of these expensive vortices could be created above the small $`T_c`$ in underdoped samples. Then we would expect that the transport properties of the pseudogap state resemble those of a flux-flow (FF) phase. Transport resembling the normal state only occurs at high temperatures of order of the core energy when the vortices proliferate and overlap. However this picture is inconsistent with experiments. These show that above a limited fluctuation regime close to $`T_c`$ the in-plane transport looks rather normal and signs of FF conductivity do not extend far above $`T_c`$ – contrary to the Nernst signal. Apparently the conductivity $`\sigma _n`$ due to a significant number of quasiparticles dominates over the FF conductivity $`\sigma _{\mathrm{FF}}`$ in the total conductivity $`\sigma =\sigma _n+\sigma _{\mathrm{FF}}`$. Thus one has to explain two things: where the normal excitations come from and why the FF contribution is small. $`\sigma _{\mathrm{FF}}`$ can be estimated to be $`\eta /n_V`$, where $`n_V`$ is the density of vortices either forced in by a magnetic field in the mixed state or generated thermally above the BKT transition. $`\eta `$ is the friction coefficient for the vortex motion. The FF conductivity is small if $`\eta `$ is small and $`n_V`$ is not. In the following we show that our model produces a moderately large $`n_V`$ even for underdoped samples with low $`T_c`$. We will also present an approximate calculation that yields a finite density of normal excitations. We will not attempt to calculate $`\eta `$ in this work. An effect that may reduce $`\eta `$ is the experimentally observedvortexcore  small low-energy density of states in the vortex cores. This translates into small dissipation due to vortex motion and thus small $`\sigma _{\mathrm{FF}}`$. Note that $`\eta `$ cannot become arbitrarily small if we want $`\sigma _n`$ to dominate the conductivity, as the the quasiparticles responsible for $`\sigma _n`$ will cause some dissipation. Thus a small FF contribution requires a moderately large $`n_V`$ at temperatures above the limited fluctuation regime near $`T_c`$.

An extreme way to obtain normal transport properties just above $`T_c`$ is to make the vortices very dense, such that they overlap just above the transition. Then it is hard to understand why the vortex Nernst signal persists to temperatures so high above $`T_c`$. A theory considering purely Gaussian SC fluctuationsussishkin  gave good agreement for the Nernst effect in overdoped and optimally doped samples. However the description of underdoped samples becomes problematic and an additional suppression of $`T_c`$ had to be invoked.

Thus neither very few and expensive nor too many and too cheap vortices seem to match the experimental picture. Basically what is needed is a theory which produces a core energy of order $`T_c`$ rather than $`J`$. In other words we need to decrease the core energy by placing in the vortex core another non-superconducting state that is nearby in energy. There are several proposals in the literaturearovas  for such a cheap vortex core, mainly emphasizing the vicinity of the $`d`$-wave superconductor to other ordered states. Here we study the possibility of a staggered flux (SF) state inside the vortexsu2v ; dhlee . This scenario has the advantage that it emerges naturally from the SU(2) invariance, i.e. the Mott insulating nature of the undoped state. Note that the vicinity to the Mott state is also responsible for the small superfluid stiffness $`\rho _sx`$.

The idea that vortices in the underdoped system have SF cores ties in with a more general picture of the pseudogap state. This is derived from the SU(2) gauge theory for the $`t`$-$`J`$ model and views the pseudogap regime as a thermally disordered state, where the system fluctuates between various types of short range order corresponding to mean-field states that would all become identical at zero doping. The two most prominent correlations are $`d`$-wave superconductivity - determining the ground state as soon as the other fluctuations freeze out at low $`T`$ – and SF correlations. The latter represent - in addition to conventional phase fluctuations of the SC order parameter - the lowest lying fluctuations around the $`d`$-wave SC state with the largest spectral weightcollmodes . We have arguedhonlee  that the scattering of quasiparticles with these SF fluctuations may be related to the partial loss of the quasiparticle peaks in the pseudogap state.

Some of the ideas presented here carry over to other types of cheap vortices. Indeed there are experimental indicationskaku  for antiferromagnetic ordering at low temperatures in the vortex cores of optimally doped Tl-compounds, and it is quite likely that the vortex core may contain both SF and AF correlationsogata . The focus on SC and SF correlations is an attempt to concentrate on the main tendencies suggested from the SU(2) approach. The SF state has the advantage that it has a gap structure identical to the $`d`$-wave superconductor and naturally explains the gap in the vortex core. Other types of correlations, such as antiferromagnetic tendencies, may be viewed as additional instabilities, which naturally coexist with the SF order for small doping.

Recently Ivanov and Leedima  calculated the energy differences between SF and $`d`$-wave SC states using the Gutzwiller projection technique. Together with the computed superfluid stiffness this provides an estimate for the energy of a SF vortex. Note that in Ref. dima  the vortex core turns out to very small, but on the other hand the core state, taken to be a pure SF state, is likely to be somewhat too high in energy, and a better core state will increase the size of the core again. Keeping in mind these uncertainties and for the lack of better parameters we will use the numbers of Ref. dima  as input for a generalized XY-model.

Let us begin with the SU(2) mean-field (MF) theorycollmodes . Here we are interested in low temperatures. Hence we assume that the bosons carrying the electronic charge are condensed. The Hamiltonian for the fermionic spin degrees of freedom $`f_i`$, $`f_i`$ on the lattice sites $`i`$ is given by

$$H_f=\frac{J}{2}\underset{ij}{}\left(\begin{array}{c}f_i\\ f_i^{}\end{array}\right)^{}\left(\begin{array}{cc}\chi +W_{ij}& \mathrm{\Delta }_{ij}\\ \mathrm{\Delta }_{ij}^{}& \chi +W_{ij}\end{array}\right)\left(\begin{array}{c}f_j\\ f_j^{}\end{array}\right).$$
(1)
$`\chi `$ contains the hopping contributions and the constraint field $`a_{0,3}`$, which is a function of the doping. Next we allow for local admixtures of SF amplitude in exchange for SC pairing amplitude, described by an angle $`\theta _i`$collmodes , and fluctuations of the SC phase, $`\alpha _i`$. Then the SF amplitude on the bond $`ij`$ is given by $`W_{ij}=i\mathrm{\Delta }_0(1)^{i_x+j_y}\mathrm{cos}\frac{\theta _i+\theta _j}{2}`$ and the pairing amplitude is $`\mathrm{\Delta }_{ij}=(1)^{i_y+j_y}\mathrm{\Delta }_0\mathrm{sin}\frac{\theta _i+\theta _j}{2}\mathrm{exp}\left[i(\alpha _i+\alpha _j)/2\right]`$. The pure superconductor has $`\theta =\pi /2`$, while the two degenerate SF states have $`\theta =0`$ and $`\theta =\pi `$.

Now consider an effective Hamiltonian for $`\alpha _i`$ and $`\theta _i`$,

$`H`$ $`=`$ $`{\displaystyle \underset{ij}{}}\rho _s(x,\theta _i,\theta _j)\mathrm{cos}(\alpha _i\alpha _j)+{\displaystyle \underset{i\text{VP}}{}}m_\theta \mathrm{cos}^2\theta _i`$ (2)
$`+K{\displaystyle \underset{ij}{}}\left(\theta _i\theta _j\right)^2+{\displaystyle \underset{\text{VP}}{}}H_V.`$
The first term is the phase stiffness of the SC phase that depends on doping $`x`$ and the local $`\theta `$. The second term is a mass term for $`\theta `$ that takes into account the energy difference between SF and SC state outside the vortex plaquettes (VP). The third term is a gradient term for the $`\theta `$-variation. The last term is the vortex core energy, $`H_V=_{i\text{VP}}(m_n\mathrm{sin}^6\overline{\theta }+m_\theta )`$. The sum is over the four sites on each of the VP. Here $`\mathrm{\Delta }_{ij}`$ is assumed to vanish and according to Eq. 1, $`\theta `$ now describes an interpolation between the SF state ($`\theta =0`$,$`\pi `$) and the zero flux state ($`\theta =\pi /2`$) which is a Fermi liquid. As discussed earlier, the Fermi liquid core is expected to be costly and the energy costs of the SF state and the Fermi liquid state are $`m_\theta `$ and $`m_\theta +m_n`$, respectively. The specific $`\theta `$-dependence chosen for $`H_V`$ comes from the MF theory for uniform $`\theta `$ and $`\overline{\theta }`$ is the average $`\theta _i`$ over the four sites of the VP. Requiring that the vortex core area $`4a^2`$ equals the vortex size $`\pi \xi ^2`$ of a pure SF vortex in the microscopic calculation in Ref. dima , we obtain that the lattice constant of Eq. 2 is slightly larger than the one of the underlying $`t`$-$`J`$ model with a scale factor of $`1.2`$. We assume that the variation of the vortex core size can be neglected in the small doping and low temperature region we are interested in.

In principle, all parameters for (2) can be obtained from the MF theory by integrating out the fermions as in Ref.collmodes , but where possible we use parameters obtained by the Gutzwiller variational treatment of Ivanov and Leedima . With $`t/J=3`$, this gives $`\rho _s(x,\theta =\pi /2)0.75xJ`$, $`m_\theta 0.33xJ`$ and $`m_n(0.25x)J/2`$. The $`\theta `$-dependence of $`\rho _s`$ is obtained numerically from MF theory which shows a rapid linear increase from zero when $`\theta `$ deviates from $`0`$ or $`\pi `$. The coefficient $`K`$ for the gradient of $`\theta `$ is difficult to extract from the collective mode spectrumcollmodes , as the $`\theta `$-mode energy is not simply quadratic around $`(\pi ,\pi )`$. Nevertheless the $`q`$-dependence of the $`\theta `$-variations is quite weak and less than that of the SC phase $`\alpha `$ and goes to zero for $`x0`$. Thus we estimate $`K\rho _s/2`$. Other choices give similar results.

In the scaling theory for the BKT transitionchailub  the vortex core energy is assumed to be large compared to the temperature such that the fugacity $`y=\mathrm{exp}(E_c/T)`$ can be used as small parameter. In the limit $`y0`$ the transition occurs at $`T_c=\pi \rho _s/2`$. A nonzero fugacity $`y>0`$ leads to a reduction of $`T_c`$ from this upper bound by roughly $`\pi ^2\rho _sy`$. With the parameters above the core energy of a single ideal SF vortex is $`E_c=0.75\pi xJ`$ or twice the maximal $`T_c^{max}=\pi \rho _s/2`$. This leads to a reduction of the critical temperature down to $`T_c1.06\rho _s`$.

The model described by Eq. 2 can be simulated with Monte Carlo methods. To estimate $`T_c`$ we calculate the helicity modulus $`\mathrm{{\rm Y}}`$ whichteitel  measures the rigidity with respect to a phase gradient in the system. $`\mathrm{{\rm Y}}`$ vanishes above the SC transition: in BKT theory it jumps to zero at $`T_c`$, the height of the jump being $`2T_c/\pi `$, independent of the core energynelson . We also measure the average $`\theta `$-variation inside and outside the vortex cores and how the number of vortices depends on doping and temperature.

A snapshot from the Monte Carlo is shown in Fig. 1 for a sample above $`T_c`$. We observe two relatively well isolated and other less separated vortices. The SC phase $`\alpha _i`$ is disordered but exhibits remnants of short range order. For $`T=0.09J`$ and $`x=0.06`$, $`\theta `$ varies rapidly in space due to its light mass $`m_\theta <T`$ and small gradient terms. Hence there is a significant amount of SF admixture reducing the SC pairing amplitude locally even outside the vortex cores (see shaded areas in Fig. 1). Outside the vortices, the average SF amplitude is $`|\mathrm{cos}\theta |0.3`$ (compared to 1 for the pure SF state). Inside the vortex cores it is strongly enhanced, but not maximal ($`|\mathrm{cos}\theta |0.7`$ for the sample shown, see also left plot in Fig. 3).

In the left plot of Fig. 2 we show the helicity modulus $`\mathrm{{\rm Y}}`$ for two dopings $`x`$ and a 100$`\times `$100 system. $`\mathrm{{\rm Y}}`$ goes to zero above a doping-dependent temperature, but does not exhibit the universal jump of the XY-modelnelson ; teitel . This is clearly a finite-size effect which is exacerbated by the small vortex density and which does not occur in the XY-model where the core energy is zero. An estimate for the true $`T_c`$ is the intersection of the data with the line $`2T/\pi `$, that is based on the the jump criterion $`\left(\mathrm{\Delta }\mathrm{{\rm Y}}\right)_{T_c}=\pi /2`$nelson . With that we arrive at $`T_c^e0.75xJ\rho _s`$. A numerical bound is $`T_cxJ`$. Note that in a $`d`$-wave superconductor thermally excited nodal quasiparticles lead to an additional reduction of the superfluid density that is not included here. We expect however that this does not affect the nature of the SC transition.

The temperature and doping dependence of the vortex density is summarized in the right plot of Fig. 2. The onset temperature for a finite vortex density is an increasing function of $`x`$, approximately given by $`T_V\rho _s`$. In the underdoped system there is wider temperature range, starting at the BKT $`T_c`$ and extending up to the mean-field $`T_c0.18J`$ where the vortex density continues to increase and does not saturate. Thus the vortices do not overlap and there is the possibility that SC phase coherence is still well defined locally in a range above $`T_c`$. For $`x=0.06`$ the phase correlation length $`\xi _\alpha `$ is $`7`$ lattice spacings $`a`$ at $`T2T_c^e`$. This is in contrast with the normal XY-model where at $`2T_c`$, $`\xi _\alpha a`$, and on average there is a vortex on every 5th site. We emphasize that in our model the phase fluctuations are not the only source of disorder. This can be clearly seen in the right panel of Fig. 1. The SF fluctuations are not limited to the vortex cores and lead to sizable amplitude fluctuations as well. Thus the proximity to the Mott state gives rise to the low energy scale for SF fluctuations and is responsible for both phase and amplitude fluctuations proliferating at comparable temperature scales.

We now address the electronic excitations. At finite doping the bulk SF state has small Fermi pockets, and correspondingly it is natural to expect that a state which fluctuates between a superconductor with gap nodes and a SF state with small Fermi arcs will produce a finite density of states (DOS) at low energies. When we restrict the considerations to static configurations, we can calculate the quasiparticle spectrum from Hamiltonian (1) for a given configuration of fluctuations on a finite system and average over many samples.

Results are shown in the right plot of Fig. 3. The DOS exhibits a suppression for all temperatures below the mean-field transition at $`T0.18J`$, but the gap fills in when $`T`$ is increased through the SC transition. The local DOS is inhomogeneous, but is not simply correlated with or confined to the positions of vortices or regions of higher SF amplitude. These results share some aspects with a recent work by Eckl et al.eckl  who considered disordering the $`d`$-wave superconductor by phase fluctuations modeled by a simple XY model. In their model the vortex cores are conventional and the core energy is zero by construction. Our generalized model accomodates both cheap vortices with SF core and energetically more expensive vortices with Fermi liquid core. Our calculation shows that already at the BKT $`T_c`$ there is a finite number of quasiparticle excitations at low energies, and it is likely that the conductivity will be dominated by these normal excitations. The quasiparticle contribution to the Nernst signal however is smallwang  and the vortex contribution will dominate the Nernst signal at low $`T`$.

In conclusion, we have presented a simple model, motivated by the SU(2) approach for the high-$`T_c`$ cuprates, that describes the superconducting transition as an unbinding transition of vortices with staggered flux core. Using parameter estimates from projected wave functions and the SU(2) mean-field theory it allows us to understand the occurrence of a moderate vortex density even for underdoped systems with low $`T_c`$. The vortex density is determined by an energy scale that is closely related to the energy difference between the SF and the $`d`$-wave SC state and that disappears towards zero doping, making the vortices relatively cheap. However there is wider temperature range above $`T_c`$ where the vortices are sufficiently dilute and do not overlap such that we expect that the superconducting phase coherence stays intact locally. This is in accordance with the interpretations of the Nernst effect measurementswang . Phase and SF (i.e. amplitude) fluctuations lead to a filling in of the gap in the density of states at small energies already near the superconducting $`T_c`$. This could account for the normal-looking in-plane transport of the pseudogap phase which onsets just above the superconducting transition. A full analysis of this issue requires a calculation of the flux-flow resistivity in the fluctuating state. Although our approach involves rather strong simplifications and approximations we believe that it describes a way to understand the simultaneous occurrence of normal transport behavior (e.g. in the resistivity) and strong SC fluctuations, as witnessed by the vortex Nernst effect.

We thank D.A. Ivanov, O. Motrunich and A. Vishwanath for useful discussions. CH acknowledges support by the German research foundation (DFG) and PAL by NSF grant DMR-0201069.


text:
# Major Mergers of Galaxy Haloes: Cuspy or Cored Inner Density Profile?

## 1 Introduction

In the currently favored cosmological constant plus cold dark matter ($`\mathrm{\Lambda }`$CDM) model of cosmology, structure forms hierarchically: small dark matter haloes collapse at high redshifts, while larger and more massive objects form through a series of minor mergers that accrete smaller mass haloes and major mergers with comparable mass haloes. Since dark matter haloes provide the gravitational potential wells necessary for galaxy formation (White & Rees, 1978), a detailed understanding of the effects of mergers is necessary in order to test models of galaxy formation and evolution within the $`\mathrm{\Lambda }`$CDM paradigm.

Recent numerical studies indicate that the hierarchical build-up process results in haloes with significant substructure populations: hundreds to thousands of subhaloes to the smallest resolvable mass scales appear to reside within galaxy sized haloes, comprising approximately 10% of a halo’s mass (Ghigna et al., 1998; Klypin et al., 1999; Moore et al., 1999b). A number of recent papers have investigated the evolution and detailed properties of the substructure using both semi-analytic models and numerical simulations. Depending on factors such as the relative mass, pericentric distance, and halo concentration, the orbits and masses of the subhaloes can be strongly affected by dynamical friction and tidal stripping (e.g. Velázquez & White 1999; Taylor & Babul 2001; Hayashi et al. 2003; Taffoni et al. 2003; Benson et al. 2004). Dense subhaloes that survive the tidal process can sink to the centre of the parent halo, impacting its inner structure (e.g. Dekel et al. 2003; Dekel, Devor, & Hetzroni 2003).

In contrast to minor mergers, where large haloes swallow up smaller ones in a relatively gentle fashion, major mergers are more violent events involving systems of approximately equal mass that can be studied in detail with numerical simulations only. The interest in major mergers has a long history. Toomre & Toomre (1972) noted that many observed features of galaxy pairs can be explained by tidal interactions that could promote merging. They further speculated that if such interactions were common enough, elliptical galaxies might naturally be explained as the merger remnants of spiral galaxies. This suggestion has led to numerous studies of galaxy merging and has been quite successful in helping explain the origin of ellipticals. In some of the earliest simulations, White (1978, 1979) studied several types of mergers using 250 simulation particles per halo and concluded that products of mergers of cored galaxies were significantly more centrally concentrated than their progenitors, with more extended haloes having an approximate power-law density profiles $`\rho r^3`$ that resembled the de Vaucouleurs’ empirical profile for the surface-brightness of ellipticals (de Vaucouleurs, 1948). The higher resolution simulations by Pearce, Thomas, & Couchman (1993) with $`8000`$ particles per halo also exhibited the $`r^3`$ behavior with evidence for homology. Simulations with multiple components representing dark matter haloes, stellar disks, and stellar bulges have been a source of much theoretical progress (e.g. Barnes & Hernquist 1992 and references therein) and are essential for elucidating the dynamical coupling between the dark matter and stellar components in galaxies and for comparing simulation results with observations. Pure dark halo simulations, however, are required for a controlled study of the non-baryonic dark matter that dominates the gravitational potential of seemingly all galactic systems.

In this paper we perform high-resolution $`N`$-body simulations of mergers of two equal-mass dark matter haloes and study how the phase-space distribution of dark matter particles is affected by such violent events and how the merger product relaxes into equilibrium. We pay particular attention to the central regions of the haloes and examine whether the initial inner structure of a dark matter halo, either with a flat core or a sharp cusp, is disrupted or preserved by major collisions. Our study complements the numerous recent cosmological $`N`$-body studies of the properties and evolution of dark matter haloes, most of which have reached the somewhat surprising conclusion that over a wide range of mass scales, these haloes are universally cuspy and have radial density profiles with $`\rho (r)r^\gamma `$ (with $`1.0<\gamma <1.5`$) in the central region and $`\rho r^3`$ near the virial radius (Navarro, Frenk, & White, 1997; Fukushige & Makino, 1997; Moore et al., 1999a; Ghigna et al., 2000; Fukushige, Kawai, & Makino, 2004). The seemingly generic prediction of an inner cusp agrees with earlier simulations (e.g. Dubinski & Carlberg 1991) but is contrary to the inferred dark matter distribution in some dwarf and low surface brightness galaxies (e.g. McGaugh & de Blok 1998; Simon et al. 2003, though also see Swaters et al. 2003) as well as in spiral galaxies (Salucci, 2001) and has lead to speculation that the $`\mathrm{\Lambda }`$CDM model is in crisis.

The simulations of major mergers performed in this paper provide a set of controlled experiments for us to address questions such as: How does the phase space density of dark matter haloes evolve and re-equilibrate during major mergers? Is the cuspy profile an attractor region of the phase space so that merger remnants of two non-cuspy (i.e. cored) haloes will end up cuspy? Do cuspy haloes, once formed (by any process), retain their inner structure after inevitable major mergers with other haloes? Unlike cosmological simulations that are designed to incorporate many complicated processes at once, the simulations performed here allow us to separate out the effects of minor mergers and to focus on how major collisions change the spatial and velocity structures of any initial haloes of our choice. Additionally, cosmological parameters are unimportant for the simulations themselves since major mergers are local processes of two fully collapsed objects. In this case, all we need is the initial phase space distribution of each galaxy halo and the initial relative positions and velocities of the haloes themselves. By simulating collisions of two galaxy haloes with high spatial resolution and up to 5 million particles (for each halo), we are able to follow their phase-space evolution and the density profiles of the remnants down to $`0.004`$ $`r_{200}`$ in our highest resolution simulation.

Section 2 of this paper describes the initial setup and the parameters used in the numerical simulations. We also discuss in some detail the tests performed to quantify the numerical effects of two-body relaxation on the central regions of the haloes and the technical issue of how to determine the centre of a merged halo. Section 3 contains our results from the ten production runs for the energy, density profiles, velocity distributions, and phase space structure of merger remnants. Section 4 discusses the energy distribution and the relaxation of a merged halo to an equilibrium state, with an emphasis on violent relaxation and phase mixing. Section 5 contains a summary of our results.

## 2 Procedure

### 2.1 Initial Conditions

We study two general types of inner density profiles for galaxy haloes: cored and cuspy. For the cored profile, we use an isothermal sphere with a core radius and an exponential cut-off (Hernquist, 1993):

$$\rho (r)=\frac{1}{2\pi ^{3/2}}\frac{\alpha M}{r_{\mathrm{cut}}}\frac{\mathrm{exp}(r^2/r_{\mathrm{cut}}^2)}{r^2+r_{\mathrm{core}}^2}$$
(1)
where $`M`$ is the halo’s mass, $`\alpha `$ is a normalization constant given in terms of $`xr_{\mathrm{cut}}/r_{\mathrm{core}}`$ as $`\alpha ^1=1\sqrt{\pi }xe^{x^2}\mathrm{erfc}(x)`$, $`r_{\mathrm{cut}}`$ is the outer exponential fall-off radius and $`r_{\mathrm{core}}`$ is the core radius. For the cuspy haloes, we use the Navarro, Frenk, and White (NFW) profile (Navarro, Frenk, & White, 1997):

$$\rho (r)=\frac{\rho _c\overline{\delta }}{(r/r_s)(1+r/r_s)^2}$$
(2)
with scale radius $`r_s`$ and characteristic density $`\overline{\delta }`$ given by

$`r_s={\displaystyle \frac{r_{200}}{c}}={\displaystyle \frac{1}{c}}\left({\displaystyle \frac{3M_{200}}{800\pi \rho _c}}\right)^{1/3}`$ (3)
$`\overline{\delta }={\displaystyle \frac{200c^3}{3[\mathrm{ln}(1+c)c/(1+c)]}}`$ (4)
where $`c`$ is the halo concentration parameter and $`r_{200}`$ is the radius at which the average interior density is 200 times the critical density of the universe $`\rho _c=3H^2/8\pi G`$. (We also adopt $`r_{200}`$ as a halo’s virial radius and therefore $`M_{200}800\pi \rho _cr_{200}^3/3`$ as its virial mass.)

We generate our initial conditions using the methods described by Hernquist (1993). The initial positions of the particles in a halo are drawn randomly from its density profile. All of our initial haloes are chosen to be $`10^{12}M_{}`$. In the case of the NFW profile the total mass is formally infinite, so it is necessary to introduce a truncation radius for the halo. We choose to truncate our haloes at $`r=r_{200}`$. For our choice of $`c=10`$, the virial and truncation radius is 162.6 $`h^1`$ kpc and the scale radius is $`r_s=r_{200}/c=16.26`$ $`h^1`$ kpc. For the cored profile, we choose $`r_{\mathrm{core}}=8`$ $`h^1`$ kpc, $`r_{\mathrm{cut}}=81`$ $`h^1`$ kpc, and truncate the halo at $`2r_{\mathrm{cut}}`$. This choice of parameters allows us to start with two types of galaxies, equal in mass and radius, but with very different inner density profiles (core vs. cusp).

For the initial velocities of particles in the frame of a given halo, we first use the Jeans equation (the first moment of the collisionless Boltzmann equation) to relate the density profile $`\rho (r)`$ and potential $`\mathrm{\Phi }(r)`$ to the velocity dispersions (Binney & Tremaine, 1987):

$$\frac{d(\rho \sigma _r^2)}{dr}+2\frac{\beta }{r}\rho \sigma _r^2=\rho \frac{d\mathrm{\Phi }}{dr},$$
(5)
where the velocity anisotropy is parameterized by

$$\beta =\beta (r)1\frac{\sigma _\theta ^2+\sigma _\varphi ^2}{2\sigma _r^2},$$
(6)
and $`\sigma _i^2`$ is the dispersion of velocity component $`i`$. If we further assume initial velocity isotropy, $`\beta (r)=0`$, we can integrate Eqn. (5) to yield the dispersion of the radial velocity as a function of radius,

$$\sigma _r^2(r)=\frac{1}{\rho (r)}_r^{\mathrm{}}𝑑R\rho (R)\frac{d\mathrm{\Phi }}{dR}.$$
(7)
From there, we randomly draw individual velocity components from a Gaussian distribution with variance $`\sigma _r^2(r)`$. Following Hernquist (1993), we impose the additional constraint $`v0.95v_{\mathrm{esc}}`$, which insures that all particles are initially bound to the halo.

We set up the two colliding haloes to be just touching at the truncation radii initially. The initial relative speed of the two haloes is taken to be the speed one of the haloes would have if the two were in circular orbit about their common centre of mass, i.e., $`v_{\mathrm{rel}}=2\sqrt{G\mu /r_{\mathrm{rel}}}`$, where $`\mu `$ is the reduced mass and $`r_{\mathrm{rel}}`$ is the distance between the halo centres. If the two haloes were point particles, this configuration would give $`v_{\mathrm{rel}}=v_{\mathrm{esc}}/\sqrt{2}`$; all our simulations are therefore relatively high-speed encounters for bound orbits. In this set-up, and with our choice of halo parameters, the two haloes first cross centres of mass at approximately 1.2 Gyr, then again at 2.2 Gyr. This leaves ample time after the first encounter for the haloes to interact and evolve within the 5.0 Gyr run.

Most of our simulations are head-on collisions with an impact parameter $`b=0`$. We have also performed two simulations in which the two halo centres are initially offset by $`b=0.5r_{200}`$ in order to test the robustness of our conclusions. In terms of the dimensionless parameter

$$\lambda =\frac{J|E|^{1/2}}{GM^{5/2}}$$
(8)
for the total mass $`M`$, energy $`E`$, and orbital angular momentum $`J`$, this corresponds to $`\lambda =0.042`$. In terms of the parameter $`ϵ`$, defined as the ratio of the angular momentum to that needed for circular orbit at the same velocity and separation, this set-up yields $`ϵ=0.24`$. As a final test of the generality of our findings, we performed two simulations of collisions involving unequal mass haloes with cores.

### 2.2 Simulation Parameters and Numerical Effects

All of our simulations are done using GADGET, a publicly available N-body tree code (Springel, Yoshida, & White, 2001). GADGET is available in serial and massively parallel forms; we have tested both versions on a local linux cluster and the parallel version on the IBM SP2 computers at the National Energy Research Scientific Computing Center (NERSC). In these test runs, we use either $`10^5`$ or $`10^6`$ particles and run for 5.0 Gyr. An important objective of the test runs is to check the level of stability against two-body relaxation for an isolated halo initially in equilibrium. Two-body relaxation is a numerical artifact of running simulations where each simulation particle is on the order of $`10^{60}10^{75}`$ times more massive than a typical elementary particle dark matter candidate. This leads to significant particle encounters that one does not expect to find in a realistic dark matter halo. As a result, the two-body relaxation timescale effectively puts a limit on the region of the halo that faithfully represents the phase space distribution in the collisionless limit.

Though there has been some disagreement about what regions can be considered reliable beyond the level of two-body relaxation (e.g. Moore et al. 1998; Klypin et al. 2001; Fukushige & Makino 2001; Power et al. 2003; Diemand et al. 2003), a conservative criterion is that the local two-body relaxation timescale $`t_r`$, which we define in terms of the circular velocity $`V_{\mathrm{circ}}=V_{\mathrm{circ}}(r)`$, period $`T(r)=2\pi r/V_{\mathrm{circ}}`$, and number of particles N (or equivalently mass M) interior to a radius r as

$`t_r(r)`$ $`=`$ $`T(r_{200}){\displaystyle \frac{N}{8\mathrm{ln}N}}{\displaystyle \frac{r}{r_{200}}}{\displaystyle \frac{V_{\mathrm{circ}}}{V_{\mathrm{circ}}(r_{200})}}`$ (9)
$`=`$ $`{\displaystyle \frac{\pi }{4}}{\displaystyle \frac{N}{\mathrm{ln}N}}\sqrt{{\displaystyle \frac{r^3}{GM(r)}}},`$
be at least as large as the age of the universe (Fukushige & Makino, 2001; Power et al., 2003). Since our simulations are purely gravitational, with no cosmological expansion, the current age of the universe has no direct meaning for the simulations themselves. As a result, we use the slightly different convergence criterion that the two-body relaxation time be longer than the length of time the simulation has run. We can invert Eqn. 9 to find $`r(t_r)`$, the minimum radius which satisfies $`t_r>t_f`$ for the final time of the simulation $`t_f`$; we refer to this as our minimum converged radius (see Table 1).

Our test runs of single isolated haloes show that the density profiles are stable at all radii $`r(t_r)<r0.8r_{200}`$ for the full 5.0 Gyr of simulation: any change in the profile is under 1%. (Near $`r_{200}`$ there are small deviations due to truncating the halo.) Inside $`r(t_r)`$, the behavior is dictated by the structure of the halo. NFW haloes have “temperature inversions” where the maximum velocity dispersion is reached not at the centre but at a finite radius. As a result, two-body relaxation does not cause a contraction of the central regions as is the case with the singular isothermal sphere (the so-called gravothermal catastrophe (Binney & Tremaine, 1987)) but rather an expansion as energy is transfered from regions with higher velocity dispersion to lower. Any model with a cusp shallower than $`\gamma =2`$ starts with a temperature inversion (Quinlan, 1996) and therefore undergoes an expansion of the central region and a concurrent reduction in central density. The central expansion is then followed by a collapsing stage, initiated once the temperature inversion is eliminated (Quinlan, 1996; Hayashi et al., 2003). This process occurs on the evaporation timescale, which is much longer than the length of our simulations at all radii. We do see evolution as a result of two-body interactions well within $`r(t_r)`$ for our isolated haloes with $`10^5`$ particles, but seems to be minimal even at $`0.5r(t_r)`$. We conclude that our simulations will be converged on all scales greater than $`r(t_r)`$ and are likely reliable to even smaller radii.

After obtaining stable isolated haloes, we performed production runs (all using the massively parallized version of GADGET) of halo collisions with different particle numbers to check for convergence in the inner structure of merger remnants. We did production runs with zero impact parameter using $`2\times 10^5`$ and $`2\times 10^6`$ total particles for the core-core collisions, and $`10^7`$ particles (in addition to both lower particle number simulations) for an cusp-cusp collision; for a core-cusp collision, we did one run with $`2\times 10^6`$ total particles. For the two runs with non-zero impact parameter, we used $`2\times 10^6`$ total particles. Our two unequal mass mergers were both core-core collisions with mass ratios of $`3:1`$. Both runs used the parameters described in Section 2.1 for the more massive galaxy. In one run (core4), we reduced all characteristic radii by a factor of $`3^{1/3}`$, yielding a core radius of $`5.55`$ kpc and a cut-off radius of $`55.5`$ kpc. In the other run (core5), the less massive galaxy had the same length scales as the more massive halo: $`\mathrm{r}_{\mathrm{core}}=8`$ kpc and $`\mathrm{r}_{\mathrm{cut}}=81`$ kpc, giving a halo that is spatially identical to the heavier one but with the density uniformly lower by a factor of 3. By using these parameters, we can identify how a merging halo having either a smaller core or a lower central density affects the merger process.

For all runs, we use a force softening of $`0.175`$ $`h^1`$ kpc for our particles, the new cell-opening accuracy parameter (errTolForceAcc) of 0.001, and timesteps that are inversely proportional to the acceleration (timestep criterion 1 of GADGET). (Our choice of errTolForceAcc is somewhat smaller than the maximum allowable value suggested by Power et al. (2003) in their recent thorough investigation of numerical parameters in cosmological simulations, but our tests have shown this value is necessary to obtain results that are independent of machine architechture and whether or not the run is parallel.) A typical run involved between 7000 (core-core runs) and 10000 (cusp-cusp runs) timesteps. Since the crossing time $`t_c=\sqrt{r_{200}^3/GM}`$ is approximately 1 Gyr for our haloes, an average timestep corresponds to approximately $`5.9\times 10^4`$ $`t_\mathrm{c}`$. The parameters for all of the production runs, as well as some results, are summarized in Table 1.

In all of our simulations, at least 65% of the total initial mass ends up within the virial radius of the final halo, meaning there are over $`6\times 10^6`$ particles within $`r_{200}`$ in our highest resolution simulation, a large enough number to reliably resolve density contrasts exceeding $`10^6`$. In all runs, at least 97% of the particles stays bound to the final halo. All mergers appear to be mostly completed by 3.5 Gyr, with the virial mass of the haloes increasing by only $`10`$% between 3.5 to 5.0 Gyr. Fig. 1 shows the density profiles at the final timestep for the three cusp-cusp collisions. The profiles of the remnants are extremely similar at all radii greater than $`2`$ $`h^1`$ kpc, and runs “Cusp2” and “Cusp3” have indistinguishable density profiles for radii greater than $`0.7`$ $`h^1`$ kpc. These results indicate two-body relaxation has not significantly affected any of our haloes at radii larger than $`2`$ $`h^1`$ kpc and that our two highest resolution cusp-cusp collisions yield density profiles that are reliable on sub-kiloparsec scales.

### 2.3 Definition of Halo Centres

The choice of how to define a halo’s centre is crucial in analyzing the final state of the merger remnant since quantities such as the spherically averaged density profile and radial velocity dispersion depend on the location of the centre. Previous studies have mainly used two methods for defining the centre of a halo: the most bound particle(s) (MBP) and the centre of mass (COM). For example, NFW (Navarro, Frenk, & White, 1996) use the COM; Moore et al. (1999a) use the MBP and note that it gives very similar results to the COM; Pearce et al. (1993) use the MBP and find that using the COM gives a consistently shallower inner profile.

We adopt an iterative version of the COM technique, which we find to be largely independent of initial assumptions. We have tested this method of computing a halo’s centre against using either the MBP or the COM of the five hundred most bound particles. The resulting spherically averaged density profiles are nearly identical at all radii that we consider to be converged – we find no evidence that one method systematically yields steeper cusps than the other. It is also worth noting that in all cases, the density profiles using different criteria for convergence of the centre of mass calculation differ only at radii that we already consider to be unreliable due to two-body relaxation and Poisson noise. As noted earlier, the density profiles are quite similar over a factor of 50 in particle number (see Fig. 1).

## 3 Results

### 3.1 Scatter Plots; Energy

Fig. 2 shows five snapshots at $`t=0`$, 1, 1.5, 2.5, and 5 Gyr from the merger simulation of a cored and a cuspy halo (run “Mixed”). (The general sequence of events is the same for all three types of mergers listed in Table 1.) To illustrate the evolution of the individual haloes, we also plot separately the cored halo (middle panels) and cuspy halo (right panels). During the encounter between 1 and 1.5 Gyr, the outer regions of the haloes pass through each other mostly unimpeded; the main effect in the outer portions is for the halo to become more extended. The panels at 2.5 Gyr show that the cuspy halo’s outer region is more easily stripped away than that of the cored profile; this is as expected, since the NFW profile is more dense in the centre but less tightly bound towards the outer region.

The merger remnants of the head-on collisions are noticeably prolate (extended along the collision axis) with axis ratios of $``$ 1.6:1:1 at half-mass radius. This result is quite similar regardless of the initial density profile. The two offset simulations “Core3” and “Cusp4” with orbital angular momentum result in mildly triaxial remnants.

Fig. 3 shows the evolution of the kinetic, potential and total energies as well as the virial ratio for the same run in Fig. 2. The initial crossing of halo centres is seen to occur at $``$ 1.2 Gyr, followed by a second crossing at $``$ 2.2 Gyr. Subsequent oscillations are strongly damped, and the system moves towards a steady state in virial equilibrium. The major features in Fig. 3 are common to all production runs listed in Table 1. All runs conserve energy to better than 0.4% and have virial ratios $`2T/|W|`$ within 2% of 1.0 by the final timestep.

### 3.2 Density Profile

In Fig. 4, we plot the spherically-averaged density profiles $`\rho (r)`$ of the merger products at 3.5 Gyr (dashed curve) and 5.0 Gyr (solid curve) for each of the three types of collisions: core-core (left), cusp-cusp (middle), and core-cusp (right). The density profiles in all cases are quite stable and constant from 3.5 to 5.0 Gyr, indicating that the merging process is essentially complete by 3.5 Gyr (also evident in Fig. 3) and that our simulations with $`2\times 10^6`$ particles do not suffer much from relaxation effects for the radii shown in Fig. 4. The dotted curves in Fig. 4 show our fitting results to be discussed below. For comparison, we overlay in a separate plot (Fig. 5) the final profiles of the merged haloes at 5 Gyr from the three collisions.

The final product of the core-core collision does not show signs of moving towards an NFW profile: the remnant has a core of approximately the same size as in the original haloes. In order to quantify our results, we have tried to fit the final profile to the same functional form as used for the initial cored profile (Eqn. 1). We find the inner region of the remnants well fit by a core radius of $`r_{\mathrm{core}}=7.8`$ $`h^1`$ kpc (for run “Core2”), which is nearly identical to the initial $`r_{\mathrm{core}}=8`$ $`h^1`$ kpc, but the density in the outer region of the remnants falls off more gradually, as $`r^3`$ instead of the initial exponential cutoff in Eqn. (1). Overall, we find the merged halo of core-core collisions better fit by

$$\rho (r)=\frac{A}{(r^2+B^2)^{3/2}},$$
(10)
with $`A5\times 10^{10}M_{}`$ and $`B11`$ $`h^1`$ kpc, than by the original profile (1). Since the remnants have virial radii about 20% larger than their progenitors, this means the relative size of the core in the final product is smaller, a result consistent with those of earlier simulations (White, 1978; Pearce, Thomas, & Couchman, 1993). The high spatial and mass resolution used here allows us to simulate scales more than five times below the core radius and to conclude that the core is preserved during one head-on major merger.

In the cusp-cusp collision, the product is essentially a new cuspy halo with an NFW-like form, which we find well fit by

$$\rho (r)=A\frac{B^3}{r^C(B+r)^{3C}}$$
(11)
for radius down to the relaxation scale $`r=0.7`$ $`h^1`$ kpc for the $`N_p=10^7`$ particle run. However, the best-fitting value for the logarithmic slope of the cusp is shallower than 1, with $`C0.7`$, perhaps hinting that collisions of galaxies serve to reduce the cusp at very small radii. The best-fitting parameters for the density profiles of our two highest resolution cusp-cusp collisions, as well as for the core-cusp collision, are summarized in Table 2.

In the case of the core-cusp merger, any hint of the core is entirely eliminated. The remnant has an overall density profile much closer to NFW than to the original core model and again well fit by Eqn. (11) with an inner slope of $`C0.65`$.

To understand the core-cusp merger in more detail, we calculate the spherically-averaged density profile for particles in each halo separately, first at $`t=0.0`$ Gyr and again at $`t=5`$ Gyr; the results are shown in Fig. 6. Interestingly, the particles originating from the cuspy halo maintain a cuspy, NFW-like profile (left panel), while the particles originating from the cored halo retain their original density profile to a remarkable extent as well (right panel). The final density profile is dominated by the cusp particles in the inner $`10`$ kpc and by the core particles from $`10100`$ kpc; outside of $`200`$ kpc, both types of particles seem to contribute approximately equally (also see Fig. 10).

The results above are all for head-on collisions with a zero impact parameter. In Fig. 7 we compare the density profiles of the non-zero impact parameter runs with the corresponding zero impact parameter simulations. The profiles are almost exactly the same, leading us to the conclusion that introducing moderate amounts of orbital angular momentum into a merging system will not influence our results.

To ensure that the results are not specific to equal-mass mergers, we plot the density profiles of our two $`3:1`$ mass mergers in Fig. 8, along with the density profile of the remnant in an equal mass merger for reference. The remnants of both unequal mass mergers retain cores of $`8`$kpc, the same result we found in the equal mass case. In fact, the final density profiles in all three are remarkably similar; this shows that core-core major mergers produce cored remnants.

### 3.3 Structure of Central Region

Although our simulations are collisionless, resulting in a phase space density that must be conserved due to Liouville’s theorem, we have no a priori reason to believe the configuration and momentum space distributions will be independently conserved, i.e. that a merger of two cored (cuspy) haloes will necessarily result in a cored (cuspy) remnant. A system’s coarse-grained distribution function $`\overline{f}`$ is constrained to be non-increasing as a function of time (Binney & Tremaine, 1987), and thus $`\overline{f}`$ in the final state is required at every point to be less than the maximal initial $`\overline{f}`$, but this requirement does not prevent cusps from becoming cores or vice-versa in collisionless simulations. In fact, violent collisions provide strongly fluctuating potentials that cause particles to rapidly exchange energy with the background, a process that could quickly affect the phase space structure of both initial haloes. Since the main objective of our study is to investigate the effects of major mergers on the phase-space distributions of dark matter in the inner regions of galaxy haloes, we examine this region in more detail in this section.

In the core-core and cusp-cusp collisions, we find that as expected, both progenitors contribute essentially equally to the particles at each radius of the merged halo. As Figs. 4 and 5 show, the inner density profile of the merger remnants is similar to their progenitors, i.e., cored haloes remain cored; cuspy haloes remain cuspy.

The situation is more intriguing for the mixed core-cusp collision. In order to elucidate the inner structure of this type of mergers, we first examine where the initial inner particles in each halo end up in the merged halo after 5 Gyr. (The results at 3.5 Gyr are very similar.) Fig. 9 shows the fraction of particles in each halo that start out within radius $`r`$ initially and also remain within $`r`$ at 5 Gyr. (By construction, this fraction approaches 1 at large $`r`$.) It shows that in the innermost several kiloparsecs of the merged halo, about 60% of the cusp particles (solid curve) that originated from within a given radius stay within that same radius, while very few of the core particles (dashed curve) remain. We believe it is this “lack of mobility” of the cusp particles that gives the cuspy halo a stronger influence at the final remnant’s centre and forces the final density profile to be cuspy rather than cored.

Next we examine the relative contribution of cusp vs core particles to the total density (or potential) in the core-cusp collision. Fig. 10 shows the ratio of core to cusp particles initially (dashed curve) and at 5.0 Gyr (solid curve). The figure shows that the cusp particles dominate the core particles number-wise in the centre initially (if the two galaxies were superimposed with a common centre), and that at 5 Gyr, the central region of the merger remnant becomes even more dominated by cusp particles. This result implies that the tightly bound particles in a cusp tend to stay bound together, while a core is much easier to disrupt. This complements the information from Fig. 9 by showing that for all particles of a given type, regardless of origin, the cuspy particles tend to take over the central region. Moving outward, the core particles begin to contribute more to the remnant. Note that the radius where the number of core particles is equal to the number of cusp particles stays almost constant at $`r20h^1`$ kpc. This is one indication that there is not a major reconfiguration of the particles occuring during the course of the simulation. We have also calculated the quantities in Figs. 9 and 10 at 3.5 Gyr and find them to be very similar to those found at 5 Gyr. This corroborates the evidence in Fig. 3 that the merger is essentially completed by 3.5 Gyr.

### 3.4 Velocity Profile

The particle velocities in our initial haloes are isotropic by construction, but those in our final haloes all exhibit some degree of anisotropy. For the cusp-cusp collision, Fig. 11 shows the velocity dispersions (top panels) and velocity anisotropy (bottom panels) as a function of halo radius $`r`$ at $`t=0`$, 3.5, and 5 Gyr for the simulation with $`10^7`$ particles (run “Cusp3” in Table 1).

The velocities remain mostly isotropic (i.e. $`\beta =0`$) near the centre of the merger remnant, but become mildly anisotropic ($`\beta <0.40`$) with increasing radius due to radial infalls.

We have also examined the velocity distribution function for the same cusp-cusp run.

We find that in the central region of the merged haloes, the velocities are quite Gaussian in addition to being isotropic; Fig. 12 shows $`f(v_i)`$ vs. $`v_i`$ for $`i=r,\theta ,\varphi `$ for $`10^4`$ particles in a radial bin arbitrarily centred on 10.106 kpc. The distribution function of $`v_r`$ (left panel) is well fit by a Gaussian (dashed line) with width $`\sigma _r=177`$ km s<sup>-1</sup>, computed from the particles in that bin. The distributions of $`v_\theta `$ (middle panel) and $`v_\varphi `$ (right panel) are also reasonably well fit by the corresponding Gaussians: $`\sigma _\theta =153`$ km s<sup>-1</sup> and $`\sigma _\varphi =161`$ km s<sup>-1</sup>. Moving away from the centre, we find the radial distribution becomes substantially less Gaussian and the angular distributions remain mostly Gaussian in the intermediate portion of the halo. Near the virial radius, however, none of the distributions are Gaussian.

The core-core collision exhibits similar trends as in Figs. 11 and 12, so we will not show the results here.

For the core-cusp collision, we plot in Fig. 13 the radial velocity dispersion $`\sigma _r`$ as a function of radius $`r`$ at $`t=0,2.5`$, and 3.5 Gyr for the core (crosses) and cusp (diamonds) particles separately. (Little evolution occurs after 3.5 Gyr.) The dispersion of the particles in the cuspy halo changes moderately between 0.0 and 3.5 Gyr, with the peak in $`\sigma _r`$ rising and moving outward but the overall shape maintained. For the particles in the cored halo, however, the velocity dispersion near the centre has almost doubled to $`\sigma _r230`$ km s<sup>-1</sup> after the merger.

To understand why the velocity dispersion of the core particles increases toward the centre of the merged halo upon merging with a cuspy halo, it is instructive to consider a situation where the two haloes are spherically symmetric and share a common centre. In this case, the potential $`\mathrm{\Phi }`$ is a function of the radial coordinate only and can be written as $`\mathrm{\Phi }(r)=\mathrm{\Phi }_1(r)+\mathrm{\Phi }_2(r)`$, a sum of the individual potentials $`\mathrm{\Phi }_1(r)=\mathrm{\Phi }_{\mathrm{core}}`$ and $`\mathrm{\Phi }_2(r)=\mathrm{\Phi }_{\mathrm{cusp}}`$. For particles in either the cored or cuspy halo, we then have

$`\sigma _r^2(r)`$ $`=`$ $`{\displaystyle \frac{1}{\rho (r)}}{\displaystyle _r^{\mathrm{}}}𝑑R\rho (R)\left({\displaystyle \frac{d\mathrm{\Phi }_1}{dR}}+{\displaystyle \frac{d\mathrm{\Phi }_2}{dR}}\right)`$ (12)
$`=`$ $`{\displaystyle \frac{1}{\rho (r)}}{\displaystyle _r^{\mathrm{}}}𝑑R{\displaystyle \frac{G\rho (R)}{R^2}}M_{\mathrm{tot}}(R),`$
where $`M_{\mathrm{tot}}(R)=M_1(R)+M_2(R)`$ is the total mass interior to radius $`R`$. From the point of view of the core particles, adding the cuspy particles has the effect of adding considerable extra mass and thus deepening the potential.

If the density profile of the core particles remains roughly constant, then adding extra mass in the form of cuspy particles essentially heats the core particles to higher a velocity dispersion, as seen in Fig. 13. The reverse is also true: the core particles also heat the cusp particles. This effect is actually seen in Fig. 13, as the cuspy particles’ dispersion increases relative to the original value, especially noticeably near the original peak of $`10`$ $`h^1`$ kpc. Since the cuspy particles dominate the total mass at radii less than this, however, the core particles are a small perturbation and do not change the cuspy particles’ dispersion significantly. The change in the peak of the NFW particles’ dispersion from $`10`$ to $`30`$ $`h^1`$ kpc can also be understood in this manner, since the core particles actually contribute more of the interior mass at this radius, as shown in Fig. 10.

### 3.5 Phase Space Profiles

The previous subsections discussed the density and velocity profiles separately. It is also instructive to examine the local phase space density $`\rho (r)/\sigma ^3(r)`$, which is shown in Fig. 14 for the three types of mergers studied in this paper. We find that for both the cusp-cusp (middle panel) and core-cusp (right panel) mergers, although the individual $`\rho (r)`$ and $`\sigma (r)`$ are not a pure power-law in $`r`$ (see Figs. 4, 11, and 13), the ratio is interestingly well-approximated by a power law $`\rho /\sigma ^3r^\alpha `$, where the best-fitting slope is $`\alpha 1.75`$ for both mergers. The core-core merger (left panel), on the other hand, shows a core in the phase space density since both $`\rho (r)`$ and $`\sigma (r)`$ are constant at $`r<r_{\mathrm{core}}`$.

Our power-law results can be compared with Taylor & Navarro (2001), who have suggested that the phase space structure of haloes is intimately related to that of self-similar collapse onto spherical density perturbations in an expanding universe. They solved the Jeans equation under the assumption of isotropy and power-law $`\rho /\sigma ^3r^\alpha `$ and found that haloes with the NFW profile are well-approximated by the same power law, $`\alpha =1.875`$, as the self-similar solutions in secondary spherical infall models (Bertschinger, 1985; Fillmore & Goldreich, 1984; Gott, 1975; Gunn, 1977). Since the phase space density is inversely proportional to the local entropy, Taylor & Navarro were led to the conclusion that CDM halo structure may be driven toward the most sharply peaked phase space distribution allowable for a monotonically decreasing density profile. For a power law phase space density, this solution yields density profiles that are very similar to the NFW profile but with a slightly shallower inner cusp, $`\rho r^{0.75}`$ rather than $`r^1`$.

Using the power-law result $`\rho /\sigma ^3r^{1.75}`$ in Fig. 14, we obtain from the Jeans equation a density profile $`\rho r^{0.72}`$, which is consistent with the best-fitting inner slope shown in Sec. 3.2. One can also look at the question in reverse, i.e., given a power-law density $`\rho r^\gamma `$, what solutions are admitted for the phase space density? Calculations show that $`\gamma 0.7`$ cusps lead to phase space densities well approximated by $`\rho /\sigma ^3r^{1.75}`$ as $`r0`$, but the slope of $`\rho /\sigma ^3`$ is not strictly a constant beyond the central region (Singh & Ma, in preparation).

## 4 Discussion

Our simulations offer insight into the mechanisms that drive a merging system toward an equilibrium structure, e.g., violent relaxation (Lynden-Bell, 1967) and phase mixing. Violent relaxation is due to a time-varying potential which allows the non-conservation of individual particles’ energies. The relative motion of the two haloes at the start of our simulations provides such a fluctuating potential. From Fig. 3, it is clear our mergers involve a rapid transfer of orbital energy to internal energy of the individual haloes (and subsequently the merger remnants). The decay of orbital motion occurs approximately on the crossing timescale, with essentially all energy in the form of the remnants’ internal energy by 2.5 Gyr. That the orbital motion is so heavily damped and that the remnant haloes do not expand or collapse significantly indicate that the time-varying potential necessary for violent relaxation is eliminated relatively quickly. One possible exception is that density waves could enable violent relaxation to persist even once the haloes have merged. Phase mixing tends to damp these waves strongly, however, with a decay timescale of $`\tau \lambda /\sigma `$ for a perturbation of size $`\lambda `$. Taking $`\sigma _r=150`$ km s<sup>-1</sup> to be typical for our merged haloes, this means that phase mixing will damp out all perturbations smaller than the virial radius within $`1`$ Gyr of the haloes merging. In our simulations, this seemingly short timescale for violent relaxation is sufficient to drive the velocity distribution function in the inner region of the merged halo to a Gaussian form, as predicted by statistical reasoning (e.g. Nakamura 2000). That the distribution function becomes less Gaussian with increasing radius is also in line with expectations because much less mixing is possible in the outer regions of the halo on the timescales we are studying.

Following the energy distribution through the collision can help in understanding the structure of the final product. A particularly interesting quantity is the change in energy $`\mathrm{\Delta }E`$ for a particle from initial to final state as a function of position. Funato, Makino, & Ebisuzaki (1992) describe an energy segregation in their collision simulations: particles with higher (less negative) energy tend to gain more energy than particles with lower energy. This means a particle near the centre of an initial halo should, on average, lose energy (i.e. become more tightly bound) relative to a particle near the virial radius of a halo. We follow the energy change of each particle, and plot in Fig. 15 the change in energy as a function of impact parameter for the core-cusp run “Mixed”. As Funato et al. describe, the particles with small impact parameters tend to lose energy, while the particles with large impact parameters tend to gain it. The effect is more dramatic for the cuspy particles than for the core particles, perhaps because the inner portion of the cuspy halo is initially more tightly bound but with an essentially equal velocity dispersion profile. As noted by Merral & Henriksen (2003), however, this energy segregation does not prevent the system from equilibrating via violent relaxation, at least in the central region.

Another interesting quantity is the differential energy distribution $`dN/dE`$, which can be written as a product of the phase-space distribution function $`f(E)`$ and the density of states $`g(E)`$ (Binney & Tremaine, 1987):

$$\frac{dN}{dE}=f(E)g(E).$$
(13)
The shape of $`dN/dE`$ is therefore *not* the same as the shape of $`f(E)`$. In Fig. 16 we plot $`dN/dE`$ for the core-cusp merger, for both types of particle, at the start of the simulation (thin curves) and at 5 Gyr (thick curves). It shows that $`dN/dE`$ is broadened, a natural outcome of mergers where the range of energies accessible to an individual particle can change drastically. The wave-like features in $`dN/dE`$ at 5 Gyr are indications of shells of particles formed via phase-mixing (Spergel & Hernquist, 1992). We note that unlike $`f(E)`$, which generally decreases with increasing $`E`$, the shape of $`dN/dE`$ is skewed towards higher energies because of the density of states. Both $`f(E)`$ and $`dN/dE`$ have been studied analytically for isolated NFW (Widrow, 2000) and King haloes (Binney & Tremaine, 1987) in equilibrium. The overall shape of $`dN/dE`$ resembles that shown in Fig. 16.

## 5 Summary

We have simulated major mergers of galaxy haloes with both cuspy and cored inner density profiles using GADGET, a versatile N-body tree code. Our findings can be summarized as follows:

1. Mergers of two cored haloes of equal mass result in a cored halo. The core radius of the merged halo is nearly identical to those of the original haloes, but the core to virial radius ratio is reduced, with the outer density profile of the merger remnant well approximately by $`\rho r^3`$. This result also holds for mergers with parent halo mass ratios of 3:1, implying all major mergers of two cored haloes result in cored haloes.
2. Mergers of a cuspy halo with either a cored halo or a second cuspy halo result in a cuspy halo. If the cuspy halo starts out with an NFW form, the merged halo also has an NFW-like form, with a slightly shallower inner cusp of $`\rho r^{0.7}`$ instead of $`r^1`$.
3. The particles from a given halo retain memory of their progenitor’s density profile. In major mergers of a cored halo with a cuspy halo, the particles originating from the cuspy halo end up with a cuspy profile, and those from the cored halo end up with a cored profile. (The cuspy particles dominate the total potential at the centre of the merged halo, so the resulting profile is cuspy.) The core particles in the inner 10 $`h^1`$ kpc of the remnant are significantly heated.
4. The merger remnants relax from the inside outwards, attaining Gaussian velocity distributions near their centres but not throughout the halo. Violent relaxation due to time-varying potential is effective only during the initial phase of the mergers; phase mixing is likely the dominant relaxation process at late times.
5. Including a moderate amount of orbital angular momentum in the form of a non-zero impact parameter has little discernible effect on our findings.

The set of simulations presented in this paper has been performed for mergers of non-spinning haloes in a set bound orbit. This clearly only probes a restricted range of parameter space. A possible area for future investigations is to include non-zero initial spins, which can influence tidal interactions significantly but may only affect the final phase space structure weakly (White, 1979). We have also chosen not to convolute the issues in this paper by ignoring minor mergers with satellite haloes and continuous accretion of diffuse material. The latter mechanisms may be responsible for generating the first generation of cuspy haloes from initially cored haloes, a process not achieved by major mergers of cored haloes according to our study.

Although the length and time scales quoted in this paper are for Milky-Way sized haloes of $`10^{12}M_{}`$, there is nothing in the numerical code or in our analysis that is special to this particular mass scale (except a mildly mass-dependent concentration parameter for NFW haloes). Our conclusions in this paper are thus likely to be applicable to major mergers of comparable-mass haloes over a wide mass range. Our finding that cuspy haloes are resilient to major mergers hints that the existence of cuspy haloes is natural in the hierarchical structure formation scenario: the build-up to larger mass haloes via repeated merging may be sufficient to produce cuspy haloes from cored haloes and preserve the general structure of cuspy haloes.

## Acknowledgments

We thank Jon Arons, Andrew Benson, Andrey Kravtsov, and Volker Springel for useful discussions and Jeff Filippini for help with the energy analysis. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098. C.-P. M is partially supported by an Alfred P. Sloan Fellowship, a Cottrell Scholars Award from the Research Corporation, and NASA grant NAG5-12173.


text:
# Modeling the Tidal Tails of the Sagittarius Dwarf Galaxy

## 1. Introduction

Since the discovery of the Sgr dwarf by Ibata et al. (1994) many groups (e.g., Johnston, Hernquist, & Bolte 1996, Ibata et al. 1997, Ibata & Lewis 1998, Gómez-Flechoso, Fux, & Martinet 1999, Johnston et al. 1999, Helmi & White 2001) have sought to model the Sgr - Milky Way interaction with respect to a modest patchwork of observational constraints. Recently, Majewski et al. (2003a, hereafter “Paper I”) have shown that the extensive length of the Sgr tidal tails can be traced by M giant stars visible in the all-sky view of the system provided by the 2MASS database. Spectroscopy of Sgr candidate stars has allowed determination of radial velocities throughout the trailing tail (Majewski et al. 2003b, hereafter “Paper II”), and these substantial new constraints can be used to develop more refined models of the Sgr system.

In this contribution, we briefly describe some of the major results of such modeling. A comprehensive description of this new Sgr disruption model can be found in Law, Johnston, & Majewski (2003, hereafter “Paper III”).

## 2. Modeling the Sgr System

Following previous work by Johnston et al. (1996, 1999) the Milky Way potential is represented numerically by a Miyamoto-Nagai (1975) disk, Hernquist spheroid, and a logarithmic halo. The total mass and radial profile are fixed by requiring that the rotation curve of this model Galaxy be consistent with HI & CO tangent point observations (e.g., Honma & Sofue 1997).

The Sgr dwarf itself is represented by $`10^5`$ self-gravitating particles (representing both the dark and light matter components of the satellite), which are initially distributed according to a Plummer (1911) model. This satellite is evolved through the simulated Galactic potential for five orbital periods using a self-consistent field code (Hernquist & Ostriker 1992). The present-day simulated dwarf is constrained to be located at $`(l,b)=(5.6^{},14.2^{})`$ at a solar distance of $`D_{\mathrm{Sgr}}=24`$ kpc (Paper I, Ibata et al. 1995) and have a radial velocity of $`v_{\mathrm{LOS},\mathrm{Sgr}}=171`$ km s-1 (Ibata et al. 1997). The direction of the dwarf’s space velocity vector is determined by requiring that the dwarf orbit in the orbital plane observed in Paper I.

Subject to these requirements, test-particle orbits (i.e. orbits calculated for a test particle with the observed kinematical characteristics of Sgr) and N-body simulations are performed for simulated satellites with a variety of orbital speeds. These simulations can be additionally constrained using the 2MASS M giant distance and radial velocity data presented in Papers I and II. Fig. 1 compares the M giant data (Panels a-b, filled squares) with the model Sgr dwarf whose tidal tails best reproduce the observations (Panels c-d). Note the close agreement between model and observed debris distances and radial velocities along the trailing debris tail ($`\mathrm{\Lambda }_{}=0^{}`$ \- $`100^{}`$)111We use the orbital longitude coordinate system in the Sgr orbital plane defined in Paper I.. This best-fit model is characterized by a period of 0.75 Gyr with apoGalacticon 52 kpc, periGalacticon 14 kpc, and a present space velocity of $`(U,V,W)=(237.2,43.4,218.9)`$ km s-1.

Although we do not attempt to model the Sgr core in detail, it is nonetheless possible to use the width of the Sgr debris stream to estimate such global characteristics as the bound mass of the dwarf. The simulated dwarf which appears to best fit the width of streams shown in Fig. 1 has a present mass of $`M_{\mathrm{Sgr}}=3\times 10^8M_{}`$ and a mass-to-light ratio $`M_{\mathrm{Sgr}}/L_{\mathrm{Sgr}}=21`$.

## 3. Discussion

As demonstrated in the previous section, the tidal tails of this model provide a good fit to the all-sky view of M giants presented in Papers I and II. It is therefore possible to use this model to determine what range of Milky Way models permit simulated satellites to reproduce observations. Particularly, N-body simulations can be used to constrain the flattening of the Galactic halo (e.g., Ibata et al. 2001). Fitting an orbital plane to leading and trailing M giant debris separately, we determine that the orbital pole of Sgr debris has precessed by $`1.7^{}\pm 2.4^{}`$ over about $`300^{}`$ of orbital longitude. Repeating this calculation for N-body simulations in model dark halos with a variety of flattenings, we calculate pole precessions of $`2.2^{}\pm 1.6^{}`$, $`3.5^{}\pm 1.7^{}`$, and $`5.6^{}\pm 1.4^{}`$ for flattenings in the halo potential of $`q=1,0.95`$, and $`0.90`$ respectively. It therefore appears likely that the halo of the Milky Way can be described by an almost spherical potential.

Although this model provides a good match to the distances and velocities of trailing Sgr debris given in Papers I and II, it does not fit recent data obtained by Majewski et al. (2003c, hereafter “Paper IV”) in the region of the Sgr leading arm. Fig. 1 (Panel a, filled triangles) plots these new data, which has velocities slower than that of the model by up to 200 km s-1 in the range $`\mathrm{\Lambda }_{}=300^{}`$ \- $`200^{}`$. There is no simple modification of the velocity of the model satellite that serves to reproduce this new trend, and this may be an indication of such other effects as dynamical friction. However, simulations suggest that including corrections from Chandrasekhar’s formulation of dynamical friction should not have a substantial effect on the observed velocities of leading tidal debris for model satellites with mass $`M_{\mathrm{Sgr},0}10^{10}M_{}`$, and we find that accurately reproducing the observed trend is difficult even for satellites with initial masses greater than this.

This inconsistency and implications of the best-fit model for the size and shape of the Milky Way are discussed at greater length in Paper III.

### Acknowledgments.

DRL acknowledges support from the Local Organizing Committee, a U.Va. Small Research Fellowship, and the U.Va. Echols Program. The authors also acknowledge NASA/JPL contracts 1228235 (SRM/KVJ) and 1234021 (MFS).

## References

Gómez-Flechoso, M. A., Fux, R., & Martinet, L. 1999, A&A, 347, 77

Helmi, A. & White, S. D. M. 2001, MNRAS, 323, 529

Hernquist, L. & Ostriker, J. P. 1992, ApJ, 386, 375

Honma, M. & Sofue, Y. 1997, PASJ, 49, 453

Ibata, R. A., Gilmore, G., & Irwin, M. J. 1994, Nature, 370, 194

Ibata, R. A., Gilmore, G. & Irwin, M. J. 1995, MNRAS, 277, 781

Ibata, R. A. & Lewis, G. F. 1998, ApJ, 500, 575

Ibata, R.A., Lewis, G. F., Irwin, M.J., Totten, E., & Quinn, T. 2001, ApJ, 551, 294

Ibata, R. A., Wyse, R. F. G., Gilmore, G., Irwin, M. J., & Suntzeff, N. B. 1997, AJ, 113, 634

Johnston, K. V., Hernquist, L., & Bolte, M. 1996, ApJ, 465, 278

Johnston, K. V., Majewski, S. R., Siegel, M. H., Reid, I. N., & Kunkel, W. E. 1999, AJ, 118, 1719

Law, D. R., Johnston, K. V., & Majewski, S. R. 2003, in prep. (“Paper III”)

Majewski, S. R., Skrutskie, M. F., Weinberg, M. D. & Ostheimer, J. C. 2003a, ApJ submitted (“Paper I”)

Majewski, S. R. et al. 2003b, ApJ submitted (“Paper II”)

Majewski, S. R. et al. 2003c, in prep. (“Paper IV”)

Miyamoto, M. & Nagai, R. 1975, PASJ, 27, 533

Plummer, H. C. 1911, MNRAS, 71, 460