text:
# Enhancement of 𝓑(𝑩̄_𝒅→𝝁⁺𝝁⁻)/𝓑(𝑩̄_𝒔→𝝁⁺𝝁⁻) in the MSSM with Modified Minimal Flavour Violation and Large 𝐭𝐚𝐧𝜷

## I Introduction

In the Standard Model (SM), the dominant contributions to the decays $`\overline{B}_ql^+l^{}`$, where $`q=d,s`$ and $`l=e,\mu ,\tau `$, arise from $`Z^0`$-penguin diagrams and box diagrams involving $`W^\pm `$ bosons. The contributions due to neutral Higgs boson exchange, on the other hand, are suppressed by a factor of $`m_lm_{b,q}/M_W^210^3`$, and therefore are completely negligible.

The SM predicts the $`\overline{B}_s\mu ^+\mu ^{}`$ branching ratio to be

$$(\overline{B}_s\mu ^+\mu ^{})=(3.2\pm 1.5)\times 10^9,$$
(1)
and the ratio of branching fractions

$$R_{\mathrm{SM}}\frac{(\overline{B}_d\mu ^+\mu ^{})}{(\overline{B}_s\mu ^+\mu ^{})}|_{\mathrm{SM}}\frac{\tau _{B_d}}{\tau _{B_s}}\frac{M_{B_d}}{M_{B_s}}\frac{f_{B_d}^2}{f_{B_s}^2}\frac{|V_{td}|^2}{|V_{ts}|^2}O(10^2),$$
(2)
where $`\tau _{B_q}`$ is the lifetime of the $`B_q`$ meson, $`M_{B_q}`$ and $`f_{B_q}`$ are the corresponding mass and decay constant. We note that the main uncertainty on the branching ratio in Eq. (1) arises from $`f_{B_s}`$, with a typical value of $`(230\pm 30)\text{MeV}`$ obtained from lattice QCD calculations . On the other hand, the relative rates of $`\overline{B}_d`$ and $`\overline{B}_s`$ decays \[Eq. (2)\] have a smaller theoretical uncertainty due to the appearance of the ratio $`f_{B_d}/f_{B_s}`$, which can be more precisely determined than $`f_{B_s}`$ alone. A determination of $`R`$ can therefore potentially provide valuable information on $`|V_{td}/V_{ts}|`$ in the SM. However, given the SM prediction of $`(\overline{B}_d\mu ^+\mu ^{})O(10^{10})`$, the $`\overline{B}_d\mu ^+\mu ^{}`$ decay is experimentally remote unless it is significantly enhanced by new physics. Thus, the purely leptonic decays of neutral $`B`$ mesons provide an ideal testing ground for physics outside the SM, with the current experimental upper bounds

$$(\overline{B}_d\mu ^+\mu ^{})<2.8\times 10^7(90\%\mathrm{C}.\mathrm{L}.),$$
(4)
$$(\overline{B}_s\mu ^+\mu ^{})<2.0\times 10^6(90\%\mathrm{C}.\mathrm{L}.).$$
(5)
Our main interest in this paper is in a qualitative comparison of the $`\overline{B}_{d,s}\mu ^+\mu ^{}`$ branching fractions in the presence of non-standard interactions, which can be made by using the ratio

$$R\frac{(\overline{B}_d\mu ^+\mu ^{})}{(\overline{B}_s\mu ^+\mu ^{})}.$$
(6)
Referring to Eq. (2), it is important to note that the suppression of $`R`$ in the SM is largely due to the ratio of the Cabibbo-Kobayashi-Maskawa (CKM) elements. This dependence on the CKM factors allegedly pertains to all models in which the quark mixing matrix is the only source of flavour violation. It is therefore interesting to ask if $`R`$ could be of the order unity in some non-standard models where flavour violation is governed exclusively by the CKM matrix. Working in the framework of the Minimal Supersymmetric Standard Model (MSSM) with a large ratio of Higgs vacuum expectation values, $`\mathrm{tan}\beta `$ (ranging from $`40`$ to $`60`$), we show that such a scenario does exist, and study its consequences for the $`\overline{B}_{d,s}\mu ^+\mu ^{}`$ branching ratios.***A number of authors  have noted that in the large $`\mathrm{tan}\beta `$ regime neutral Higgs-boson contributions can increase the SM decay rates of $`\overline{B}_{d,s}\mu ^+\mu ^{}`$ by up to several orders of magnitude. Moreover, there is a correlation between $`(\overline{B}_s\mu ^+\mu ^{})`$ and $`B_s^0`$$`\overline{B_s^0}`$ mixing , and, in mSUGRA scenarios, between $`(\overline{B}_s\mu ^+\mu ^{})`$ and $`(g2)_\mu `$  (see also Ref. ). To this end, we extend our previous analysis  to include the effects of gluino and neutralino contributions, in addition to those coming from diagrams with $`W^\pm ,\stackrel{~}{\chi }^\pm `$, charged and neutral Higgs boson exchange.

The outline of this paper is as follows. In Sec. II, we define modified minimal flavour violation ($`\overline{\mathrm{MFV}}`$) and discuss briefly three distinct scenarios within the MSSM. The effective Hamiltonian describing the decays $`\overline{B}_{d,s}\mu ^+\mu ^{}`$ in the presence of non-SM interactions is given in Sec. III, together with the branching ratios. In Sec. IV, we present analytic formulae for the gluino and neutralino contributions to scalar and pseudoscalar Wilson coefficients governing the $`bql^+l^{}`$ transition. Numerical results for the branching fractions $`(\overline{B}_{d,s}\mu ^+\mu ^{})`$ and the ratio $`R`$ are presented in Sec. V. In Sec. VI, we summarize and conclude.

## II The MSSM with Modified Minimal Flavour Violation

There exists no unique definition of minimal flavour violation (MFV) in the literature (see, e.g., Refs. ). The common feature of these MFV definitions is that flavour violation and/or flavour-changing neutral current (FCNC) processes are entirely governed by the CKM matrix. On the other hand, they differ, for example, by the following additional assumptions:

* there are no new operators present, in addition to those of the SM ,
* FCNC processes are proportional to the same combination of CKM elements as in the SM ,
* flavour transitions occur only in charged currents at tree level .

While these ad hoc assumptions are useful for certain considerations, such as the construction of the universal unitarity triangle , they cannot be justified by symmetry arguments on the level of the Lagrangian. For example, the number of operators with a certain dimension is always fixed by the symmetry of the low-energy effective theory. Whether the Wilson coefficients are negligible or not, depends crucially on the model considered and on the part of the parameter space. Furthermore, the requirement that FCNC processes are proportional to the same combination of CKM elements as in the SM fails, for example, in the MSSM and can be retained only after further simplifying assumptions. The last statement in iii) is of pure phenomenological relevance in order to avoid huge contributions to FCNC processes.

Using symmetry arguments, we propose an approach that relies only on the key ingredient of the MFV definitions in Refs. , without considering the above mentioned additional assumptions i)–iii). We call an extension of the SM a modified minimal flavour-violating ($`\overline{\mathrm{MFV}}`$) model if and only if FCNC processes or flavour violation are entirely ruled by the CKM matrix; that is, we require that FCNC processes vanish to all orders in perturbation theory in the limit $`V_{\text{CKM}}𝟙`$.

Let us briefly motivate our definition of $`\overline{\mathrm{MFV}}`$. The structure of the SM is such that the conservation of lepton numbers $`L_e`$, $`L_\mu `$ and $`L_\tau `$ is exact. This is reflected by three global $`U(1)`$ symmetries (one for each lepton flavour) of the SM Lagrangian. As a consequence, there are no FCNC processes in the lepton sector to all orders in perturbation theory. We regard the breaking of these lepton flavour $`U(1)`$ symmetries as flavour violation.This is different from Ref. , where flavour violation is regarded as the breaking of the flavour symmetry $`G_F`$. As far as the quark sector of the SM is concerned, the CKM matrix explicitly breaks the quark flavour $`U(1)`$ symmetries. Thus, in the quark sector, FCNC processes can be generated in electroweak interactions. In the limit $`V_{\text{CKM}}𝟙`$ the three $`U(1)`$ quark symmetries are restored and therefore FCNC transitions vanish to all orders in perturbation theory. This observation leads us to the definition of $`\overline{\mathrm{MFV}}`$ given above. As will become clear, the advantage of $`\overline{\mathrm{MFV}}`$ is that it is less restrictive than MFV, while the CKM matrix remains the only source of FCNC transitions.

Our definition of $`\overline{\mathrm{MFV}}`$ is manifest basis independent. However, in order to find a useful classification of different $`\overline{\mathrm{MFV}}`$ scenarios within the MSSM, we will work in the super-CKM basis (see Ref.  for details). In this basis the quark mass matrices are diagonal, and both quarks and squarks are rotated simultaneously. The scalar quark mass-squared matrices in this basis have the structure

$$_U^2=\left(\begin{array}{cc}_{U_{LL}}^2& _{U_{LR}}^2\\ _{U_{LR}}^2& _{U_{RR}}^2\end{array}\right),_D^2=\left(\begin{array}{cc}_{D_{LL}}^2& _{D_{LR}}^2\\ _{D_{LR}}^2& _{D_{RR}}^2\end{array}\right),$$
(7)
where the $`3\times 3`$ submatrices are given by

$$_{U_{LL}}^2=M_{\stackrel{~}{U}_L}^2+M_U^2+\frac{1}{6}M_Z^2\mathrm{cos}2\beta (34\mathrm{sin}^2\theta _W)𝟙,$$
(9)
$$_{U_{LR}}^2=M_U(A_U^{}\mu \mathrm{cot}\beta 𝟙),$$
(10)
$$_{U_{RR}}^2=M_{\stackrel{~}{U}_R}^2+M_U^2+\frac{2}{3}M_Z^2\mathrm{cos}2\beta \mathrm{sin}^2\theta _W𝟙,$$
(11)
$$_{D_{LL}}^2=M_{\stackrel{~}{D}_L}^2+M_D^2\frac{1}{6}M_Z^2\mathrm{cos}2\beta (32\mathrm{sin}^2\theta _W)𝟙,$$
(13)
$$_{D_{LR}}^2=M_D(A_D^{}\mu \mathrm{tan}\beta 𝟙),$$
(14)
$$_{D_{RR}}^2=M_{\stackrel{~}{D}_R}^2+M_D^2\frac{1}{3}M_Z^2\mathrm{cos}2\beta \mathrm{sin}^2\theta _W𝟙.$$
(15)
Here, $`M_{\stackrel{~}{U}_{L,R}}^2`$ and $`M_{\stackrel{~}{D}_{L,R}}^2`$ are the soft SUSY breaking squark mass-squared matrices, $`\mu `$ is the Higgsino mixing parameter, $`A_U`$ and $`A_D`$ are soft SUSY breaking trilinear couplings, $`𝟙`$ denotes the $`3\times 3`$ unit matrix, and

$$M_U\mathrm{diag}(m_u,m_c,m_t),M_D\mathrm{diag}(m_d,m_s,m_b).$$
(16)
Because of SU(2) gauge invariance, the mass matrix $`M_{\stackrel{~}{D}_L}^2`$ is intimately connected to $`M_{\stackrel{~}{U}_L}^2`$ via

$$M_{\stackrel{~}{D}_L}^2=V_{\text{CKM}}^{}M_{\stackrel{~}{U}_L}^2V_{\text{CKM}},$$
(17)
which is important for our subsequent discussion. In the case of $`\overline{\mathrm{MFV}}`$, the matrices $`A_U,A_D`$, $`M_{\stackrel{~}{U}_R}^2,M_{\stackrel{~}{D}_R}^2`$ must be diagonal:

$$A_U\mathrm{diag}(A_u,A_c,A_t),A_D\mathrm{diag}(A_d,A_s,A_b),$$
(18)
$$M_{\stackrel{~}{U}_R}^2\mathrm{diag}(m_{\stackrel{~}{u}_R}^2,m_{\stackrel{~}{c}_R}^2,m_{\stackrel{~}{t}_R}^2),M_{\stackrel{~}{D}_R}^2\mathrm{diag}(m_{\stackrel{~}{d}_R}^2,m_{\stackrel{~}{s}_R}^2,m_{\stackrel{~}{b}_R}^2).$$
(19)
Furthermore, assuming that there are no new CP-violating phases, in addition to the single CKM phase, the matrices $`A_U`$ and $`A_D`$, as well as $`\mu `$, are real.

Taking into account the relation in Eq. (17), one encounters three cases of $`\overline{\mathrm{MFV}}`$.

* Scenario (A):
$`M_{\stackrel{~}{U}_L}^2`$ is proportional to the unit matrix, and so $`M_{\stackrel{~}{U}_L}^2=M_{\stackrel{~}{D}_L}^2`$. As a result, there are no gluino and neutralino contributions to flavour-changing transitions at one-loop level. The implications of this scenario for rare $`B`$ decays have been discussed, e.g., in Refs. . This scenario of $`\overline{\mathrm{MFV}}`$ coincides with the MFV scenario at low $`\mathrm{tan}\beta `$, as defined in Refs. .
* Scenario (B):
$`M_{\stackrel{~}{D}_L}^2`$ is diagonal but not proportional to the unit matrix and, in consequence, $`M_{\stackrel{~}{U}_L}^2`$ has non-diagonal entries. In such a case, there are again no gluino and neutralino contributions to flavour-changing one-loop transitions involving only external down-type quarks and leptons. However, additional chargino contributions show up, due to non-diagonal entries of $`M_{\stackrel{~}{U}_L}^2`$. This scenario has recently been investigated in Ref. .
* Scenario (C):
$`M_{\stackrel{~}{U}_L}^2`$ is diagonal but not proportional to the unit matrix, which gives rise to off-diagonal entries in $`M_{\stackrel{~}{D}_L}^2`$. Accordingly, gluino and neutralino exchange diagrams (in addition to those involving $`W^\pm ,\stackrel{~}{\chi }^\pm `$, charged and neutral Higgs bosons) contribute to flavour-changing transitions at one-loop level that involve external down-type quarks.

The common feature of all these scenarios is that the CKM matrix is the only source of flavour violation. For the remainder of this paper we will concentrate mainly on scenarios (B) and (C). To our knowledge, the consequences of scenario (C) for the decays $`\overline{B}_{d,s}\mu ^+\mu ^{}`$ have not yet been discussed in the literature.

## III Effective Hamiltonian and branching ratio

The effective Hamiltonian responsible for the processes $`\overline{B}_ql^+l^{}`$, with $`q=d,s`$ and $`l=e,\mu ,\tau `$, in the presence of non-standard interactions is given byWe omit the index $`q`$ on the operators and the corresponding short-distance coefficients.

$$H_{\mathrm{eff}}=\frac{G_F\alpha }{\sqrt{2}\pi }V_{tb}V_{tq}^{}\underset{i=10,S,P}{}[c_i(\mu )𝒪_i(\mu )+c_i^{}(\mu )𝒪_i^{}(\mu )],$$
(20)
with the short-distance coefficients $`c_i^{()}(\mu )`$ and the local operators

$$𝒪_{10}=(\overline{q}\gamma ^\mu P_Lb)(\overline{l}\gamma _\mu \gamma _5l),𝒪_{10}^{}=(\overline{q}\gamma ^\mu P_Rb)(\overline{l}\gamma _\mu \gamma _5l),$$
(22)
$$𝒪_S=m_b(\overline{q}P_Rb)(\overline{l}l),𝒪_S^{}=m_q(\overline{q}P_Lb)(\overline{l}l),$$
(23)
$$𝒪_P=m_b(\overline{q}P_Rb)(\overline{l}\gamma _5l),𝒪_P^{}=m_q(\overline{q}P_Lb)(\overline{l}\gamma _5l),$$
(24)
where $`P_{L,R}=(1\gamma _5)/2`$. In addition to the operators in Eqs. (20), there are tensor and vector operators, $`(\overline{q}\sigma ^{\mu \nu }P_{L,R}b)(\overline{l}\sigma _{\mu \nu }P_{L,R}l)`$ and $`(\overline{q}\gamma ^\mu P_{L,R}b)(\overline{l}\gamma _\mu l)`$. However, they do not contribute to the $`\overline{B}_ql^+l^{}`$ decays and do not mix with those appearing in Eq. (20). As a matter of fact, the matrix element that involves the antisymmetric tensor $`\sigma ^{\mu \nu }`$ must vanish since $`p^\mu p_{B_q}^\mu `$ is the only four-momentum vector available. Similarly, the matrix element $`l^+l^{}|\overline{l}\gamma _\mu l|0`$ does not contribute when contracted with $`0|\overline{q}\gamma ^\mu P_{L,R}b|\overline{B}_q(p)p^\mu `$ by the equation of motion.

Evolution of the short-distance coefficients $`c_i^{()}(\mu )`$ from the matching scale $`\mu _t=m_t^{\mathrm{pole}}`$ down to the low-energy scale $`\mu _b=m_b^{\mathrm{pole}}`$ can be performed by means of the renormalization group equation. We note that the anomalous dimensions of the operators in Eqs. (20) vanish, and thus the renormalization group evolution is trivial.

Because of the pseudoscalar nature of the decaying $`B_q`$ meson, the hadronic matrix elements are non-zero only for axial-vector and pseudoscalar operators; namely,

$$0|\overline{q}\gamma _\mu \gamma _5b|\overline{B}_q(p)=ip_\mu f_{B_q},$$
(25)
and, employing the equation of motion,

$$0|\overline{q}\gamma _5b|\overline{B}_q(p)=if_{B_q}\frac{M_{B_q}^2}{m_b+m_q}.$$
(26)
Here $`f_{B_q}`$ is the $`B_q`$ meson decay constant, which can be obtained from lattice QCD computations :

$$f_{B_d}=(200\pm 30)\text{MeV},f_{B_s}=(230\pm 30)\text{MeV},\frac{f_{B_s}}{f_{B_d}}=1.16\pm 0.04.$$
(27)
(Similar results have recently been obtained from a QCD sum rule analysis .)

If the lepton spins are not measured, the branching ratio for the case $`l=\mu `$ takes the general form

$$(\overline{B}_q\mu ^+\mu ^{})=\frac{G_F^2\alpha ^2M_{B_q}^3f_{B_q}^2\tau _{B_q}}{64\pi ^3}|V_{tb}V_{tq}^{}|^2\sqrt{14\widehat{m}_\mu ^2}\left\{(14\widehat{m}_\mu ^2)|F_S|^2+|F_P+2\widehat{m}_\mu F_A|^2\right\},$$
(28)
with the notation $`\widehat{m}_\mu m_\mu /M_{B_q}`$ and the dimensionless form factors

$$F_{S,P}=M_{B_q}\left[\frac{c_{S,P}m_bc_{S,P}^{}m_q}{m_b+m_q}\right],F_A=c_{10}c_{10}^{}.$$
(29)
In the SM, the contributions involving the neutral Higgs boson are completely negligible, and so $`(\overline{B}_q\mu ^+\mu ^{})\widehat{m}_\mu ^2`$, which is a consequence of helicity suppression.§§§Because the $`B`$ meson has spin zero, both $`\mu ^+`$ and $`\mu ^{}`$ must have the same helicity. Finally, allowing the input parameters in Eq. (28) to vary over the interval in Eq. (27), and using the ranges for the CKM factors given in Ref. , the ratio of decay rates of $`\overline{B}_{d,s}\mu ^+\mu ^{}`$ within the SM is expected to be in the range

$$0.02R_{\mathrm{SM}}0.05,$$
(30)
which is largely due to the imprecisely known CKM elements.

## IV Higgs-boson contributions to the decays $`\overline{𝑩}_{𝒅\mathbf{,}𝒔}\mathbf{}𝒍^\mathbf{+}𝒍^{\mathbf{}}`$

We now turn to the computation of the scalar and pseudoscalar Wilson coefficients in the $`bql^+l^{}`$ transition arising from gluino and neutralino exchange diagrams within the general MSSM. As mentioned earlier, we perform our calculation in the large $`\mathrm{tan}\beta `$ regime (i.e. $`40\mathrm{tan}\beta 60`$). For the remaining contributions ($`W^\pm ,H^\pm ,\stackrel{~}{\chi }^\pm `$) to these short-distance coefficients, we refer to Refs. .

The relevant box and penguin diagrams are depicted in Fig. 1, where $`H^0,h^0,A^0`$ and $`G^0`$ are the neutral Higgs and would-be-Goldstone bosons respectively, $`\stackrel{~}{l}_a`$ are the charged sleptons, $`\stackrel{~}{d}_a`$ denote the down-type squarks, $`\stackrel{~}{\chi }_k^0`$ are the neutralinos, and $`\stackrel{~}{g}`$ represents the gluino. We perform the calculation in the ’t Hooft-Feynman gauge, using the Feynman rules of Ref. , and adopting the on-shell renormalization prescription described in Ref. .

In our subsequent calculation, we will exploit the tree-level relations

$$M_{A^0}^2=M_H^2M_W^2,M_{H^0}^2+M_{h^0}^2=M_{A^0}^2+M_Z^2,$$
(31)
$$\frac{\mathrm{sin}2\alpha }{\mathrm{sin}2\beta }=\left(\frac{M_{H^0}^2+M_{h^0}^2}{M_{H^0}^2M_{h^0}^2}\right),\frac{\mathrm{cos}2\alpha }{\mathrm{cos}2\beta }=\left(\frac{M_{A^0}^2M_Z^2}{M_{H^0}^2M_{h^0}^2}\right),$$
(32)
where $`M_{A^0}`$ and $`M_H`$ are the masses of the $`\mathrm{CP}`$-odd and charged Higgs boson respectively. $`M_{h^0,H^0}`$ and $`\alpha `$ are the masses and mixing angle in the $`\mathrm{CP}`$-even Higgs sector. This leaves two free parameters in the Higgs sector that we choose to be $`M_H`$ and $`\mathrm{tan}\beta `$.

Introducing, for convenience, the dimensionless variables

$$x_a=\frac{m_{\stackrel{~}{d}_a}^2}{M_{\stackrel{~}{g}}^2},x_{ak}=\frac{m_{\stackrel{~}{d}_a}^2}{M_{\stackrel{~}{\chi }_k^0}^2},$$
(33)
our results for the gluino and neutralino contributions at large $`\mathrm{tan}\beta `$ can be summarized as follows.Our results for the gluino and neutralino contributions differ somewhat from those given in Ref. , which contains typographical errors in the formulae for the various Wilson coefficients .

(i) Gluino:

$`c_{S,P}^{\stackrel{~}{g}}`$ $`=`$ $`\pm {\displaystyle \frac{1}{V_{tb}V_{tq}^{}}}{\displaystyle \frac{4\alpha _s}{3\alpha }}{\displaystyle \frac{m_l\mathrm{tan}^2\beta }{m_b(M_H^2M_W^2)M_{\stackrel{~}{g}}}}{\displaystyle \underset{a,b=1}{\overset{6}{}}}(\mathrm{\Gamma }^{D_L})_{qb}(\mathrm{\Gamma }^{D_R})_{a3}`$ (34)
$`\times `$ $`\{D_2(x_a,x_b)[\mathrm{\Gamma }^{D_L}M_DA_D^{}\mathrm{\Gamma }^{D_R}\pm \mathrm{H}.\mathrm{c}.]_{ba}M_{\stackrel{~}{g}}^2D_3(x_a)\delta _{ab}]\},`$ (35)
$`c_{S,P}^{\stackrel{~}{g}}`$ $`=`$ $`\pm {\displaystyle \frac{1}{V_{tb}V_{tq}^{}}}{\displaystyle \frac{4\alpha _s}{3\alpha }}{\displaystyle \frac{m_l\mathrm{tan}^2\beta }{m_q(M_H^2M_W^2)M_{\stackrel{~}{g}}}}{\displaystyle \underset{a,b=1}{\overset{6}{}}}(\mathrm{\Gamma }^{D_R})_{qb}(\mathrm{\Gamma }^{D_L})_{a3}`$ (36)
$`\times `$ $`\{D_2(x_a,x_b)[\mathrm{\Gamma }^{D_L}M_DA_D^{}\mathrm{\Gamma }^{D_R}\pm \mathrm{H}.\mathrm{c}.]_{ba}M_{\stackrel{~}{g}}^2D_3(x_a)\delta _{ab}\},`$ (37)
where the subscript $`q`$ on the $`\mathrm{\Gamma }`$ matrices (see Appendix B for details) is equal to $`1`$ $`(2)`$ for $`\overline{B}_d`$ $`(\overline{B}_s)`$ decays. Furthermore, $`\theta _W`$ is the Weinberg angle and the functions $`D_{2,3}`$ are given in Appendix A. Observe that in Eqs. (34) and (36) we have retained the leading and subleading terms in $`\mathrm{tan}\beta `$, which might become comparable in size in some part of the MSSM parameter space. A term is called leading or subleading if it has the following structure:

$$C_i𝒪_i|_{\mathrm{leading}}\mathrm{tan}^{n+1}\beta \left(\frac{m}{M}\right)^n,C_i𝒪_i|_{\mathrm{subleading}}\mathrm{tan}^n\beta \left(\frac{m}{M}\right)^n,$$
(38)
$`n=0,1,2`$,…. Further, $`m`$ denotes lepton and light quark masses while $`M`$ stands for masses of particles that have been integrated out.

(ii) Neutralino:

$`c_{S,P}^{\stackrel{~}{\chi }^0}`$ $`=`$ $`\pm {\displaystyle \frac{1}{V_{tb}V_{tq}^{}}}{\displaystyle \frac{m_l}{12M_W(M_H^2M_W^2)\mathrm{sin}^2\theta _W}}{\displaystyle \underset{k=1}{\overset{4}{}}}{\displaystyle \underset{a=1}{\overset{6}{}}}M_{\stackrel{~}{\chi }_k^0}D_3(x_{ak})`$ (39)
$`\times `$ $`\{\mathrm{tan}^4\beta {\displaystyle \frac{3m_q}{M_W}}N_{k3}^2(\mathrm{\Gamma }^{D_R})_{qa}(\mathrm{\Gamma }^{D_L})_{a3}`$ (42)
$`+\mathrm{tan}^3\beta N_{k3}\left[(\mathrm{tan}\theta _WN_{k1}3N_{k2})(\mathrm{\Gamma }^{D_L})_{qa}(\mathrm{\Gamma }^{D_L})_{a3}+{\displaystyle \frac{2m_q}{m_b}}\mathrm{tan}\theta _WN_{k1}(\mathrm{\Gamma }^{D_R})_{qa}(\mathrm{\Gamma }^{D_R})_{a3}\right]`$
$`+\mathrm{tan}^2\beta {\displaystyle \frac{2M_W}{3m_b}}\mathrm{tan}\theta _WN_{k1}(\mathrm{tan}\theta _WN_{k1}3N_{k2})(\mathrm{\Gamma }^{D_L})_{qa}(\mathrm{\Gamma }^{D_R})_{a3}\},`$
$`c_{S,P}^{\stackrel{~}{\chi }^0}`$ $`=`$ $`{\displaystyle \frac{1}{V_{tb}V_{tq}^{}}}{\displaystyle \frac{m_l}{12M_W(M_H^2M_W^2)\mathrm{sin}^2\theta _W}}{\displaystyle \underset{k=1}{\overset{4}{}}}{\displaystyle \underset{a=1}{\overset{6}{}}}\{M_{\stackrel{~}{\chi }_k^0}D_3(x_{ak})`$ (43)
$`\times `$ $`\{\mathrm{tan}^4\beta {\displaystyle \frac{3m_b}{M_W}}N_{k3}^2(\mathrm{\Gamma }^{D_L})_{qa}(\mathrm{\Gamma }^{D_R})_{a3}`$ (46)
$`+\mathrm{tan}^3\beta N_{k3}^{}\left[(\mathrm{tan}\theta _WN_{k1}^{}3N_{k2}^{})(\mathrm{\Gamma }^{D_L})_{qa}(\mathrm{\Gamma }^{D_L})_{a3}+{\displaystyle \frac{2m_b}{m_q}}\mathrm{tan}\theta _WN_{k1}^{}(\mathrm{\Gamma }^{D_R})_{qa}(\mathrm{\Gamma }^{D_R})_{a3}\right]`$
$`+\mathrm{tan}^2\beta {\displaystyle \frac{2M_W}{3m_q}}\mathrm{tan}\theta _WN_{k1}^{}(\mathrm{tan}\theta _WN_{k1}^{}3N_{k2}^{})(\mathrm{\Gamma }^{D_R})_{qa}(\mathrm{\Gamma }^{D_L})_{a3}\},`$
where $`N`$ is the neutralino mixing matrix, defined in Appendix B. Unlike the gluino contributions, we have kept only the leading term in $`\mathrm{tan}\beta `$ \[cf. Eq. (38)\], which stems from the counterterm of the electroweak wave function renormalization (see Fig. 1). However, we have checked numerically whether the subleading term is important.

## V Numerical analysis

### A Experimental constraints

Currently available data on rare $`B`$ decays already constrain new-physics contributions to the various Wilson coefficients governing the $`bs`$ transitions.A recent analysis of exclusive and inclusive rare $`B`$ decays within SUSY has been performed in Ref.  by using the SM operator basis. Recall that in the present work we utilize an extended operator basis. We take into account the following experimental upper limits :

$$(\overline{B}X_se^+e^{})<10.1\times 10^6(90\%\mathrm{C}.\mathrm{L}.),$$
(47)
$$(\overline{B}X_s\mu ^+\mu ^{})<19.1\times 10^6(90\%\mathrm{C}.\mathrm{L}.),$$
(48)
$$(\overline{B}K^{}\mu ^+\mu ^{})<3.1\times 10^6(90\%\mathrm{C}.\mathrm{L}.),$$
(49)
besides those on $`\overline{B}_{d,s}\mu ^+\mu ^{}`$ decays given in Eqs. (2). As far as the measured branching fractions of $`\overline{B}K\mu ^+\mu ^{}`$  and $`\overline{B}X_s\gamma `$  are concerned, we will allow the following ranges:

$$0.5\times 10^6(\overline{B}K\mu ^+\mu ^{})1.5\times 10^6,$$
(50)
$$2.0\times 10^4(\overline{B}X_s\gamma )4.5\times 10^4.$$
(51)
The constraints on the flavour-changing entries in the matrices $`M_{\stackrel{~}{U}_L}^2`$ and $`M_{\stackrel{~}{D}_L}^2`$ will be discussed below.

### B Implications of $`\overline{\mathrm{𝐌𝐅𝐕}}`$ for the decays $`\overline{𝑩}_{𝒅\mathbf{,}𝒔}\mathbf{}𝝁^\mathbf{+}𝝁^{\mathbf{}}`$

In our analysis, we will perform a scan over the following ranges of MSSM parameters:

$$125\text{GeV}M_H500\text{GeV},$$
(53)
$$100\text{GeV}M_{1,2}1\text{TeV},250\text{GeV}M_{\stackrel{~}{g}}1\text{TeV},$$
(54)
$$1\text{TeV}\{\mu ,(A_U)_{ii},(A_D)_{ii}\}1\text{TeV},$$
(55)
$$300\text{GeV}\{(M_{\stackrel{~}{U}_L})_{ii},(M_{\stackrel{~}{D}_L})_{ii},(M_{\stackrel{~}{U}_R})_{ii},(M_{\stackrel{~}{D}_R})_{ii}\}1\text{TeV}.$$
(56)
Further, we fix $`\mathrm{tan}\beta =50`$ and $`m_{\stackrel{~}{l},\stackrel{~}{\nu }}=100\text{GeV}`$, and take the lower bounds on the sparticle masses from . As far as the decay constants $`f_{B_q}`$ are concerned, we will take the central values as given in Eq. (27). The remaining input parameters are fixed as given in Tables 1 and 2 of Ref. .

We have analyzed more than 100,000 parameter sets, which have been produced randomly. The so-called unconstrained plots contain those parameter sets that are consistent with the lower bounds on the sparticle masses. On the other hand, in the constrained plots we have also taken into account the bounds from Eqs. (47)–(51) as well as the constraints on the flavour-changing entries in the matrices $`M_{\stackrel{~}{U}_L}^2`$ and $`M_{\stackrel{~}{D}_L}^2`$ (see discussion below).

As for the various branching ratios, it turns out that the primed Wilson coefficients are negligibly small, and hence can be safely neglected. \[Comments on this issue can be found at the end of the discussion of scenario (C).\]

### Scenario (A)

Recall that in scenario (A) the matrices $`M_{\stackrel{~}{U}_L}^2`$ and $`M_{\stackrel{~}{D}_L}^2`$ are equal and proportional to the unit matrix. Therefore, the gluinos and neutralinos do not contribute at one-loop level. In Fig. 2, we have plotted $`(\overline{B}_d\mu ^+\mu ^{})`$ versus $`(\overline{B}_s\mu ^+\mu ^{})`$, where the dashed line represents our reference curve with the slope $`R_{\mathrm{SM}}`$ while the solid line is our result obtained within scenario (A). The scan over the parameter region given in Eqs. (V B) shows that the ratio $`R`$ is approximately constant and close to $`R_{\mathrm{SM}}0.03`$.

### Scenario (B)

In this scenario, the matrix $`M_{\stackrel{~}{D}_L}^2`$ is diagonal,

$$M_{\stackrel{~}{D}_L}^2=\mathrm{diag}(m_{\stackrel{~}{d}_L}^2,m_{\stackrel{~}{s}_L}^2,m_{\stackrel{~}{b}_L}^2),$$
(57)
with at least two different entries; hence there are no gluino and neutralino contributions. Employing the relation in Eq. (17), the matrix $`M_{\stackrel{~}{U}_L}^2`$ becomes

$$M_{\stackrel{~}{U}_L}^2=V_{\text{CKM}}M_{\stackrel{~}{D}_L}^2V_{\text{CKM}}^{}=\left(\begin{array}{ccc}m_{\stackrel{~}{u}_L}^2& \mathrm{\Delta }_{\stackrel{~}{U}_L}^{12}& \mathrm{\Delta }_{\stackrel{~}{U}_L}^{13}\\ \mathrm{\Delta }_{\stackrel{~}{U}_L}^{12}& m_{\stackrel{~}{c}_L}^2& \mathrm{\Delta }_{\stackrel{~}{U}_L}^{23}\\ \mathrm{\Delta }_{\stackrel{~}{U}_L}^{13}& \mathrm{\Delta }_{\stackrel{~}{U}_L}^{23}& m_{\stackrel{~}{t}_L}^2\end{array}\right),$$
(58)
with the off-diagonal elements

$$\mathrm{\Delta }_{\stackrel{~}{U}_L}^{12}=(m_{\stackrel{~}{d}_L}^2m_{\stackrel{~}{s}_L}^2)V_{ud}V_{cd}^{}+(m_{\stackrel{~}{b}_L}^2m_{\stackrel{~}{s}_L}^2)V_{ub}V_{cb}^{},$$
(60)
$$\mathrm{\Delta }_{\stackrel{~}{U}_L}^{13}=(m_{\stackrel{~}{d}_L}^2m_{\stackrel{~}{s}_L}^2)V_{ud}V_{td}^{}+(m_{\stackrel{~}{b}_L}^2m_{\stackrel{~}{s}_L}^2)V_{ub}V_{tb}^{},$$
(61)
$$\mathrm{\Delta }_{\stackrel{~}{U}_L}^{23}=(m_{\stackrel{~}{d}_L}^2m_{\stackrel{~}{s}_L}^2)V_{cd}V_{td}^{}+(m_{\stackrel{~}{b}_L}^2m_{\stackrel{~}{s}_L}^2)V_{cb}V_{tb}^{}.$$
(62)
In writing these equations, we have used the unitarity of the CKM matrix. The flavour-changing entries given above, together with the corresponding elements in the down squark sector \[see Eqs. (64) below\], are constrained by experimental data on $`K^0`$$`\overline{K}^0`$, $`B^0`$$`\overline{B}^0`$, $`D^0`$$`\overline{D}^0`$ oscillations, and the $`bs\gamma `$ decay .\**\**\**Strictly speaking, these constraints are only valid if all squark masses are close in size, but the order of magnitude should also be valid for non-degenerate masses .

In Fig. 3, we have plotted $`(\overline{B}_d\mu ^+\mu ^{})`$ versus $`(\overline{B}_s\mu ^+\mu ^{})`$ for unconstrained and constrained parameter sets. In case of the former, the deviation of $`R`$ from $`R_{\mathrm{SM}}`$ can be one order of magnitude, with the range $`0.002R0.115`$. Applying the constraints on the flavour-changing entries in the matrix $`M_{\stackrel{~}{U}_L}^2`$, as well as the bounds given in Eqs. (47)–(51), our numerical predictions for $`R`$ in scenario (B) are similar to those in scenario (A). In fact, $`R`$ varies between $`0.026`$ and $`0.030`$. It is important to note that the bounds on $`\mathrm{\Delta }_{\stackrel{~}{U}_L}^{ij}`$  severely constrain the additional chargino contributions in scenario (B) in contrast to scenario (A).

### C Scenario (C)

We now repeat the analysis of the $`\overline{B}_{d,s}\mu ^+\mu ^{}`$ decays within the framework of scenario (C). In this case, the matrix $`M_{\stackrel{~}{U}_L}^2`$ is diagonal,

$$M_{\stackrel{~}{U}_L}^2=\mathrm{diag}(m_{\stackrel{~}{u}_L}^2,m_{\stackrel{~}{c}_L}^2,m_{\stackrel{~}{t}_L}^2),$$
(63)
with at least two different entries. According to the relation in Eq. (17), this implies that $`M_{\stackrel{~}{D}_L}^2`$ has non-diagonal entries, so that gluinos and neutralinos contribute to the $`bql^+l^{}`$ transition already at one-loop level. In this case, $`M_{\stackrel{~}{D}_L}^2`$ can be written as

$$M_{\stackrel{~}{D}_L}^2=V_{\text{CKM}}^{}M_{\stackrel{~}{U}_L}^2V_{\text{CKM}}=\left(\begin{array}{ccc}m_{\stackrel{~}{d}_L}^2& \mathrm{\Delta }_{\stackrel{~}{D}_L}^{12}& \mathrm{\Delta }_{\stackrel{~}{D}_L}^{13}\\ \mathrm{\Delta }_{\stackrel{~}{D}_L}^{12}& m_{\stackrel{~}{s}_L}^2& \mathrm{\Delta }_{\stackrel{~}{D}_L}^{23}\\ \mathrm{\Delta }_{\stackrel{~}{D}_L}^{13}& \mathrm{\Delta }_{\stackrel{~}{D}_L}^{23}& m_{\stackrel{~}{b}_L}^2\end{array}\right),$$
(64)
where the flavour-changing off-diagonal entries are given by

$$\mathrm{\Delta }_{\stackrel{~}{D}_L}^{12}=(m_{\stackrel{~}{u}_L}^2m_{\stackrel{~}{c}_L}^2)V_{us}V_{ud}^{}+(m_{\stackrel{~}{t}_L}^2m_{\stackrel{~}{c}_L}^2)V_{ts}V_{td}^{},$$
(66)
$$\mathrm{\Delta }_{\stackrel{~}{D}_L}^{13}=(m_{\stackrel{~}{u}_L}^2m_{\stackrel{~}{c}_L}^2)V_{ub}V_{ud}^{}+(m_{\stackrel{~}{t}_L}^2m_{\stackrel{~}{c}_L}^2)V_{tb}V_{td}^{},$$
(67)
$$\mathrm{\Delta }_{\stackrel{~}{D}_L}^{23}=(m_{\stackrel{~}{u}_L}^2m_{\stackrel{~}{c}_L}^2)V_{ub}V_{us}^{}+(m_{\stackrel{~}{t}_L}^2m_{\stackrel{~}{c}_L}^2)V_{tb}V_{ts}^{}.$$
(68)
As before, we take the constraints of Refs.  on these off-diagonal elements.

The scatter plots in Fig. 4 exhibit an order-of-magnitude deviation from $`R_{\mathrm{SM}}0.03`$. In the unconstrained case (left plots) $`R`$ is predicted to be in the range $`0.026R20.636`$, while for the constrained parameter sets (right plots) we find $`0.026R2.863`$. A noticeable feature of scenario (C) is that there exists a lower bound on $`R`$, i.e. $`R0.95R_{\mathrm{SM}}`$ (see upper plots of Fig. 4), which is due to the structure of the CKM matrix. We stress that this bound is valid only within scenario (C) and does not apply to scenario (B) or scenarios with new sources of flavour violation (see Sec. II).

In Fig. 5, the branching ratios of $`\overline{B}_d\mu ^+\mu ^{}`$ and $`\overline{B}_s\mu ^+\mu ^{}`$ are displayed as functions of the charged Higgs boson mass, $`M_H`$, for $`\mathrm{tan}\beta `$ equal to $`50`$ (left plot) and $`60`$ (right plot). The remaining parameters are given by

$$\mu =800\text{GeV},M_1=500\text{GeV},M_2=200\text{GeV},M_{\stackrel{~}{g}}=250\text{GeV},$$
(70)
$$A_U=A_D=\mathrm{diag}(100\text{GeV},100\text{GeV},75\text{GeV}),$$
(71)
$$M_{\stackrel{~}{U}_R}=M_{\stackrel{~}{D}_R}=\mathrm{diag}(500\text{GeV},500\text{GeV},500\text{GeV}),$$
(72)
$$M_{\stackrel{~}{U}_L}=\mathrm{diag}(460\text{GeV},500\text{GeV},517\text{GeV}).$$
(73)
Both plots are consistent with the constraints on rare $`B`$ decays in Eqs. (47)–(51), and with those on the flavour-changing entries $`\mathrm{\Delta }_{\stackrel{~}{D}_L}^{ij}`$ . Interestingly, in the left (right) plot, the ratio $`R`$ ranges between $`0.15R0.81`$ ($`0.44R1.78`$), while the magnitude of the individual branching fractions decreases drastically with increasing charged Higgs boson mass, $`M_H`$. Note that for small $`M_H`$ and $`\mathrm{tan}\beta `$ close to $`60`$ both branching ratios are in a region that can be probed experimentally, in Run II of the Fermilab Tevatron, BABAR and Belle. The large enhancement of $`R`$ is due to $`\text{Im}[c_{S,P}^{\stackrel{~}{g}}(\overline{B}_d\mu ^+\mu ^{})]`$ being by an order-of-magnitude larger than $`\text{Im}[c_{S,P}^{\stackrel{~}{g}}(\overline{B}_s\mu ^+\mu ^{})]`$ and a cancellation between the real parts of $`c_{S,P}^{\stackrel{~}{\chi }^\pm }`$ and $`c_{S,P}^{\stackrel{~}{g}}`$. As a consequence, the neutralino contributions to the scalar and pseudoscalar coefficients, although negligible compared to the corresponding contributions in the chargino and gluino sector, cannot be neglected in the $`\overline{B}_s\mu ^+\mu ^{}`$ decay. The neutralinos contribute significantly to the ratio $`R`$; for illustration, the ratio $`R0.94`$ in the right plot ($`\mathrm{tan}\beta =60`$) at $`M_H=130\text{GeV}`$ increases to $`R1.55`$ after neglecting the contributions from neutralinos.

As already mentioned, we have found that the primed Wilson coefficients, $`c_{S,P}^{}`$, are numerically suppressed. In fact, applying the mass insertion approximation as outlined in Appendix C, one can convince oneself that the enhancement factor $`1/m_q`$ (recall $`q=s,d`$) appearing in the scalar and pseudoscalar coefficients cancels out. Therefore, the primed Wilson coefficients are suppressed by a factor $`m_q/m_b1`$, compared to the unprimed coefficients.

## VI Summary and conclusions

In the framework of the minimal supersymmetric extension of the SM with large $`\mathrm{tan}\beta `$, we have computed the gluino and neutralino exchange diagrams contributing to the purely leptonic decays $`\overline{B}_{d,s}\mu ^+\mu ^{}`$. Together with our previous analysis , the present paper provides a complete study of the decays $`\overline{B}_{d,s}\mu ^+\mu ^{}`$ at one-loop level in the MSSM with modified minimal flavour violation and large $`\mathrm{tan}\beta `$.

We have defined $`\overline{\mathrm{MFV}}`$ using symmetry arguments and have shown that $`\overline{\mathrm{MFV}}`$ is less restrictive than MFV, while the CKM matrix remains the only source of flavour violation. We have given a criterion for testing whether a theory belongs to the class of $`\overline{\mathrm{MFV}}`$. Within the MSSM we have investigated three scenarios that are possible within the context of $`\overline{\mathrm{MFV}}`$. In particular, we have studied the case where the gluino and neutralino exchange diagrams contribute besides $`W^\pm ,H^\pm ,\stackrel{~}{\chi }^\pm `$ \[scenario (C)\]. The neutralino Wilson coefficients are numerically smaller than those coming from the chargino and gluino contributions. However, we have found that in certain regions of the MSSM parameter space cancellations between the chargino and gluino coefficients occur, in which case the neutralino contributions become important. As a matter of fact, for the SUSY parameter sets examined, we found that a large value of $`R(\overline{B}_d\mu ^+\mu ^{})/(\overline{B}_s\mu ^+\mu ^{})`$ always involves such a cancellation.

Including current experimental data on rare $`B`$ decays, as well as on $`K,B,D`$ meson mixing, we found that in certain regions of the SUSY parameter space the branching ratios $`(\overline{B}_d\mu ^+\mu ^{})`$ and $`(\overline{B}_s\mu ^+\mu ^{})`$ can be up to the order of $`10^7`$ and $`10^6`$ respectively. Specifically, we showed that there exist regions in which the branching fractions of both decay modes are comparable in size, and may well be accessible to Run II of the Fermilab Tevatron.

We wish to stress that a measurement of the branching ratios $`(\overline{B}_{d,s}\mu ^+\mu ^{})`$, or equivalently, a ratio $`R`$ of $`O(1)`$, does not necessarily imply the existence of new flavour violation outside the CKM matrix. Nevertheless, any observation of these decay modes in ongoing and forthcoming experiments would be an unambiguous signal of new physics.

###### Acknowledgements.

We would like to thank Andrzej J. Buras, Manuel Drees, Gino Isidori, and Janusz Rosiek for useful discussions. We are grateful to Andrzej J. Buras for his comments on the manuscript. We would also like to thank Martin Gorbahn for suggesting the name Modified Minimal Flavour Violation. This work was supported in part by the German ‘Bundesministerium für Bildung und Forschung’ under contract 05HT1WOA3 and by the ‘Deutsche Forschungsgemeinschaft’ (DFG) under contract Bu.706/1-1.

## A Auxiliary functions

The loop functions appearing in the formulae of Sec. IV are given by

$$D_2(x,y)=\frac{x\mathrm{ln}x}{(1x)(xy)}+(xy),$$
(A1)
$$D_3(x)=\frac{x\mathrm{ln}x}{1x}.$$
(A2)
## B Mass matrices

### 1 Squark mass matrices

The squark mass-squared matrices given in Eq. (7) are diagonalized according to

$$\mathrm{\Gamma }^U_U^2\mathrm{\Gamma }^U=\mathrm{diag}(m_{\stackrel{~}{u}_1}^2,\mathrm{},m_{\stackrel{~}{u}_6}^2),\mathrm{\Gamma }^D_D^2\mathrm{\Gamma }^D=\mathrm{diag}(m_{\stackrel{~}{d}_1}^2,\mathrm{},m_{\stackrel{~}{d}_6}^2).$$
(B1)
It is convenient to split the $`6\times 6`$ squark mixing matrices, $`\mathrm{\Gamma }^{U,D}`$, into two $`6\times 3`$ submatrices:

$$(\mathrm{\Gamma }^Q)_{ai}=(\mathrm{\Gamma }^{Q_L})_{ai},(\mathrm{\Gamma }^Q)_{a,i+3}=(\mathrm{\Gamma }^{Q_R})_{ai},Q=U,D,$$
(B2)
where $`a=1,\mathrm{},6`$ and $`i=1,2,3`$.

### 2 Neutralino mass matrix

The neutralino mass matrix has the structure

$$_{\stackrel{~}{\chi }^0}=\left(\begin{array}{cccc}M_1& 0& M_Z\mathrm{sin}\theta _W\mathrm{cos}\beta & M_Z\mathrm{sin}\theta _W\mathrm{sin}\beta \\ 0& M_2& M_Z\mathrm{cos}\theta _W\mathrm{cos}\beta & M_Z\mathrm{cos}\theta _W\mathrm{sin}\beta \\ M_Z\mathrm{sin}\theta _W\mathrm{cos}\beta & M_Z\mathrm{cos}\theta _W\mathrm{cos}\beta & 0& \mu \\ M_Z\mathrm{sin}\theta _W\mathrm{sin}\beta & M_Z\mathrm{cos}\theta _W\mathrm{sin}\beta & \mu & 0\end{array}\right),$$
(B3)
which is diagonalized by a unitary matrix $`N`$ such that

$$N^{}_{\stackrel{~}{\chi }^0}N^{}=\mathrm{diag}(M_{\stackrel{~}{\chi }_1^0},M_{\stackrel{~}{\chi }_2^0},M_{\stackrel{~}{\chi }_3^0},M_{\stackrel{~}{\chi }_4^0}).$$
(B4)
## C Perturbative diagonalization of $`𝓜_𝑫^\mathrm{𝟐}`$

Within the framework of scenario (C), all flavour-changing interactions and CP violation are due exclusively to the CKM matrix. In this case, the matrices $`_{D_{LR}}^2`$ and $`_{D_{RR}}^2`$ in Eq. (7) are diagonal and real whereas $`_{D_{LL}}^2`$ contains complex off-diagonal entries \[cf. Eq. (64)\]. In order to diagonalize the down squark mass-squared matrix $`_D^2`$ perturbatively, we rewrite it as

$$_D^2=\stackrel{~}{}_D^2+\stackrel{~}{\mathrm{\Delta }}_D=\left(\begin{array}{cc}\stackrel{~}{}_{D_{LL}}^2& _{D_{LR}}^2\\ _{D_{LR}}^2& _{D_{RR}}^2\end{array}\right)+\left(\begin{array}{cc}\mathrm{\Delta }_{D_{LL}}& 0\\ 0& 0\end{array}\right),$$
(C1)
where $`\stackrel{~}{}_{D_{LL}}^2`$ ($`\mathrm{\Delta }_{D_{LL}}`$) contains the diagonal (off-diagonal) elements of $`_{D_{LL}}^2`$. We remark parenthetically that for large values of $`\mathrm{tan}\beta `$ the LR elements in $`_D^2`$ cannot in general be treated as small perturbations. \[This is the reason for choosing the decomposition of $`_D^2`$ given in Eq. (C1).\]

Writing $`\mathrm{\Gamma }^D`$ as a product of two unitary matrices, $`\mathrm{\Gamma }^DXY`$, we have

$$(_D^2)_{\mathrm{diag}}=\mathrm{\Gamma }^D_D^2\mathrm{\Gamma }^D=X\left\{D+\mathrm{\Delta }_D\right\}X^{},$$
(C2)
where $`D=Y\stackrel{~}{}_D^2Y^{}\mathrm{diag}(D_1,\mathrm{},D_6)`$ and $`\mathrm{\Delta }_D=Y\stackrel{~}{\mathrm{\Delta }}_DY^{}`$. Then, if $`|(\mathrm{\Delta }_D)_{ij}|D_a`$, we can make the ansatz

$`X`$ $`=`$ $`𝟙+\delta 𝕏^{(\mathrm{𝟙})}+\delta 𝕏^{(\mathrm{𝟚})}+\mathrm{},`$ (C3)
$`(_D^2)_{\mathrm{diag}}`$ $`=`$ $`D+\delta D^{(1)}+\delta D^{(2)}+\mathrm{}`$ (C4)
to solve Eq. (C2) perturbatively in terms of $`\mathrm{\Delta }_D`$. As a result, we obtain

$`\mathrm{\Gamma }_{ac}^Df(m_{\stackrel{~}{d}_c}^2)\mathrm{\Gamma }_{cb}^D`$ $`=`$ $`Y_{ac}^{}f(D_c)Y_{cb}+Y_{ac}^{}\mathrm{\Delta }_{cd}f_1(D_c,D_d)Y_{db}`$ (C5)
$`+`$ $`Y_{ac}^{}\mathrm{\Delta }_{ce}\mathrm{\Delta }_{ed}f_2(D_c,D_d,D_e)Y_{db}+\mathrm{},`$ (C6)
where $`f`$ denotes an arbitrary loop function, and

$`f_1(x,y)={\displaystyle \frac{f(x)f(y)}{xy}},f_2(x,y,z)={\displaystyle \frac{f_1(x,z)f_1(y,z)}{xy}},`$ (C7)
for $`xyxzyz`$. In all other cases the corresponding limit of the functions $`f_{1,2}`$ has to be taken. Setting $`Y𝟙`$ in Eq. (C5) reproduces the result given in Refs. .


format: text

text:
# Critical dynamics of superconducting BSCCO films

## I Introduction

Direct observation of the scaling properties of the superconducting transition has at last become feasible due to the discovery of the cuprate superconductors. Microwave measurements in ultrapure single crystal YBa2Cu3O7-x (YBCO) have provided strong evidence that for accessible temperature ranges, the effective static universality class is the three dimensional XY model (3D-XY)kamal , where the penetration depth displayed scaling over three decades in reduced temperature. Measurements on thin films have yielded varying results: microwave measurements booth  report $`\nu =1.2`$, contactless ac conductivity measurements naki  yield $`\nu =1.7`$ at frequencies up to $`2`$ GHz and dc conductivity in finite magnetic fields robe  give $`\nu =0.91.0`$. While none of these agree with the 3D-XY scaling behavior, it is widely believed that the results on pure single crystals are indicative of the intrinsic critical fluctuations of the superconducting transition.

Transport measurements have probed the dynamics of the critical fluctuations, which are characterized by the value of the dynamic critical exponent $`z`$ describing how the relaxation time $`\tau `$ scales with the correlation length $`\xi `$: $`\tau \xi ^z`$. However for YBCO, experiments have not yielded a consistent picture. For example, longitudinal dc-resistivity measurements JTK  yield $`z=1.5\pm 0.1`$ in magnetic fields up to $`6`$T, while microwave conductivity measurements booth  are consistent with $`z=2.33.0`$. Still larger values for the dynamic exponent ranging from $`z=5.6`$ to $`z=8.3`$, reminiscent of glassy dynamics, have also been seen in resistivity measurements naki ; robe  in thin films. These large values of $`z`$ strongly suggest that disorder plays a dominant role in these thin films, and this may explain why their observed static critical behavior differs from that of the best single crystals.

Critical fluctuation results on Bi2Sr2CaCu2O8+δ (BSCCO) are even more dissimilar, partially due to highly anisotropic fluctuations. The complex conductivity of underdoped BSCCO films measured at terahertz frequencies was found to agree with a Kosterlitz-Thouless-Berezinskii (KTB) model with diffusive dynamics.corson ; minnhagen  In previous work on BSCCO crystals a non-linear transport study found 2D critical behavior with an anomalous dynamic exponent of $`z5.6`$ and no crossover to 3D behavior.pierson  The 3D universal phase angle of the complex conductivity deduced from the magnetic susceptibility of BSCCO epitaxial films and crystals has even yielded larger values of $`z`$.kotzler  Only dc fluctuation conductivity measurements in BSCCO rapp  find relaxational dynamics, $`z2.0`$, assuming a 3D-XY static exponent, $`\nu 2/3`$, in agreement with recent theoretical work that predicts that $`z2.0`$ in 3D.aji  These results indicate the lack of consensus on the true nature of the dynamic fluctuations representative of the superconducting transition, which is the primary motivation of the present investigation. While a large value of the dynamic exponent $`(z>2)`$ can be understood as disorder effects, only the presence of a previously undetected collective mode coupled to the superfluid density would lead to $`z1.5`$ in 3D, just as occurs in superfluid He4.HH

In this paper we provide results on the critical fluctuations in high quality epitaxially-grown BSCCO films. We have measured the ab-plane complex conductivity with a two-coil inductive technique at zero applied magnetic field ($`H=0`$), described in Section II. Our measurements, presented in section III, indicate that the critical fluctuations in BSCCO films are consistent with the 3D-XY critical fluctuation model, rather than a KTB transition within the layers. This contrasts with earlier studies using this technique on YBCO, which have reported KTB critical fio  and mean-field lin  fluctuations. In section IV, we estimate the temperature intervals in which possible crossovers to 2D fluctuations can occur, due either to decoupling of the bilayers or to the thin film behaving as a 2D system with a large $`c`$-axis correlation length. We find that our measurements should indeed be well-described by anisotropic 3D-XY critical fluctuations. In section V, we provide a crossover analysis of the approach to the critical point, and examine the approach to a universal phase angle. In particular, we predict the qualitative forms of the variations with temperature and frequency. These are in agreement with our measurements, but contrast earlier results which report very different values of the static and dynamic critical exponents than those reported here. Analysis of the phase angle of the film response using linear response theory in the critical region, which yields a dynamic exponent of $`z=2.0`$, is presented in section VI. However, the breakdown of linear response theory becomes important near the transition temperature $`T_c`$, so that the analysis is complicated by the need to extract the superfluid density from the measured mutual inductance by an approximate inversion technique. Our analysis shows that the determination of $`z`$ purely from linear response theory is contaminated by non-linear effects, and at probe currents that are too large the effective value of $`z`$ inferred can be as low as 1.5. To circumvent this problem, we derive in section VII a non-linear scaling law obeyed by the raw measured mutual inductance at $`T_c`$ (i.e. without the need for a data inversion) and extract the dynamic exponent directly from the mutual inductance. In agreement with the linear response analysis, we again obtain $`z2.0`$. Section VIII summarizes our conclusions.

## II Samples and Experimental Setup

We have measured several high quality oxygen-doped BSCCO films, grown by molecular-beam-epitaxy on SrTiO3 substrates and analyzed in-situ by RHEED ECK . The BSCCO films are nominally oriented along the c-axis, with unit cell thickness of $`d=15.4`$ Å. AFM images of similar BSCCO films reveal an rms roughness of 5 Angstroms primarily due to single unit cell step height variation over the surface. The rms roughness of the substrates are approximately 1.5 Angstroms.

The inductively-measured $`T_c`$ from several films of different dopings was fit to the empirical formula $`T_c=T_{c,max}(182.6(p0.16))^2`$, to determine the hole doping per Cu, p.presland  The two branches of the formula are distinguished by the relative low temperature superfluid density, which is known to increase monotonically through optimal doping (p=0.16). Films A, B, and C have $`n=21,40`$ and $`60`$ unit cells respectively and are grown on $`10`$ and $`14`$ mm square substrates. Film A is overdoped at $`p0.19`$, while film B and C are slightly overdoped and underdoped, respectively. We estimate that a doping gradient in the samples introduces a experimental temperature broadening of 0.15 K near $`T_c`$.

To obtain the ab-plane complex conductivity, we employ a two-coil inductive technique. The mutual inductance of two coaxial coils is measured with a film placed in between and normal to the axis of the coils. An ac current is applied through the drive coil; a second coil, attached to a large input impedance lock-in amplifier, measures the fields produced by the drive coil and the screening currents in the film. The coils have $`135`$ turns and an average radius of $`1.5`$ mm. The drive coil is placed $`0.35`$ mm above the film so that at a current of $`I=50\mu `$A rms, the applied ac magnetic field is $`<0.1`$ Gauss normal to the film. Measurements are performed in a He4 continuous-flow dewar, which allows the temperature of the film to be controlled without significant heating from the drive coil. The sample temperature is monitored by a Si diode attached to the back side of the substrate.

The complex conductivity $`\sigma =\sigma _1i\sigma _2`$ of each film is calculated from the mutual inductance using the exact analytical expression for an infinite diameter film clem  and a numerical inversion technique similar to that of Fiory et al. fiory . The method of Turneaure et al. turneaure  is employed to correct for film diameter effects prior to calculating the in-plane complex conductivity. The measured penetration depth $`\lambda _{ab}=(\mu _0\omega \sigma _2)^{1/2}`$, is independent of the measurement frequency from $`f=10100`$ kHz below the transition temperature, where the dissipation is small ($`\sigma _1\sigma _2)`$. To analyze the critical fluctuations, we study the complex superfluid density per CuO bilayer, $`\rho =i\sigma \omega d\mathrm{\Phi }_0^2/(4\pi ^2k_B)`$, where $`\mathrm{\Phi }_0`$ is the superconducting flux quantum.FFH  The real part of this quantity, $`Re[\rho ]`$ has units of temperature and sets the energy scale for critical fluctuations. In Fig. 1(a), $`Re[\rho ]`$ obtained at $`f=80`$ kHz and $`I=40\mu `$A is plotted for films A through C. At $`5`$ K, films A, B, and C, have penetration depths of $`235`$, $`265`$, and $`285`$ nm, respectively. Note that film A, which is overdoped, has a smaller low temperature penetration depth than film B and C. All films exhibit a $`T^2`$ component to the superfluid density at low temperatures, however the zero-temperature intercept of the films near optimal doping is approximately zero, indicating the crystalline films have very few impurities.

## III Static critical behavior

In Figs. 1(a) and 1(b) a power of $`\rho =Re[\rho ]`$ and $`2T/\pi `$ is shown to locate where we would expect the KTB transition for a noninteracting bilayer, $`\rho =2T/\pi `$.bea  We observe no such drop in the superfluid density at this temperature in our films. The films instead show that coupling among the bilayers induces an anisotropic 3D transition, precluding a KTB transition in the bilayers. At the transition temperature of an isolated bilayer $`T=\rho \pi /2`$, however, the measured superfluid density has been renormalized from the mean field value by quasi-2D fluctuations.

The mean field behavior, expected at low temperatures, crosses over to the 3D-XY behavior as the critical temperature is approached from below. In layered systems, the 3D-XY behavior is exhibited once the c-axis correlation length exceeds the interlayer spacing, $`d`$; i.e. when $`\rho (T)=c_sT_c`$, where $`c_s0.5`$ and $`\rho =Re[\rho ]`$.kamal  In Fig. 1(b), the data from Fig. 1(a) are plotted as $`Re[\rho ]^{3/2}`$ near the critical regime. Also shown in Figs. 1(a) and 1(b) is the corresponding power of $`c_sT_c`$ to locate the onset of 3D-XY behavior. The static 3D-XY model gives $`\rho \xi ^1|TT_c|^\nu `$, where $`\nu 2/3`$. We observe that $`Re[\rho ]^{3/2}`$ varies linearly in temperature for $`\rho (T)c_sT_c`$ indicating a 3D-XY, rather than mean-field ($`\rho |TT_c|`$), behavior. At our measurement current ($`40\mu `$A), we observed an absence of frequency dependence below the dissipation peak over the temperature range of the fit. In film C, the linear regime extends as far as 7 K below $`82`$ K. The linear fits extrapolate to $`T_c=79.1`$ K , $`84.9`$ K and $`84.7`$ K for films A, B and C respectively.

Although our static scaling analysis is carried out at temperatures and at a drive current for which there is no discernible frequency dependence, at higher temperatures we observe variations with both drive current and frequency. This is apparent in Fig. 2, for which larger currents are used on film B with 40 layers. Film C, with 60 layers, exhibited a weaker dependence over the same frequencies and currents. Since the superfluid density depends on frequency and current at higher temperatures, the superfluid density scales as $`\rho \xi ^{2D}f(\omega \xi ^z,E\xi ^{1+z})`$, where E is the amplitude of the electric field. Only in the critical regime where the data are independent of both frequency and electric field can one expect to see the linear response scaling behavior and extract the exponent $`\nu `$. At higher temperatures where we observe non-linear behavior as well as frequency dependence, the data are still consistent with the 3D-XY fluctuations, and we can extract $`z`$. In section VI we analyze the phase angle in terms of linear response theory, but a nonlinear scaling analysis is required for completeness and is given in section VII.

## IV Dimensionality of the static fluctuations

When the superfluid density is a fraction of the 3D crossover ($`\rho (T)<(1/2)c_sT_c`$), we observe a drop in the superfluid density from the static scaling value. Although this could be interpreted as a crossover to 2D behavior, we do not believe this to be the case in our films since the drop is observed to be current as well as frequency dependent (see Fig. 2). Notice that the superfluid density does drop more rapidly than expected from the 3D-XY model as the frequency is decreased for a given drive current, but does so more slowly as the current is decreased at fixed frequency. While the former might lead one to conclude that the anisotropic 3D-XY behavior is no longer valid in the static limit, the latter would imply the opposite. More accurate measurements are needed at low frequency and currents to establish unambiguously the true nature of the static limit. Nevertheless the trend with decreasing current, coupled with the phase angle dependence discussed in the next section, is consistent with 3D-XY behavior. Another putative crossover to 2D behavior is expected when $`\xi _c`$ becomes as large as the film thickness, $`h`$. The 3D-XY scaling implies that this occurs at a temperature, $`T^{}`$, which is $`0.05`$ K, $`0.01`$ K and $`0.02`$ K below the transition temperature for films A, B and C respectively. This implies that the region of pure 2D fluctuations is too small to be resolved in our data and therefore we expect 3D fluctuations to dominate.

## V Crossover phenomena in the linear phase angle near $`T_c`$

Since the experiments are performed at finite frequency, it is important to consider whether or not the correct asymptotic behavior is being probed and what conclusions can be drawn from the observed behavior. We begin with the superfluid density, within linear response, written as

$$\rho =(t^\nu )^{2D}f_1(\omega t^{\nu z})$$
(1)
where $`t=\left|TT_c\right|/T_c`$ and $`f_1`$ is a scaling function. In the limit of $`t0`$ ($`TT_c`$), the $`t`$-dependence must cancel out of the expression for $`\rho `$, which implies that the asymptotic behavior of the scaling function $`f_1(x)x^\theta `$ as $`x\mathrm{}`$ where $`\theta =(2D)/z`$. For $`t0`$, we can write

$$\rho =(i\omega )^{(D2)/z}f_2(t\omega ^{1/\nu z})$$
(2)
where the scaling function $`f_2(y)=`$ const as $`y0`$. We now ask what behavior we should observe as $`\omega 0`$ for $`t0`$. For $`\omega 0`$ at $`t=0`$, we should see the critical point behavior. For $`t/\omega ^{1/\nu z}1`$ we are probing the $`y0`$ regime of the scaling function $`f_2(y)`$ and expect to see asymptotic scaling behavior. If $`t/\omega ^{1/\nu z}1`$ we are probing the $`y\mathrm{}`$ regime of the argument of $`f_2`$ and do not expect to see the scaling behavior. Thus at a finite $`t_0`$, for $`\omega t_0^{4/3}`$ we would observe critical behavior while for $`\omega t_0^{4/3}`$ we would not. We have used the critical exponents $`\nu =2/3`$ and $`z=2`$ for specificity, but the argument does not rely on these precise values. Note that, as is generally the case with crossover, the behavior of the superfluid density away from the critical point at finite frequencies is not universal except in the large frequency limit - a result that is sometimes found to be surprising, but which follows from a generic crossover analysis of relevant variables. goldenfeldbook  This predicted low frequency behavior is consistent with our data where we observe the critical phase angle in Fig. 3(b) and yet we do not observe scaling below $`T_c`$. More importantly, this is also consistent with the non-linear scaling observed at $`T_c`$ in Section VII. Away from $`T_c`$ an asymptotic scaling regime is demonstrated in the microwave measurements of Booth et al.booth

## VI Dynamic critical behavior and the breakdown of linear response theory

Next we study the phase of $`\rho `$, $`\varphi _\rho `$, close to the transition temperature. The frequency dependence of $`2\varphi _\rho /\pi `$ for film B measured at $`I=`$100 $`\mu A`$ is shown in Fig. 3(a). Within linear response theory, the phase angle is independent of frequency at the critical temperature and is given by $`\varphi _\rho (T_c)=\pi /(2z)`$ in three dimensions, where $`z`$ is the dynamic critical exponent.DOR  In Fig. 3(a) we notice that the curves for different frequencies seem to approach a phase angle consistent with a dynamic exponent of $`z1.5`$ and $`T_c=84.9`$ K. When repeated at a smaller value of the current, $`I=50\mu A`$, $`\varphi _\rho (T_c)`$ is smaller, as shown. Since the constant phase angle is a result from linear response theory, the estimate of the dynamic exponent from the data is expected to be more accurate at smaller currents.

This is indeed seen in the current dependence of the phase angle at different frequencies for film C. In Fig. 3(b), $`2\varphi _\rho /\pi `$ is shown for a set of frequencies and currents. As the current is lowered, the frequency independent phase angle shifts to lower values. At $`20\mu `$A, $`10\mu `$A and $`5\mu `$A, this phase angle is $`0.57\pi /2`$, $`0.54\pi /2`$ and $`0.51\pi /2`$ respectively. This suggests that $`z2.0`$ in film C and the result of $`z1.5`$ from film B is an artifact of the large current used in the measurement. With this analysis, the thickest film (C) exhibits the best estimate of the critical phase angle, the film with intermediate thickness (B) provides evidence for similar behavior, but the thinnest film did not show a critical phase angle. Due to the crossover phenomena discussed in the previous section, the behavior of the phase angle away from $`T_c`$ is not universal at low frequencies, but we can still infer the dynamic critical exponent from the critical phase angle at $`T_c`$ which is indeed universal.

Both the phase angle and the superfluid density are found to be dependent on the drive current near $`T_c`$. In the limit of smaller currents, the frequency independent critical phase angle does indeed yield a dynamic exponent of $`z2.0`$. This is a result consistent with the presence of $`3D`$ fluctuations. This is also consistent with the behavior of the superfluid density with decreasing current at fixed frequency, where it is seen to approach the asymptotic static scaling law. Furthermore, the smooth evolution of the phase angle suggests that the crossover to the 2D regime near $`T^{}`$ is not resolved in our experiment. Nevertheless, considering that nonlinear scaling and dimensional crossover is expected close to $`T_c`$, we now proceed to analyze the nonlinear response to get a more accurate determination of the dynamic exponent.

## VII Nonlinear scaling at the superconducting transition

The field equation governing the vector potential, $`\stackrel{}{A}`$, is

$$^2\stackrel{}{A}/\mu _0=\stackrel{}{J_d}h\delta (z)\sigma \stackrel{}{E}$$
(3)
where $`\stackrel{}{J_d}`$ is the current density in the drive coil, $`\sigma `$ is the conductivity of the superconducting film, and $`\stackrel{}{E}`$ is the electric field. For the geometry of the setup,

$$\stackrel{}{J_d}=I_d\delta (rr_c)\delta (zz_c)\widehat{\varphi }$$
(4)
in cylindrical coordinates, where $`I_d`$ is the drive current, $`r_c`$ is the radius of the drive coil and $`z_c`$ is the distance from the film to the drive coil.clem  For an anisotropic system, in the critical regime, the superfluid density, $`\rho `$, scales as,

$$\rho i\omega \sigma \xi _{ab}^{3D}\xi _c^1f(i\omega \xi _{ab}^z,E\xi _{ab}^{1+z})$$
(5)
where $`\xi _{ab}`$ and $`\xi _c`$ are the correlation lengths parallel and perpendicular to the CuO bilayers, and $`f`$ is a scaling function. We can rewrite the field equation as

$`^2A=`$ $``$ $`\mu _0I_d\delta (rr_c)\delta (zz_c)`$
$``$ $`h\xi _{ab}^{3D}\xi _c^1f(\omega \xi _{ab}^z,A\xi _{ab})A\delta (z)`$
where we have used the scaling form for the complex conductivity. Since the scaling function is written in terms of dimensionless variables we can non dimensionalize the equation by recasting it in terms of $`A\xi _{ab},\omega \xi _{ab}^z,r/h`$ and $`I_d\xi _{ab}`$. The solution then takes the form,

$$Ai\omega \xi _{ab}^{1+z}=G(i\omega \xi _{ab}^z,i\omega \xi _{ab}^{1+z}I_d,h^2\xi _{ab}^{3D}\xi _c^1\mathrm{\Lambda }^1,r/h)$$
(7)
where $`G`$ is a function that can be determined by solving the full non-linear field equation and $`\mathrm{\Lambda }`$ is the thermal length given by $`\mathrm{\Phi }_0^2/(4\pi \mu _0k_BT)`$. So far all we have done is to use the scaling form for the superfluid density and looked for a generic solution of a differential equation by identifying all the non-dimensional variables. For a given geometry of the experiment the right hand side of Eq. 7 is a function of three variables. One can further simplify the expression if we realize that there exists a regime $`\left|TT_c\right|<T^{}`$ where we are governed by the limit where $`\xi _c=h`$. This allows us to eliminate one of the three terms as follows. Close to $`T_c`$, using the 3D-XY behavior of the correlation length and the definition of $`T^{}`$, we can rewrite $`h^2\xi _{ab}^{3D}/\xi _c`$ as,

$$h^2\xi _{ab}^{3D}/\xi _c=h\left|1+(TT^{})/(T^{}T_c)\right|^\nu .$$
(8)
For a temperature range $`\left|TT^{}\right|\left|T^{}T_c\right|`$, $`h^2/\xi _ch`$. Thus, for temperatures within $`0.05K`$, $`0.02K`$, and $`0.01`$K of $`T^{}`$ for films a, b, and c respectively, this approximation is valid, and the field in Eq. 7 is a function of only two variables for a given geometry. Thus one can now perform the standard data collapse at $`T_c`$ of the measured mutual inductance to obtain the dynamic exponent. In this regime we can eliminate $`\xi _{ab}`$ from Eqn.7 so that

$$A=\omega ^{1/z}\stackrel{~}{G}(I_d\omega ^{1/z},h\mathrm{\Lambda }^1,r/h).$$
(9)
The mutual inductance, $`M`$, is related to the vector potential at the pick up coil and scales as

$$MA\omega ^{1/z}\stackrel{~}{G}(I_di\omega ^{1/z}).$$
(10)
Within the resolution of our data, any data collapse observed is a 3D phenomena starting to crossover to the asymptotic 2D scaling. The reason for this is that the solution (Eq. 7) depends on a number of variables, but sufficiently close to $`T_c`$ there exists a regime wherein the mutual inductance is only a function of one scaling variable. If we do not observe any data collapse we would conclude that we are never in the regime where the approximation $`h^2/\xi _ch`$ is valid. It is indeed true that the elimination of the dependence of $`\xi _{ab}`$ in Eq. 7 is possible only at $`T_c`$ and need not hold in the entire regime where the crossover occurs. It is tempting to interpret the experimental results in terms of a 2D system. Our data on the other hand are unable to resolve the crossover regime, and yet do exhibit data collapse (see Fig 4). This suggests that the temperature at which we observe data collapse lies in the crossover regime, but is not necessarily the true $`T_c`$. From our analysis we conclude that the $`T_c`$ is within $`0.02`$ K of the temperature at which we observe data collapse for film B and C. We now look for data collapse of our data measured by varying frequency, current and temperature. Notice that the value of the exponent $`\nu `$ is required only to establish the regime of validity of the approximation and is not necessary to extract the dynamic exponent from data collapse at $`T_c`$.

For the scaling analysis at $`T_c`$, the mutual inductance at fixed frequency, $`M_f`$, was compared to the mutual inductance at fixed current, $`M_I`$. For convenience, $`|M(T)|`$ is normalized to unity above the $`T_c`$. Values of $`M_f`$ ($`M_I`$), were measured at 12.5 $`\mu `$A (10 kHz) intervals for a set of fixed temperatures. Then values of $`|M_f|`$ were compared with values of $`|M_I|`$ taken at the same temperature by fitting the raw data of $`|M_f|`$ to a polynomial and selecting points with the values of $`I\omega ^{1/z}`$ equivalent to those of $`|M_I|`$. In Fig. 4(a) $`|M_f|`$ and $`|M_I|`$ for film B are plotted versus $`I\omega ^{1/2}`$ for a series of temperatures. The best data collapse of the curves appear at $`z=2.0`$ and $`T=85.22`$ K. The error in the scaling, $`\mathrm{\Delta }_M`$, is shown in Fig. 4(b) for different temperatures and values of $`z`$. This is obtained from $`\mathrm{\Delta }_M^2=\left(\left|M_f\right|\left|M_I\right|\right)^2`$, where the sum runs over the 10 values of $`I\omega ^{1/z}`$ measured for $`|M_I|`$. The lowest error value indicates $`z=2.0\pm 0.1`$ and $`T=85.2\pm 0.1`$ K. A similar analysis for the film C yields $`z=1.8\pm 0.2`$ and $`T=85.9\pm 0.2`$ K, which is also consistent with $`z=2.0\pm 0.1`$. The transition temperature of film B is $`T_c=85.2\pm 0.1`$K, while for film C it is $`T_c=85.9\pm 0.2`$ K. These values of $`T_c`$ are close to those obtained with the phase angle measurements.

In Fig. 4(b), the evolution of the error for a fixed trial $`z`$ also suggests that the dynamic exponent observed is indeed due to $`3D`$ fluctuations. Notice that as the temperature is increased for $`z=2.0`$, the error decreases, approaching the best fit (i.e. data collapse), and then increases again. Within the resolution of the data, we never observe any change to the $`2D`$ regime. In other words, the approximation that allowed us to look for the data collapse is indeed supported by the behavior of the mutual inductance.

## VIII Conclusions

The superfluid density in BSCCO films with different thicknesses has been measured and shown to be consistent with the predictions of the 3D-XY model. The dynamic scaling exponent has been obtained by performing both a linear and non-linear scaling analysis. While the linear scaling analysis did yield a 3D dynamic exponent of $`z2.0`$, non-linear effects needed to be studied to get an accurate determination. This is evident from the current dependence of the critical phase angle of the superfluid density. A 3D non-linear scaling analysis also yields a dynamic exponent of $`z=2.0\pm 0.1`$.

###### Acknowledgements.

We thank H. Westfahl and D. Sheehy for many helpful discussions. This work was supported by the U. S. Department of Energy, Division of Materials Sciences, grant DEFG02-91ER45439 through the Frederick Seitz Materials Research Laboratory at the University of Illinois at Urbana-Champaign, and by the National Science Foundation under grant number NSF-DMR-99-70690.


text:
# Evolution of percolating force chains in compressed granular media

## Abstract

The evolution of effective force chains percolating through a compressed granular system is investigated. We performed experiments by compressing an ensemble of spherical particles in a cylindrical container monitoring the macroscopic constitutive behavior and the acoustic signals emitted by microscopic rearrangements of particles. As a novel approach, we applied the continuous damage model of fiber bundles to describe the evolution of the array of force chains during the loading process. The model provides a nonlinear constitutive behavior in good quantitative agreement with the experimental results. For a system of hard particles the model predicts a universal power law divergence of stress when approaching a critical deformation. The amplitude distribution of acoustic signals was found experimentally to follow a power law with exponent $`\delta =1.15\pm 0.05`$ which is in a good agreement with the analytic solution of the model.

Recently, the behavior of granular materials has been extensively studied under various conditions due to their scientific and technological importance. Huge experimental and theoretical efforts have been devoted to obtain a better understanding of the global behavior of granular media in terms of microscopic phenomena which occur at the level of discrete particles . Subjecting a confined granular packing to an uniaxial compression a rather peculiar constitutive behavior can be observed: for small strains a strong deviation from the linear elastic response can be found implying that the system drastically hardens in this regime . Linear elastic behavior can only be achieved asymptotically at larger deformations when the system gets highly compacted. When the external load is decreased again the system shows an irreversible increase in its effective stiffness, furthermore, under cyclic loading hysteretic behavior is obtained.

Microscopically, inside a compressed granular packing, stresses are transferred by the contact of particles. Under gradual loading conditions the particles get slightly displaced changing their contacts and the local load supported by them, which can be experimentally visualized using photoelastic materials for grains. These experiments revealed that in a compressed granular system the stresses are transmitted along the direction of the external load by force chains which can branch at the grains and form a complex network . Particles lying between lines of the force network do not support any load and can even be removed from the packing without changing its mechanical properties. Increasing the external load, more and more force lines appear and they all undergo erratic changes until the system reaches a saturated state when all the particles hold typically the same load and the system behaves as a bulk material. The creation and restructuring of percolating force chains implies relative displacements of particles which can be followed experimentally by recording the acoustic waves emitted, however, up to now no such experiments have been performed systematically. Theoretically, this problem has been mainly studied by means contact dynamics simulations using spherical or cylindrical particles, and cellular automata . Computer simulations also revealed the generation and evolution of force chains in compressed granular materials, however, the statistics of microscopic restructuring events, the emergence of the array of force chains and their relation to the macroscopic constitutive behavior remained unclarified.

In the present letter the generation and evolution of percolating force chains is studied experimentally and theoretically in granular packings subjected to an uniaxial external load. We measured experimentally the macroscopic constitutive behavior and the acoustic signals emitted by microscopic restructuring events compressing an ensemble of spherical glass bead confined in a cylinder. Based on the analogy of force lines percolating through the system and fibers of a fiber composite we propose a novel theoretical approach, namely, an inversion of the Continuous Damage Model (CDM) of fiber bundles  to describe the stress transmission through granular assemblies. The model naturally captures the emergence and gradual hardening of force chains and provides analytic solutions for the constitutive behavior and acoustic activity as well. We also propose an efficient Monte Carlo simulation technique.

In our experiments a cylindrical container made of PMMA was filled with glass beads of 5 mm diameter and water. The cylinder has a thickness of 5 mm and a diameter of 140 mm. A punch test was carried out applying monotonically increasing displacements at the top level of the glass beads. Eight acoustic sensors were placed at the container wall to record the signals emitted during the compression of the beads, as can be seen in Fig. 1. The actual force and the displacement, measured at the traverse of the loading machine, and the acoustic signals were recorded simultaneously. An eight channel transient recorder was used as an analogue-digital converter to enable the storing of the acoustic emission waveforms and a signal-based data. Experiments were performed under strain controlled conditions at a fixed strain rate, i.e. moving the traverse at a constant speed of 1 mm/minute.

The nonlinear elastic response of the system can be observed in Fig. 2 where the measured force $`F`$ is presented as a function of strain $`\epsilon `$ imposed. A power law of an exponent 2.6 was obtained as a best fit to the measured data in reasonable agreement with former experiments of Ref. .

The emergence and gradual hardening of force chains is responsible on the microscopic level for the strong non-linearity observed macroscopically. To obtain information about microscopic processes, the acoustic waves emitted due to sudden relative displacements of particles were monitored. Typically several hundred signals were recorded during the experiment. The inset of Fig. 3 shows the automatically extracted peak amplitudes of the burst signals versus time. The energy is defined as the integral of the acoustic emission signal amplitude following the onset time. The energy values of the acoustic emissions are summed up in intervals of 30 seconds to elucidate the time dependent evaluation of acoustic emission activity. More details of acoustic emission data analysis and especially signal-based techniques can be found in . The statistics of restructuring events is characterized by the distribution $`D(s)`$ of the height of peaks $`s`$, which is presented in Fig. 3 on a double logarithmic plot. It can be seen that $`D(s)`$ shows a power law behavior over two orders of magnitude, the exponent of the fitted straight line is $`\delta =1.15\pm 0.05`$.

Based on the analogy of force chains in granular packings and fibers of fibrous materials or fiber reinforced composites, we propose a novel theoretical approach, namely, an inverted fiber bundle model, to describe the gradual emergence and hardening of force chains occurring under uniaxial compression. Fiber bundles are composed of parallel fibers of identical elastic properties but stochastically distributed breaking thresholds. A fiber fails during the loading process when the local load on it exceeds its breaking threshold. Fiber failures are followed by a redistribution of load on the remaining intact fibers according to the range of interaction in the system. The so-called Continuous Damage Model introduced recently , is particularly suited to model granular materials since it captures gradual stiffness changes of elements of the model. In our model of compressed packings, force lines formed by particles are represented by an array of lines organized in a square lattice as illustrated in the inset of Fig. 1. A randomly distributed rearrangement thresholds $`d`$ is assigned to each line of the array from a cumulative probability distribution $`P_0(\frac{d}{d_c})`$, where $`d_c`$ denotes the characteristic strength of force lines. During the compression process, when the local load on a line exceeds its threshold value $`d`$ the line undergoes a sudden restructuring as a result of which it becomes stiffer and straighter. The lines’ stiffness increases in a multiplicative manner, i.e. the stiffness is multiplied by a factor $`a>1`$ at each restructuring so that the constitutive equation of a single line after suffering $`k`$ restructurings reads as $`\sigma =E_\mathrm{o}a^k\epsilon ^\alpha `$. Here $`E_\mathrm{o}`$ denotes the stiffness characterizing single particle contacts, and the exponent $`\alpha `$ takes into account possible non-linearities of a single contact like for Hertz law $`\alpha =1.5`$. After each rearrangement the force chain gets a new threshold value (annealed disorder) from a distribution of the same functional form, but the characteristic strength $`d_c`$ of the distribution is increased in a multiplicative way so that after $`k`$ rearrangement events the disorder distribution takes the form $`P_k(\frac{d}{d_c})=P_0(\frac{d}{d_oq^k})`$. The maximum value of possible restructurings $`k_{max}`$ is proportional to the number of contacts and therefore to $`\frac{h}{l_o}`$ in the nomenclature of the inset of Fig. 1. The ratio $`\tau =a/q`$ is a very important parameter of the model, it decides whether the force chain becomes more fragile ($`\tau >1`$) or more ductile ($`\tau <1`$) as a result of restructuring. Essentially, $`\tau `$ is material dependent and determined by the stability of contacts. Calculations showed that the most interesting behavior of the model can be found when $`\tau =1`$, hence, in the present letter our investigation is restricted to this case. The full complexity of the model will be explored in a forthcoming publication .

Experiments and discrete element simulations  have revealed that the number of effective force chains increases during the compression process until it reaches a saturation value. To capture this effect in our model, for the number of elements we prescribe the form $`N(\epsilon )=N_oG(\epsilon )`$, where $`N_o`$ denotes the saturation number of chains, and the profile $`G(\epsilon )`$ has the property $`G(\epsilon )1`$ with increasing $`\epsilon `$. Hence, the number of force lines $`dN`$ emerging due to an infinitesimal deformation increment from $`\epsilon `$ to $`\epsilon +d\epsilon `$ is $`dN=N_og(\epsilon )d\epsilon `$, where $`g(\epsilon )=dG(\epsilon )/d\epsilon `$. Following the derivation of the constitutive behavior of the continuous damage model of fiber bundles , the macroscopic constitutive equation of the compressed granular system can be cast into the form

$`\sigma (\epsilon )`$ $`=`$ $`E_o{\displaystyle \underset{i=0}{\overset{k_{max}1}{}}}{\displaystyle _0^\epsilon }a^i(\epsilon \epsilon ^{})^\alpha g(\epsilon ^{})P_0^i(\epsilon \epsilon ^{})\times `$ (2)
$`\left[1P_0(\epsilon \epsilon ^{})\right]d\epsilon ^{}`$
$`+`$ $`{\displaystyle _0^\epsilon }a^{k_{max}}(\epsilon \epsilon ^{})^\alpha g(\epsilon ^{})P_0^{k_{max}}(\epsilon \epsilon ^{})𝑑\epsilon ^{}.`$ (3)
Eq. (2) takes also into account that the local strain of force lines $`\epsilon \epsilon ^{}`$ is different from the externally imposed strain value $`\epsilon `$ since it also depends on the initial strain $`\epsilon ^{}`$. The integral is performed over the whole loading history to take into account all the lines generated. For explicit calculations we imposed an exponential form $`N(\epsilon )=N_o(1e^{\frac{\epsilon }{\beta }})`$ for the number of chains. The best fit obtained to the experimental data is presented in Fig. 2 where the force $`F=N_o\sigma `$ is plotted against deformation $`\epsilon `$. It can be seen that a good fit was achieved with physically reasonable parameter values. The maximum possible number of percolating force chains $`N_o`$ that can emerge in the system was estimated as the ratio of the total area of the container $`A_o`$ to the cross section of a single particle $`A_p`$, i. e. $`N_o=\frac{A_o}{A_p}=784`$. The value of the other parameters are $`E_o=4600\frac{N}{m^2},d_o=3N`$, $`\beta =0.01`$ and $`a=q=1.01`$. The value of $`a`$ falls close to one indicating that a single restructuring gives rise only to a slight increase of stiffness of a force chain. Model calculations revealed that the zero derivative at the starting part of the constitutive curve is due to the gradual creation of load bearing force chains. The small value of $`\beta `$ implies that the generation of new force chains stops at a relatively small strain value, and hence, the later rapid increase of $`F`$ as function of $`\epsilon `$ is mainly caused by the hardening of the existing force lines. Further information can be gained about the constitutive behavior for larger strains by simplifying Eq. (2) assuming a fixed number $`N_o`$ of force chains from the beginning of the process. Under this assumption Eq. (2) can be reformulated as

$`\sigma (\epsilon )`$ $`=`$ $`E_o\epsilon \left[1P_0(\epsilon )\right]{\displaystyle \underset{i=0}{\overset{k_{max}1}{}}}a^iP_0(\epsilon )^i`$ (4)
$`+`$ $`a^{k_{max}}E_o\epsilon P_0(\epsilon )^{k_{max}}.`$ (5)
It can be seen from Eq. (4) that if the maximum number $`k_{max}`$ of possible restructuring events goes to infinity the stress $`\sigma `$ has finite values only for $`aP_0(\epsilon )<1`$. In this case the summation can be performed in the first term, while the second term tends to zero, and the constitutive equation takes the form

$$\sigma (\epsilon )=E_o\epsilon \left[1P_0(\epsilon )\right]\frac{1}{1aP_0(\epsilon )}.$$
(6)
It follows that the stress $`\sigma `$ diverges when $`\epsilon `$ approaches a critical value $`\epsilon _c`$, where $`\epsilon _c`$ satisfies the equation $`P_0(\epsilon _c)=1/a`$. Expanding $`P_0(\epsilon )`$ into a Taylor series at $`\epsilon _c`$ as $`P_0(\epsilon )=P_0(\epsilon _c)+p_0(\epsilon _c)(\epsilon _c\epsilon )+\mathrm{}`$, and substituting it into Eq. (6) the behavior of $`\sigma `$ in the vicinity of $`\epsilon _c`$ reads as

$$\sigma (\epsilon )\frac{1}{ap_0(\epsilon _c)(\epsilon _c\epsilon )}(\epsilon _c\epsilon )^1.$$
(7)
It means that the stress $`\sigma `$ shows a power law divergence when $`\epsilon `$ approaches the critical value $`\epsilon _c`$. The value of the exponent is universal; it does not depend on the form of disorder distribution $`P_0`$, while the value of $`\epsilon _c`$ depends on it. It is interesting to note that in Ref.  the same power law divergence was found in large scale molecular dynamic simulations of a hard sphere system.

When the number of force lines is fixed it is possible to obtain analytic results also for the statistics of restructuring events. Restructuring occurs during the compression process when the local load on a force line exceeds its threshold value. Since loading is performed under strain controlled conditions, there is no load redistribution among existing force lines, i.e. restructuring of a force line does not affect other elements of the system. If the new threshold value, assigned to the force line after rearrangement, is smaller than the local load, the force line undergoes successive restructurings until it gets stabilized. The number of steps to reach the stable state defines the size $`s`$ of the restructuring event, which is the analog of the acoustic signals measured experimentally. The number of restructuring events $`N_{k,s}`$ of size $`s`$ starting in force chains which have already suffered $`k`$ restructurings can be deduced as

$`{\displaystyle \frac{N_{k,s}(\epsilon )}{N_o}}=p_0(\epsilon )P_0^{s+k1}(\epsilon )\left[1P_0(\epsilon )\right],`$ (8)
for $`s+kk_{max}1`$, and

$`{\displaystyle \frac{N_{k,s}(\epsilon )}{N_o}}=p_0(\epsilon )P_0^{k_{max}1}(\epsilon ),`$ (9)
for $`s+k=k_{max}`$ (see also Ref. ). The number of events $`D(s)`$ of size $`s`$ can be determined by integrating over the entire loading history and summing over all possible $`k`$ values

$`D(s)={\displaystyle \underset{k=0}{\overset{k_{max}s}{}}}{\displaystyle _0^{\epsilon _c}}{\displaystyle \frac{N_{k,s}(\epsilon )}{N_o}}𝑑\epsilon +{\displaystyle _0^{\epsilon _c}}{\displaystyle \frac{N_{k_{max}s,s}(\epsilon )}{N_o}}𝑑\epsilon .`$ (10)
Finally, substituting Eqs. (8,9) into Eq. (10) and performing the calculations yields

$`D(s)=s^1,\text{where}1sk_{max},`$ (11)
i.e. the distribution of microscopic restructuring events exhibits an universal power law behavior with an exponent 1, which is completely independent on the disorder distribution. Numerical simulations revealed that the universal power law behavior also holds when the gradual creation of force chains is taken into account, i.e. when the system is described by the full Eq. (2). The statistics of restructuring events obtained by Monte Carlo simulations is presented in Fig. 4, where simulations were performed under the assumption $`N(\epsilon )=N_o(1e^{\frac{\epsilon }{\beta }})`$. The inset shows local events of different sizes that occurred during the loading process, and their distribution is presented in the main figure. The power law behavior of the analytic prediction of eq. 11 is verified. It is important to emphasize that the theoretical results on event statistics (Fig. 4) are in a very good quantitative agreement with the experimental findings (Fig. 3).

Note, that the functional form and the value of the exponent of $`D(s)`$ in the analytic calculations is mainly the consequence of the locality of restructurings due to the absence of load redistribution. The excellent agreement observed indicates that this is likely the microscopic mechanism responsible for the power law statistics of acoustic signals observed experimentally.

This work was supported by the project SFB381, and by the NATO grant PST.CLG.977311. F. Kun acknowledges financial support of the Bólyai János Fellowship of the Hungarian Academy of Sciences and of the Research Contract FKFP 0118/2001.