text:
# Quantum-Critical Dynamics of the Skyrmion Lattice.

## Abstract

Near to filling fraction $`\nu =1`$, the quantum Hall ferromagnet contains multiple Skyrmion spin excitations. This multi-Skyrmion system has a tremendously rich quantum-critical structure. This is simplified when Skyrmions are pinned by disorder. We calculate the nuclear relaxation rate in this case and compare the result with experiment. We discus how such measurements may be used to further probe the quantum-critical structure of the multi-Skyrmion system.

At exact filling of a single Landau level the quantized Hall state forms an almost perfect ferromagnet. This quantum Hall ferromagnet (QHF) has some novel features due to the phenomenology of the underlying quantized Hall state. Magnetic vortices, or Skyrmions, in the QHF carry quantized electrical charge. These Skyrmions are stabilised by a chemical potential so that the ground state slightly away from filling fraction $`\nu =1`$ contains a finite density of them. This original proposal of Sondhi et al. has been confirmed in a number of experiments.

The $`T=0`$ phase diagram of the multi-Skyrmion system has been thoroughly investigated. In the absence of disorder, crystalline arrangements are expected. For filling fractions very close to $`\nu =1`$ a triangular lattice is formed with a transition to a square lattice as the deviation from $`\nu =1`$ is increased. The statistical mechanics of the possible melting transitions has been considered. At the highest Skyrmion densities, zero point fluctuations are expected to give rise to a quantum-melted state.

Despite this wealth of study, a complete account of the experimental observables has not been achieved. For example, nuclear magnetic resonance provides one of the clearest probes of the spin polarisation in the QHF. Although it is understood in general terms how low-energy spin-fluctuations of the Skyrmion system may enhance the relaxation of nuclear spins, attempts to calculate relaxation rates have been flawed. The fundamental physics, missed in other considerations, is the quantum-critical nature of fluctuations of the Skyrmion lattice. One immediate consequence of this quantum-criticallity is that the limits of temperature/frequency$`0`$ and frequency/temperature$`0`$ are very different. Typically, experimental probes are at frequencies much less than temperature and the latter limit is appropriate. This means that zero-temperature calculations cannot model experiments correctly.

Let us consider these points a little further. In our analysis below, we will find an underlying gappless XY-model governing orientational fluctuations of the multi-Skyrmion system. If a spinwave expansion is attempted for a gappless magnet at or below its critical dimension, the occupation of low-frequency modes is found to diverge. The constraint, fixing the magnitude of the local spin, restricts this divergence. Interplay between divergence and constraint gives rise to a finite (temperature dependent) correlation length, $`\xi (T)`$, beyond which correlations of the magnet decay exponentially. The dynamics of the critical magnet are very different on length scales greater or less than $`\xi `$. On lengthscales less than $`\xi `$, the groundstate is ordered (albeit in a quantum superposition of all possible orientations due to long wavelength spin fluctuations). Fluctuations with wavelength less than $`\xi (T)`$ may, therefore, be described by a modified spinwave expansion. On lengthscales greater than $`\xi `$, the groundstate is disordered and fluctuations are overdamped. This quantum relaxational dynamics is a striking feature of quantum-critical systems and leads to interesting universalities. Regimes of renormalized classical and quantum activated behaviour at low-temperature, cross over to universal behaviour in the high-temperature quantum-critical regime.

The Skyrmion spin-configuration consists of a vortex-like arrangement of in-plane components of spin with the z-component reversed in the centre of the Skyrmion and gradually increasing to match the ferromagnetic background at infinity. At large distances, the spin distribution decays exponentially to the ferromagnetic background on a length scale determined by the ratio of spin stiffness to Zeeman energy. An individual Skyrmion may be characterized completely by its position (i.e. the point at which the spin points in the opposite direction to the ferromagnetic background) its size (i.e. the number of flipped spins) and the orientation of the in-plane components of spin. The equilibrium size of the Skyrmion is determined by a balance between its coulomb and Zeeman energies (In the presence of a disorder potential, the potential energy of the Skyrmion also enters this balancing act).

Consider a ferromagnet with a dilute distribution of Skyrmions. The normal modes of this system are relatively easy to identify. Firstly, ferromagnetic spinwaves propagate in-between the Skyrmions. The spectrum of these is gapped by the Zeeman energy and will be ignored from now on. Positional fluctuations, or phonon modes, of the Skyrmions are gapless in a pure system, but gapped when the lattice is pinned by disorder. Finally, fluctuations in the in-plane orientation and size must be considered. These two types of fluctuation are intimately connected; rotating a Skyrmion changes its size. This follows from the commutation relations of quantum angular momentum operators.

The orientation, $`\theta (𝐱_i,t)`$, of a Skyrmion centred at a point $`𝐱_i`$ is described by the following effective action:

$$S=\frac{1}{2}𝑑t\left[\underset{i}{}I_i\theta (𝐱_i,t)_t^2\theta (𝐱_i,t)\underset{<i,j>}{}J_{ij}\mathrm{cos}(\theta (𝐱_i,t)\theta (𝐱_j,t))\right].$$
(1)
$`I_i`$ is the moment of inertia of the $`i^{th}`$ Skyrmion and $`J_{ij}`$ is the stiffness to relative rotations of neighbouring Skyrmions. The first term in Eq.(1) arises due to the change in energy of a Skyrmion when its size fluctuates; $`\mathrm{\Delta }E=I^1\delta s^2/8`$. Since the z-component of spin and orientation are conjugate coordinates, a cross-term $`i\delta s(𝐱_i,t)_t\theta (𝐱_i,t)/2`$ appears in their joint effective action. Integrating out $`\delta s(𝐱_i,t)`$ gives the first term in Eq.(1). Clearly, $`1/4I`$ is the second derivative of the Skyrmion energy with respect to its spin. $`I`$ is related to the Skyrmion size (in fact $`I=24\mu _BgB/s`$) and, in the absence of disorder, is the same for all Skyrmions. Correlation functions involving $`\delta s`$ may be calculated using Eq.(1) by making the replacement $`\delta s(𝐱_i,t)2iI_t\theta (𝐱,t)`$\- the result of a simple Gaussian integration over $`\delta s`$. The second term in Eq.(1) is an effective dipole-interaction of Skyrmions due to the energetics of overlapping Skyrmion tails. In a square Skyrmion lattice, $`J_{ij}`$ is independent of lattice site. A continuum limit may be taken where $`\theta _i`$ is replaced by a staggered field, $`\theta _i\theta _i+\eta _i\pi `$, with $`\eta _i=0,1`$ on adjacent sites, and $`\theta _i\theta _j`$ is replaced by a derivative;

$$S=\frac{1}{2}\frac{d\omega d^2k}{(2\pi )^3}\theta (𝐤,\omega )\left(\delta \nu \overline{\rho }I\omega ^2+J𝐤^2\right)\theta (𝐤,\omega ).$$
(2)
The frequency integral in this expression is shorthand for a Matsubara summation at finite temperature and the momentum integral is over the Brillouin zone. A factor of the Skyrmion density $`\delta \nu \overline{\rho }`$ has been introduced, where $`\delta \nu `$ is the deviation from filling fraction $`\nu =1`$ and $`\overline{\rho }`$ is the electron density. There are a few caveats to the use of Eq.(2). We defer discusion of these until later. This model is perhaps most familiar as an effective theory of the Josephson junction array. In this case, $`\theta `$ is the phase of the superconducting order parameter and its conjugate coordinate is the charge of the superconducting junction.

In order to calculate properties of the Skyrmion lattice, we must relate fluctuations in the Skyrmion orientation to fluctuations in the orientation of local spin. We use a coherent-state representation of the polarization of the local spin, via an O(3)-vector field $`𝐧(𝐱,t)`$. The static spin distribution at a point $`𝐱`$ relative to the centre of a single Skyrmion is denoted by $`𝐧(𝐱)`$ and its in-plane components by $`n_x+in_y=n_re^{i\varphi _0}`$. The in-plane components of local spin at a point $`𝐱`$, in response to rotational fluctuations of a Skyrmion centred at $`𝐱_i`$, are given by

$$n_re^{i\theta }(𝐱,t)=()^{\eta _i}\left[n_r(𝐱𝐱_i)+2iI\frac{n_r(𝐱𝐱_i)}{s}_t\theta (𝐱_i,t)\right]e^{i\theta +i\varphi _0}(𝐱_i,t).$$
(3)
We have used the conjugate relationship between Skyrmion spin and orientation in writing down this expression. In a distribution of many Skyrmions, one must in principle sum the contributions of all Skyrmions to the fluctuation in spin at the point $`𝐱`$. However, in the dilute limit in which we are performing our explicit calculation, Skyrmions are exponentially localized. The dominant spin fluctuations occur near to the centre of Skyrmions and so, to logarithmic accuracy, the local fluctuations at a point $`𝐱`$ are due only to the nearest Skyrmion.

We will calculate the nuclear relaxation rate due to low energy quantum fluctuations of the Skyrmion lattice;

$$\frac{1}{T_1}=T\gamma \underset{\omega 0}{lim}\frac{d^2k}{(2\pi )^2}\frac{mS_+(𝐤,\omega )S_{}(𝐤,\omega )}{\omega },$$
(4)
where $`\gamma `$ is the hyperfine coupling constant. Other physical observables, such as the temperature dependence of magnetization, $`M=𝑑td^2xn_z(𝐱,t),`$ may be calculated similarly. Our first task is to replace the expectation of the spin raising and lowering operators in Eq.(4) by correlators of the Skyrmion orientation. Substituting from Eq.(3) into Eq.(4) and ignoring terms higher order in frequency, the nuclear relaxation rate at a point $`𝐱`$ is

$$\frac{1}{T_1}(𝐱)=T\gamma \underset{\omega 0}{lim}n_r^2(𝐱)\frac{d^2k}{(2\pi )^2}\frac{me^{i\theta }(𝐤,\omega )e^{i\theta }(𝐤,\omega )}{\omega }.$$
(5)
This takes the form of a correlation function of the Skyrmion orientation, multiplied by a profile function characteristic of the Skyrmion groundstate. The average rate is given by integrating this over an area containing a single Skyrmion and multiplying by the Skyrmion density $`\delta \nu \overline{\rho }`$. The result is identical to Eq.(5) with the replacement $`n_r^2(𝐱)\delta \nu \overline{n_r^2}`$. $`\overline{n_r^2}=\overline{\rho }d^2xn_r^2(𝐱)`$ is a number characteristic of a single Skyrmion. For a pure Skyrmion spin distribution, $`\overline{n_r^2}=2s`$ to logarithmic accuracy in the Skyrmion spin, $`s=\overline{\rho }d^2x(1n_z(𝐱))`$. Notice that radial fluctuations of local spin contribute only to higher order in frequency and have been neglected in writing down Eq.(5).

The problem of finding the nuclear relaxation rate has now been reduced to evaluating the correlation function in Eq.(5) using the effective action Eq.(2). This is rather tricky. The O(2)-quantum rotor, described by Eq.(2), is quantum-critical. It is necessary to employ a non-perturbative scheme, such as 1/N or epsilon expansions, to calculate in the quantum-critical regime of this model. Here we use the result of the 1/N expansion of Chubukov et al.

An important feature of the effective action, Eq.(2), is that it displays a zero-temperature phase transition. For $`IJ<1`$, the Skyrmion moment of inertia is sufficiently small that quantum fluctuations destroy long range order even at zero temperatures. For $`IJ>1`$, the $`T=0`$ groundstate has an infinite correlation length and is ordered. Notice that arbitrarily small temperatures destroy this long-range order even when $`IJ>1`$. This arises due to the interplay between fluctuations and constraint and can be seen in a simple mean-field calculation. In the O(2) representation of Eq.(2), $`S=𝑑td^2x\left[𝐧(I_t^2+J_𝐱^2)𝐧+\lambda (𝐱,t)(𝐧^21)\right]`$, where $`𝐧(𝐱,t)`$ is an O(2)-vector field and $`\lambda (𝐱,t)`$ is an auxiliary field that imposes the constraint $`𝐧^2=1`$. Imposing the constraint at mean-field level, $`𝐧^2=1`$, determines a temperature dependent gap, $`\lambda (T)`$. The spin correlations decay exponentially on a length scale $`\xi (T)=\sqrt{J/\lambda (T)}`$. The results of such a calculation are sketched in Fig.1.

Above a temperature of about $`T_{QC}=|2\pi JE_{\text{max}}|`$ (The cut-off, $`E_{\text{max}}=2\pi \sqrt{J/I}`$ corresponds to fluctuations with momentum at the Brillouin zone boundary, $`k_{\text{max}}=\pi \sqrt{\overline{\rho }\delta \nu }`$), the gap/correlation length develops a universal temperature dependence. In this region, thermal and quantum fluctuations are of similar importance and are very difficult to disentangle. This crossover to universal high-temperature behaviour from distinct low-temperature behaviours is a feature of all correlation functions of Eq.(2). It is usual to summarise this behaviour by the phase diagram sketched in Fig.2. In this figure, $`J(T)`$ is the renormalized stiffness in the ordered phase and $`\lambda (T)`$ is the gap in the paramagnetic phase.

The correlation function required in Eq.(5) has been calculated in Ref. by means of a 1/N expansion, with the result

$`\underset{\omega 0}{lim}{\displaystyle \frac{1}{\omega }}{\displaystyle \frac{d^2k}{(2\pi )^2}me^{i\theta }(𝐤,\omega )e^{i\theta }(𝐤,\omega )}`$
$`\begin{array}{cccc}=\hfill & \frac{1}{T}\frac{0.015}{\sqrt{IJ}}\left(\frac{kT}{\pi J}\right)^\eta ,\hfill & TT_{QC}\hfill & \\ =\hfill & \frac{1}{T}\frac{0.085}{\sqrt{IJ}}\frac{IT^2}{\lambda (0)}e^{2\sqrt{\lambda (0)/IT^2}},\hfill & TT_{QC}\text{,}\hfill & IJ<4,\hfill \\ =\hfill & \frac{0.18}{T}\left(\frac{kT}{2\pi \lambda (T)}\right)^{1/2},\hfill & TT_{QC}\text{,}\hfill & IJ>4,\hfill \end{array}`$
where $`\eta `$ is a number close to zero. Substituting these results into Eq.(5) we obtain

$$\frac{1}{T_1}=\gamma \frac{0.03s}{\sqrt{IJ}}\left(\frac{kT}{\pi J}\right)^\eta $$
(7)
at high temperature. The full behaviour is sketched in Fig.3. The kink at $`T_{KT}`$ is due to the discontinuous change in spin stiffness seen at the Kosterlitz-Thouless transition. This effect is not seen in the 1/N expansion of Ref. and must be calculated by some other means.

The nuclear relaxation rate obtained here is very different from that obtained in Ref.. As pointed out in the introduction, this is due to the unphysical limit $`T/\omega 0`$ used in Ref.. Nevertheless, it is instructive to see how the results of Ref. relate to the present formalism. Since $`T\omega `$ in their work, Côté et al consider fluctuations on length scales much less than the correlation length. The groundstate is ordered and a spinwave expansion may be used. Long wavelength fluctuations lead to a superposition of orientations- rotational averaging- the immediate consequence of which is that $`e^{i\theta }(𝐤,\omega )e^{i\theta }(𝐤,\omega )=0`$. Returning to the substitution of Eq.(3) into Eq.(4), we must retain terms to next order in frequency. This is a cross-term between radial and transverse fluctuations and involves a correlator $`_t\theta e^{i\theta }(𝐤,\omega )e^{i\theta }(𝐤,\omega )`$. Evaluating this correlation function via a $`T=0`$ spinwave expansion of the effective action, Eq.(2), reproduces the result of Ref. (up to a numerical factor due to our estimate of the Skyrmion profile function, $`n_rn_r/s`$). The zero-temperature phonon contribution may be calculated similarly.

Up to now, we have assumed that the Skyrmion lattice is pinned by disorder and that phonons may be ignored as a consequence. The situation is rather subtle and a fuller discussion is appropriate at this juncture. A quadratic effective action for phonons of the Skyrmion lattice is known. It is identical to that of an electronic Wigner crystal in a magnetic field with a vanishing effective mass. Under this effective action, Skyrmions move in small ellipses with a frequency $`\omega _𝐤|𝐤|^{3/2}/B`$ and major axes orientated transverse to the phonon momentum.

At finite temperature, the occupation of transverse phonons is infra-red divergent. This divergence is restricted by interactions between phonons arising from the non-harmonicity of the Skyrmion interaction. Unlike fluctuations in orientation, where the spinwave interaction is due to a topological constraint, these phonon interactions are non-universal (at best the universality is hidden in the details of the groundstate spin distribution and effective Skyrmion interaction potential). The resulting physics is very similar to that discussed for the rotation mode above; the phonon system is quantum-critical and has low-temperature ordered and quantum-melted phases and a high-temperature quantum-critical regime.

Even this is not the full story. Although the low energy dispersions of phonons and Skyrmion rotations are independent, non-linear interactions exist between these modes. The lattice stiffness is, in part, due to the dipole interaction of Skyrmions and is affected by fluctuations in orientation. Similarly, the dipole interaction between Skyrmions is strongly dependent upon the separation of Skyrmions and is affected by phonons. These non-linearities occur on the same footing as the phonon-phonon interactions and rotation-rotation interactions. The full quantum-critical structure of the multi-Skyrmion system is tremendously complicated.

The position taken here in neglecting this wealth of structure is that the Skyrmion lattice is pinned by disorder and the phonon spectrum gapped. Phononic fluctuations are suppressed at low temperatures and the associated critical structure occurs at higher temperature. The residual effect of phonons is a slight thermal renormalization of the rotational stiffness. The slight distortion in static positions of Skyrmions, in response to the disorder potential, gives a small random contribution to the stiffness $`J`$. This randomness produces a small region of Bose-glass phase at low temperatures, intervening between the paramagnet and renormalized classical regimes. For weak disorder, this phase only affects the physics very close to the critical point and does not affect our conclusions.

We now turn to a discussion of the experimental implication of the above calculations. Detailed measurements of $`1/T_1`$ have been carried out by Bayot et al. Above $`40mK`$, $`T_1`$ is independent of temperature. This is consistent with the rotational degrees of freedom being in their quantum-critical regime. Values of $`I`$ and $`J`$ for this system extracted from TDHFA calculations put $`IJ`$ very close to $`1`$. The system is close to criticality and the crossover to the quantum-critical regime occurs at correspondingly low temperature. At $`40mK`$ there is an abrupt step in $`1/T_1`$(and attendant peak in heat capacity). This is consistent with a Kosterlitz-Thouless transition in the orientational order (notice that the crossover temperature $`|2\pi JE_{\text{max}}|`$ may be much less than $`T_{KT}=2\pi J`$ and so the behaviour may be quantum-critical either side of the transisiton). There are a number of other candidate transitions, however, and it is not easy to discriminate between them. Considerations along the lines of those presented here allow some elaboration, but this is necessarily rather speculative and we refrain from its discussion at present. We may make some firm predictions for nuclear relaxation measurements below $`40mK`$. Changing the deviation in filling fraction or using tilted filed measurements both change the parameter $`IJ`$ and allow exploration of the phase diagram shown in Fig.1. The divergence or otherwise of the nuclear relaxation rate as temperature is reduced to zero should give a clear indication of the quantum-critical structure.

I would like to thank N. R. Cooper, J. R. Chalker, S. M. Girvin and N. Read for enlightening discusions, comments and suggestions. This work was supported by Trinity College Cambridge.


format: text

text:
# Tidal Disruption of a Solar Type Star by a Super-Massive Black Hole

## 1 Introduction

It has long been suggested that supermassive black holes (with a black hole mass $`M_{\mathrm{bh}}10^6`$ M) in relatively low-luminosity active galactic nuclei (AGNs) can be fed by tidal disruptions of stars that are on nearly radial orbits (e.g. Frank, 1978; Lacy et al., 1982; Carter & Luminet, 1985; Rees, 1988, 1990). The frequency of such events is expected to be of the order of $`10^4`$ yr<sup>-1</sup> in a galaxy like M31. In AGNs with more massive central black holes ($`M_{\mathrm{bh}}10^8`$ M), stars are “swallowed” whole, since the ratio of the event horizon radius to the tidal disruption radius increases with the black hole mass $`\left(M_{\mathrm{bh}}^{2/3}\right)`$.

Although direct observational evidence for the disruption process is lacking, a few observations have been tentatively identified as being associated with disruption events. For example, an ultraviolet flare at the center of the elliptical galaxy NGC 4552, although rather weak, has been suggested to result from such an event (Renzini et al., 1995). Other sudden eruptions that have been interpreted as potentially resulting from tidal disruptions include outbursts in IC 3599, NGC 5905, RXJ 1242.6-1119, RXJ 1624.9+7554, and RXJ 1331.9-3243 (e.g. Komossa & Greiner, 1999; Komossa & Bade, 1999; Grupe et al., 1999, and references therein). It has also been suggested that the sudden appearance and variability of the double-peaked Balmer lines in the active nucleus NGC 1097 is at least broadly consistent with being produced in a ring resulting from tidal disruption (Storchi-Bergmann et al., 1995). Similarly, it has been proposed that an outburst observed in the Seyfert galaxy NGC 5548 was due to a single star falling into a $`10^7`$ M black hole (Peterson & Ferland 1986; although other interpretations exist, e.g. Terlevich & Melnick 1988; Kallman & Elitzur 1988). More recent work on the broad emission line region of this galaxy, although interesting in its own right, did not shed any new light on the cause of the variability (Goad & Koratkar, 1998).

Several aspects of the problem of stellar encounters with supermassive black holes and of tidal disruption have been examined both analytically and numerically (e.g. Nolthenius & Katz, 1982; Bicknell & Gingold, 1983; Carter & Luminet, 1983; Hills, 1988; Kochanek, 1992, 1994; Syer & Clarke, 1992, 1993; Laguna et al., 1993; Marck et al., 1996). In addition, the interaction of a white dwarf with a massive black hole has been studied both generally (e.g. Frolov et al., 1994) and in the context of gamma-ray bursts (e.g. Fryer et al., 1999).

In the present work, we follow the evolution of a star that is tidally disrupted. In particular, we calculate the properties of the disruption debris to longer times than in some of the previous studies, using a Post-Newtonian, Smooth-Particle-Hydrodynamics (SPH) code.

The numerical method is described in $`\mathrm{\S }2`$, the results are presented in $`\mathrm{\S }3`$, and a discussion and conclusions follow.

## 2 The Numerical Method

We first introduce some notation that will be used throughout the paper. The mass of the black hole is denoted by $`M_{\mathrm{bh}}`$. The gravitational radius of the black hole is $`R_g=2GM_{\mathrm{bh}}/c^2`$. To measure the strength of the tidal encounter we use the dimensionless parameter (e.g. Press & Teukolsky, 1977)

$$\eta _t=\left(\frac{M_{}}{M_{\mathrm{bh}}}\frac{R_p^3}{R_{}^3}\right)^{1/2},$$
(1)
and to measure the magnitude of the relativistic effects we use

$$\eta _r=\frac{R_p}{R_g},$$
(2)
where $`M_{}`$ and $`R_{}`$ are the mass and radius of the star respectively and $`R_p`$ is the radius at the pericenter.

Simulating the full evolution of a star that is tidally disrupted by a black hole is a challenging numerical task. Once passing the pericenter, the star is tidally disrupted into a very long and dilute gas stream. As noted already by Kochanek (1994), each fluid element in this stream follows an almost test particle orbit. The fraction of the gas that is bound to the black hole eventually returns to pericenter and moves on to start another orbit. Relativistic orbital precession can cause this outgoing gas stream to collide with the incoming stream. Clearly, this problem consists of highly varying spatial configurations of matter, with much “empty space.” These attributes led to our selection of SPH as the numerical scheme used to tackle this problem. Indeed, previous works that used grid based codes (e.g. Khokhlov et al. (1993b, a); Frolov et al. (1994); Diener et al. (1997)) only followed the star close to the pericenter. This is also true of previous works by Evans & Kochanek (1989), and Laguna et al. (1993) that used SPH.

We use the (0+1+2.5) Post-Newtonian (PN) SPH code described in Ayal et al. (2000). This code implements the formalism of Blanchet, Damour & Schäfer (1990; hereafter BDS) and features full Newtonian gravity and hydrodynamics (0PN), the first order effects of general relativity on the gravity and hydrodynamics \[known as the first post-Newtonian, (1PN) approximation\], and gravitational wave damping (2.5PN). The independent matter variables used in this formalism consist of the following set: $`\rho _{}`$ the coordinate rest mass density, $`\epsilon _{}`$ the coordinate specific internal energy and $`𝐰`$ the specific linear momentum. In fully relativistic terms these are defined as:

$`\rho _{}`$ $`=`$ $`\sqrt{g}u^0\rho ,`$ (3)
$`\epsilon _{}`$ $`=`$ $`\epsilon (\rho _{}),`$ (4)
$`w_i`$ $`=`$ $`\left(c^2+\epsilon +p/\rho \right){\displaystyle \frac{u_i}{c}},`$ (5)
where $`\rho `$ is the rest mass density, $`\epsilon (\rho )`$ is the specific energy, $`p(\epsilon ,\rho )`$ is the pressure and $`u^\mu `$ is the four-velocity (Greek indices run from 0 to 4, latin indices from 1 to 3). The corresponding BDS variables are the above quantities neglecting all terms except 0PN, 1PN and 2.5PN. Using these variables the formalism yields an evolution system which consists of 9 Poisson equations and 4 hyperbolic equations which we solve as explained in Ayal et al. (2000).

We model the black hole using a massive point particle that has no hydrodynamical interactions. In order to be consistent with the 1PN approximation we must ensure that all the relativistic effects are small (of the order of 10%). This excludes simulating the strong field regions near the black hole’s horizon. We enforce this limit using a fully absorbing boundary condition at $`10R_g`$. Every particle crossing into this region is assumed to be accreted and taken out of the simulation. Since the mass of the star is negligible compared the mass of the black hole, we do not need to increase the mass of the black hole for every particle crossing the $`10R_g`$ boundary.

The statistical nature of SPH does not handle well single, separate particles. Indeed, the entire formalism is built on the assumption that each particle interacts hydrodynamically with about 60 (in 3D) other particles at all times. This number of interactions is needed to form a good sample of the fluid properties for an accurate calculation of gradients. In order to maintain this number of interactions under varying conditions, most SPH codes employ adaptive smoothing lengths (e.g. Benz, 1990; Monaghan, 1992), in effect changing the resolution at each point. In our problem, we have widely varying length scales when the gas stream expands, where varying the smoothing length is advantageous. On the other hand, the difference in the particle eccentricities causes them to separate when approaching pericenter, and they arrive there almost one by one. At this stage, maintaining hydrodynamic interaction with 60 other particles would require huge smoothing lengths, where a small change in the smoothing length would lead to a large change in the number of interactions. This, coupled with an algorithm that tries to maintain a fixed number of interactions, can cause large oscillations in the smoothing length, which in turn introduce a highly varying number of hydrodynamic interactions. These huge smoothing lengths for particles approaching pericenter also mean that we have a very coarse resolution at a crucial stage.

In order to overcome this numerical difficulty we introduce a “particle splitting” (PS) scheme. Whenever a particle satisfies some splitting criterion we split it into new particles each having a smaller mass and smoothing length so that the overall mass is conserved. Thus the PS algorithm consists of two parts—the splitting criterion and the splitting method. In the method we use, maximal splitting, we split each particle into 13 particles, giving the original particle 12 new neighbors (12 is the maximum number of spheres that can “touch” a sphere with the same radius). Another possible method is minimal splitting (while still maintaining a quasi spherical symmetry)—we split each particle into 5, adding 4 particles at the edges of a tetrahedron centered on the original particle. We found that the minimal splitting method tends to produce lumpier particle distributions, we therefore used maximal splitting in both the PS runs presented in this work. We split each particle into 13 particles, each having half the smoothing length of the original particle, spaced by one original smoothing length. This splitting method conserves the particle’s interaction radius which is twice it’s smoothing length. As a splitting criterion we used the ratio between each particle’s smoothing length and the average smoothing length. Whenever the smoothing length of a bound particle exceeds twice the average smoothing length we split this particle. We split only bound particles so as not to waste computational time on the unbound debris, in which we are not interested in the present work. In order for the unbound debris smoothing length not to dominate the average smoothing length which we use for the splitting criterion, we enforce a maximum smoothing length of twice the average smoothing length for unbound particles. This causes unbound particles of gas to have fewer and fewer hydrodynamical interactions, and their motion to be dominated by gravity, as is expected. Another computational time reducing technique that we use is to delete any particle that is both unbound and is farther away than $`2500R_g`$. The latter criterion ensures that the deletion of these particles does not affect the dynamics near pericenter.

In order to estimate the errors in using the particle splitting method, we compared the results of three runs. The first run was a conventional SPH run with a fixed number of 4295 particles, denoted by F1. The second and third runs were PS runs denoted by PS1 and PS2 and they differ only in the number of particles at the initial time. These parameters are summarized in Table 1.

## 3 Results

As initial conditions we took a $`\mathrm{\Gamma }=5/3`$ polytrope with a solar radius and mass. This star was then put at a distance of $`100R_g`$ and given the velocity of an appropriate parabolic Keplerian orbit with a pericenter at $`20R_g`$. Giving the initial conditions at such a distance ensures that any relativistic effects (at that position) are negligible so that using Newtonian expressions for various quantities (such as energy) is justified. As can be seen in Fig. 1, the star’s center of mass (CM) follows an almost relativistic orbit, which validates the use of the 1PN approximations for this problem. For this run we obtained the values of $`\eta _r18`$ and $`\eta _t0.65`$ at the pericenter. These values ensure that on one hand the 1PN approximation is valid, and on other, that the tidal interaction is sufficiently strong to lead to disruption.

We can roughly divide the disruption process into 2 qualitatively different regimes. In the first, the tidal forces dominate and the star is destroyed. This is followed, at about $`t=8`$ hours, by the post disruption regime. This latter stage consists of the gas stream phase and the accretions phase. During the gas stream phase the pressure is negligible and the particles move on almost Keplerian orbits. The accretion phase starts at $`t=15`$ days when the first particles return to pericenter.

### 3.1 Disruption

The disruption phase for a solar type star was studied, using Newtonian physics, by Evans & Kochanek (1989). Later it was studied in greater detail by Khokhlov et al. (1993b, a); and by Diener et al. (1997) who included a general relativistic treatment of the tidal potential of a Kerr black hole and Newtonian hydrodynamics for the star. This phase was also studied by Laguna et al. (1993) using Newtonian hydrodynamics on a fixed Kerr background. We compare our results for this phase with these previous studies.

The results shown are taken from run F1 which has a higher resolution at this stage. The rest-mass density contours presented in Fig. 2 show the effect of the tidal forces on the star. In Fig. 3 we show the central coordinate rest-mass density $`\rho _{}`$, the angular momentum $`|𝐉|`$, and total energy $`E_{}`$ of the star, during the disruption. The value of $`\rho _{}`$ does not increase beyond the initial value, and it falls rapidly after the disruption. Comparing to Laguna et al. (1993), we find that our encounter is somewhere between a non-relativistic encounter and a single-compression relativistic encounter, as expected. The values of $`\rho _{}`$ are also close to those in Diener et al. (1997) and Khokhlov et al. (1993a). The total angular momentum relative to the CM of the star increases by about 5 orders of magnitude at the disruption, and a more gentle increase occurs afterwards, caused by the gravitational torques acting on the debris. The total energy of the star also has a sharp increase at the disruption as the star becomes unbound, followed by another rise as the star moves away from the BH.

The differential mass distribution in specific energies $`dM/d\epsilon _{}`$ is shown in Fig. 4. We use $`\mathrm{\Delta }ϵ=GM_{\mathrm{bh}}R_{}/R_p^2`$, the change in the BH potential across the star (Lacy et al., 1982) as our energy scale. The mass distribution is almost constant as predicted by Rees (1988). A comparison with Evans & Kochanek (1989), who calculated this quantity for the Newtonian case, shows that in our calculation the width of the distribution is smaller by approximately $`0.5\mathrm{\Delta }ϵ`$ (FWHM).

In general, we find that the 1PN energy $`E_{\mathrm{pn}}`$ is conserved to better than 2%. The energy radiated by gravitational radiation during the disruption is negligible (compared to the total), amounting to $`1.6\times 10^{46}`$ erg.

### 3.2 Post Disruption

The post disruption phase begins at about $`t=8`$ hours. At this stage the gas is sufficiently dilute and cold to make the pressure very small. Consequently, the gas elements in the stream follow almost exact geodesics. A previous study on the gas stream phase, by Kochanek (1994), used the thinness of the stream to decouple the transverse properties of the stream from the variations along its length. This same thinness introduces difficulties into the numerical approach. The spherical nature of the SPH particles causes the code to overestimate the stream width. This effect can be overcome by using more particles but the increase in computational time renders this approach impractical with current hardware. As the SPH particles approach pericenter for the second time, the PS method therefore becomes essential. Without PS, the SPH particles approach pericenter almost one by one and the strong compressional effects are manifested through two- and three-particle interactions near the pericenter. This small number of interactions could reduce the reliability of the results. By using PS however, we increase the number of particles approaching pericenter, thus increasing the number of interactions and thereby ensuring that the hydrodynamical interactions remain adequate. The effects of PS can be seen in Fig. 5 where the conventional run has a very non-uniform particle density even when compared to the PS1 run (with only 1350 particles at this stage).

In Fig. 6 we show the distribution of eccentricities for the SPH particles in the gas stream phase. The mean eccentricity is about 0.994, as the star was initially marginally bound. The distribution does not change by much until the end of the stream phase, when the first particle approaches pericenter. The differential mass distribution in specific internal energy as a function of time is shown in Fig. 7. The earliest time shown is close to that of Fig. 4. As can be seen, at the gas stream phase the internal energy distribution shifts towards lower energies, while the second passage through pericenter heats the particles up. This heating is caused by the strong compression characterizing the second passage of the gas through pericenter. The compression is caused by the gas orbits converging to the space occupied by the star at the initial pericenter (e.g. Kochanek, 1994). This compression can be seen in Fig. 8 where we show the pressure at some fixed time after the first return to pericenter in the PS2 run. The pressure is shown along a path as a function of $`l`$, the length of the path. The path was chosen so as to pass through the gas stream and pericenter. As can be seen, the pressure rises sharply at negative $`l`$, corresponding to the passage through pericenter, and there is a large pressure gradient in the $`z`$ direction. We find that the bounce that follows this compression is sufficiently strong to impart a significant fraction of the gas with the escape velocity. In Fig. 9 we show the amount of the star’s mass that is bound, unbound, and accreted (we note again that accreted here means that it has come closer then $`10R_g`$ to the BH). As the star leaves the vicinity of the BH for the first time, 65% of its mass is bound. This stays constant during the gas stream phase since the pressures are small. Following the second passage through pericenter however, mass gets unbound, until, at the time we stop our simulation, 50% of the mass is unbound, 40% is bound and the remaining 10% is accreted. The relative difference between the two runs is about 10% in all of these quantities.

In Fig. 10 we show the differential mass distribution in orbital periods in the stream. This has been used in previous works (e.g. Rees, 1988; Evans & Kochanek, 1989; Laguna et al., 1993) to estimate the mass infall rate onto the BH, under the assumption that all the mass returning to pericenter is accreted. In Fig. 11 we show the estimated accretion rate according to Fig. 10 together with the actual accretion rate we get. Our calculations show that since some of the mass becomes unbound when reaching pericenter, the actual accretion rate is a factor of 3 lower than that inferred previously. The total accreted mass up to 60 days (Fig. 9) was overestimated by a factor of 4 when assuming that all returning mass is accreted. Another important consequence of the strong compression near pericenter is that the expected self-intersection of the gas stream in fact does not occur, since the debris is given a high velocity perpendicular to the orbital plane.

The inner $`10^4`$ R around the BH is composed of the infalling debris and a low density cloud with a density of $`10^{11}\mathrm{g}\mathrm{cm}^3`$ and with a temperature range of higher than $`10^6`$ K. In Figs. 12 and 13 we show the density, specific internal energy and velocity in the inner 2000 R around the BH. Even using PS is is clear that we cannot resolve structures well within this radius (e.g. an accretion disk, as proposed by Cannizzo et al. 1990). Nevertheless, there is some evidence of circularization of the gas flow in the velocity plots.

## 4 Discussion

We have performed a 1PN simulation of a solar mass star being disrupted by a super massive ($`10^6`$ M) black hole. The disruption process itself causes about 1/2 of the star’s matter to become unbound. We follow this matter up to and beyond the time when it returns to pericenter and starts accreting onto the BH. This in turn enables us to determine the amount of returning mass that is actually accreted. Contrary to previous, more heuristic estimates, we find that only about 25% of the returning mass actually gets accreted. The rest becomes unbound, following being heated by the strong compression accompanying the approach to pericenter.

The main consequence of this process is that the maximum accretion rate is a factor of 3 lower than expected from a simple examination of the rate of mass return to pericenter. The accretion rate into the volume resolved by our simulation is also quite constant, at about 1 $`\mathrm{M}_{}\mathrm{yr}^1`$, for the last 20 days of the simulation as opposed to a power law decay which is expected from the rate of return. When integrated over time, these results show that only about 10% of the original star’s mass actually gets accreted as opposed to the 65% expected if all the bound mass would be finally accreted. This means that the mass involved in the accretion flow, be it in the form of an accretion disk (e.g. Cannizzo et al., 1990), or spherical accretion (e.g. Loeb & Ulmer, 1997) is considerably smaller than previously estimated. Consequently, the duration of the expected ‘flare’ can be shorter by a factor 4–5 compared to these early estimates, which makes the detectability of these events much harder (the rather high temperature of some of the infalling debris also makes the bolometric correction high). Indeed, supernova searches in distant galaxies have failed so far to identify any such event unambiguously (Filippenko 2000, private communication).

The mass of the debris that becomes unbound because of the pericenter compression is comparable to the mass of the unbound debris resulting from the original disruption event. The main difference is that this new debris component has a different orbital distribution. Most notably, the velocity has a much larger component in the direction perpendicular to the orbital plane. This additional component of unbound high velocity debris could produce interesting consequences (e.g. components like Sgr A East in the Galactic center) when colliding with the interstellar medium surrounding the BH (e.g. Khokhlov & Melia, 1996). In particular, the fact that the mass ejection is more spherically symmetric, makes the event more similar to a normal supernova. Thus, tidal disruption events may produce a quite unique signature (in terms of their impact on surrounding gas), in which two supernova-type events are separated by a few weeks to a few months, with the first one being very anisotropic (mass being ejected within a solid angle $`\mathrm{\Omega }\stackrel{>}{}16(R_{}/R_\rho )^{1/2}(M_{}/M_{\mathrm{bh}})^{1/2}`$ rad<sup>2</sup>) and the second more spherically symmetric.

ML acknowledges support from NASA Grant NAG5-6857.


text:
# Modelling galaxy clustering at high redshift

## 1. Introduction

There has been significant recent progress in the study of galaxy formation within a cosmological context, mainly due to a phenomenological approach to this problem. The idea is to start with a structure formation model that describes where and when galactic dark haloes form. A simple description of gas dynamics and star formation provides a means to calculate the amount of stars forming in these haloes. Stellar population synthesis models then provide the spectral evolution, i.e. luminosities and colours, of these galaxies.

Many physical processes are modelled as simple functions of the circular velocity of the galaxy halo. Therefore, the Tully-Fisher relation is the most obvious observational relation to try and predict, as it relates the total luminosity of a galaxy to its halo circular velocity. However, most phenomenological galaxy formation models do not simultaneously fit the I-band Tully-Fisher relation and the B or K band luminosity function. When one sets the model parameters such that the Tully-Fisher relation has the right normalization, the luminosity functions generally overshoot (e.g. Kauffmann, White & Guiderdoni 1993; Kauffmann, Colberg, Diaferio & White 1999), certainly for the $`\mathrm{\Omega }=1`$, $`H_0=50`$ km s-1 Mpc-1 standard CDM cosmology (in the form given by Davis et al. 1985) that we consider in this paper. Alternatively, when making sure that the luminosity functions matches by changing some of the model parameters, the Tully-Fisher relation ends up significantly shifted with respect to the observed relation (e.g. Cole et al. 1994; Heyl et al. 1995).

In order to keep the modelling as analytical as possible, an extension to the Press & Schechter (1974) prescription for the evolution of galaxy haloes (e.g. Bond et al. 1991; Bower 1991; Lacey & Cole 1993; Kauffmann & White 1993) has been a popular ingredient for implementations of a phenomenological theory of galaxy formation. However, the EPS formalism is designed to identify collapsed systems, irrespective of whether these contain surviving subsystems. This ‘overmerging’ of subhaloes into larger embedding haloes is relevant to the problem of matching both the galaxy luminosity function and the Tully-Fisher relation, as the central galaxy in an overmerged halo is the focus of a much larger cooling gas reservoir than the reservoir that galaxy is to focus of in case its parent subhalo survives. Traditional N-body simulations suffer from a similar overmerging problem (e.g. White 1976), which is of a purely numerical nature, caused by two-body heating in dense environments when the mass resolution is too low (Carlberg 1994; van Kampen 1995).

In order to circumvent these problems, we use an N-body simulation technique that includes a built-in recipe for galaxy halo formation, designed to prevent overmerging (van Kampen 1995, 1997), to generate the halo population and its formation and merger history. This resolves most of the discrepancy sketched above, and allows us to make the modelling more realistic by adding chemical evolution and a merger-driven bursting mode of star formation to the modelling. Once stars are formed, we apply the stellar population synthesis models of Jimenez et al. (1998) to follow their evolution. We have enhanced these models with a model for the evolution of the average metallicity of the population, which depends on the starting metallicity. Feedback to the surrounding material means that cooling properties of that material will change with time, affecting the star formation rate, and thus various other properties of the parent galaxy.

## 2. Overview of the phenomenological model

The key ingredients of the model are described below. We refer to van Kampen et al. (1999) for a much more detailed description and discussion of the model, and a list of the choices for the various parameters involved.

### The merging history of dark-matter haloes.

This is often treated by Monte-Carlo realizations of the analytic ‘extended Press-Schechter’ formalism, which ignores substructure. We use a special N-body technique to prevent galaxy-scale haloes undergoing ‘overmerging’ owing to inadequate numerical resolution.

### The merging of galaxies within dark-matter haloes.

Each halo contains a single galaxy at formation. When haloes merge, a criterion based on dynamical friction is used to decide how many galaxies exist in the newly merged halo. The most massive of those galaxies becomes the single central galaxy to which gas can cool, while the others become its satellites.

### The history of gas within dark-matter haloes.

When a halo first forms, it is assumed to have an isothermal-sphere density profile. A fraction $`\mathrm{\Omega }_b/\mathrm{\Omega }`$ of this is in the form of gas at the virial temperature, which can cool to form stars within a single galaxy at the centre of the halo. Application of the standard radiative cooling curve shows the rate at which this hot gas cools below $`10^4`$ K, and is able to form stars. Energy output from supernovae reheats some of the cooled gas back to the hot phase. When haloes merge, all hot gas is stripped and ends up in the new halo.

### Quiescent star formation.

The star formation rate is equal the ratio of the amount of cold gas available and the star-formation timescale. The amount of cold gas available depends on the merger history of the halo, the star formation history, and the how much cold gas has been reheated by feedback processes.

### Starbursts.

We also model star bursts, i.e. the star-formation rate may suffer a sharp spike following a major merger event.

### Feedback from star-formation.

The energy released from young stars heats cold gas in proportion to the amount of star-formation, returning it to the reservoir of hot gas.

### Stellar evolution and populations.

Our work assumes the spectral models of Jimenez et al. (1998); for solar metallicity, the results are not greatly different from those of other workers. The IMF is generally taken to be Salpeter, but any choice is possible. Unlike other workers, we take it as established that the population of brown dwarfs makes a negligible contribution to the total stellar mass density, and we do not allow an adjustable $`M/L`$ ratio, $`\mathrm{{\rm Y}}`$, for the stellar population.

### Chemical evolution.

The evolution of the metals must be followed, for two reasons: (i) the cooling of the hot gas depends on metal content; (ii) for a given age, a stellar population of high metallicity will be much redder. The models of Jimenez et al. (1998) allow synthetic stellar populations of any metallicity to be constructed.

## 3. Low-redshift results

With the set-up described above we match both the B and K band luminosity function and the I-band Tully-Fisher relation, for an $`\mathrm{\Omega }=1`$ standard CDM structure formation scenario. Resolving the overmerging problem is the major contributor to this result, but the inclusion of chemical evolution and starbursts are also important ingredients.

The new ingredients we have added to the modelling of galaxy formation are needed in order to make the models more realistic, and are not introduced simply in order to give yet more free parameters. Nevertheless, our resolution to the Tully-Fisher / luminosity function discrepancy may well not be unique, and various other changes to the ingredients of the phenomenological galaxy formation recipe might produce similar results. For example, we have not studied the influence cosmological parameters have on the model galaxy populations, where $`\mathrm{\Omega }`$, $`\mathrm{\Lambda }`$, and $`\sigma _8`$ are likely to be the important parameters. Other types of ingredients are possible as well: Somerville & Primack (1998) resolve some of the discrepancy using a dust extinction model plus a halo-disk approach to feedback.

## 4. High-redshift clustering

One way of resolving the worries about degeneracies in the cosmological/physical parameter space will be to include data at intermediate and high redshifts, which is being gathered with increasing speed and ease, and at increasingly higher redshifts. In this contribution we show a preliminary comparison of the correlation properties of galaxies at redshift $`z=3`$. Recently, Giavalisco et al. (1998) gave an estimate for the galaxy-galaxy correlation function $`\xi (r)=(r_0/r)^\gamma `$ for a sample of Lyman-break galaxies at this redshift. They found $`r_0=2.1h^1`$Mpc and $`\gamma =2.0`$.

We selected our model galaxies in exactly the same way as Giavalisco et al. (1998) did, and compared two of the models produced by van Kampen et al. (1999), models $`n`$ and $`b`$, to the observational data. The first model ($`n`$), which is as close as possible to the model by Cole et al. (1994), but with the mass-to-light parameter $`\mathrm{{\rm Y}}=1`$, gives $`r_0=3.5h^1`$Mpc and $`\gamma =1.72`$. The second model ($`b`$), which includes starbursts and chemical evolution, gives $`r_0=4.4h^1`$Mpc and $`\gamma =2.1`$. Both models fit the correlation function at $`z=0`$ very well, and cannot be distinguished from each other.

As the observed correlation data are still relatively uncertain at this moment in time, it is premature to rule out models on the basis of this data. The two models discussed above have similar predictions for low redshifts, but predict different clustering properties at high-redshift. However, the differences are not large, so one needs either really good data, or a much larger variety of observational characteristics of the high-redshift galaxy population.

#### Acknowledgments.

Many thanks for the loud vocal support from outside the conference building during my presentation. I like to think that the people of Marseille just wanted to show how much they supported everything I said …

## References

Bond J.R., Cole S., Efstathiou G., Kaiser N., 1991, ApJ, 379, 440

Bower R.J., 1991, MNRAS, 248, 332

Carlberg R.G., 1994, ApJ, 433, 468

Cole S., Aragón-Salamanca A., Frenk C.S., Navarro J.F., Zepf S.E., 1994, MNRAS, 271, 781

Davis M., Efstathiou G., Frenk C.S., White S.D.M., 1985, ApJ, 292, 371

Heyl J.S., Cole S., Frenk C.S., Navarro J.F., 1995, MNRAS, 274, 755

Giavalisco, M., Steidel C.C., Adelberger, K.L., Dickinson, M.E., Pettini, M., Kellogg, M., 1998, ApJ, 503, 543

Jimenez R., Padoan P., Matteucci F., Heavens A., 1998, MNRAS, 299, 123

Kauffman G., White S.D.M., 1993, MNRAS, 261, 921

Kauffmann G., Colberg J.M., Diaferio A., White S.D.M., 1999, MNRAS, 303, 188

Kauffmann G., White S.D.M., Guiderdoni, 1993, MNRAS, 264, 201

Lacey C.G., Cole S., 1993, MNRAS, 262, 627

Press W.H., Schechter P., 1974, ApJ, 187, 425

Somerville R.S., Primack J.R., 1998, astro-ph/9802268

van Kampen E., 1995, MNRAS, 273, 295

van Kampen E., 1997, in Clarke D.A., West M.J., eds., Proc. 12th ‘Kingston meeting’ on Theoretical Astrophysics: Computational Astrophysics, ASP Conf. Ser. Vol. 123. Astron. Soc. Pac., San Francisco, p. 231, astro-ph/9904270

van Kampen E., Jimenez R., Peacock J.A, 1999, MNRAS, in press, astro-ph/9904274

White S.D.M., 1976, MNRAS, 177, 717


text:
# Formation and evolution of disk galaxies within cold dark matter halos

## 1. Introduction

The hierarchical cosmic structure formation picture based on the inflationary cold dark matter (CDM) provides a solid framework for models of galaxy formation and evolution. On the other hand, the unprecedented observations of galaxies at different redshifts make it possible to probe and constrain these models. Here we discuss some of the results obtained with a self-consistent scenario of disk galaxy formation and evolution within the context of the hierarchical picture (the extended collapse scenario).

## 2. Dark matter halos

Using the extended Press-Schechter approximation, we generate the mass aggregation histories (MAHs) of the dark matter (DM) halos. Collapse and virialization of these halos are then calculated assuming spherical symmetry and adiabatic invariance, using a method based on a generalization of the secondary infall model (Avila-Reese, Firmani, & Hernández 1998). These halos will mainly correspond to isolated systems. The diversity of MAHs results in diversity of density profiles which, in our model, mainly depend on the MAH. The density profile corresponding to the average MAH is well described by the Navarro et al. (1996) profile. In Avila-Reese et al. (1999) we have compared the outer density profiles, concentrations and structural relations of thousands of halos identified as isolated in a cosmological ($`\mathrm{\Lambda }`$CDM) N-body simulation with those obtained with our seminumerical method. We found a good agreement between the model and simulation results.

We have found that $`13\%`$ of the halos in the numerical simulation at $`z=0`$ are contained within larger halos and $`17\%`$ have significant companions within three virial radii. The remaining 70% of the halos are isolated objects. The slope $`\beta `$ of the outer density profile ($`\rho r^\beta `$) and the halo concentration defined as $`c_{1/5}=r_h/r(M_h/5)`$, where $`r_h`$ and $`M_h`$ are the virial radius and mass, depend on the halo environment. For a given $`M_h`$, halos in clusters have typically steeper outer profiles and are more concentrated than the isolated halos (for the latter $`\beta 2.9`$ in average and $`\beta `$ between 2.5 and 3.8 for 68% of the halos). Contrary to naive expectations, halos in galaxy and group systems as well as the halos with significant companions, systematically have flatter and less concentrated density profiles than isolated halos. A tight correlation between $`M_h`$ and the maximum circular velocity $`V_m`$ is observed: $`M_hV_m^n`$, $`n3.2`$. This is roughly the slope of the infrared Tully-Fisher relations (TFR). Thus, it seems that there is no room for the mass dependence of the infrared $`M_h/L`$ ratio.

## 3. Galaxy evolutionary models

We model the formation and evolution of baryon disks in centrifugal equilibrium within the growing CDM halos formed as described in §2. We assume that halos acquire angular momentum from large-scale torques with the spin parameter $`\lambda `$ distributed log-normally and constant in time. The disks are built inside-out with the gas infall rate (no mergers) proportional to the cosmological mass aggregation rate and assuming detailed angular momentum conservation. The gravitational drag of the disk on the DM halo is calculated. The local SF is assumed to be induced by disk instabilities and regulated by energy balance within the disk turbulent ISM (no SF feedback and self-regulation at the level of the interhalo medium is allowed). We also calculate the secular formation of a bulge. This way, at each epoch and at each radius, the growing disk is characterized by the infall rate of fresh gas, the gas and stellar surface density profiles, the total rotation curve (including the DM component), the local SF rate, and the size of the inner region transformed into bulge component.

## 4. Highlights of the model results

Results on the structure and dynamics of our model disk galaxies were discussed in Firmani & Avila-Reese (2000); the luminosity properties and topics related to the disk Hubble sequence were treated in Avila-Reese & Firmani (2000), while some evolutionary aspects of the galaxies were presented in Firmani & Avila-Reese (1999). In the following, we highlight some of the results.

Local properties. The (stellar) surface density and brightness profiles are exponential, the sequence of high to low surface brightness (SB) being mainly determined by $`\lambda `$. The gas profiles at $`z=0`$ are also exponential although much lower in density and with a scale radius $`24`$ times larger than the stellar profiles. There is a negative radial gradient of the color index: stars in the outer regions of the disk form later than stars in the inner regions. We find that the local SF rate per unit area correlates with the gas surface density as $`\mathrm{\Sigma }_{\mathrm{SFR}}(r)\mathrm{\Sigma }_g^n(r)`$ with $`n2`$ for most of the models and over a major portion of the disks. The shape of the rotation curves correlates with the SB ($`\lambda `$) and in most cases is approximately flat. The dark halo dominates in the rotation curve decomposition down to very central regions.

The infrared Tully-Fisher relations (TFR). The slope of the $`M_hV_m`$ relation of the CDM halos remains imprinted in the TFR and agrees with observations. This slope is almost independent of the assumed disk mass fraction $`f_d`$ when the disk component in the rotation curve decomposition is gravitationally important ($`f_d\genfrac{}{}{0pt}{}{_>}{^{}}0.03`$ for the $`\mathrm{\Lambda }`$CDM model used here). The zero point of the model TFR is only slightly larger than the observed zero point. The rms scatter in our TFR slightly decreases with mass; from $`V_m=70`$ to 300 km/s the scatter is between 0.38 and 0.31 mag. We have found that a major contribution to this scatter is from the scatter in the DM halo structures due to the dispersion of the MAHs; a minor contribution to the scatter is due to the dispersion of $`\lambda `$. The TFR for high and low SB models is approximately the same. The slope of the correlation among the residuals of the TF and luminosity-radius relations is small and non-monotonic, although the shape of the rotation curves of our models correlates with the SB. For a given total (star+gas) disk mass, the $`V_m`$ decreases with decreasing SB. However, owing to the dependence of the SF efficiency on the disk surface density, the stellar mass $`M_s`$ (luminosity) also decreases. This combined influence of the SB ($`\lambda `$) on $`V_m`$ and $`M_s`$ puts models of different SB on the same $`M_sV_m`$ relation. As a result, high and low SB models follow similar TFRs.

The Hubble sequence.  The main properties of the high and low SB disk galaxies and their correlations are determined by the combination of three fundamental physical factors and their dispersions: the halo virial mass, the MAH and the angular momentum given through $`\lambda `$. The MAH determines mainly the halo structure, the integral color index, and the gas fraction $`f_g`$ while $`\lambda `$ determines mainly the disk SB, the bulge-to-disk (b/d) ratio and the shape of the rotation curve. Our models show that the redder and more concentrated (higher SB) is the disk, the smaller is $`f_g`$ and the larger is the b/d ratio (disk Hubble sequence). The values of all these magnitudes are in good agreement with observations.

Evolutionary features. In the inside-out hierarchical disk formation scenario galaxies undergo not only luminosity but also structural (size, SB, b/d ratio) evolution. For an Einstein- de Sitter universe we find that the scale radius for normal disk galaxies decreases roughly as $`(1+z)^{0.5}`$ up to $`z1.5`$, while the central $`B`$band SB from $`z=0`$ to $`z=1`$ increases by $`1.2`$ mag.

The SF history in the models is driven both by the MAH and the disk gas surface density. For the average MAH and $`\lambda =0.05`$, the SF rate reaches a maximum at $`z1.52.5`$ which is a factor of 2.5-4.0 higher than the rate at $`z=0`$. In the same way, $`L_B`$ increases towards the past by factors slightly smaller than the SF rate. The less massive galaxies present a slightly more active luminosity evolution than the massive galaxies. The model galaxies are somewhat bluer in the past; from $`z=0`$ to $`z=1`$ the $`BV`$ decreases on average 0.25-0.30 magnitudes. The total mass-to-$`L_B`$ ratio also decreases towards higher redshifts: from $`z=0`$ to $`z=1`$ it decreases on average by a factor $`3.3`$, i.e. a galaxy at $`z=1`$ is more luminous in the $`B`$band and less massive than at $`z=0`$. Again, this is a result related to the hierarchical MAHs of the protogalaxies.

Owing to the mass (size) evolution, for a fixed $`V_m`$, the $`H`$band luminosity is a factor $``$2.2 less at $`z=1`$ than at $`z=0`$; however, owing to the luminosity evolution, $`L_B`$ is a factor $``$ 2.1 larger. Therefore, while the zero-point of the $`H`$band TFR increases towards the past, in the case of the $`B`$band TFR, compensation due to the $`L_B`$ evolution results in the zero-point remaining approximately constant with time. The slopes in both cases also remain constant.

## 5. Potential difficulties of the hierarchical scenario

Although several main properties, correlations, and evolutionary features of normal disk galaxies have been successfully predicted by our models, it is important to remark on their problems. We find the following potential conflicts with the observations: 1) the size and SB evolution of the disks is too pronounced, 2) the radial color index gradients are too steep and the $`f_g`$ is slightly over-abundant, 3) the DM component dominates in the rotation curve decompositions almost down to the center and the halos are too cuspy.

Regarding item 1), if selection effects in the deep field are not so significant as Simard et al. (1999) have claimed, then probably it is not so serious. In fact, some physical ingredients not considered in our models (e.g., merging, angular momentum transfer, and non-stationary SF) all work in the direction to improve models regarding problems 1) and 2). The problem 3) can probably be solved if the inner density profile of the CDM halos can be shallower than predicted (several solutions such as self-interacting CDM, warm DM, non-Gaussian fluctuations, have been proposed). Nevertheless, it is possible that all these problems together with those of the dearth of satellites and the high frequency of disk disruptive mergers, are in general pointing out to serious problems for the Gaussian CDM-based hierarchical picture of structure formation. More observational tests regarding the problems mentioned above and more theoretical effort in modeling galaxy formation and evolution are urgently required.

## References

Avila-Reese, V., Firmani, C., & Hernández X. 1998, ApJ, 505, 37

Avila-Reese, V, Firmani, C., Klypin A., & Kravtsov A. 1999, MNRAS, 309, 527

Avila-Reese, V, & Firmani, C. 2000, RevMexA&A, v. 36, in press

Firmani, C., & Avila-Reese, V. 1999, ASP Conf. Series 176, 406

Firmani, C., & Avila-Reese, V. 2000, MNRAS, in press

Navarro, J., Frenk, C.S. & White, S.D.M. 1997, ApJ, 462, 563

Simard, L. et al. 1999, ApJ, 519, 563