text:
# The shape of the CMB power spectrum.

## 1 Introduction

In recent years high resolution measurements of cosmic microwave background (CMB) anisotropies have measured the CMB power spectrum beyond the first acoustic peak. More recently the all-sky WMAP results represent a new step towards high precision cosmological parameter estimation. When combining these measurements with large-scale structure data and distant supernovae an increasingly consistent picture of the universe arises.

In order to investigate cosmological parameters in detail the experiments themselves have to be thoroughly understood. CMB anisotropy experiments are essentially of two kinds. High frequency experiments use bolometers and are usually balloon-born. The calibration, pointing and beam shape are the main sources of uncertainty of such instruments and tend to smear features in the power spectrum. Low frequency experiments use interferometers and direct;y measure the CMB in Fourier space. The calibration of such instruments is usually well understood and they do not suffer from beam uncertainty but the band-powers inferred from the observations tend to be anti-correlated. This can lead to increased contrast in the bands and is usually well encoded in the correlation matrix of the band-powers.

This work is an extension to a previous project presented in odman . We use our previously obtained results as an input to a Markov Chain Monte-Carlo (MCMC) program to estimate the shape of the CMB anisotropy power spectrum using a phenomenological curve. We use the WMAP data along with other CMB data separated in three subsets according to the type of experiment.

## 2 Method

In recent cosmological parameter estimations based on CMB data (see melchpars  and references therein) it has been shown that the simplest theoretical models based on $`6`$ fundamental parameters are sufficient to describe the universe from currently available data. Common assumptions are a flat universe in which adiabatic initial perturbations evolve with no contribution from background gravity waves or massive neutrinos.

With the recently released WMAP data a more complicated cosmological picture can be tested spergel  including a varying spectral index, an equation of state of dark energy etc. It is has now become common to test a parameter space of $`11`$ dimensions or more. The shape of the CMB power spectrum allowed in such analyses is limited by physical models. It is therefore timely to investigate what model-independent shape is preferred by the data in order to test the influence of experimental systematics which might not be fully under control.

In this analysis, we parametrise the CMB power spectrum using a sum of five consecutive Gaussians

$$\mathrm{}(\mathrm{}+1)C_{\mathrm{}}/2\pi =\underset{i=1}{\overset{N}{}}\mathrm{\Delta }T_i^2\mathrm{exp}((\mathrm{}\mathrm{}_i)^2/2\sigma _i^2),N=5.$$
(1)
We then use an MCMC routine to constrain the fifteen parameters of our phenomenological curve; $`\mathrm{\Delta }T_i`$, $`l_i`$ and $`\sigma _i`$. MCMC is an iterative process. After a certain burn-in time during which the Markov chains reach equilibrium, the parameters are sampled such that their number density is proportional to their posterior distribution. As explained in verde  such a routine is sensitive to the prior probability assumed for the parameters. Its efficiency also depends on the sampling strategy. If degeneracies exist among parameters it is useful to sample along degeneracy directions.

We make use of previously obtained results as an input to our MCMC routine for both the prior and the sampling strategy. We assume our previous results represent a $`15`$-dimensional Gaussian which we use a prior. For the sampling strategy, we diagonalise the covariance matrices obtained from our analysis before the WMAP data release. We then sample the parameter space as follows:

$$P_i=P_i^0+\underset{j}{}C_{ij}\lambda _jp_j,i,j=1,\mathrm{},15,$$
(2)
where $`P_i`$ are the parameters of our Gaussians, $`\mathrm{\Delta }T_i,l_i`$ and $`\sigma _i`$. $`P_i^0`$ are the previously obtained best-fit values of those parameters, $`\lambda _i`$ are the eigenvalues of the covariance matrices, $`C_{ij}`$ is the matrix of eigenvectors and $`p_j`$ are sampled from a Gaussian of mean zero and unit variance. This is equivalent to following the axes of our Gaussian prior. We further check that the samples are within the previously assumed priors.

We then run our MCMC code for three subsets of data. We first separate high (HF) and low frequency (LF) data as listed in table 1 and then combine them. We always use the WMAP data. For each sample we calculate the phenomenological power spectrum and return its likelihood to the MCMC routine. The WMAP likelihood is calculated using the method described in verde . For the other data the likelihood for a given phenomenological model is defined by $`2\mathrm{l}\mathrm{n}=(C_B^{ph}C_B^{exp})M_{BB^{}}(C_B^{}^{ph}C_B^{}^{exp})`$ where $`M_{BB^{}}`$ is the Gaussian curvature of the likelihood matrix at the peak. When available we use the lognormal approximation to the band-powers. We marginalize over the reported Gaussian distributed calibration error for each experiment and we include the beam uncertainties by the analytical marginalization method presented in sara .

## 3 Results and discussion

We run our code until we obtain around $`30,000`$ samples for each run. The number density of each parameter is proportional to its marginalised posterior distribution. We summarize our resuls in table 2 along with the pre-WMAP results for the fifteen Gaussian parameters of our phenomenological power spectra. The first peak is consistenly higher when including WMAP. It is shifted to slightly higher multipoles. The second peak is more pronounced than pre-WMAP when including all data. The third peak is very high for the LF data. This is mainly due to the VSAE band power at $`\mathrm{}800`$. At high multipoles, no distinct fourth or fifth peaks are detected, but the damping of the power is present in all cases.

In order to get the shape of the power spectrum we divide the $`\mathrm{}`$ \- $`\mu K^2`$ space into pixels of size $`\mathrm{\Delta }l=10`$, $`\mathrm{\Delta }T=10`$ and count the number of phenomenological spectra that go through each pixel. The resulting allowed shapes are shown in figure 2 along with pre-WMAP results. This confirms the consistency and accuracy of the constraints obtained on the first two peaks. At smaller angular scales, as we had obtained before, the discrepancy remains between high and low frequency experiments. Even if we include correlations and marginalise over calibration and beam uncertainties in the data the high frequency data show less contrasted features in the power spectrum than low frequency experiments.

In particular the third peak is higher than the second peak for low frequency experiments. This has important consequences in terms of the estimate of cold dark matter in the universe as shown in figure 3. The yellow region in figure 3 which corresponds to the low frequency experiments intersects the orange regions for higher values of the CDM physical density than from the WMAP parameter estimation (vertical gray region). The ratio of the first to the third peak from high-frequency data allows for a lower CDM density owing to the less marked secondary features in the power spectrum.

## 4 Conclusions

In this work, we take advantage of the power of the WMAP data to constrain the shape of the CMB power spectrum. We find that when including data from both high and low frequency experiments the preferred shape of the power spectrum is very close to the shape of the physical ‘concordance model’. Low frequency data show very contrasted second and third peaks. High frequency data prefer more shallow features. In all cases a damping of the power on small scales is detected.

Such a model-independent approach sheds light on discrepancies between high and low frequency experiments which would be less marked in traditional parameter estimation techniques. Our best fit shapes are not physical but the main observables of the spectrum are consistent with the WMAP cosmological parameter estimations.

Further investigation of the experiment-related discrepancies would be useful and WMAP provides a strong starting point. In combination with higher multipole data the shape of the CMB power spectrum is now strongly constrained. There is no doubt that a second and third peaks have been detected. This illustrates the new era cosmology has entered: the question now is not whether secondary peaks are present, thereby strengthening the already compelling agreement with inflationary predictions but how strong those secondary features are.

We would like to thank Alessandro Melchiorri, Anthony Lasenby, Licia Verde, Hiranya Peiris and Rob Izzard for useful comments. CJÖ is supported by an Isaac Newton Studentship and a Girton College Scholarship.


format: text

text:
# Goodness–of–fit Statistics and CMB Data Sets

## 1 Introduction

Measurement of the cosmic microwave background (CMB) temperature anisotropies has proven to be one of the most powerful tools for estimating important cosmological parameters (Netterfield et al. 2002 , Pryke et al. 2002 , Rubino-Martin et al. 2002 , Sievers et al. 2002 , Wang et al. 2002 ). The observed angular power spectrum shows the coherent peak structure expected in inflationary models, and fitting model curves to the data111See http://webast.ast.obs-mip.fr/cosmo/CMB for an up–to–date compilation. yields constraints on many parameters. This leads in particular to the conclusion that the geometry of space is flat (Lineweaver et al. 1997 , de Bernardis et al. 2000 , Hanany et al. 2000 , Lange et al. 2000 , Balbi et al. 2000 ). In terms of statistics, the procedure just described is one of parameter estimation.

Parameter estimation proceeds via the identification of a best model (set of parameters) within a family of models, an evaluation of the quality of the fit and the construction of parameter constraints. The method of maximum likelihood (ML), for example, is a useful, general procedure for finding a best–fit model. As a general rule, one must judge the quality of the fit before any serious consideration of parameter constraints. This requires the application of a Goodness–of–fit (GOF) statistic. Such a statistic is, usually, some scalar function of the data whose distribution may be calculated once given an underlying physical model and a model of the statistical fluctuations in the data. It is generally a function $`gof(d,T)`$ of both the data $`d`$ and theory $`T`$, such that $`gof`$ attains, for example, a minimum when $`d`$ is generated by the theory $`T`$. It is defined in a ’monotonic’ way, in the sense that $`gof`$ becomes larger as $`d`$ gets ’further’ from a realization of $`T`$. The ’significance’ may then be defined as the probability of obtaining $`gof>gof_{\mathrm{obs}}`$. On this basis, it permits a quantitative evaluation of the quality of the best model’s fit to the data: if the probability of obtaining the observed value of the GOF statistic (from the actual data set) is low (low significance), then the model should be rejected. Without such a statistic, one does not know if the best model is a good model, or simply the “least bad” of the family.

In this paper, we examine in some detail the issue of GOF when analysing anisotropy data on the cosmic microwave background. The vast majority of present analyses of the power spectrum data do not include proper GOF evaluations. The problem is particularly complicated by the fact that approximate likelihood methods must be employed in order to process the large volume of data and to explore a significant part of parameter space. These methods usually rely on power spectrum estimates, such as flat band–powers, extracted either from scan data, or from reconstructed sky maps. Because the power is quadratic in the temperature fluctuations, it is clear that these estimates are not Gaussian distributed. The traditional approach of $`\chi ^2`$ minimisation incorrectly assumes that power estimates are Gaussian distributed, something that can lead to a bias in determining the best model (e.g., Douspis et al. 2001a). For the same reason, the value of the reduced $`\chi ^2`$ at the best model does not retain its usual statistical meaning and may therefore not be simply used as a GOF statistic.

Approximations to the band–power likelihood function that permit more rigorous analyses have been proposed (Bond, Jaffe & Knox 2000; Bartlett et al. 2001). The question remains, however, of how to correctly evaluate the GOF of the best model. Such an evaluation requires knowledge of the distribution of the power estimates, which is not necessarily the same as the likelihood function. Using the same approach as Bartlett et al. (2001; hereafter paper 1), we propose an ansatz for the distribution of band–power estimates and test it against Monte Carlo simulations of certain MAX and Saskatoon data sets. The ansatz requires only minimal experimental information, and it appears to work well. We therefore use it to construct two GOF statistics, which we then apply to various ensembles of the present CMB data set.

## 2 Likelihood Method

It is useful to begin with a discussion of GOF in the context of a complete likelihood analysis. Although computationally challenging (in fact, impossible for large data sets: Bond et al. 2000; Borrill 1999ab), a likelihood approach is conceptually straightforward and our discussion serves to highlight certain important points. Such an analysis is in any case required for a small subset of data in order to test approximate methods (see, for example, Douspis et al. 2001a, hereafter paper 2).

Following the notation of papers 1 and 2, we write the likelihood function as (we consider only Gaussian perturbations)

$$(\stackrel{}{\mathrm{\Theta }})\mathrm{Prob}(\stackrel{}{d}|\stackrel{}{\mathrm{\Theta }})=\frac{1}{(2\pi )^{N_{\mathrm{pix}}/2}|C|^{1/2}}e^{\frac{1}{2}\stackrel{}{d}^tC^1\stackrel{}{d}}$$
(1)
where $`C(\stackrel{}{\mathrm{\Theta }})`$ is the correlation matrix (a function of the model parameters $`\stackrel{}{\mathrm{\Theta }}`$ and including a contribution from instrumental noise), and $`\stackrel{}{d}`$ is column vector listing the pixel values222These ‘pixels’ may either be the simple pixels of a map, or temperature differences, as given by, for example, MAX.. The elements of $`\stackrel{}{\mathrm{\Theta }}`$ may be either the cosmological parameters, or a set of band–powers. Maximising the likelihood function over the parameters defines the “best model” corresponding to the parameters $`\stackrel{}{\mathrm{\Theta }}_{best}`$.

In the present situation, we are greatly aided by the Gaussian form of Eq. (1) in the data vector, $`\stackrel{}{d}`$. Given the best model, the most obvious GOF statistic is then clearly

$$gof=\stackrel{}{d}^t\stackrel{~}{C}^1\stackrel{}{d}$$
(2)
where $`\stackrel{~}{C}C(\stackrel{}{\mathrm{\Theta }}_{best})`$ is the correlation matrix evaluated at the best model. For the Gaussian fluctuations we have assumed, this quantity follows a $`\chi ^2`$ distribution, with a number of degrees–of–freedom (DOF) approximately equal to the number of pixels minus the number of parameters333This recipe does not strictly apply in the present case, because the parameters are non–linear functions of the data; it is nevertheless standard practice. In any case, the number of pixels is in practice much larger than the number of parameters..

## 3 $`\chi ^2`$ method

For a variety of reasons (e.g., increased computational speed or inaccessible pixel data) most parameter estimations use power estimates, $`\delta T^2`$, as their starting point, such as those shown in Figure 1. A classic minimisation of $`\chi ^2`$

$$\chi ^2(\stackrel{}{\mathrm{\Theta }})=\underset{n=1}{\overset{N_{exp}}{}}\left(\frac{\delta T_n^{obs}\delta T_n(\stackrel{}{\mathrm{\Theta }})}{\sigma _n}\right)^2$$
(3)
is commonly used to find $`\stackrel{}{\mathrm{\Theta }}_{best}`$ and the best model, where $`\sigma _n=\sigma _+`$ ($`\sigma _{}`$) if the model passes above (below) the data point. The obvious GOF statistic would then be the value of the $`\chi ^2`$ evaluated at the minimum: $`gof=\chi ^2(\stackrel{}{\mathrm{\Theta }}_{best})`$. As already noted, this whole procedure is inappropriate because power estimates do not follow a Gaussian distribution. It is of course true that if the number of contributing effective degrees–of–freedom 444less than the number of original pixels by a factor depending on the pixel–pixel correlations; see paper 1 is large, a power estimate will closely follow a Gaussian; this, however, is never the case on the largest scales probed by a survey. We shall see in the following that, for actual CMB data, the $`\chi ^2`$ approach leads to quantitatively different results than other, more appropriate GOF statistics. For future reference, we show the value of this classic $`\chi ^2`$ in Table 1.

## 4 Proposed approximation

To improve the $`\chi ^2`$ analysis, several authors have proposed approximations to the band–power likelihood function $`(\delta T)`$ that may be constructed based on only minimal information about the experimental set–up (Bond, Jaffe & Knox 2000; Paper 1). One then arrives at the likelihood as a function of cosmological parameters $`\stackrel{}{\mathrm{\Theta }}`$ with $`[\delta T(\stackrel{}{\mathrm{\Theta }})]`$. Unfortunately, these approximate likelihood functions do not retain the normalisation of the full likelihood over pixels (Eq.1). This is a crucial point for GOF: we cannot deduce the quantity in Eq. (2) from the value of the approximate likelihood at its maximum.

An alternative way to build a GOF statistic would be from the expected distribution of power estimates, i.e., the distribution of points in Figure 4 around the model curve. Testing the observed dispersion of actual power points around the best–fit model against this expectation amounts to a GOF. The main difficulty in this approach is that we do not have an expression for the distribution of ML power estimates. It is important to understand that this distribution is not the same as the band–power likelihood, whose maximum is used to find the estimated power. In this section, we first motivate and then test an approximation to the distribution of ML power estimates.

### 4.1 Motivating an ansatz

Our approach will be the same as in Paper 1, and the following results thus apply when using the approximate band–power likelihood introduced therein. We motivated our likelihood approximation with an unrealistically simplified situation of $`N_{\mathrm{pix}}`$ uncorrelated pixels and uniform noise (refered to hereafter as the simple picture). This suggested a functional form depending on two parameters, an effective number of degrees–of–freedom $`\nu `$ and a noise parameter $`\beta `$; in the simple picture, $`\nu =N_{\mathrm{pix}}`$ and $`\beta ^2`$ is the noise variance. These two parameters could be found in realistic situations by adjusting to published flat–band confidence intervals (“errors”). The particular advantage of such a technique is that it permits an approximate likelihood analysis based on rather rudimentary information often found in the literature; this is an important advantage for many first generation experiments. In this same spirit, we now propose an ansatz for the ML band–power estimators.

For the simple picture ($`\nu =N_{\mathrm{pix}}`$ and $`\beta ^2=`$ noise variance), we showed in Paper 1 that the ML band–power estimator, $`\delta T^2`$, was a linear transform of a $`\chi _{N_{\mathrm{pix}}}^2`$ random variable:

$$\chi _\nu ^2=\nu \frac{([\delta T]^2+\beta ^2)}{([\delta T(\stackrel{}{\mathrm{\Theta }})]^2+\beta ^2)}$$
(4)
where $`\delta T^2(\stackrel{}{\mathrm{\Theta }})`$ is the band–power of the underlying model. In a realistic situation where $`\nu `$ and $`\beta `$ are found from published power estimates, there is no a priori guarantee that this formula applies with the same values of $`\nu `$ and $`\beta `$. One is, of course, tempted to suppose that the same values may in fact be used, at least approximatively. This hope forms the basis of our proposed ansatz for the band–power estimator distribution:

$`𝒫(\delta T^2|\stackrel{}{\mathrm{\Theta }})`$ $``$ $`Y^{(\nu /21)}e^{Y/2}`$ (5)
$`Y[\delta T^2]`$ $``$ $`\nu .{\displaystyle \frac{([\delta T]^2+\beta ^2)}{([\delta T(\stackrel{}{\mathrm{\Theta }})]^2+\beta ^2)}}`$
The underlying model band–power $`\delta T^2(\stackrel{}{\mathrm{\Theta }})`$ is in practice taken to be the ML estimate. The essential spirit of our approach is that, knowing the flat–band estimates and the 68 and 95% confidence levels, one is able to reconstruct the entire likelihood function and (now) the probability distribution of the estimate $`\delta T_{\mathrm{fb}}`$.

The only way to be sure that this proposed method actually works is by testing it against Monte Carlo simulations of some experiments before generalised it. We mention at least one reason for caution: the quantity $`\nu `$ represents an effective number of DOF, reduced from $`N_{\mathrm{pix}}`$ by inter–pixel correlations, applicable to the likelihood function; it is not at all clear that this same effective DOF applies equally well to the power estimator distribution (as it does in the simple picture). In particular, note that since the same data where used to find the best–fit model, we might expect a reduction in DOF, something familiar from the classic reduced $`\chi ^2`$ test. Here, however, we have no clear idea of the reduction. Fortunately, the proposed method nevertheless appears valid, as the following Monte Carlo simulations demonstrate it.

### 4.2 Testing the ansatz

We simulated many different data realizations of the MAX ID (Clapp et al. 1994) and Saskatoon (Netterfield et al. 1996) experiments in order to reconstruct the corresponding ML power estimator distribution. For example, we ran 30000 realizations of MAX ID at a frequency of $`3.5\mathrm{cm}^1`$ in the following manner: we first compute the flat band–power and the one dimensional likelihood function for the actual observational data. Knowledge of the latter provided the value of the pair $`(\nu ,\beta )`$. The maximum of the likelihood function gave us the “best model”, which was used to simulate pixels on the sky. In order to take into account all correlations, we simulated our pixels using the full pixel–pixel correlation matrix. We first computed the theoretical part of the correlation matrix evaluated for our “best model”. After diagonalization, we drew 30000 realisations of 21 pseudo–pixels from a Gaussian distribution centered on 0 and with the variances given by the eigenvalues. We reconstructed the “true” sky pixels using the transformation matrix (eigenvector matrix) and adding realizations of Gaussian noise (given by the known noise correlation matrix). We thus obtained 30000 sets of 21 pixels, correlated and drawn according to the best model ($`\delta T_{\mathrm{fb}}=57.3\mu K`$). For each realization, we derived the ML power estimate and build a histogram of its distribution.

Figure 2 shows the resulting distribution for $`MAXID3.5cm^1`$. Overplotted in red as the smooth curve is the ansatz Eq. (5) with the same values of $`(\nu ,\beta )`$ as found from the likelihood function. We see that the proposed approximate distribution is indeed a good representation of the true power estimator distribution.

The same kind of analysis was performed for the Saskatoon K–band 3–point data, an altogether different observing strategy. Once again, the approximation fitted the distribution to high accuracy. On the basis of these agreements, we will now adopt the proposed form in Eq. (5) as a good representation of the distribution of ML band–power estimators.

### 4.3 From probability function to GOF

On the basis of the distribution Eq. (5), we now construct two GOF statistics. The goal is to define a scalar quantity, $`gof`$, that measures the scatter of points around a given model and whose distribution is known under the hypothesis that this model represents the “truth” (the null hypothesis). An improbable value of $`gof`$ would indicate that there is a problem.

Both constructions assume that the band–powers are independent. This of course is not strictly true, but generally speaking published band–powers do not have strong statistical correlations; for example, the residual correlation between the Saskatoon bands is at a level of $`10\%`$. Calibrations errors, on the other hand, do induce important band–band correlations. As already mentioned, the present work does not include calibration errors, and any “bad fit” indicated by our GOF tests could indicate either a false model, or that calibration errors are important. Our aim here is to show the ability of a proper GOF to identify problems with CMB power data fits, and to demonstrate the advantage of the two proposed GOF statistics on the naive and inappropriate classic $`\chi ^2`$.

#### 4.3.1 Generalized $`\chi ^2`$

For each band–power $`i`$, consider the variables $`\alpha _i`$ defined as follows:

$$_{\mathrm{}}^{\alpha _i}\frac{1}{\sqrt{\pi }}\mathrm{exp}(x^2/2)𝑑x=p_i$$
(6)
where $`p_i_{\beta _i^2}^{\delta T_i^2}𝒫_i(\delta T^2|\stackrel{}{\mathrm{\Theta }})𝑑\delta T^2`$ is calculated using Eq. (5). The $`\alpha _i`$ are thus Gaussian random variables with zero mean and a variance of unity. Hence, the sum $`gof=_1^{N_{exp}}\alpha _i^2`$ follows a $`\chi ^2`$ distribution with $`N_{\mathrm{band}}`$ DOF and provides a handy GOF statistic.

#### 4.3.2 Characteristic functions

Another way to define a GOF statistic for a fit to $`N_{\mathrm{band}}`$ power points relies on the following property of characteristic functions: the characteristic function for the sum of independent random variables is given by the product of the individual characteristic functions. Given the $`N_{\mathrm{band}}`$ random variables $`Y_i`$ and their probability distributions $`𝒫_i`$ (Eq. 5), we calculate the distribution of the random variable $`z_iY_i`$, which will represent the goodness of fit, as follows: for each $`Y_i`$ we can compute the corresponding characteristic function $`\mathrm{\Phi }_i(k)`$. Then using the property cited above, we can construct the characteristic function $`\mathrm{\Phi }_z(k)`$ of the variable $`z`$ by $`\mathrm{\Phi }_z(k)=\mathrm{\Phi }_1(k)\mathrm{}\mathrm{\Phi }_{N_{\mathrm{band}}}(k)`$. The probability distribution function of $`z`$, $`(z)`$, is then given by the inverse Fourier transform of $`\mathrm{\Phi }_z(k)`$. This approach is particularly straightforward in our case because the probability function given in Eq. (5) is just a $`\chi ^2`$ law with $`\nu _i`$ DOF, whose characteristic function $`\mathrm{\Phi }_i`$ has an analytic form. Multiplication of the individual characteristic functions thus gives an analytical expression whose inverse Fourier transform is itself a $`\chi ^2`$ distribution in $`z`$, with $`\nu =\nu _i`$ DOF:

$`(z)=z^{\nu /2}e^{z/2}`$ (7)
$`\mathrm{with}z={\displaystyle \underset{i}{}}Y_i\mathrm{and}\nu ={\displaystyle \nu _i}`$
The variable $`gof=z`$ is thus (another) $`\chi ^2`$–distributed quantity that provides a useful GOF statistic.

## 5 “The good, the bad and the GOF” or Are CMB fluctuations consistent with a Gaussian distribution?

### 5.1 Application

In this section we apply each of the above GOF statistics to the CMB data set shown in Figure 1; note that this does note include the most recent BOOMERanG, MAXIMA and DASI results. Adding these new data will essentially results in reducing the “$`\chi ^2`$” distributed $`gof`$ values without changing drastically the results presented in this section. Our overall approach is as described in Le Dour et al. 2000 (hereafter paper 3) and Douspis et al. 2001a, where we used the likelihood approximation given in paper 1 to find the best model. We consider three combinations of data: Data set 1 contains all points (ALL)555Actually, we noticed that our approximation fails to recover the MC simulations for upper limits. For this reason we do not include them in our analysis; set 2 consists of all the data minus the Python 5 results (Coble et al. 1999) (ALL-5); and set 3 combines just COBE (Tegmark & Hamilton 1997), MAXIMA (Hanany et al. 2000) and BOOMERanG (de Bernardis et al. 2000) (CMB).

The best model for each data set will be referred to as $`\mathrm{BM}_{\mathrm{ALL}}`$, $`\mathrm{BM}_{\mathrm{ALL}5}`$,$`\mathrm{BM}_{\mathrm{CMB}}`$. A summary of the various GOF statistics for these models is given in Table 1; the lines are labelled by GC for “generalised $`\chi ^2`$”, CF for “characteristic functions”, and “$`\chi ^2`$” for the classic $`\chi ^2`$ of Eq. (3)666We noticed that $`\sigma _i`$ is given different definitions in the literature. When considering the evaluation of the GOF using each definition, we found that the value of the GOF is quite sensitive to the definition of $`\sigma _i`$. We consider in this paper the technique giving the best value of the GOF.

For GC technique, the $`gof`$ is directly equivalent to the absolute value of a $`\chi ^2`$. To convert this into percentage, we need to know the DOF. The latter is given in our case by the number of experiments taken into account in each set minus the number of free cosmological parameters.

For the CF technique, the percentage given in table 1 is obtained by integrating the probability distribution function $`f`$ of $`Z_Y`$ from infinity to $`Z_{obs}=_i^{N_{set}}Z_i^2`$ where $`N_{set}`$ is the number of experiments in each set.

The figure 4 summarises the results given by our CF test on both data sets 1 and 2. The line gives the function to integrate and the shaded part is the integrated part corresponding to the numbers given in Table 1. The solid (blue) line and shaded part correspond to data set 1, and the dashed (red) line and arrow to data set 2.

### 5.2 Discussion

The first remark to be made based on Table 1 is that the complete data set (set 1) is inconsistent with a Gaussian sky fluctuations, according to all three techniques; the GC method, for example, excludes this hypothesis at more than 99.99 %. This means in particular that it is not appropriate to search cosmological constraints, because the whole class of models considered is ruled out. This could be due to several effects, in particular the fact that we do not include calibration uncertainties in our analysis.

The situation is different if we remove Python 5 (set 2) from the analysis. In this case, our two evaluations of the GOF (GC and CF) both accept the hypothesis of Gaussian sky fluctuations. In contrast, the classic (but inappropriate) $`\chi ^2`$ statistic marginally excludes such hypothesis. Figure 3 illustrates the difference between our GC method and the classic $`\chi ^2`$, data point by data point (for a subset of data set 1). Triangles show individual $`\chi ^2`$ values, while boxes correspond to the $`\nu _i`$ defined in section 3. We see that the classic $`\chi ^2`$ overpenalizes the fit for outliers, a conclusion already noted in paper 2.

Finally, we can see that all three methods accept the Gaussian hypothesis as a good representation to the COBE, MAXIMA and BOOMERanG data (set 3).

## 6 Conclusion

We have discussed three different ways of estimating the GOF to CMB band–powers. A GOF statistic is a key element of any parameter estimation study, and a good fit must be insured before considering parameter constraints. The classic $`\chi ^2`$ GOF statistic is not rigorously applicable to power spectrum data, because power estimates are not Gaussian distributed quantities. We propose instead two alternative GOF statistics based on an approximation to the distribution of power estimators. This approximation was motivated by the same kind of arguments presented in paper 1 for the likelihood function. The distribution of a power estimator is a different quantity than the likelihood function used to define the estimator. We tested the approximation presented here against Monte Carlos simulations of CMB observations and found that it reproduced well the distribution of the maximum likelihood band–power estimator.

We then constructed two different GOF statistics, whose distributions were found using the approximate power estimator distribution. With the same, rather minimal information required to build the likelihood approximation (paper 1), we are now also able to develop a GOF statistic to test the quality of the maximum likelihood model to a set of band–power data, thereby allowing a complete statistical analysis of anisotropy data from diverse observations. The method is limited by the fact that we are unable to account for correlations between band–powers; this, however, is not a serious restriction, as these correlations are usually rather unimportant for the final results based on current data sets.

In applying this approach to a set of band–power data of Figure 1 we found that the “best model” obtained is in fact a bad fit. In other words, the data are unlikely to have been drawn from a Gaussian distribution represented by such a model. The fit becomes acceptable if we exclude the Python 5 points from the analysis, according to our GOF statistics. This is most likely due to the fact that we do not account for calibration errors, and so the bad fit probably just indicates that the adopted calibration is incorrect. It is interesting to note that, even with Python 5 removed, the classic $`\chi ^2`$ still marginally rejects the best fit. We traced this behaviour to the fact that this method over weights the importance of “outliers”.

The important cosmological conclusion is that this CMB data set (excluding Python 5, due to our inability to account for calibration errors) is consistent with Gaussian sky fluctuations drawn from the best–fit inflationary model.

A final remark concerns the possibility offered by the development of an approximated distribution function of the estimators. In the application of current Monte Carlo methods for $`C_{\mathrm{}}`$’s extraction (eg. Szapudi et al. 2000, MASTER: Hivon et al. 2001), the estimator distribution is a natural output. The likelihood function needed in parameter estimations is however unknown. The present study suggests that we could reconstruct the likelihood function directly from the estimator distribution. The two parameters ($`\nu `$ and $`\beta `$) can be fitted on the estimator distribution and then used in the approximated likelihood function of Bartlett et al. 2001. Consequently one is then able to reconstruct all the likelihood function and to perform a proper parameter estimation.

###### Acknowledgements.

M. D. would like to thank Nabila Aghanim for usefull comments and corrections.


text:
# Active galactic nuclei and the minor merger hypothesis

## 1 Introduction

In this paper we investigate the hypothesis that activity in galactic nuclei is induced by the acquisition and disruption of a small gas-rich companion galaxy. There has been considerable discussion of this ‘minor merger’ hypothesis in the literature (see, for example, recent discussions in Taniguchi, 1999; Chatzichristou, 2000,a,b; 2001 a,b, and papers therein).

In elliptical galaxies, a study of the effects of minor mergers is motivated by the observational evidence (van Dokkum & Franx, 1995) that dust discs are present in the cores of a high fraction of ellipticals, with a higher detection rate in radio galaxies (Verdoes Kleijn et al., 2000). This result is consistent with the results of surveys of 3CR radio galaxies by de Koff et al. (2000) and by Martel et al. (1999). The simplest interpretation of these is that the dust structure seen in these galaxies is the debris trail of a small merger incident. The analysis of de Koff et al. (2000), and the finding by Sparks et al. (2000) that galaxies with optical jets tend to display face-on (i.e. round) dust discs, both suggest that the radio jets in radio galaxies show a tendency to be aligned perpendicular to the observed dust disc structure. In this picture a recent merger would have stimulated the current nuclear activity. The surrounding dust structure would indicate the plane of the orbit of the disrupted satellite, and the resulting radio jets would naturally be expected to be perpendicular to this plane.

In view of this, we investigate first the conditions required for a small gas-rich satellite galaxy to be accreted by a larger galaxy in such a way that it gives rise the a dust disc similar to those observed. To focus the discussion we concentrate on modelling the formation of the dust disc in the otherwise unremarkable E1 galaxy NGC 3379. Our results are straightforwardly applicable to other gas-poor ellipticals. The dust disc in NGC 3379 (van Dokkum & Franx, 1995; Ferrari et al., 1999) has a radius $`100\mathrm{pc}`$, and a mass estimated to be around $`150\mathrm{M}_{}`$. In addition Ferrari et al. (1999) show that there is more patchy dust extending out to about $`1\mathrm{kpc}`$ from the nucleus. We adopt a simplified procedure to model the interaction. We use spherical models for both the major galaxy and the minor mergee and model the motion of the orbit of the satellite using simple dynamical friction, along with a straightforward model for the tidal disruption of the satellite (Section 2). We use a simple model for the subsequent evolution of the stripped gas (Section 3). Our basic conclusion (Sections 5.1 and 6) is that in order to generate a dust disc of the kind observed it is necessary for the initial trajectory of the merging satellite to be accurately directed towards the nucleus of the major galaxy.

However, reality turns out to be more complicated than this simple picture might suggest. First there is the recent analysis by Schmitt et al. (2002) who consider only those E/S0 galaxies with active nuclei which have $`\mathrm{𝑤𝑒𝑙𝑙}\mathrm{𝑑𝑒𝑓𝑖𝑛𝑒𝑑}`$ dust discs in their central regions, rather than galaxies whose dust distribution is more dispersed. Schmitt et al. conclude that the radio jets from the centres of these discs are not in fact oriented perpendicular to the dust discs’ axes, and the angles the jets make with the axes are consistent with being distributed evenly in the range $`0\mathrm{°}75\mathrm{°}`$. Second there is the finding by Kinney et al (2000) that in Seyfert galaxies the directions of the central jets are consistent with being oriented randomly in space, independent of the plane of the gas disc of the host spiral galaxy. If the minor merger hypothesis is to hold, then the apparent misalignments between between the angular momentum of the accreting material (as evidenced by the observed discs in early type galaxies and by the disc plane in Seyferts) and the observed direction of the jets needs some explanation. One possible explanation is the idea that the jet direction is governed by the spin of the central black hole (through the Bardeen-Petterson effect (e.g. Scheuer and Feiler, 1996 and references therein) and that for some reason the spin of the black hole is misaligned with the current accretion flow. While this is a reasonable explanation if the galaxy is spherically symmetric it would seem difficult to achieve if the spin of the black hole has been generated by the gas accretion process and if the potential of the galaxy is such the the accreted gas would (in the long term) be expected to have a preferred direction for its angular momentum vector. One way round this problem, following the ideas of Wilson & Colbert (1995; see also Merritt, 2002), is if the black hole spin is generated not by accretion of gas, but rather by accretion of black holes presumably contained in the nuclei of the small merging satellites. Even a small black hole (say a tenth of the mass of that in the host galaxy) can in its final approach add substantial orbital angular momentum to the spin of the final coalesced black hole. For this process to work, what is required is that the central (nuclear) black hole of a merging satellite be able to reach the nuclear regions of the major galaxy with the angle between its orbital angular momentum and the axis of symmetry of the gravitational potential of the major galaxy distributed over a wide range of angles.

To investigate this possibility we extend our analysis to the case where the major host galaxy is a spiral galaxy. Here we are not concerned with the fate of any gas that might have been in the merging satellite, but are, rather, concerned with the fate of its nuclear black hole. We make use of similar simplified dynamics using the concepts of gravitational drag and a similar approximation to take account of tidal stripping to those we used before, but we now apply these concepts to a a galaxy with a more complicated internal structure (Section 4). As before we find that the cross-section for an incoming satellite to reach the nucleus of the host galaxy in a Hubble time is relatively low. However here we pose the additional question: of those satellite nuclei which do reach the central regions of the host galaxy in a reasonable time, what is the distribution of the inclinations of their orbits relative to normal to the disc of the galaxy? In order for the merger hypothesis to remain viable we shall require the distribution we find to very broad. We summarise our findings and give discussion of our results in Sections 5.2 and 6 respectively.

## 2 Merging with a spherical galaxy

In this Section we describe our simplified galaxy dynamics and merging procedure.

### 2.1 Galactic models

As mentioned in Section 1, we first investigate the particular case in which the major galaxy is spherically symmetric as a reasonable approximate model for the merging of a satellite with the galaxy NGC 3379.

#### 2.1.1 Major galaxy (NGC 3379)

To model the density profile of NGC 3379 we have made use of Lima Neto, Gerbal and Márquez’s \[Lima Neto et al. 1999\] analytic approximation to the deprojection of its surface brightness profile obtained with HST \[Pastoriza et al. 2000\]111Note that there is a typographical error in Pastoriza et al. (2000): in equation 2 ‘exp’ should be ‘10 to the power of’ (Ferrari, private communication). Lima Neto et al. \[Lima Neto et al. 1999\] assume that the galaxy is spherical, with a constant mass-to-light ratio $`\mathrm{{\rm Y}}`$, and find that

$$\rho (r)=\rho _0\left(\frac{r}{a}\right)^p\mathrm{exp}\left[\left(\frac{r}{a}\right)^\nu \right],$$
(1)
where $`\nu =0.424`$ and $`p=0.572`$. Using the normalization to the surface brightness given by Pastoriza et al. \[Pastoriza et al. 2000\], together with an assumed distance to NGC 3379 of $`9.9\mathrm{Mpc}`$ and a mass to light ratio of $`\mathrm{{\rm Y}}=5.3\mathrm{{\rm Y}}_{}`$ (Magorrian et al. 1998), gives $`\rho _\mathrm{o}=3.96\times 10^3\mathrm{M}_{}/\mathrm{pc}^3`$, and $`a=18.6\mathrm{pc}`$. The resulting visible mass is $`2.9\times 10^{10}\mathrm{M}_{}`$, with a half-mass radius of about $`800\mathrm{pc}`$, excluding any dark halo. We find that the inclusion of a dark halo or a central $`10^8\mathrm{M}_{}`$ black hole (Magorrian et al. 1998; Gebhardt et al. 2000) has negligible effect on our results below.

We calculate the mass distribution, potential and velocity dispersion of the galaxy through the usual equations

$`M(r)`$ $`=`$ $`{\displaystyle _0^r}\rho (r^{})4\pi r^2\text{d}r^{}`$ (2)
$`\mathrm{\Phi }(r)`$ $`=`$ $`{\displaystyle _r^{\mathrm{}}}{\displaystyle \frac{GM(r^{})}{r^2}}\text{d}r^{}`$ (3)
$`\sigma ^2(r)`$ $`=`$ $`{\displaystyle \frac{1}{\rho (r)}}{\displaystyle _r^{\mathrm{}}}{\displaystyle \frac{GM(r^{})\rho (r^{})}{r^2}}\text{d}r^{},`$ (4)
where we have assumed for simplicity that the galaxy’s velocity dispersion tensor is isotropic. We evaluate the integrals adaptively using Simpson’s rule until an accuracy of one part in $`10^6`$ is reached and store the calculated values of $`M(r)`$, $`\mathrm{\Phi }(r)`$ and $`\sigma (r)`$ on grids with 30 points spaced logarithmically in radius between $`100\mathrm{pc}`$ and $`100\mathrm{kpc}`$. Values between tabulated points are obtained via cubic spline interpolation.

#### 2.1.2 Satellite

We have considered two different models for the small infalling satellite galaxy. One is a rigid Plummer sphere with density distribution given by

$$\rho (r)=\frac{3M}{4\pi b^3}\left(1+\frac{r^2}{b^2}\right)^{5/2},$$
(5)
and the other is a rigid Hernquist sphere with density distribution given by

$$\rho (r)=\frac{M}{2\pi }\frac{a}{r}\frac{1}{\left(a+r\right)^3}.$$
(6)
For both of these cases we calculate the mass and velocity dispersion as for the host galaxy, except with a grid spaced logarithmically between $`0.1\mathrm{kpc}`$ and $`10\mathrm{kpc}`$. The total mass ($`M=10^9\mathrm{M}_{}`$) and scale lengths ($`b=0.4\mathrm{kpc}`$ and $`a=0.23\mathrm{kpc}`$) were chosen to match the estimated mass and half-mass radius ($`0.55\mathrm{kpc}`$) of Sagittarius \[Helmi & White 2000\], as a typical dwarf galaxy. We will see in Section 2.3 that the choice of model for the satellite does not make a significant difference to our results.

### 2.2 The orbit of the infalling satellite

In computing the orbit of the infalling satellite, we neglect changes which might be induced by the satellite in the host galaxy. This approximation has been the subject of debate in the literature \[Zaritsky & White 1988, Hernquist & Weinberg 1989\], with the conclusion being reached that it is reasonable to ignore the galaxy response for low mass satellites \[Velázquez & White 1999\]. Physically, this is what is expected, as the response induced in the galaxy is proportional to the mass ratio of the satellite to the galaxy $`M_s/M_g`$, and the effect of this response is again proportional to the mass ratio, so the overall effect will be proportional to $`\left(M_s/M_g\right)^2`$.

We model the motion of the satellite using Chandrasekhar’s dynamical friction formula \[Binney & Tremaine 1987\]. Thus

$$\frac{\text{d}𝒗_𝒔}{\text{d}t}=\frac{kM_s\rho }{\left|𝒗_𝒔\right|^3}𝒗_𝒔,$$
(7)
where $`𝒗_𝒔`$ is the velocity of the satellite, $`\rho `$ is the density of the medium, $`k`$ is given by

$$k=4\pi \mathrm{ln}\mathrm{\Lambda }G^2\left(\text{erf}(X)\frac{2X}{\sqrt{\pi }}\mathrm{exp}\left(X^2\right)\right),$$
(8)
$`G`$ is Newton’s gravitational constant, $`X`$, the dimensionless speed of the satellite, is given by

$$X=\frac{\left|𝒗_𝒔\right|}{\sqrt{2}\sigma (r)},$$
(9)
and $`\mathrm{\Lambda }`$ is the Coulomb logarithm discussed below.

Equation (7) applies to the case of a point mass moving through an infinite homogeneous medium with a Maxwellian velocity distribution. However, numerical experiments \[Velázquez & White 1999\] show that it also provides a reasonable approximation to the gravitational drag experienced in more realistic situations, provided one chooses the Columb logarithm $`\mathrm{ln}\mathrm{\Lambda }`$ appropriately. To determine whether it is adequate for the situation we are considering, we use it to reproduce the results from the infall of a satellite into a dark halo when the galaxy was modelled using a full $`N`$-body treatment \[van den Bosch et al. 1999\]. We note that this procedure has been investigated in some detail by Taylor & Babul (2001), who demonstrate that it does provide a reasonable approximation.

To perform the integration, we used the algorithm from the fifth-order Cash-Karp Runge-Kutta integrator ‘odeint\[Press et al. 1992\]. The timestep used is automatically adjusted by the routine to control errors, and in no case did the timestep reach a significant fraction of the dynamical time of the satellite.

As a test of our procedure, we use our method to model van den Bosch et al.’s \[van den Bosch et al. 1999\] Plummer sphere sinking into a truncated isothermal halo. We find that we are able to reproduce their results over a wide range of conditions if we choose $`\mathrm{ln}\mathrm{\Lambda }2`$. In Figure 1 we show a comparison between the full $`N`$-body calculations of van den Bosch et al (1999) and our approximate procedure for a typical case. As can be seen, there is reasonable agreement between the two. Using these results we shall adopt the same choice $`\mathrm{ln}\mathrm{\Lambda }=2`$ for the orbital computations contained in this paper. We should note that, for the spherical case, our treatment of dynamical friction also ignores the effects of possible galaxy rotation and velocity anisotropy, but unless these are extreme we expect that this approximation will be adequate in the current case.

### 2.3 Tidal stripping

In the model used by van den Bosch et al. \[van den Bosch et al. 1999\], the satellite maintained a constant mass all the way through the simulation. Due to tidal effects, the outer regions of the satellite will not remain bound to the satellite. We are interested, therefore, in the manner in which the infalling satellite is tidally stripped as it passes through the host galaxy. This has an effect upon the orbital dynamics, because the deceleration due to Chandrasekhar friction is proportional to the mass of the satellite.

We use a simple algorithm to simulate the tidal stripping of the satellite. At each instant, we define the tidal radius, $`r_t`$, of the satellite as being the radius at which the mean density of the satellite is equal to the mean density of the galaxy, viz.,

$$\frac{M_s(r_t)}{\frac{4}{3}\pi r_t^3}=\frac{M_g(r_s)}{\frac{4}{3}\pi r_s^3},$$
(10)
where $`M_s(r)`$ is the mass of the satellite contained within a radius $`r`$, $`M_g(r)`$ is the same for the galaxy and $`r_s`$ is the galactocentric radius of the centre of the satellite. We then truncate the satellite at $`r=r_t`$, keeping the density structure internal to this point unchanged. This reduces the mass of the satellite to $`M_s(r_t)`$. In reality the satellite will react to this density truncation on an internal dynamical timescale, but since by construction222The tidal radius is defined by equating the mean density of the satellite to the mean density of the main galaxy, and dynamical timescales just depend on mean densities. this timescale is similar to the orbital timescale at this point, neglect of this effect is likely to be a reasonable approximation for current purposes.

In Figure 2, we show the outcome of this stripping algorithm by plotting the fraction of the satellite galaxy which remains bound as a function of its distance from the centre of the host galaxy. We show the results for both models of the satellite galaxy. The Plummer sphere is completely disrupted by the time it comes within 1 kpc of the centre of the host galaxy. The Hernquist sphere loses 99% of its initial mass by the time it comes within is $`400\mathrm{pc}`$ from the centre of the galaxy. These distances are much greater than the size of the observed dust disc, which has a radius of around $`100\mathrm{pc}`$, but are not that dissimilar to the overall size of the dust spread within the galaxy. We note that these result depends only upon the density profiles of the satellite and the galaxy. They are just the result of solving Equation (10), and do not involve any consideration of the dynamics.

## 3 Deposition and location of the stripped gas

### 3.1 Low energy orbits

To make a preliminary investigation of the initial conditions required for a merging satellite to be able to deposit significant gas into the nucleus of the host, we start by considering a series of low-energy orbits for the satellite galaxy. In each one the satellite starts on an orbit at an apocentre of $`10\mathrm{kpc}`$ from the galactic centre (that is, it starts with zero radial velocity at that point), and we vary the initial varying angular momentum. We characterize the initial angular momentum in terms of an eccentricity $`e`$, where $`e=(r_+r_{})/(r_++r_{})`$ and $`r_+(r_{})`$ are the apo-(peri-)centres of the orbit of a test particle with the same initial location and velocities. We use the models of NGC 3779 and the satellite as described in Section 2.

To simplify the discussion we shall assume that the gas in the satellite galaxy is distributed spatially in the same way as the stellar matter. We make this assumption for two reasons: First, such an assumption probably overestimates the concentration of the gas within the dwarf. This will lead to the gas being stripped stripped from the dwarf later, and thus at a lower angular momentum than would otherwise occur. Hence this assumption will be an upper bound on the magnitude of any effect seen. As we show later, even with this assumption, we cannot replicate the observed effects. Second, by assuming a spherical dwarf, we can dramatically reduce the computational power required for our simulations and thus investigate a larger region of orbital parameter space.

As an initial estimate of where the stripped gas ends up, the simplest assumption is that stripped material conserves its angular momentum about the centre of the host galaxy (since the host is spherically symmetric), and that each element of gas, by losing energy through dissipative processes ends up in a circular orbit corresponding to its initial angular momentum $`L`$. The radius $`r_c`$ of this orbit is given by

$$L^2=GM_g(r_c)r_c.$$
(11)
The results of making this simple approximation are shown in Figure 3. The outermost spike in each of these mass distributions is caused by the stripping which occurs at the first pericentre passage. For a high angular momentum orbit ($`e=0.5`$), the stripping occurs over a range of radii as the satellite’s orbit gradually decreases in size due to gravitational drag. For the intermediate angular momentum orbit ($`e=0.75`$) the stripping occurs more rapidly and thus the range of angular momentum in the stripped gas is reduced.

In order to acheive the required radial penetration, we also investigated a low angular momentum orbit ($`e=0.98`$) for which the satellite’s orbit is almost radial. In this case, the gravitational drag plays almost no role, and the satellite is totally disrupted on the first periastron passage, but the gas still ends up at radii exceeding $`250\mathrm{pc}`$. Thus we conclude that, using these approximations, the gas cannot end up within around $`100\mathrm{pc}`$ unless the satellite was initially on a very eccentric, indeed essentially radial, orbit.

### 3.2 Gas dynamics

The assumption we made above, that an element of stripped gas, as it falls in, loses energy by interactions with other elements, but does not exchange angular momentum with them is not fully realistic. Since the angular momentum is crucial in determining the ultimate location of the gas, and since the discs we are obtaining are still much larger than observed discs we now consider a more realistic model for the radial evolution of the stripped gas.

In the spirit of the approximations made in this paper, rather than undertaking a full hydrodynamical treatment (for example using SPH), we opt for computational simplicity and have followed the approach of Lin and Pringle \[Lin & Pringle 1976\]. In such a scheme, the gas particles move freely in the combined potential of the host galaxy and satellite, except that a form of dissipation is introduced which simulates thermally radiative shocks, by preventing the crossing of particle orbits while exactly conserving mass and angular momentum. In this scheme, if two particles are close to each other, they interact, losing a fraction of their energy, but retaining their total momentum. One effect of this is to redistribute the angular momentum between different particles. For example, a circular ring of particles acts as an accretion disc and gradually spreads out to form a wider ring. In this paper, we consider inelastic collisions between particles in order to produce the greatest dissipation. This is clearly the most extreme situation, and therefore the results of this section will provide an upper bound on the magnitude of this effect.

The scheme used to determine when collisions occur is simple: we divide the space in the host galaxy into cells, and, after every time step of the integration, if two (or more) particles occupy the same cell, they interact. We use a cylindrical polar grid $`(R,\varphi ,z)`$ to describe the cells, with the $`z=0`$ plane aligned with the orbit of the satellite galaxy. The cell boundaries are spaced logarithmically in $`R`$, from $`100\mathrm{pc}`$ to $`20\mathrm{kpc}`$ with $`R_{n+1}/R_n=1.02`$, and linearly in the $`z`$ from $`z=5.05\mathrm{kpc}`$ to $`z=5.05\mathrm{kpc}`$ with $`z_{n+1}z_n=0.1\mathrm{kpc}`$. The choice of $`R_{n+1}/R_n`$ is governed by a balance between the need for computational efficiency (if the zones are too small then there are two many cells with too few particles in them which permits interpenetration of fluid flows), and a need to keep the zones small, so that viscous effects do not become too rapid, resulting in unphysically rapid spreading of the gas discs which form. We performed runs with $`R_{n+1}/R_n`$ of $`1.05`$, $`1.02`$ and $`1.01`$. The results from the runs at $`1.02`$ and $`1.01`$ were very similar, implying that convergence had been achieved. Hence we choose $`R_{n+1}/R_n=1.02`$ to reduce the computational time required for our simulations. We choose 314 equispaced points in $`\varphi `$. This means that the projection of the cells onto the $`z=0`$ plane is approximately square and that the effective cross section of our gas particles varies from $`2\mathrm{pc}`$ at $`R=100\mathrm{pc}`$ to $`400\mathrm{pc}`$ at $`R=20\mathrm{kpc}`$.

We create one gas particle, of mass $`10^5\mathrm{M}_{}`$ for every $`10^5\mathrm{M}_{}`$ stripped from the satellite. The particle is started with a velocity

$$𝒗_𝒈=𝐯_s+\sigma _s\left(\left|𝒓_𝒈𝒓_𝒔\right|\right)\widehat{𝒓},$$
(12)
where $`𝒗_𝒈`$ and $`𝒗_𝒔`$ are the velocity of the gas particle and the satellite respectively, $`𝒓_𝒈`$ and $`𝒓_𝒔`$ are the positions of the centres of the gas particle and the satellite, $`\widehat{𝒓}`$ is a unit vector in a random direction and $`\sigma _s(r)`$ is the velocity dispersion of the satellite at a radius $`r`$. A particle is removed from the computation when one of the following conditions is met: the particle is within $`100\mathrm{pc}`$ of the centre of the galaxy (at which point, the gas is considered to have become part of the nuclear gas disc); the particle is within 90% of the current tidal radius of the satellite (when it is considered to have been reabsorbed by the satellite); or, the particle moves beyond the outermost edges of the grid at $`R=20\mathrm{kpc}`$ and $`\left|z\right|=5.05\mathrm{kpc}`$.

The radii at which the gas settles initially using this model are shown in Fig. 4. For comparison with the previous estimates (Figure 3) we show the result for initial eccentricities of $`0.5`$ and $`0.75`$. We do not show the results for an initial eccentricity of $`0.98`$, because the majority of the stripped mass passes within $`100\mathrm{pc}`$ of the galaxy centre and is therefore removed from the simulation. By comparing these profiles with those in Figure 3, it can be seen that the result of allowing the gas particles to interact with each other in a more realistic manner is to broaden the distributions, and to move the innermost gas particles inwards by about a factor of two in radius. Some of this difference is due to the more realistic treatment of the stripping process, some is due to treating the gas dynamics more realistically and some is due to the artifical viscous effects introduced by our approximate scheme for treating the gas dynamics. We believe that the viscous effects are not significant here. The mean radius of the final gas distribution is within 20% of the results obtained using the simple prescription of the previous section, which implies the total energy loss is similar. This approximate similarity between the two approaches serves to indicate that any viscous effects introduced due by our grid-based dissipation scheme are reasonably small.

### 3.3 High energy orbits

The above results make it clear that in order to form a dust disc of the size of a few hundred parsecs of the kind which is seen in NGC 3379, the orbit of the hypothesized incoming satellite galaxy is not of the form of a low energy, low eccentricity orbit which goes through multiple pericentre passages before the satellite is disrupted. Rather the initial orbit of the incoming satellite needs to have a sufficiently high eccentricity (i.e. low angular momentum) that the satellite is likely to be completely disrupted on its first pass.

For a more detailed investigation of the high eccentricity regime, it is more useful to parameterize the problem in terms of an impact parameter $`b`$ and an initial velocity $`V_0`$. In addition, to model the approach of the satellite galaxy to the host’s core, we now start the satellite’s orbit at a galactocentric radius of $`80\mathrm{kpc}`$. We use the same procedure as outlined above (Section 3.2) to treat the gas stripping and subsequent gas dynamics, but, since the pericentric distances are now considerably reduced, we now decrease the inner radius at which particles are removed from the computation from $`100`$ to $`10\mathrm{pc}`$. This is necessary because of the change in the nature of the orbits of the gas particles. When we were considering low energy orbits, the particles which reached $`100\mathrm{pc}`$ were on essentially circular orbits and are assumed to become part of the nuclear gas disc, but in the high eccentricity case, particles reach $`100\mathrm{pc}`$ with velocities much greater than circular and as such do not necessarily get incorporated into the disc.

With these initial conditions, we find that a disc of similar size to that observed in NGC 3379 can be reproduced if the specific angular momentum of the incoming satellite, $`bV_0`$ is around $`50\text{kpc km s}^1`$. The surface density profile of the disc obtained using this value for the orbital angular momentum is shown in Figure 5. In this situation, around 1% of the gas particles fall within the radius $`10\mathrm{pc}`$ and so are removed from the calculation.

We conclude that in order to obtain a dust disc the approximate size of the one observed it is necessary for the satellite galaxy to be on an orbit which is almost precisely oriented towards the galactic nucleus. In order to estimate of the likelihood of such an occurrence, it is of course necessary to have some knowledge of the space and velocity distribution of the prospective mergee population. This we do not have. Numerical computations of galaxy formation through mergers should in principle be able to give some insight into this (e.g. Steinmetz 1999), but current simulations do not currently have high enough mass resolution for our purposes.

We may, however, use the above results to estimate the fraction of those merging satellites disrupted on their first pericentre passage which give rise to a dust disc of the kind observed in NGC 3379. If we assume that the satellite galaxy, as it falls into the main galaxy, has a velocity approximately equal to that of the dispersion velocity in the main galaxy, which is about 200 km s-1, then for this satellite give rise to the observed gas disc, we have seen that we require $`bV_050\text{kpc km s}^1`$, and thus that the pericentre of the orbit to be $`b250\mathrm{pc}`$. For the satellite to be disrupted on the first pass through the main galaxy, we have seen (Section 2.3) that we require the pericentre distance to be less than around $`b1\mathrm{kpc}`$. If we assume that incoming satellites have their impact parameters evenly distributed over this range then this would imply that only of the order of $`4\%`$ of those satellites disrupted on the first pass are on appropriate orbits to form a disc such as observed in NGC 3379. Of course, the majority of satellites have larger impact parameters. They would be captured by the main galaxy, and would then slowly work their way inwards in the manner described in Section 2.3, distributing their gas and dust widely through the host galaxy. This should lead to observable effects (Section 5). This in any case we expect that for every merging satellite which could produce the observed disc, there are many more which could not.

## 4 Merging with a spiral galaxy

For the reasons discussed in the Introduction, we now extend the foregoing analysis to consider the merging of a small satellite galaxy with a spiral host galaxy. Our interest in this case centres not so much on where the stripped gas may end up (since it likely in any case to collide with and merge with the gas already in the gas disc of the spiral), but rather on the fate of the nucleus of the satellite, or more especially, on the fate of the black hole presumably contained within it.

### 4.1 Galactic models

#### 4.1.1 The spiral

We choose to use a simple model of a typical disc galaxy, in which the galaxy is modelled using four axisymmetric components (representing the central black hole/star cluster, the bulge, the disc and the dark halo), each of which is represented by a Miyamoto-Nagai potential (Binney and Tremaine 1987, Equation 2-50a), which takes the form

$$\mathrm{\Phi }(R,z)=\frac{GM}{\sqrt{R^2+\left(a+\sqrt{z^2+b^2}\right)^2}},$$
(13)
where we use a cylindrical coordinate system $`(R,\varphi ,z)`$.

The density due to this potential is

$$\rho (R,z)=\frac{b^2M}{4\pi }\frac{aR^2+\left(a+3\sqrt{b^2+z^2}\right)\left(a+\sqrt{b^2+z^2}\right)^2}{\left(R^2+\left(a+\sqrt{b^2+z^2}\right)^2\right)^{5/2}\left(b^2+z^2\right)^{3/2}}.$$
(14)
We note that this reduces to a Plummer sphere in the case of $`a=0`$. From these expressions, and using the assumption that the internal galactic dynamics consists of an isotropic velocity dispersion coupled with a mean azimuthal streaming motion, we can calculate the velocity dispersion and mean streaming velocity of each component, via (Binney and Tremaine, 1987, Equations 4-65 and 4-66):

$`\sigma ^2(R,z)`$ $`=`$ $`{\displaystyle \frac{1}{\rho (R,z)}}{\displaystyle _z^{\mathrm{}}}\rho (R,z){\displaystyle \frac{\mathrm{\Phi }(R,z^{})}{z^{}}}\text{d}z^{}`$ (15)
$`\overline{v}_\varphi ^2(R,z)`$ $`=`$ $`R{\displaystyle \frac{\mathrm{\Phi }(R,z)}{R}}`$ (16)
$`+`$ $`{\displaystyle \frac{R}{\rho (R,z)}}{\displaystyle \frac{}{R}}{\displaystyle _z^{\mathrm{}}}\rho (R,z^{}){\displaystyle \frac{\mathrm{\Phi }(R,z^{})}{z^{}}}\text{d}z^{}`$
As before, enlarging the parameter space to make more complicated assumptions about the stellar distribution function is not warranted given the other simplifications we adopt here. In any case, the results are not likely to be strongly affected unless the stellar distribution function is extremely anisotropic.

For our disc galaxy, we take the parameters used by Sofue (1996) for the Milky Way which assumes: (i) a spherically symmetric ($`a=0`$) nuclear star cluster of mass $`5\times 10^9\mathrm{M}_{}`$, and scale $`b=120\mathrm{pc}`$ (ii) a spherically symmetric bulge of mass $`10^{10}\mathrm{M}_{}`$ and scale $`b=750\mathrm{pc}`$; (iii) an axisymmetric disc of mass $`1.6\times 10^{11}\mathrm{M}_{}`$, radial scale $`a=6\mathrm{kpc}`$ and vertical scale $`b=500\mathrm{pc}`$; and (iv) a spherically symmetric dark halo of mass $`3\times 10^{11}\mathrm{M}_{}`$ and scale $`b=15\mathrm{kpc}`$.

#### 4.1.2 The satellite

As above, we model the satellite galaxy in this case as a Plummer sphere, subject to tidal stripping but with the addition of a central point mass to represent a nuclear black hole and surrounding star cluster. We increase the mass of the satellite to be $`4\times 10^9\mathrm{M}_{}`$ but retain the same scale length of $`b=0.4\mathrm{kpc}`$. This is consistent with the models of Sagittarius of Helmi and White (2001). This represents the inclusion of dark matter in the model for the dwarf. This could have be added into the previous simulations, but would have had little effect as the satellite was very rapidly destroyed by tidal stripping. We set the central point mass to be $`1\times 10^7\mathrm{M}_{}`$, although we shall see later that our simulations are basically independent of this.

### 4.2 Orbital dynamics

We model the orbital dynamics of the satellite in a manner using basically the same ideas as used above (Section 2.3). The satellite galaxy is subject to the gravitational influence of the main galaxy, to dynamical friction and to tidal stripping.

The gravitational force due to the main galaxy is straightforward to calculate. Due to the axisymmetry of our model, there are only two independent components to consider. The radial force component (in a cylindrical sense, that parallel to the plane of the disc) is

$$F_R=\frac{\mathrm{\Phi }}{R}=\frac{GMR}{\left(R^2+\left(a+\sqrt{b^2+z^2}\right)^2\right)^{3/2}}$$
(17)
and the vertical component (perpendicular to the disc),

$$F_z=\frac{\mathrm{\Phi }}{z}=\frac{GMz\left(a+\sqrt{b^2+z^2}\right)}{\left(R^2+\left(a+\sqrt{b^2+z^2}\right)^2\right)^{3/2}\left(a+\sqrt{b^2+z^2}\right)}.$$
(18)
The force in both directions is calculated for each of the four components of our model galaxy and these forces summed to produce the resultant force on the satellite due to the main galaxy.

To calculate the dynamical friction on the satellite, we first calculate the mean streaming velocity, $`\overline{v}_\varphi `$, of the disc component of the galaxy at the satellite’s current location using Equation 16 (note that all other components have zero mean streaming velocity since they are spherically symmetric) and from this we calculate the relative velocity of the satellite which respect to the disc. Second we compute the local velocity dispersion of each component of the galaxy. Then, finally, we use the standard gravitational drag formula (Equation 7) to calculate the dynamical friction force on the satellite. Once again, the contribution from all four components of the galaxy is then summed, and this added to the gravitational force to produce the total force on the satellite.

After each timestep in our integration, we apply simple tidal stripping to the satellite. We calculate the mass internal to the satellite’s current position, via

$$M(R,z)=2\pi _0^{\sqrt{R^2+z^2}}_0^\pi \rho (r,z)r^2\mathrm{sin}\theta \text{d}\theta \text{d}r$$
(19)
and then use this mass to calculate the current tidal radius of the satellite, as given in Equation 10. Whilst this is not strictly accurate, it is an adequate approximation for our current purposes. As we shall see below, the stripping of the satellite is not the most important factor in these simulations.

Again, we use a variant of the ‘odeint’ routine from Numerical Recipes \[Press et al. 1992\] to control the timestep in the simulation, so that we use small timesteps when necessary, but larger ones when possible to speed up our simulation.

### 4.3 Results

Since we know from our previous set of calculations that the only satellite trajectories which are likely to have a chance of reaching the centre of the host galaxy are those which are almost radial, we assume that the satellite galaxy is initially approaching the host galaxy from a large distance ($`100\mathrm{kpc}`$).

To be specific, we define the centre of the Seyfert as the origin of our coordinate system, and then consider the trajectories of small, incoming galaxies which start on the plane which has a normal vector at an angle $`90\alpha `$ from the axis of rotation of the Seyfert and has a point of closest approach of $`100\mathrm{kpc}`$ to the origin. We consider satellite orbits which start in this plane, and which have their initial velocity being perpendicular to the plane. We shall call the angle $`\alpha `$ the ‘approach angle’, as it corresponds to the angle between the initial velocity vector and the plane of the host galaxy. Thus $`\alpha =90\mathrm{°}`$ corresponds to an orbit which is initially directed at right angles to the disc of the host galaxy. For this paper, we investigated orbits with approach angles of $`0\mathrm{°}`$, $`15\mathrm{°}`$, $`30\mathrm{°}`$, $`45\mathrm{°}`$, $`60\mathrm{°}`$, $`75\mathrm{°}`$, $`85\mathrm{°}`$, and $`90\mathrm{°}`$ and with two sets of initial velocities: $`100\mathrm{km}\mathrm{s}^1`$ and $`200\mathrm{km}\mathrm{s}^1`$.

We continue to track the satellite’s orbit either until the satellite has settled into the plane of the disc, after which orbital evolution is very slow, or until the satellite reaches a radius, $`r_{\mathrm{sph}}`$ at which the potential of the host galaxy is spherical and so there will be no further changes in the direction of satellite’s orbital angular momentum. For our model, we take $`r_{\mathrm{sph}}`$ to be $`200\mathrm{pc}`$.

The main outcome of our orbital computations is that most satellites settle fairly quickly into the orbital plane of the disc galaxy. This comes about simply because the gravitational drag due to the plane is relatively large because the dispersion velocity of the disc component is relatively small. The fraction of the initial conditions which result in the satellite reaching the central regions of the main galaxy *before* it settles into the plane of the disc is relatively small, but not vanishingly so. Those that do reach the centre are generally those for which their first dynamically significant impact with the disc is in the centre of the galaxy. Three orbits, chosen such that two of them (the solid and dotted lines) are those which reach the centre of the galaxy, are shown in Figure 6. Once the satellite has settled into the plane its orbital evolution becomes very slow, and does not in general reach the centre within a Hubble time. We should note that for our approximations of a static potential for the disc to be valid, we require that the satellite not be more massive than the mass it interacts with when passing through the disc. The satellite will interact with material in a cylinder of radius of order of $`GM/v^2`$, where $`v`$ is the velocity of the satellite when it passes through the disc. This mass is always greater than the current mass of the satellite when passing through the disc, except for the earliest passes through the disc when there is no significant drag.

For each orbit which reached $`r_{\mathrm{sph}}`$ before settling into the plane of the disc, (i.e. which reached the centre in less than a Hubble time) we recorded the angle, $`\beta `$, which the satellite’s final orbital angular momentum vector made with the normal to the plane of the disc. Thus $`\beta =0\mathrm{°}`$ means the satellite’s angular momentum is parallel to that of the rotational angular momentum of the disc, i.e the final orbit of the satellite lies in the plane of the galaxy disc.

In general, the plane ($`\alpha =const.`$), from which we launch our initial trajectories contains a small number of regions (typically around 5) from which the satellite reaches the centre of the host galaxy, but these are spread out over a relatively large area. To find the cross-sectional area $`A\left(\alpha \right)`$ at each approach angle, we first undertook a low-resolution scan (out to a radius where satellite simply flew past the galaxy), to find the regions in which the satellite reached the centre before reaching the plane of the galaxy. We picked out the areas of interest by eye, and then investigated those regions at higher resolution. The cross-sectional area was then calculated simply by summing the area of these regions. The exception to this is the face-on case, which is rotationally symmetric in the plane of the disc. Here we simply undertook a high resolution scan along one line parallel to the plane of the disc and then calculated the cross-sectional area from the circular annuli produced. In Figure 7, we plot the contribution $`A\left(\alpha \right)\mathrm{cos}\alpha `$ to the total cross-sectional area made by each approach angle $`\alpha `$. The factor of $`\mathrm{cos}\alpha `$ accounts for the fact that we expect fewer satellites to have approach angles close to $`90\mathrm{°}`$, assuming the incoming flux of satellites is isotropic.

We have also determined the angle $`\beta `$ between the orbital angular momentum of the satellite’s orbit when it reaches the centre ($`r<r_{\mathrm{sph}}`$) and the normal to the plane of the galaxy disc. Figure 8 shows the distribution of $`\beta `$ for different approach angles $`\alpha `$. We note that there is little change with $`\alpha `$. The absence of orbits which end with a low value of $`\beta `$ is due to the fact that these orbits are close to the plane of the disc, and so quickly settle into the plane of the disc and are never reach the radius $`r_{\mathrm{sph}}`$.

We also compute an overall distribution for the angle $`\beta `$ for those mergees which reach the centre of the host galaxy. In doing so we encounter the problem that the distribution of the initial orbits of the satellites’ orbits is unknown. In the absence of any strong evidence to the contrary, we assume that the incoming merging satellites approach from random directions. With this assumption, we have computed the distribution for $`\beta `$, and this is as shown in Figure 9. We note that this distribution is reasonably uniform in $`\beta `$, except for a zone of avoidance for angles $`\beta \mathrm{}<20\mathrm{°}`$.

We have repeated the above experiments with the black hole in the centre of the dwarf galaxy having a mass from $`10^5\mathrm{M}_{}`$ to $`10^8\mathrm{M}_{}`$ (as opposed to the $`10^7\mathrm{M}_{}`$ used previously), and obtain no significant differences. This may have been expected from the fact that the satellite is not stripped down to a small mass by the time it reaches a radius of $`200\mathrm{pc}`$ from the centre of the Seyfert.

### 4.4 Comparision with observations

A comparision of our model with the available observational data is not entirely trivial. In order to obtain a quantitative result, we follow the method of Kinney et al (2000, and references therein). A brief description of the method is given here.

Two parameters are observed for every galaxy used in Kinney et al.: $`i`$, the inclination of the galaxy’s disc to the plane of the sky, and $`\delta `$, the difference between the position angles of the galaxy’s major axis and that of the jet when projected onto the plane of the sky. Each $`(i,\delta )`$ pair constrains the jet to lie on a pair of great circles on a sphere about the nucleus of the galaxy. The angle $`\beta `$ between the jet and the rotation axis of the disc varies around these great circles. We assume that the chance of a jet being at any point on either of these great circles is given by the probability of our simulations producing the value of $`\beta `$ (as shown in Figure 9) for that point on the great circle. At a fixed inclination, different values of $`\delta `$ give great circles which have different orientations in space, and hence we calculate the probability of obtaining a given value of delta given the inclination, $`p\left(\delta |i\right)`$. We then integrate this distribution from $`\delta =0`$ to $`\delta =\delta _k`$ where $`\delta _k`$ is the observed value of $`\delta `$ for the galaxy $`k`$ to obtain the centile point $`c_k`$. The distribution of centile points across the data set should then be uniform over the range $`[0,1]`$. Full details of this technique can be found in Kinney et al. (2000).

We perform a K-S test to compare our model with the the $`60\mu \text{m}`$ data of Kinney et al. (2000). We find that the data is consistent with being drawn from our model at the 97% level.

## 5 Discussion

We have undertaken an investigation into the possibility that the activity of, and jet direction in, at least some active galactic nuclei might be fuelled by the infall and merger of small satellite galaxies.

### 5.1 Cross-sections for mergers

We considered first the evidence that many of the host galaxies active nuclei in the 3CR catalogue are found to have discs of dust around their nuclei. Taking this as prima facie evidence that the nuclear activity might be related to infall, we have modelled the disruption of a small satellite galaxy by a larger elliptical, taking the particular case of NGC 3379, which while not in the 3CR catalogue, does show a radio bright nucleus and does display a dust disc.

Our modelling procedure is of necessity somewhat crude. We have adopted a realistic model for the inner regions of the host elliptical where the disruption of the infalling satellite takes place, but have not taken proper account of a possible massive, but low density, dark halo. We considered idealised models for the internal structure of the satellite, one with and one without a central density cusp, but found that the results did not depend significantly upon which models we used. Gravitational drag on the orbit of the satellite was modelled approximately using the standard Chandrasekhar formula, suitably parameterized by comparison with full numerical computations taken from the literature.

Probably the least realistic part of our modelling procedure is the description of the disruption of the satellite and the subsequent dynamics of the gas it contained. For example, we have assumed that the main galaxy is sufficiently devoid of gas that ram pressure stripping of gas from the satellite can be ignored, and that the effect of the interactions between the stripped gas and the gas already present may be ignored. This may not always be the case. Even so, because our description does contain much of the relevant physics it is unlikely that our results will be grossly in error. The infalling satellite is disrupted at the point where tidal effects of the host come into play using simple truncation of the underlying model. This ignores the dynamical response of the satellite to the truncation process, and so, if anything, tends to overestimate the radius at which disruption occurs. On the other hand, we have implicitly assumed that the gas and stars in the satellite occupy the same volume, and thus strip gas particles in proportion to the mass removed. If the gas were more centrally located than the bulk of the mass, then our procedure would underestimate the radius at which the gas from the satellite is deposited. However, as can be seen from Figure 2, the radius at which tidal disruption takes place is fairly well-defined, with the radius by which 90 per cent of the satellite has been tidally removed being almost independent of the model for the internal structure. The only way to bring about a significant change in the radius at which the gas is deposited is to use the cusped model and to assume that all the gas lies within the radius occupied by the central one percent of the total galaxy mass. This would imply for a gas to mass ratio of a percent that the gas was dynamically important in the nuclear region of the satellite. The modelling of the dynamics of the gas once it has been stripped conserves mass and angular momentum, but uses a simple algorithm to allow intersecting streams of gas to shock and lose energy by radiation. This is probably an adequate approximation for the current paper.

Despite the relative crudeness of our modelling procedure, we are able to draw some conclusions which we believe are relatively robust. We find that while it is possible to reproduce the dust disc structure seen in NGC 3379 via the infall of a satellite galaxy, it seems likely that this type of event is very rare. Our initial simple model for the satellite, which started within the galactic halo on an orbit with a given eccentricity, was not able to reproduce the structure seen in NGC 3379. A more realistic model, in which the satellite starts well outside the galaxy with a small impact parameter was found to be capable of reproducing the observed structure, provided that the satellite galaxy has an initial angular momentum about the galactic nucleus of around $`50\text{kpc km s}^1`$. Given the lack of observational or theoretical evidence as to the distribution of the orbits of dwarf satellites, we cannot make a quantitative statement about how likely this is, but we do not believe that this event is common. If we assume a typical velocity of around $`200\mathrm{km}\mathrm{s}^1`$ for the satellite when at infinity, this means the initial orbit must have an impact parameter of the order of $`250\mathrm{pc}`$, which would represent a small fraction of phase space if the initial orbits were uniformly distributed. Even if the orbital distribution were moderately biased towards low angular momentum orbits, this would still imply that the formation of a structure similar to that in NGC 3379 being formed by this method is a rare event.

With regard to powering a single AGN, this model produces results which are not inconsistent with current theories. To order of magnitude level, if an event forming a disc similar to that observed in NGC 3379 occurs, we find that a few percent of the gas particles reach a distance of $`10\mathrm{pc}`$ from the nucleus. At this point it is reasonable to assume that they will be swallowed by the central black hole within a reasonably short time span \[Shlosman et al. 1990, Bekki 2000\]. If we assume that a few percent of the mass of the original dwarf galaxy was gas, this gives around $`10^6\mathrm{M}_{}`$ being swallowed by the black hole, which, at an accretion rate of $`0.1\text{M}_{}/\text{yr}`$, will power a low level AGN for around $`10^7\mathrm{yr}`$. This is not unreasonable.

However, if we further assume that all galaxy activity is triggered only by events of the above type, then the model becomes harder to sustain. If we assume that every merger triggers nuclear activity which lasts around $`10^7`$ years, and if around $`1\%`$ of elliptical galaxies are active, then this implies that each elliptical galaxy has undergone about 10 such events over the lifetime of the universe. However, if only a few percent of mergers occur at sufficiently low angular momentum to give rise to nuclear activity, then the total number of satellites swallowed by the host over its lifetime rises to around 300 to 1000. If each satellite has a mass of around $`10^9\mathrm{M}_{}`$ as we assume here, this implies the addition of a substantial amount of mass. Even if our estimates are too pessimistic with regard to the efficiency with which merger process is able to fuel the nucleus, it is clear that we would nevertheless expect that essentially all elliptical galaxies should show some evidence of such a merger having occurred, and perhaps more than one. In this regard, we note that whilst van Dokkum and Franx \[van Dokkum & Franx 1995\] found dust structures in approximately half of their sample, more recent more sensitive surveys \[Pastoriza et al. 2000\] are finding dust in a higher proportion of elliptical galaxies, including some in which dust remained undetected in the van Dokkum and Franx survey.

### 5.2 Directionality

Having established that the likelihood of any particular merging satellite galaxy reaching the host nucleus is relatively small, we then turn our attention to incoming orbits for which the successful mergees actually reach the nucleus. The motivation for this part of the investigation is the discovery (Kinney et al. 2000) that the jet directions in the centres of Seyfert galaxies appear to bear no relation to the orientation of the discs of the host galaxies. The possibility we are investigating here is whether the mergees reaching the centre do so with an angular momentum distribution which could produce the observational data observed by Kinney et al.

For this case we have taken the host galaxy to be a spiral galaxy, and have employed the same procedures for orbital decay and tidal stripping as before. Again we find that the cross-sectional areas for merging of the satellite with the host nucleus are relatively small, but we now concentrate our analysis on the angular momentum of the orbits of those satellites (or what is left of them) when they reach the nuclear regions of the host. Those satellite nuclei which do reach the centre in a finite time (that is, before settling into the disc) have the angle between their angular momentum vector and the disc normal ($`\beta `$) distributed approximately randomly, as is shown in Figure 9. The deficit at low values of $`\beta `$ comes about because these satellites are close to, or in, the disc of the Seyfert galaxy, and as such rapidly settle down into the plane. We have noted that is distribution is in good agreement with the data presented in Kinney et al. (2000).

The simulations we have performed here give us a distribution for the angle $`\beta `$. However, in order to relate these results to the observed jet distribution, we need the resultant spin of the black hole formed between the merger of the two holes. For the general case, this is an unsolved problem in general relativity and hence we again adopt a simple prescription. We assume the mass of the dwarf’s black hole, $`M_D`$, is much less than the mass of the Seyfert’s black hole, $`M_{\mathrm{AGN}}`$ and can then use the prescription from Colbert and Wilson (1995) to first order to obtain the spin of the daughter of the merger as

$$L_S=\frac{2\sqrt{3}GM_DM_{\mathrm{AGN}}}{c}$$
(20)
where $`L_S`$ is the spin of the resultant hole. This applies only for a merger between two black holes which were not originally spinning, which we note is not the case here. If $`M_D/M_{\mathrm{AGN}}=0.1`$, then the resultant hole is a Kerr black hole with spin parameter $`a=0.35`$, where $`a`$ is the angular momentum of the hole in units of $`GM^2/c`$ and lies between 0 and 1. Thus the accumulation of a few satellite nuclear black holes by the black hole in the AGN nucleus, from orbital directions which are more of less randomly distributed can lead to random orientations for the spins of the AGN nuclear black holes. A more detailed discussion of this is given by Merritt (2002).

## 6 Conclusions

We have investigated the dynamics of the merging process in the minor merger hypothesis for active galactic nuclei. Our analysis has resulted in two main conclusions.

First, for the satellite galaxy to be able to merge directly with the nucleus of the host galaxy (for example, to give rise to the compact dust discs which are seen in early type active galaxies) requires the initial orbit of the satellite to be well aimed. The corollary of this is that for each merging satellite which gives rise to a compact nuclear dust disc, there must be many which do not. Similarly, in disc galaxies, for every minor merger which causes an active nucleus, there must be many mergers which simply give rise to a merger of the incoming satellite with the galaxy disc.

Second, in the case of the host galaxy being a disc galaxy, if the initial orbits of the satellites are randomly oriented with respect to the host galaxy, then the orbits of those which reach the host nuclear regions in a reasonable time, are also fairly randomly oriented once they reach the nucleus. We note that this result might be able to provide an explanation of why the jet directions in the nuclei of Seyfert galaxies are apparently unrelated to the plane of the galaxy disc. The distribution of angles $`\beta `$ between the final orbit of the successful merger candidates and the Seyfert disc which we find once we have averaged over an assumed isotropic incoming satellite distribution is consistent with the observational data presented by Kinney et al (2000).

## 7 Acknowledgments

We thank Frank van den Bosch for providing us with the diagrams in Figure 1. We would also like to thank an anonymous referee for constructive comments made on earlier versions of this paper.