text:
# Relaxation of finite perturbations: beyond the Fluctuation-Response relation

## Lead Paragraph

Understanding the behavior of a dynamical system out of its equilibrium is a crucial issue of statistical physics. In the case of an infinitesimal perturbation that shifts the system out of equilibrium, the classical Fluctuation-Response theorem allows to determine the linear response of the system in term of its equilibrium properties, i.e. correlation functions. While the behavior of infinitesimal perturbations gives relevant information for problems of statistical mechanics, for climate and geophysical models the main goal is to characterize the relaxation of large perturbations, which can not be obtained from the linear response theorem. We present here a generalization of Fluctuation-Response relation, which holds for finite amplitude perturbations, providing a tool for extracting non equilibrium behavior out of equilibrium features of the system. We also discuss the non trivial role of the amplitude of perturbations in systems where many characteristic time scales are present.

## I Introduction

The Fluctuation-Response (F/R) relation has a deep relevance in statistical physics and more generally in systems with chaotic dynamics (in particular in hydrodynamics ). The relevance of a connection between “non equilibrium” features (i.e. response to an external perturbation) and “equilibrium” properties (i.e. time correlations computed according to the invariant measure) is well known in statistical mechanics. We can mention the important Green - Kubo formulas in the linear response theory . Beyond statistical physics, an other field where the F/R problem has an obvious relevance is climate research . One of the key problems is the possibility to understand the response of the present climate to some violent changes (e.g. a volcanic eruption). The essential point is the possibility that the recovery of the climate system from a perturbation (response) can be estimated from its time history (correlations time of the unperturbed system).

Assuming that the system is mixing and has invariant probability density function (pdf) $`\rho (𝒙)`$, it is possible to derive the following F/R relation. Let us denote by $`𝒙(t)=(x_1(t),\mathrm{},x_N(t))`$ the state of the system at time $`t`$. If at the initial time $`t=0`$ the system is perturbed by $`\delta 𝒙(0)=(\delta x_1(0),\mathrm{},\delta x_N(0))`$, the average evolution of the perturbation $`\delta x_i(t)`$ with respect the unperturbed trajectory is

$$\delta x_i(t)=\underset{i}{}R_{i,j}(t)\delta x_j(0)$$
(1)
where

$$R_{i,j}(t)=\frac{\delta x_i(t)}{\delta x_j(0)}=x_i(t)f_j(𝒙(0))$$
(2)
and the function $`f_j`$ depends on $`\rho (𝒙)`$ as

$$f_j(𝒙)=\frac{\mathrm{ln}\rho (𝒙)}{x_j}$$
(3)
In Section II we will give a complete derivation of the above formulas.

As far as we know, the F/R problem had been studied only for infinitesimal perturbations. For statistical mechanics problems it is relevant to deal with infinitesimal perturbations on the microscopic variables. In a similar way this problem has importance in many analytical approaches to the statistical description of hydrodynamics where Green functions are naturally involved both in perturbative theory and closure schemes .

On the other hand in geophysical or climate problems the interest for infinitesimal perturbation seems to be rather academic, while the interesting problem is the behavior of relaxation of large fluctuations in the system due to fast changes of the parameters.

In this paper we want to address the problem of the F/R relation for non infinitesimal perturbations. In Section II we will show that it is possible to generalize the F/R relation to large perturbations, involving rare events of the invariant measure. Section III is devoted to a discussion on the connections, and differences, between our approach and well known results in dynamical system theory. In Section IV we will discuss the application to systems involving a single characteristic time, while Section V is devoted to system with many characteristic times. Section VI is devoted to conclusions and the Appendix VII contains some technical remarks.

## II Theoretical background

In the following we will consider a dynamical system with evolution $`𝒙(t)=\varphi ^t𝒙(0)`$ of the $`N`$-dimensional vector $`𝒙`$. For generality, we will explicitly consider the case in which the time evolution can also be not completely deterministic (e.g. stochastic differential equations). We will assume the existence of an invariant probability distribution $`\rho (𝒙)`$ and the ergodicity of the system so that

$$A\underset{T\mathrm{}}{lim}\frac{1}{T}_0^TA(𝒙(t))𝑑t=\mu (A)A(𝒙)\rho (𝒙)𝑑𝒙$$
(4)
for any (smooth enough) observable $`A`$.

Our aim is the understanding of the mean response $`\delta A(t)`$ of a generic observable $`A`$ initially perturbed with $`\delta A(0)`$. The first step is the study of one component of $`𝒙`$ i.e. $`\delta x_i(t)`$ with an initial (non-random) perturbation $`\delta 𝒙(0)=\delta 𝒙_0`$. Introducing the probability of transition from $`(𝒙_0,0)`$ to $`(𝒙,t)`$, $`W(𝒙_0,0𝒙,t)`$ (for a deterministic system we have $`W(𝒙_0,0𝒙,t)=\delta (𝒙(t)\varphi ^t𝒙_0)`$), we can easily write an expression for the mean value of the variable computed along the perturbed trajectory $`x_i^{}(t)=x_i(t)+\delta x_i(t)`$:

$$x_i^{}(t)=x_i\rho ^{}(𝒙_0)W(𝒙_0,0𝒙,t)d𝒙d𝒙_0$$
(5)
where $`\rho ^{}(𝒙)`$ is the initial distribution of perturbed system, which is related to the invariant distribution by $`\rho ^{}(𝒙_0)=\rho (𝒙_0\delta 𝒙_0)`$. Noting that the mean value of $`x_i(t)`$ can be written in a similar way:

$$x_i(t)=x_i\rho (𝒙_0)W(𝒙_0,0𝒙,t)d𝒙d𝒙_0$$
(6)
one has:

$`\delta x_i(t)`$ $`=`$ $`{\displaystyle x_i\frac{\rho (𝒙_0\delta 𝒙_0)\rho (𝒙_0)}{\rho (𝒙_0)}\rho (𝒙_0)}`$ (7)
$`\times `$ $`W(𝒙_0,0𝒙,t)d𝒙d𝒙_0`$ (8)
$`=`$ $`x_i(t)F(𝒙_0,\delta 𝒙_0)`$ (9)
where

$$F(𝒙_0,\delta 𝒙_0)=\left[\frac{\rho (𝒙_0\delta 𝒙_0)\rho (𝒙_0)}{\rho (𝒙_0)}\right]$$
(10)
For an infinitesimal perturbation $`\delta 𝒙(0)=(\delta x_1(0)\mathrm{}\delta x_N(0))`$ expanding (10) to first order one ends with the expression

$`\delta x_i(t)`$ $`=`$ $`{\displaystyle \underset{j}{}}x_i(t){\displaystyle \frac{\mathrm{ln}\rho (𝒙)}{x_j}}|_{t=0}\delta x_j(0)`$ (11)
$``$ $`{\displaystyle \underset{j}{}}R_{i,j}(t)\delta x_j(0)`$ (12)
which defines the linear response

$$R_{i,j}(t)=x_i(t)\frac{\mathrm{ln}\rho (𝒙)}{x_j}|_{t=0}$$
(13)
of the variable $`x_i`$ with respect to a perturbation of $`x_j`$. Relation (13) is the generalization for non Hamiltonian systems of the well known fluctuation/response (F/R) relation .

Let us note that in the general case the invariant measure $`\rho (𝒙)`$ is not known, so the equation (13) gives just a qualitative information. In the case of Gaussian distribution, $`\rho (𝒙)`$ factorizes and the linear response recovers the correlator

$$R_{i,j}(t)=\frac{x_i(t)x_j(0)x_ix_j}{x_jx_jx_jx_j}$$
(14)
In the case of finite perturbations, the F/R relation (9) is typically non-linear in the perturbation $`\delta 𝒙_0`$ and thus no simple relations analogous to (13) exist. Nevertheless we can disentangle the different contributions in the response (9) by studying an initial perturbation whose only non-zero component is the $`jth`$ one,

$$\delta ^{(j)}𝒙(0)=(0,\mathrm{},0,\delta x_j(0),0,\mathrm{},0)$$
(15)
We therefore generalize the F/R relation (13) to non-linear response of $`x_i`$ to a perturbation on the $`j`$ variable as

$$R_{i,j}(t)=x_i(t)f_j(0)$$
(16)
where $`f_j`$ is given by

$$f_j(𝒙_0)=\frac{\rho (𝒙_0\delta ^{(j)}𝒙(0))\rho (𝒙_0)}{\rho (𝒙_0)\delta x_j(0)}$$
(17)
The explicit prediction of the response from (16) requires the analytic expression of the invariant pdf, which is in general not known. Nevertheless (9) guarantees the existence of a link between equilibrium properties of the system and the response to finite perturbations. This fact has a relevant consequence for systems with one single characteristic time: a generic correlation (e.g. the correlation (14)) in principle gives informations on the relaxation time of finite size perturbations, even when the invariant measure $`\rho `$ is not known .

## III Remarks on the connections between F/R relation, dynamical system theory and statistical mechanics

Since the F/R relation involves the evolution of differences between variables computed on two different realizations of the system, it is natural to conclude that this issue is closely related to the predictability problem and, more in general, to chaotic behavior. Actually, a detailed analysis shows that the two problems, i.e. F/R relation and predictability, have only a very weak connection. For the sake of completeness, we shortly discuss here the analogies and differences between these two issues.

The typical problem for the characterization of predictability is the evolution of the trajectory difference $`\delta 𝒙(t)`$, in particular of $`\mathrm{ln}|\delta 𝒙(t)|`$ which defines the leading Lyapunov exponent $`\lambda `$. For small $`|\delta 𝒙(0)|`$ and large enough $`t`$ one has

$$\mathrm{ln}|\delta 𝒙(t)|\mathrm{ln}|\delta 𝒙(0)|\lambda t$$
(18)
On the other hand, in F/R issue one deals with averages of quantities with sign, such as $`\delta 𝒙(t)`$. This apparently marginal difference is very important and it is at the basis of the famous objection by van Kampen related to the standard derivation of the linear response theory . In a nutshell, using the modern dynamical systems terminology, the van Kampen’s argument is as follows. Since in presence of chaos $`|\delta 𝒙(t)|`$ grows exponentially in time, it is not possible to linearize (10) for time larger than $`(1/\lambda )\mathrm{ln}(\mathrm{\Delta }/|\delta 𝒙(0)|)`$, where $`\mathrm{\Delta }`$ is the typical fluctuation of the variable $`𝒙`$. As a consequence, the linear response theory is expected to be valid only for extremely small and unphysical perturbations, in clear disagreement with the experience. A solution of this apparent paradox was proposed by Kubo who suggested that “instability \[of the trajectories\] instead favors the stability of distribution functions, working as the cause of the mixing. More recent works have demonstrated the constructive role of chaos in F/R relation and the non relevance of van Kampen’s criticism . The objection by van Kampen remains nevertheless relevant for numerical computations of F/R relation (see Appendix).

Fluctuation/response relation was developed in the context of statistical mechanics of Hamiltonian systems, but it also holds for non conservative systems, and even non deterministic systems (e.g. Langevin equations) and has no general relation with “chaotic quantities” such as Lyapunov exponents or Kolmogorov-Sinai entropy. This generated in the past some confusion about the applicability of F/R/ relation. For example, some authors claimed (with qualitative arguments) that in fully developed turbulence there is no relation between equilibrium fluctuations and relaxation to equilibrium  while the correct statement concerns the non validity of the simplified relation (15) which holds only for systems with Gaussian statistics.

Thanks to its general validity and robustness, the F/R relation has been also used to obtain informations on the unknown invariant measure $`\rho (𝒙)`$ on the basis of the linear response $`R_{i,j}(t)`$. An important example comes from the field of disordered systems where the F/R had been applied to the study of aging phenomena .

Concluding this short discussion on the connections between F/R relation, dynamical system theory and statistical mechanics, we mention recent results about rigorous derivation of the Onsager reciprocity relations  and the macroscopic fluctuation theory for stationary non-equilibrium states  in a class of stochastic models describing interacting particles systems.

## IV Systems with a single characteristic time

Let us start by studying two examples of systems with a single characteristic time: a deterministic chaotic system (the Lorenz model) and a nonlinear Langevin process.

We first consider the Lorenz model

$`{\displaystyle \frac{dx}{dt}}`$ $`=`$ $`\sigma (yx)`$ (19)
$`{\displaystyle \frac{dy}{dt}}`$ $`=`$ $`xz+rxy`$ (20)
$`{\displaystyle \frac{dz}{dt}}`$ $`=`$ $`xybz`$ (21)
with standard parameters for chaotic behavior: $`b=8/3`$, $`\sigma =10`$ and $`r=28`$. The correlation function (14) for the variable $`z`$, shown in Fig. 1, qualitatively reproduces the behavior of the response to different sizes of the perturbation of the $`z`$ variable, ranging from infinitesimal ones up to the size of the attractor. The accuracy does not increase when decreasing the perturbation because the invariant distribution is not Gaussian (see Fig. 1) and thus the general correlation (13) should be used. We observe that the use of (13) instead of (14) is in general much more difficult because the invariant distribution is in general non factorable.

To better illustrate this point, let us now consider a system whose invariant probability distribution is known. In this case we can quantitatively compare the differences between the responses to infinitesimal and finite perturbations. Our example is provided by the stochastic process $`x(t)`$ determined by

$$\frac{dx}{dt}=\frac{dU(x)}{dx}+\sqrt{2D}\xi (t)$$
(22)
where $`\xi (t)`$ is a white noise, i.e. a Gaussian process with $`\xi (t)=0`$ and $`\xi (t)\xi (t^{})=\delta (tt^{})`$. The invariant probability distribution is :

$$\rho (x)=𝒩e^{U(x)/D}$$
(23)
where $`𝒩`$ is fixed by normalization.

A Gaussian pdf is obtained using $`U(x)=x^2/2`$ which corresponds to the linear Ornstein-Uhlenbeck process $`dx/dt=x+\sqrt{2D}\xi (t)`$. Our example uses a modified version of the Gaussian case,

$$U=\{\begin{array}{cc}\frac{1}{2}x^2\hfill & ,|x|<1\hfill \\ |x|\frac{1}{2}\hfill & ,|x|>1.\hfill \end{array}$$
(24)
The resulting pdf, shown in the inset of Fig. 2, has a Gaussian core, with exponential tails. Figure 2 also shows the response function for an infinitesimal and for a finite size perturbation. For both perturbations, the response function measured from the perturbed trajectories is exactly predicted by statistics of the unperturbed system according to (16), while the Gaussian correlation $`C(t)=x(t)x(0)/\sigma ^2`$ gives only an estimate of the relaxation time. By construction, the pdf of this system has larger tails than in the Gaussian case, thus large fluctuations decay slower than small ones. In the linear case the mean response is simply $`R(t)=\mathrm{exp}(t)`$ and does not depend on the amplitude of the initial perturbation $`\delta x(0)`$.

The results obtained for the Lorenz model and for the nonlinear Langevin equations suggest that if only one characteristic time is present, the existence of the F/R relation allows for some qualitative results even in the absence of precise knowledge of $`\rho `$, both for infinitesimal and finite perturbation.

## V Systems with many characteristic times

In systems with many characteristic times, different correlation functions do not show the same behavior, i.e. depending on the observable one can observe very different time scales, corresponding to the different decay times of the correlation functions $`C_{j,j}=x_j(t)x_j(0)`$ . In addition, at variance with systems with one single time scale, here the amplitude of the perturbation can play a major role in determining the response, because different amplitudes may affect features with different time properties.

The link between equilibrium and relaxation properties established by the F/R relation (16) suggests that it is possible to relate different relaxation rates with the time scales measured by means of correlations. Consider the case of an observable $`A`$ which depends on all the variables of the system $`\{x_1,\mathrm{},x_N\}`$. For infinitesimal perturbations, a straightforward generalization of (1,2) gives:

$$\delta A(t)=A\left(𝒙(t)\right)f_j\left(𝒙(0)\right)\delta x_j(0)$$
(25)
In the case of finite perturbations, as stressed in Sect. (II), it is possible to write a F/R relation:

$$\delta A(t)=A\left(𝒙(t)\right)F(𝒙(0),\delta 𝒙(0))$$
(26)
in which, at variance with (25), all the variables are mixed. In (26) the relaxation properties depend explicitly on the initial perturbation $`\delta 𝒙(0)`$.

Depending on the choice of $`A(𝒙)`$, different perturbations on $`A`$ correspond to different amplitudes of the perturbations on each variable $`x_j`$. Consequently, one can think that it is possible to associate each perturbation to a certain subset of variables which are mainly perturbed. The relaxation of $`\delta A(t)`$ will be ruled by the characteristic time of that particular subset.

In order to illustrate this issue we consider a shell model for turbulence . Shell models are a simplified model for turbulent energy cascade, that describe the dynamics of velocity fluctuations at a certain scale $`\mathrm{}_n=k_n^1`$ with a single shell-variable $`u_n`$. Wave-numbers $`k_n`$ are geometrically spaced as $`k_n=k_0\lambda ^n`$, allowing to cover a large range of scales with relatively few variables. A quadratic interaction between neighbor shell reproduces the main features of three-dimensional turbulence. The specific model we will use is

$`\left({\displaystyle \frac{d}{dt}}+\nu k_n^2\right)u_n`$ $`=`$ $`i[k_{n+1}u_{n+1}^{}u_{n+2}ϵk_nu_{n+1}u_{n1}^{}`$ (27)
$`+`$ $`(1ϵ)k_{n1}u_{n2}u_{n1}]+f_n`$ (28)
where $`\nu `$ is the molecular viscosity, $`f_n`$ is an external forcing which injects energy at large scale, and $`ϵ`$ is a free parameter. In order to have the correct conservation laws (energy and helicity) in the inviscid unforced case one has to fix $`ϵ=1/2`$. The observable considered is the total energy $`E(t)=\frac{1}{2}_{n=1}^N|u_n(t)|^2`$ which is the conserved quantity in the inviscid, unforced limit .

In order to study the response to perturbations with different amplitude on $`E`$, we consider the following perturbed systems labeled with $`i=1,..,N`$: $`u_n^{(i)}(t)=u_n(t)+\delta u_n^{(i)}(t)`$ where the initial perturbations $`\delta u_n^{(i)}(0)`$ are set in the following way:

$$\delta u_n^{(i)}(0)=\{\begin{array}{cc}0\hfill & ,\mathrm{\hspace{0.17em}1}ni1\hfill \\ \sqrt{|u_n|^2}\hfill & ,inN\hfill \end{array}$$
(29)
This corresponds to a set of initial perturbations of the energy

$$\delta E_i(0)=\frac{1}{2}\underset{n=i}{\overset{N}{}}|u_n|^2$$
(30)
Such a perturbation is motivated by the fact that in the unperturbed system the energy is distributed among the shells according to the Kolmogorov scaling $`|u_n|^2k_n^{2/3}`$, and the smaller scales give smaller contributions to the energy $`E(t)`$. Thus it is natural to assume that a small perturbation of the energy will affect mainly the small scales.

For each perturbation $`\delta E_i`$, the average response of energy

$$\frac{\delta E_i(t)}{\delta E_i(0)}=\frac{\underset{n=1}{\overset{N}{}}|u_n^{(i)}(t)|^2|u_n(t)|^2}{_{n=1}^N|u_n^{(i)}(0)|^2|u_n(0)|^2}$$
(31)
reveals a close relation with the time correlation of the corresponding largest perturbed shell $`u_i(t)`$, as shown in Fig. 3. A measure of the relaxation time can be provided by the halving times $`T_{1/2}`$ of the mean response, at which $`\delta E_i(T_{1/2})=1/2\delta E_i(0)`$. The dependence of response times on the amplitude of the initial perturbation, shown in Fig.4, reflects Kolmogorov scaling for characteristic times $`\tau _nk_n^{2/3}u_n^2\delta E_n`$

$$T_{1/2}\delta E.$$
(32)
The above results on the shell model show that the response to a finite size perturbation of a system with many characteristic times may depend on the amplitude of the perturbation. Thanks to the existence of F/R relation is possible to establish a link between relaxation times of different perturbation and characteristic times of the system.

## VI Conclusions

Starting from the seminal works of Leith , who proposed the use of F/R relation for understanding the response of the climatic system to changes in the external forcing, many authors tried to apply this relation to different geophysical problems, ranging from simplified models , to general circulation models  and to the covariance of satellite radiance spectra . In most of the applications it has not been taken into account the limits of applicability of the F/R relation which has been used as a kind of approximation. We have shown that a F/R relation holds under very general conditions. The derivation in Section II clearly shows the limits of applicability in its simplest form (i.e. the Gaussian approximation (14)).

Our main result is the demonstration that an exact fluctuation/response relation holds also for non infinitesimal perturbation. This relation involves the detailed form of the invariant probability distribution. In particular, in order to predict the mean response to large perturbations, one needs a precise knowledge of the tails of the pdf.

We believe that this generalization of the usual linear response theory can be relevant in many applications. As an example, we can mention climate research, where our results imply the possibility, at least in principle, to understand the behavior of the system after a large impulsive perturbation (e.g. a volcanic eruption) in terms of the knowledge obtained from its time history. Of course one has to take into account the strong limitations due to the need to have a good statistics of rare events.

###### Acknowledgements.

This work has partially supported by MIUR (Cofinanziamento Fisica Statistica di Sistemi Complessi Classici e Quantistici). We acknowledge the allocation of computer resources from INFM Progetto Calcolo Parallelo. Angelo Vulpiani acknowledges support from the INFM Center for Statistical Mechanics and Complexity (SMC). We thank M. Falcioni for useful remarks.

## VII Appendix

In this appendix we want to discuss how van Kampen criticism is relevant for the numerical evaluation of infinitesimal response function.

In numerical simulations, $`R_{i,j}(t)`$ is computed perturbing the variable $`x_i`$ at time $`t=t_0`$ with a small perturbation of amplitude $`\delta x_i(0)`$ and then evaluating the separation $`\delta x_i(t)`$ between the two trajectories $`𝒙(t)`$ and $`𝒙^{}(t)`$ which are integrated up to a prescribed time $`t_1=t_0+\mathrm{\Delta }t`$. At time $`t=t_1`$ the variable $`x_i`$ of the reference trajectory is again perturbed with the same $`\delta x_i(0)`$, and a new sample $`\delta 𝒙(t)`$ is computed and so forth. The procedure is repeated $`M1`$ times and the mean response is then evaluated according to (14).

In presence of chaos, the two trajectories $`𝒙(t)`$ and $`𝒙^{}(t)`$ typically separate exponentially in time and the perturbed system relaxes to the unperturbed one only in average, therefore the mean response is the result of a delicate balance of terms which grow in time in different directions. The average error in the computation of $`R_{i,j}(t)`$ typically increases in time as $`e^{L(2)t/2}/\sqrt{M}`$, where $`L(2)`$ is the generalized Lyapunov exponent . Thus very high statistics is needed in order to compute $`R_{i,j}(t)`$ for large $`t`$ .

We remark that the exponential growth is generally valid only for infinitesimal perturbation. When the perturbation reaches the typical size of the system, the difference between the perturbed and the unperturbed trajectory tends to saturate. Thus, for finite amplitude perturbations the mean response is the average of terms that remain of order $`O(1)`$, and less statistics is required to obtain convergence. In this sense the mean response to finite perturbation is more representative of the behavior of a single perturbation than in the infinitesimal case.

On the other hand, even if equation (16) is formally valid for arbitrary large perturbations, for practical use an upper limit exist due to finiteness of statistics. To predict the relaxation of a perturbation $`\delta 𝒙(0)`$, one needs sufficient statistics for the convergence of $`\rho (𝒙(0)\delta 𝒙(0))`$. This request is more severe in systems where large fluctuations are suppressed. An example is provided by the stochastic model (22) with

$$U(x)=\frac{1}{2}x^2+\frac{1}{4}x^4$$
(33)
Here the pdf has sub-Gaussian tails, and we observe the opposite behavior of the system (24), as shown in Fig. 5. While in the case with exponential tails we have a good statistical convergence for a perturbation greater than $`2\sigma `$ in the second system this perturbation is too large to obtain convergence even with huge statistics ($`10^9`$ runs).


format: text

text:
# On the central black hole mass in Mkn 501

## 1 Introduction

According to the commonly accepted paradigm, Active Galactic Nuclei (AGN) are thought to harbour supermassive black holes (BH) surrounded by geometrically thin accretion disks, the latter triggering the formation of relativistic jets. Increasing evidence now indicates, however, that this picture probably has to be expanded: First, hierarchical galaxy evolution schemes suggest, that as a result of mergers between galaxies, binary BH systems (BBHSs) should be generally expected in the center of elliptical galaxies, and hence, for example, in the typical hosts of BL Lac objects (cf. Begelman et al. 1980; Kauffmann 1996; Richstone et al. richstoneetal1998 (1998); Milosavljevic & Merritt milo01 (2001); Yu yu2002 (2002)). Secondly, from a phenomenological point of view, a multitude of observational evidence has been plausibly related to the presence of BBHSs in AGN, including the observed misalignment, precession and wiggling of extragalactic jets (e.g. Conway & Wrobel conwaywrobel1995 (1995); Kaastra & Roos kaastraroos1992 (1992)), periodic optical outburst events as observed in the quasar OJ 287 (e.g. Sillanpää et al. sillanetal1988 (1988); Lehto & Valtonen letho1996 (1996); Villata et al. villataetal1998 (1998)) and the helical motion of knots (e.g. Abraham & Carrara abrahamcarrara1998 (1998); Britzen et al. britzenetal2001 (2001)).
One expects that a binary framework for AGN may be particularly relevant for our understanding of BL Lac objects such as Mkn 501 (e.g. De Paolis et al. depaolisetal2002 (2002)).
Mkn 501 ($`z=0.034`$) is one of at least four AGN which have been reliably detected at TeV energies (see Catanese & Weekes cataneseweekes99 (1999) for a review). From optical observations its host is known to be the elliptical galaxy UGC 10599 (cf. Stickel et al. stickeletal93 (1993)). As a BL Lac object, Mkn 501 belongs to the blazar class of AGN which are thought to have relativistic jets oriented at a small angle to the line of sight, so that the observed flux is strongly enhanced by relativistic beaming effects. Detailed observations since 1997 have confirmed the picture of powerful activity in Mkn 501 and revealed a wide range of variability patterns during outbursts (e.g. Protheroe et al. protheroeetal1998 (1998)) including flaring episodes of several days and rapid variability on subhour timescale (e.g. Sambruna et al sambrunaetal2000 (2000)). The evidences for a $`23`$-day periodicity in both the TeV and the X-ray lightcurves during the 1997 high state (cf. Hayashida et al. hayashidaetal1998 (1998); Kranich et al. kranichetal1999 (1999); Nishikawa et al. nishikawaetal1999 (1999)) may be counted among the most fascinating features (see Kranich kranichdiss2001 (2001) and Kranich et al. kranichetal2001 (2001) for a recent assessment of significance) and possibly indicates the presence of a BBHS in Mkn 501 (cf. Rieger & Mannheim riegermannheim2000 (2000)). In this paper we analyse the implications of the apparent disagreement of independent central mass estimates for Mkn 501. Starting with an investigation of the constraints given by host galaxy observations (Sect. 2), we proceed with an analysis of the robustness of the HE emission estimates (Sect. 3) and finally consider the consequences expected from a binary BH modelling (Sect. 4).

## 2 Central mass estimates derived from host galaxy observations

Recent mass estimates using host galaxy observations mainly rely on the following two statistical correlations:
(1) Dynamical studies of nearby elliptical galaxies have revealed an apparent, almost linear correlation (albeit with significant intrinsic scatter) between the central BH mass and the B-band luminosity of the bulge part of the host galaxy, which is given by $`M_H0.78\times 10^8M_{}(L_{\mathrm{B},\mathrm{bulge}}/10^{10}L_{})^{1.08}`$ (Magorrian et al. magorrianetal1998 (1998); Kormendy & Gebhardt KG2001 (2001), hereafter KG01).
(2) A much tighter correlation $`M_H\sigma ^\alpha `$ seems to exist between the BH mass $`M_H`$ in nearby inactive galaxies and the stellar velocity dispersion $`\sigma `$ of their host bulge (Gebhardt et al. gebhardtetal2000 (2000); Ferrarese & Merritt FM00 (2000)). However, up to now there is considerable debate over the true slope $`\alpha `$. Using different samples, Gebhardt et al. (gebhardtetal2000 (2000)) found $`\alpha =3.75\pm 0.3`$, while Ferrarese & Merritt (FM00 (2000)) obtained $`\alpha =4.8\pm 0.54`$ and $`\alpha =4.72\pm 0.36`$ (Merritt & Ferrarese 2001a ). Currently, further research is required to settle the question whether this difference is mostly caused by lower quality data and a less precise regression algorithm (cf. Merritt & Ferrarese 2001a ) or by systematic differences in the velocity dispersions used by the different groups for the same galaxies (cf. Tremaine et al. tremaineetal2002 (2002)).
The results using reverberation mapping (RM) (e.g. Kaspi et al. kaspietal2000 (2000); Nelson nelson2000 (2000); Wandel wandel2002 (2002)) indicate that the $`M_H\sigma `$ correlation may also hold for nearby AGN. However, a critical test of this conclusion depends on both a secure measure of the BH mass and an accurate determination of the stellar velocity dispersion. So far, the quality of BH mass estimates from stellar or gas kinematical data (whether ground- or HST-based), which require the BH sphere of influence $`r_H=GM_H/\sigma ^2`$ to be well-resolved, seems to increase only modestly (Merritt & Ferrarese 2001a ), so that over-estimation may be quite possible. Further progress has been expected using RM methods (Ferrarese et al. ferrareseetal2001 (2001)). Yet, the accuracy of RM may be strongly affected by systematic errors, e.g. due to uncertainties in the geometry and kinematics of the BLR or due to an unknown angular radiation pattern of the line emission, which may result in a systematic error up to at least a factor of $`3`$ (cf. Krolik krolik2001 (2001)). Moreover, only few accurate measurements of $`\sigma `$ seem to exist for AGN. Ferrarese et al. (ferrareseetal2001 (2001)) have recently analysed six AGN with well-determined RM BH masses by a careful measurement of their velocity dispersions and found a general consistency with the $`M_H\sigma `$ relation for quiescent galaxies. However, only BH masses below $`10^8M_{}`$ have been included so far, thus leaving out the high mass end of the correlation, and in addition, a large scatter is indicated. Besides providing a promising tool for the determination of BH masses in AGN, the current uncertainties in the correlations should be considered, if one tries to assess its implication for individual sources such as Mkn 501.

In the case of Mkn 501, Barth et al. (2002a ) have recently determined the stellar velocity dispersion from the calcium triplet lines to be $`\sigma =(372\pm 18)`$ km/s (cf. also Barth et al. 2002b ). Applying the $`M_H\sigma `$ relations of KG01 and Merritt & Ferrarese (2001a ), they derived a BH mass for Mkn 501 of $`M_H(0.93.4)\times 10^9M_{}`$. This mass estimate was supported by the study of Wu et al. (wuetal2002 (2002)), who estimated the velocity dispersions and BH masses from the fundamental plane for ellipticals for a large AGN sample including $`63`$ BL Lac objects (but not Mkn 501). They derived BH masses up to $`10^9M_{}`$, but with a potential error up to a factor of two. In particular, inspection of the fit in their Fig. 1 indicates a possible BH mass for Mkn 501 of $`(47)\times 10^8M_{}`$ for $`M_R(\mathrm{host})=24.2`$ mag (Pursimo et al. pursimoetal2002 (2002)). The general challenge of determining $`\sigma `$ accurately may be illustrated in more detail with reference to the recent work by Falomo et al. (falomoetal2002 (2002)), who provided a systematical study of the stellar velocity dispersion in seven BL Lacs. Using measurements in two spectral ranges, they found a velocity dispersion of $`\sigma =(291\pm 13)`$ km/s for Mkn 501, which is significantly lower than the one derived by Barth et al. (2002a ). Hence, if this value is used instead, the BH mass estimated by Barth et al. (2002a ) is reduced by up to a factor of three, i.e. one obtains $`M_H(3.610.7)\times 10^8M_{}`$. Additional support for such a low $`\sigma `$-value in Mkn 501 seems to be indicated by the original Faber & Jackson relation, which yields $`\sigma 270`$ km/s (see Fig. 2 in Falomo et al. falomoetal2002 (2002)) for $`M_R(\mathrm{host})=24.2`$ mag. Future research is needed to test whether the discrepancy in $`\sigma `$ is mainly induced by the difference in the method deriving $`\sigma `$ (direct fitting versus Fourier quotient routine).

As noted above, an additional mass estimate for Mkn 501 can also be derived from the $`M_HL_{\mathrm{bulge}}`$ correlation. The reported large uncertainties in this relation have recently been examined by McLure & Dunlop (mclure2002 (2002)) using R-band luminosities, which are less sensitive to extinction. By analysing the virial BH masses for a sample of $`72`$ AGN, they found the scatter to be quite smaller than previously estimated and stressed its usefulness. For application to Mkn 501, we may exploit the absolute R-band luminosity of its host galaxy recently derived by Pursimo et al. (pursimoetal2002 (2002)) (see also Nilsson et al. nilssonetal1999 (1999)). Assuming $`H_0=50`$ km s-1 Mpc-1, they obtained $`M_R(\mathrm{host})=24.2`$ mag. If we convert R- to B-band luminosity assuming $`BR=1.56`$ (e.g. Goudfrooij et al. goudfrooijetal1994 (1994); Fukugita et al. fukugitaetal1995 (1995); Urry et al. urryetal2000 (2000)), we have $`M_B(\mathrm{host})=22.46`$ mag, which results in $`L(M_B(\mathrm{host}))=1.406\times 10^{11}L_{}`$. Using the KG01-relation for the B-band luminosity, the expected BH mass in Mkn 501 is $`M_H1.3\times 10^9M_{}`$, but with a potential error of up to at least a factor of three. Using the more recent McLure & Dunlop (mclure2002 (2002))-relation $`\mathrm{log}(M_H/M_{})=0.5M_R2.96(\pm 0.48)`$, one finds $`M_H(0.464.2)\times 10^9M_{}`$. Uncertainties in the determination of $`M_R`$ may further reduce the expected BH mass. For example, values from the literature presented in Table 4 of Nilsson et al. (nilssonetal1999 (1999)) indicate that $`M_R`$ might be up to $`0.4`$ mag higher and therefore $`M_H`$ correspondingly smaller. More importantly, if a Hubble constant $`H_0=70`$ km s-1Mpc-1 is assumed, one finds $`M_R=23.47`$, which results in $`M_H7.8\times 10^8M_{}`$ (KG01), again with substantial scatter of up to at least a factor of three, or $`M_H(218)\times 10^8M_{}`$ (McLure & Dunlop mclure2002 (2002)), thus allowing for a central mass as low as $`2\times 10^8M_{}`$.

## 3 Mass estimates derived in the context of high energy emission models

Following a quite different approach, the BH mass of Mkn 501 could also be estimated in the framework of high energy emission models:
(1) With respect to the high energy emission, Fan et al. (fan99 (1999)) have recently determined the central black hole masses for several $`\gamma `$-ray loud blazars (including Mkn 501) by assuming that the observed $`\gamma `$-rays are produced at $`200R_g`$. Accordingly, the central BH mass obeys the relation $`M_H/M_{}500\delta \mathrm{\Delta }t`$, where $`\mathrm{\Delta }t`$ denotes the doubling timescale in seconds and $`\delta `$ the Doppler factor. For Mkn 501 they obtained a central BH mass of $`M_H0.9\times 10^7M_{}`$, using $`\mathrm{\Delta }t=6`$ hrs and $`\delta 0.9`$. However, current evidence indicates that the relevant timescale $`\mathrm{\Delta }t`$ might be substantially smaller. Recent observations reveal considerable sub-hour variability (e.g. Ghosh et al. ghoshetal2000 (2000); Catanese & Sambruna catanesesambruna (2000)) on a timescale as low as $`\mathrm{\Delta }t=1200`$ s (Sambruna et al. sambrunaetal2000 (2000)). Such a low value for the observed timescale could possibly be accommodated by assuming a high doppler boosting factor $`\delta 18`$. Hence, it appears that the more crucial point in this derivation is the assumption that the $`\gamma `$-rays dominating the emission are produced at $`200`$ Schwarzschild radii. Indeed, at least in the case of the TeV-blazars, the variable, high energy emission is usually regarded as produced by moving knots or shocks in the jet far from the accretion disk (for a review, cf. Mannheim mannheim97 (1997); Aharonian & Völk ahavoelk2001 (2001) and references therein). While instructive, the derived estimate should thus not be considered as a robust, general limit.
(2) A further mass estimate for Mkn 501 has been derived by DeJager et al. (dejager1999 (1999)) following an approach developed by Hayashida & Miyamoto et al. (hayashidaetal1998k (1998)). Assuming the variation in the accretion process to drive the X-ray and TeV variation in the jet via the dynamo effect, their result yields a central BH mass of $`M_H=(16)\times 10^7M_{}`$ for $`\delta =(1015)`$. However, due to the absence of a physical basis for the required scaling of the Fourier spectrum and due to the assumption of a linear proportionality between variability timescale and BH mass (cf. Hayashida & Miyamoto et al. hayashidaetal1998k (1998)), which is probably not valid for the blazar class (Kataoka et al. kataokaetal2001 (2001)), this estimate again does not appear to be robust.
(3) The estimates (1) and (2) which suggest a BH mass less than $`6\times 10^7M_{}`$, are strongly model-dependent as shown above. We may illustrate this in more detail by comparing them with results derived in the context of another, high energy emission model. Bednarek, Kirk & Mastichiadis (bednareketal1996 (1996)) for example, have developed a special model for the origin of the high energy particles in TeV blazars like Mkn 421 and Mkn 501, assuming the electrons responsible for the high energy emission to be accelerated rectilinearly in an electric field. In this model, the mass of the central BH is expected to be limited by $`M_H\stackrel{>}{}10^8M_{}(E_{\mathrm{max}}/34\mathrm{TeV})^{2.5}/l_{\mathrm{Edd}}^{1/2}`$, where $`E_{\mathrm{max}}`$ denotes the maximum photon energy and $`l_{\mathrm{Edd}}`$ the disk luminosity in units of the Eddington luminosity. There is strong evidence for a sub-Eddington accretion mode in BL Lacs in general (e.g. Cavaliere & D’Elia cavalieredelia2002 (2002)) and particularly for the TeV emitting blazars (e.g. Celotti et al. celottietal1998 (1998)). Thus, using characteristic values, i.e. $`l_{\mathrm{Edd}}=(0.010.001)`$ and $`E_{\mathrm{max}}=20`$ TeV (cf. Samuelson et al. samuelsonetal98 (1998); Konopelko konopelko99 (1999)), we arrive at a mass $`M_H\stackrel{>}{}(2.4614.66)\times 10^8M_{}`$, which is up to ten times larger than the estimates (1) and (2).

## 4 Mass estimates derived in the context of a binary black hole scenario

In a recent contribution (Rieger & Mannheim riegermannheim2000 (2000), hereafter RMI; riegermannheim2001 (2001)), we have shown that the periodicity of $`P_{\mathrm{obs}}=23`$ days, observed during the 1997 high state of Mkn 501, could be plausibly related to the orbital motion in a BBHS, provided the jet, which dominates the observed emission, emerges from the less massive (secondary) BH. If such an interpretation (henceforth called the standard scenario) is appropriate, we may derive a third estimate for the central mass in Mkn 501. We have demonstrated in RMI that, due to relativistic effects, the observed period appears drastically shortened, so that for the intrinsic keplerian orbital period one finds $`P_k=(614)\mathrm{yrs}`$. Taking into account that the observed emission is periodically modulated by differential doppler boosting due to the orbital motion, one may derive a simple equation for the required mass dependence in the standard scenario (cf. RMI, Eq. (8)):

$`{\displaystyle \frac{M}{(m+M)^{2/3}}}`$ $`=`$ $`{\displaystyle \frac{P_{\mathrm{obs}}^{1/3}}{(2\pi [1+z]G)^{1/3}}}{\displaystyle \frac{c}{\mathrm{sin}i}}`$ (1)
$`\times {\displaystyle \frac{f^{1/(3+\alpha )}1}{f^{1/(3+\alpha )}+1}}(1{\displaystyle \frac{v_z}{c}}\mathrm{cos}i)^{2/3}.`$
Here $`\alpha `$ denotes the spectral flux index, $`z`$ the redshift, $`v_z`$ the outflow velocity, $`i1/\mathrm{\Gamma }_b`$ the angle of the jet axis to the line of sight, $`\mathrm{\Gamma }_b1015`$ (e.g. Mannheim et al. mannheimetal1996 (1996); Spada et al. spadaetal1999 (1999)) the bulk Lorentz factor and $`f`$ the observed flux ratio between maximum and minimum (for the TeV range $`f8`$), while $`m`$ and $`M`$ denote the masses of the smaller and larger BH, respectively.
In order to break the degeneracy in this mass ratio, we may utilize an additional constraint by assuming that the current binary separation $`d`$ corresponds to the separation at which gravitational radiation becomes dominant (cf. RMI). Such a constraint yields an upper limit for the allowed binary masses and might be associated with the key aspect that BL Lac objects are old, more evolved and underluminous sources, i.e. they might be close binaries, probably settled above or near the critical gravitational separation, because the possibility of removing further angular momentum has been almost terminated as a result of declining gas accretion rates. We can specify the corresponding gravitational separation $`d_g`$ by equating the timescale $`\tau _{\mathrm{grav}}=|d/\dot{d}|=(5c^5/64G^3)d^4/(Mm[m+M])`$ on which gravitational radiation shrinks the binary orbit, with the dynamical timescale $`\tau _{\mathrm{gas}}`$ for gas accretion (cf. Begelman et al. begeletal80 (1980); note that compared with RMI, this estimate for $`\tau _{grav}`$ is a factor $`2.5`$ more precise, cf. Rieger & Mannheim riegermannheim2001 (2001)). A characteristic measure for $`\tau _{\mathrm{gas}}`$ is given by the Eddington limit $`\tau _{\mathrm{gas}}3.77\times 10^7(\eta /0.1)`$ yrs, assuming a canonical $`10\%`$ conversion efficiency (cf. Krolik krolik99 (1999)). Using $`\tau _{\mathrm{grav}}`$ and $`\tau _{\mathrm{gas}}`$, one finally arrives at

$$d_g3.50\times 10^{16}M_8^{1/4}m_8^{1/4}(m_8+M_8)^{1/4}\mathrm{cm}.$$
(2)
If we combine this expression with the relevant expression for the binary separation (i.e. Eq. (7) of RMI), we obtain a second constraint on the allowed mass ratio given by

$`{\displaystyle \frac{(m_8+M_8)^{3/4}}{M_8^{5/4}m_8^{1/4}}}=7.3\times 10^6{\displaystyle \frac{f^{1/(3+\alpha )}+1}{f^{1/(3+\alpha )}1}}{\displaystyle \frac{(1+z)}{P_{\mathrm{obs}}}}\mathrm{sin}i,`$ (3)
assuming the jet to arise from the less massive BH (for the reverse case the BH masses should be interchanged in Eq. (1) and Eq. (3)). By using Eq. (1) and Eq. (3), we may determine the binary masses permitted for a BBHS with separation near that for which gravitational radiation becomes dominant. In Fig. 1 we have plotted the allowed range for the case of $`\mathrm{\Gamma }_b=15`$ and $`\alpha (1.21.7)`$ (e.g. Aharonian et al. aharonianetal1999 (1999)). The run of the curves is determined by Eq. (1) with the upper bounds given by Eq. (3).

We note that the upper bounds depend on the presumed timescale for the binary evolution, i.e. on the estimate of the gravitational separation. In particular, the derived upper bounds may be larger, if, for example, disk-driven migration (e.g. Armitage & Natarajan armitage2002 (2002)) might be relevant. Fig. 1 reveals that (a) the typical combined mass for a BBHS is expected to be less than $`3\times 10^8M_{}`$, and that (b) a BBHS with, e.g. $`m6\times 10^7M_{}`$ and combined central mass $`(M+m)(1.52)\times 10^8M_{}`$ appears well-conceivable given the current limits.
On the other hand, if a high central mass of $`10^9M_{}`$ will be established by further research, the proposed binary scenario appears to be ruled out. We note, however, that even in this case a binary scenario may be still possible provided that the jet, which dominates the emission, is produced by the primary BH. To illustrate the implications in this case, let us consider a (combined) central mass of $`0.9\times \times 10^9M_{}`$ by demanding the primary to be in the range $`(56)\times 10^8M_{}`$ (see Table 1). The mass of the secondary BH then is determined by Eq. (1) with the masses interchanged, the current separation $`d`$ by Eq. (3) of RMI (again with the masses interchanged) and the gravitational separation $`d_g`$ by Eq. (2). The results as shown in Table 1 imply a close binary system with $`d\stackrel{<}{}d_g`$. If the optically bright QSO stage thus occurs during the binary evolution and the applied doppler factors are considered as typical, phases of super-Eddington accretion and/or with decreased conversion efficiency seem to be necessary for the binary to be above its gravitational separation.

## 5 Discussion and conclusions

In this paper, we have considered the apparent discrepancies between independent central mass estimates for the paradigmatic object Mkn 501. Our main conclusions can be summarized as follows:
(1) Estimates inferred from host galaxy observations using the reported, statistical correlations for the central velocity dispersion and the bulge luminosity reveal a tendency for the central BH mass in Mkn 501 to be larger than $`5\times 10^8M_{}`$. If verified by further investigations, this would exclude both the standard binary scenario, where the jet is emitted from the less massive BH, and several high energy emission models proposed so far. However, as the correlations still show substantial intrinsic scatter and large uncertainties due to the use of limited samples – especially with respect to the high mass end – the accuracy of this estimate should be treated with caution. In particular, we have shown that the central mass limit could be easily up to a factor of $`(23)`$ smaller, thus allowing for the standard binary scenario, in which case a simple explanation for the divergent mass estimates from emission models and host observations appears possible.
(2) Independent mass estimates from high energy emission models suggest a mass for the jet emitting BH in Mkn 501 which is smaller than $`6\times 10^7M_{}`$. However, this estimate is quite model-dependent and so cannot be used as a universally valid constraint. In particular, emission models have been developed for Mkn 501 where a central BH mass up to ten times larger appears quite possible.
(3) If a binary BH interpretation for the periodical variability (cf. RMI) is appropriate in the case of Mkn 501, the jet dominating the emission has to be produced by the less massive (secondary) BH. Using characteristic jet parameters and assuming both Eddington-limited accretion and a binary separation comparable to the gravitational one, we have shown that the combined central (primary $`+`$ secondary BH) mass should be smaller than $`3\times 10^8M_{}`$. The binary model particularly permits a system with $`m=6\times 10^7M_{}`$ and $`(m+M)(1.52)\times 10^8M_{}`$, possibly allowing for a convergence of constraints from emission models and host observations.
(4) If further research clearly establishes a high central BH mass $`10^9M_{}`$ in Mkn 501, this would call for a change of fundamental assumptions and/or a modification of the parameter space regarded to be typical in the context of several high energy emission models. The only way to avoid this conclusion seems to be the assumption of a highly unequal BBHS with the jet dominating the emission produced by the less massive and the central mass dominated by the primary BH. However, it then appears no longer be possible to explain the observed periodicity via the orbital motion. On the other hand, even in the case of a high central mass, a binary model for the observed periodicity may be still possible provided the jet is produced by the primary BH. Yet, for the binary to be near or above the separation at which gravitational radiation becomes dominant, one then requires the earlier binary evolution to pass through phases of super-Eddington accretion and/or decreased conversion efficiency $`\eta <0.1`$. Such conditions need not necessarily be ad hoc but have already been considered in the context of galaxy evolution and the growth of massive BH (e.g. Collin et al. collinetal2002 (2002); Yu & Tremaine yutremaine2002 (2002)).
In view of Mkn 501 as a paradigmatic object, bringing its mass estimates to convergence continues to be an important task. Further research on the parameter space of emission models and the careful observations of nearby galaxies will be particularly valuable in order to minimise the uncertainties in the statistical correlations, to evaluate the impact of host observations on high energy emission models, and to assess the plausibility of a BBH model.

###### Acknowledgements.

We are grateful to J. Heidt for fruitful discussions about the host of Mkn 501, and C. Hettlage and R. Hessman for a helpful reading of the manuscript. F.M.R. gratefully acknowledges support under DFG MA 1545/8-1.