format: text

text:
# The Magnetohydrodynamics of Convection-Dominated Accretion Flows

## 1 Introduction

Models of radiatively inefficient accretion flows provide a useful framework for interpreting observations of low-luminosity black hole X-ray binaries and active galactic nuclei (see, e.g., Narayan, Mahadevan & Quataert 1998; Quataert 2001; Narayan 2002 for reviews). In the past few years, there has been rapid theoretical progress in understanding the dynamics of such flows. Much of this advance has been driven by numerical simulations. In particular, a number of hydrodynamic simulations have been reported in which angular momentum transport is put in “by hand” using an $`\alpha `$ prescription (e.g., two-dimensional simulations by Igumenshchev, Chen, & Abramowicz 1996; Igumenshchev & Abramowicz 1999, 2000; Stone, Pringle, & Begelman 1999; and three-dimensional simulations by Igumenshchev, Abramowicz & Narayan 2000). Such simulations reveal a flow structure very different from advection-dominated accretion flow models (ADAFs) that were proposed to describe the structure of radiatively inefficient accretion flows (Ichimaru 1977; Narayan & Yi 1994; Abramowicz et al. 1995). In an ADAF, the gas accretes rapidly, and the gravitational potential energy of the accreting gas is stored as thermal energy and advected into the central black hole. By contrast, in the hydrodynamic simulations, the rate of mass accretion is much smaller and strong radial convection efficiently transports energy outwards. The simulations have been interpreted in terms of a “convection-dominated accretion flow” model (CDAF; Narayan, Igumenshchev & Abramowicz 2000; Quataert & Gruzinov 2000; Abramowicz et al. 2002). In a CDAF, convection simultaneously transports energy outwards and angular momentum inwards, strongly suppressing the accretion rate onto the central black hole.

The relevance of the hydrodynamic simulations, and the CDAF model derived from them, is unclear. This is primarily because of the ad hoc treatment of angular momentum transport. It is believed that angular momentum transport in accretion flows is primarily due to MHD turbulence initiated by the magnetorotational instability (MRI; Balbus & Hawley 1991; BH91). Balbus & Hawley have argued that, because of the fundamental role played by magnetic fields, hydrodynamic models typically cannot describe the structure of the accretion flow (or differentially rotating systems more generally; e.g., Balbus & Hawley 1998; Balbus 2000; 2001). They have applied this criticism in detail to the CDAF model (Hawley, Balbus, & Stone 2001; Balbus & Hawley 2002; hereafter BH02). In an independent argument, BH02 also suggest that CDAF models violate the second law of thermodynamics. We do not consider this issue here, but intend to deal with it in a separate paper.

In this paper we discuss the physics of the CDAF model within the framework of MHD. Apart from clarifying the theoretical underpinnings of the CDAF model, we believe that our analysis also provides a useful illustration of the conditions under which a hydrodynamic model is applicable to an intrinsically magnetohydrodynamic situation; this is important in other astrophysical contexts, e.g., the excitation of density waves at Lindblad resonances in a magnetized accretion flow, or diskoseismology models of QPOs.

The paper is organized as follows. In §2 we discuss the linear stability of a differentially rotating and thermally stratified plasma in the presence of a weak vertical magnetic field. We identify and explain the difference between “convective” modes and “MRI” modes and show that the former transport angular momentum inwards while the latter usually transport angular momentum outwards. Then, in §2.1 we generalize this result to an arbitrary field orientation and arbitrary stratification. In §3 we discuss the implications of our results. In particular, we use the linear analysis to speculate about the conditions under which the saturated nonlinear turbulence in an accretion flow might give rise to a CDAF-like structure. In §4 we summarize our results.

## 2 Convective Modes in a Magnetized Differentially Rotating Plasma

We consider the linear stability of a differentially rotating and thermally stratified plasma in the presence of a weak magnetic field. We specialize to the case of a vertical magnetic field, $`B_zB`$, and modes with purely vertical wavevectors, $`k_zk`$. We also take the rotation rate, pressure, and density to be constant on cylinders, i.e., $`\mathrm{\Omega }(R)`$, $`P(R)`$, and $`\rho (R)`$. We thus employ the same equations as those in BH02 (which are a simplification of the more general analysis in BH91). Although the problem we analyze is identical to that considered by BH02, we differ in the interpretation of the results. In particular, we identify a set of modes that have all the features of growing convective modes in hydrodynamics, including inward transport of angular momentum.

The thermal stratification of the flow can be described by the Brunt-Väisälä frequency,

$$N^2=N_R^2=\frac{3}{5\rho }\frac{P}{R}\frac{\mathrm{ln}P\rho ^{5/3}}{R}.$$
(1)
In the absence of rotation, the system is convectively unstable if $`N^2>0`$. In a rotating medium, but ignoring magnetic fields, convection is present only if $`N^2>\kappa ^2`$, where $`\kappa ^2=4\mathrm{\Omega }^2+d\mathrm{\Omega }^2/d\mathrm{ln}R`$ is the epicyclic frequency.111The requirements for convection in an unmagnetized rotating medium are the Høiland criteria (e.g., Tassoul 1978). These reduce to the condition $`N^2>\kappa ^2`$ for the case considered here. For simplicity, we will refer to this condition as the Høiland criterion in this section. §2.1 treats the general case. Note that $`N^2>\kappa ^2`$ requires a sound speed comparable to the rotation speed.

For a weakly magnetized medium, BH91 derived the behavior of linear perturbations including the effects of thermal stratification. For fluid displacements of the form $`\xi \mathrm{exp}[ikz+\gamma t]`$, and for the particular geometry considered here, the growth rate $`\gamma `$ is given by

$$\gamma ^2=(kv_A)^2+\frac{1}{2}\left[N^2\kappa ^2\pm \sqrt{(N^2\kappa ^2)^2+16\mathrm{\Omega }^2(kv_A)^2}\right],$$
(2)
where $`v_A=B/\sqrt{4\pi \rho }`$ is the Alfvén speed. For the remainder of this section we focus on the growing mode in equation (2); this has $`\gamma ^2>0`$ and corresponds to the positive sign of the square root term. In the Appendix we present a more general analysis of both growing and oscillatory modes that helps clarify some of the issues discussed in this section.

To understand the physics behind the dispersion relation in equation (2) consider the equations for the radial and azimuthal displacements of fluid elements (BH91):

$$\frac{^2\xi _R}{t^2}2\mathrm{\Omega }\frac{\xi _\varphi }{t}=\left(\frac{d\mathrm{\Omega }^2}{d\mathrm{ln}R}+(kv_A)^2\right)\xi _R+N^2\xi _R,$$
(3)
$$\frac{^2\xi _\varphi }{t^2}+2\mathrm{\Omega }\frac{\xi _R}{t}=(kv_A)^2\xi _\varphi .$$
(4)
Solving equation (4) for $`\xi _\varphi `$ we can substitute into equation (3) to find

$$\frac{^2\xi _R}{t^2}=a_H+a_{M,s}+a_{M,d},$$
(5)
where the three acceleration terms are of the form

$$a_H=(N^2\kappa ^2)\xi _R,a_{M,s}=(kv_A)^2\xi _R,a_{M,d}=\left[\frac{4\mathrm{\Omega }^2k^2v_A^2}{\gamma ^2+k^2v_A^2}\right]\xi _R.$$
(6)
Since $`^2\xi _R/^2t=\gamma ^2\xi _R`$ it is straightforward to solve equation (5) and obtain the dispersion relation in equation (2).

The first acceleration term in equation (5), $`a_H`$, is a pure “hydrodynamic” or “Høiland” term (hence the subscript $`H`$) which is present even in the absence of magnetic fields. This term can be either stabilizing or destabilizing; it is stabilizing when $`N^2<\kappa ^2`$ (the medium is convectively stable by the Høiland criteria) and destabilizing when $`N^2>\kappa ^2`$ (convectively unstable). The second and third terms in equation (6) owe their existence to the magnetic field (hence subscript $`M`$). One of them is always stabilizing (subscript $`s`$) and the other is always destabilizing (subscript $`d`$).

Equations (3)–(6) have a very clean physical interpretation, as we explain in the rest of this section. For the simple geometry considered here (vertical field and vertical wavevector), gas and magnetic pressure forces are unimportant. If $`N^2=0`$, the only forces in the problem are magnetic tension and rotating frame dynamics (e.g. the Coriolis force). Equations (3) and (4) then describe the behavior of fluid elements that are coupled by the tension of magnetic field lines. As Balbus & Hawley (1992, 1998) have discussed, these equations are identical to those describing fluid elements coupled by a spring of spring constant $`(kv_A)^2`$. For a stratified medium with $`N^20`$, the radial acceleration of a fluid element has an additional component due to the radial buoyancy force. This is the $`N^2\xi _R`$ term in equations (3) and (5). When $`N^2<\kappa ^2`$, this term does not have a large effect and introduces modest changes to the values of numerical factors. However, when $`N^2>\kappa ^2`$, which corresponds to a convectively unstable system by the Høiland criterion, the acceleration term $`a_H`$ becomes positive and destabilizing. This introduces qualitatively new effects.

The angular momentum flux of a linear mode is given by $`R\mathrm{\Omega }T_{R\varphi }`$, where $`T_{R\varphi }`$ is the $`R\varphi `$ component of the fluid stress tensor,

$$T_{R\varphi }=\rho (\delta v_R\delta v_\varphi \delta v_{AR}\delta v_{A\varphi }),$$
(7)
$`\delta 𝐯`$ is the perturbed fluid velocity, and $`\delta 𝐯_𝐀\delta 𝐁/\sqrt{4\pi \rho }`$. Using the eigenfunctions in BH91 we confirm BH02’s expression for the stress tensor, namely

$$\frac{T_{R\varphi }}{\rho |\delta v_R^2|}=\frac{\mathrm{\Omega }}{D\gamma }\left[\frac{k^2v_A^2}{\gamma ^2}\left(\left|\frac{d\mathrm{ln}\mathrm{\Omega }}{d\mathrm{ln}R}\right|+2\right)\frac{\kappa ^2}{2\mathrm{\Omega }^2}\right]\left(\frac{\mathrm{\Omega }}{\gamma }\right)t_{R\varphi },$$
(8)
where $`D=1+k^2v_A^2/\gamma ^2`$ and $`t_{R\varphi }`$ is a useful dimensionless stress whose sign gives the sign of $`T_{R\varphi }`$.

We are now in a position to elucidate the physics of the linear instability calculation. Figure 1a shows a plot of the growth rate $`\gamma /\mathrm{\Omega }`$ as a function of the dimensionless wavevector $`kv_A/\mathrm{\Omega }`$ for 5 values of the Brunt-Väisälä frequency: $`N^2/\mathrm{\Omega }^2=0,0.5,1,1.5,`$ and $`2`$. We have assumed a Keplerian rotation profile, for which $`\kappa ^2=\mathrm{\Omega }^2`$. Thus, the first two choices of $`N^2`$ correspond to a stable entropy gradient by the Høiland criterion, the third is neutrally stable, and the last two choices correspond to an unstable entropy gradient. Figure 1b shows the dimensionless stress $`t_{R\varphi }`$ as a function of $`kv_A/\mathrm{\Omega }`$ for the same five values of $`N^2`$. It is useful to note that the self-similar ADAF model gives $`N^2/\mathrm{\Omega }^2=15(\gamma _{ad}1)/4\gamma _{ad}`$, where $`\gamma _{ad}`$ is the adiabatic index of the gas (Narayan & Yi 1994).222This expression is not valid if $`\gamma _{ad}=5/3`$ since then $`\mathrm{\Omega }^2=0`$ in the analytical models. For example, if $`\gamma _{ad}=1.5`$, $`N^2/\mathrm{\Omega }^2=1.25`$.

For $`N^2=0`$, the results reproduce the usual MRI. The growth rate is very small for $`kv_A\mathrm{\Omega }`$ and peaks when $`kv_A`$ is $`\mathrm{\Omega }`$, with a peak value of $`\gamma =3\mathrm{\Omega }/4`$; the angular momentum flux is always outwards ($`T_{R\varphi }>0`$). The instability is clearly triggered by the magnetic field via the destabilizing term $`a_{M,d}`$, since the other two terms in equation (5) are both stabilizing (negative). The MRI survives with qualitatively the same behavior even for non-zero $`N^2`$, so long as the Høiland criterion for convective stability ($`N^2<\kappa ^2`$) is satisfied.

Consider, however, a flow with $`N^2>\kappa ^2`$. Such a flow is unstable to convection according to hydrodynamics, and is therefore the case of interest for CDAFs. Figure 1 shows that, for $`N^2>\kappa ^2`$, modes with $`kv_A\mathrm{\Omega }`$ are not very different from their $`N^2=0`$ counterparts. The angular momentum flux is still outwards and the growth rate of the mode is comparable (though somewhat larger). By contrast, long-wavelength modes are very different. As Figure 1 shows, long-wavelength modes are strongly unstable and they transport angular momentum inwards. In fact, as can be seen by inspection, the $`kv_A\mathrm{\Omega }`$ limit of the dispersion relation in equation (2) (with the + root) gives $`\gamma ^2=N^2\kappa ^2`$, and the corresponding stress tensor is $`T_{R\varphi }/(\rho |\delta v_R|^2)=\kappa ^2/[2\mathrm{\Omega }(N^2\kappa ^2)^{1/2}]<0`$. The long-wavelength modes are clearly independent of the magnetic field and behave very differently from MRI modes.

The physics of the $`kv_A\mathrm{\Omega }`$ modes is as follows. In the pure MRI problem, namely $`N^2=0`$ (or more generally $`N^2<\kappa ^2`$), long-wavelength modes are only weakly unstable because the magnetic tension forces that are central to the MRI are weak; the field lines are hardly bent by a long-wavelength perturbation. By contrast, when $`N^2>\kappa ^2`$, the fluid has two destabilizing forces, the buoyancy term ($`a_H`$) due to an unstable entropy gradient and the standard MRI term ($`a_{M,d}`$). The buoyancy force that drives convection is independent of the wavelength of the perturbation, whereas the MRI depends on $`k`$. For long wavelengths, buoyancy is much more important and controls all the physics of the mode, both the growth rate and the angular momentum transport. Moreover, we show in the Appendix that the “convective” mode driven by buoyancy is the only unstable long-wavelength mode in MHD for a Høiland-unstable medium ($`N^2>\kappa ^2`$). In particular, long-wavelength perturbations for which magnetic tension dominates (referred to as $``$ modes in the Appendix) are stable oscillatory waves rather than unstable modes (see eq. A3 and Fig. 3).

Figure 2 explicitly compares the relative importance of the two destabilizing terms, the buoyancy term $`a_H`$ and the magnetic tension term $`a_{M,d}`$. Let us identify a mode as an “MRI” mode if aH
<aM,d
<subscript𝑎𝐻subscript𝑎𝑀𝑑a_{H}\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{$\m@th\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr} }}}a_{M,d} and as a “convective” mode if $`a_H>a_{M,d}`$. With this natural identification, we see that, for $`N^2<\kappa ^2`$, all unstable modes are MRI modes, regardless of the wavelength. Even when $`N^2>\kappa ^2`$, short-wavelength modes with kvA
>Ω
>𝑘subscript𝑣𝐴Ωkv_{A}\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{$\m@th\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}}\Omega are still MRI modes. However, unstable long-wavelength modes with $`kv_A\mathrm{\Omega }`$ are clearly convective modes. The growth of these modes is due entirely to the unstable buoyancy force; correspondingly, the modes transport angular momentum inwards (Fig. 1b), just as indicated by hydrodynamic studies (e.g., Ryu & Goodman 1992). Indeed these modes are virtually indistinguishable from their hydrodynamic counterparts. All of the properties of hydrodynamic convective instabilities therefore carry over to a magnetohydrodynamic analysis.

### 2.1 The Hydrodynamic Limit

Magnetic tension should be weak for fluctuations with $`kv_A\mathrm{\Omega }`$ regardless of the details of the magnetic field geometry and thermal stratification. In this subsection, we show this explicitly by considering the long-wavelength limit of the most general axisymmetric MHD dispersion relation.

Balbus (1995) considered the linear, adiabatic, and axisymmetric stability of a differentially rotating weakly magnetized plasma in the presence of both vertical and radial stratification. He showed that the dispersion relation is given by

$$\stackrel{~}{\omega }^4\frac{k^2}{k_z^2}\stackrel{~}{\omega }^2\left[\frac{3}{5\rho }(DP)(D\mathrm{ln}[P\rho ^{5/3}])+\frac{1}{R^3}D(R^4\mathrm{\Omega }^2)\right]4\mathrm{\Omega }^2(𝐤𝐯_𝐀)^2=0,$$
(9)
where

$$D=\left(\frac{k_R}{k_z}\frac{}{z}\frac{}{R}\right),$$
(10)
and $`\stackrel{~}{\omega }^2=\gamma ^2(𝐤𝐯_𝐀)^2`$ (for fluctuations $`\mathrm{exp}[\gamma t+i𝐤𝐫]`$). Unlike in the previous section, $`\mathrm{\Omega }`$, $`P`$, and $`\rho `$ are now allowed to be functions of both $`z`$ and $`R`$, the magnetic field has an arbitrary direction (i.e., $`B_R`$, $`B_\varphi `$, and $`B_z`$ components), and the wavevector is given by $`𝐤=k_R\widehat{𝐑}+k_z\widehat{𝐳}`$.

For $`𝐤𝐯_𝐀=0`$ (but $`k/k_z`$, $`k_R/k_z`$, etc., finite), equation (9) reduces to

$$\gamma ^2\frac{k^2}{k_z^2}+\left[\frac{3}{5\rho }(DP)(D\mathrm{ln}[P\rho ^{5/3}])+\frac{1}{R^3}D(R^4\mathrm{\Omega }^2)\right]=0.$$
(11)
Equation (11) is independent of the magnetic field. Indeed, as noted by Balbus (1995), it describes the linear, adiabatic, axisymmetric, hydrodynamic stability of a differentially rotating and thermally stratified medium (e.g., Goldreich & Shubert 1967). As is readily confirmed, equation (11) has unstable solutions if and only if the Høiland criteria are satisfied (e.g., Tassoul 1978).

The hydrodynamic limit in equation (11) formally requires setting $`𝐤𝐯_𝐀=0`$ in the MHD dispersion relation (eq. ). If $`𝐤𝐯_𝐀0`$, the full dispersion relation (eq. ), which has twice the number of modes (because it is fourth order in $`\gamma `$ rather than second order), must be used. However, if the medium is unstable by the Høiland criteria, and if we focus on long-wavelength MHD instabilities with $`𝐤𝐯_𝐀\mathrm{\Omega }`$, then the MHD system will behave like its hydrodynamic counterpart; in particular, the only unstable mode is effectively hydrodynamical and its growth rate is given by equation (11) with small corrections $`(𝐤𝐯_𝐀/\mathrm{\Omega })^21`$. Note that equation (11) describes the unstable growing mode only when the medium is Høiland unstable. If the medium is Høiland stable, the only possible instabilities are intrinsically MHD in nature and equation (9) must be used even for $`𝐤𝐯_𝐀\mathrm{\Omega }`$. This point is explained in more detail in the Appendix.

The analysis in this section shows that all of the properties of axisymmetric hydrodynamic instabilities carry over to MHD in the long-wavelength limit. In particular, the results of the previous section are general, and are not an artifact of the simplifying assumptions made there. Although our quantitative analysis is restricted to axisymmetry, we suspect that long-wavelength non-axisymmetric modes will behave hydrodynamically as well (since they also do not significantly bend the magnetic field lines).

## 3 Implications for CDAFs

To assess the implications of our linear analysis it is useful to consider wavelengths relative to the scale height of the disk, $`Hc_s/\mathrm{\Omega }`$, where $`c_s`$ is the sound speed. The maximally growing MRI mode has $`kH\beta ^{1/2}`$, where $`\beta c_s^2/v_A^2`$ is the ratio of the thermal energy to the magnetic energy in the disk. By contrast, the longer wavelength convective modes have $`kH\beta ^{1/2}`$; more quantitatively, the linear analysis suggests that buoyancy forces become important for kH
<0.3β1/2
<𝑘𝐻0.3superscript𝛽12kH\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{$\m@th\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr} }}}0.3\beta^{1/2} (Fig. 2). For these modes to be of interest they must fit in the disk and so must have kH
>1
>𝑘𝐻1kH\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{$\m@th\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}}1. This in turn requires $`\beta 1`$, i.e., that the magnetic field must be sub-equipartition (β
>10
>𝛽10\beta\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{$\m@th\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}}10 may suffice).

Given the above analysis, we propose the following identification between the hydrodynamic CDAF model and the present MHD results. If the magnetic field saturates at a sub-equipartition value, the most unstable MRI mode operates on a scale $`\beta ^{1/2}HH`$ and transports angular momentum outwards. Moreover, if the medium is Høiland unstable, only such short-wavelength fluctuations transport angular momentum outwards. Longer wavelength convective fluctuations are dominated by scales $`H`$ and transport angular momentum inwards (Fig. 1).

Which fluctuations are more important? Or, more precisely, which part of the power spectrum of fluctuations is more important, those with $`kH1`$ or those with $`kH\beta ^{1/2}`$? This cannot really be addressed analytically because we do not understand the saturation of the linear instabilities (and it is clear that they are coupled; e.g., convection, by stirring the fluid, can help to build up the magnetic field). Nonetheless, it is typically the case that the energy in a turbulent plasma is dominated by the largest unstable scales in the medium. This suggests that, if $`\beta 1`$ (perhaps β
>10
>𝛽10\beta\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{$\m@th\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}}10 is sufficient), the larger-scale convective fluctuations will dominate the dynamics, as proposed in the hydrodynamic models. Moreover, the inward angular momentum transport by convective modes depends on the magnitude of $`N^2\kappa ^2`$ (Fig. 1b). Thus there is a natural way for the inward convective transport to self-adjust to counteract outward transport by small scale MRI fluctuations. In the limit of efficient convection and $`\beta 1`$ we would expect marginal stability to convection according to the Høiland criteria, just as in the hydrodynamic models (Quataert & Gruzinov 2000).

It is important to stress several points:

(1) Although we believe that our linear stability analysis sheds important light on the physics of radiatively inefficient accretion flows, it does not directly address the fully developed nonlinear turbulence that is of primary importance to the accretion flow structure. This requires numerical simulations.

(2) The picture proposed here is very similar to the hydrodynamic models of CDAFs; the primary change is that we have theoretical support from a magnetohydrodynamic analysis.

(3) Balbus (2001) has shown that magnetized dilute plasmas can be convectively unstable in the presence of an outwardly decreasing temperature, rather than just an outwardly decreasing entropy. It would be interesting to incorporate this into future work.

(4) The extent to which one can usefully distinguish scales in the power spectrum that are convective from those that are magnetorotational depends on $`\beta `$, and becomes more and more useful if $`\beta 1`$. It would be straightforward if the convective and MRI modes occupied distinct and well-separated regions of wavevector space. Then even the nonlinear physics of the instabilities could be distinct. To cite an example from plasma physics, heat transport in fusion devices is due to instabilities on both proton and electron Larmor radii scales and one can often distinguish the physics of each mechanism separately even in the nonlinear regime (e.g., Dorland et al. 2000; Rogers, Dorland & Kotschenreuther 2000). As Figure 1 shows, such a separation in wavevector space is not present for our problem. On the other hand, Figs. 1 and 2 do show quite clearly that the physics of the convective modes is very different from that of the MRI modes.

(5) For a Høiland-unstable medium ($`N^2>\kappa ^2`$), there is a nice continuity between the hydrodynamic limit and the strong field MHD limit. In wavevector space, the field-affected modes (the MRI modes) are restricted to larger values of $`k`$ where magnetic tension is important. For very weak fields, these modes are at very large $`k`$, and most of the long-wavelength modes that probably dominate the dynamics are convective. As the field strength increases, however, the MRI modes expand downwards in $`k`$ space and become more and more important. Another way to look at this is as follows. When the fluid is unstable according to the Høiland criterion, the switch that determines whether or not convection is important is not the presence or absence of a magnetic field, but rather the magnitude of the dimensionless parameter $`kv_A/\mathrm{\Omega }kH/\beta ^{1/2}`$. When $`kv_A/\mathrm{\Omega }`$ is small, convection dominates the physics, whereas when it is
>1
>absent1\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{$\m@th\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}}1, the magnetic field dominates.

## 4 Conclusion

In this paper we have used linear stability theory to assess the conditions under which a hydrodynamic analysis is appropriate for an intrinsically magnetohydrodynamic accretion flow. We have shown that, when the medium is unstable by the Høiland criteria, long-wavelength instabilities with $`kv_A\mathrm{\Omega }`$ are effectively hydrodynamical. Physically, this is because magnetic tension is unimportant for sufficiently long-wavelength perturbations. Formally, the dispersion relation for a magnetohydrodynamic plasma reduces in this limit to that of hydrodynamics (§2). This shows that all of the properties of hydrodynamic instabilities carry over to an MHD analysis in the long-wavelength limit.

For the particular case of radiatively inefficient accretion flows that are unstable to both convection and the MRI, the relevant parameter that distinguishes convective behavior from magnetorotational behavior is the quantity $`kv_A/\mathrm{\Omega }H/(\lambda \beta ^{1/2})`$. Long-wavelength modes for which $`\lambda /H\beta ^{1/2}`$ are convective, while modes for which $`\lambda /H\beta ^{1/2}`$ are magnetorotational. The convective modes are primarily driven by buoyancy, not magnetic tension (Fig. 2). Moreover, they transport angular momentum inwards rather than outwards (Fig. 1b), just as indicated by hydrodynamic studies. For the long-wavelength modes to be of interest the accretion flow must have $`\beta 1`$ so that the convective fluctuations fit inside the disk. Under these conditions we believe that the hydrodynamic CDAF models can accurately describe the structure of the accretion flow. The precise $`\beta `$ that allows a hydrodynamical analysis is uncertain and can be determined only by magnetohydrodynamic simulations (the linear analysis suggests that β
>10
>𝛽10\beta\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{$\m@th\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}{\vbox{\offinterlineskip\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr} }}}10 may suffice).

It is interesting to note that the saturation $`\beta `$ of the MRI does not appear to be universal. Given the plausible dependence on $`\beta `$ highlighted above, this may have important implications for the structure of radiatively inefficient accretion flows. Three-dimensional MHD simulations find that that if the field is initially toroidal a radiatively inefficient accretion flow develops into a CDAF-like state (Igumenshchev, Abramowicz & Narayan 2002), while if the field has a vertical component the structure is very different (Hawley, Balbus & Stone 2001; Hawley & Balbus 2002; see also Stone & Pringle’s 2001 2D simulations).

To conclude, it is important to emphasize that the primary conclusion of the analytic arguments in this paper is not that CDAFs actually do describe the structure of radiatively inefficient accretion flows, but rather that they could. In particular, our linear MHD analysis supports all of the key CDAF assumptions so long as $`\beta 1`$.

###### Acknowledgements.

RN, IVI, and EQ were supported in part by NSF grant AST-9820686, RFBR grant 00-02-16135, and NASA grant NAG 5-12043, respectively.

## Appendix A Appendix

In the main text we focused exclusively on unstable modes and therefore considered only the positive sign of the root in equation (2). Here we generalize the discussion and consider both unstable and stable modes. To this end, we rewrite equation (2) as follows

$$\gamma ^2=(kv_A)^2+\frac{(N^2\kappa ^2)}{2}\left[1\pm \sqrt{1+16\mathrm{\Omega }^2(kv_A)^2/(N^2\kappa ^2)^2}\right].$$
(A1)
The + and $``$ signs in equation (A1) map directly onto the + and $``$ signs in equation (2) when $`N^2\kappa ^2>0`$, but the mapping is reversed when $`N^2\kappa ^2<0`$. Equation (A1) is, of course, completely equivalent to equation (2), but for this Appendix we find it more convenient to work with the $`\pm `$ convention in equation (A1) because of its behavior when $`N^2\kappa ^2`$ changes sign. In what follows, we refer to the + sign in equation (A1) as the “+ mode” and the $``$ sign as the “$``$ mode.”

Let us focus on the long-wavelength limit $`kv_A/\mathrm{\Omega }1`$. The growth rates of the + and $``$ modes are then given by

$$\gamma _+^2(N^2\kappa ^2)+\left(\frac{4\mathrm{\Omega }^2}{N^2\kappa ^2}1\right)(kv_A)^2,$$
(A2)
$$\gamma _{}^2\left(\frac{4\mathrm{\Omega }^2}{N^2\kappa ^2}+1\right)(kv_A)^2.$$
(A3)
In the limit $`kv_A/\mathrm{\Omega }1`$, the growth rate of the + mode is the same as in hydrodynamics: $`\gamma _+^2=N^2\kappa ^2`$. If $`N^2>\kappa ^2`$ (Høiland unstable), $`\gamma _+^2`$ is positive and the mode is unstable, while if $`N^2<\kappa ^2`$, $`\gamma _+^2`$ is negative and the mode is oscillatory in nature. In either case, the mode is completely insensitive to the magnetic field.

The $``$ mode, on the other hand, has a growth rate that vanishes as $`kv_A0`$, which is a clear indication that this mode is strongly influenced by the magnetic field. Indeed, the mode owes its very existence to the presence of the magnetic field. The $``$ mode is associated with the classic MRI of BH91.

An inspection of equations (A2) and (A3) reveals the following interesting result: assuming that $`|N^2\kappa ^2|<4\mathrm{\Omega }^2`$ (which is almost always the case), only one of the two modes can be unstable. When the medium is convectively stable ($`N^2<\kappa ^2`$), the + mode (which we argued above is hydrodynamic in nature) is obviously stable. In this case, the $``$ mode is unstable and gives the long-wavelength limit of the MRI. However, when the medium is convectively unstable ($`N^2>\kappa ^2`$), then it is the + mode that grows, while the $``$ mode is stable. Thus, in the long wavelength limit, the MRI actually disappears if the medium is unstable by the Høiland criterion. Fortunately, this is precisely when the medium is convectively unstable so there are always unstable modes available, regardless of the sign of $`N^2\kappa ^2`$.

Figure (3) shows $`\gamma ^2`$ as a function of $`(kv_A)^2`$ for the + and $``$ modes for representative values of $`N^2/\mathrm{\Omega }^2`$. We see that modes that are stable ($`\gamma ^2<0`$) in the limit $`kv_A0`$ remain stable for all values of $`k`$, while modes that are unstable at small values of $`kv_A`$ reach a maximum growth rate for $`kv_A\mathrm{\Omega }`$ and then become stable with increasing $`kv_A`$.


text:
# Milli-arcsecond scale Rotation Measure in the CSS Quasars 0548+165 and 1524-136

## 1 Introduction

Compact Steep-spectrum Sources (CSSs, e.g. Pearson et al. 1985; Fanti et al. 1990) are physically small objects with sub-galactic dimensions. Their structures and sizes are possibly affected by the ambient gas in the central regions of the parent optical objects. At sub-arcsecond resolution these sources often show strongly distorted structures with recognizable jet-like features (Fanti et al. 1986; Spencer et al. 1991), consistent with strong dynamical interactions between the jets and the ambient media. Most CSSs show low ($``$ 1$`\%`$) percentage polarization at or below 5 GHz (Saikia 1991; Saikia, Singal & Cornwell 1987). However, the median polarization increases with frequency (van Breugel et al. 1992; Saikia 1995, Saikia et al. 2001) suggesting that the observations at lower frequencies are affected by Faraday depolarization. A number of CSSs have very high RMs, of the order of several thousand rad m-2. To explain both the small linear sizes and the low polarizations at centimeter wavelengths, it has been postulated that CSSs are cocooned in dense gaseous envelopes (e.g. Mantovani et al. 1994). Supporting evidence for such envelopes comes from optical spectroscopy (O’Dea 1998).

The largest structural distortions in CSSs are often seen in sources dominated by jets, though this may be due in part to a selection effect caused by projection effects (e.g. Spencer 1994). However, images of CSSs do indicate that intrinsic source distortions are caused by interactions with a dense, inhomogeneous gaseous environment. This view is supported by the fact that the most distorted and complex structures are found in objects with very weak cores. Additionally, apparent superluminal motions are rare among CSSs, the only known cases to exhibit such motions being 3C147 (Alef et al. 1990) and 3C138 (Cotton et al. 1997).

In addition to distorted structures, increased asymmetry in the location of the outer components also suggests that these components are evolving through a dense asymmetric environment in the central regions of the galaxies (Sanghera et al. 1995; Saikia 1995). Sub-arcsecond polarimetry has provided evidence in favour of the interaction of these components with dense clouds of gas. For example, the southern component of 3C147 has a much higher RM than the northern component, is brighter, closer to the nucleus, and has the expected signature of a jet colliding with a cloud of gas on the southern side of the galaxy (Junor et al. 1999).

As noted by Murgia et al. (1999), about 10-15% of CSSs, mainly quasars, are jet-dominated with complex or highly asymmetric structures and may represent a different class of “frustrated” objects where strong jet-cloud interactions are at work.

In most of the core-jet structures, the derived magnetic field configurations are parallel to the source major axis, as is found for VLBI-scale jets in quasars (Cawthorne et al. 1993). In core-dominated QSOs, high rest-frame Faraday rotation measures ($`>10^3`$ rad m-2) are found in the inner 20 pc of the nuclear region (Taylor 1999).

Polarimetric observations of CSSs provide a useful probe of the physical conditions in the gaseous environment of these young sources. However, the number of CSSs for which detailed polarization information is available is still rather small. To improve this situation, we are conducting a series of observations to image several sources from a list of CSSs having reasonable degrees of polarized emission and clear signatures of interaction between radio emission and the environment, namely fractional polarizations that decrease with increasing wavelength, and high values of rotation measure (RM $`>`$450 rad m-2 in the source rest frame).

### 1.1 The two target sources

So far, most of the sources observed in our candidate list are QSOs. The resulting images will be useful for a direct comparison with the polarization structures found in core-dominated QSOs. Here we will present and discuss new VLBA polarization observations of the two CSS quasars, 0548$`+`$165 and 1524$``$136.

In arcsecond-scale images, 0548$`+`$165 shows an asymmetric structure with a strong unresolved component coincident with a quasar at redshift 0.474, and a faint secondary component 3 arcsec northward. Most of the polarized emission comes from the main component. This emission is strongly depolarized going from 15 to 5 GHz, with a rest-frame rotation measure of 1934 rad m-2 (Mantovani et al., 1994). The magnetic field direction derived from the derotated electric vectors is aligned almost east-west, which is also the direction of the inner part of the milli-arcsecond jet. Images obtained with the European VLBI Network at 1.6 and 5 GHz (Mantovani et al., 1998) show a collimated thin jet with a noticeable bend at 80 mas ($`253`$ pc) from the core. However, the jet does not lose its collimation after the bend and can be tracked for an additional 75 mas. It is worth mentioning that while the faint component is located 3 arcsec northward, the jet points South.

1524$``$136 has a steep spectrum classified as C- (i.e. inverted at low frequency) by Steppe et al. (1995) from observations made in the range 90 MHz–230 GHz. The Einstein 0.2–4.5 keV X-ray luminosity is 2.7 $`10^{46}`$ erg sec-1 (Wilkes et al. 1994). Images from VLA observations of 1524$``$136 (Mantovani et al. 1994) show a rather compact object which is resolved into a double source, at 15 GHz. The separation between the two components is $``$330 mas. 1524$``$136 has a rest-frame RM=$``$840 rad m-2 and does not depolarise between 15 and 5 GHz.

Table 1 summarizes the physical parameters derived for the two sources from these previous observations. The source name is followed by the Optical Identification (OI), Redshift, Apparent Visual Magnitude, log(Spectral Power) at 5 GHz in W Hz-1, Spectral Index ($`S_\nu \nu ^\alpha `$), Linear Size (LS) in pc (H0 = 100 km s-1 Mpc-1 and q0 = 0.5).

## 2 Observations and data analysis

The observations were carried out with the VLBA 111The Very Long Baseline Array and the Very Large Array are facilities of the National Radio Astronomy Observatory, USA, operated by Associated Universities Inc., for the National Science Foundation. and one antenna of the VLA, recording both right- and left-circular polarization (RCP, LCP) and using 1-bit sampling. The observations are summarized in Table 2. Amplitude calibration was derived using both measurements of the system temperatures made during the observations and knowledge of the antenna gains of each element of the array. Complex correlation coefficients were recovered at the Array Operations Center, Socorro (NM, USA) correlator. The recovered complex correlation coefficients are corrupted principally by phase gradients in frequency and time. These are corrected by the use of global fringe-fitting (Schwab & Cotton, 1983). Polarization calibration was performed following Cotton (1993). Calibration and subsequent imaging were done in NRAO’s $`𝒜𝒫𝒮`$ analysis package. The RCP$``$LCP delay difference was derived by fringe-fitting a short segment of the cross-hand data from a strong calibrator. Strong calibrators (DA193, OQ208, 3C84, 3C138, 3C286, 3C273, 3C345) were used to determine the instrumental polarizations of the antennas. The polarization angle calibration used observations of 3C286 at 8 GHz and of 3C273 at the four frequencies in the 5-GHz band. Low resolution images of 3C286 were made from the calibrated 8.4-GHz data using the shortest baselines in the array. The polarization position angle (PA) of the electric vector of the linearly polarized emission was then compared with that measured in VLA images of 3C286. The 8.4-GHz data for each source, calibrated in complex amplitude and phase, were edited and averaged in frequency and time to improve the signal-to-noise ratio. These data were subsequently used in the standard iterative self-calibration and imaging process. Polarized images of un-polarized calibrator sources were made to check our residual instrumental polarization. These images showed no systematic polarization miscalibration and the peak polarization signal was $`<`$0.3% of the peak total intensity signal.

Due to logistical problems, a slightly different PA calibration strategy was used at 5 GHz. Here, four images, one for each of the observing frequencies, were made of the source 3C273. This source has a small rotation measure (RM$`=80`$ rad m-2; Taylor et al. 1998) for the strongest polarized component along the jet. We then derived the PA values for each frequency from a small box (5$`\times `$5 pixels) around the peak of the polarized emission along the jet of 3C273. The LCP-RCP phase differences were corrected such as to place the PAs on a slope of 80 rad m-2 in a plot of PA against $`\lambda ^2`$. The raw results from our analysis are shown in Fig. 1. The required corrections applied to the PA for each IF were as follows: $`\mathrm{\Delta }`$PA$`{}_{4615}{}^{}=6.9^{}`$, $`\mathrm{\Delta }`$PA$`{}_{4653}{}^{}=45.5^{}`$, $`\mathrm{\Delta }`$PA$`{}_{4850}{}^{}=3.5^{}`$, $`\mathrm{\Delta }`$PA$`{}_{5090}{}^{}=0.8^{}`$. The calibrated data were then imaged and the polarization PA recovered. These agree well with the nearly contemporaneous integrated values in the UMRAO database 222This research has made use of data from the University of Michigan Radio Astronomy Observatory which is supported by funds from the University of Michigan. (UMRAO 2001).

## 3 Results

### 3.1 Total intensity and polarization images

Naturally-weighted images were made separately at the five observing frequencies. Natural weighting gives the best sensitivity at the expense of some resolution. To allow direct comparison between the images at different frequencies, the images were then convolved to a common effective resolution via a two dimensional, circularly-symmetric Gaussian restoring beam with a half-width of 6 mas. The r.m.s. noises in the linearly polarized intensity images are 0.15 and 0.4 mJy beam-1 at 5 and 8.4 GHz respectively. These values are close to the theoretically expected noises.

As examples, Figs. 2 and 3 show respectively the total intensity and linearly polarized intensity images obtained at 4619 MHz for 0548$`+`$165, while Figs. 4 and 5 display the images for 1524$``$136 at 4854 MHz. For the latter source, we concentrate here on the bright polarized jet-like component to the North. The wide field images at 8.4 GHz will be presented elsewhere.

The polarized and total flux density for five regions along the jet of 0548+165 (see Fig. 2) have been measured and are presented in Table 3 together with the Spol / Stot ratio for each region. The Q map and the U map were integrated over a region separately and then the polarized flux density, S$`{}_{pol}{}^{}=(Q^2+U^2)^{0.5}`$, derived. Stot and Spol are in mJy, the frequency $`\nu `$ in GHz. In regions 2 and 3 the polarized flux density is below the 3-$`\sigma `$ noise level. The same information for 1524$``$136 has been collected in Table 4 for the regions marked A-F on Fig. 5.

### 3.2 Rotation Measures

Three regions of linearly polarized emission, labelled 5A, 5B and 5C on Fig. 3, have been taken into consideration for the source 0548+165. The PA values for each of the IFs have been measured considering the distribution of values in a small box (5$`\times `$5 pixels) around the peak of the polarized emission. The standard deviation for each value was also derived and added in quadrature to the error in the PA correction applied by the calibration procedure (3). The PA values and their standard deviations from the mean are given in Table 5, together with the total flux densities, the linearly polarized flux densities and the polarization percentages of the emission from the selected regions.

In Table 6 the same information have been collected for six regions, A-F, within 1524$``$136. It should be noted that two of the four images made at the different C-band frequencies had better uv coverages, and consequently a higher dynamic ranges, due to the presence in the array of a VLA antenna. This causes the differences in the total flux densities for each region of the two sources seen between the four C-band frequencies listed in Tables 3 and 4.

The PA vs $`\lambda ^2`$ plots for the selected regions are presented in Fig. 6 for the source 0548+165, and in Fig. 7 for the source 1524$``$136. In Table 7, the derived RMs, the intrinsic PAs at the point of emission plus the RM in the rest frame ($`rf`$, i.e. corrected for the redshift) are presented for the two sources. In each plot the best fit to the data is shown. The RMs are estimated by fitting the points with a linear least-squares fit. The correlation coefficients are $`0.98`$ in all cases but for regions B and D in 1524$``$136 where we found 0.95. Subtracting $`180^{}`$ from the PAs of the 8.4 GHz point for the six regions selected in 1524$``$136, the correlation coefficients drop to values $`<0.61`$ in four cases. We found 0.82 for region B and 0.94 for region D. In conclusion, only in region D do the RM values reported in Table 7 present an ambiguity.

## 4 Discussion

### 4.1 0548$`+`$165

The nucleus of 0548+165 is almost certainly component 1 (see Fig. 2). Using the flux density values derived from the 1.6-GHz observations by Mantovani et al. (1998) and from the present observations at 5 and 8.4 GHz (Table 3), a spectral index of $`\alpha 0.2`$ is found for component 1. The remaining components have steeper spectral indices. The projected linear size for 0548$`+`$165 from the nucleus to the region where the jet shows an $`90^{}`$ bend is $``$330 pc. Regions, 5A, 5B and 5C all show large RMs, much larger than the typical RM$``$100 rad m-2 due to the ISM of our Galaxy in this direction (Simard-Normandin, Kronberg & Button 1981).

The magnetic field direction derived from the observations lies parallel to the tangent of the total intensity contours as we move around the $`90^{}`$ bend. This suggests that the field lines are compressed by an interaction of the jet with a dense cloud. It is also of interest that this interaction does not appear to significantly disrupt or perturb the jet from a smooth flow. This is also supported by sensitive (Mk3, Mode A) $`\lambda 6`$-cm VLBI observations (Mantovani et al. 1998).

The largest value of RM is found in region 5C, where the jet is most bent. The value of RM then decreases going outwards along the jet from regions 5B to 5A. In the three regions, the PAs are rotated by angles larger than $`\pi `$/2 by 6-cm band. In addition, in the region where the source structure suggests that an interaction between the jet and the ISM is taking place, the depolarization between 8.4 and 4.6 GHz is p4.6/p8.4=0.26. This indicates (Burn 1966) that the rotation is due to a ‘thin slab’ plus something else, such as a two-phase foreground medium, i.e. narrow line clouds and a hot medium are present.

The polarized signal lies predominantly on the western edge of the bent jet i.e. on the side where the putative collision occurs. This may be explained simply by compression of magnetic field lines by the collision, with consequent enhancement of both the total intensity and the linearly polarised emission.

The depolarization produced by a Faraday screen is given by:

$$p(\lambda )=(1+8\sigma ^2\lambda ^4t^2/s_{}^2)^{1/2}$$
(1)
(Tribble 1992) where $`p`$ is the fractional polarization, $`\lambda `$ the wavelength, $`\sigma `$ the RM dispersion assuming a stochastic Gaussian distribution, $`t`$ the resolution, and $`s_{}`$ the RM scale fluctuation.

Assuming $`s_{}`$ to be of the order of the separation between the peaks of polarized emission ($`7.5`$ mas), $`\sigma =`$ 1372 rad m-2, i.e. the dispersion in the rest frame, and $`t=`$ 7 mas, a value p4.6/p8.4 = 0.32 for the depolarization is obtained, in good agreement with the value measured in the area of the bend (see Table 3). It follows that the observed depolarization might be due to beam depolarization since the RM scale fluctuation is of the order of the resolution achieved by the present observations.

The RM values in regions 5A and 5B are similar (at the 3-$`\sigma `$ level) and can be understood as being due to an external screen. In region 5C, the RM excess might be due to the presence of a NLR, which could also be responsible for the large bend in the jet. In fact, 3-D simulations of jet-cloud interactions (Norman & Balzara 1993) have shown that a jet can be deflected by angles up to 90 preserving its collimation after the bend. De Young (1991) showed that a jet can be deflected by interacting with a cloud in the NLR on a time scale of $`10^6`$ yr at a distance of 2–4 kpc from the nucleus.

Making the assumptions that the external screen is responsible for the RM of regions 5A and 5B, and partially in region 5C, and that the magnetic field, $`B`$, along the line of sight and the thermal electron density, $`n_e`$, are almost constant, the rotation measure is RM = $`0.8n_eB_{}l`$. With RM = 1200 rad m-2, and a $`l`$ = 1 kpc thick screen, the product $`n_eB_{}`$ is 1.5 cm$`{}_{}{}^{3}\mu `$G. The typical range for the density of thermal electron is $`n_e`$ 0.1 – 1 cm-3, suggesting that $`B_{}1.515\mu `$G.

Similar values for the magnetic field were derived by Venturi & Taylor (1999) for the bending jet of 3C216.

### 4.2 1524$``$136

The region of jet-cloud interaction in the case of 1524$``$136 is situated 475 pc from the nucleus. The RM value changes going from component A to E making the astrophysical analysis more intriguing. The peak of polarized emission is found in region A. The jet changes direction in region B where, however, the RM is lower. Concerning depolarization, this is seen for regions A$`+`$B and C, while region F shows no depolarization.

Following Tribble (1992), 1524$``$136 has a depolarization p4.6/p8.4 = 0.33 assuming $`s_{}=24`$ mas (mean distance between components), $`t=`$ 6 mas and $`\sigma =`$ 2577 rad m-2. Integrating the total flux and the polarized flux over the whole source, a similar value for the depolarization is found, i.e. p4.6/p8.4 = 0.28. This result might suggest that the depolarization can be due to inhomogenities in the structure of the magnetic field responsible for the changes in the RM, caused by the external screen, on a comparable scale to that of the distances between the components. Since the differences between PAs is $`>90^{}`$, we can say that the RM is due to an external screen in which the orientation of the magnetic field along the line of sight and/or the electron density, do vary along the jet in such a way that the observed variation in RMs can be reproduced.

The intensity of the magnetic field along the line of sight for each region, assuming $`n_e0.11`$ cm-3 and a 1 kpc thick screen, is summarized in Table 8. The values obtained for $`B_{}`$ in 1524$``$136 are larger than those found for the screen in 0548$`+`$165.

## 5 Conclusions

In this paper we have presented results obtained for the CSSs 0548+165 and 1524$``$136 from VLBA+VLA1 observations at five frequencies. These sources were the first observed in a larger project planned to study CSSs which have polarization percentages that decrease with decreasing frequency and high values of rotation measure (RM$`>`$450 rad m-2) on the arcsecond scale.

Both sources show polarized emission along their mas-scale jets, while weak polarized emission is detected in their cores. We were able to estimate rotation measures by a five-point linear fit to PA vs $`\lambda ^2`$ plots in three and six regions along the jets of 0548+165 and 1524$``$136 respectively. The RMs obtained show very high values in both sources (RM$`>`$ 4000 and up to 104 rad m-2 respectively). We suggest that the high RMs are produced by a foreground screen in 1524$``$136 and by a foreground screen plus a RM excess due to the effect of a NLR cloud for 0548+165. In recent HST optical images of CSSs, a central compact component of $`<`$300 mas, dominated by forbidden \[OIII\]$`\lambda \lambda 4959,5007`$ Å and \[OIII\]$`\lambda 3727`$ Å emission lines, is detected. Compact emission line regions with dimensions typical of the NLR of Seyfert galaxies are also observed (Axon et al. 2000).

Estimates of the magnetic field intensity along the line of sight in the foreground screen responsible for Faraday rotation have been made on the basis of thermal electron density and screen depth. Assuming $`n_e0.1\text{cm}^3`$ for the thermal electron density and $`l1\text{kpc}`$ for the Faraday screen depth, we find a magnetic field intensity value of $``$ 15 $`\mu `$G in the foreground screen. A magnetic field in the range 40–125 $`\mu `$G was found in 1524$``$136. The electron density and Faraday-screen depth estimates used, are uncertain (by a factor of ten) because of lack of information on X–ray emission and on the position angle between the jet axis and the observer’s line of sight.

The observed depolarization is consistent with that derived using the Tribble (1992) model which analyzes the effects of a Faraday screen. This suggests that the depolarization is probably due to beam depolarization in 0548+165 and most likely due to inhomogeneities in the magnetic field structure of the foreground screen in 1524$``$136.

The derived magnetic field directions follow the jets in both sources, even when the source is strongly bent due to the interaction with the dense external medium. However, the jets maintain their collimation despite a possible strong interaction with a NLR.

Finally, we note the difference between the distribution of polarised emission in this pair of two jet-dominated QSOs and that found in core-dominated QSOs. For the latter class, the highest RMs are found in the inner 20 pc of the nuclear region (Taylor 1999). In the CSSs 0548+165 and 1524$``$136 the most polarised regions are along the jet while the nuclei are weakly polarised.

###### Acknowledgements.

The authors whish to thank the referee Dr. G.B. Taylor for valuable comments on the manuscript.


text:
# Search for Charmless Two-body Baryonic Decays of 𝐵 Mesons

## Abstract

We report the results of a search for the rare baryonic decays $`B^0p\overline{p}`$, $`\mathrm{\Lambda }\overline{\mathrm{\Lambda }}`$, and $`B^+p\overline{\mathrm{\Lambda }}`$. The analysis is based on a data set of $`31.7\times 10^6B\overline{B}`$ events collected by the Belle detector at the KEKB $`e^+e^{}`$ collider. No statistically significant signals are found, and we set branching fraction upper limits $`(B^0p\overline{p})<1.2\times 10^6`$, $`(B^0\mathrm{\Lambda }\overline{\mathrm{\Lambda }})<1.0\times 10^6`$, and $`(B^+p\overline{\mathrm{\Lambda }})<2.2\times 10^6`$ at the 90% confidence level.

K. Abe7, K. Abe40, T. Abe41, I. Adachi7, Byoung Sup Ahn15, H. Aihara42, M. Akatsu22, Y. Asano47, T. Aso46, V. Aulchenko1, T. Aushev12, A. M. Bakich37, Y. Ban33, E. Banas27, S. Behari7, P. K. Behera48, A. Bondar1, A. Bozek27, M. Bračko20,13, J. Brodzicka27, T. E. Browder6, P. Chang26, Y. Chao26, B. G. Cheon36, K.-F. Chen26, R. Chistov12, Y. Choi36, L. Y. Dong10, J. Dragic21, A. Drutskoy12, S. Eidelman1, V. Eiges12, Y. Enari22, C. W. Everton21, F. Fang6, H. Fujii7, C. Fukunaga44, M. Fukushima9, N. Gabyshev7, A. Garmash1,7, T. Gershon7, A. Gordon21, R. Guo24, J. Haba7, H. Hamasaki7, F. Handa41, K. Hara31, T. Hara31, N. C. Hastings21, H. Hayashii23, M. Hazumi7, E. M. Heenan21, T. Higuchi42, T. Hojo31, T. Hokuue22, Y. Hoshi40, S. R. Hou26, W.-S. Hou26, S.-C. Hsu26, H.-C. Huang26, Y. Igarashi7, T. Iijima7, H. Ikeda7, K. Inami22, A. Ishikawa22, R. Itoh7, H. Iwasaki7, Y. Iwasaki7, P. Jalocha27, H. K. Jang35, J. H. Kang51, J. S. Kang15, N. Katayama7, H. Kawai2, H. Kawai42, Y. Kawakami22, H. Kichimi7, D. W. Kim36, Heejong Kim51, H. J. Kim51, H. O. Kim36, Hyunwoo Kim15, S. K. Kim35, T. H. Kim51, K. Kinoshita4, S. Korpar20,13, P. Križan19,13, P. Krokovny1, R. Kulasiri4, S. Kumar32, A. Kuzmin1, Y.-J. Kwon51, J. S. Lange5,34, G. Leder11, S. H. Lee35, D. Liventsev12, R.-S. Lu26, J. MacNaughton11, G. Majumder38, F. Mandl11, T. Matsuishi22, S. Matsumoto3, Y. Mikami41, W. Mitaroff11, K. Miyabayashi23, Y. Miyabayashi22, H. Miyake31, H. Miyata29, G. R. Moloney21, S. Mori47, T. Mori3, T. Nagamine41, Y. Nagasaka8, Y. Nagashima31, T. Nakadaira42, E. Nakano30, M. Nakao7, J. W. Nam36, Z. Natkaniec27, K. Neichi40, S. Nishida16, O. Nitoh45, S. Noguchi23, T. Nozaki7, S. Ogawa39, F. Ohno43, T. Ohshima22, T. Okabe22, S. Okuno14, S. L. Olsen6, W. Ostrowicz27, H. Ozaki7, P. Pakhlov12, H. Palka27, C. S. Park35, C. W. Park15, H. Park17, L. S. Peak37, J.-P. Perroud18, M. Peters6, L. E. Piilonen49, N. Root1, M. Rozanska27, K. Rybicki27, H. Sagawa7, Y. Sakai7, H. Sakamoto16, M. Satapathy48, A. Satpathy7,4, O. Schneider18, S. Schrenk4, C. Schwanda7,11, S. Semenov12, K. Senyo22, M. E. Sevior21, H. Shibuya39, B. Shwartz1, V. Sidorov1, J. B. Singh32, S. Stanič47,⋆, A. Sugi22, K. Sumisawa7, T. Sumiyoshi7, K. Suzuki7, S. Suzuki50, S. Y. Suzuki7, S. K. Swain6, T. Takahashi30, F. Takasaki7, M. Takita31, K. Tamai7, N. Tamura29, J. Tanaka42, M. Tanaka7, G. N. Taylor21, Y. Teramoto30, S. Tokuda22, M. Tomoto7, T. Tomura42, S. N. Tovey21, K. Trabelsi6, T. Tsukamoto7, S. Uehara7, K. Ueno26, Y. Unno2, S. Uno7, K. E. Varvell37, C. C. Wang26, C. H. Wang25, J. G. Wang49, M.-Z. Wang26, Y. Watanabe43, E. Won35, B. D. Yabsley7, Y. Yamada7, M. Yamaga41, A. Yamaguchi41, H. Yamamoto41, Y. Yamashita28, M. Yamauchi7, S. Yanaka43, J. Yashima7, P. Yeh26, K. Yoshida22, Y. Yuan10, Y. Yusa41, C. C. Zhang10, J. Zhang47, Y. Zheng6, V. Zhilich1, and D. Žontar47

1Budker Institute of Nuclear Physics, Novosibirsk

2Chiba University, Chiba

3Chuo University, Tokyo

4University of Cincinnati, Cincinnati OH

5University of Frankfurt, Frankfurt

6University of Hawaii, Honolulu HI

7High Energy Accelerator Research Organization (KEK), Tsukuba

8Hiroshima Institute of Technology, Hiroshima

9Institute for Cosmic Ray Research, University of Tokyo, Tokyo

10Institute of High Energy Physics, Chinese Academy of Sciences, Beijing

11Institute of High Energy Physics, Vienna

12Institute for Theoretical and Experimental Physics, Moscow

13J. Stefan Institute, Ljubljana

14Kanagawa University, Yokohama

15Korea University, Seoul

16Kyoto University, Kyoto

17Kyungpook National University, Taegu

18IPHE, University of Lausanne, Lausanne

19University of Ljubljana, Ljubljana

20University of Maribor, Maribor

21University of Melbourne, Victoria

22Nagoya University, Nagoya

23Nara Women’s University, Nara

24National Kaohsiung Normal University, Kaohsiung

25National Lien-Ho Institute of Technology, Miao Li

26National Taiwan University, Taipei

27H. Niewodniczanski Institute of Nuclear Physics, Krakow

28Nihon Dental College, Niigata

29Niigata University, Niigata

30Osaka City University, Osaka

31Osaka University, Osaka

32Panjab University, Chandigarh

33Peking University, Beijing

34RIKEN BNL Research Center, Brookhaven NY

35Seoul National University, Seoul

36Sungkyunkwan University, Suwon

37University of Sydney, Sydney NSW

38Tata Institute of Fundamental Research, Bombay

39Toho University, Funabashi

40Tohoku Gakuin University, Tagajo

41Tohoku University, Sendai

42University of Tokyo, Tokyo

43Tokyo Institute of Technology, Tokyo

44Tokyo Metropolitan University, Tokyo

45Tokyo University of Agriculture and Technology, Tokyo

46Toyama National College of Maritime Technology, Toyama

47University of Tsukuba, Tsukuba

48Utkal University, Bhubaneswer

49Virginia Polytechnic Institute and State University, Blacksburg VA

50Yokkaichi University, Yokkaichi

51Yonsei University, Seoul

on leave from Nova Gorica Polytechnic, Slovenia

\[

\]

The Belle collaboration recently reported the observation of the decay process $`B^+p\overline{p}K^+`$, which is the first known example of a $`B`$ meson decay to a charmless final state containing baryons. In this paper we report the results of a search for the related two-body modes $`B^0p\overline{p}`$, $`\mathrm{\Lambda }\overline{\mathrm{\Lambda }}`$ and $`B^+p\overline{\mathrm{\Lambda }}`$. In the Standard Model, these decays are expected to proceed via color-suppressed $`bu`$ tree diagrams (Figs. 1(a) and (c)) and $`bs,d`$ penguin diagrams (Figs. 1(b) and (d)). The search is based on a 29.4 fb-1 sample of $`e^+e^{}`$ data accumulated at the $`\mathrm{{\rm Y}}(4S)`$ resonance, which contains 31.7 million $`B\overline{B}`$ pairs. A previous search for these decays by the CLEO collaboration using a 5.41 fb-1 sample of $`\mathrm{{\rm Y}}(4S)`$ data yielded 90% confidence-level (C.L.) upper limits.

Belle is a general purpose detector operating at the KEKB asymmetric $`e^+e^{}`$ collider. Tracking information is provided by a silicon vertex detector and a central drift chamber (CDC) in a 1.5 Tesla magnetic field. Hadron identification (PID) for $`\pi /K/p`$ discrimination is obtained from CDC $`dE/dx`$ measurements, aerogel $`\stackrel{ˇ}{\mathrm{C}}`$erenkov counter pulse heights, and timing information from the time-of-flight system. Electron identification is based mainly on CsI(Tl) electromagnetic calorimeter and CDC $`dE/dx`$ information. $`K_L`$ and muons are identified by a system of resistive plate counters interleaved with the iron plates of the flux-return iron yoke.

The event selection criteria are based on tracking and PID requirements, and are optimized using Monte Carlo (MC) simulated event samples.

All primary charged tracks are required to satisfy the following track quality criteria based on the track impact parameters relative to the interaction point (IP), which is determined run-by-run. The $`z`$ axis is defined by the positron beam line. The deviations from the IP position are required to be within $`\pm `$0.05 cm in the transverse ($`x`$-$`y`$) plane, and within $`\pm `$2 cm in the $`z`$ direction. Tracks that satisfy the muon or electron identification requirements are rejected.

Primary proton candidates are selected based on $`p/K/\pi `$ likelihood functions obtained from the hadron identification system. The selection criteria are $`\frac{L_p}{L_p+L_K}>0.6`$ and $`\frac{L_p}{L_p+L_\pi }>0.6`$, where $`L_{p/K/\pi }`$ stands for the proton/kaon/pion likelihood.

$`\mathrm{\Lambda }`$ candidates are reconstructed via the $`p\pi ^{}`$ decay channel and are selected using cuts on four parameters: the angular difference between the $`\mathrm{\Lambda }`$ flight direction and the direction pointing from IP to the decay vertex in the transverse plane; the distance between each track and the IP in the transverse plane; the distance between the decay vertex and the IP in the transverse plane; and the displacement in $`z`$ of the closest approach points of the two tracks to the beam axis. The secondary protons are required to have $`\frac{L_p}{L_p+L_\pi }>0.6`$. The $`p\pi ^{}`$ mass spectrum after the application of the above selection criteria is shown in Fig. 2 for a typical run period. The peak position is consistent with the nominal $`\mathrm{\Lambda }`$ mass  and the mass resolution is about 0.9 MeV/$`c^2`$. Finally, we require the invariant mass of the $`\mathrm{\Lambda }`$ candidate to be within $`\pm `$5 MeV/$`c^2`$ of the nominal $`\mathrm{\Lambda }`$ mass.

Due to the high momentum of the primary particles from these two-body decay modes, the background from generic $`B`$ decays is negligible. This is checked with MC samples of $`B^+B^{}`$ and $`B^0\overline{B^0}`$ pairs where the $`B`$-mesons decay dominantly via $`bc`$ processes. We also checked backgrounds using MC samples of charmless $`B`$ meson decays that include low multiplicity $`B`$ decays into final states with $`\pi `$, $`K`$, $`K^{}`$, $`\rho `$, $`\omega `$, $`\varphi `$, $`\eta `$, and $`\eta ^{}`$ mesons. Only the $`p\overline{p}`$ mode is found to have any background contamination. However, the level of this contamination is negligible, corresponding to less than one event for the current data set.

The main background is from continuum $`q\overline{q}`$ processes. This is confirmed using an off-resonance data set (2.3 fb-1 taken 60 MeV below the $`\mathrm{{\rm Y}}(4S)`$) and a MC sample of 65 million continuum events. These continuum events have a jet-like topology while $`B\overline{B}`$ events are more spherical in the $`\mathrm{{\rm Y}}(4S)`$ center of mass (CM) frame. For continuum event rejection we use $`\mathrm{cos}\theta _T`$, the cosine of the angle between the direction of one primary decay particle and the thrust axis of the non-candidate tracks and showers. This distribution is nearly flat for signal events and is strongly peaked at $`\pm `$1 for continuum background. We also use $`\mathrm{cos}\theta _B`$, the cosine of the angle between the $`B`$ candidate flight direction and the positron beam direction. The signal has a $`\mathrm{sin}^2\theta _B`$ distribution while the background is uniformly distributed. We require the absolute values of $`\mathrm{cos}\theta _T`$ and $`\mathrm{cos}\theta _B`$ to be less than 0.9 for $`\mathrm{\Lambda }\overline{\mathrm{\Lambda }}`$ decays and less than 0.8 for the other modes. For the latter case, the background reduction factor is more than five, while $``$70% of the signal is retained.

We use the following two kinematic variables to identify the reconstructed $`B`$ meson candidates: the beam constrained mass, $`m_{\mathrm{bc}}=\sqrt{E_{\mathrm{beam}}^2p_B^2}`$, and the energy difference, $`\mathrm{\Delta }E=E_BE_{\mathrm{beam}}`$, where $`E_{\mathrm{beam}}`$, $`p_B`$ and $`E_B`$ are the beam energy, the momentum and energy of the reconstructed $`B`$ meson in the $`\mathrm{{\rm Y}}(4S)`$ CM frame, respectively. We retain events with 5.20 GeV/$`c^2`$ $`<m_{\mathrm{bc}}<`$ 5.29 GeV/$`c^2`$ and $`0.2`$ GeV $`<\mathrm{\Delta }E<`$ 0.2 GeV. The signal yield is extracted by maximizing the likelihood function

$$L=e^{(s+b)}\underset{i=1}{\overset{N}{}}[sP_s(m_{\mathrm{bc}_i},\mathrm{\Delta }E_i)+bP_b(m_{\mathrm{bc}_i},\mathrm{\Delta }E_i)],$$

where $`N`$ is the total number of candidate events, $`s(b)`$ denotes the signal (background) yield and $`P_{s(b)}`$ denotes the signal (background) probability density function (pdf).

The $`m_{\mathrm{bc}}`$ and $`\mathrm{\Delta }E`$ signal pdf’s are determined by MC. We use a Gaussian function as the signal pdf for the $`m_{\mathrm{bc}}`$ distribution and a sum of two Gaussians for the $`\mathrm{\Delta }E`$ distribution. The Gaussian parameters (mean and $`\sigma `$) are determined separately for each mode. Background shapes are determined from events in sideband regions of $`\mathrm{\Delta }E`$ and $`m_{\mathrm{bc}}`$ separately. We adopt an empirical function to model the $`m_{\mathrm{bc}}`$ background shape (for events with 0.1 GeV $`<|\mathrm{\Delta }E|<`$ 0.2 GeV) and a first-order polynomial for the $`\mathrm{\Delta }E`$ background shape (for events with 5.20 GeV/$`c^2`$ $`<m_{\mathrm{bc}}<`$ 5.26 GeV/$`c^2`$).

Table 1 summarizes the results. The efficiencies for the $`\mathrm{\Lambda }\overline{\mathrm{\Lambda }}`$ and $`p\overline{\mathrm{\Lambda }}`$ modes include the $`\mathrm{\Lambda }p\pi ^{}`$ branching fraction (64%). The efficiencies are determined from signal MC with the identical event selection and fitting procedure as for the data. Figures 3 and 4 show the $`m_{\mathrm{bc}}`$ (with $`|\mathrm{\Delta }E|<0.05`$ GeV) and $`\mathrm{\Delta }E`$ (with 5.27 GeV/$`c^2`$ $`<m_{\mathrm{bc}}<`$ 5.29 GeV/$`c^2`$) projections for these three modes, respectively. Projections of the fits are shown as smooth curves. No statistically significant signals are found and we determine 90% C.L. upper limits on the signal yields by integrating the likelihood function. We also compute limits using a counting method. We define a signal region by 5.27 GeV/$`c^2`$ $`<m_{\mathrm{bc}}<`$ 5.29 GeV/$`c^2`$ and $`|\mathrm{\Delta }E|<0.05`$ GeV, and treat the number of background events from the maximum likelihood fit as a prediction for the background in this region. We then count the number of events actually observed, and apply the method of  to obtain 90% C.L. intervals.

To estimate the possible influence of fluctuations on the determination of the upper limits, we vary the parameters of the pdf’s by one standard deviation and change the form of the $`\mathrm{\Delta }E`$ signal pdf to a single Gaussian. The relative change of the upper limit is 25%, which is mainly due to the changes in the background shape. The upper limit determination by the counting method is checked by redefining the signal region (varying the $`\mathrm{\Delta }E`$ range from 2$`\sigma `$ to 4$`\sigma `$) and comparing the outcomes. The upper limits from both methods listed in column 3 of Table 1 include these fluctuations. We quote the values from the counting method as the most conservative upper limits. Fit projections with the signal yield fixed at the upper limit for each mode are shown as the superimposed dashed curves in Figs. 3 and 4 for comparison.

The systematic error due to the efficiency of the proton identification ($`p`$-$`K`$ and $`p`$-$`\pi `$) criteria is studied using $`\mathrm{\Lambda }`$ samples. We vary the likelihood ratio requirement for protons and compare the ratio of reconstructed $`\mathrm{\Lambda }`$ yields in data and MC. The overall error is about 3%. We include a 2% error per track to account for the uncertainty in tracking efficiency. The $`\mathrm{\Lambda }`$ reconstruction efficiency is checked by comparing the flight distance distributions of data and MC. They agree very well and no additional error is assigned. The correlated parts of the errors are added linearly to obtain the overall uncertainty in the tracking efficiency and the uncertainty in the PID efficiencies (for $`p`$ and $`\overline{p}`$), then the resulting errors are combined in quadrature. When determining the upper limit for the branching fraction, the efficiency was reduced by one standard deviation. The efficiencies and upper limits for all three decay modes are listed in Table 1.

In summary, we have performed a search for the rare baryonic decays $`B^0p\overline{p}`$, $`\mathrm{\Lambda }\overline{\mathrm{\Lambda }}`$, and $`B^+p\overline{\mathrm{\Lambda }}`$ with 31.7 million $`B\overline{B}`$ events collected by the Belle detector at the KEKB $`e^+e^{}`$ collider. No statistically significant signals are found for these modes, and we set upper limits on their branching fractions at the 90% C.L.. The upper limits are:

$`(B^0`$ $`p\overline{p}`$ $`)<1.2\times 10^6;`$
$`(B^0`$ $`\mathrm{\Lambda }\overline{\mathrm{\Lambda }}`$ $`)<1.0\times 10^6;`$
$`(B^+`$ $`p\overline{\mathrm{\Lambda }}`$ $`)<2.2\times 10^6.`$
These are currently the most stringent limits for these decays. The limit on $`B^0p\overline{p}`$ has been improved by a factor of six compared to the existing bound. These limits are considerably below corresponding charmless mesonic branching fractions which are typically at the level of $`10^5`$. In contrast, the recent Belle measurement of $`B^+p\overline{p}K^+`$ indicates a relatively large branching fraction of $`4.3\times 10^6`$. The underlying physics is now being probed with rapidly accumulating data and new experimental results.

We wish to thank the KEKB accelerator group for the excellent operation of the KEKB accelerator. We acknowledge support from the Ministry of Education, Culture, Sports, Science, and Technology of Japan and the Japan Society for the Promotion of Science; the Australian Research Council and the Australian Department of Industry, Science and Resources; the National Science Foundation of China under contract No. 10175071; the Department of Science and Technology of India; the BK21 program of the Ministry of Education of Korea and the CHEP SRC program of the Korea Science and Engineering Foundation; the Polish State Committee for Scientific Research under contract No. 2P03B 17017; the Ministry of Science and Technology of the Russian Federation; the Ministry of Education, Science and Sport of Slovenia; the National Science Council and the Ministry of Education of Taiwan; and the U.S. Department of Energy.