text:
# Dependence of the Inner DM Profile on the Halo Mass

## 1 Introduction

According to the currently favoured galaxy formation scenarios, galaxies are formed from the gravitational growth of tiny density perturbations imprinted on a uniform universe during inflation. Small mass perturbations grow faster and constitute the building blocks for the assembly of larger galaxies and clusters. In the nonlinear phase of the gravitational collapse the dark matter (DM) particles are shock heated to the virial temperature and settle into a dark halo with mean overdensity about $`\mathrm{\Delta }=178`$ (according to the simple spherical collapse model in a flat $`\mathrm{\Omega }_0=1`$ universe). Analysing N-body simulations of hierarchical structure formation Navarro, Frenk & White (1996, 1997) (hereafter, NFW) have proposed that DM halos develop a universal density profile, valid for virialized halos at any redshift and with any mass, from dwarf galaxies to clusters. Further works have shown that an approximatively universal profile develops regardless of the details of the adopted cosmology and initial power spectrum of density perturbations (Huss et al., 1999; Eke et al., 2001). The NFW density profile has a core cusp $`\rho r^1`$ and in the outer regions the density decreases as $`\rho r^3`$. The details of the gravitational collapse affect only the relative size of the core with respect to the virial radius (the halo concentration parameter, $`c_\mathrm{\Delta }`$) but not the slope of the profile. This result implies that during the virialization process the memory of the mass function and profile shapes of the building blocks that formed the halo is lost.

Recent observations of dwarf galaxies (spheroidal and irregular) and low surface brightness (LSB) galaxies seem to show that a flat core for the DM density profile is more compatible with the measurements of the rotation curves (van den Bosch et al., 2000; de Blok et al., 2001; de Blok & Bosma, 2002; Borriello & Salucci, 2001; Salucci & Burkert, 2000). These observations have renewed interest in the subject and a number of solutions have been proposed to solve this possible discrepancy. The proposed solutions range from feedback effects of the baryons and stars on the DM density profile to modifications of the physical properties (e.g., scattering cross section, temperature) of the DM particles. Some examples of alternatives to CDM are warm dark matter (Bode et al., 2001) and self-interacting DM (Spergel & Steinhardt, 2000). The potential discrepancy between the steep DM cores found in N-body simulations with the observed flat core of dwarf galaxies is perhaps the most serious problem faced in CDM cosmologies today.

High resolution simulations and careful convergence tests have been done or are underway to understand if the discrepancies found by some authors (e.g., Moore et al. (1998, 1999) find steeper core profiles) are due to numerical artifacts. In the core of DM halos the particle crossing time is much shorter than the Hubble time and the integration of trajectories is subject to subtle numerical effects that are difficult to keep under control. Power et al. (2002) estimated that about 1 million particles are needed in order to have convergent results in the inner 1% of the halo virial radius. To have such a high number of particles per halo the use of a special technique is required. This technique consists of simulating with high mass resolution only the particles that will end up in the halo of interest at $`z=0`$, while having less mass resolution for the particles that end up far away from the halo of interest. This has the drawback that each simulation can resolve only one halo at a time, therefore selection effect biases and a poor statistical sample could affect the reliability of the final result even if the simulation is very accurate and has high resolution. Note that most of the detailed and computationally expensive work on the subject has focused on simulating halos at $`z=0`$ with typical mass similar to the Milky Way.

The mass function of the building blocks of dwarf galaxies is different in comparison to the one for larger galaxies or clusters. On small scales (or masses) the mass fraction of virialized halos is about constant as a function of the logarithmic DM mass of the halos because the power spectrum has a slope $`n3`$. This means that small and large mass satellites contribute equally to the mass of the accreting halo. In massive halos as in the Milky Way or in clusters of galaxies, the contribution from large mass accreting satellites is instead dominant. For this reason it is more crucial to resolve very small mass satellites (i.e., to have high mass resolution) in simulating dwarf galaxies than in simulating massive galaxies or clusters.

In this paper I compare the density profiles of DM halos in four CDM N-body simulations that are identical apart from the box sizes that are 1 h-1 Mpc, 32 h-1 Mpc, 256 h-1 Mpc and 1024 h-1 Mpc. The simulation method is a two-level particle-particle particle-mesh (P3M) with $`256^3`$ particles and $`\mathrm{\Lambda }`$CDM cosmology. In dimensionless units the only difference between the simulations is the initial power spectrum of density perturbations. I compare the profiles when the most massive halos are composed of about $`10^5`$ DM particles and the clustering in the simulations is similar (see Fig. 2). This is done in order to minimise possible systematic errors in the simulations and to be more confident that any variation in the profile shapes is produced by the different initial power spectrum. The masses of the halos studied in the four simulations, from the small to the large box, are typical of dwarf galaxies, Milky Way mass galaxies, clusters and superclusters of galaxies, respectively.

If the cores of small mass halos are dense enough to resist disruption and survive “undigested” when incorporated into a larger object, they should determine the halo structure in the inner regions. Using scaling arguments Subramanian, Cen & Ostriker (2000) have shown that, in a flat universe with an initial power spectrum of density perturbation $`P(k)k^n`$, the core density profile is $`\rho _{dm}r^\alpha `$ with $`\alpha =(9+3n)/(5+n)`$. The numerical results that they report approximately confirm this expectation. In $`\mathrm{\Lambda }`$CDM cosmologies the logarithmic slope of the initial power spectrum ranges between $`3<n<1`$, depending on the mass scale (see Fig. 1). From dwarfs to superclusters the logarithmic slope of the core density profile should range between $`0\stackrel{<}{}\alpha \stackrel{<}{}2`$. As shown in Fig. 1, the power spectrum in the 1 h-1 Mpc box is a power law with $`n2.7`$. According to the aforementioned scaling relationship $`\alpha 0.4`$. This slope is in agreement with what I find in the 1 h-1 Mpc box. In the 32 h-1 Mpc box the power spectrum is not a single power law but the slope is close to $`n2`$ and therefore $`\alpha 1`$ is expected, consistent with the NFW universal profile and with the result found for this simulation. In cluster and supercluster size halos $`\alpha `$ should be larger than in galaxy size halos ($`\alpha 1.5`$). Indeed, in the 256 h-1 Mpc simulation I find slopes of the inner profile $`\alpha 1.4`$. But this is not a strong result because the number of particles in each halo is barely sufficient to start seeing a deviation with respect to the NFW profile. Moreover, the introduction of the parameter $`\alpha `$ for the slope of the inner part of the halo profile (which can be easily calculated knowing the mass of the galaxy halo) does not increase the number of adjustable parameters in the profile fitting formula. This happens because the concentration parameter, $`c_\mathrm{\Delta }^{NFW}(M_{dm})`$, needed to fit the profiles with the NFW formula, is in this case a constant. The changing slope of the inner profile mimics the variation of the concentration parameter in the NFW profile, shifting the value of the radius where the circular velocity reaches its maximum value (see appendix A).

This paper is organised as follows: in § 2 I explain the methods for the simulations and in § 3 I show the results for dwarf galaxies, normal galaxies and clusters of galaxies. In § 4 I demonstrate how a density profile with changing inner slope $`\alpha `$ can fit all the profiles from dwarfs to clusters with a constant concentration parameter. The discussion and conclusions are presented in § 6. The equations for the circular velocity and integrated mass for the NFW profile and for a halo with arbitrary inner and outer slopes of the density profile are contained in appendix A.

## 2 The Method

In the literature most efforts have focused on studying galaxies at $`z=0`$ with masses typical of the Milky Way. Ricotti et al. (2002a, b) have performed high-resolution simulations of the formation of the first galaxies using a cosmological code that solves the equations for the DM particles, baryons, stars and radiative transfer. These simulations achieve a mass resolution of $`M_p=4.9\times 10^3`$ M for the DM using $`256^3`$ particles and 1 h-1 Mpc box size. In these simulations at $`z10`$ or, expressing the redshift in terms of the scale factor, at $`a=(1+z)^10.09`$, the profiles of the most massive halos ($`M_{dm}2\times 10^8`$ M) appear flatter than the NFW profile. I therefore perform the same simulation for only the DM particles, finding again systematically flatter halo cores than found by NFW. Since this result could be produced by subtle numerical problems of the N-body integration or the small number of particles in the halo (typically between $`10^4`$ and $`10^5`$), I repeat the simulation for a 32 h-1 Mpc box down to $`a=0.3`$ when the most massive halos have the same number of DM particles as in the 1 h-1 Mpc box at $`a=0.09`$. In this case the profiles are well fitted by the NFW universal profile. The differences of the profiles in this two simulations are likely to be caused by the different initial power spectrum of density perturbations, since the they have identical smoothing length in dimensionless units, identical Fourier modes of the initial density perturbations, similar number of integration time steps and comparable number of particles per halo. According to this paradigm and the results of Subramanian et al. (2000), the slope of the DM profiles in simulations with box sizes larger than 32 h-1 Mpc should be steeper than the NFW profiles. I perform two additional simulations with 256 h-1 Mpc and 1024 h-1 Mpc box sizes to test this hypothesis. The large box simulations present the problem that in order to achieve the same degree of clustering as in the smaller boxes it is necessary to evolve the simulation in the future when $`a\stackrel{>}{}1`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }>0.7`$. For the 256 h-1 Mpc simulation the required level of clustering is reached when $`0.9<a<1.5`$ and the number of particles per halo is barely sufficient to start noticing a steepening of the inner profile with respect to the NFW profile. In the case of the 1024 h-1 Mpc box the level of clustering reached in the smaller boxes can never be achieved because the universe begins to inflate (the scale factor becomes $`a>1000`$ in a few time steps) and the event horizon starts decreasing, preventing any further clustering. I present the results for this simulation because they are interesting for understanding what determines the outer slope of the density profile. The number of particles in these halos is too small for measuring the slope of the inner profile but it is evident that, in an accelerating universe, the outer slope of the density profile is steeper. Analytically, Subramanian et al. (2000) show that the slope of the outer profile is $`\gamma 4`$ in a low density universe and in an accelerating $`\mathrm{\Lambda }`$CDM universe. Table 1 lists the scale factors when the simulations with different box sizes, $`L_{box}`$, have run for a comparable number of time steps and the most massive halos are composed of the same number of DM particles.

In the following paragraphs I explain in greater detail the numerical techniques and the sanity checks adopted to derive the density profiles and the rotation curves.

### 2.1 N-body Simulations

I use an N-body code based on the P3M method with two levels of mesh refinement (Bertschinger & Gelb, 1991; Gnedin & Bertschinger, 1996). In all the simulations the number of particles is $`256^3`$ and the Plummer softening parameter is $`ϵ=0.02`$ in cell units $`\mathrm{\Delta }L`$, where $`\mathrm{\Delta }L=L_{box}/256`$ and $`L_{box}`$ is the box size. Unless stated otherwise, I use scale-free (dimensionless) units in this paper. The cosmology adopted is a flat $`\mathrm{\Lambda }`$CDM model with $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_m=0.3`$ and $`h=0.7`$. The initial conditions for the density and velocity fields are calculated using the COSMIC package (Bertschinger, 1995) assuming a scale-invariant (tilt $`n=1`$) initial spectrum of density perturbations and $`\mathrm{\Lambda }`$CDM transfer function. The spectrum, shown in Fig. 1, is normalised imposing a variance $`\sigma _8=0.91`$ in spheres with radius of 8 h-1 Mpc at $`z=0`$. The initial conditions are calculated using the Zel’dovich approximation until $`z=166,91,48`$ and $`26`$ for the $`L_{box}=1,32,256`$ and $`1024`$ h-1 Mpc boxes, respectively. The N-body simulations start at those redshifts. I use the same random realization for the wave phases and directions in all simulations.

The halos are identified using DENMAX (Bertschinger & Gelb, 1991) with smoothing parameter $`G=1000`$. The DENMAX algorithm identifies halos as maxima of the smoothed density field with smoothing length $`L_G=L_{box}/G`$. The algorithm assigns each particle in the simulation to a group (halo) by moving each particle along the gradient of the density field until it reaches a local maximum. Unbound particles are then removed from the group. The results of DENMAX depend on the degree of smoothing used to define the density field. A finer resolution in the density field will split large groups into smaller subunits and vice versa. This arbitrariness in results is a common problem of any group-finding algorithm, and it is not only a numerical problem but often a real physical ambiguity. Especially at high redshift, since the merger rate is high, it is difficult to identify or define a single galaxy halo. In appendix B in Ricotti et al. (2002a) the mass functions obtained using DENMAX with different smoothing lengths are compared with the analytical expectation using the Press-Schechter formalism. It was found that the smoothing parameter $`G=1000`$ gives the best results for boxes with $`256^3`$ particles.

In each simulation with a different box size the profiles of the halos are compared when the level of clustering is similar and the positions of the most massive halos coincide. This happens after a comparable number of time steps, but the number of steps required becomes larger as the box size increases since the halos are more concentrated. The most massive halos in the simulations are composed of about $`N_p=10^410^5`$ DM particles and the typical number of time steps to reach this level of clustering is a few thousands (see Table 1). Fig. 2 shows the sizes and positions of bound halos in the 1 h-1 Mpc box at $`a=0.085`$ and in the 32 h-1 Mpc box at $`a=0.230`$, when the clustering in the two simulations is similar.

### 2.2 Density Profiles

The density profiles are calculated in three different ways:
(i) Using the list of particles identified as halo members by DENMAX. This means all the bound particles.
(ii) Using, additionally, the particles identified by DENMAX as unbound particles.
(iii) Using all the particles. This includes all the halo satellites and neighbouring halos.
Usually profiles using (i) and (ii) are almost indistinguishable. But the profiles using (iii) look different in the outer parts if the halo is disturbed by accreting satellites. The profiles obtained with method (iii) are used for the analysis as in previous studies. The other profiles are useful to discriminate between isolated halos and halos that are undergoing accretion of massive satellites. One can choose whether to exclude or include in the analysis halos that are undergoing major mergers; either choice does not change the result significantly.

The profiles are shown in dimensionless units: the radii are expressed in number of cells $`x`$ (the cell size is $`\mathrm{\Delta }L=L_{box}/256`$ in comoving coordinates) and the masses in number of DM particles $`N_p`$. The mass and density of each halo are sampled as a function of the radius, $`x`$, in shells with constant logarithmic spacing $`dx/x=0.01`$. The dimensionless circular velocity is more appropriate to compare the profiles since the virial radii are difficult to determine when the halos are not isolated. Moreover, the relative error on the density profile is smaller at radii corresponding to the maximum of the circular velocity since the number of particles per shell at intermediate radii is maximum. Therefore, if the halo profile is not disturbed by accreting satellites, the location and the value of the maximum circular velocity provides the best method to compare halos of different masses. The dimensionless circular velocity is defined as,

$$v_c=\frac{V_c(x)}{V_c^0}=\sqrt{\frac{N_{dm}(x)}{x}}$$

where $`V_c^0=(GM_p/\mathrm{\Delta }L)^{1/2}=7.365\times 10^2(L_{box}/1h^1\mathrm{Mpc})`$ km s-1. If the slope of the inner density profile is $`\rho (r)r^\varphi `$, the circular velocity has a slope $`V_cr^{(1\varphi /2)}`$. The Poisson error for the density is $`\delta \rho /\rho =dN_p(x)^{1/2}`$, where $`dN_p(x)`$ is the number of particles in each shell of radius $`x`$. The error on the mass as a function of the radius is $`\delta M_{dm}/M_{dm}(x)=N_p(x)^{1/2}`$, where $`N_p(x)`$ is the number of particles in each sphere of radius $`x`$. The error on the circular velocity is $`\delta v_c/v_c=0.5\delta M_{dm}/M_{dm}`$. Table 2 lists the values of useful constants for converting dimensionless units into physical units.

## 3 Results

In Figs. 3-5 (top) the dimensionless circular velocities are shown for the nine most massive halos in the simulations with box sizes $`L_{box}=1,32`$, and 256 h-1 Mpc, respectively. The statistical errors, $`\delta v_c`$, on each point are about the sizes of the symbols and therefore are not shown to avoid excessive crowding. I show the profiles obtained using methods (ii) (small dots) and (iii) (large dots). I do not show those obtained using method (i) since they are almost indistinguishable from (ii). I show the fit with a NFW profile obtained by matching the position and value of the maximum of $`v_c`$ and the virial radius (arrow) for the halos that are not severely perturbed by mergers. The points are connected in the range of radii that are reliable: $`x_{mi}<x<x_\mathrm{\Delta }`$. Here $`x_\mathrm{\Delta }=2.37(\mathrm{\Delta }/178)^{1/3}(N_p/10^4)^{1/3}`$ is the dimensionless virial radius and $`x_{mi}`$ is the radius that contains a number of particles $`N_p(x_{mi})300(\rho (x_{mi})/10^4)^{1/2}`$. This condition ensures that two-body relaxation is not important in flattening the inner profile (see Power et al. (2002) for details on convergence studies). Note that $`x_{mi}`$ is always larger than twice the Plummer softening length $`2ϵ=0.04`$, therefore an incorrect force calculation does not affect the results at those radii. In summary, the sanity checks adopted are as follows.
1) Only radii well above twice the softening length are considered.
2) Only radii that enclose $`N_p300(\rho /10^4)^{1/2}`$ particles are considered.
3) No merging halos: profiles (ii) and (iii) give similar results for the value and radius of the maximum circular velocity.
4) The same criteria are adopted for all the simulations with different box sizes, $`L_{box}`$, and scale factor $`a`$.

Figs. 3-5 (bottom) show the mean profiles obtained by scaling (shifting in logarithmic scale) the radii to the same $`x_{max}=1`$ and normalising the maximum of $`v_c`$ to $`v_c^{max}=1`$. Each radial bin, logarithmically spaced, is weighted equally so long as $`N_p300(\rho /10^4)^{1/2}`$. Only the halos that are not undergoing a major merger are considered. If all the ten profiles or a different subset of profiles are considered, I find that the resulting mean profile does not change significantly. Also, using the arithmetic or geometric mean does not change the mean profiles. A possible criticism to this method is that the halos analysed are not isolated as in works where a single galaxy is re-simulated with higher resolution. Time-dependent effects could thus produce flatter or steeper cores. For this reason, I present the mean profile at different redshifts (listed in the figure captions). If the halos are not relaxed to their final configuration the mean profiles should differ at different redshifts. Moreover, the increased statistical sample of halos makes the result more robust.

In the 1 h-1 Mpc simulation (see Fig. 3) the slope of the inner profiles is, on average, flatter (i.e., the circular velocity rises more steeply) than the NFW profile. The typical masses of the halos in this simulation are $`M_{dm}10^8`$ M. Even at radii close to $`v_{max}`$ the circular velocity starts to deviate from the NFW fit. The slope of the outer profile is instead consistent with $`\gamma =3`$ as in the NFW formula. Note that method (ii) does not attribute several DM particles to the halo, producing an outer profile that decreases more steeply. It is not clear if DENMAX fails not to include those particles or if they effectively do not belong to the halo. The definition of halo boundary (or virial radius) is somewhat arbitrary and time dependent since halos are constantly growing. Especially for small mass halos the flat mass function of the accreting satellites and the short accretion timescale makes it more difficult to define the outer edge of the galaxy. This is also true for dwarf galaxies at $`z=0`$, even if the accreting satellites are mostly invisible because they do not contain stars. For this reason an “isolated” dwarf galaxy cannot be found or defined in the same way as for massive galaxies for which the mass function of accreting satellites is dominated by halos with masses similar to the host halo. There are halos for which the use of either method (iii) or (ii) gives a different radius where the circular velocity is maximum. Those halos are excluded from the analysis since they are severely disturbed by accreting satellites (but if they were to be included in the analysis, the mean slope would remain flatter than $`\alpha =1`$). The mean density profiles are best fitted by inner slopes in the range $`\alpha 0.40.5`$, equivalent to mean slopes of the inner circular velocity $`\beta 0.750.8`$.

In the 32 h-1 Mpc simulation (see Fig. 4) the NFW profile provides a good fit to all the halos that are not severely disturbed by accreting satellites. The typical masses of the halos in this simulation are $`M_{dm}5\times 10^{12}`$ M. The mean density profiles are best fitted assuming inner slopes in the range $`\alpha 0.91`$, equivalent to mean slopes of the inner circular velocity $`\beta 0.50.55`$.

In the 256 h-1 Mpc simulation (see Fig. 5) the inner density profile slopes appear slightly steeper than predicted by NFW (i.e., flatter circular velocities). Unfortunately, the number of particles in the inner region is barely sufficient to start seeing a discrepancy. A larger number of particles per halo is needed to construct reliable rotation curves at small radii where the discrepancy, expected by extrapolating the reliable part of the rotation curve, should become evident. The typical masses of the halos in this simulation are $`M_{dm}2\times 10^{15}`$ M. The mean density profiles are best fitted by inner slopes in the range $`\alpha 1.31.4`$, equivalent to mean slopes of the inner circular velocity $`\beta 0.30.35`$. In Fig. 5 (bottom), it appears that the slope of the outer profile becomes steeper when $`a>1`$. For this reason, included in the figure are the best fits for two cases: (i) when the outer slope of the profile is $`\gamma 3`$ at $`a=1.0,1.2`$, and (ii) when the slope is $`\gamma 4.0`$ at $`a=1.5,1.7`$. This dependence of the outer profile slope on the scale factor becomes more evident in the 1024 h-1 Mpc simulation and is discussed in the following paragraph.

In Fig. 6 (top), the dimensionless circular velocities are shown for the nine most massive halos in the simulation with box size $`L_{box}=1024`$ h-1 Mpc at $`a=10`$. The number of particles in these halos is too small for measuring the slope of the inner profile but it is evident that the outer slope is steeper than the NFW profile. Fig. 6 (bottom) shows the density profiles (solid lines) and slope of the density profiles (dashed lines) for the nine more massive halos in the same simulation. The outer density profile has a slope $`\gamma 4`$. At $`a=2.3`$ the slope is $`\gamma 3.5`$ and increases to $`\gamma 4`$ when the scale factor is $`a=4,10,100`$ and 1000. The physical explanation of this result is likely to be related to the expansion rate of the universe. When $`\mathrm{\Omega }_\mathrm{\Lambda }(z)0`$ at $`a\stackrel{<}{}0.5`$ the universe is decelerating and $`\gamma 3`$. As $`a`$ increases to $`a\stackrel{>}{}1`$ the universe starts inflating (at $`a=2`$, $`\mathrm{\Omega }_\mathrm{\Lambda }0.95`$ in a flat universe) and the slope of the outer profiles becomes $`\gamma 4`$. It is possible that even at $`z=0`$ the outer profile of recently virialized clusters is slightly steeper than $`\gamma =3`$. According to the aforementioned interpretation, the observation of a steep outer profile in nearby clusters would indicate that the universe is accelerating (see also the results of Subramanian et al. (2000) for low density universes). In § 4 I show that for the halo profiles in this simulation the radii, $`r_{max}`$, where the circular velocity is maximum, decrease with increasing scale factor, consistent with the predictions assuming an inner slope of the density profile $`\alpha 1.5`$.

### 3.1 Radiative Feedback Effects on the DM profile

Fig. 7 shows the mean $`v_c`$ for the nine most massive halos in the simulation 256L1p3 presented by Ricotti et al. (2002b) at $`a=0.075`$ and $`a=0.090`$. This simulation attempts to simulate realistically the formation of the first galaxies in the universe modeling radiative feedback effects in the early universe. The parameters of the simulation are the same as in the 1 h-1 Mpc simulation presented in this work. But in addition to DM particles it includes gas dynamics, star formation using a Schmidt-Law, metal enrichment from star formation and radiative transfer. Radiative feedback effects produced by UV radiation emitted by the first stars are calculated following the molecular processes involving H2 formation, destruction and cooling. From the results presented in this work it appears that radiative feedback processes are not responsible for the flattening of the DM profile in small-halos. The slopes of the inner profile of the halos in this simulation are similar to the corresponding N-body simulation, or perhaps slightly steeper.

## 4 Constant concentration DM profiles

In comoving coordinates, $`r`$, all of the simulated halos in this work can be fitted by a profile with the form (cf., Zhao, 1996)

$$\rho _{dm}(X)\frac{1}{X^\alpha (1+X)^{\gamma \alpha }},$$
(1)
where $`X=c_\mathrm{\Delta }r/r_\mathrm{\Delta }`$, $`r_\mathrm{\Delta }`$ is the comoving virial radius, $`\alpha =\alpha (M_{dm})`$ is a function of the halo mass, $`M_{dm}`$, and $`c_\mathrm{\Delta }7`$ is the concentration parameter that, as I will show, is a constant. In physical coordinates, $`R`$, it is simply $`X=c_\mathrm{\Delta }R/R_\mathrm{\Delta }`$, where $`R_\mathrm{\Delta }=(1+z)^1r_\mathrm{\Delta }`$ is the virial radius. Given the initial power spectrum of density perturbations, the parameter $`\alpha `$ is a function of the halo mass only. The values found in this work agree with the predictions of Subramanian et al. (2000) (see Fig. 8). An acceptable fit to the theoretical value of $`\alpha `$ for the case of the $`\mathrm{\Lambda }`$CDM cosmology, adopted in this work, has the form $`\alpha 1.3+(M_{dm,14}^{1/6}1)/(M_{dm,14}^{1/6}+1)`$, where $`M_{dm,14}=M_{dm}/(3\times 10^{14}`$ M). This fit is five percent accurate for masses $`10^6\stackrel{<}{}M_{dm}/M_{}\stackrel{<}{}10^{18}`$.

The logarithmic slope, $`\gamma `$, of the outer parts of the halo appears to depend on the acceleration of the universe: when the scale parameter is $`a\stackrel{<}{}1`$, the slope is $`\gamma 3`$ as in the NFW profile, but is $`\gamma 4`$ at $`a>1`$ when $`\mathrm{\Omega }_\mathrm{\Lambda }1`$ and the universe is inflating. Equation (1) is normalised imposing the condition that the integrated mass inside the virial radius is $`m_{dm}(X=c_\mathrm{\Delta })=M_{dm}`$ (see appendix A). The comoving virial radius is defined as

$$r_\mathrm{\Delta }(32\mathrm{kpc})\left(\frac{\mathrm{\Delta }}{178}\right)\left(\frac{M_{dm}}{10^6M_{}}\right)^{1/3},$$

to satisfy the relation $`M_{dm}=(4\pi /3)\rho _0\mathrm{\Delta }r_\mathrm{\Delta }`$, where $`\mathrm{\Delta }`$ is the mean overdensity of the halo and $`\rho _0`$ is the mean DM density at $`z=0`$. $`\mathrm{\Delta }`$ depends on the cosmology since it is related to the mean overdensity and the time it takes for a density perturbation to stop expanding with the universe and “turn around,” starting the collapse. According to the simple top-hat spherical collapse approximation, $`\mathrm{\Delta }178\mathrm{\Omega }_0(z)^{0.45}`$ for a flat universe with a cosmological constant (Eke et al., 1998).

I now explain why this profile has the property of a constant concentration parameter. The changing slopes of the inner and outer parts of the profile mimics the variation of the concentration parameter in the NFW profile, shifting the value of the radius, $`r_{max}`$, where the circular velocity reaches its maximum value. For the NFW profile (i.e., for $`\alpha =1`$) $`r_\mathrm{\Delta }/r_{max}=c_\mathrm{\Delta }^{NFW}/2.16`$, therefore the location of $`r_{max}`$ with respect to the virial radius is inversely proportional to the concentration parameter. In the general case

$$\frac{r_\mathrm{\Delta }}{r_{max}}=\frac{c_\mathrm{\Delta }}{f(\alpha ,\gamma )}.$$
(2)
The location of the maximum circular velocity depends on the inner and outer profile slopes and on the concentration parameter. In appendix A I show that a good approximation for the function $`f(\alpha ,\gamma )`$ is $`f(4/\gamma )^{2.68}(2\alpha )`$. In Fig. 8 it can be seen that the values of $`r_\mathrm{\Delta }/r_{max}`$ (circles), for all the simulations in this work (listed in Table 1), agree with equation (2), assuming constant concentration parameter $`c_\mathrm{\Delta }=7`$. Moreover the values of $`\alpha `$ (triangles) are in good agreement with the theoretical relationship proposed by Subramanian et al. (2000), $`\alpha =(9+3n)/(5+n)`$ (dashed line), where $`n`$ is the effective slope of the power spectrum. The solid and open circles show $`r_\mathrm{\Delta }/r_{max}`$ measured in profiles with $`\gamma =3`$ (i.e., at $`a<1.2`$) and with $`\gamma >3`$, respectively. The solid and dash-dotted lines show $`r_\mathrm{\Delta }/r_{max}`$ given by equation (2) using the theoretical value of $`\alpha `$ (shown by the dashed line), $`\gamma =3`$ and $`\gamma =3.54`$, respectively. In the 1024 h-1 Mpc box simulation I find that $`\gamma 4`$. The inner slope of the profile cannot be measured because of the small number of particles in the inner regions. But the measured values of $`r_\mathrm{\Delta }/r_{max}`$ are consistent with values of the inner slope $`\alpha 1.5`$, predicted by the theory.

## 5 Discussion

The results presented in this work and their interpretation seem to disagree with most previously published work, although the simulations presented here are not of substantially higher resolution than those in the literature. In this section I analyse the discrepancies with previous works and point out that the results do not disagree strongly with most previous works, save for the 1 h-1 Mpc simulation. The result of the 1 h-1 Mpc simulation cannot be compared with previous works since a similar simulation has not been performed before. The joint analysis of the three simulations with $`L_{box}=1,32`$ and 256 h-1 Mpc suggests that the slope of the inner part of the DM profile depends on the power spectrum of linear perturbation. This interpretation of the results disagrees with the common view of universality of DM profiles, although several groups have found results that contrast with this idea, both in the context of galactic dynamics (Nipoti et al., 2003, e.g.,) and cosmology (Syer & White, 1998; Kravtsov et al., 1998; Jing & Suto, 2000; Jing, 2000). Here I identify and summarise three possible reasons why the present results differ from a large number of published work on N-body simulations.

1. The profiles of halos with typical masses resolved by the 1 h-1 Mpc simulation have not been analysed before. This simulation is equivalent to a scale free simulation with linear power spectrum $`P(k)k^n`$ with $`n=2.7`$ and Einstein-de Sitter cosmology, since at $`z>2`$ the cosmological constant contribution is negligible. Previous works have studied scale free simulations with power spectrum index $`n2`$. For instance Cole & Lacey (1996) studied the $`n=0,1,2`$ cases and Navarro et al. (1997) the $`n=0,0.5,1,1.5`$ cases. The simulations in those works have the same resolution as in the present study. Their results do not contrast with the findings of this work because when $`n\stackrel{>}{}2`$ the inner slope of DM profiles is also $`\alpha 11.5`$. But, as mentioned, the resolution of the simulations (about $`10^410^5`$ particles per halo) is marginally sufficient to discern between $`\alpha 1`$ and $`\alpha 1.5`$. The $`n=2.7`$ case is markedly different from the others since it corresponds to a scenario in which all mass perturbations become non-linear at about the same time (an index $`n<3`$ would correspond to the top-down formation scenario, in which small mass halos form after the large ones). For this reason, I believe that the result found analysing the 1 h-1 Mpc simulation is new and complements previous studies.
2. Interpretation of the results. Theoretically it is unclear what determines the slope of the inner DM profiles. Some works argue that the slope depends on the power spectrum of initial perturbations (Syer & White, 1998; Subramanian et al., 2000), others argue that the slope tends asymptotically to a steep slope $`\alpha 12`$ (Alvarez et al., 2003; Dekel et al., 2003). Recent still unpublished high-resolution N-body results from Navarro and collaborators show that the inner slope does not converge to $`\alpha =1`$ but the slope slowly flattens to $`\alpha \stackrel{<}{}1`$. This result contrasts with most aforementioned semi-analytical results indicating that physical processes that determine the density profile of DM halos and their time evolution are not well understood.

The halos analysed in the present work are the most massive halos at a given redshift. For this reason they are recently virialized halos with different masses, compared after about the same number dynamical times from their formation. The halos analysed in high-resolution N-body studies are usually chosen to be “isolated”: halos that did not experience “recent” major mergers. Therefore, depending on their mass, their age from the virialization differs greatly at $`z=0`$. Semi-analytical works study the equilibrium profile of halos given an accretion history. The different evolutionary states of the halos studied in literature could explain the discrepancies in the results if the time scale to reach the equilibrium configuration is long compared to the accretion time. I find that the halo profile depends on the accretion history. But it is possible that this is true only during a transient period of time, before the halo relaxes to a different equilibrium configuration. This concern is especially valid for the 1 h-1 Mpc simulation since in a scale free simulation with $`n2.7`$ the halos continuously accrete satellites. The result is therefore statistically significant at redshift about 10 but the shallow slope of the profiles might not survive to redshift $`z=0`$. Even if the shallow profiles are transient, depending on the time scale to reach the equilibrium configuration, the result at $`z=10`$ might suggest a reason for why not all dwarf galaxy halos at $`z=0`$ present shallow cores. According to the work of Dekel et al. (2003), substructure is not tidally striped but instead is compressed when the halo profile is shallow. Therefore substructure could survive undigested for quite a long time. Moreover even at $`z=0`$, small mass halos, due to their larger number density, are expected to capture dark satellites of non-negligible mass with higher frequency than larger mass halos. Further work is needed, and is currently in progress (Ricotti & Wilkinson, 2003), to extend the result of the present work from redshift $`z=10`$ to redshift $`z=0`$. The aim of this study is to understand if the flat profiles at $`z=10`$ are in a stable configuration that would survive for a Hubble time and if the profiles are consistent with the observed flat cores of dSph galaxies in the local group. Preliminary results show that the dynamical state of the halos is stable and the inner profiles remain flat for $`20`$ Grys if the halo is evolved in isolation.

To further support the case for a dependence of the inner profile on the mass of the halos I show that flat inner profiles are associated with an almost isotropic velocity dispersion of DM particles while steeper profiles have a velocity dispersion slightly biased in the radial direction. This result is shown in Fig. 9. The lines show the parameters $`\beta =1\sigma _\theta ^2/\sigma _r^2`$ (dashed lines) and $`\beta ^{}=1v_\theta ^2/v_r^2`$ (solid lines) as a function of the scaled radius $`x=r/r_{max}`$, where $`r_{max}`$ is the radius where the circular velocity is maximum. Here, $`\sigma _\theta ^2=v_\theta ^2v_\theta ^2`$ and $`\sigma _r^2=v_r^2v_r^2`$, where $`v_r`$ is the radial velocity and $`v_\theta `$ is the tangential velocity. The $`\beta `$ parameter measures the degree of anisotropy of the velocity dispersion: $`\beta =0`$ corresponds to an isotropic velocity dispersion, $`\beta <0`$ to a tangentially biased velocity dispersion and $`0<\beta <1`$ to a radially biased velocity dispersion. The lines from top to bottom show the mean value of $`\beta `$ and $`\beta ^{}`$ for the 9 most massive halos in the 256 h-1 Mpc, 32 h-1 Mpc and 1 h-1 Mpc simulations. The shallow inner DM profiles in the 1 h-1 Mpc simulation have almost isotropic velocity dispersion. Steeper DM profiles in the 32 h-1 Mpc and 256 h-1 Mpc simulations have a velocity dispersion slightly biased in the radial direction. This result is consistent with the analytical expectations based on the solution of the Jeans equations in spherical systems.
3. Method and analysis of the results In this paper I present statistical results without preselecting the halos to re-simulate with higher resolution as has been done in previous studies. In order to compare simulations with observations a statistically significant sample should be analysed. Halos that did not reach an asymptotic equilibrium configuration (if they exist) or halos presenting substructures produced by recent mergers should be included in the analysis if they appear to be statistically significant. Especially if there is not an easy way to determine the dynamical stage of a halo because the accreting substructure is invisible (small halos at $`z=0`$ should accrete mostly dark halos), understanding what fraction of dark halos of a given mass is in a transient configuration is relevant to solve the problem of the DM cores observed in a fraction of galaxies at low redshift.

As expected the halos studied in this work show substantial substructure. 3D images of the halos show that numerous satellites are present, independently of the halo mass. Perhaps in small mass halos substructure can survive “undigested” closer to the halo cores. Nevertheless, in most cases, this does not affect the approximately spherical symmetry of the halos and the determination of the halo centres. The halo centres are calculated, using DENMAX, as the minimum of the gravitational potential. I also calculate the centres as the maximum of halo density and as the centre of mass of bound particles. The halo centres calculated using those alternative definitions are practically coincident and choosing either of them does not change the profiles. Exceptions are profile 8 in the 1 h-1 Mpc simulation and profile 4 in the 256 h-1 Mpc simulation that have maximum density centres offset by 2$`ϵ`$, where $`ϵ`$ is the Plummer softening parameter.

In this work the inner slope of the profiles and the core radii are determined fitting the circular velocities instead of the density profiles. This procedure partially remove the degeneracy between the slope $`\alpha `$ of the inner profile and the core radius $`r_s`$. This is illustrated in Fig. 10, which shows the maximum likelihood estimate of the parameters $`\alpha `$ and $`r_s`$ fitting the density profile (left) or the circular velocity (right). The $`1\sigma ,2\sigma `$ and $`3\sigma `$ confidence contours are shown assuming constant relative errors on the density and circular velocity in each logarithmically spaced radial bin.

## 6 Conclusions

Two main problems faced by CDM cosmologies are related to the properties of small mass galaxies. Namely, (i) the number of visible galactic satellites (dwarf galaxies) in the Local Group is smaller than predicted by N-body simulations (Moore et al., 1999; Klypin et al., 1999) and (ii) the flatness of the DM cores in dwarf and LSB galaxies is not reproduced by N-body simulations. The first problem can be solved by including radiative feedback effects (Chiu et al., 2001; Ricotti et al., 2002b) in cosmological simulations. Galaxies with masses $`M_{dm}\stackrel{<}{}10^810^9`$ M are too small to retain the gas that is photoevaporated by stars or by the ionising background after reionization. Consequently, the luminosity of most dwarf galaxies is predicted to be too low to be detected (i.e., the mass-to-light ratio should increase for halos with smaller masses) or zero for very small mass halos. Here I have shown results of N-body simulations suggesting that the second problem might not contradict the predictions of CDM cosmologies. But further investigation is needed to understand if the dwarf galaxies showing shallow DM profiles at $`z10`$ (e.g., just after the time when most of the dwarf galaxies formed) are in a configuration that will survive unchanged until redshift $`z=0`$. In summary the main findings of this work are:

* The slope, $`\alpha `$ of the inner profile of DM halos is determined by the mass function of the accreting substructure (i.e., the initial power spectrum of initial perturbations). Dwarf galaxies have on average flatter DM cores, and clusters, steeper cores than galaxies similar to the Milky Way for which the NFW profile is a good fit to the halo density profile.
* The logarithmic slope, $`\gamma `$, of the outer parts of the halo appears to depend on the acceleration of the universe: when the scale parameter is $`a\stackrel{<}{}1`$, $`\gamma 3`$ as in the NFW profile, but $`\gamma 4`$ at $`a>1`$ when $`\mathrm{\Omega }_\mathrm{\Lambda }1`$ and the universe is inflating.
* A density profile in the form $`\rho _{dm}X^\alpha (1+X)^{\gamma +\alpha }`$, where $`X=c_\mathrm{\Delta }r/r_\mathrm{\Delta }`$, $`\alpha `$ is a function of the halo mass ($`\alpha 1.3+[M_{dm,14}^{1/6}1]/[M_{dm,14}^{1/6}+1]`$, where $`M_{dm,14}=M_{dm}/[3\times 10^{14}`$ M\]) and $`\gamma 3`$ (but $`\gamma 4`$ if $`\mathrm{\Omega }_\mathrm{\Lambda }1`$), can fit all the profiles from dwarfs to superclusters with constant concentration parameter $`c_\mathrm{\Delta }7`$. This provides a physical explanation for the variation of the concentration parameter in the NFW profile as a function of the halo mass: the changing slopes mimic the variation of the concentration parameter by shifting the radius where the circular velocity reaches its maximum value.
* The dependence of $`\alpha `$ on the mass of the halo, and therefore on the slope of the power spectrum, $`n`$, agrees with the theoretical expectation $`\alpha =(9+3n)/(5+n)`$, derived assuming that “undigested” satellites determine the halo structure in the inner regions (Subramanian et al., 2000). Note that if the DM halo develops a flat core, the chance for accreting satellites to survive tidal disruption is larger (Dekel et al., 2003).

The shape of the profiles for a given mass and redshift present significant statistical variations. This can be attributed to two factors. Firstly, galaxies are not isolated entities and especially at high redshift are severely disturbed by neighbouring satellites or by undergoing major mergers. This produces irregularities in the azimuthally averaged density profile. Secondly, the accretion history of satellites has statistical variation and depends on the local environment. A dwarf galaxy forming at high redshift from a large $`\sigma `$ peak of the power spectrum of initial perturbations (the type of galaxy simulated in this work) could have an inner profile different from the same mass galaxy forming at lower redshift from a smaller $`\sigma `$ perturbation. Also, the slope of the profile could differ if a galaxy forms in a void or in an overdense region, even for galaxies with similar masses and radii.

I have identified three main reasons to explain why previous studies (e.g., Eke et al., 2001; Moore et al., 1999) found that the inner slope of the DM profiles are not determined by the power spectrum or accretion history (but see also Syer & White (1998); Kravtsov et al. (1998); Jing & Suto (2000); Jing (2000) for a different point of view). (i) The profiles of halos with typical masses resolved by the 1 h-1 Mpc simulation have not been analysed before. (ii) The dynamical state of the the halos at $`z=10`$ might be different at $`z=0`$. But there are not published high-resolution simulations for dwarf-size halos ($`M_{DM}10^8`$ M) at $`z=0`$. The classic results of Cole & Lacey (1996); Navarro et al. (1997) are based on the same type of simulations presented here. (iii) Most of the computational efforts have focused on simulating a single galaxy-size or cluster-size halo with mass $`M_{dm}\stackrel{>}{}10^{12}`$ M at $`z=0`$ with higher resolution (about $`10^6`$ particles per halo) than in this study. Such a large number of particles is required to study the profile in the very inner regions of the halo. But the task of achieving reliable profiles with such a high spatial resolution is sensitive to numerical integration errors and requires careful resolution studies (Power et al., 2002); perhaps this is a reason for the disagreement between groups, using different codes, on the slope of the inner profile. In this work I find that the shape of the DM profiles present significant scatter around the mean. It is, therefore, important to analyse a significant statistical sample of halos in order to measure the mean profile. Works where a single halo is re-simulated with higher resolution are affected by a bias, difficult to control, determined by the criteria for picking the halos to re-simulate. Finally, the fitting method could be important as well for the results. The degeneracy between the slope of the inner profile and the value of the core radius is more easily broken fitting the circular velocities instead of the density profiles. This is because the radius where the circular velocity is maximum is univocally determined by the value of the core radius, $`r_s`$ (for the NFW profile is $`r_{max}=2.16r_s`$). The halo profiles at any redshift and mass, when renormalised matching the values and radii of the maximum circular velocities, have to be identical if their profiles can be all fitted by the NFW profile.

Moore et al. (1999) find steep profiles in cluster size halos ($`\alpha 1.5`$), in agreement with the results of this work. The results of this work also agree with the classic NFW result for Milky Way size galaxies but not with Moore et al. (1999) for halos with the same mass. In this work the spatial resolution is not as large as in the aforementioned studies, therefore I cannot study the slope of the DM profile in the very centre of the halos. Nevertheless, I find evidence for disagreement with the NFW predictions even at relatively large radii. Perhaps this new result has been found because it is the first time that a simulation of halos with masses typical of dwarf galaxies and mass resolution $`M_p10^3`$ M has been performed. The result for the cluster size halos could also be understood with similar reasoning (i.e., the cluster masses are larger than in previous studies: $`M_{dm}10^{15}10^{16}`$ M). The steeper core density profile found in clusters may not be as evident as the flatter core found in dwarf galaxies but, combined with the result for normal galaxies, the findings from the three simulations strengthen the case for an inner slope that changes with the mass of the halos.

I conclude this work with two final remarks:

* According to the scaling relation $`\alpha =(9+3n)/(5+n)`$, the inner slope of DM profile of dwarf galaxies is more sensitive to the power spectrum index, $`n`$, than more massive galaxies: $`\delta \alpha =1.5\delta n`$ if $`n3`$. A tilt $`\delta n0.2`$ in the initial power spectrum should produce a variation $`\delta \alpha 0.3`$ of the inner slope of the DM profile.
* The best place to look for the effects of the cosmological constant is the outer profile of clusters and superclusters at low redshift. The prediction is for a slope steeper than $`\gamma =3`$. Perhaps gravitational lensing studies could be useful in this regard.

### ACKNOWLEDGEMENTS

M.R. is supported by a PPARC theory grant. Research conducted in cooperation with Silicon Graphics/Cray Research utilising the Origin 3800 supercomputer (COSMOS) at DAMTP, Cambridge. COSMOS is a UK-CCC facility which is supported by HEFCE and PPARC. I would like to thank Jerry Ostriker for his valuable comments on the draft and exciting discussions.

## Appendix A Generalised Profile

The comoving virial radius, $`r_\mathrm{\Delta }`$, here is defined as the radius where the halo mean overdensity is $`\mathrm{\Delta }`$, i.e., $`M_{dm}(4\pi /3)\rho _0\mathrm{\Delta }r_\mathrm{\Delta }^3`$. Here $`M_{dm}`$ is the mass of the halo, $`\rho _0`$ is the mean DM density at $`z=0`$ and $`\mathrm{\Delta }178`$ if $`\mathrm{\Omega }(z)=1`$. A density profile with inner slope $`\alpha `$, and outer slope $`\gamma `$, can be written in the form

$$\rho _{dm}(r)=\frac{\rho _0\mathrm{\Delta }\delta _{(c,\alpha ,\gamma )}}{X^\alpha (1+X)^{\gamma \alpha }}$$

where $`X=r/r_s`$, $`r_s=r_\mathrm{\Delta }/c_\mathrm{\Delta }`$ is the core radius and $`\delta _{(c,\alpha ,\gamma )}`$ is a normalisation constant. The NFW profile has $`\alpha =1`$ and $`\gamma =3`$. The mass enclosed inside the radius, $`X`$, is

$$\begin{array}{c}\hfill m_{dm}(X)=4\pi \rho _0\mathrm{\Delta }\delta _{(c,\alpha ,\gamma )}r_s^3_0^X\frac{y^g}{(1+y)^{\gamma \alpha }}\frac{dy}{y}=\\ \hfill 4\pi \rho _0\mathrm{\Delta }\delta _{(c,\alpha ,\gamma )}r_s^3\left(\frac{X^g}{g}\right){}_{2}{}^{}F_{1}^{}(g;\gamma \alpha ;g+1;X)\end{array}$$
(3)
where $`g=(3\alpha )`$ and $`{}_{2}{}^{}F_{1}^{}`$ is the hypergeometric function. If $`\gamma =3`$, the hypergeometric function simplifies as

$${}_{2}{}^{}F_{1}^{}(g;g;g+1;X)=g\underset{n=0}{\overset{\mathrm{}}{}}\frac{(n;g)}{(g+n)n}(X)^n,$$

where $``$ is the beta function. When $`\alpha =1`$ and $`\gamma =3`$ (NFW profile) equation (3) can be integrated analytically and the solution is

$$m_{dm}(X)=\frac{4\pi }{3}\rho _0\mathrm{\Delta }\delta _c^{NFW}r_s^3\left[\mathrm{ln}(1+X)\frac{X}{1+X}\right].$$

The normalisation constant $`\delta _{(c,\alpha ,\gamma )}`$ is such that $`m_{dm}(X=c_\mathrm{\Delta })=M_{dm}`$. For the NFW profile, $`\delta _c^{NFW}=c_\mathrm{\Delta }^3/[\mathrm{ln}(1+c_\mathrm{\Delta })c_\mathrm{\Delta }/(1+c_\mathrm{\Delta })]`$. In the general case equation (3) has to be evaluated numerically. If $`\gamma =3`$, analytical solutions exist for $`\alpha =0,0.5,1,1.5`$ and 2:

$${}_{2}{}^{}F_{1}^{}(X)=\{\begin{array}{cc}\frac{3}{X^3}\left[\frac{(1+3X/2)X}{(1+X)^2}+\mathrm{ln}(1+X)\right],\text{if}\alpha =0\hfill & \\ \frac{5}{X^{5/2}}\left[\frac{(1+4X/3)X^{1/2}}{(1+X)^{3/2}}+\mathrm{ArcSinh}(X^{1/2})\right],\text{if}\alpha =0.5\hfill & \\ \frac{3}{X^{3/2}}\left[\left(\frac{X}{1+X}\right)^{1/2}+\mathrm{ArcSinh}(X^{1/2})\right]\text{if}\alpha =1.5\hfill & \\ \frac{\mathrm{ln}(1+X)}{X},\text{if}\alpha =2\hfill & \end{array}$$

If $`\gamma =4`$ an analytical solution exists for a generic $`\alpha `$:

$${}_{2}{}^{}F_{1}^{}(X)=\left(\frac{1+X}{X}\right)^\alpha [13(1+X)^1+3(1+X)^2(1+X)^3].$$

The circular velocity is defined as

$$V_c(X)=\sqrt{\frac{m_{dm}(X)}{Xr_s}}.$$
(4)
For the NFW profile equation (4) has a maximum $`V_c^{max}r_\mathrm{\Delta }`$ at $`r_{max}=2.16r_s`$. For a fixed value of $`r_\mathrm{\Delta }`$, $`V_c^{max}`$ depends weakly on $`c_\mathrm{\Delta }`$, and $`r_{max}`$ is inversely proportional to $`c_\mathrm{\Delta }`$. In the general case, $`r_{max}=f(\alpha ,\gamma )r_s`$ and must be calculated numerically. In Fig. 11(a) I show the numerical solution and a simple fit to $`f(\alpha ,\gamma )=r_{max}/r_s`$ as a function of $`\alpha `$ for $`\gamma =3,3.5`$ and 4. If $`\alpha =0`$ the maximum, $`f`$, of $`v_c`$ can be found solving the equation

$$2[1(1+f)^\gamma ]+(\gamma 1)f\left[(\gamma 2)^2f^2+\gamma f+\frac{2\gamma }{\gamma 1}\right]=0.$$

If $`\gamma =4`$ the solution is $`f=2`$, but in general the solution has to be evaluated numerically. The simplest fit to the numerical solution, shown in Fig. 11(b), is given by $`f=(4/\gamma )^{2.68}`$. Using this fit in conjunction with the previous for $`f(\alpha ,3)`$, a general fitting formula is found to be

$$\frac{r_{max}}{r_s}=f(\alpha ,\gamma )(4/\gamma )^{2.68}(2\alpha ),$$
(5)
which is at least 5 percent accurate in the parameter interval $`0<\alpha <2`$ and $`3<\gamma <4`$.


format: text

text:
# 1 Introductory remarks

## 1 Introductory remarks

At present time there are certain indications that formation and evolution of a significant part of black holes occur in active galactic nuclei (see, e .g. Ref.? and references therein). Under the circumstances the notion of black hole as an isolated object being described, for example, by the Schwarzschild (SW) metric evidently becomes inadequate to the situation. In active galactic nuclei the regions filled with a nonzero energy density or strong magnetic fields are the rule rather than exception. As a result, the properties of black holes may be changed under these conditions and to take into account this change one should modify the description of black holes. One can note that solutions of the Einstein equations which can model the mentioned conditions are known for a long time. In particular, the Kottler metric which can be written down in the form

$$ds^2=g_{\mu \nu }dx^\mu dx^\nu adt^2a^1dr^2r^2(d\vartheta ^2+\mathrm{sin}^2\vartheta d\phi ^2)$$
$`(1)`$
with $`a=12M/r+\mathrm{\Lambda }r^2/3`$ is the solution of the Einstein equations with cosmological constant $`\mathrm{\Lambda }`$ ($`M`$ is the black hole mass) and is known since 1918 (see, e. g., Refs.?). Inasmuch as the cosmological constant $`\mathrm{\Lambda }`$ can be associated with the nonzero energy density (in usual units) $`\rho =\mathrm{\Lambda }c^4/(8\pi G)`$ (see e. g. Ref. 3) then the metric (1) may be interpreted as the one describing the SW black hole immersed in a medium with the nonzero energy density. It should be noted, however, that metric (1) is not asymptotically flat one and the given circumstance can lead to difficulties of formulating the scattering problem necessary for the Hawking radiation analysis (see below).

On the other hand, the nontrivial topological properties of spacetimes may play essential role while studying quantum geometry of fields on them (see, e. g., our reviews?). The black hole spacetimes are not exceptions. Really, as was shown in a recent series of the papers ?,?,?,?, the black holes can carry the whole spectrum of topologically inequivalent configurations (TICs) for miscellaneous fields, in the first turn, complex scalar and spinor ones. The mentioned TICs can essentially modify the Hawking radiation from black holes? and also can give some possible statistical explanation of the black hole entropy?. Physically, the existence of TICs should be obliged to the natural presence of magnetic U(N)-monopoles (with $`N1`$) on black holes though the total (internal) magnetic charge (abelian or nonabelian) of black hole remains equal to zero.

In present paper we would like to clear up to what extent the mentioned features can be inherent to the black holes described by metric (1). For to simplify calculations our considerations will restrict themselves to the TICs of complex scalar field and, accordingly, to the Dirac U(1)-monopoles. Sec. 2 serves for obtaining the auxiliary relations necessary for the rest of paper. Sec. 3 is devoted to description of Dirac monopoles on the black holes under exploration and also to separation of variables for the corresponding wave equations while Sec. 4 contains the correct statement of the scattering problem for scalar particles interacting with the mentioned Dirac monopoles and description of the conforming $`S`$-matrix necessary to analyse the Hawking radiation process. The latter is discussed in Sec. 5 and Sec. 6 gives concluding remarks.

Throughout the paper we employ the system of units with $`\mathrm{}=c=G=1`$, unless explicitly stated. Finally, we shall denote $`L_2(F)`$ the set of the modulo square integrable complex functions on any manifold $`F`$ furnished with an integration measure.

## 2 Preliminary considerations

First of all it should be noted that the standard spacetime topology on which the metric (1) with arbitrary $`a=a(r)`$ can be realized in a natural way is of $`\mathrm{}^2\times \mathrm{𝕊}^2`$-form. In the given paper we shall, however, deal only with the Kottler metric where the form of $`a`$ was fixed above in Sec. 1. Let us introduce a number of quantities that for metric (1) with $`\mathrm{\Lambda }0`$ will generalize the corresponding ones when $`\mathrm{\Lambda }=0`$, i. e., in the SW black hole case. Besides we shall consider $`\mathrm{\Lambda }>0`$ within the paper.

Let us consider the equation $`a(r)=12M/r+\mathrm{\Lambda }r^2/3=0`$ that can be rewritten as $`r^3+3r/\mathrm{\Lambda }6M/\mathrm{\Lambda }=0`$. Discriminant of the latter equation is equal to $`D=(1/\mathrm{\Lambda })^3+(3M/\mathrm{\Lambda })^2>0`$ which means the mentioned equation to have one real and two complex conjugated roots. The only real root is

$$r_+=\frac{1}{\sqrt{\mathrm{\Lambda }}}(\sqrt[3]{\alpha }+\sqrt[3]{\beta })$$
$`(2)`$
with $`\alpha =3M\sqrt{\mathrm{\Lambda }}\sqrt{1+9M^2\mathrm{\Lambda }}`$, $`\beta =3M\sqrt{\mathrm{\Lambda }}+\sqrt{1+9M^2\mathrm{\Lambda }}`$ and we have the expansion

$$a(r)=\frac{\mathrm{\Lambda }}{3r}(rr_+)\left(r^2+r_+r+\frac{6M}{\mathrm{\Lambda }r_+}\right).$$
$`(3)`$
It is clear that $`r_+2M`$ when $`\mathrm{\Lambda }0`$ and we can see that horizon of the SW black hole is shifted to $`r_+`$ at $`\mathrm{\Lambda }>0`$. Under this situation when defining the black hole entropy as $`S=\pi r_+^2`$ we shall find the black hole temperature $`T`$ from the relation $`1/T=S/M`$ which entails

$$T=\frac{\mathrm{\Lambda }\sqrt{1+9M^2\mathrm{\Lambda }}(\alpha \beta )^{2/3}}{2\pi (\alpha ^{1/3}+\beta ^{1/3})[3M\mathrm{\Lambda }(\alpha ^{2/3}\beta ^{2/3})+\sqrt{\mathrm{\Lambda }(1+9M\mathrm{\Lambda }^2}(\alpha ^{2/3}+\beta ^{2/3})]},$$
$`(4)`$
so when $`\mathrm{\Lambda }0`$, $`T1/8\pi M`$, the SW black hole temperature.

Further, we define the variable $`r_{}`$ by the relation

$$r_{}=\frac{dr}{a}=\frac{r_+^2}{2(3Mr_+)}\mathrm{ln}\left(\frac{r}{r_+}1\right)\frac{r_+^2}{4(3Mr_+)}\mathrm{ln}\left[1+\frac{\mathrm{\Lambda }r_+}{6M}(r^2+rr_+)\right]$$

$$+\left[\frac{6M}{\mathrm{\Lambda }(3Mr_+)}+\frac{r_+^3}{2(3Mr_+)}\right]\frac{1}{\sqrt{\frac{24M}{\mathrm{\Lambda }r_+}r_+^2}}\mathrm{arctan}\frac{2r+r_+}{\sqrt{\frac{24M}{\mathrm{\Lambda }r_+}r_+^2}}\frac{r_+}{2},$$
$`(5)`$
where we choose the integration constant for the expression (5) to pass on to the standard Regge-Wheeler variable $`r_{}=r+2M\mathrm{ln}[r/(2M)1]`$ at $`\mathrm{\Lambda }0`$. At last, let us introduce the dimensionless quantities $`x=r_{}/M`$, $`y=r/M`$, $`y_+=r_+/M`$ so that for them we shall have

$$x=\frac{y_+^2}{2(3y_+)}\mathrm{ln}\left(\frac{y}{y_+}1\right)\frac{y_+^2}{4(3y_+)}\mathrm{ln}\left[1+\frac{\mathrm{\Lambda }M^2y_+}{6}(y^2+yy_+)\right]$$

$$+\left[\frac{6}{\mathrm{\Lambda }M^2(3y_+)}+\frac{y_+^3}{2(3y_+)}\right]\frac{1}{\sqrt{\frac{24}{\mathrm{\Lambda }M^2y_+}y_+^2}}\mathrm{arctan}\frac{2y+y_+}{\sqrt{\frac{24}{\mathrm{\Lambda }M^2y_+}y_+^2}}\frac{y_+}{2},$$
$`(6)`$
i. e, $`y=y(x)`$ is a reverse function to (6) and again at $`\mathrm{\Lambda }0`$ the relation is simplified to $`x=y+2\mathrm{ln}(0.5y1)`$. At this point we face the difficulty that at $`yy_+`$, $`x\mathrm{}`$ while when $`y+\mathrm{}`$, $`x`$ some finite value which is different from the SW black hole case, where $`y_+=2y+\mathrm{}`$, whereas $`\mathrm{}<x<+\mathrm{}`$. Under this situation the correct statement of the scattering problem necessary for the Hawking radiation analysis is impossible so long as the latter problem should be considered on the whole $`x`$-axis?.

To overcome this obstacle it should be noted that the region in universe with the nonzero energy density (and where the SW black hole can be located) has, generally speaking, some finite sizes. Under the circumstances, if $`r_0`$ is such a characteristic size then the metric (1) should describe the black hole at $`r_+rr_0`$, while at $`r_0r+\mathrm{}`$ the black hole should be described by metric (1) with $`\mathrm{\Lambda }=0`$, i. e., by the standard SW metric. As a result, the relation (6) should be applied when $`\mathrm{}<xx_0`$, otherwise, at $`x_0x<+\mathrm{}`$ the corresponding $`y(x)`$ should be determined from the equation $`x=y+2\mathrm{ln}(0.5y1)`$ and one should put $`y_0=r_0/M`$ while $`x_0=y_0+2\mathrm{ln}(0.5y_01)`$.

Now we shall have a function $`y(x)`$ defined on the whole $`x`$-axis and the conforming scattering problem can be posed correctly (see Sec. 4). It is clear the size $`r_0`$ should be specified under the concrete conditions so in what follows we shall use it as a parameter in further considerations.

## 3 Dirac monopoles and separation of variables

As was disscussed in Refs.?,?,?, the black holes can carry TICs of a complex scalar field that are conditioned by the availability of a countable number of complex line bundles over the $`\mathrm{}^2\times \mathrm{𝕊}^2`$-topology underlying the 4D black hole physics. Each TIC corresponds to sections of a complex line bundle $`E`$ while the latter can be characterized by its Chern number $`n\mathrm{}`$ (the set of integers). TIC with $`n=0`$ can be called untwisted, while the rest of the TICs with $`n0`$ should be referred to as twisted. Using the fact that all the mentioned line bundles can be trivilized over the chart of local coordinates $`(t,r,\vartheta ,\phi )`$ covering almost the whole manifold $`\mathrm{}^2\times \mathrm{𝕊}^2`$one can obtain a suitable wave equation on the given chart for TIC $`\varphi `$ with mass $`\mu _0`$ and Chern number $`n\mathrm{}`$ in the form

$$|g|^{1/2}(_\mu ieA_\mu )[g^{\mu \nu }(|g|)^{1/2}(_\nu ieA_\nu )\varphi ]=\mu _0^2\varphi ,$$
$`(7)`$
where the specification depends on the gauge choice for connection $`A_\mu `$ (vector-potential for the corresponding Dirac monopole) in the line bundle $`E`$ with the Chern number $`n`$. We shall here consider the gauge where one can put the above connection 1-form equal to $`A=A_\mu dx^\mu =\frac{n}{e}\mathrm{cos}\vartheta d\phi `$. Then the curvature of the bundle $`E`$ is $`F=dA=\frac{n}{e}\mathrm{sin}\vartheta d\vartheta d\phi `$. We can further introduce the Hodge star operator on 2-forms $`F`$ of any $`k`$-dimensional (pseudo)riemannian manifold $`B`$ provided with a (pseudo)riemannian metric $`g_{\nu \tau }`$ by the relation (see, e. g., Ref.?)

$$FF=(g^{\nu \alpha }g^{\tau \beta }g^{\nu \beta }g^{\tau \alpha })F_{\nu \tau }F_{\alpha \beta }\sqrt{|g|}dx^1dx^2\mathrm{}dx^k$$
$`(8)`$
in local coordinates $`x^\nu `$. In the case of the metric (1) this yields

$$(dtdr)=\sqrt{|g|}g^{tt}g^{rr}d\vartheta d\phi =r^2\mathrm{sin}\vartheta d\vartheta d\phi ,$$

$$(d\vartheta d\phi )=\sqrt{|g|}g^{\vartheta \vartheta }g^{\phi \phi }dtdr=\frac{1}{r^2\mathrm{sin}\vartheta }dtdr,$$
$`(9)`$
so that $`^2==1`$, as should be for the manifolds with lorentzian signature?. After this we find $`F=\frac{n}{er^2}dtdr,`$ and integrating over the surface $`t=const`$, $`r=const`$ with topology $`\mathrm{𝕊}^2`$ gives rise to the Gauss theorem

$$_{S^2}F=0$$
$`(10)`$
and the Dirac charge quantization condition

$$_{S^2}F=4\pi \frac{n}{e}=4\pi q$$
$`(11)`$
with magnetic charge $`q`$. Besides, the Maxwell equations $`dF=0`$, $`dF=0`$ are clearly fulfilled with the exterior differential $`d=_tdt+_rdr+_\vartheta d\vartheta +_\phi d\phi `$ in coordinates $`t,r,\vartheta ,\phi `$. Also it should be emphasized that the total (internal) magnetic charge $`Q_m`$ of black hole which should be considered as the one summed up over all the monopoles remains equal to zero because

$$Q_m=\frac{1}{e}\underset{n\mathrm{}}{}n=0,$$
$`(12)`$
so the external observer does not see any magnetic charge of the black hole though the monopoles are present on black hole in the sense described above.

### 3.1 Monopole masses

We can estimate the corresponding monopole masses so long as we are effectively in the asymptotically flat spacetime due to the remarks made late in Sec. 2 and one can use the $`T_{00}`$-component of the energy-momentum tensor for electromagnetic field

$$T_{\mu \nu }=\frac{1}{4\pi }(F_{\mu \alpha }F_{\nu \beta }g^{\alpha \beta }+\frac{1}{4}F_{\beta \gamma }F_{\alpha \delta }g^{\alpha \beta }g^{\gamma \delta }g_{\mu \nu })$$
$`(13)`$
to find the monopole masses according to

$$m_{\mathrm{mon}}(n)=\underset{\mathrm{}\times \mathrm{𝕊}^2}{}T_{00}\sqrt{\gamma }d^3x=\underset{\mathrm{}\times \mathrm{𝕊}^2}{}T_{00}r^2\mathrm{sin}\vartheta a^{1/2}d^3x$$
$`(14)`$
with $`a(r)`$ of (1) and $`\sqrt{\gamma }d^3x=r^2\mathrm{sin}\vartheta a^{1/2}d^3x`$ is the volume element for the spatial part (with topology $`\mathrm{}\times \mathrm{𝕊}^2`$) of the metric (1) while $`T_{00}`$-component is proved to be equal to $`\frac{1}{16\pi }aF_{\vartheta \phi }^2g^{\vartheta \vartheta }g^{\phi \phi }`$ with $`F_{\vartheta \phi }=n\mathrm{sin}\vartheta /e`$. This entails

$$m_{\mathrm{mon}}(n)=\frac{n^2}{4e^2}\left[\underset{r_+}{\overset{r_0}{}}\sqrt{1\frac{2M}{r}+\frac{\mathrm{\Lambda }r^2}{3}}\frac{dr}{r^2}+\underset{r_0}{\overset{\mathrm{}}{}}\sqrt{1\frac{2M}{r}}\frac{dr}{r^2}\right]$$

$$=\frac{n^2}{4e^2}\left[\underset{r_+}{\overset{r_0}{}}\sqrt{1\frac{2M}{r}+\frac{\mathrm{\Lambda }r^2}{3}}\frac{dr}{r^2}+\frac{1}{3M}\frac{1}{3M}\left(1\frac{2M}{r_0}\right)^{3/2}\right].$$
$`(15)`$
It is clear that at $`r_0+\mathrm{}`$, $`m_{\mathrm{mon}}(n)\mathrm{}`$ and the conforming Compton wavelength $`\lambda _{\mathrm{mon}}=1/m_{\mathrm{mon}}0`$, i. e., $`\lambda _{\mathrm{mon}}<<r_+`$ and the given monopoles might reside in black hole under discussion as quantum objects.

### 3.2 Separation of variables

If returning to the Eq. (7), then, in the way analogous to that of Refs.?,?,?, one can show the equation (7) to have in $`L_2(`$$`\mathrm{}^2\times \mathrm{𝕊}^2`$$`)`$ a complete set of solutions of the form

$$f_{\omega nml}=\frac{(2\pi \omega )^{1/2}}{r}e^{i\omega t}Y_{nlm}(\vartheta ,\phi )R_{\omega nl}(r),l=|n|,|n|+1,\mathrm{},|m|l,$$
$`(16)`$
where the explicit form of the monopole (spherical) harmonics $`Y_{nlm}(\vartheta ,\phi )`$ can be found in Refs. ?,?,?. As to the function $`R_{\omega nl}(r)`$, then it obeys the equation

$$a_rar^2\frac{R_{\omega nl}}{r}+\omega ^2rR_{\omega nl}=a\left[\mu _0^2+\frac{l(l+1)n^2}{r^2}\right]rR_{\omega nl},$$
$`(17)`$
which can be rewritten as

$$\frac{d^2R_{\omega nl}}{dr_{}^2}+(\omega ^2\mu _0^2)R_{\omega nl}=\left[\mu _0^2(a1)+\frac{a}{r}\frac{da}{dr}+a\frac{l(l+1)n^2}{r^2}\right]R_{\omega nl}$$
$`(18)`$
with the help of (5). At last, when denoting $`k=M\omega `$ and introducing the functions $`\psi (x,k,n,l)=R_{\omega nl}[My(x)]`$, where $`y(x)`$ is a function reverse to (6), we shall obtain that $`\psi (x,k,n,l)`$ obeys the Schrödinger-like equation

$$\left[\frac{d^2}{dx^2}+(k^2\mu _0^2M^2)\right]\psi (x,k,n,l)=V(x,M,\mathrm{\Lambda },n,l)\psi (x,k,n,l),$$
$`(19)`$
where

$$V(x,M,\mathrm{\Lambda },n,l)=\left[1\frac{2}{y(x)}+\frac{\mathrm{\Lambda }M^2}{3}y^2(x)\right]\left[\frac{2}{y^3(x)}+\frac{2\mathrm{\Lambda }M^2}{3}+\frac{l(l+1)n^2}{y^2(x)}\right]$$

$$+(\mu _0M)^2\left[\frac{\mathrm{\Lambda }M^2}{3}y^2(x)\frac{2}{y(x)}\right],\mathrm{}xx_0,$$
$`(20)`$
and

$$V(x,M,\mathrm{\Lambda },n,l)=\left[1\frac{2}{y(x)}\right]\left[\frac{2}{y^3(x)}+\frac{l(l+1)n^2}{y^2(x)}\right](\mu _0M)^2\frac{2}{y(x)},x_0x+\mathrm{}.$$
$`(21)`$
Now we can discuss how to pose the correct scattering problem for Eq. (19) because the corresponding solutions of this problem will be just necessary to analyse the Hawking radiation from black hole.

## 4 The scattering problem

One can notice that potential $`V(x,M,\mathrm{\Lambda },n,l)`$ of (20)–(21) satisfies the conditions explored in Refs.? and referring for more details to those works we may here use their results to correctly formulate the scattering problem for Eq. (19). Namely, the correct statement of the scattering problem will consist in searching for two solutions $`\psi ^+(x,k^+,n,l)`$, $`\psi ^{}(x,k,n,l)`$ of the equation (19) obeying the following conditions

$$\psi ^+(x,k^+,n,l)=\{\begin{array}{cc}e^{ikx}+s_{12}(k,l,n)e^{ikx}+o(1),\hfill & x\mathrm{}\text{,}\hfill \\ s_{11}(k,l,n)w_{i\mu ,\frac{1}{2}}(2ik^+x)+o(1),\hfill & x+\mathrm{}\text{,}\hfill \end{array}$$

$$\psi ^{}(x,k,n,l)=\{\begin{array}{cc}s_{22}(k,l,n)e^{ikx}+o(1),\hfill & x\mathrm{}\text{,}\hfill \\ w_{i\mu ,\frac{1}{2}}(2ik^+x)+s_{21}(k,l,n)w_{i\mu ,\frac{1}{2}}(2ik^+x)+o(1),\hfill & x+\mathrm{}\text{,}\hfill \end{array}$$
$`(22)`$
where $`k^+(k)=\sqrt{k^2(\mu _0M)^2}`$, $`\mu =\frac{(\mu _0M)^2}{(k^+)^2}=\frac{(\mu _0M)^2}{k^2(\mu _0M)^2}`$ and the functions $`w_{\pm i\mu ,\frac{1}{2}}(\pm z)`$ are related to the Whittaker functions $`W_{\pm i\mu ,\frac{1}{2}}(\pm z)`$ (concerning the latter ones see e. g. Ref.?) by the relation

$$w_{\pm i\mu ,\frac{1}{2}}(\pm z)=W_{\pm i\mu ,\frac{1}{2}}(\pm z)e^{\pi \mu /2},$$

so that one can easily gain asymptotics (using the corresponding ones for Whittaker functions?)

$$w_{i\mu ,\frac{1}{2}}(2ik^+x)=e^{ik^+x}e^{i\mu \mathrm{ln}|2k^+x|}[1+O(|k^+x|^1)],x+\mathrm{},$$

$$w_{i\mu ,\frac{1}{2}}(2ik^+x)=e^{ik^+x}e^{i\mu \mathrm{ln}|2k^+x|}[1+O(|k^+x|^1)],x+\mathrm{}.$$
$`(23)`$
We can see that there arises some $`S`$-matrix with elements $`s_{ij},i,j=1,2`$ and it satisfies certain unitarity relations (for more details see Refs.?) As will be seen below, for calculating the Hawking radiation we need the coefficient $`s_{11}`$, consequently, we need to have some algorithm for numerical computation of it inasmuch as the latter cannot be evaluated in exact form. The given algorithm can be extracted from the results of Refs.? To be more precise

$$s_{11}(k,l,n)=2ik/[f^{}(x,k,l,n),f^+(x,k^+,l,n)],$$
$`(24)`$
where \[,\] signifies the Wronskian of functions $`f^{},f^+`$, the so-called Jost type solutions of Eq. (19). In their turn, these functions and their derivatives obey the certain integral equations. Since the Wronskian does not depend on $`x`$ one can take the following form of the mentioned integral equations

$$f^{}(x_0,k,l,n)=e^{ikx_0}+\frac{1}{k}\underset{\mathrm{}}{\overset{x_0}{}}\mathrm{sin}[k(x_0t)]V^{}(t,l,n)f^{}(t,k,l,n)𝑑t,$$
$`(25)`$
$$(f^{})_x^{}(x_0,k,l,n)=ike^{ikx_0}+\underset{\mathrm{}}{\overset{x_0}{}}\mathrm{cos}[k(x_0t)]V^{}(t,l,n)f^{}(t,k,l,n)𝑑t,$$
$`(26)`$
$$f^+(x_0,k^+,l,n)=w_{i\mu ,\frac{1}{2}}(2ik^+x_0)+$$

$$\frac{1}{k^+}\underset{x_0}{\overset{+\mathrm{}}{}}\mathrm{Im}[w_{i\mu ,\frac{1}{2}}(2ik^+x_0)w_{i\mu ,\frac{1}{2}}(2ik^+t)]V^+(t,l,n)f^+(t,k^+,l,n)𝑑t,$$
$`(27)`$
$$(f^+)_x^{}(x_0,k^+,l,n)=\frac{d}{dx}w_{i\mu ,\frac{1}{2}}(2ik^+x_0)+$$

$$\frac{1}{k^+}\underset{x_0}{\overset{+\mathrm{}}{}}\mathrm{Im}[\frac{d}{dx}w_{i\mu ,\frac{1}{2}}(2ik^+x_0)w_{i\mu ,\frac{1}{2}}(2ik^+t)]V^+(t,l,n)f^+(t,k^+,l,n)𝑑t,$$
$`(28)`$
where the potentials

$$V^{}(x,l,n)=\left[1\frac{2}{y(x)}+\frac{\mathrm{\Lambda }M^2}{3}y^2(x)\right]\left[(\mu _0M)^2+\frac{2}{y^3(x)}+\frac{2\mathrm{\Lambda }M^2}{3}+\frac{l(l+1)n^2}{y^2(x)}\right],$$
$`(29)`$
$$V^+(x,l,n)=2(\mu _0M)^2\left[\frac{1}{x}\frac{1}{y(x)}\right]+\left[1\frac{2}{y(x)}\right]\left[\frac{l(l+1)n^2}{y^2(x)}+\frac{2}{y^3(x)}\right]$$
$`(30)`$
tend to 0 when, respectively, $`x\mathrm{}`$ or $`x+\mathrm{}`$. The relations (24)–(30) can be employed for numerical calculation of $`s_{11}`$ ( see Refs.? for the case $`\mathrm{\Lambda }=0`$).

To summarize, we can now obtain the general solution of (19) in the class of functions restricted on the whole $`x`$-axis in the linear combination form

$$\psi (x,k,n,l)=C_{nl}^1(k)\psi ^+(x,k^+,n,l)+C_{nl}^2(k)\psi ^{}(x,k,n,l),$$
$`(31)`$
since a couple of the functions $`\psi ^+(x,k^+,n,l)`$, $`\psi ^{}(x,k,n,l)`$ forms a fundamental system of solutions for the equation (19).

Having obtained all the above, we can discuss the Hawking radiation process for any TIC of complex scalar field.

## 5 Modification of Hawking radiation

As can be seen, the equation (7) corresponds to the lagrangian

$$=(g)^{1/2}(g^{\mu \nu }\overline{𝒟_\mu \varphi }𝒟_\nu \varphi \mu _{0}^{}{}_{}{}^{2}\overline{\varphi }\varphi ),$$
$`(32)`$
where the overbar signifies complex conjugation, $`𝒟_\mu =_\mu ieA_\mu `$ and, as a result, we have the conforming energy-momentum tensor for TIC with the Chern number $`n`$

$$T_{\mu \nu }=\mathrm{Re}[(\overline{𝒟_\mu \varphi })(𝒟_\nu \varphi )\frac{1}{2}g_{\mu \nu }g^{\alpha \beta }(\overline{𝒟_\alpha \varphi })(𝒟_\beta \varphi )\mu _0^2g_{\mu \nu }\overline{\varphi }\varphi ].$$
$`(33)`$
When quantizing twisted TIC with the Chern number $`n`$ we shall take the set of functions $`f_{\omega nlm}`$ of (16) as a basis in $`L_2`$($`\mathrm{}^2\times \mathrm{𝕊}^2`$) so that $`|m|l`$, $`l=|n|,|n|+1`$,… and we shall normalize the monopole spherical harmonics $`Y_{nlm}(\vartheta ,\phi )`$ to 1, i. e., by the condition

$$\underset{0}{\overset{\pi }{}}\underset{0}{\overset{2\pi }{}}\overline{Y_{nlm}}(\vartheta ,\phi )Y_{nl^{}m^{}}(\vartheta ,\phi )\mathrm{sin}\vartheta d\vartheta d\phi =\delta _{ll^{}}\delta _{mm^{}}.$$
$`(34)`$
After this we can evidently realize the procedure of quantizing TIC with the Chern number $`n`$, as usual, by expanding this TIC in the modes $`f_{\omega nlm}`$

$$\varphi =\underset{l=|n|}{\overset{\mathrm{}}{}}\underset{|m|l}{}\underset{\mu _0}{\overset{\mathrm{}}{}}𝑑\omega (a_{\omega nlm}^{}f_{\omega nlm}+b_{\omega nlm}^+\overline{f_{\omega nlm}}),$$

$$\varphi ^+=\underset{l=|n|}{\overset{\mathrm{}}{}}\underset{|m|l}{}\underset{\mu _0}{\overset{\mathrm{}}{}}𝑑\omega (b_{\omega nlm}^{}f_{\omega nlm}+a_{\omega nlm}^+\overline{f_{\omega nlm}}),$$
$`(35)`$
so that $`a_{\omega nlm}^\pm `$, $`b_{\omega nlm}^\pm `$ should be interpreted as the corresponding creation and annihilation operators for charged scalar particle in both the gravitational field of black hole and the field of the conforming monopole with the Chern number $`n`$ and we have the standard commutation relations

$$[a_i^{},a_j^+]=\delta _{ij},[b_i^{},b_j^+]=\delta _{ij}$$
$`(36)`$
and zero for all other commutators, where $`i=\{\omega nlm\}`$ is a generalized index. Under the circumstances one can speak about the Hawking radiation process for any TIC of complex scalar field. Actually, one should use the energy-momentum tensor of (33) for TIC with the Chern number $`n`$ in quantum form

$$T_{\mu \nu }=\mathrm{Re}[(\overline{𝒟_\mu }\varphi ^+)(𝒟_\nu \varphi )\frac{1}{2}g_{\mu \nu }g^{\alpha \beta }(\overline{𝒟_\alpha }\varphi ^+)(𝒟_\beta \varphi )].$$
$`(37)`$
To get the luminosity $`L(n)`$ with respect to the Hawking radiation for TIC with the Chern number $`n`$, we define a vacuum state $`|0>`$ by the conditions

$$a_i^{}|0>=0,b_i^{}|0>=0,$$
$`(38)`$
and then

$$L(n)=\underset{r\mathrm{}}{lim}\underset{S^2}{}<0|T_{tr}|0>𝑑\sigma $$
$`(39)`$
with the vacuum expectation value $`<0|T_{tr}|0>`$ and the surface element $`d\sigma =r^2\mathrm{sin}\vartheta d\vartheta d\phi `$. Now with the help of (34), (35), (37), (38) and passing on to the quantities $`x`$, $`k`$ introduced in Sec. 3, we shall find

$$L(n)=M^2\underset{x\mathrm{}}{lim}\mathrm{Re}\left[\frac{i}{2\pi }\underset{l=|n|}{\overset{\mathrm{}}{}}(2l+1)_{\mu _0M}^{\mathrm{}}\overline{(_x\psi (x,k,n,l)}\psi (x,k,n,l)𝑑k\right],$$
$`(40)`$
where the relations

$$\frac{d}{dr}=\frac{1}{a}\frac{d}{dr^{}}=\frac{1}{aM}\frac{d}{dx}$$
$`(41)`$
were used.

As a consequence, the choice of vacuum state can be defined by the choice of the suitable linear combination of (31) for $`\psi (x,k,n,l)`$. Under this situation the Hawking radiation corresponds to the mentioned choice in the form

$$C_{nl}^1(k)=\frac{1}{\sqrt{\mathrm{exp}(\frac{k^+}{TM})1}},C_{nl}^2(k)=0,$$
$`(42)`$
which at last with using asymptotics (22)–(23) entails (in usual units)

$$L(n)=A\underset{l=|n|}{\overset{\mathrm{}}{}}(2l+1)\underset{\stackrel{~}{\mu }}{\overset{\mathrm{}}{}}\frac{|s_{11}(k,n,l)|^2k^+dk}{\mathrm{exp}(\frac{k^+}{TM})1},$$
$`(43)`$
where the black hole temperature $`T`$ is given by (4), $`A=\frac{1}{2\pi \mathrm{}}\left(\frac{\mathrm{}c^3}{GM}\right)^20.27367310^{50}\mathrm{erg}\mathrm{s}^1M^2`$, $`\stackrel{~}{\mu }=(\mu _0M)/\mu _{pl}^2`$, $`\mu _0`$ and $`M`$ in g while $`\mu _{pl}^2=c\mathrm{}/G`$ is the Planck mass square in g2.

We can interpret $`L(n)`$ as an additional contribution to the Hawking radiation due to the additional charged scalar particles leaving the black hole because of the interaction with monopoles and the conforming radiation can be called the monopole Hawking radiation?,?. Under this situation, for the all configurations luminosity $`L`$ of black hole with respect to the Hawking radiation concerning the complex scalar field to be obtained, one should sum up over all $`n`$, i. e.

$$L=\underset{n\mathrm{}}{}L(n)=L(0)+2\underset{n=1}{\overset{\mathrm{}}{}}L(n)$$
$`(44)`$
since $`L(n)=L(n)`$.

As a result, we can expect a marked increase of Hawking radiation from black holes under consideration. But for to get an exact value of this increase one should apply numerical methods. In the case of the pure SW black hole ($`\mathrm{\Lambda }=0`$), for example, it was found that the contribution due to monopoles can be of order 11 % of the total pion-kaon luminosity?. So that it would be interesting enough to evaluate similar increase in the case $`\mathrm{\Lambda }0`$.

## 6 Concluding remarks

The results of the present paper show that nontrivial topological properties of black holes may have essential influence on their physical properties even when the black holes cannot be considered as the isolated objects. Though we dwelt only upon the scalar particle case it is beyond doubt that similar considerations will hold true for spinor particles too if employing the results of Refs.? The correct statement of the corresponding scattering problem for spinor particles on the SW black hole has been only recently explored? and one can extend the results obtained here to the spinor case. Also the important case is that of black hole immersed in strong magnetic field. The conforming metric is the Schwarzschild-Ernst one? which is defined on the topology $`\mathrm{}^2\times \mathrm{𝕊}^2`$ as well, consequently, one should investigate how the (internal) Dirac monopoles will behave under these conditions. We hope to continue the work along this line in future explorations.

## 7 Acknowlegments

The work of authors was supported in part by the Russian Foundation for Basic Research (grant no. 01-02-17157).

References