text:
# 1 Introduction

## 1 Introduction

Solar neutrino experiments, including the measurement of the neutral current rate for solar neutrinos by the SNO collaboration Ahmad:2002jz  provide a solid evidence for solar neutrino conversions Maltoni:2002ni . This has been recently confirmed by the first results from the KamLAND experiment using reactor (anti)-neutrinos :2002dm , Pakvasa:2003zv . Combining the information from reactors with all of the solar neutrino data leads to the best fit point Maltoni:2002aw :

$$\mathrm{tan}^2\theta sol=0.46,\mathrm{\Delta }m^2sol=6.9\times 10^5\mathrm{eV}^2,$$
(1)
confirming that the solar neutrino mixing angle is large, but significantly non-maximal. The $`3\sigma `$ region for $`\theta `$ is:

$$0.29\mathrm{tan}^2\theta sol0.86,$$
(2)
based on a combination of all experimental data. However, one finds a significant reduction of the allowed $`\mathrm{\Delta }m^2sol`$ range. As shown in Ref. Maltoni:2002aw , the pre-KamLAND LMA-MSW region is now split into two sub-regions. At $`3\sigma `$ (1 dof.) one obtains

$$5.1\times 10^5\mathrm{eV}^2\mathrm{\Delta }\mathrm{m}^2\mathrm{sol}9.7\times 10^5\mathrm{eV}^2,1.2\times 10^4\mathrm{eV}^2\mathrm{\Delta }\mathrm{m}^2\mathrm{sol}1.9\times 10^4\mathrm{eV}^2.$$
(3)
Altogether, KamLAND results exclude all oscillation solutions except for the large mixing angle MSW solution (LMA-MSW) to the solar neutrino problem Gonzalez-Garcia:1999aj .

On the other hand, current atmospheric neutrino data require oscillations involving $`\nu _\mu \nu _\tau `$ Fukuda:1998mi . The most recent global analysis gives Maltoni:2002ni ,

$$\mathrm{sin}^2\theta atm=0.5,\mathrm{\Delta }m^2atm=2.5\times 10^3\mathrm{eV}^2$$
(4)
with the 3$`\sigma `$ ranges (1 d.o.f.)

$`0.3\mathrm{sin}^2\theta atm0.7`$ (5)
$`1.2\times 10^3\mathrm{eV}^2\mathrm{\Delta }\mathrm{m}^2\mathrm{atm}4.8\times 10^3\mathrm{eV}^2.`$ (6)
These data have triggered a rush of theoretical and phenomenological papers on models of neutrino masses and mixings, most of which introduce a large mass scale in order to implement various variants of the see-saw mechanism seesaw , Mohapatra:1980yp , Schechter:1981cv . Broken R-parity supersymmetry provides a theoretically interesting and phenomenologically viable alternative to the origin of neutrino mass and mixing Ross:1985yg . Here we focus on the simplest case of supersymmetry with bilinear R-parity breaking Diaz:1997xc . In contrast to the seesaw mechanism, here neutrino masses are generated at the electro-weak scale. Such low-scale schemes for neutrino masses have the advantage of being testable also in accelerator experiments Bartl:2000yh -Hirsch:2002ys  through the decay properties of the lightest supersymmetric particle if the LSP is a neutralino Porod:2000hv , Mukhopadhyaya:1998xj , Choi:1999tq , a slepton Hirsch:2002ys  or a stop Restrepo:2001me , Allanach:1999bf .

Supersymmetric models with explicit bilinear breaking of R-parity (BRpV) deCampos:1995av , Banks:1995by , deCarlos:1996du , Akeroyd:1997iq , BRpVrecent , BRpVother , BRpV\_tau , BRpVmore  provide a simple and calculable framework for neutrino masses and mixing angles in agreement with the experimental data Hirsch:2000ef . In this model the atmospheric neutrino mass scale is generated at tree-level, through an effective ‘low-scale” variant of the seesaw mechanism Ross:1985yg . In contrast, the solar mass and mixings are generated radiatively Hirsch:2000ef . Tree-level neutrino masses within BRpV have been treated extensively in the literature.

This paper is mainly devoted to the solar neutrino masses and mixing. An accurate and reliable calculational method is now necessary in order to confront the model with the new experimental data from KamLAND and other neutrino experiments. A complete one-loop calculation of the neutrino-neutralino mass matrix has been given Hirsch:2000ef  but is rather complex. On the other hand, approximations to the full 1-loop calculation which exist in the literature Chun:2002vp  have not been tested yet against the full calculation. Especially in view of future experimental sensitivities we think such a “benchmark” is important.

In this paper we give an accurate determination of neutrino mass and mixing within an analytical approximation and obtain formulae which can be rather simple, in some cases. For definiteness we will stick to the case of explicit BRpV only. This is the simplest of all R parity violating models. It can be considered either as a minimal three–parameter extension of the MSSM (with no new particles) valid up to some very high unification energy scale, or as the effective description of a more fundamental theory in which the breaking of R-parity is spontaneous Masiero:1990uj , Romao:vu , spo:2 . The latter implies the absence of trilinear R-parity breaking parameters in the superpotential 111Alternatively, such absence may arise from suitable symmetries Mira:2000gg .

This paper is organized as follows. In Sec. 2 we introduce the main features of the model and the relevant mass matrices and corresponding diagonalization matrices. In particular we identify the relevant Feynman graph topologies and rules, and derive approximate formulae for the couplings relevant for the determination of radiatively induced solar neutrino mass scale. We give approximate formulas for the bottom quark/squark loop as well as for the charged scalar loop. In Sec. 3 we check the accuracy of our approximation formulas by a comparison with a full numerical calculation, studying first the role of the simplest bottom-sbottom loop, then the charged scalar loop, before comparing the sum of the two to the full numerical result. In 4 we give simplified approximation formulas for the solar mass and solar mixing angle and conclude and summarize our results in 5.

## 2 BRpV Formalism

In this section we introduce the main features of the model and the relevant mass matrices, and develop approximate formulas, first for couplings and then for the radiative contributions to the neutrino masses due to the exchange of bottom and sbottom quarks, and due to charged scalars and charged fermion loops.

### 2.1 BRpV Model

The minimal BRpV model we are working with is characterized by the presence of three extra bilinear terms in the superpotential analogous to the $`\mu `$ term present in the MSSM

$$W=W_{Yuk}+\epsilon _{ab}\left(\mu \widehat{H}_d^a\widehat{H}_u^b+ϵ_i\widehat{L}_i^a\widehat{H}_u^b\right)$$
(7)
where $`W_{Yuk}`$ includes the usual MSSM Yukawa terms, $`\mu `$ is the Higgsino mass term of the MSSM, and $`ϵ_i`$ are the three new terms which violate R-Parity and lepton number. The smallness of $`ϵ_i`$ may arise dynamically (the product of a Yukawa coupling times a singlet sneutrino vacuum expectation value) in models with spontaneous breaking of R parity Masiero:1990uj .

Alternatively, the smallness of the $`ϵ_i`$ may arise from suitable family symmetries Mira:2000gg . In fact any solution to the $`\mu `$ problem Giudice:1988yz  potentially explains also the “$`ϵ_i`$-problem” Nilles:1996ij . In fact a common origin for the $`ϵ_i`$ terms responsible for the explanation of the neutrino anomalies, and the $`\mu `$ term accounting for electroweak symmetry breaking can be ascribed to a suitable horizontal symmetry that may also predict their ratio, as in Mira:2000gg .

In addition we have the corresponding soft supersymmetry breaking terms in the scalar potential,

$$V_{soft}=V_{soft}^{}+\epsilon _{ab}\left(B\mu H_d^aH_u^b+B_iϵ_i\stackrel{~}{L}_i^aH_u^b\right)$$
(8)
where $`B`$ and the three $`B_i`$ have units of mass and in $`V_{soft}^{}`$ we include all the usual mass and trilinear supersymmetry breaking terms of the MSSM.

### 2.2 Rotation Matrices

If the effective RpV parameters are smaller than the weak scale, we can work in a perturbative expansion defined by $`\xi 1`$, where $`\xi `$ denotes a $`3\times 4`$ matrix given as Hirsch:1998kc

$`\xi _{i1}`$ $`=`$ $`{\displaystyle \frac{g^{}M_2\mu }{2\mathrm{\Delta }_0}}\mathrm{\Lambda }_i`$ (9)
$`\xi _{i2}`$ $`=`$ $`{\displaystyle \frac{gM_1\mu }{2\mathrm{\Delta }_0}}\mathrm{\Lambda }_i`$ (10)
$`\xi _{i3}`$ $`=`$ $`{\displaystyle \frac{ϵ_i}{\mu }}+{\displaystyle \frac{M_{\stackrel{~}{\gamma }}v_u}{4\mathrm{\Delta }_0}}\mathrm{\Lambda }_i`$ (11)
$`\xi _{i4}`$ $`=`$ $`{\displaystyle \frac{M_{\stackrel{~}{\gamma }}v_d}{4\mathrm{\Delta }_0}}\mathrm{\Lambda }_i`$ (12)
where $`\mathrm{\Delta }_0`$ is the determinant of the $`4\times 4`$ neutralino mass matrix, $`M_{\stackrel{~}{\gamma }}=g^2M_1+g^2M_2`$ and

$$\mathrm{\Lambda }_i=\mu v_i+v_dϵ_i$$
(13)
The neutralino/neutrino mass matrix is diagonalized by a $`7\times 7`$ rotation matrix $`𝒩`$ according to

$$𝒩^{}𝐌_{F^0}𝒩^1=𝐌_{F^0}^{\mathrm{d}iag}$$
(14)
and the eigenvectors are given by

$$F_i^0=𝒩_{ij}\psi _j$$
(15)
using the basis $`\psi =(i\lambda ^{},i\lambda ^3,\stackrel{~}{H}_d^1,\stackrel{~}{H}_u^2,\nu _e,\nu _\mu ,\nu _\tau )`$. In this approximation, the rotation matrix can be written as

$$𝒩^{}\left(\begin{array}{cc}N^{}& N^{}\xi ^{}\\ V_\nu ^T\xi & V_\nu ^T\end{array}\right)$$
(16)
Here, $`N`$ is the rotation matrix that diagonalizes the $`4\times 4`$ MSSM neutralino mass matrix, $`V_\nu `$ is the rotation matrix that diagonalizes the tree level neutrino $`3\times 3`$ mass matrix, and $`\xi _{ij}1`$ are the expansion parameters Hirsch:1998kc , Now96 . The terms we need are

$$V_\nu ^T\xi =\left(\begin{array}{cccc}0& 0& b\stackrel{~}{ϵ}_1& 0\\ 0& 0& b\stackrel{~}{ϵ}_2& 0\\ a_1|\stackrel{}{\mathrm{\Lambda }}|& a_2|\stackrel{}{\mathrm{\Lambda }}|& a_3|\stackrel{}{\mathrm{\Lambda }}|+b\stackrel{~}{ϵ}_3& a_4|\stackrel{}{\mathrm{\Lambda }}|\end{array}\right)$$
(17)
where $`b=1/\mu `$,

$`a_1={\displaystyle \frac{g^{}M_2\mu }{2\mathrm{\Delta }_0}},a_2={\displaystyle \frac{gM_1\mu }{2\mathrm{\Delta }_0}},a_3={\displaystyle \frac{M_{\stackrel{~}{\gamma }}v_u}{4\mathrm{\Delta }_0}},a_4={\displaystyle \frac{M_{\stackrel{~}{\gamma }}v_d}{4\mathrm{\Delta }_0}},`$ (18)
The $`\stackrel{~}{ϵ}`$ parameters in eq. (17) are defined as $`\stackrel{~}{ϵ}_i=\left(V_\nu ^T\right)^{ij}ϵ_j`$, and are given by

$`\stackrel{~}{ϵ}_1`$ $`=`$ $`{\displaystyle \frac{ϵ_e(\mathrm{\Lambda }_\mu ^2+\mathrm{\Lambda }_\tau ^2)\mathrm{\Lambda }_e(\mathrm{\Lambda }_\mu ϵ_\mu +\mathrm{\Lambda }_\tau ϵ_\tau )}{\sqrt{\mathrm{\Lambda }_\mu ^2+\mathrm{\Lambda }_\tau ^2}\sqrt{\mathrm{\Lambda }_e^2+\mathrm{\Lambda }_\mu ^2+\mathrm{\Lambda }_\tau ^2}}}`$
$`\stackrel{~}{ϵ}_2`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Lambda }_\tau ϵ_\mu \mathrm{\Lambda }_\mu ϵ_\tau }{\sqrt{\mathrm{\Lambda }_\mu ^2+\mathrm{\Lambda }_\tau ^2}}}`$ (19)
$`\stackrel{~}{ϵ}_3`$ $`=`$ $`{\displaystyle \frac{\stackrel{}{\mathrm{\Lambda }}\stackrel{}{ϵ}}{\sqrt{\mathrm{\Lambda }_e^2+\mathrm{\Lambda }_\mu ^2+\mathrm{\Lambda }_\tau ^2}}}`$
On the other hand the chargino/charged slepton mass matrix is diagonalized with two different $`5\times 5`$ mass matrices,

$$𝒰^{}𝐌_{F^+}𝒱^1=𝐌_{F^+}^{\mathrm{d}iag}$$
(20)
with the eigenvectors satisfying

$$F_{Ri}^+=𝒱_{ij}\psi _j^+,F_{Li}^{}=𝒰_{ij}\psi _j^{}$$
(21)
in the basis $`\psi ^+=(i\lambda ^+,\stackrel{~}{H}_2^1,e_R^+,\mu _R^+,\tau _R^+)`$ and $`\psi ^{}=(i\lambda ^{},\stackrel{~}{H}_1^2,e_L^{},\mu _L^{},\tau _L^{})`$, and with the Dirac fermions being

$$F_i^+=\left(\begin{array}{c}F_{Ri}^+\\ \\ \overline{F_{Li}^{}}\end{array}\right)$$
(22)
To first order in the R-Parity violating parameters we have

$$𝒱\left(\begin{array}{cc}V& V\xi _R^T\\ V_R^{\mathrm{}}\xi _R^{}& V_R^{\mathrm{}}\end{array}\right),𝒰\left(\begin{array}{cc}U& U\xi _L^{}\\ V_L^{\mathrm{}}\xi _L& V_L^{\mathrm{}}\end{array}\right)$$
(23)
where $`V_L^{\mathrm{}}`$ and $`V_R^{\mathrm{}}`$ diagonalize the charged lepton mass matrix according to $`V_L^{\mathrm{}}𝐌^{\mathrm{}}V_R^{\mathrm{}}=𝐌_{\mathrm{d}iag}^{\mathrm{}}`$. For the purposes of our approximate formula, it is sufficient to take $`\xi _R=\mathrm{𝟎}_{2\times 3}`$, because the mixing between right-handed leptons and the charginos is supressed with respect to $`\xi _L`$ by a factor of $`m_l/M_{SUSY}`$ Hirsch:1998kc , Now96 . Note, that we can choose $`V_L^{\mathrm{}}=V_R^{\mathrm{}}=\mathrm{𝟏}_{3\times 3}`$. We then have

$$\xi _L^{i1}=a_1^L\mathrm{\Lambda }_i,\xi _L^{i2}=a_2^L\mathrm{\Lambda }_i+bϵ_i$$
(24)
and

$$a_1^L=\frac{g}{\sqrt{2}\mathrm{\Delta }_+},a_2^L=\frac{g^2v_u}{2\mu \mathrm{\Delta }_+}$$
(25)
where $`\mathrm{\Delta }_+`$ is the determinant of the $`2\times 2`$ chargino mass matrix.

In the BRpV model the charged Higgs fields mix with the charged sleptons forming an $`8\times 8`$ mass matrix Hirsch:2000ef , which is diagonalized by a rotation matrix $`𝐑_{S^\pm }`$. The construction of $`𝐑_{S^\pm }`$ to first order in small (RpV) parameters is quite straightforward but lengthy. The interested reader can find the details in Appendix A.

### 2.3 Approximate Couplings

The relevant Feynman rules for the bottom-sbottom loops are, in the case of left sbottoms:

with

$`O_{Lij}^{bn\stackrel{~}{b}}`$ $`=`$ $`R_{j1}^{\stackrel{~}{b}}h_b𝒩_{i3}^{}R_{j2}^{\stackrel{~}{b}}{\displaystyle \frac{2g}{3\sqrt{2}}}t_W𝒩_{i1}^{}`$
$`O_{Rij}^{bn\stackrel{~}{b}}`$ $`=`$ $`R_{j1}^{\stackrel{~}{b}}{\displaystyle \frac{g}{\sqrt{2}}}\left(𝒩_{i2}\frac{1}{3}t_W𝒩_{i1}\right)R_{j2}^{\stackrel{~}{b}}h_b𝒩_{i3}^{}`$ (26)
where $`t_W=\mathrm{tan}\theta _W`$. After approximating the rotation matrix $`𝒩`$ we find that expressions similar to eq. (26) with the replacement $`𝒩\mathrm{N}`$ are valid when the neutral fermion is a neutralino. When the neutral fermion $`F^0`$ is a neutrino, the following expressions hold

$`O_{Lij}^{bn\stackrel{~}{b}}`$ $``$ $`R_{j1}^{\stackrel{~}{b}}h_b\left(a_3|\stackrel{}{\mathrm{\Lambda }}|\delta _{i^{}3}+b\stackrel{~}{ϵ}_i^{}\right)+R_{j2}^{\stackrel{~}{b}}{\displaystyle \frac{2g}{3\sqrt{2}}}t_Wa_1|\stackrel{}{\mathrm{\Lambda }}|\delta _{i^{}3}`$
$`O_{Rij}^{bn\stackrel{~}{b}}`$ $``$ $`R_{j1}^{\stackrel{~}{b}}{\displaystyle \frac{g}{\sqrt{2}}}\left(\frac{1}{3}t_Wa_1a_2\right)|\stackrel{}{\mathrm{\Lambda }}|\delta _{i^{}3}+R_{j2}^{\stackrel{~}{b}}h_b\left(a_3|\stackrel{}{\mathrm{\Lambda }}|\delta _{i^{}3}+b\stackrel{~}{ϵ}_i^{}\right)`$ (27)
where $`i^{}=i4`$ label one of the neutrinos. $`R_{jk}^{\stackrel{~}{b}}`$ are the rotation matrices connecting weak and mass eigenstate basis for the scalar bottom quarks. In case of no intergenerational mixing in the squark sector $`R_{jk}^{\stackrel{~}{b}}`$ can be parameterized by just one diagonalizing angle $`\theta _{\stackrel{~}{b}}`$.

The relevant Feynman rule for the charged Higgs/slepton loops is

where the $`O_{Lijk}^{cns}`$ and $`O_{Rijk}^{cns}`$ couplings are given in Appendix B in Eqs. (B, B).

After approximating the rotation matrices $`𝒰`$, and $`𝒱`$ in the chargino sector, and $`𝒩`$ in the neutralino sector we find approximate expressions for these couplings that we will use below. These formulae are collected in Eqs. (69, 70, 71, B) of Appendix B.

### 2.4 Relevant Topologies

We now give the structure of the mass matrices relevant for the determination of solar neutrino masses and mixings. While in the BRpV model the atmospheric anomaly is explained at the tree-level, the solar neutrino masses and mixings are both generated radiatively. In particular, the “solar angle” has no meaning at the tree level due to the degeneracy of the two lightest neutrinos in this limit.

Diagonalizing the tree-level neutrino mass matrix first and adding then the 1-loop corrections before re-diagonalization the resulting neutrino/neutralino mass matrix has non-zero entries in the neutrino/neutrino, the neutrino/neutralino and in the neutralino/neutralino sectors. We have found that the most important part of the 1-loop neutrino masses derives from the neutrino/neutrino sector and that the 1-loop induced neutrino/neutralino mixing is usually subdominant.

The relevant topologies for the one loop calculation of neutrino masses are then illustrated in Fig.1.

Here our conventions are as follows: open circles with a cross inside indicate genuine mass insertions which flip chirality. On the other hand open circles without a cross correspond to small R-Parity violating projections, indicating how much of an Rp-even/odd mass eigenstate is present in a given Rp-odd/even weak eigenstate. Strictly speaking these projections are really coupling matrices attached to the vertices, and this is what appears in the numerical code. However, given the smallness of Rp-violating effects, the “insertion-method” proves to be a rather useful tool to develop an analytical perturbative expansion and to acquire some simple understanding of the results.

### 2.5 Bottom-sbottom loops

The simplest contribution to the radiatively induced neutrino mass arises from loops involving bottom quarks and squarks is given by Hirsch:2000ef

$$\stackrel{~}{\mathrm{\Pi }}_{ij}(0)=\frac{N_c}{16\pi ^2}\underset{r}{}\left(O_{Rjr}^{bn\stackrel{~}{b}}O_{Lir}^{bn\stackrel{~}{b}}+O_{Ljr}^{bn\stackrel{~}{b}}O_{Rir}^{bn\stackrel{~}{b}}\right)m_bB_0(0,m_b^2,m_r^2)$$
(28)
$`B_0(0,m_b^2,m_r^2)`$ is the usual Passarino-Veltman function Passarino:1978jh , FFLP . This contribution can be expressed as being proportional to the difference of two $`B_0`$ functions,

$$\mathrm{\Delta }B_0^{\stackrel{~}{b}_1\stackrel{~}{b}_2}=B_0(0,m_b^2,m_{\stackrel{~}{b}_1}^2)B_0(0,m_b^2,m_{\stackrel{~}{b}_2}^2)$$
(29)
as follows

$$\mathrm{\Delta }\stackrel{~}{\mathrm{\Pi }}_{ij}=\frac{N_cm_b}{16\pi ^2}2s_{\stackrel{~}{b}}c_{\stackrel{~}{b}}h_b^2\mathrm{\Delta }B_0^{\stackrel{~}{b}_1\stackrel{~}{b}_2}\left[\frac{\stackrel{~}{ϵ}_i\stackrel{~}{ϵ}_j}{\mu ^2}+a_3b\left(\stackrel{~}{ϵ}_i\delta _{j3}+\stackrel{~}{ϵ}_j\delta _{i3}\right)|\stackrel{}{\mathrm{\Lambda }}|+\left(a_3^2+\frac{a_La_R}{h_b^2}\right)\delta _{i3}\delta _{j3}|\stackrel{}{\mathrm{\Lambda }}|^2\right]$$
(30)
where we have defined

$$a_R=\frac{g}{\sqrt{2}}\left(\frac{1}{3}t_Wa_1a_2\right),a_L=\frac{g}{\sqrt{2}}\frac{2}{3}t_Wa_1$$
(31)
The different contributions can be understood as coming from the graphs corresponding to the first topology of Fig. 1. They have been depicted in more detail in Fig. 2,

where we have adopted the following conventions: a) as before, open circles correspond to small R-parity violating projections, indicating how much of a weak eigenstate is present in a given mass eigenstate, (b) full circles correspond to R-parity conserving projections and (c) open circles with a cross inside indicate genuine mass insertions which flip chirality.

The open and full circles should really appear at the vertices since the particles propagating in the loop are the mass eigenstates. We have however separated them to better identify the origin of the various terms. There is another set of graphs analogous to the previous ones which corresponds to the heavy sbottom. They are obtained from the previous graphs making the replacement $`\stackrel{~}{b}_1\stackrel{~}{b}_2`$, $`s_{\stackrel{~}{b}}c_{\stackrel{~}{b}}`$ and $`c_{\stackrel{~}{b}}s_{\stackrel{~}{b}}`$. Note that for all contributions to the $`2\times 2`$ submatrix corresponding to the light neutrinos the divergence from $`B_0(0,m_b^2,m_{\stackrel{~}{b}_1}^2)`$ is canceled by the divergence from $`B_0(0,m_b^2,m_{\stackrel{~}{b}_2}^2)`$, making finite the contribution from bottom-sbottom loops to this submatrix, as it should be, since the mass is fully “calculable”.

### 2.6 Charged Scalar-Charged Fermion Loops

Another contribution to the radiatively induced neutrino mass comes from charged-scalar/charged-fermion loops, given as Hirsch:2000ef

$$\stackrel{~}{\mathrm{\Pi }}_{ij}(0)=\frac{1}{16\pi ^2}\underset{k,r}{}\left(O_{Rkjr}^{cns}O_{Lkir}^{cns}+O_{Lkjr}^{cns}O_{Rkir}^{cns}\right)m_kB_0(0,m_k^2,m_r^2)$$
(32)
The structure of the contribution from charged Higgs/slepton loops is substantially more complex than that of the bottom-sbottom loop considered above. It can be expressed as

$`\mathrm{\Delta }\stackrel{~}{\mathrm{\Pi }}_{ij}`$ $`=`$ $`{\displaystyle \frac{m_\tau }{16\pi ^2}}[C_{ij}^{\stackrel{~}{\tau }_2\stackrel{~}{\tau }_1}\mathrm{\Delta }B_0^{\stackrel{~}{\tau }_2\stackrel{~}{\tau }_1}+C_{ij}^{H^\pm \stackrel{~}{\tau }_1}\mathrm{\Delta }B_0^{H^\pm \stackrel{~}{\tau }_1}+C_{ij}^{H^\pm \stackrel{~}{\tau }_2}\mathrm{\Delta }B_0^{H^\pm \stackrel{~}{\tau }_2}`$ (33)
$`+C_{ij}^{H^\pm \stackrel{~}{L}_1}\mathrm{\Delta }B_0^{H^\pm L_1}+C_{ij}^{H^\pm \stackrel{~}{L}_2}\mathrm{\Delta }B_0^{H^\pm L_2}+C_{ij}^{G^\pm \stackrel{~}{L}_1}\mathrm{\Delta }B_0^{G^\pm L_1}`$
$`+C_{ij}^{G^\pm \stackrel{~}{L}_2}\mathrm{\Delta }B_0^{G^\pm L_2}+C_{ij}^{G^\pm \stackrel{~}{\tau }_1\stackrel{~}{\tau }_2}\mathrm{\Delta }B_0^{G^\pm \stackrel{~}{\tau }_1\stackrel{~}{\tau }_2}+C_{ij}^{G^\pm H^\pm \stackrel{~}{\tau }_1\stackrel{~}{\tau }_2}\mathrm{\Delta }B_0^{G^\pm H^\pm \stackrel{~}{\tau }_1\stackrel{~}{\tau }_2}`$
$`+(ij)]`$
where

$`\mathrm{\Delta }B_0^{XY}B_0(0,m_\tau ^2,m_X^2)B_0(0,m_\tau ^2,m_Y^2);X,Y=(G^\pm ,H^\pm ,L_1,L_2,\stackrel{~}{\tau }_1,\stackrel{~}{\tau }_2)`$
$`\mathrm{\Delta }B_0^{G^\pm \stackrel{~}{\tau }_1\stackrel{~}{\tau }_2}c_{\stackrel{~}{\tau }}^2B_0(0,m_\tau ^2,m_{\stackrel{~}{\tau }_1}^2)+s_{\stackrel{~}{\tau }}^2B_0(0,m_\tau ^2,m_{\stackrel{~}{\tau }_2}^2)B_0(0,m_\tau ^2,m_{G^\pm }^2)`$
$`\mathrm{\Delta }B_0^{G^\pm H^\pm \stackrel{~}{\tau }_1\stackrel{~}{\tau }_2}c_\beta ^2B_0(0,m_\tau ^2,m_{G^\pm }^2)+s_\beta ^2B_0(0,m_\tau ^2,m_{H^\pm }^2)`$ (34)
$`c_{\stackrel{~}{\tau }}^2B_0(0,m_\tau ^2,m_{\stackrel{~}{\tau }_1}^2)s_{\stackrel{~}{\tau }}^2B_0(0,m_\tau ^2,m_{\stackrel{~}{\tau }_2}^2)`$
and

$`C_{ij}^{\stackrel{~}{\tau }_2\stackrel{~}{\tau }_1}=s_{\stackrel{~}{\tau }}c_{\stackrel{~}{\tau }}\{\sqrt{2}g^{}a_1|\stackrel{}{\mathrm{\Lambda }}|[gV_{\nu ,j3}^Ta_1^L\mathrm{\Lambda }_3\frac{1}{\sqrt{2}}(ga_2+g^{}a_1)|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3}]\delta _{i3}`$
$`+h_\tau ^2(b\stackrel{~}{ϵ}_i+a_3|\stackrel{}{\mathrm{\Lambda }}|\delta _{i3}c_\beta \frac{v_3}{v}V_{\nu ,i3}^T)[b\stackrel{~}{ϵ}_j+a_3|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3}V_{\nu ,j3}^T(a_2^L\mathrm{\Lambda }_3+bϵ_3)]\}`$
$`C_{ij}^{H^\pm \stackrel{~}{\tau }_1}=s_\beta \mathrm{\Theta }_{HL_3}\{c_{\stackrel{~}{\tau }}h_\tau V_{\nu ,i3}^T[gV_{\nu ,j3}^Ta_1^L\mathrm{\Lambda }_3\frac{1}{\sqrt{2}}(ga_2+g^{}a_1)|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3}]`$
$`+s_{\stackrel{~}{\tau }}h_\tau ^2V_{\nu ,i3}^T[b\stackrel{~}{ϵ}_j+a_3|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3}V_{\nu ,j3}^T(a_2^L\mathrm{\Lambda }_3+bϵ_3)]\}`$
$`C_{ij}^{H^\pm \stackrel{~}{\tau }_2}=s_\beta \mathrm{\Theta }_{HR_3}\{s_{\stackrel{~}{\tau }}h_\tau V_{\nu ,i3}^T[gV_{\nu ,j3}^Ta_1^L\mathrm{\Lambda }_3\frac{1}{\sqrt{2}}(ga_2+g^{}a_1)|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3}]`$
$`c_{\stackrel{~}{\tau }}h_\tau ^2V_{\nu ,i3}^T[b\stackrel{~}{ϵ}_j+a_3|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3}V_{\nu ,j3}^T(a_2^L\mathrm{\Lambda }_3+bϵ_3)]\}`$
$`C_{ij}^{H^\pm L_1}=s_\beta \stackrel{~}{\mathrm{\Theta }}_{HL_1}h_\tau gV_{\nu ,i3}^TV_{\nu ,j1}^Ta_1^L\mathrm{\Lambda }_3`$
$`C_{ij}^{H^\pm L_2}=s_\beta \stackrel{~}{\mathrm{\Theta }}_{HL_2}h_\tau gV_{\nu ,i3}^TV_{\nu ,j2}^Ta_1^L\mathrm{\Lambda }_3`$
$`C_{ij}^{G^\pm L_1}=c_\beta {\displaystyle \frac{v_1}{v}}h_\tau gV_{\nu ,i3}^TV_{\nu ,j1}^Ta_1^L\mathrm{\Lambda }_3`$ (35)
$`C_{ij}^{G^\pm L_2}=c_\beta {\displaystyle \frac{v_2}{v}}h_\tau gV_{\nu ,i3}^TV_{\nu ,j2}^Ta_1^L\mathrm{\Lambda }_3`$
$`C_{ij}^{G^\pm \stackrel{~}{\tau }_1\stackrel{~}{\tau }_2}=c_\beta {\displaystyle \frac{v_3}{v}}h_\tau V_{\nu ,i3}^T\left[gV_{\nu ,j3}^Ta_1^L\mathrm{\Lambda }_3{\displaystyle \frac{1}{\sqrt{2}}}(ga_2+g^{}a_1)|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3}\right]`$
$`C_{ij}^{G^\pm H^\pm \stackrel{~}{\tau }_1\stackrel{~}{\tau }_2}=h_\tau g\stackrel{~}{ϵ}_iV_{\nu ,i3}^Tba_1^L\mathrm{\Lambda }_3`$
The result of Eq. (33) can be represented graphically for better understanding. The terms proportional to $`\mathrm{\Delta }B_0^{\stackrel{~}{\tau }_2\stackrel{~}{\tau }_1}`$ come from the graphs of Fig. 3. There is another set of four graphs corresponding to $`\stackrel{~}{\tau }_2`$. These are found after making the replacements $`\stackrel{~}{\tau }_1\stackrel{~}{\tau }_2`$, $`s_{\stackrel{~}{\tau }}c_{\stackrel{~}{\tau }}`$, and $`c_{\stackrel{~}{\tau }}s_{\stackrel{~}{\tau }}`$.

The diagrams in the first row are the ones that are equivalent to those in the bottom-sbottom loop. They have as characteristic feature the presence of two Rp-violating insertions (open circles) in the external legs. However, in contrast to the quark sector, R-parity violation can also appear in the charged internal lines running in the loops, since it occurs in the charged fermion sector. This explains the origin of the second row in Fig. 3. The presence of R-parity violating insertions in the internal lines of the second row in Fig. 3 correspond to the second topology in Fig. 1. The full diagrammatic explanation of the rest of the terms appearing in Eq. (33) is given in detail in Appendix C.

## 3 Analytical versus Numerical results

In this section we will check the accuracy of the approximation formulas given in Sec. 2.5 and 2.6. We do this by comparing the results obtained with their use with a full numerical calculation of the one-loop contributions to the neutrino mass, whose details can be found in Ref. Hirsch:2000ef .

As will be explained in more detail below, the relative importance of the various loops depends on the - currently unknown - supersymmetric parameters. In order to reduce the number of free parameters in the following we will adopt the minimal constrained supergravity (mSUGRA) version of the MSSM. As a rule of thumb it can be said that the bottom-sbottom loop usually gives the main contribution to the neutrino mass matrix when the neutralino is the LSP. On the other hand, if the scalar tau is the LSP, both bottom-sbottom and charged scalar loops are of approximately comparable magnitudes.

We have therefore constructed two different random scans over SUSY parameter space. Both sets start with the following rather generous parameter ranges: $`M_2`$ from \[0,1.2\] TeV, $`|\mu |`$ from \[0,2.5\] TeV, $`m_0`$ in the range \[0,1.0\] TeV, $`A_0/m_0`$ and $`B_0/m_0`$ \[-3,3\] and $`\mathrm{tan}\beta `$ \[2.5,10\]. All randomly generated points were subsequently tested for consistency with the minimization (tadpole) conditions of the Higgs potential, as well as for phenomenological constraints from supersymmetric particle searches. We then selected points in which a) the lightest neutralino is the LSP (called set “Ntrl” in the following) or b) at least one of the charged sleptons was the LSP (called set “Stau” in the following). Note that in the Stau set $`m_0<<M_2`$ and large $`\mu `$ values are strongly preferred.

R-parity violating parameters are chosen in such a way that neutrino oscillation data are reproduced approximately. As discussed in the introduction, atmospheric neutrino experiments require a near-to-maximal atmospheric mixing angle $`\theta atm`$, with $`\mathrm{\Delta }m^2atm`$ in the range given in Eq. (5). On the other hand reactor data constrain the electron-neutrino component in the third mass eigenstate to be small. And, finally, in combination with solar neutrino data, the KamLAND data require a $`\theta sol`$ in the range given in Eq. (2) with $`\mathrm{\Delta }m^2sol`$ as given in Eq. (3). The latter ranges belong to the LMA-MSW region indicated by a solar-only global analysis of neutrino data given in Ref. Maltoni:2002ni . For completeness we also include the (pre-KamLAND) LOW and VAC-type solutions of the solar neutrino anomaly. In the following we will first discuss the bottom-sbottom and the charged scalar loops separately, before considering a calculation taking into account both loops in comparison to the full calculation.

### 3.1 Bottom-sbottom loop

In Fig. 4 we show the ratio of the approximate-over-exact solar neutrino mass parameter $`m_{\nu _2}^{Appr}/m_{\nu _2}^{exact}`$ versus $`\mathrm{\Delta }m^2sol`$ for the case in which only the bottom-sbottom loop is taken into account, both in the approximate and in the exact calculation. The horizontal bands indicate attainable neutrino mass values when the parameters are scanned as indicated previously. As can be seen from the figure the approximate formula works quite well for points in both Ntrl and Stau sets, as long as the neutrino masses fall in the LMA-MSW range indicated by the right vertical bands. Note that the LMA-MSW and LOW bands indicated in the figure correspond to the full analysis of solar data only, presented in Ref. Maltoni:2002ni . The recent KamLAND reactor neutrino data rule out the LOW solution and restricts the LMA-MSW to somewhat narrower ranges indicated in Eq. (3). One finds that the mass values inferred from our present analytical approximation are always within 10 % or less of the exact numerical calculation of the bottom-sbottom loop. Larger deviations show up only in the Ntrl set, for very small neutrino masses, which we trace to the neglection of the 1-loop neutrino/neutralino mixing terms in our approximate treatment. Although not strictly ruled out by a solar-only global neutrino data analysis Maltoni:2002ni , these LOW and VAC-type solutions are now strongly disfavored by the latest KamLAND reactor neutrino data.

### 3.2 Charged scalar loop

In Fig. 5 we show the ratio of the approximate-over-exact solar neutrino mass parameter $`m_{\nu _2}^{Appr}/m_{\nu _2}^{exact}`$ plotted versus $`\mathrm{\Delta }m^2sol`$, for a calculation which takes into account only the charged scalar loop in both the approximate and the exact calculation. As can be seen from the figure the approximate formula is accurate for all points in the LMA-MSW region, indicated by the right vertical bands Maltoni:2002ni , both for the Stau and for the Ntrl sets. The only case where our analytic results gives a poorer approximation (to better than a factor-of-2) of the full numerical result is for the Ntrl set, when the neutrino mass falls in the LOW or VAC ranges, now strongly disfavored by the KamLAND results. We have checked numerically that for these very small neutrino masses all terms in eq. (33) are of approximately equal importance and there are significant cancellations among terms, which leads to a less reliable final result.

### 3.3 Comparison with full calculation

In supersymmetric models with mSUGRA-like boundary conditions the bottom-sbottom and the charged scalar loop usually give the most important contribution to the neutrino mass matrix. This is demonstrated in Fig. 6 (left) for the set Ntrl and in Fig. 6 (right) for the set Stau. In both figures we show the ratio of the approximate-over-exact solar neutrino mass parameter $`m_{\nu _2}^{Appr}/m_{\nu _2}^{exact}`$ versus $`\mathrm{\Delta }m^2sol`$ in eV2, where $`m_{\nu _2}^{Appr}`$ is the approximate loop calculation involving the bottom-sbottom and the charged scalar loop, while $`m_{\nu _2}^{exact}`$ is the exact numerical computation taking into account all loops.

In the region of $`\mathrm{\Delta }m^2sol`$ appropriate for the currently preferred LMA-MSW solution to the solar neutrino problem one finds that the approximate calculation reproduces the exact result better than 10 %. Only in the set Ntrl one finds larger deviations, up to 25 % in extreme cases, when $`\mathrm{\Delta }m^2sol`$ lies in the LOW region, strongly disfavored by KamLAND. This is due to the larger errors in the bottom-sbottom calculation in this set for small neutrino masses as discussed above.

## 4 Simplified approximation formulas

### 4.1 The solar mass

First we note that for nearly all points in our random sets we find that $`m_{\nu _2}m_{\nu _3}`$. In other words, bilinear R-parity breaking favors a hierarchical neutrino spectrum. Moreover, we have found numerically that the terms proportional to $`\stackrel{~}{ϵ}_i\times \stackrel{~}{ϵ}_j`$ in the self energies in Eq. (30) give the most important contribution to $`m_{\nu _2}`$ in the bottom-sbottom loop calculation in most points of our sets. If these terms are dominant one can find a very simple approximation for the bottom-sbottom loop contribution to $`m_{\nu _2}`$. It is given by

$$m_{\nu _2}\frac{3}{16\pi ^2}\mathrm{sin}(2\theta _{\stackrel{~}{b}})m_b\mathrm{\Delta }B_0^{\stackrel{~}{\tau }_2\stackrel{~}{\tau }_1}\frac{(\stackrel{~}{ϵ}_1^2+\stackrel{~}{ϵ}_2^2)}{\mu ^2}$$
(36)
We have checked numerically that Eq. (36) reproduces the result of the full approximative formula to high accuracy if $`m_{\nu _2}0.3m_{\nu _3}`$ . Note also that Eq. (36) holds only if the 1-loop contributions to the neutrino mass matrix are smaller than the tree-level one. This condition requires that $`|\stackrel{}{ϵ}|^2/|\mathrm{\Lambda }|1`$ approximately, i.e. the bilinear parameters $`ϵ_i`$ must be suppressed with respect to $`\mu `$. Note that such a suppression could, in principle, be motivated by suitable flavour symmetries Mira:2000gg .

Due to the more complicated structure of the charged scalar loop it is not possible to give a simple equation for $`m_{\nu _2}`$ similar to Eq. (36) for the bottom-sbottom loop. However, for $`m_{\nu _2}`$ larger than (few) $`\times 10^4`$ we have found that the most important contributions to the charged scalar loop are the terms proportional to $`\mathrm{\Delta }B_0^{\stackrel{~}{\tau }_2\stackrel{~}{\tau }_1}`$, $`\mathrm{\Delta }B_0^{H^\pm \stackrel{~}{\tau }_1}`$ and $`\mathrm{\Delta }B_0^{H^\pm \stackrel{~}{\tau }_2}`$ in Eq. (33). We note in passing that Eq. (36), with appropriate replacements, allows us to estimate the typical contributions to the charged scalar loop within a factor of $`3`$. However, such an estimate will be biased toward too small (large) $`m_{\nu _2}`$ for scalar tau (neutralino) LSPs.

In Fig. (7) we show a comparison of our simplified approximation formula, including the simple form of the sbottom-bottom loop and the three most important coefficients for the charged scalar loop, as discussed above, to the full numerical calculation including all loops. As one can see, even the simplified version of our formula works surprisingly well in the LMA-MSW regime, although the agreement with the full calculation is now less good for the LOW region, as could have been expected from the results discussed previously.

### 4.2 The solar mixing angle

In the basis where the tree-level neutrino mass matrix is diagonal the mass matrix at one–loop level can be written as

$`\stackrel{~}{m}_\nu =V_\nu ^{(0)T}m_\nu V_\nu ^{(0)}=\left(\begin{array}{ccc}c_1\stackrel{~}{ϵ}_1\stackrel{~}{ϵ}_1& c_1\stackrel{~}{ϵ}_1\stackrel{~}{ϵ}_2& c_1\stackrel{~}{ϵ}_1\stackrel{~}{ϵ}_3\\ c_1\stackrel{~}{ϵ}_2\stackrel{~}{ϵ}_1& c_1\stackrel{~}{ϵ}_2\stackrel{~}{ϵ}_2& c_1\stackrel{~}{ϵ}_2\stackrel{~}{ϵ}_3\\ c_1\stackrel{~}{ϵ}_3\stackrel{~}{ϵ}_1& c_1\stackrel{~}{ϵ}_3\stackrel{~}{ϵ}_2& c_0|\stackrel{}{\mathrm{\Lambda }}|^2+c_1\stackrel{~}{ϵ}_3\stackrel{~}{ϵ}_3\end{array}\right)+\mathrm{}`$ (40)
where the $`\stackrel{~}{ϵ}_i`$ were defined before in Eq. (2.2). Coefficients $`c_0`$ and $`c_1`$ contain couplings and supersymmetric masses. Since they cancel in the final expression for the angle their exact definition is not necessary in the following. Dots stand for other terms which we will assume to be less important in the following, see the discussion at the end of this subsection. This matrix can be diagonalized approximately under the condition

$`x{\displaystyle \frac{c_1|\stackrel{}{\stackrel{~}{ϵ}}|^2}{c_0|\stackrel{}{\mathrm{\Lambda }}|^2}}1`$ (41)
i.e. if the 1-loop contribution to the neutrino mass matrix is smaller than the tree-level contribution, as also discussed above for Eq. (36). Then

$`\stackrel{~}{m}_\nu =c_0|\stackrel{}{\mathrm{\Lambda }}|^2\left(\begin{array}{ccc}x{\displaystyle \frac{\stackrel{~}{ϵ}_1\stackrel{~}{ϵ}_1}{|\stackrel{}{\stackrel{~}{ϵ}|^2}}}& x{\displaystyle \frac{\stackrel{~}{ϵ}_1\stackrel{~}{ϵ}_2}{|\stackrel{}{\stackrel{~}{ϵ}|^2}}}& x{\displaystyle \frac{\stackrel{~}{ϵ}_1\stackrel{~}{ϵ}_3}{|\stackrel{}{\stackrel{~}{ϵ}|^2}}}\\ x{\displaystyle \frac{\stackrel{~}{ϵ}_2\stackrel{~}{ϵ}_1}{|\stackrel{}{\stackrel{~}{ϵ}|^2}}}& x{\displaystyle \frac{\stackrel{~}{ϵ}_2\stackrel{~}{ϵ}_2}{|\stackrel{}{\stackrel{~}{ϵ}|^2}}}& x{\displaystyle \frac{\stackrel{~}{ϵ}_2\stackrel{~}{ϵ}_3}{|\stackrel{}{\stackrel{~}{ϵ}|^2}}}\\ x{\displaystyle \frac{\stackrel{~}{ϵ}_3\stackrel{~}{ϵ}_1}{|\stackrel{}{\stackrel{~}{ϵ}|^2}}}& x{\displaystyle \frac{\stackrel{~}{ϵ}_3\stackrel{~}{ϵ}_2}{|\stackrel{}{\stackrel{~}{ϵ}|^2}}}& 1+x{\displaystyle \frac{\stackrel{~}{ϵ}_3\stackrel{~}{ϵ}_3}{|\stackrel{}{\stackrel{~}{ϵ}|^2}}}\end{array}\right)`$ (45)
We now calculate the eigenvalues and eigenvectors of this matrix as series expansions in the small $`x`$ parameter. For the eigenvalues we get

$`m_1`$ $`=`$ $`0`$
$`m_2`$ $`=`$ $`xc_0{\displaystyle \frac{|\stackrel{}{\mathrm{\Lambda }}|^2}{|\stackrel{}{\stackrel{~}{ϵ}|^2}}}+𝒪(x^2)=c_1\left(\stackrel{~}{ϵ}_1^2+\stackrel{~}{ϵ}_2^2\right)+𝒪(x^2)`$ (46)
$`m_3`$ $`=`$ $`c_0|\stackrel{}{\mathrm{\Lambda }}|^2+c_1\stackrel{~}{ϵ}_3^2+𝒪(x^2)`$
and for the first two eigenvalues (the third can also be easily obtained but it will not be necessary for the discussion of the solar mixing angle),

$`e_1`$ $`=`$ $`({\displaystyle \frac{\stackrel{~}{ϵ}_2}{\stackrel{~}{ϵ}_1}}\sqrt{{\displaystyle \frac{\stackrel{~}{ϵ}_1^2}{\stackrel{~}{ϵ}_1^2+\stackrel{~}{ϵ}_2^2}}},\sqrt{{\displaystyle \frac{\stackrel{~}{ϵ}_1^2}{\stackrel{~}{ϵ}_1^2+\stackrel{~}{ϵ}_2^2}}},0)`$
$`e_2`$ $`=`$ $`(e_{2,1},e_{2,2},e_{2,3})`$ (47)
where up to $`𝒪(x^2)`$ we have,

$`e_{2,1}`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{ϵ}_1\stackrel{~}{ϵ}_3}{\sqrt{\stackrel{~}{ϵ}_3^2(\stackrel{~}{ϵ}_1^2+\stackrel{~}{ϵ}_2^2)}}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\stackrel{~}{ϵ}_1\stackrel{~}{ϵ}_3\sqrt{\stackrel{~}{ϵ}_3^2(\stackrel{~}{ϵ}_1^2+\stackrel{~}{ϵ}_2^2)}}{|\stackrel{}{\stackrel{~}{ϵ}|^4}}}x^2+𝒪(x^3)`$
$`e_{2,2}`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{ϵ}_2\stackrel{~}{ϵ}_3}{\sqrt{\stackrel{~}{ϵ}_3^2(\stackrel{~}{ϵ}_1^2+\stackrel{~}{ϵ}_2^2)}}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\stackrel{~}{ϵ}_2\stackrel{~}{ϵ}_3\sqrt{\stackrel{~}{ϵ}_3^2(\stackrel{~}{ϵ}_1^2+\stackrel{~}{ϵ}_2^2)}}{|\stackrel{}{\stackrel{~}{ϵ}|^4}}}x^2+𝒪(x^3)`$ (48)
$`e_{2,3}`$ $`=`$ $`{\displaystyle \frac{\sqrt{\stackrel{~}{ϵ}_3^2(\stackrel{~}{ϵ}_1^2+\stackrel{~}{ϵ}_2^2)}}{|\stackrel{}{\stackrel{~}{ϵ}|^2}}}x+{\displaystyle \frac{\left(\stackrel{~}{ϵ}_1^2+\stackrel{~}{ϵ}_2^2\stackrel{~}{ϵ}_3^2\right)\sqrt{\stackrel{~}{ϵ}_3^2(\stackrel{~}{ϵ}_1^2+\stackrel{~}{ϵ}_2^2)}}{|\stackrel{}{\stackrel{~}{ϵ}|^4}}}x^2+𝒪(x^3)`$
Knowing the eigenvectors we can write down the rotation matrix that diagonalizes $`\stackrel{~}{m}_\nu `$,

$`\stackrel{~}{V}_\nu ^T\stackrel{~}{m}_\nu \stackrel{~}{V}_\nu =\text{diag}(m_1,m_2,m_3)`$ (49)
where

$`\stackrel{~}{V}_\nu ^T=\left(\begin{array}{ccc}e_{1,1}& e_{1,2}& e_{1,3}\\ e_{2,1}& e_{2,2}& e_{2,3}\\ e_{3,1}& e_{3,2}& e_{3,3}\end{array}\right)`$ (53)
The neutrino mixing matrix is then given by

$`U=\left(V_\nu ^T\stackrel{~}{V}_\nu ^T\right)^T`$ (54)
Using the fact that $`U_{e3}`$ has to be small one can get the following expression for the solar mixing angle:

$`\mathrm{tan}^2\theta sol={\displaystyle \frac{U_{e2}^2}{U_{e1}^2}}`$ (55)
Now using the Eqs. (53), (4.2) and substituting in Eq. (54) we obtain the very simple expression for the solar mixing angle,

$`\mathrm{tan}^2\theta sol={\displaystyle \frac{\stackrel{~}{ϵ}_1^2}{\stackrel{~}{ϵ}_2^2}}`$ (56)
This formula is a very good approximation if the one–loop matrix has the structure $`ϵ_i\times ϵ_j`$, as is the case of the bottom-sbottom loop (and, to a lesser extend also for the charged scalar loop, which has one coefficient with the same index structure), and if $`m_{\nu _3}m_{\nu _2}`$ . This is illustrated in Fig. 8.

In the left panel we show a calculation comparing for all points in the set Ntrl the approximate to the exact solar angle, while the right panel shows a subset of points using the cut $`\mathrm{sin}(2\theta _{\stackrel{~}{b}})\mathrm{\Delta }B_0^{\stackrel{~}{\tau }_2\stackrel{~}{\tau }_1}>0.02`$. Note that this cut is designed such as to prefer points in which there is a sizeable contribution to the full 1-loop neutrino mass due to the bottom-sbottom loop. For points in which the charged scalar loop dominates eq. (56) gives only a factor-of-two estimate of the true solar angle.

Note finally that eq. (56) will fail completely, if $`\mathrm{\Lambda }_\mu \mathrm{\Lambda }_\tau `$ and $`ϵ_\mu ϵ_\tau `$, since then $`\stackrel{~}{ϵ}_2^2=0`$, see Eq. (2.2). This is the origin of the “sign condition” discussed in Hirsch:2000ef .

## 5 Discussion and conclusions

We have presented an approximate calculation of the neutrino mass matrix at one-loop in supersymmetry with bilinearly broken R-parity. The method is based on a systematic perturbative expansion of R-parity violating vertices to leading order. We have identified the bottom-sbottom and the charged scalar loop as the most important ones, at least in supersymmetric models with mSUGRA-like boundary conditions. Taking into account only these loops, we have given explicit formulas and discussed their validity as well as the accuracy with which they describe solar neutrino mass and mixing parameters. This was done by comparing our analytical results to the exact numerical calculation. We have found that for the case of the large mixing MSW solution our formulas - even within the simplified form Eq. (36) and Eq. (56) - yield good agreement with the full numerical calculation, but are much simpler to implement than the full numerical one-loop calculation. The only solar neutrino “solutions” for which our analytical approximation is less accurate are those that are now ruled out by the recent reactor neutrino data from KamLAND.

Let us finally discuss some possible caveats to the success of our approximate treatment. One is the assumption that supersymmetry breaking mass terms are flavour diagonal, which we have adopted, motivated by constraints from flavour changing processes. Although such terms could be included into our approximate treatment, we have not done so, mainly due to the fact that the resulting formulas would be much more complicated and, therefore, of very limited practical use. A second concern is that our sample points were all generated using mSugra assumptions for the soft breaking masses. Clearly there are other possibilities to break supersymmetry and even though we expect that the bottom-sbottom loop and the charged scalar loop will still be well described by our approximation formulas, other loops, which we didn’t take into account, might be more important than what we have found in our data sets.

## Acknowledgements

This work was supported by Spanish grant BFM2002-00345, by the European Commission RTN grant HPRN-CT-2000-00148 and the ESF *Neutrino Astrophysics Network*. M. H. is supported by a Spanish MCyT Ramon y Cajal contract. W. P. is supported by the ’Erwin Schrödinger fellowship No. J2095’ of the ‘Fonds zur Förderung der wissenschaftlichen Forschung’ of Austria FWF and partly by the Swiss ‘Nationalfonds’. M.A.D. was supported by Conicyt grant No. 1010974. J.R. was partly supported by the fellowship SFRH/BSAB/269/2002 from FCT (Portugal) with funds from the EU.

## Appendix A Rotation Matrices

In the basis $`(H_d^+,H_u^+,\stackrel{~}{e}_L^+,\stackrel{~}{\mu }_L^+,\stackrel{~}{\tau }_L^+,\stackrel{~}{e}_R^+,\stackrel{~}{\mu }_R^+,\stackrel{~}{\tau }_R^+)`$, one can write, to first order in R-parity violating parameters, the Goldstone rotation matrix as

$$𝐑_G=\left[\begin{array}{cccccccc}c_\beta & s_\beta & v_1/v& v_2/v& v_3/v& 0& 0& 0\\ s_\beta & c_\beta & 0& 0& 0& 0& 0& 0\\ c_\beta v_1/v& s_\beta v_1/v& 1& 0& 0& 0& 0& 0\\ c_\beta v_2/v& s_\beta v_2/v& 0& 1& 0& 0& 0& 0\\ c_\beta v_3/v& s_\beta v_3/v& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]$$
(57)
where $`v^2=v_d^2+v_u^2`$ to this order, and $`\mathrm{tan}\beta =v_u/v_d`$, as usual. We have also used the shorthand notation $`c_\beta (s_\beta )=\mathrm{cos}\beta (\mathrm{sin}\beta )`$.

Neglecting the electron and muon Yukawa couplings, the rotation that diagonalizes the sleptons at tree level is given by (in the same basis as above)

$$𝐑_{\stackrel{~}{\tau }}=\left[\begin{array}{cccccccc}1& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& c_{\stackrel{~}{\tau }}& 0& 0& s_{\stackrel{~}{\tau }}\\ 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& s_{\stackrel{~}{\tau }}& 0& 0& c_{\stackrel{~}{\tau }}\end{array}\right]$$
(58)
After the rotations $`𝐑_{\stackrel{~}{\tau }}𝐑_G`$ are performed, the charged scalar mass matrix is diagonalized up to small R-parity violating entries. In the approximation where there is no intergenerational mixing and $`h_\mu h_e0`$, these are

$$𝚫𝐌_{𝐒^\pm }^{\mathrm{𝟐}}{}_{\stackrel{~}{\tau }}{}^{}=\left[\begin{array}{cccccccc}0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& \stackrel{~}{X}_{HL_1}& \stackrel{~}{X}_{HL_2}& \stackrel{~}{X}_{HL_3}& 0& 0& \stackrel{~}{X}_{HR_3}\\ 0& \stackrel{~}{X}_{HL_1}& 0& 0& 0& 0& 0& 0\\ 0& \stackrel{~}{X}_{HL_2}& 0& 0& 0& 0& 0& 0\\ 0& \stackrel{~}{X}_{HL_3}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& \stackrel{~}{X}_{HR_3}& 0& 0& 0& 0& 0& 0\end{array}\right]$$
(59)
where

| $`\stackrel{~}{X}_{HL_i}=X_{HL_i}`$ | $`\stackrel{~}{X}_{HR_i}=X_{HR_i},(i=1,2)`$ |
| --- | --- |
| $`\stackrel{~}{X}_{HL_3}=c_{\stackrel{~}{\tau }}X_{uL_3}+s_{\stackrel{~}{\tau }}X_{dR_3}`$ | $`\stackrel{~}{X}_{HR_3}=s_{\stackrel{~}{\tau }}X_{HL_3}+c_{\stackrel{~}{\tau }}X_{HR_3}`$ |

(60)
with

| $`X_{HL_i}=s_\beta X_{uL_i}+c_\beta X_{dL_i}`$ | $`X_{HR_3}=s_\beta X_{uR_3}+c_\beta X_{dR_3}`$, (i=1,3) |
| --- | --- |

(61)
and

| $`X_{uL_i}=\frac{1}{4}g^2v_dv_i\mu ϵ_i\frac{1}{2}h_\tau ^2v_dv_i\delta _{i3},`$ | $`X_{dL_i}=\frac{v_i}{v_d}\frac{c_\beta }{s_\beta }m_{\stackrel{~}{\nu }}^2\mu ϵ_i\frac{c_\beta }{s_\beta }+\frac{1}{4}g^2v_uv_i,`$ |
| --- | --- |
| $`X_{uR_3}=\frac{1}{\sqrt{2}}h_\tau (A_\tau v_3+ϵ_3v_u),`$ | $`X_{dR_3}=\frac{1}{\sqrt{2}}h_\tau (\mu v_3+ϵ_3v_d).`$ |

(62)
These mixings are removed with the rotation matrix $`𝐑_𝐗`$ given by

$$𝐑_𝐗=\left[\begin{array}{cccccccc}1& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& \mathrm{\Theta }_{HL_1}& \mathrm{\Theta }_{HL_2}& \mathrm{\Theta }_{HL_3}& 0& 0& \mathrm{\Theta }_{HR_3}\\ 0& \mathrm{\Theta }_{HL_1}& 1& 0& 0& 0& 0& 0\\ 0& \mathrm{\Theta }_{HL_2}& 0& 1& 0& 0& 0& 0\\ 0& \mathrm{\Theta }_{HL_3}& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& \mathrm{\Theta }_{HR_3}& 0& 0& 0& 0& 0& 1\end{array}\right]$$
(63)
in the small mixing approximation $`\mathrm{sin}\mathrm{\Theta }\mathrm{\Theta }`$. Note that here we have defined

$$\mathrm{\Theta }_{HL_i}\frac{\stackrel{~}{X}_{HL_i}}{m_{H^\pm }^2m_{\stackrel{~}{\mathrm{}}_{Li}}^2},\mathrm{\Theta }_{HR_i}\frac{\stackrel{~}{X}_{HR_i}}{m_{H^\pm }^2m_{\stackrel{~}{\mathrm{}}_{Ri}}^2}$$
(64)
Putting everything together we get the final form of the charged scalar diagonalization matrix $`𝐑_X𝐑_{\stackrel{~}{\tau }}𝐑_G`$ which can be expressed as

$`𝐑_X𝐑_{\stackrel{~}{\tau }}𝐑_G`$ $`=`$
$`\left[\begin{array}{cccccccc}c_\beta & s_\beta & v_1/v& v_2/v& v_3/v& 0& 0& 0\\ s_\beta & c_\beta & \mathrm{\Theta }_{HL_1}& \mathrm{\Theta }_{HL_2}& \stackrel{~}{\mathrm{\Theta }}_{HL_3}& 0& 0& \stackrel{~}{\mathrm{\Theta }}_{HR_3}\\ s_\beta \mathrm{\Theta }_{HL_1}c_\beta \frac{v_1}{v}& c_\beta \mathrm{\Theta }_{HL_1}+s_\beta \frac{v_1}{v}& 1& 0& 0& 0& 0& 0\\ s_\beta \mathrm{\Theta }_{HL_2}c_\beta \frac{v_2}{v}& c_\beta \mathrm{\Theta }_{HL_2}+s_\beta \frac{v_2}{v}& 0& 1& 0& 0& 0& 0\\ s_\beta \mathrm{\Theta }_{HL_3}c_{\stackrel{~}{\tau }}c_\beta \frac{v_3}{v}& c_\beta \mathrm{\Theta }_{HL_3}+c_{\stackrel{~}{\tau }}s_\beta \frac{v_3}{v}& 0& 0& c_{\stackrel{~}{\tau }}& 0& 0& s_{\stackrel{~}{\tau }}\\ 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ s_\beta \mathrm{\Theta }_{HR_3}+s_{\stackrel{~}{\tau }}c_\beta \frac{v_3}{v}& c_\beta \mathrm{\Theta }_{HR_3}s_{\stackrel{~}{\tau }}s_\beta \frac{v_3}{v}& 0& 0& s_{\stackrel{~}{\tau }}& 0& 0& c_{\stackrel{~}{\tau }}\end{array}\right]`$
where we have defined,

$$\stackrel{~}{\mathrm{\Theta }}_{HL_3}=c_{\stackrel{~}{\tau }}\mathrm{\Theta }_{HL_3}s_{\stackrel{~}{\tau }}\mathrm{\Theta }_{HR_3},\stackrel{~}{\mathrm{\Theta }}_{HR_3}=s_{\stackrel{~}{\tau }}\mathrm{\Theta }_{HL_3}+c_{\stackrel{~}{\tau }}\mathrm{\Theta }_{HR_3}.$$
(66)
## Appendix B Charged Higgs/slepton couplings

The couplings of the five (generalized to include also the three charged leptons) charginos to the eight charged scalars (including Higgs bosons and sleptons of both chiralities) and seven neutralinos (generalized to include also the three neutrinos) are given by Hirsch:2000ef

$`O_{Lijk}^{cns}`$ $`=`$ $`R_{k1}^{S^\pm }h_\tau 𝒩_{j7}𝒱_{i5}R_{k2}^{S^\pm }\left({\displaystyle \frac{g}{\sqrt{2}}}𝒩_{j2}𝒱_{i2}+{\displaystyle \frac{g^{}}{\sqrt{2}}}𝒩_{j1}𝒱_{i2}+g𝒩_{j4}𝒱_{i1}\right)`$
$`R_{k5}^{S^\pm }h_\tau 𝒩_{j3}𝒱_{i5}g^{}\sqrt{2}\left(R_{k6}^{S^\pm }𝒩_{j1}𝒱_{i3}+R_{k7}^{S^\pm }𝒩_{j1}𝒱_{i4}+R_{k8}^{S^\pm }𝒩_{j1}𝒱_{i5}\right)`$
where $`i`$ labels the charginos, $`j`$ labels neutralinos, and $`k`$ labels the charged scalars, respectively. For the the right-handed couplings the corresponding couplings are given by

$`O_{Rijk}^{cns}`$ $`=`$ $`R_{k1}^{S^\pm }\left({\displaystyle \frac{g}{\sqrt{2}}}𝒩_{j2}𝒰_{i2}+{\displaystyle \frac{g^{}}{\sqrt{2}}}𝒩_{j1}𝒰_{i2}g𝒩_{j3}𝒰_{i1}\right)`$
$`+R_{k3}^{S^\pm }\left({\displaystyle \frac{g}{\sqrt{2}}}𝒩_{j2}𝒰_{i3}+{\displaystyle \frac{g^{}}{\sqrt{2}}}𝒩_{j1}𝒰_{i3}g𝒩_{j5}𝒰_{i1}\right)`$
$`+R_{k4}^{S^\pm }\left({\displaystyle \frac{g}{\sqrt{2}}}𝒩_{j2}𝒰_{i4}+{\displaystyle \frac{g^{}}{\sqrt{2}}}𝒩_{j1}𝒰_{i4}g𝒩_{j6}𝒰_{i1}\right)`$
$`+R_{k5}^{S^\pm }\left({\displaystyle \frac{g}{\sqrt{2}}}𝒩_{j2}𝒰_{i5}+{\displaystyle \frac{g^{}}{\sqrt{2}}}𝒩_{j1}𝒰_{i5}g𝒩_{j7}𝒰_{i1}\right)`$
$`+R_{k8}^{S^\pm }h_\tau \left(𝒩_{j7}𝒰_{i2}𝒩_{j3}𝒰_{i5}\right)`$
After approximating the rotation matrices $`𝒰`$, and $`𝒱`$ in the chargino sector, and $`𝒩`$ in the neutralino sector we find the expressions given in Eqs. (69)-(B). Note that we have divided them into cases where the charged fermion is a lepton or a chargino. For the left couplings when the charged fermion is a chargino we have,

$$O_{Lijk}^{cns}=R_{k2}^{S^\pm }\left[\frac{g}{\sqrt{2}}a_2V_{i^{}2}+\frac{g^{}}{\sqrt{2}}a_1V_{i^{}2}+ga_4V_{i^{}1}\right]|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3}$$
(69)
where $`V`$ is the reduced $`2\times 2`$ chargino diagonalization matrix of the MSSM, and $`i^{}=1,2`$. If the charged fermion is a lepton we have

$`O_{Lijk}^{cns}`$ $`=`$ $`R_{k1}^{S^\pm }h_\tau V_{\nu ,j3}^T\delta _{i3}+R_{k5}^{S^\pm }h_\tau \left(b\stackrel{~}{ϵ}_j+a_3|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3}\right)\delta _{i3}`$ (70)
$`+\left[R_{k6}^{S^\pm }\delta _{i1}+R_{k7}^{S^\pm }\delta _{i2}+R_{k8}^{S^\pm }\delta _{i3}\right]\sqrt{2}g^{}a_1|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3}`$
For the right-handed couplings when the charged fermion is a chargino we get

$`O_{Rijk}^{cns}`$ $`=`$ $`R_{k1}^{S^\pm }\left[\frac{1}{\sqrt{2}}(ga_2+g^{}a_1)|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3}U_{i^{}2}+g(b\stackrel{~}{ϵ}_j+a_3|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3})U_{i^{}1}\right]`$ (71)
$`R_{k3}^{S^\pm }gV_{\nu ,j1}^TU_{i^{}1}R_{k4}^{S^\pm }gV_{\nu ,j2}^TU_{i^{}1}R_{k5}^{S^\pm }gV_{\nu ,j3}^TU_{i^{}1}+R_{k8}^{S^\pm }h_\tau V_{\nu ,j3}^TU_{i^{}2}`$
where $`U`$ is the second $`2\times 2`$ chargino rotation matrix of the MSSM. Finally, if the charged fermion is a lepton one has

$`O_{Rijk}^{cns}`$ $`=`$ $`R_{k3}^{S^\pm }\left[\frac{1}{\sqrt{2}}(ga_2+g^{}a_1)|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3}\delta _{i1}gV_{\nu ,j1}^Ta_1^L\mathrm{\Lambda }_i\right]`$
$`R_{k4}^{S^\pm }\left[\frac{1}{\sqrt{2}}(ga_2+g^{}a_1)|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3}\delta _{i2}gV_{\nu ,j2}^Ta_1^L\mathrm{\Lambda }_i\right]`$
$`R_{k5}^{S^\pm }\left[\frac{1}{\sqrt{2}}(ga_2+g^{}a_1)|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3}\delta _{i3}gV_{\nu ,j3}^Ta_1^L\mathrm{\Lambda }_i\right]`$
$`R_{k8}^{S^\pm }h_\tau \left[V_{\nu ,j3}^T(a_2^L\mathrm{\Lambda }_i+bϵ_i)(b\stackrel{~}{ϵ}_j+a_3|\stackrel{}{\mathrm{\Lambda }}|\delta _{j3})\delta _{i3}\right]`$
## Appendix C Charged Scalar–Charged Fermion Loops

There are nine different terms contributing to the Charged Scalar-Charged Fermion loop, as it was shown in Eq. (33). All these terms give a finite contribution to the $`2\times 2`$ submatrix corresponding to the light neutrinos. In this Appendix we will explain with graphs the origin of the different terms. The conventions used were explained in section 2.6.

### C.1 $`\mathrm{\Delta }B_0^{\stackrel{~}{\tau }_2\stackrel{~}{\tau }_1}`$

The terms proportional to $`\mathrm{\Delta }B_0^{\stackrel{~}{\tau }_2\stackrel{~}{\tau }_1}`$ come from the graphs of Fig. 3 as explained in section 2.6.

### C.2 $`\mathrm{\Delta }B_0^{H^+\stackrel{~}{\tau }_1}`$ and $`\mathrm{\Delta }B_0^{H^+\stackrel{~}{\tau }_2}`$

Now consider the terms proportional to $`\mathrm{\Delta }B_0^{H^+\stackrel{~}{\tau }_1}`$ and $`\mathrm{\Delta }B_0^{H^+\stackrel{~}{\tau }_2}`$ in Eq. (33). Of these terms, the ones which are related to the charged Higgs mixing with staus, can be understood as coming from the four graphs of Fig. 9. Associated to these charged Higgs graphs are those related to the $`\stackrel{~}{\tau }_1`$ mixing with charged Higgs. These are given in Fig. 10.

There is another set of four graphs corresponding to $`\stackrel{~}{\tau }_2`$ that are obtained from those in Fig. 10 by replacing $`\stackrel{~}{\tau }_1\stackrel{~}{\tau }_2`$, $`s_{\stackrel{~}{\tau }}c_{\stackrel{~}{\tau }}`$ and $`c_{\stackrel{~}{\tau }}s_{\stackrel{~}{\tau }}`$. These three groups of four graphs, when combined, form a set which is ultraviolet finite and account for the terms in Eq. (33) proportional to $`\mathrm{\Delta }B_0^{H^+\stackrel{~}{\tau }_1}`$ and $`\mathrm{\Delta }B_0^{H^+\stackrel{~}{\tau }_2}`$.

### C.3 $`\mathrm{\Delta }B_0^{H^\pm L_1}`$ and $`\mathrm{\Delta }B_0^{H^\pm L_2}`$

We now turn our attention to the terms proportional to $`\mathrm{\Delta }B_0^{H^\pm L_1}`$ and $`B_0^{H^\pm L_2}`$ which are related to the mixing between charged Higgs with selectrons and smuons. The terms proportional to $`\mathrm{\Delta }B_0^{H^\pm L_1}`$ come from the diagrams of Fig. 11. The terms proportional to $`\mathrm{\Delta }B_0^{H^\pm L_2}`$ are easily obtained from these by replacing the corresponding slepton lines and couplings. Notice that in Fig. 11 b) there is a contribution proportional to $`v_1/v`$ that does not belong to this term. We will show below that it will contribute to the $`\mathrm{\Delta }B_0^{G^\pm L_1}`$ term.

### C.4 $`\mathrm{\Delta }B_0^{G^\pm L_1}`$ and $`\mathrm{\Delta }B_0^{G^\pm L_2}`$

The graphs contributing to the $`\mathrm{\Delta }B_0^{G^\pm L_1}`$ and $`\mathrm{\Delta }B_0^{G^\pm L_2}`$ terms are related to those of Fig. 11. They are given by Fig. 12 and by the term proportional to $`v_1/v`$ in Fig. 11 b), for the case of the selectron. The terms proportional to $`\mathrm{\Delta }B_0^{G^\pm L_2}`$ are easily obtained from these by replacing the corresponding slepton lines and couplings.

### C.5 $`\mathrm{\Delta }B_0^{G^\pm \stackrel{~}{\tau }_1\stackrel{~}{\tau }_2}`$

We now consider a more complicated term, the one proportional to $`\mathrm{\Delta }B_0^{G^\pm \stackrel{~}{\tau }_1\stackrel{~}{\tau }_2}`$. This term gives a finite ultraviolet contribution and comes from the diagrams of Fig. 13, together with the parts of the diagrams of Fig. 10 that are proportional to $`v_3/v`$. Corresponding to the diagrams in Fig. 13 proportional to $`v_3/v`$, there is another set with $`\stackrel{~}{\tau }_1`$ and $`\stackrel{~}{\tau }_2`$ interchanged in the usual way.

### C.6 $`\mathrm{\Delta }B_0^{G^\pm H^\pm \stackrel{~}{\tau }_1\stackrel{~}{\tau }_2}`$

Let us consider finally the last term in Eq. (33), the one proportional to $`\mathrm{\Delta }B_0^{G^\pm H^\pm \stackrel{~}{\tau }_1\stackrel{~}{\tau }_2}`$. This term gives an ultraviolet finite contribution and comes from four diagrams. The first two are those represented in Fig. 14 corresponding to a $`H^\pm `$ and $`\stackrel{~}{\tau }_1`$ propagating in the loop. The other two are obtained from these with the replacements,

| $`H^\pm G^\pm `$ | , | $`s_\beta c_\beta `$ |
| --- | --- | --- |
| $`\stackrel{~}{\tau }_1\stackrel{~}{\tau }_2`$ | , | $`s_{\stackrel{~}{\tau }}c_{\stackrel{~}{\tau }}`$ |

(73)


format: text

text:
# The variability of accretion onto Schwarzschild black holes from turbulent magnetized discs

## 1 Introduction

A major aim of the astrophysical study of black holes is to determine the fundamental parameters of a rotating black hole – the mass $`M`$ and spin parameter $`a`$ – from astronomical observations. For Galactic black hole candidates such as Cygnus X-1, GRO J1655-40 and GRS 1915+105, the temporal power spectrum holds promise as a route to achieving this goal (van der Klis 2000). The variability of these sources is conventionally described as the superposition of a broadband noise component with one or more quasi-periodic oscillations (QPOs), whose frequencies presumably encode information about the inner accretion flow and black hole properties. For AGN, where the long dynamical times associated with supermassive black holes hamper similar studies, the main focus is on inferring black hole properties from the profiles of relativistically broadened iron emission lines arising from near the marginally stable orbit (Fabian et al. 1989; Tanaka et al. 1995; Nandra et al. 1997; Wilms et al. 2001; Lee et al. 2002). Iron line profiles themselves are highly variable (Iwasawa et al. 1996).

For accreting black holes, and for other relativistically compact objects, the observed lightcurve arises from the modulation of rest-frame variability (i.e. variability that would be detected by an observer orbiting with the disc gas) by photon propagation effects. Since there is now a consensus that magnetorotational instabilities (MRI, Balbus & Hawley 1991; Balbus & Hawley 1998) provide a source of angular momentum transport in well ionized discs, an unavoidable source of rest-frame variability is the fluctuating dissipation inherent to magnetohydrodynamic (MHD) disc turbulence. Simulations of the MRI show that the high frequency temporal power spectra of the mass accretion rate and magnetic stress display steep power-law spectra (Kawaguchi et al. 2000; Hawley & Krolik 2001, 2002), and suggest that the broadband noise component of observed power density spectra may be identified with the temporal fluctuations of MHD disc turbulence. Spatial fluctuations in the dissipation may also play an important role. For a disc inclined to the line of sight, beaming enhances the emission from the approaching side of the disc relative to that from the receeding side. The combination of this photon propagation effect with any source of persistent non-axisymmetric dissipation in the inner regions of the disc could lead to high frequency QPOs in the lightcurve (e.g. Abramowicz et al. 1991; Karas 1999).

In this paper, we use non-relativistic global MHD simulations of the MRI (Armitage 1998; Hawley 2000; Arlt & Rudiger 2001) to model the rest-frame variability of accretion onto a Schwarzschild black hole. A pseudo-Newtonian potential (Paczynski & Wiita 1980) is used to describe the gravitational potential in the vicinity of such a black hole. We extend recent three-dimensional numerical studies of black hole accretion (Hawley 2000; Hawley & Krolik 2001; Armitage, Reynolds & Chiang 2001; Hawley 2001; Reynolds & Armitage 2001; Hawley, Balbus & Stone 2001; Machida, Matsumoto & Mineshige 2001; ; Hawley & Krolik 2002; Igumenshchev & Narayan 2002; Machida & Matsumoto 2003; Igumenshchev, Narayan & Abramowicz 2003) by taking full account of the distinction between the rest-frame disc emission and that seen by a distant observer. We use these simulations to study both the predicted power spectra from black holes accreting via a geometrically thin disc, and the expected variability of the iron line emission that has been observed in the X-ray spectra of Seyfert galaxies (Tanaka et al. 1995; Wilms et al 2001; Lee et al. 2002) and Galactic Black Hole Candidates (Balucinska-Church & Church 2000; Miller et al. 2002). The details of the MHD simulations are presented in Section 2, while Section 3 describes our computation of the observed continuum and line emission taking into account relativistic photon propagation. Predicted temporal power spectra are presented in Section 4. In Section 5, we explore our predictions for variability of the broad iron line. It is suggested that the non-axisymmetric structure discussed in this paper may already have been observed in Chandra High Energy Transmission Grating (HETG) data for the Seyfert galaxy NGC 3516, and we proceed to quantify the nature of the variability that might be observed in the future, using higher signal-to-noise Constellation-X data. Our results, predictions and speculations are summarized in Section 6.

## 2 The MHD simulations

We model the accretion disc around a non-rotating black hole using ideal, Newtonian MHD. The fluid is assumed to be isothermal, with sound speed $`c_s`$, so the equations to be solved are,

$`{\displaystyle \frac{\rho }{t}}+(\rho 𝐯)`$ $`=`$ $`0`$ (1)
$`{\displaystyle \frac{𝐒}{t}}+(\mathrm{𝐒𝐯})`$ $`=`$ $`P\rho \mathrm{\Phi }+𝐉\times 𝐁`$ (2)
$`{\displaystyle \frac{𝐁}{t}}`$ $`=`$ $`\times (𝐯\times 𝐁)`$ (3)
$`P`$ $`=`$ $`\rho c_s^2,`$ (4)
where $`𝐒=\rho 𝐯`$ is the momentum field and the other symbols have their usual meanings. These equations are solved using the ZEUS MHD code (Stone & Norman 1992a, 1992b), operating in cylindrical polar co-ordinates $`(z,r,\varphi )`$. ZEUS is a fixed-grid, time-explicit Eulerian code which uses an artificial viscosity to handle shocks, and a combination of constrained transport (Evans & Hawley 1988) and the method of characteristics (Stone & Norman 1992b) to evolve the magnetic fields. When operated using van Leer advection, as in this work, it is formally of second order accuracy.

A pseudo-Newtonian gravitational potential (Paczynski & Wiita 1980) is used to mimic the effects of general relativity near a non-rotating black hole of mass $`M`$,

$$\mathrm{\Phi }=\frac{GM}{rr_g},$$
(5)
where $`r_g=2GM/c^2`$. This potential reproduces the presence of a marginally stable circular orbit at $`r_{\mathrm{ms}}=6GM/c^2`$. For computational reasons, we also omit the vertical component of the gravitational force. This ‘cylindrical disc’ approximation has been used in several previous studies (Hawley, Gammie & Balbus 1995; Armitage 1998; Reynolds & Armitage 2001; Hawley 2001; Steinacker & Papaloizou 2002), and greatly reduces the computational cost of global disc simulations. Of course, it also precludes study of physical effects such as the Parker instability and the development of a disc corona (e.g. Miller & Stone 2000).

To investigate variability caused by the combination of a turbulent disc with relativistic beaming effects, it is necessary to model all $`2\pi `$ of the disc in azimuth. Apart from this (trivial) extension, our numerical setup is the same as that which we have used previously (Armitage, Reynolds & Chiang 2001; Reynolds & Armitage 2001). A hydrodynamically stable disc flow, threaded by a weak vertical magnetic field, is set up as the initial state of the simulations. We begin by defining a Gaussian surface density profile, centred at 4 $`r_{\mathrm{ms}}`$, where $`r_{\mathrm{ms}}`$ is the radius of the marginally stable circular orbit, with a width of 2.66 $`r_{\mathrm{ms}}`$. This disc, initialized with Keplerian angular velocity and zero magnetic fields, is run with the code in 1D mode until a numerical equilibrium is attained. We then add a weak vertical magnetic field (ratio of gas to magnetic pressure $`\beta =10^4`$) at all radii where the density exceeds a threshold value (set to be approximately 10 percent of the maximum density), and continue the run in 3D. The boundary conditions are periodic in both $`z`$ and $`\varphi `$, and set to outflow at the radial boundaries. ‘Outflow’ boundary conditions are implemented in ZEUS by setting the values of all physical quantities in boundary zones to be equal to their values in the active zones at the edge of the grid.

Two simulations were run, differing only in the adopted sound speed. One run had a ratio of the sound speed to the Keplerian velocity at $`r_{\mathrm{ms}}`$ of $`c_s/v_\varphi =0.061`$, while a second run had $`c_s/v_\varphi =0.041`$ at $`r_{\mathrm{ms}}`$. We will refer to these as the ‘hot’ and ‘cold’ run respectively, though the actual temperature differences are modest. Both simulations were evolved for a total of 140 orbital periods at $`r_{\mathrm{ms}}`$, in identical computational domains,

$`0.7r_\mathrm{g}<`$ $`z`$ $`<0.7r_\mathrm{g}`$ (6)
$`0.7r_{\mathrm{ms}}<`$ $`r`$ $`<8r_{\mathrm{ms}}`$ (7)
$`0<`$ $`\varphi `$ $`<2\pi ,`$ (8)
with 48 grid points in $`z`$, 192 grid points in $`r`$, and 288 grid points in $`\varphi `$. A non-uniform radial grid, with mesh points spaced such that $`r_{i+1}=(1+\delta )r_i`$, with $`\delta `$ a constant, is employed, so that $`\mathrm{\Delta }r/\mathrm{\Delta }\varphi `$ is constant at all radii. The inner radial boundary was placed as close to the marginally stable orbit as possible (to maximise the timestep which is almost always determined by the azimuthal velocity at $`r_{\mathrm{in}}`$), but far enough within $`r_{\mathrm{ms}}`$ that the flow across the boundary was both supersonic and super-Alfvenic. As has been noted by Igumenshchev et al. (2003), this inner boundary condition becomes inappropriate if large net amounts of magnetic flux are able to be dragged inward by the accretion flow. Here, we make the assumption – justified in part by the results of Lubow, Papaloizou & Pringle (1994) – that extensive field dragging does not occur for thin discs.

### 2.1 Properties of the simulations

The basic properties of the MHD disc turbulence seen in the simulations are, unsurprisingly, very similar to those reported previously (i.e. although we need a full $`2\pi `$ in azimuth to study variability, smaller azimuthal domains suffice for most other purposes). Here, we summarize a few features that are pertinent to the subsequent analysis.

Figure 1 shows the total mass in the computational domain as a function of time for the hot and cool runs respectively. Significant accretion commences after around 30 orbits of evolution at the marginally stable orbit. For the hot run, the accretion rate is very roughly constant after around 50 orbits of evolution, and it is this portion of the simulation which we analyse for variability. Note, however, that there remains a secular trend in the accretion rate even over this latter portion (initially increasing, then decreasing once a significant fraction of the initial mass has been accreted). Since this is an artefact of the computational setup, we can only look at variability on substantially shorter timescales than the run length111Because we have only one realization of each simulation, there are also large uncertainties in our estimation of the temporal power spectrum at the lowest frequencies. These purely statistical uncertainties would preclude us from making statements about variability on timescales comparable to the run length, even if the mean accretion rate showed no secular trend.

A single persistent hotspot, orbiting near the innermost stable orbit where relativistic effects are most pronounced, would lead to a strong periodic feature in the lightcurve (Karas 1999). Such QPOs will be diluted in more realistic models, first by the existence of patchy dissipation at all radii, and second by the limited lifetime of regions of enhanced dissipation. Figure 2 illustrates the latter point. We plot the magnetic stress, vertically integrated through the disc, at a single point corotating with the local Keplerian orbital velocity. Results for the hot simulation are shown, although very similar results are obtained for the cold case. Large fluctuations are observed, with the peak stress exceeding the local mean by a factor that can be as large as $`4`$. However, these are also transient features, with a temporal coherence time of only a few orbital periods. If, as we assume later, the magnetic stress provides a proxy for the disc dissipation, it is clear that ‘hotspots’ arising from MHD disc turbulence will lead to at most broad QPOs in power spectra of the disc emission.

Recent analytic (Krolik 1999; Gammie 1999; Agol & Krolik 2000; Afshordi & Paczynski 2003) and numerical work (Hawley 2000; Armitage, Reynolds & Chiang 2001; Reynolds & Armitage 2001; Hawley & Krolik 2001, 2002) has investigated the possibility that the structure of the inner accretion disc might be modified by the existence of magnetic coupling to gas in the plunging region at $`r<r_{\mathrm{ms}}`$. The presence of stress at the last stable orbit is potentially important for studies of variability, as a disc with non-zero stress at $`r_{\mathrm{ms}}`$ has a dissipation profile that is more centrally concentrated than a standard disc (Agol & Krolik 2000). Figure 3 shows how important this effect is for our simulations. We plot the mean stress at $`r_{\mathrm{ms}}`$ as a fraction of the stress at $`2r_{\mathrm{ms}}`$ (this ratio is adopted because the changing accretion rate means that the stress at $`r_{\mathrm{ms}}`$ does not tend to a constant value). At late times, we find that this ratio fluctuates around a value of around 0.3 for the hot disc run. There are, however, large fluctuations, which persist over timescales that can be as long as 10 orbits at $`r_{\mathrm{ms}}`$. As discussed by Papaloizou & Nelson (2003), and by Winters, Balbus & Hawley (2003), these persistent fluctuations are a generic feature of disc MHD turbulence. They are present in local (e.g. Brandenburg et al. 1995; Winters, Balbus & Hawley 2003) as well as global simulations (Papaloizou & Nelson 2003), and their existence means that averaging over long timescales is needed to recover the mean properties of the flow.

As in previous simulations (Reynolds & Armitage 2001), there is evidence that the contribution of stress at $`r_{\mathrm{ms}}`$ becomes relatively less important as the sound speed (or, equivalently, the implied disc thickness $`h/r`$) is reduced. In our subsequent analysis, we compute both raw lightcurves (which include disc dissipation arising from stress at $`r_{\mathrm{ms}}`$), and normalized lightcurves which explicitly fix the form of the radial dissipation profile to match the predictions of zero-torque boundary condition disc models. Comparing these lightcurves allows us to investigate whether the presence of stress at the marginally stable orbit has a significant (and potentially observable) effect on the predicted temporal power spectrum.

## 3 Relativistic transfer function

With a number of approximations, described in Section 4.1, the MHD simulations can be used to define a plausible time-dependent pattern of emission across the disc plane, $`F(r,\varphi ,t)`$. When the region of interest is close to the black hole, as it is in this study, one must also account for general relativistic effects (i.e., beaming, light bending and time-delays) in order to determine the appearance of the simulated disc to an observer at infinity. In this section, we briefly discuss how we incorporate these effects.

Suppose the observer is viewing the disc from a large distance at an inclination of $`i`$ from the normal to the disc plane. We define cartesian coordinates $`(x,y)`$ on the observer’s image plane such that the origin $`x=y=0`$ lies at the center of the event horizon’s image and the $`y`$-axis is parallel with the apparent (i.e. projected) normal to the disc. Let us consider a photon that is emitted from the disc at position $`(r,\varphi )`$ with frequency $`\nu _{\mathrm{em}}`$ and is subsequently observed, a Schwarzschild time $`t`$ later, at position $`(x,y)`$ on the image plane with frequency $`\nu _{\mathrm{obs}}`$. Our task is to compute the functions $`x(r,\varphi ;i)`$, $`y(r,\varphi ;i)`$, $`t(x,y;i)`$, and $`g(x,y;i)`$, where we have defined the frequency shift factor as

$$g(x,y)=\frac{\nu _{\mathrm{obs}}}{\nu _{\mathrm{em}}}.$$
(9)
For a number of inclinations ($`i=1^{},10^{},20^{},\mathrm{},80^{}`$), we integrated null geodesics (photon paths) through the Schwarzschild metric

$`ds^2`$ $`=`$ $`\left(1{\displaystyle \frac{2GM}{c^2r}}\right)dt^2+\left(1{\displaystyle \frac{2GM}{c^2r}}\right)^1dr^2`$ (10)
$`+`$ $`r^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2),`$
which exactly describes the space-time geometry around a non-rotating and uncharged black hole. Starting from a particular point $`(x,y)`$ on the image plane, we integrate photon paths down to the disc plane $`\theta =\pi /2`$ using the first integrals of the geodesic equation (conservation of angular momentum and energy, arising from the $`\varphi `$ and $`t`$ independence of the metric, respectively), the Carter constant, and the light-like nature of the path. For these details, we refer the reader to Appendix A of Reynolds et al. (1999). Turning points in the photon orbits are identified and treated using the method of Rauch & Blandford (1994). These integrations allow us to determine the maps $`x(r,\varphi ;i)`$, $`y(r,\varphi ;i)`$ and $`t(x,y;i)`$.

The frequency shift factor for a given photon path, $`g`$, is given by

$$g=\frac{(p_\mu u^\mu )_{\mathrm{disc}}}{(p_\mu u^\mu )_{\mathrm{obs}}},$$
(11)
where $`p_\mu `$ is the 4-momentum of the photon at the disc/observer, and $`u^\mu `$ is the 4-velocity of the disc/observer. The photon 4-momentum is determined by the photon integration, and we suppose that the observer is stationary (in Schwarzschild coordinates) and very distant from the black hole; thus $`u_{\mathrm{obs}}^\mu =(1,0,0,0)`$. Hence, we can compute $`g`$ for any given photon path provided we know the 4-velocity field on the disc plane. Initially, one might be tempted to take the 3-velocity field directly from our MHD simulations (which are approximately ‘Keplerian’ with respect to the modified potential) and construct the appropriate 4-velocity field as input for our computation of $`g`$. However, this approach is problematic for the following reason. Since our MHD simulation is intrinsically non-relativistic, and based upon a potential that diverges at the event horizon, the 3-velocities within the simulated plunging region can significantly exceed those that would be found in a relativistic calculation. This would lead to artificially high values of $`|g1|`$, and hence artificially strong relativistic beaming.

Motivated by the fact that the velocity field in our simulated disc is very close to Keplerian with respect to the pseudo-Newtonian potential, we choose instead to use the correct relativistic 4-velocity for test-particle circular motion in our computation of $`g`$. This is well defined down to the photon circular orbit at $`3GM/c^2`$, which lies inside the inner edge of our computational grid. Of course, disc material within some radius close to the radius of marginal stability (at $`6GM/c^2`$) will no longer be on circular orbits, but the deviations in the relevant parts of the disc are small. Employing this assumption, we compute the function $`g(x,y;i)`$.

## 4 Temporal variability of the simulated flux

### 4.1 The simulated lightcurve

A direct calculation of the emission from the disc would require, first, that we solve an energy equation, and second that we attain sufficient resolution in a global simulation to model a stratified disc and corona (Miller & Stone 2000) where a significant fraction of the dissipation probably occurs. Since our actual calculation is more limited than this, we must instead choose some proxy for the dissipation. Out of several possible proxies, we have opted to follow Hawley & Krolik (2001), and use the magnetic stress as a local tracer of where in the flow dissipation is expected to occur 222Data storage limitations prompted us to make this basic choice prior to running full resolution simulations. During the runs, vertically integrated slices of $`B_rB_\varphi `$ were the only quantities output at high time resolution.. We then compute lightcurves and power spectra by making the further assumption that the dissipated energy is promptly radiated from the disc.

The most important assumption is obviously the aforementioned one – that at a given radius, the azimuthal dependence of the emission is proportional to the vertically integrated magnetic stress. Even given this assumption, however, there is some ambiguity as to the best way to calculate the lightcurve. As described below, we calculate both ‘raw’ lightcurves, and lightcurves in which the average disc dissipation has been normalized to match an analytic form.

#### 4.1.1 The raw lightcurve

To compute the raw lightcurve, we first calculate the emission map $`F_{\mathrm{raw}}`$ that would be seen by observers orbiting with the disc flow. $`F_{\mathrm{raw}}`$ is given by (e.g. Hubeny & Hubeny 1998).

$$F_{\mathrm{raw}}(r,\varphi ,t)=\sqrt{\frac{GM}{r^3}}\left(\frac{A}{B}\right)B_rB_\varphi dz,$$
(12)
where the relativistic correction factors are,

$`A`$ $`=`$ $`1{\displaystyle \frac{2GM}{rc^2}}`$
$`B`$ $`=`$ $`1{\displaystyle \frac{3GM}{rc^2}}.`$ (13)
Using the transfer function, which gives the mapping between $`(r,\varphi )`$ in the disc plane and $`(x,y)`$ in the image plane, we then compute the image that would be seen by a distant observer viewing the disc at a particular inclination angle, $`i`$. We assume that the observer images the disc in a fixed spectral band, which introduces a ‘K-correction’ that depends upon the spectrum of the disc emission. If at each radius the energy spectrum is $`\nu ^1`$ (fairly typically for accreting black holes), the resulting image is,

$$I(x,y,t)g^4(x,y)F_{\mathrm{raw}}(r,\varphi ,t_{\mathrm{emit}}),$$
(14)
where one factor of $`g`$ comes from the K-correction. We allow for the varying flight time of photons from different parts of the accretion disc by distinguishing between the time $`t_{\mathrm{emit}}`$, at which photons are emitted from the disc, and time $`t`$ when the photons are observed. Examples of these images are shown in Figure 4. The raw lightcurve at time $`t_i`$, $`(t_i)`$ is then simply the sum over the image of the intensity at time $`t_i`$ at each pixel $`(x_j,y_k)`$,

$$(t_i)=\underset{j,k}{}I(x_j,y_k,t_i).$$
(15)
We emphasize that the predicted variability properties depend upon the assumed emission spectrum from the disc as well as on the statistical fluctuations of MHD disc turbulence. In particular, modeling blackbody emission from the disc would require a different radial weighting of the simulated disc dissipation.

The resulting light curves are shown in Figure 5 for several different inclinations. As already noted, the shape of these lightcurves over the longest timescales reflects only the viscous evolution of the ring of gas used for our initial conditions. There is a secular increase in the inner accretion rate (and accompanying dissipation) as the MRI saturates at increasingly large radii, followed (for the hot simulation) by an eventual decline once a substantial fraction of the gas has been accreted. These slow changes in accretion rate are artefacts of the initial conditions used in the simulation — real accretion discs are long lived objects with very large reservoirs of material. Of more relevance to real systems is the inclination dependence of these light curves. One observes a marked increase in rapid variability for the more edge on systems. This is readily interpreted as the effects of relativistic beaming coupled with the ‘blobby’ nature of the simulated emission. We shall explore this point in more detail below when we discuss temporal power spectra.

#### 4.1.2 Novikov-Thorne normalized lightcurve

The raw lightcurve computed from the simulation self-consistently includes the contribution from torque at the marginally stable orbit to the emission of the inner disc. As shown in Figure 3, this contribution can be important, especially for hotter (and therefore geometrically thicker) discs (Reynolds & Armitage 2001). However, because the domain of the simulation is limited, only the inner region of the simulated disc is in a steady-state. The raw lightcurve therefore underestimates the contribution to the flux from larger radii in the disc.

Cognisant of this drawback, we also compute lightcurves by combining the simulated fluctuations in the stress (in $`\varphi `$ and $`t`$) with an analytic dissipation profile due to Novikov & Thorne (1973). Specializing to the case of a Schwarzschild black hole, the radial dependence of the disc flux measured by an observer orbiting with the disc gas is (Page & Thorne 1974),

$`F_{\mathrm{NT}}(r){\displaystyle \frac{1}{x^4(x^33x)}}[x\sqrt{6}{\displaystyle \frac{\sqrt{3}}{2}}\mathrm{ln}\left({\displaystyle \frac{x\sqrt{3}}{\sqrt{6}\sqrt{3}}}\right)`$
$`+{\displaystyle \frac{\sqrt{3}}{2}}\mathrm{ln}({\displaystyle \frac{x+\sqrt{3}}{\sqrt{6}+\sqrt{3}}},)]`$ (16)
where,

$$x\sqrt{\frac{r}{(GM/c^2)}}.$$
(17)
Given this dissipation profile, we proceed to construct a lightcurve as before, except that we replace the emission map $`F_{\mathrm{raw}}(r,\varphi ,t)`$ with $`F_{\mathrm{NT}}(r,\varphi ,t)`$, defined as,

$$F_{\mathrm{NT}}(r,\varphi ,t)=\frac{B_rB_\varphi (r,\varphi ,t)}{<B_rB_\varphi >(r)}F_{\mathrm{NT}}(r).$$
(18)
In this expression $`B_rB_\varphi `$ in both the numerator and the denominator is understood to be vertically integrated, and the angle brackets denote averaging over both $`\varphi `$ and $`t`$. By construction, the resulting dissipation map $`F_{\mathrm{NT}}(r,\varphi ,t)`$ has the same time-averaged dissipation profile as a Novikov-Thorne disc, on top of which are superimposed fluctuations whose statistical properties are computed directly from the MHD simulation. Using this normalized dissipation map (which we shall refer to as a ‘NT-normalization’), one can construct the observed emission map and a lightcurve in exactly the manner discussed above.

### 4.2 The temporal power spectrum

For a given lightcurve, we compute the temporal power spectrum $`G(f)`$, where $`f`$ is the frequency, via a straightforward Fourier transform together with the “Leahy normalization”,

$$G(f)=\frac{2|_i(t_i)\mathrm{exp}(2\pi jft_i)|^2}{_i(t_i)}.$$
(19)
We truncate the resulting power spectra at the high-frequency end at the Nyquist frequency (i.e., half of the sampling frequency of the lightcurve) to account for the finite sampling time. In addition, frequencies are rescaled such that the orbital frequency at the radius of marginal stability is unity.

Figure 6 shows power spectra computed using emission that originates from several narrow annuli in the inner disc of the hot simulation. The power spectra are described well by broken power laws. Below a break frequency, which is comparable to the local orbital frequency, $`G(f)`$ is flat or modestly declining ($`G(f)f^1`$). At higher frequencies, above the break, we obtain a steep spectrum with $`G(f)f^{3.5}`$.

Power spectra computed using the integrated disc emission (ie directly from the lightcurves in Figure 5) are shown in Figures 7 and 8 for the cold and hot runs respectively. For the cold run, $`G(f)`$ calculated from the raw lightcurve is accurately approximated by a power law with an index of -2 for all frequencies which are accessible using the simulation. There are no significant changes in the slope of $`G(f)`$ with increasing system inclination. Figure 7 also shows power spectra computed from the NT-normalized lightcurves. Although not identical, these are almost indistinguishable from power spectra calculated using the raw lightcurve. We conclude that the presence of torque at the marginally stable orbit – at least at the rather modest levels obtained in our calculations – does not significantly affect the predicted power spectrum.

The power spectra shown in Figure 7 appear to be accurately described by single power laws over the range of frequencies accessible using the simulation. Observationally, high frequency QPOs in black hole systems are typically weak features. For example, the broad QPO at 300 Hz reported by Remillard et al. (1999) in GRO J1655-40 had an amplitude of less that 1%. Substantially longer simulations would be needed before we could measure the theoretical power spectrum accurately enough to detect broad QPOs of such small amplitude. We can, however, place limits on the possible existence of stronger QPOs that could arise from the ‘lighthouse’ mechanism discussed by Abramowicz et al. (1991) and Karas (1999). Such QPOs, since they originate from the combination of hotspots with relativistic beaming, should become stronger at higher system inclinations. To investigate this possibility, the lower panels of Figure 7 show the relative power for inclined systems as compared to face-on systems. For $`i=50^{}`$, any QPOs at frequencies comparable to the orbital frequency contribute at most $`20`$% of the broadband power.

Figure 8 shows the corresponding power spectra extracted from the hot simulation. The results are generally consistent, though the deviations from a power law form for $`G(f)`$ are larger. At higher inclinations there are some indications of the presence of excess power at a frequency of around 0.2–0.3 times the orbital frequency at $`r_{\mathrm{ms}}`$. Excess power is also seen in the cold simulation at similar frequencies, though at a lower level. Longer simulations are required to investigate the reality of such weak, broad features.

Since the emission from a particular annulus has a power-spectrum that breaks at the orbital frequency of that annulus, we expect that the disk-integrated power-spectrum should possess a high-frequency break corresponding to the innermost emitting annulus. The overall power-spectra should therefore steepen appreciably above $`f1`$. In fact, we do not see such a high-frequency break in Figures 7 and 8. We have investigated the source of this additional high frequency noise in our lightcurves, and find that it originates from emission from the outermost parts of the simulated disc. Our suspicion is that the spatial discretization of the disc into cells – which are physically largest in the outer disc – introduces a low level of high frequency flicker into the lightcurve. This hampers our ability to see high-frequency cutoffs in the disk-integrated power-spectra.

## 5 Iron line variability

The broad iron line is one of the best understood observational probes of black hole accretion discs at the present time (e.g., see reviews by Fabian et al. 2000; Reynolds & Nowak 2003). Thus, it is useful to examine the predictions of our simulation for temporal variability of the iron line.

The predicted temporal variability of the iron line obviously depends on how that line is excited. Different authors have computed the iron line response to flares in the inner regions of the flow, using a variety of different geometries and kinematics for the source of exciting photons (Fabian et al. 1989; Stella 1990, Matt & Perola 1992; Reynolds et al. 1999; Ruszkowski 2000). In this work, we assume that the local iron line intensity coming from the surface of the disc is proportional to the local continuum intensity (and, hence, is proportional to the local vertically integrated azimuthal magnetic stress). The line intensity at image-plane position $`(x,y)`$, time $`t`$ and energy $`E`$ is then given by

$$I_{\mathrm{line}}(x,y,t;E)g^4(x,y)F(r,\varphi ,t_{\mathrm{emit}})\delta (EgE_0),$$
(20)
where $`E_0`$ is the rest-frame energy of the emission line and $`F(r,\varphi ,t)`$ is taken to be either the raw of NT-normalized disc flux. Note that the K-correction is not relevant for computing line fluxes and, in this expression, one power of $`g`$ results from the transformation of the frequency-element. The total simulated observed line is obtained from integrating this expression over the image plane. Figure 9 shows instantaneous iron line profile resulting from the hot simulation for an inclination of $`i=40^{}`$. The skewed and redshifted profile is caused by a combination of the line-of-sight Doppler effect, the transverse Dopper effect (i.e. special relativistic time dilation) and gravitational redshifts (see Fabian et al. 2000 and Reynolds & Nowak 2003 for a detailed discussion of the physics determining these line profiles).

Under the assumptions employed above, the simulated iron line at any particular instant is a weighted redshift map of the turbulent emissivity pattern across the disc. Thus, the temporal and azimuthal variability caused by the MHD turbulence will produce a level of inevitable and irreducible temporal variability in the iron line profile. This can be seen by comparing the instantaneous line profiles in Figure 9. From this figure, it can be seen that the iron line from the hot disc displays dramatic (factor of $`2`$) variability between independent time slices, whereas the line profile from the cold disc is more stable, showing deviations at only the $`30`$ percent level.

The nature of the iron line variability is demonstrated more clearly in Figure 10, which shows the intensity of the simulated iron line emission as a function of energy and time for a representative interval during our simulations, assuming an inclination of $`i=20^{}`$. This is a plausible inclination for many Seyfert-1 nuclei, the systems in which broad iron lines have been best studied to date. One can see a tremendous degree of structure in the iron line profile as a function of both energy and time. Of particular note are the distinctive arcs that correspond to the orbital motion of especially bright regions within the turbulent flow, visible in Figure 4. These arcs essentially trace out (in the energy-time plane) the paths of test particle with orbits close to the radius of marginally stability. If they could be detected, they would provide a quantitative test of the flow extremely close to the black hole. High temporal resolution is required to observe such features, but the required capabilities should be readily within reach of the planned Constellation-X and XEUS missions.

In fact, observational evidence for transient non-axisymmetric structure already exists. Chandra-HETG observations of the Seyfert-1 galaxy NGC 3516 have found substructure within the well known broad iron line profile of this object (Turner et al. 2002). The substructure has the form of narrow line-line features at rest-frame energies of 5.6 keV, 6.2 keV and 7.1 keV. Plausible explanations exist for the latter two features; the 7.1 keV feature can be readily attributed to the fluorescent K$`\beta `$ emission of cold iron, and the 6.2 keV feature might be the ‘Compton shoulder’ of the prominent 6.4 keV K$`\alpha `$ fluorescent line of cold iron. However, the 5.6 keV feature cannot be attributed to any expected features. Turner et al. suggest that this feature is, in fact, iron K$`\alpha `$ fluorescence with a rest-frame energy of 6.4 keV from an anomalously bright annulus or spot on the accretion disc surface. The flux contained in the 5.6 keV feature is approximately 5% of that contained in the overall broad iron line.

Is it possible to relate the substructure observed in the NGC 3516 iron line with the turbulence induced substructure seen in our simulations? While large amplitude substructure is present in virtually every instantaneous simulated line profile, it is much diminished once the line profiles are averaged over more than a few orbits. The Chandra-HETG observation of Turner et al. that yielded the clear detection of the iron line substructure was $`8\times 10^4`$ s in duration. Assuming a black hole mass of $`2.3\times 10^7\mathrm{M}_{}`$ (estimated using the velocity dispersion from di Nella et al. 1995, in conjunction with the scaling relations of Gebhardt et al. 2000 or Ferrarese & Merritt 2000)333Independent estimates of the black hole mass are comparable. For example, Onken et al. (2003) estimate $`M_{BH}=1.7\times 10^7M_{}`$ from reverberation mapping., this corresponds to approximately 20 orbital periods of the accretion disc at the radius of marginal stability. Thus, to compare our simulations to these data, we averaged our simulated line profiles over 20 orbits of the inner disc. We assumed a disc inclination of $`i=40^{}`$ (Wu & Han 2001). Given the useful length of our simulation, five independent 20-orbit line profiles could be formed. We found that 4 out of 5 of these profiles possessed no strong deviations from the average. In one profile, however, there was a slight excess centered at an energy of 5.5 keV, possessing 2–3% of the total flux. Therefore, while the true interpretation of the Chandra data for NGC 3516 remains unclear, turbulence-induced substructure is a viable possibility even when one considers that the disc has likely undergone several orbits during the observation.

## 6 Summary

In this paper, we have explored the predicted variability of accretion discs around Schwarzschild black holes. The temporal and spatial fluctuations in the magnetic stress which drives accretion have been calculated using three-dimensional global simulations in a pseudo-Newtonian potential, which ought to provide a reasonable first approximation to the disc dynamics outside the marginally stable orbit. By making the further critical assumption that the radiation from the disc follows the local stress, we have calculated the predicted emission that would be seen by an observer using a fully relativistic ray-tracing code that accounts for the relativistic effects of beaming, gravitational redshift, and light bending.

Our results for the temporal power spectrum of the predicted emission are consistent with those of previous authors (Kawaguchi et al. 2000; Hawley & Krolik 2001, 2002). Fluctuations in the magnetic stress, when integrated across the disc, produce a high frequency power spectrum described by a power law with slope $`f^2`$. To a first approximation, observers viewing the system at a wide range of inclinations – from face-on to almost edge-on – are predicted to derive similar power spectra. In more detail, high inclinations are found to boost the relative power at frequencies a factor of a few below the orbital frequency at the marginally stable orbit. There is little change in the predicted power spectra due to the presence or absence of modest levels of magnetic stress at the marginally stable orbit. No strong QPOs are seen in the simulations, though weak features could be present and would not be clearly observed given the limited statistics available from $`10^2`$ orbit simulations.

The broad iron line provides an alternate probe of conditions in the inner disc, and has the advantage of retaining velocity information. We have used the simulations to quantify the predicted iron line variability within a model in which the local iron line emissivity tracks the locally generated continuum emission. This model represents the opposite extreme from reverberation models, in which the source of line exciting flux is often assumed to be widely spatially separated from the disc. We find that order unity fluctuations in the line flux at fixed energy occur for discs with $`h/r0.1`$, with the fluctuation amplitude increasing rapidly with increasing disc thickness. Averaging over at least several orbits at the radii of interest is needed before meaningful averages can be extracted (see also Papaloizou & Nelson 2003; Winters, Balbus & Hawley 2003). Observationally, this implies that long integrations, extending over tens of orbital timescales, are required in order to reduce the fluctuation amplitude to percent levels and thereby recover the mean line profile.

## Acknowledgements

We thank Cole Miller, Eve Ostriker, Daniel Proga and Mateusz Ruszkowski for extensive discussions, and Pete Ruprecht for assisting with computational support. We thank the referee for providing an extremely prompt report on the paper. CSR acknowledges support from the National Science Foundation under grant AST0205990.