text:
# Evidence of self-interacting cold dark matter from galactic to galaxy cluster scales

## 1 Introduction

Constraining the nature of dark matter is presently one of the most relevant problems in cosmology and particle physics. The current most popular scenarios for structure formation in the universe are based on the inflationary CDM theory, according to which cosmic structures arise from small Gaussian density fluctuations composed mostly of non-relativistic collisionless particles. Luminous galaxies are thought to form by gas cooling and condensing into the dark matter haloes which grow by gravitational accretion and merging in a hierarchical fashion.

The question on the inner density profiles of the virialized dark matter haloes is at present controversial. In the last few years much observational and theoretical effort has been employed into investigating the inner structure of dark haloes. On galactic scales, the rotation curves of dwarf galaxies offer a way to study the inner mass distribution of their dark haloes directly since these galaxies are dominated by dark matter. By analysing the rotation curves of some near dwarf galaxies, Moore (1994), Flores $`\&`$ Primack (1994) and Burkert (1995) have shown that the central mass distribution of their dark haloes is soft, i.e. the haloes have a constant density core. A similar result concerns low surface brightness galaxies (LSB, hereafter) (de Blok $`\&`$ McGaugh 1997) even though the uncertainty in the observational data is larger than in the case of dwarf spirals. Hernández $`\&`$ Gilmore (1998) showed that the observed rotation curves of both LSB and normal large galaxies are consistent with a fixed initial halo shape, characterized by a significant soft core inner region. On scales of clusters of galaxies, unfortunately there is not much information available. Recently, from strong gravitational lensing observations, Tyson, Kochanski, & Dell’Antonio (1998) have obtained an unprecedent high-resolution mass map for the cluster CL0024+1654, which has not a central cD galaxy, and found the existence of a soft core. Taken together these studies suggest the univesality of constant density cores across both large mass scales and galactic types.

On the theoretical side, the structure of the CDM haloes was studied over a wide range of masses by means of high-resolution N-body cosmological simulations (e.g., Navarro, Frenk, & White 1997; NFW hereafter) and semi-analytical approaches (e.g., Avila-Reese, Firmani, & Hernández 1998). It was found that the universal density profile firstly introduced by NFW describes very well the mass distribution of most of the CDM haloes. This profile is univocally determined by the mass, and in the centre diverges as $`\rho r^1`$ producing a cusp in the core. Recent high-resolution N-body simulations have shown that, as the numerical resolution is increased, the inner profiles result even steeper than $`r^1`$ (e.g., Moore et al. 1999b), making the CDM haloes more cuspy than in the case of the NFW profile.

So far, the predicted inner density profile of the CDM haloes seems to be in conflict with the observations. Another potential difficulty for the CDM models was recently reported: the N-body simulations predict an overly large number of haloes within group-like systems compared to observations (Klypin et al. 1999; Moore et al. 1999a). In light of these difficulties, the current stance of the hierarchical CDM-based scenario of structure formation remains somewhat confusing because, in fact, this scenario successfully accounts for: the distribution of matter at large scales (Bahcall et al. 1999), the uniformity of the cosmic microwave radiation and its small temperature anisotropies, and the observationally inferred cosmological parameters.

The aim of this letter is to analyse the halo core properties inferred from observations which might suggest explanations of the origin as to the soft halo cores and clarify the discrepancies that appear on small scales with the hierarchical scenario of structure formation. We investigate whether some modifications on the initial conditions of this scenario are able to improve the results with respect to the observations. We demonstrate that the introduction of self-interaction in the CDM particles as was suggested by Spergel & Steinhardt (1999) offers the most viable solution to the core problem in a context that preserves the hierarchical CDM-based scenario.

## 2 Halo central density from observations

We select from the literature dwarf and LSB galaxies with accurately measured rotation curves and clearly dominated by dark matter. These restrictions considerably reduce ambiguities in the estimates of the dark matter mass distribution due to uncertain stellar mass-to-light ($`M/L`$) ratios and modifications of the original halo profile produced by the gravitational drag of baryons during disc formation. Hence the dark haloes of these galaxies can be rightly assumed almost “virgin”. These constraints reduce the sample to six dwarf galaxies: DDO154 (Carignan et al. 1998), DDO170 (Lake et al. 1990), DDO105 (Schramm 1992, quoted by Moore 1994), NGC3109 (Jobin et al. 1990), IC2574 (Martimbeau et al. 1994), NGC5585 (Cóte et al. 1991). Six LSB galaxies are selected with the same criterion from a published sample: F568-v1, F571-8, F574-1; F583-1, F583-4, UGC5999 (de Blok $`\&`$ McGaugh 1997). The rotation curves measured for all these galaxies were used by the different authors to estimate the halo parameters, particularly the central density.

Our analysis also includes the density profile obtained for the cluster CL0024+1654 from a high resolution mass map derived using strong lensing techniques (Tyson et al. 1998). Because of the lack of a massive cD galaxy in the core, this cluster can be assumed to be dark matter dominated at the centre. Two clusters of galaxies, CL1455+22 and CL0016+16, with evident shallow mass profiles in the inner regions obtained by weak gravitational lensing studies (Smail et al. 1995) have also been considered, even though the uncertainty of the observational data is larger in these cases.

In Figure 1 we plot a very suggestive result: for a broad range of masses, the central density of the dark haloes is independent of mass (or circular velocity). Most dwarf galaxies (filled squares), LSB galaxies (open squares) and clusters (circles) indicate an average halo core density close to $`\rho _c=0.02M_{}/pc^3`$. The arrow shows a fiducial value derived from a published sample of LSB galaxies (de Blok $`\&`$ McGaugh 1997). The galaxy error bars are based on the observational uncertainty, and when possible from the range given by the maximum and minimum disc models. The cluster error bars take into account the uncertainty in observations and in a normalization factor of three in going from strong to weak lensing techniques (Wu et al. 1998).

This observational evidence makes the cosmological puzzle quite complex: how can one explain the origin of soft halo cores with roughly the same central density over the entire mass range sampled?

## 3 Shallow cores from collisionless cold dark matter

As was discussed above, observations seem to show that the inner density profile of the dark matter haloes is (i) shallow, and (ii) with a central value independent from the total halo mass (or maximum circular velocity $`V_m`$). These facts disagree with the predictions of the hierarchical CDM models. Now, we investigate some alternatives which might alleviate these difficulties within the cosmological context. For this we have performed a quantitative study of the CDM halo profiles using a semi-numerical method (Avila-Reese, Firmani, & Hernández 1998) aimed at calculating the collapse and virialization of spherically symmetric density fluctuations starting from an arbitrary mass aggregation history. Results obtained with this method are in excellent agreement with those of the N-body simulations (see Avila-Reese et al. 1999; Firmani & Avila-Reese 2000). The method is based on a generalization of the secondary infall model where non-radial motions and adiabatic invariance are taken into account. The only free parameter is the orbital parameter of particles (the perihelion to aphelion ratio) which regulates the thermal orbital energy of the system. This parameter is fixed independently of the halo mass and is constant during halo formation. Cosmological N-body simulations suggest $`r_{\mathrm{peri}}/r_{\mathrm{apo}}0.20.3`$ (see Ghigna et al. 1998).

Recently, Moore et al. (1999b) have simulated CDM haloes formed by monolithic collapse with N-body simulations introducing for this a lower cut-off at some wavelength in the power spectrum of fluctuations which suppresses substructures. The result was that the steep inner density profile of the haloes persisted. We suggest that this result might be partially a consequence of the lack of thermal orbital energy. In a monolithic collapse scenario the thermal orbital energy plays a significant role in producing soft cores: as $`r_{\mathrm{peri}}/r_{\mathrm{apo}}`$ increases a larger soft core is obtained. The density profiles of our haloes obtained for a CDM model with a lower cut-off in the variance of the power spectrum and a non zero initial thermal content, present soft cores. However, these models are unable to predict the observed central density trend shown in Figure 1 (Avila-Reese et al. 1998). In fact, the central density $`\rho _c`$ increases with $`V_m`$ in such a way that if $`\rho _c`$ is reproduced at galactic scales, for the cluster scales, $`\rho _c`$ overshoots the observed value by more than an order of magnitude. A hypothetical injection of thermal energy to the dark matter at a specific time in the life of the universe leads to a similar negative result.

An interesting way to produce soft halo cores in agreement with observations is to simply truncate the hierarchical halo mass aggregation histories at a given redshift towards the past. This may be done assuming that the halo mass fraction instantaneously collapses with some thermal energy (monolithic thermal collapse), while the rest of the mass is aggregated at the normal hierarchical rate. We have calculated the density profiles for haloes whose mass aggregation histories correspond to a hierarchical flat $`\mathrm{\Lambda }`$CDM model ($`\mathrm{\Omega }_m=0.3`$, h=0.7, $`\sigma _8=1`$) from $`z=5`$ and $`r_{peri}/r_{apo}=0.3`$; before this epoch the hierarchical aggregation was truncated. The results for this toy model are in good agreement with the observations: the haloes have a soft core, the core densities are independent from the mass and have a value similar to that what observational inferences indicate. It is interesting to note that the most distant QSOs and galaxies are at redshifts $`z5`$. Although the toy model presented here might look attractive, it is difficult to imagine a physical process capable of delaying the collapse of the central parts of the CDM haloes until $`z5`$.

## 4 Shallow cores from self-interacting cold dark matter

Self-interacting dark matter has been proposed as a possible solution for two potential conflicts of the hierarchical CDM models (Spergel & Steinhardt 1999; Hannestad 1999): the shallow core of the haloes and the dearth of dwarf galaxies in the Local Group. Astrophysical consequences of collisional dark matter have been pointed out by Ostriker (1999). It is easy to show that a configuration with the NFW density distribution is very far from thermal equilibrium: the inner velocity dispersion (temperature) has a positive gradient. Consequently, the presence of some self-interaction in the CDM particles introduces in the dark haloes a process of thermalization with heat transfer inwards, avoiding the formation of a cuspy profile. Heat capacity in the core is negative. This is a typical property of self-gravitating systems, like the interiors of the stars. For this reason, the heat transfer inwards cools the core exacerbating even more the temperature gradient. The heat transfer inwards increases causing the core to expand and cool due to gravothermal instability, leading to runaway core expansion. This physical mechanism is the key point for core expansion if self-interaction is effective. This process is similar to the post-collapse gravothermal instability well-known in dynamical studies of globular clusters (Bettwieser & Sugimoto 1983) where the minimum central density is reached roughly after a thermalization time.

The expansion of the core does not last forever. Since as the core expands the central density decreases, this would make the self-interaction less efficient and the core formation mechanism a self-limiting process. Although attractive, this mechanism is difficult to investigate because of our lack of knowledge regarding the cross section of the self-interacting dark matter particles. For this reason we start our analysis with a thermodynamical approach: we shall estimate the central density of CDM haloes assuming a thermodynamical equilibrium is reached due to strong self-interaction of the CDM particles. The final result will be the formation in the CDM halo of a central isothermal non-singular density profile established by competition between 1) mass and energy hierarchical aggregation, and 2) the thermalization due to self-interaction. The hierarchical mass and energy aggregation tends to stablish a NFW density profile (with the corresponding heat transfer inwards) while the self-interaction process tends to lead the system to a thermal equilibrium with the corresponding formation of a shallow core. For a given mass, the halo formed by a hierarchical mass aggregation identifies a gravitational binding energy (or $`V_m`$). Using this mass and binding energy to rescale a thermodynamical equilibrium configuration it is easy to find:

$$\rho _c=\alpha \frac{V_m^6}{M^2}M_{}/pc^3$$
(1)
where $`V_m`$ is in km/s, $`M`$ is the halo mass in M and $`\alpha `$ is a constant given by the detailed shape of the final equilibrium configuration. Since for the CDM haloes a tight relationship between their mass and circular velocity of the kind $`MV_m^n`$ with $`n3.2`$ is predicted (Avila-Reese et al. 1998,1999), eq. (1) implies that $`\rho _c`$ is roughly invariant with respect to the mass or $`V_m`$ as observations point out (Fig. 1). This strongly suggests that indeed a thermalization process due to dark matter self-interaction is acting in the CDM haloes.

Unfortunately, there is not a single final thermal equilibrium configuration, and as Lynden-Bell & Wood (1968) pointed out, some of configurations are even unstable. The King and Wooley configurations are examples of systems that have reached thermal equilibrium. They are characterized by a form parameter that may be related to the entropy of the system. A fiducial value for the central density of the CDM haloes with self-interaction may be estimated using a King or a Wooley profile at the state of maximum entropy (Lynden-Bell $`\&`$ Wood 1968). For these cases we derive respectively $`\alpha =\mathrm{1.3\; 10}^9`$ (short-dashed line in Fig. 1) and $`\alpha =\mathrm{2.6\; 10}^9`$ (long-dashed line in Fig. 1) in the appropriate units. The case of maximum entropy for a King profile corresponds to a value of the form parameter of $`8.5`$. A lower limit for the density may be roughly estimated from the dynamical evolution of globular clusters based on the Fokker-Planck approximation (Spitzer $`\&`$ Thuan 1972), starting from a uniform spherical distribution (this initial condition will lead to a central density lower than the density reached by the thermalization of a steep initial profile). The rescaling for this model taken at the first thermal equilibrium state gives us $`\alpha =\mathrm{1.7\; 10}^8`$ (dotted curve in Fig.1).

Global thermal equilibrium is reached when the self-interaction cross section is sufficiently large in order for the characteristic time scale of interactions across the overall halo to be shorter than the halo lifetime. An opposite situation of minimum cross section is given when self-interaction induces thermal equilibrium only in the region of the shallow core. In this case the central isothermal core appears surronded by a matter distribution characterised by a NFW profile. From the observational data it is possible now to infer an estimate of the self-interaction cross section. If $`n`$ is the dark particle number density, $`\sigma `$ the cross section and $`v`$ the dispersion velocity, assuming the collision time in the core $`\tau =1/(n\sigma v)`$ close to the Hubble time we obtain:

$$\frac{\sigma }{m_x}\mathrm{4\; 10}^{25}\left(\frac{0.02M_{}pc^3}{\rho _c}\right)\left(\frac{100kms^1}{v}\right)cm^2/GeV$$
(2)
with $`m_x`$ the mass of the dark matter particle and $`\rho _c`$ the central density. It is interesting to point out that for velocity dispersions corresponding to galaxy clusters this value is close to the upper limit estimated by Miralda-Escudé (2000) from the observationally inferred ellipticity of the cluster MS21137-23.

## 5 Summary

The discovery of a soft core in the cluster of galaxies CL0024+1254 by strong gravitational lensing measurements and the rotation curves of dark-matter dominated dwarf and LSB galaxies indicate that dark matter haloes have shallow inner density profiles from galactic to cluster scales. Studying in detail the observational data available for these cosmic objects, we found that the halo central density is nearly invariant with respect to the mass from galactic to cluster sizes.

We investigated different mechanisms and models for halo core formation within the hierarchical CDM scenario. We have shown that a lower cut-off at some wavelength in the CDM power spectrum and the assumption of high particle orbital thermal energies produce soft cores in the haloes, but the invariance of $`\rho _c`$ with respect to the mass is not reproduced. A more viable solution to the core problem is the introduction of self-interaction in the CDM particles. Being this the case, we proposed the gravothermal expansion as the mechanism responsible for the formation of soft cores in a hierarchical CDM scenario.

Using a thermodynamical approach we have estimated the central density of haloes in the case of maximum efficiency for self- interaction and found good agreement with the values inferred from observations. The central density in this case scales with the halo mass and its maximum circular velocity as $`\rho _cV_m^6/M^2`$. This result implies that $`\rho _c`$ is roughly constant because for the CDM haloes $`MV_m^n`$ with $`n3.2`$. If thermal equilibrium is restricted to the core, then the cross section given by eq.(2) may be derived consistently with observations. The cases analysed here, corresponding to a global and a local thermal equilibrium respectively, represent two limiting cases between which dark matter self-interaction may generate isothermal cores compatible with observations. We exclude from our analysis the extreme case of a very strong self-interaction which may lead the core to a gravothermal catastrophe with a central density profile steeper than NFW. Such extreme assumption of large cross section may be immediately ruled out because a singular isothermal core will be produced in contradiction with observations.

We stress the relevance confirming the existence of soft cores with scale invariant densities would have. In particular, the construction of high-resolution mass maps with gravitational lensing techniques for the inner regions of clusters is of great interest.

## Acknowledgments

ED thanks Fondazione CARIPLO for financial support.


format: text

text:
# Effect of field tilting on the vortices in irradiated Bi-2212.

## Abstract

We report on transport measurements in a Bi-$`2212`$ single crystal with columnar defects parallel to the c-axis. The tilt of the magnetic field away from the direction of the tracks is studied for filling factors $`f=B_z/B_\varphi <1`$. Near the Bose Glass transition temperature $`T_{BG}`$, the angular scaling laws are verified and we find the field independent critical exponents $`\nu ^{}=1.1`$ and $`z^{}=5.30`$. Finally, above $`H_C`$ we evidence the signature of a smetic-A like vortex phase. These experimental results provide support for the Bose Glass theory.

Columnar defects have been introduced in HTCS in order to avoid dissipation due to vortex motion. In the case of parallel tracks Nelson and Vinokur have predicted a transition between a so-called Bose Glass (BG) and vortex liquid, at a critical temperature $`T=T_{BG}`$. An important result is that field tilting does not destabilize the BG phase. Above a threshold transverse field $`H_C(T)`$ perpendicular to the columnar defects (CD), the flux lines accomodate simultaneously to CD and to the transverse field direction, defining a smectic-A like phase as proposed by Hwa and Nelson . Increasing further $`H_{}`$ leads to a vortex-liquid state . Recently, Grigera et al.  have evidenced a threshold transverse field $`H_C(T)`$ with resistivity measurements performed on a twinned Y-123 single crystal.

In this paper, we present results obtained on a Bi-2212 single crystal, irradiated with CD parallel to the c-axis with 5.8 GeV $`Pb`$ ions in GANIL (Caen- France). The CD density corresponds to a matching field $`B_\varphi =0.75`$ T. Isothermal I-V curves have been obtained varying the tilt angle between the magnetic field and the c-axis of the crystal. Measurements were performed keeping the $`B_z`$ component constant, where the z-axis is along the c-axis of the sample. In the insert of Fig. 1, the log-log plot of isothermal ohmic resistance versus $`H_{}/H_z`$ is displayed above and below $`T_{BG}`$. For $`T>T_{BG}`$, ohmic behaviour is detected, within the experimental sensitivity, even for vanishing transverse field. In contrast, for $`T<T_{BG}`$ the ohmic resistance goes to zero at some critical tilt showing the existence of a critical transverse magnetic field $`H_C(T)`$. In the BG theory, scaling functions $`f_+`$ and $`f_{}`$ are predicted to describe the resistivity in the presence of a tilted magnetic field, above and below $`T_{BG}`$, respectively. The resistance then reads:

$$R|t|^{\nu ^{}(z^{}2)}f_\pm \left((H_{}/H_z)|t|^{3\nu ^{}}\right).$$
(1)
where $`t=(TT_{BG})/T_{BG}`$ is the reduced temperature, and $`\nu ^{}`$ and $`z^{}`$ are critical exponents. Fig. 1 displays such scaling properties. We find that the scaling functions $`f_+`$ and $`f_{}`$ collapse for all the filling fractions $`f=B_z/B_\varphi <1`$ investigated in our experiment. Moreover, we obtain field independent critical exponents, $`\nu ^{}=1.1\pm 0.1`$ and $`z^{}=5.30\pm 0.05`$. These results are in good agreement with the BG theory  and other experimental results .

We shall now consider the case $`T<T_{BG}`$. For $`H_{}>H_C(T)`$, the vortex motion is mediated by kinks aligned in chains in the direction of the transverse magnetic field $`H_{}`$ , in such a way that the chain density is directly related to the linear resistance:

$$Rn_{chain}(H_{}H_C(T))^{3/2}.$$
(2)
Scaling arguments lead to a critical transverse magnetic field $`H_C(T)`$ vanishing as $`TT_{BG}^{}`$ as:

$$H_C(T)|t|^{3\nu ^{}}.$$
(3)
Figure 2 displays a plot of $`R^{\frac{2}{3}}`$ versus $`H_{}/H_z`$. The solid lines are successfull fits of Eq. 2 to data. This result does support the existence of a smectic-A behaviour. The intersection of a linear fit with the abscissa-axis directly gives the critical transverse magnetic field $`H_C(T)`$ at a given temperature. The insert of Fig. 2 shows a log-log plot of $`H_C(T)`$ thus obtained versus the reduced temperature $`t=(TT_{BG})/T_{BG}`$. The solid line represents a least-square fit of Eq. 3 to data. According to Eq. 3, we find therefrom the critical exponent value $`\nu ^{}=1.1\pm 0.1`$. Note that this value is consistent with the one we found above by another way, as predicted by Nelson and Vinokur .

In conclusion, we have investigated the Bose Glass phase transition on an irradiated Bi-2212 single crystal versus both the filling fraction $`f=B_z/B_\varphi <1`$ and the magnetic field tilt. We have shown that the Bose Glass transition in the presence of a tilted field verifies the scaling rules predicted by Nelson and Vinokur with both field independent scaling functions and critical exponents $`\nu ^{}=1.1`$ and $`z^{}=5.30`$. A smectic-A like behaviour has been evidenced. Finally, the critical tilt $`H_C/H_z`$ , separating this phase from the Bose Glass one, has been found to vary as $`H_C/H_z|t|^{3\nu ^{}}`$ in agreement with the theorical expectation of a sharp cusp in the $`TH_{}/H_z`$ phase diagram.


text:
# Vortex Pinning and Dynamics in Layered Superconductors with Periodic Pinning Arrays

## Abstract

We examine vortex dynamics and pinning in layered superconductors using three-dimensional molecular dynamics simulations of magnetically interacting pancake vortices. Our model treats the magnetic interactions of the pancakes exactly, with long-range logarithmic interactions both within and between planes. At the matching field the vortices are aligned with the pinning array. As a function of tilt angle for the pinning arrays a series of commensuration effects occur, seen as peaks in the critical current, due to pancakes finding a favorable alignment.

In superconductors with periodic pinning arrays interesting commensurability effects occur when the periodicity of the vortex lattice matches the periodicity of the pinning lattice. Experiments  and simulations  so far have been done with thin film superconductors where the vortex lattice and pinning can be considered two-dimensional. The case of vortex lattices interacting with a periodic pinning array in a layered 3D superconductor has not been studied. Such a system would correspond to an anisotropic superconductor such as BSCCO with a periodic arrangement of columnar defects. In this system the $`z`$-direction becomes important as the applied field or the pinning array is tilted. The dynamical effects of vortices moving in periodic pinning arrays in such a system have not been examined, in particular how the vortex lattice structure of the moving state differs from that of the pinned state. To study vortex pinning and dynamics in layered superconductors, we have developed a simulation containing the correct magnetic interactions between pancakes . This interaction is long range both in and between planes, and is treated using a rapidly converging summation method .

The overdamped equation of motion, for $`T=0`$, for vortex $`i`$ is given by $`𝐟_i=_{j=1}^{N_v}𝐔(\rho _{i,j},z_{i,j})+𝐟_i^{vp}+𝐟_d=𝐯_i`$, where $`N_v`$ is the number of vortices and $`\rho `$ and $`z`$ are the distance between pancakes in cylindrical coordinates. The magnetic energy between pancakes is

$`𝐔(\rho _{i,j},0)=2dϵ_0\left((1{\displaystyle \frac{d}{2\lambda }})\mathrm{ln}{\displaystyle \frac{R}{\rho }}+{\displaystyle \frac{d}{2\lambda }}E_1(\rho )\right)`$
$`𝐔(\rho _{i,j},z)={\displaystyle \frac{d^2ϵ_0}{\lambda }}\left(\mathrm{exp}(z/\lambda )\mathrm{ln}{\displaystyle \frac{R}{\rho }}E_1(R)\right)`$
where $`R=\sqrt{z^2+\rho ^2}`$, $`E_1(x)=_\rho ^{\mathrm{}}\mathrm{exp}(x/\lambda )/\rho ^{}𝑑\rho ^{}`$ and $`ϵ_0=\mathrm{\Phi }_0^2/(4\pi \xi )^2`$. The pinning is placed in a square array of parabolic traps with a radius $`r_p`$ much smaller than the distance between pins. The location of the pinning sites is the same in every layer corresponding to correlated defects. A driving force $`f_d`$ is slowly increased and the vortex velocities are measured. Here we consider the first matching field case where the number of vortices $`N_v`$ equals the number of pinning sites $`N_p`$. We conduct a series of simulations in which the pinning sites are tilted at an increasing angle with respect to the $`z`$-axis. We will only consider driving that produces vortex motion transverse to the direction of the tilt angle. We examine systems with 8 layers containing 64 vortices and pins in each layer. Work for larger systems, varied fields and coupling strength will be presented elsewhere .

In Fig. 1(a) we present the critical depinning force $`f_{dp}^c`$ as a function of tilt angle $`\theta `$. Here $`f_{dp}^c`$ peaks at $`\theta =0^{}`$ when the pancakes are aligned with pins on all layers. As $`\theta `$ is increased $`f_{dp}^c`$ drops. For small tilt angles $`\theta <5^{}`$ the vortex lines tilt with the pins. For larger angles the vortex lines realign in the $`z`$ direction. The depinning force $`f_{dp}^c`$ will then remain low as only one pancake in the straight vortex line will be sitting at a pinning site. At $`\theta =45^{}`$ $`f_{dp}^c`$ shows a peak of the same magnitude as the peak at $`\theta =0`$. At this tilt angle, and also for any angle satisfying $`\theta =\mathrm{tan}^1(n)`$ where $`n`$ is an integer, the pinning sites are again aligned in the $`z`$-direction so that a vortex line can be formed that is also aligned in the $`z`$-direction with all the pancakes in a single vortex being able to sit in a pinning site. There are also peaks in $`f_{dp}^c`$ at $`\theta =26.6^{}`$ and $`56.3^{}`$. At these angles the pancakes again sit on all the pinning sites. The individual vortex lines now consists of half the number of pancakes as at $`\theta =0.0^{}`$; however, there are now twice as many vortex lines with the pancakes from an individual vortex line being coupled in every other layer. The view from the $`z`$-direction as shown in Fig. 1 for these angles indicates that the vortex lattice is now rectangular with twice as many vortex lines as at the other angles. At $`\theta =36.9^{}`$ a smaller peak is observed. The vortex structure at this angle will be presented elsewhere .

In (b) and (c) we show the vortex structures for the pinned phase and moving phase for $`\theta =1.5^{}`$ as seen from the $`z`$-direction. In (b) the vortices can be seen to stay aligned with the pins. In (c) for $`f_d>f_{dp}^c`$ the vortices realign with the z-direction. Such a transition from a tilted to straight vortex lattice as a function of drive may be visible with neutron scattering experiments.

We acknowledge helpful discussions with L. N. Bulaevskii, A. Kolton, R.T. Scalettar, and G. T. Zimányi. This work was supported by CLC and CULAR (LANL/UC) and by the Director, Office of Adv. Scientific Comp. Res., Div. of Math., Information and Comp. Sciences, U.S. DoE contract DE-AC03-76SF00098.


text:
# Simple dynamical models of the Sagittarius dwarf galaxy

## 1 Introduction

The Sagittarius dwarf galaxy is the closest satellite of the Milky Way (Ibata, Gilmore & Irwin 1994, 1995, hereafter IGI95). Soon after its discovery, several groups carried out simulations to see if its properties are consistent with the disruption of an object similar to the other dwarf companions of the Milky Way, but none produced a model in full agreement with both the age and the structure of the observed system (Johnston, Spergel & Hernquist 1995; Velázquez & White 1995; Edelsohn & Elmegreen 1997; Ibata et al. 1997, hereafter I97; Gómez-Flechoso, Fux & Martinet 1999). All groups assumed light to trace mass and an initial system similar to observed dwarf spheroidals. All found the simulated galaxy to disrupt after one or two orbits whereas the observed system has apparently completed ten or more. Most considered this to be a problem (but cf Velázquez & White 1995). As a result, several unconventional models were proposed to explain the survival and structure of Sagittarius. In an extensive numerical study, Ibata & Lewis (1998) concluded that Sagittarius must have a stiff and extended dark matter halo if it is to survive with 25% of its initial mass still bound today. Since an extended halo cannot remain undistorted in the Galaxy’s tidal field for any conventional form of dark matter, it is unclear how this idea should be interpreted. Furthermore, it produces an uncomfortably large mass-to-light ratio ($`100`$), it cannot reproduce the observed elongation, and it suggests that little tidal debris will be liberated, in apparent conflict with the observations of Mateo, Olszewski & Morrison (1998), and Majewski et al. (1999) (see also Johnston et al. 1999). A somewhat less unorthodox model was proposed by Zhao (1998), where Sagittarius was scattered onto its current tightly bound orbit by an encounter with the Magellanic Clouds about 2 Gyr ago. This appears physically possible but requires careful tuning of the orbits of the two systems (see Ibata & Lewis 1998; and Jiang & Binney 2000). Another mechanism by which the dwarf could have moved to a short-period orbit is dynamical friction, which can be important only if Sagittarius has lost a lot of mass in the past. Jiang & Binney (2000) found a one-parameter family of initial configurations that evolve into something like the present system over a Hubble time. Their initial systems have masses $`10^{1011}\mathrm{M}_{}`$ and start from a Galactocentric radius $`200`$ kpc.

Driven by this apparent puzzle, we decided to search more thoroughly for a self-consistent model of the disruption of Sagittarius, which, after a Hubble time, has similar characteristics to those observed. (See Table 1 for a summary of the observed properties of the system.) Below we present two models which meet these requirements.

## 2 Method

In our numerical simulations, we represent the Galaxy by a fixed potential with three components: a dark logarithmic halo

$$\mathrm{\Phi }_{\mathrm{halo}}=v_{\mathrm{halo}}^2\mathrm{ln}(r^2+d^2),$$
(1)
a Miyamoto-Nagai disk

$$\mathrm{\Phi }_{\mathrm{disk}}=\frac{GM_{\mathrm{disk}}}{\sqrt{R^2+(a+\sqrt{z^2+b^2})^2}},$$
(2)
and a spherical Hernquist bulge

$$\mathrm{\Phi }_{\mathrm{bulge}}=\frac{GM_{\mathrm{bulge}}}{r+c},$$
(3)
where $`d`$=12 kpc and $`v_{\mathrm{halo}}`$ = 131.5 $`\mathrm{km}\mathrm{s}^1`$; $`M_{\mathrm{disk}}=10^{11}\mathrm{M}_{}`$, $`a`$ = 6.5 kpc and $`b`$ = 0.26 kpc; $`M_{\mathrm{bulge}}=3.4\times 10^{10}\mathrm{M}_{}`$ and $`c`$ = 0.7 kpc. This choice of parameters gives a flat rotation curve with an asymptotic circular velocity of 186 $`\mathrm{km}\mathrm{s}^1`$. The mass of the dark-matter halo within 16 kpc is 7.87 $`\times 10^{10}\mathrm{M}_{}`$ in this model.

We represent the satellite galaxy by a collection of $`10^5`$ particles and model their self-gravity by a multipole expansion of the internal potential to fourth order \[White 1983, Zaritsky & White 1988\]. This type of code has the advantage that a large number of particles can be followed in a relatively small amount of computer time. Hence a substantial parameter space can be explored while retaining considerable detail on the structure of the disrupted system. In this quadrupole expansion, higher than monopole terms are softened more strongly. We choose $`ϵ_10.20.25r_c`$ for the monopole term ($`r_c`$ is the core radius of the system) and $`ϵ_2=2ϵ_1`$ for dipole and higher terms and for the centre of expansion. The centre of expansion is a particle which, in practice, follows the density maximum of the satellite closely at all times.

For the stellar distribution of the pre-disruption dwarf we choose a King model (King 1966), since this is a good representation of the distant dwarf spheroidals. King models are defined by a combination of three parameters: $`\mathrm{\Psi }(r=0)`$ (depth of the potential well of the system), $`\sigma ^2`$ (measure of the central velocity dispersion), and $`\rho _0`$ (central density) or $`r_0`$ (King radius). The ratio $`\mathrm{\Psi }(r=0)/\sigma ^2`$ defines how centrally concentrated the system is, and for any value of this parameter, a set of homologous models with different central densities and core (or King) radii may be found. We assume that the progenitor of Sagittarius obeys the known metallicity-luminosity relation for the Local Group dSph \[Mateo 1998\]. The metallicity determinations for Sagittarius \[I97\] indicate $`[\mathrm{Fe}/\mathrm{H}]1`$, corresponding to a total luminosity in the range $`3.5\times 10^73.5\times 10^8\mathrm{L}_{}`$. To obtain an initial guess for the mass of the system, we transform this luminosity into a mass assuming a mass-to-light ratio $`2`$. The relevant initial stellar mass interval is then $`7\times 10^77\times 10^8\mathrm{M}_{}`$.

Note that our choice of a fixed potential to represent our Galaxy means that we neglect any exchange of energy between the satellite and the Galactic halo. This is an excellent approximation for the range of orbits and satellite masses that we consider, since these imply dynamical friction decay times substantially in excess of the Hubble time. The orbits are also sufficiently large that impulsive heating during disk passages can be neglected.

The orbit of Sagittarius is relatively well constrained \[I97\]. The heliocentric distance $`d25\pm 2`$ kpc and position $`(l,b)=(5.6^{},14^{})`$ of the galaxy core are well-determined; the heliocentric radial velocity $`v_r^{\mathrm{sun}}140\pm 2\mathrm{km}\mathrm{s}^1`$, and its variation across the satellite are also accurately measured. Outside the main body ($`b<20^{}`$) the radial velocity shows a small gradient $`\mathrm{d}v_r/\mathrm{d}b3\mathrm{km}\mathrm{s}^1\mathrm{deg}`$, but no gradient is detected across the main body itself. The proper motion measurements are not very accurate; $`\mu _b2.1\pm 0.7\mathrm{mas}\mathrm{yr}^1`$, and no measurement is available in the $`l`$-direction. On the other hand the strong North-South elongation of the system suggests that it has little motion in the $`l`$-direction, thus implying the orbit should be close to polar. We generate a range of possible orbits satisfying these constraints and concentrate on those with relatively long periods in order to maximise the survival chances of our satellite. We begin all our simulations half a radial period after the Big Bang to allow for the initial expansion. We place the initial satellite at apocentre, then we integrate forward until $`13`$ Gyr. The orbits are chosen so at this time the position and velocity of the satellite core correspond to those observed. We allow ourselves some slight freedom in choosing the final time in order to fit the observed data as well as possible.

## 3 Results

Figure 1 gives an example of an orbit which is consistent with all the current data on Sagittarius. It has a pericentre of $`16.3`$ kpc, an apocentre of $`68.3`$ kpc, and a radial period of $`0.85`$ Gyr. We use similar orbits for all the simulations described below. Note that the slow precession about the Galactic rotation axis is in part due to the quasi-polar nature of the orbit and in part to the fact we have assumed the Milky Way’s dark halo to be spherical.

After letting our satellite relax in isolation, we integrate each simulation for $`13`$ Gyr. In practice we needed to run a large number of simulations, and test each to see if it satisfies the observational constraints at the present time. Since it remains uncertain whether dwarf spheroidals have extended dark halos (e.g. Klessen & Kroupa 1998), we have considered both purely stellar models and models in which the initial stellar system is embedded in a more massive and more extended dark halo.

### 3.1 Constant mass-to-light ratio: A purely stellar model

Our preferred purely stellar model (Model I) initially has a core radius of $`r_c`$ = 0.44 kpc, a total velocity dispersion of 18.9 $`\mathrm{km}\mathrm{s}^1`$, and a concentration parameter $`c=\mathrm{log}_{10}(r_t/r_c)0.83`$. This implies a total mass of $`M=4.66\times 10^8\mathrm{M}_{}`$. For a satellite to survive for about 10 Gyr on an orbit with pericentre $`15`$ kpc, apocentre $`70`$ kpc, and period $`1`$ Gyr (for which the observational constraints are satisfied) its initial central density has to be $`\rho _00.360.4\mathrm{M}_{}\mathrm{pc}^3`$. Satellites with significantly smaller initial densities do not survive long enough.

In Figure 2 we plot heliocentric distance as a function of galactic latitude for stars projected near the main remnant 12.5 Gyr after infall. Streams of particles are visible at all latitudes over a broad range in distance. Sagittarius has been orbiting long enough for its debris streams to be wrapped several times around the Galaxy. (See also Figure 8.)

The remnant galaxy, i.e. the central region of the satellite’s debris, is similar to the real system. In Figure 3 we plot its mass surface density. The transformation from observed surface brightness to mass surface density (which is what the simulations give us) can be done as follows. The observed mass surface density $`\mathrm{\Sigma }`$ for an assumed mass-to-light ratio $`\mathrm{{\rm Y}}`$ is

$$\mathrm{\Sigma }=\frac{N_XL_X}{f_X}\mathrm{{\rm Y}}\left[\frac{\mathrm{M}_{}}{\mathrm{deg}^2}\right],$$
(4)
where $`N_X`$ is the number of observed stars of type $`X`$ per square degree, $`L_X`$ is their luminosity, and $`f_X`$ is the fraction of the total luminosity in stars of type $`X`$. In IGI95 the spatial structure of Sagittarius was determined from the excess of counts at the apparent magnitude of the horizontal branch. Uncertainties in the result are due primarily to contamination by sources in the Galactic bulge. Their lowest isodensity contour is at $`\mathrm{\Sigma }_{\mathrm{min}}5\times 10^5\frac{\mathrm{M}_{}}{\mathrm{deg}^2}`$, assuming $`\mathrm{{\rm Y}}2.25`$ and \[Fe/H\] $`1`$ (Bergbusch & vandenBerg 1992), and has an extent of $`7.5^{}\times 3^{}`$. This same isodensity contour is shown in Figure 3 as a thick line. It has an extent of $`8^{}\times 4.8^{}`$, in reasonable agreement with the observations given the uncertainties. In I97 isodensity contours were derived from counts of main sequence stars close to the turn-off, roughly one magnitude above the plate limit. The minimum contour in this case corresponds to $`\mathrm{\Sigma }_{\mathrm{min}}10^5\frac{\mathrm{M}_{}}{\mathrm{deg}^2}`$, and has an extent of roughly $`15^{}\times 7^{}`$. In Fig. 3 this contour is shown as a dashed-line, and has an extent of $`21^{}\times 6.5^{}`$, also in good agreement with the observations.

Note that the isophotes (or isodensity contours) become rounder towards the centre of the satellite. Its angular core radius is $`R_c1.29^{}`$, which for a distance of 26 kpc (derived from the simulations) corresponds to 0.58 kpc, again in good agreement with the observations.

The kinematic properties of the remnant galaxy are more difficult to compare with observations because a substantial amount of mass from debris streams is projected on top of the main body. Like I97, we measure the radial velocity across the system considering only particles for which $`100\mathrm{km}\mathrm{s}^1v_r^{\mathrm{sun}}180\mathrm{km}\mathrm{s}^1`$. In the left panel of Figure 4 we plot the heliocentric radial velocity, and in the right panel we plot its dispersion as a function of Galactic latitude. For comparison, we analysed the observations of I97 at CTIO in the same way (their Table 2b); these data have a precision of a few $`\mathrm{km}\mathrm{s}^1`$ (triangles in Figure 4). Our model is consistent with the observed kinematics; we obtain a heliocentric radial velocity of $`139.5\mathrm{km}\mathrm{s}^1`$ and an internal velocity dispersion in the radial direction of $`11\mathrm{km}\mathrm{s}^1`$ for the main body. However, when the radial velocity restrictions for inclusion in this calculation are relaxed, we find much larger velocity dispersions because of the contribution of stars from other streams. It is important to consider this problem when determining which stars should be considered members of Sagittarius.

### 3.2 Varying mass-to-light ratio: A model with a dark halo

The observational data for Sagittarius mainly refer to the current remnant system, which corresponds to the innermost regions of the progenitor satellite. As a consequence, models that are initially dark matter dominated in their outskirts are relatively poorly constrained.

As an example we focus on a progenitor with a mass distribution which is similar to that of Model I in its inner regions, but is considerably more extended. We take the mass distribution to be a (heavy) King model with $`r_c=0.54`$ kpc and $`r_t=10.4`$ kpc, with an initial total velocity dispersion of 25.2 $`\mathrm{km}\mathrm{s}^1`$, and total mass of $`M=1.7\times 10^9\mathrm{M}_{}`$. For an orbit like that of Model I this produces a suitable remnant after 12 Gyr. The mass distribution of this remnant satisfies many of the observational constraints of Table 1. Its core radius is slightly larger $`r_c0.65`$ kpc, and the radial velocity dispersion in the main body is 12.1 $`\mathrm{km}\mathrm{s}^1`$.

We will construct a two-component satellite with this mass distribution by solving for the dependence of mass-to-light ratio on initial binding energy that produces the initial light profile of Model I. We choose the mass-to-light ratio of satellite material to be a decreasing function of binding energy, so that the most bound particles have near “stellar” mass-to-light ratios, whereas weakly bound particles are almost entirely “dark”. From the energy distribution of the heavy King model, and that of a King model with $`r_0=0.095`$ kpc and $`\sigma =25.6\mathrm{km}\mathrm{s}^1`$, we can derive the mass-to-light ratio as a function of binding energy as

$$\mathrm{{\rm Y}}(ϵ)=\mathrm{{\rm Y}}_{}\frac{\mathrm{d}M/\mathrm{d}ϵ(ϵ)}{\mathrm{d}M_{}/\mathrm{d}ϵ(ϵ=ϵ_{}+ϵ_{\mathrm{max}}ϵ_{\mathrm{max}})}$$
(5)
where $`\mathrm{{\rm Y}}_{}`$ is the mass-to-light ratio of a stellar population. The energies $`ϵ_{}`$ of the lighter King model have been shifted by a fixed amount $`ϵ_{\mathrm{max}}ϵ_{\mathrm{max}}`$, to be on the same scale as that of the heavier King model. The resulting mass-to-light ratio is shown in Figure 5.

In Figure 6 we show the surface mass densities normalized to their central values for Model I (only stars), for the heavy King model and for the two-component model (“stars” and dark-matter). We shall refer to this two-component model as Model II, which is obtained by weighting each simulation particle by $`\mathrm{{\rm Y}}(ϵ)^1`$.

If we require that the central stellar mass surface densities of Model I and Model II be the same, we find that the total mass in stars in Model II is $`1.69\times 10^8\mathrm{M}_{}`$. To match Sagittarius surface brightness, we choose the central stellar mass-to-light ratio $`\mathrm{{\rm Y}}_{}=1.5`$. Thus, the total luminosity of Model II is then $`1.13\times 10^8\mathrm{L}_{}`$, implying a mass-to-light ratio of 15.1. Its initial velocity dispersion is 23 $`\mathrm{km}\mathrm{s}^1`$. The visible extent of the remnant has properties which are almost identical to those of Model I, and we find its velocity dispersion to be $`11.1\mathrm{km}\mathrm{s}^1`$. Both results are again in good agreement with the observations.

The two initial satellites (Models I and II) have the same stellar mass distributions in their inner regions, differing only in that one has an extended dark halo. We may thus conclude that the presence of a dark halo does not affect the final structure of the remnant, which is very similar in both models. However there is a significant difference in the properties of their debris streams. In Model I the unbound debris streams are predicted to contain 5.2 times the light in the main body of the remnant ($`M_V14.1`$), as defined by the dotted contour in Figure 3, whereas in Model II ($`M_V13.4`$) this ratio is 4.85. If we had chosen Model II to be a constant mass-to-light ratio model, we would have got an almost equally good fit to the main body of Sagittarius, but would have predicted the streams to contain 19 times the light in the main body of the remnant. In this last case, Sagittarius would have contributed $`4.56\times 10^8\mathrm{L}_{}`$ to the Galactic stellar halo in the form of debris stars (for $`\mathrm{{\rm Y}}=3.5`$). Thus we see that the observed properties of the main remnant do not usefully constrain the number of stars that may be present in the debris streams, but that the different models can be better constrained from the properties of their debris streams, as we exemplify below.

### 3.3 Discussion

#### 3.3.1 Some predictions

In this section we concentrate for simplicity on Model I. We can use it to predict star counts as a function of distance and radial velocity at different points on the sky. We focus on fields along the path defined by the orbit of Sgr, which is where we expect to find debris streams. This is illustrated in Figure 7, where the number counts are normalized to their values on the main body of our simulated Sagittarius, as shown in the first row. We assume fields which are $`1^{}\times 1^{}`$. For the distance, we use 5 kpc bins, whereas for the radial velocity we take 25 $`\mathrm{km}\mathrm{s}^1`$ bins. Note that the contrast of structures in the radial velocity counts are generally larger than in the distance counts, indicating that it should be easier to detect streams in velocity space rather than as density inhomogeneities (see also Helmi & White 1999). This is particularly true considering the much greater relative precision of the velocity measurements. Space density enhancements often occur near the orbital turning points; several are seen as sharp features in the central panel of Fig. 8.

Our model can also be used to predict where streams originating in different mass loss events should be found. This is illustrated in Figure 8 where different colours indicate material lost at different pericentric passages. Note that since the surface brightness of the unbound material decreases with time, material lost in early passages is considerably more difficult to detect than recent mass loss (for an axisymmetric potential the time dependence is $`1/t^2`$, but if the potential may be considered as nearly spherical the surface density will effectively decrease as $`1/t`$; see Helmi & White 1999). The central panel (latitude vs. heliocentric distance) explains why Sagittarius streams have been more difficult to detect above the Galactic plane than below it, even though the density contrast is higher for the northern streams (as shown in the second and third panels of Fig. 7). From the left panel,$`90^{}l90^{}`$, we see that the stream of stars lost in the previous pericentric passage (shown in blue) becomes more distant as we go north. For example, at $`b=40^{}`$, the stream is located approximately 50 kpc from the Sun. The red giant clump visual magnitude at this distance would be roughly 19.3<sup>m</sup>, compared to the 17.85<sup>m</sup> observed in the main body of Sagittarius.

#### 3.3.2 Comparison to data outside the main body of Sagittarius

Even though we have constructed our models to reproduce the properties of the main body of Sagittarius, it is nevertheless worthwhile to compare our simulations to data sets which have claimed detections of Sagittarius debris.

##### Outer Structure of Sagittarius.

Mateo et al. (1998) have traced Sagittarius material out to 30 degrees from its nucleus: the globular cluster M54. They obtained deep photometric data along the southeast extension of the major axis of Sagittarius. In Figure 9 we show the particle counts in our simulation for the strip 3 to 10 in longitude, and spanning about 30 in latitude outside the main remnant body. For comparison we plot the data by Mateo and collaborators, shifted a few degrees in latitude, and arbitrarily offset in number counts. Thus qualitatively we reproduce the break in the number counts profile. This change in slope is indicative of the transition between material which is still bound today and that lost in the last pericentric passage.

##### Star counts at $`b=40^{}`$.

Majewski et al. (1999) have claimed a detection of a possible stream from Sagittarius at $`b=40^{}`$ and $`l=11^{}`$, at a slightly smaller heliocentric distance of $`23`$ kpc and with a radial velocity of the order of $`30\mathrm{km}\mathrm{s}^1`$. As they discuss, this velocity may be strongly affected by contamination by other Galactic components. We note, however, that we would predict a stream of stars (shown in blue) going through this latitude and longitude with roughly the observed distance, and with a radial velocity of $`55\mathrm{km}\mathrm{s}^1`$. (See the central and bottom left panels of Fig. 8, $`90^{}l90^{}`$.). As mentioned above, this stream is formed mostly by material lost in the previous pericentric passage and not three passages ago, as in the model of Johnston et al. (1999). This difference reflects the different orbital timescales in the two models. The surface density of stars may be able to distinguish between them; it is predicted to be higher in our case.

Unfortunately, Majewski and collaborators could not detect the northern stream. They either did not reach the magnitude limit of $`19.3^m`$ expected for the red giant clump, or were offset by a few degrees from its expected location. Thus, for example, Majewski et al. (1999) had a limiting magnitude of $`21`$ at $`b=41^{}`$ and $`l=6^{}`$, but $`V19`$ at $`b=41^{}`$ and $`l=6^{}`$. The actual stream in our model is predicted to go through $`l1^{}`$ and to be about $`2^{}`$ wide. Note that the width prediction is more secure than the location since the motion of Sagittarius in the $`l`$-direction is poorly constrained at present, although a flattened halo would make the streams wider.

##### RR Lyrae found by the Sloan Digital Sky Survey.

The Sloan Digital Sky Survey (SDSS) commissioning data has detected 148 candidate RR Lyrae stars in about 100 deg<sup>2</sup> of sky, along the celestial equator ($`1.27^{}\delta 1.27^{}`$), and from $`\alpha =160.5^{}`$ to $`\alpha =236.5^{}`$ (Ivezic et al., 2000). Although the faint-magnitude limit of the SDSS would allow them to detect RR Lyrae stars to large Galactocentric distances, they find no candidates fainter than r\*$``$ 20, i.e., farther than 65 kpc from the Galactic center. The distribution of stars in their sample is very inhomogeneous and shows a clump of over 50 stars at about 45 kpc from the Galactic centre, which is also detected in the distribution of nonvariable objects with RR Lyrae star colors.

By studying carefully Figure 8, and from our previous discussions, we are naturally led to believe this substructure could be associated with the northern streams of Sagittarius. In the upper left panel of Figure 10 we see how, in our simulations of Model I, a stream of material intersects the area observed by SDSS. The positions of the particles in our simulations are in excellent agreement with those of the RR Lyrae candidates belonging to the reported substructure. The upper right panel shows the visual magnitude of the particles falling in the region of the sky analysed by SDSS. We note here that there are basically two substructures in this region: one at $`V19.5^m`$, and a second one, at a fainter magnitude $`V20.5^m`$ (for $`M_V=0.7^m`$ characteristic of RR Lyrae stars, e.g. Layden et al. (1996)). The first lump clearly could correspond to the substructure observed in the SDSS data. The material in this lump is mostly formed by particles that were lost in recent pericentric passages (i.e. 1 – 3 Gyr ago) as shown in the bottom left panel of Figure 10.

As Ivezic et al. (2000) discuss, they do not find any RR Lyrae stars fainter than $`V20^m`$. This would be in apparent contradiction with our results, (e.g. top right panel of Fig. 10). However, we need to estimate how much material we find in each lump, calibrate this number with respect to the number of RR Lyrae in the lump observed by SDSS, and thereby determine how many RR Lyrae SDSS could have missed. In the first lump we find 1264 particles, whereas the second has 362 particles. According to Ivezic et al. (2000) the detection efficiency decreases rapidly between $`V20^m`$, where it is fifty per cent, and $`V21^m`$ where it is zero. We here assume that for stars of $`V20.5^m`$ this efficiency is about 15%, which means that only 54 of the 362 particles could, in principle, have been observed. Therefore, we estimate that the ratio of unobserved to that of observed debris material is 0.043 in this region of the sky. Thus if SDSS found $`50`$ RR Lyrae belonging to the first substructure, it should have detected $`2.14\pm 1.46`$ RR Lyrae in the fainter magnitude range. This means that the failure to detect fainter RR Lyrae in this region of the sky is barely significant in this context. From this perspective we cannot rule out that a second stream of debris material is located at much larger distances (typically between 80 and 100 kpc from the Sun, as shown in Figure 8).

Nevertheless the absence of a visible stream may be indicating that this material could be dark-matter dominated. This second stream is formed by particles that became unbound more than 7 Gyr ago. It therefore corresponds to particles orbiting the outskirts of the progenitor of Sagittarius. If this region of the system was dark-matter dominated, such streams would remain unobservable. Fainter data ($`V2021^m`$) in this region of the sky could be crucial to constrain the initial properties of the system, e.g. size, total luminosity. This particular region of the sky should thus be explored further!

##### Carbon stars by the APM

The APM survey has detected about 75 high latitude carbon giants presumably belonging to the halo. These stars being of intermediate age, could trace streams that have recently become unbound from Sagittarius or from other Galactic satellites. Ibata et al. (2000) have proposed that a large fraction of the observed halo carbon stars belong to Sagittarius tidal debris, since they preferentially occur near the great circle of its orbit. Even though there are large uncertainties in the determination of distances to these carbon stars, and the survey is not complete, particularly in regions where we expect Sagittarius streams to be present, this proposal clearly fits within the expectations for the models we have developed here.

## 4 Conclusions

We have found viable models for the Sagittarius dwarf galaxy with a wide range of total luminosities and masses, and both with and without extended dark halos. A purely stellar progenitor could be a King model with a total velocity dispersion of 18.9 $`\mathrm{km}\mathrm{s}^1`$, a core radius of 0.44 kpc and a tidal radius of 3 kpc. For the case where the progenitor is embedded in an extended massive halo, the initial stellar distribution follows a King profile with the same core radius, a slightly larger total velocity dispersion of $`23\mathrm{km}\mathrm{s}^1`$ and similar extent. The dark-matter is more extended. The data available at present only weakly constrain the total initial extent either of the light or of the mass. The observed metallicity data, for example, are consistent with an initial galaxy similar to either of our detailed models, both of which would lie within the scatter of the luminosity–size–velocity dispersion–metallicity distribution for more distant dwarf spheroidal galaxies in the Local Group. Thus we see no indication that Sagittarius is in any way anomalous. Further work on the debris streams of Sagittarius is needed to constrain better its initial total luminosity, and to distinguish between purely stellar or dark-matter dominated progenitors.

It is certainly encouraging that our models could reproduce the data available both on the main body and on the debris streams. We wish to stress however, that this does not mean that we have found the “ultimate” model. Other models with similar characteristics may also exist. Alternatives would include progenitors with smaller stellar masses or larger dark halos; flattened systems or with anisotropic velocity distributions; or systems with a stellar disk and a spherical dark halo (as proposed for the progenitors of dSph by Mayer et al. (2000)). Moreover, our assumption of a rigid Galactic potential, which does not vary in time over 12 Gyr, is clearly simplistic in view of current models for the formation of structure in the Universe. Only when we have a better estimate of the total luminosity of Sagittarius, both in its main body, as well as on its streams, we will be able to model it in greater detail. The present interest in the debris streams of Sagittarius will help us understand not only the properties of what has turned out to be just another dwarf spheroidal, but also the formation history of our Galaxy. A complete map of the streams will, for example, allow us to derive the Galactic potential (Johnston et al. 1999). If these streams are less smooth or broader than expected, this may indicate smaller scale structure present in the halo either now or when this was assembled.

## Acknowledgments

We thank Pavel Kroupa and James Binney for comments on earlier versions of this manuscript. We have enjoyed discussions with Heather Morrison and Paul Harding. CONICET, Fundación Antorchas, DAAD–Fundación Antorchas, and EARA are acknowledged for financial support.


text:
# Critical behavior of n-vector model with quenched randomness

## 1 Introduction

The phase transitions and critical phenomena is one of the most widely investigated topic in modern physics. Nevertheless, an eliminated number of exact and rigorous results are available, and they refer mainly to the two–dimensional systems  and fractals . Rigorous results have been obtained also in four dimensions  based on an exact renormalization group (RG) technique . The RG method, obviously, provides exact results at $`d>4`$ (where $`d`$ is the spatial dimensionality), but this case is somewhat trivial in view of critical phenomena. In three dimensions, approximate methods are usually used based on perturbation theory.

Here we present a particular result obtained within the diagrammatic perturbation theory. The Ginzburg–Landau phase transition model with $`O(n)`$ symmetry (i.e., the $`n`$–vector model) is considered, which includes a quenched random temperature disorder. The usual prediction of the perturbative RG field theory  is that, in the case of the spatial dimensionality $`d<4`$ and small enough $`n`$ (at $`n=1`$ and $`n0`$, in particular), the critical behavior of the $`n`$–component vector model is essentially changed by the quenched randomness. Here we challenge this conventional point of view based on a mathematical proof. We have proven rigorously that within the diagrammatic perturbation theory the critical exponents in the actually considered model cannot be changed by the quenched randomness at $`n0`$.

## 2 The model

We consider a model with the Ginzburg–Landau Hamiltonian

$`H/T={\displaystyle \left[\left(r_0+\sqrt{u}f(𝐱)\right)\phi ^2(𝐱)+c(\phi (𝐱))^2\right]𝑑𝐱}`$
$`+uV^1{\displaystyle \underset{i,j,𝐤_1,𝐤_2,𝐤_3}{}}\phi _i(𝐤_1)\phi _i(𝐤_2)u_{𝐤_1+𝐤_2}\phi _j(𝐤_3)\phi _j(𝐤_1𝐤_2𝐤_3)`$ (1)
which includes a random temperature (or random mass) disorder represented by the term $`\sqrt{u}f(𝐱)\phi ^2(𝐱)`$. For convenience, we call this model the random model. In Eq. (2) $`\phi (𝐱)`$ is an $`n`$–component vector (the order parameter field) with components $`\phi _i(𝐱)=V^{1/2}_𝐤\phi _i(𝐤)e^{i\mathrm{𝐤𝐱}}`$, depending on the coordinate $`𝐱`$, and $`f(𝐱)=V^{1/2}_𝐤f_𝐤e^{i\mathrm{𝐤𝐱}}`$ is a random variable with the Fourier components $`f_𝐤=V^{1/2}f(𝐱)e^{i\mathrm{𝐤𝐱}}𝑑𝐱`$. Here $`V`$ is the volume of the system, $`T`$ is the temperature, and $`\phi _i(𝐤)`$ is the Fourier transform of $`\phi _i(𝐱)`$. An upper limit of the magnitude of wave vector $`k_0`$ is fixed. It means that the only allowed configurations of the order parameter field are those corresponding to $`\phi _i(𝐤)=0`$ at $`k>k_0`$. This is the limiting case $`m\mathrm{}`$ ($`m`$ is integer) of the model where all configurations of $`\phi (𝐱)`$ are allowed, but Hamiltonian (2) is completed by term $`\underset{i,𝐤}{}(k/k_0)^{2m}\phi _i(𝐤)^2`$.

The perturbation expansions of various physical quantities in powers of the coupling constant $`u`$ are of interest. In this case $`n`$ may be considered as a continuous parameter. In particular, the case $`n0`$ has a physical meaning describing the statistics of polymers .

The system is characterized by the two–point correlation function $`G_i(𝐤)`$ defined by the equation

$$\phi _i(𝐤)\phi _j(𝐤)=\delta _{i,j}G_i(𝐤).$$
(2)
Because of the $`O(n)`$ symmetry of the considered model, we have $`G_i(𝐤)G(𝐤)`$, i.e., the index $`i`$ may be omitted. It is supposed that the averaging is performed over the $`f(𝐱)`$ configurations with a Gaussian distribution for the Fourier components $`f_𝐤`$, i.e., the configuration $`\{f_𝐤\}`$ is taken with the weight function

$$P(\{f_𝐤\})=Z_1^1\mathrm{exp}\left(\underset{𝐤}{}b(𝐤)f_𝐤^2\right),$$
(3)
where

$$Z_1=\mathrm{exp}\left(\underset{𝐤}{}b(𝐤)f_𝐤^2\right)D(f),$$
(4)
and $`b(𝐤)`$ is a positively defined function of $`k`$. Eq. (2) defines the simplest random model considered in  (according to the universality hypothesis, the factor $`\sqrt{u}`$ does not make an important difference). Our random model describes a quenched randomness since the distribution over the configurations $`\{f_𝐤\}`$ of the random variable is given (by Eqs. (3) and (4)) and does not depend neither on temperature nor the configuration $`\{\phi _i(𝐤)\}`$ of the order parameter field. More precisely, the common distribution over the configurations $`\{f_𝐤\}`$ and $`\{\phi _i(𝐤)\}`$ is given by

$$P(\{f_𝐤\},\{\phi _i(𝐤)\})=P(\{f_𝐤\})\times Z_2^1(\{f_𝐤\})\mathrm{exp}(H/T),$$
(5)
where $`Z_2(\{f_𝐤\})=\mathrm{exp}(H/T)D(\phi )`$ and $`H`$ is defined by Eq. (2). For comparison, the common distribution is the equilibrium Gibbs distribution in a case of annealed randomness sometimes considered in literature.

## 3 A basic theorem

We have proven the following theorem.

Theorem. In the limit $`n0`$, the perturbation expansion of the correlation function $`G(𝐤)`$ in $`u`$ power series for the random model with the Hamiltonian (2) is identical to the perturbation expansion for the corresponding model with the Hamiltonian

$`H/T={\displaystyle \left[r_0\phi ^2(𝐱)+c(\phi (𝐱))^2\right]𝑑𝐱}`$
$`+uV^1{\displaystyle \underset{i,j,𝐤_1,𝐤_2,𝐤_3}{}}\phi _i(𝐤_1)\phi _i(𝐤_2)\stackrel{~}{u}_{𝐤_1+𝐤_2}\phi _j(𝐤_3)\phi _j(𝐤_1𝐤_2𝐤_3)`$ (6)
where $`\stackrel{~}{u}_𝐤=u_𝐤\frac{1}{2}f_𝐤^2`$.

For convenience, we call the model without the term $`\sqrt{u}f(𝐱)\phi ^2(𝐱)`$ the pure model, since this term simulates the effect of random impurities .

Proof of the theorem. According to the rules of the diagram technique, the formal expansion for $`G(𝐤)`$ involves all connected diagrams with two fixed outer solid lines. In the case of the pure model, diagrams are constructed of the vertices  , with factor $`uV^1\stackrel{~}{u}_𝐤`$ related to any zigzag line with wave vector $`𝐤`$. The solid lines are related to the correlation function in the Gaussian approximation $`G_0(𝐤)=1/\left(2r_0+2ck^2\right)`$. Summation over the components $`\phi _i(𝐤)`$ of the vector $`\phi (𝐤)`$ yields factor $`n`$ corresponding to each closed loop of solid lines in the diagrams. According to this, the formal perturbation expansion is defined at arbitrary $`n`$. In the limit $`n0`$, all diagrams of $`G(𝐤)`$ vanish except those which do not contain the closed loops. In such a way, for the pure model we obtain the expansion

$$G(𝐤)=\text{}+\text{}+\mathrm{}.$$
(7)
In the case of the random model, the diagrams are constructed of the vertices and  . Besides, it is important that only those diagrams give a nonzero contribution where each dotted line is coupled to another dotted line. The factors $`uV^1f_𝐤^2`$ correspond to the coupled dotted lines and the factors $`uV^1u_𝐤`$ correspond to the dashed lines. Thus, we have

$$G(𝐤)=\text{}+\left[\text{}+\text{}\right]+\mathrm{}.$$
(8)
In our notation the combinatorial factor corresponding to any specific diagram is not given explicitly, but is implied in the diagram itself. In the random model, first the correlation function $`G(𝐤)`$ is calculated at a fixed $`\{f_𝐤\}`$ according to the distribution $`Z_2^1\left(\{f_𝐤\}\right)\mathrm{exp}(H/T)`$ (which corresponds to diagrams where solid lines are coupled, but the dotted lines with factors $`\sqrt{u}V^{1/2}f_𝐤`$ are not coupled), performing the averaging with the weight (3) over the configurations of the random variable (i.e., the coupling of the dotted lines) afterwards. According to this procedure, the diagrams of the random model in general (not only at $`n0`$) do not contain parts like  ,  ,  , etc. Such parts would appear after the coupling of dotted lines only if unconnected (i.e., consisting of separate parts) diagrams with fixed $`\{f_𝐤\}`$ would be considered.

Thus, in the considered random model, the term of the $`l`$–th order in the perturbation expansion of $`G(𝐤)`$ in $`u`$ power series is represented by diagrams constructed of a number $`M_1`$ of vertices and an even number $`M_2`$ of vertices (i.e., $`M_2/2`$ double–vertices ), such that $`l=M_1+M_2/2`$. In the pure model, defined by Eq. (3), this term is represented by diagrams constructed of $`l`$ vertices  . It is clear and evident from Eqs. (7) and (8) that all diagrams of the random model are obtained from those of the pure model if any of the zigzag lines is replaced either by a dashed or by a dotted line, performing summation over all such possibilities. Note that such a method is valid in the limit $`n0`$, but not in general. The problem is that, except the case $`n0`$, the diagrams of the pure model contain parts like  ,  ,  , etc. If all the depicted here zigzag lines are replaced by the dotted lines, then we obtain diagrams which are not allowed in the random model, as it has been explained before. At $`n0`$, the only problem is to determine the combinatorial factors for the diagrams obtained by the above replacements. For a diagram constructed of $`M_1`$ vertices and $`M_2`$ vertices the combinatorial factor is the number of possible different couplings of lines, corresponding to the given topological picture, divided by $`M_1!M_2!`$.

It is suitable to make some systematic grouping of the diagrams of the random model. The following consideration is valid not only for the diagrams of the two–point correlation function, but also for free energy diagrams (connected diagrams without outer lines) and for the diagrams of $`2m`$–point correlation function (i.e., the diagrams with $`2m`$ fixed outer solid lines, containing no separate parts unconnected to these lines). It is supposed that at $`n0`$ the main terms are retained, which means that the free energy diagrams contain a single loop of solid lines. We define that all diagrams which can be obtained from the $`i`$–th diagram (i.e., the diagram of the $`i`$–th topology) of the pure model, belong to the $`i`$–th group. The sum of the diagrams of the $`i`$–th group can be found by the following algorithm.

1. First, the $`i`$–th diagram of pure model is depicted in an a priori defined way.

2. Each vertex is replaced either by the vertex  , or by the double-vertex  , performing the summation over all possibilities. Besides, all vertices and and all lines are numbered before coupling, and all the distributions of the numbered vertices and lines over the numbered positions (arranged according to the given picture defined in step 1 and according to the actually considered choice, defining which of the vertices must be replaced by and which of them must be replaced by ) are counted as different. Each specific realization is summed over with the weight $`1/(M_1!M_2!)`$.

3. To ensure that each specific diagram is counted with the correct combinatorial factor, the result of summation in step 2 is divided by the number of independent symmetry transformations (including the identical transformation) $`S_i`$ for the considered $`i`$–th diagram constructed of vertices  , where the symmetry transformation of a diagram is defined as any possible redistribution (such that the outer solid lines are fixed) of vertices and coupled lines not changing the given picture. Really, the coupling of lines is not changed if any of the symmetry transformations with any of the specific diagrams of the $`i`$–th group is performed, whereas, according to the algorithm of step 2, original and transformed diagrams are counted as different.

It is convenient to modify step 2 as follows. Choose any one replacement of the vertices by and  , and perform the summation over all such possibilities. For any specific choice we consider only one of the possible $`M_1!M_2!`$ distributions of the numbered $`M_1`$ vertices and $`M_2`$ vertices over the fixed numbered positions, and make the summation with weight $`1`$ instead of the summation over $`M_1!M_2!`$ equivalent (i.e., equally contributing) distributions with the weight $`1/(M_1!M_2!)`$.

Note that the location of any vertex is defined by fixing the position of dashed line, the orientation of which is not fixed. According to this, the summation over all possible distributions of lines (numbered before coupling) for one fixed location of vertices (as consistent with the modified step 2) yields factor $`8^{M_1}4^{M_2/2}`$. The $`i`$–th diagram of the pure model also can be calculated by such an algorithm. In this case the summation over all possible line distributions yields a factor of $`8^l`$, where $`l=M_1+M_2/2`$ is the total number of vertices in the $`i`$–th diagram. Obviously, the summation of diagrams of the $`i`$–th group can be performed with factors $`8^l`$ instead of $`8^{M_1}4^{M_2/2}`$, but in this case factors $`\frac{1}{2}uV^1f_𝐤^2`$ must be related to the coupled dotted lines instead of $`uV^1f_𝐤^2`$. In this case the summation over all possibilities where zigzag lines are replaced by dashed lines with factors $`uV^1u_𝐤`$ and by dotted lines with factors $`\frac{1}{2}uV^1f_𝐤^2`$, obviously, yields a factor $`uV^1\left(u_𝐤+\frac{1}{2}f_𝐤^2\right)uV^1\stackrel{~}{u}_𝐤`$ corresponding to each zigzag line with wave vector $`𝐤`$. Thus, the sum over the diagrams of the $`i`$–th group is identical to the $`i`$–th diagram of the pure model defined by Eq. (3). By this the theorem has proved.

## 4 Remarks and conclusions

The theorem has been formulated for the two–point correlation function, but the proof, in fact, is more general, as regards the relation between diagrams of random and pure models. Thus, the statement of the theorem is true also for free energy and for $`2m`$–point correlation function.

According to the proven theorem and this remark, at $`n0`$ the considered pure and random models cannot be distinguished within the diagrammatic perturbation theory. Thus, if, in principle, critical exponents can be determined from the diagrammatic perturbation theory at $`n0`$, then, in this limit, the critical exponents for the random model are the same as for the pure model. We think that in reality correct values of critical exponents can be determined from the diagrammatic perturbation theory, therefore the quenched random temperature disorder does not change the universality class at $`n0`$. This our conclusion correlates with results of some other authors. In particular, there is a good evidence that the universality class is not changed by the quenched randomness at $`n=1`$. It has been shown by extensive Monte–Carlo simulations of two–dimensional dilute Ising models  that the critical exponent of the defect magnetization is a continuous function of the defect coupling. By analyzing the stability conditions, it has been concluded in Ref.  that the critical exponent $`\nu `$ of the bulk correlation length of the random Ising model does not depend on the dilution, i.e., $`\nu =1`$ at $`d=2`$ both for diluted and not diluted Ising models. The standard (pertubative) RG method predicts the change of the universality class by the quenched randomness. We think, this is a false prediction. The fact that the standard RG method provides incorrect result is not surprising, since it has been demonstrated (in fact, proven) in Ref.  that this method is not valid at $`d<4`$.

## References