text:
# GRACE/SUSY Automatic Generation of Tree Amplitudes in the Minimal Supersymmetric Standard Model

## 1 Introduction

It is widely expected that the supersymmetric (SUSY) extension of the standard model(SM) of particle physics  is one of the most promising approaches that could resolve the problems remaining in the standard model. In particular, the minimal supersymmetric extension of the standard model (MSSM) has been extensively studied in the last decade, because it has the simplest structure containing the least number of particles, yet it is complex enough to allow the most essential features characteristic in any theory of SUSY.

Over these ten years we have been working with the automatic computation of high energy physics processes. The system for the computation of the SM, GRACE, has been published in . In order to be able to handle SUSY processes, we have made several necessary modifications and extension of the system . First we have installed the new interaction vertices according to the MSSM lagrangian allowing an extension to a more general effective Higgs potential (section 2.4). Second we have modified the program so that we can deal with the fermion number non-conserving vertices and Majorana particles, both of which are inherent to any SUSY theory.

Using this extended version of GRACE/SUSY, we have published a package of event-generator, SUSY23, which contains 23 specific SUSY processes for $`e^+e^{}`$ 2-body and 3-body .

Since then we have made extensive checks of the system. In this paper we present the detailed description of the program system which enables us to compute automatically the cross section of any MSSM tree process, and explain the structure of GRACE/SUSY together with instructions for how to use it.

Program pacakge GRACE/SUSY (GRACE v2.2.0) is available form
http://minami-home.kek.jp/.

## 2 Theoretical framework of the MSSM

### 2.1 Basic lagrangian

The basic lagrangian of the MSSM in its full form has been given in several papers. In this work we use the convention adopted in ref., where the positively (negatively) charged charginos are considered as fermion (anti-fermion). The full form of the MSSM lagrangian is too lengthy to be reproduced here. Readers should refer to the original work . We show only the soft SUSY-breaking terms in order to define the SUSY breaking parameters.

$`_{soft}`$ $`=`$ $`\mu {\displaystyle d^2\theta \mathrm{\Phi }_{𝐇_\mathrm{𝟏}}\mathrm{\Phi }_{𝐇_\mathrm{𝟐}}}+h.c.`$ (2.6)
$`{\displaystyle \frac{1}{2}}M_1\lambda \lambda {\displaystyle \frac{1}{2}}M_2\lambda ^a\lambda ^a{\displaystyle \frac{1}{2}}M_3\stackrel{~}{g}^\alpha \stackrel{~}{g}^\alpha +h.c.`$
$`\stackrel{~}{m}_1^2𝐇_1^{}𝐇_1\stackrel{~}{m}_2^2𝐇_2^{}𝐇_2(\stackrel{~}{m}_{12}^2𝐇_1\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)𝐇_2+h.c.)`$
$`{\displaystyle \frac{\sqrt{2}m_u}{v_2}}A_u𝐇_2𝐀(q_L)A(u_R)+{\displaystyle \frac{\sqrt{2}m_d}{v_1}}A_d𝐇_1𝐀(q_L)A(d_R)+h.c.`$
$`+{\displaystyle \frac{\sqrt{2}m_e}{v_1}}A_e𝐇_1𝐀(\mathrm{}_L)A(e_R)+h.c`$
$`{\displaystyle \underset{f}{}}[\stackrel{~}{m}_{\stackrel{~}{f}_L}^2\stackrel{~}{f}_L^{}\stackrel{~}{f}_L+\stackrel{~}{m}_{\stackrel{~}{f}_R}^2\stackrel{~}{f}_R^{}\stackrel{~}{f}_R].`$
In the first line $`\mathrm{\Phi }_{𝐇_1}`$ and $`\mathrm{\Phi }_{𝐇_2}`$ stand for the iso-doublet Higgs superfields with hypercharge $`Y=1`$ and $`Y=+1`$ respectively, and in the second line $`\lambda `$, $`\lambda ^a`$ and $`\stackrel{~}{g}^\alpha `$ are the gaugino fields for $`U(1)`$, $`SU(2)`$ and $`SU(3)`$ gauge symmetries, respectively. $`𝐇_i`$ appearing in the third to fifth line are the Higgs doublets which are decomposed as

$`𝐇_1\left(\begin{array}{c}H_1^0\\ H_1^{}\end{array}\right)`$ $`=`$ $`\left(\begin{array}{c}(v_1+\varphi _1^0i\chi _1^0)/\sqrt{2}\\ \varphi _1^{}\end{array}\right),`$ (2.7)
$`𝐇_2\left(\begin{array}{c}H_2^+\\ H_2^0\end{array}\right)`$ $`=`$ $`\left(\begin{array}{c}\varphi _2^+\\ (v_2+\varphi _2^0+i\chi _2^0)/\sqrt{2}\end{array}\right).`$ (2.8)
The scalar fields $`𝐀(q_L)`$, $`A(u_R)`$, etc. represent the sfermion fields, the superpartners of the corresponding fermion fields $`\left(\begin{array}{c}u_L\\ d_L\end{array}\right)`$ and $`\overline{u}_R`$, respectively,

$$𝐀(q_L)=\left(\begin{array}{c}\stackrel{~}{u}_L\\ \stackrel{~}{d}_L\end{array}\right),A(u_R)=\stackrel{~}{u}_R^{},A(d_R)=\stackrel{~}{d}_R^{}.$$
(2.9)
Although the so-called $`R`$-parity violating interactions  are not included in this lagrangian and hence in GRACE/SUSY, users can implement them in the system by coding the necessary vertex functions etc. in appropriate subprograms.

### 2.2 Particle content

The MSSM consists of superfields containing the known fermions or bosons as their component fields, and two Higgs superfields. Note that in contrast to the SM, the MSSM includes three neutral Higgs particles and a pair of charged Higgs particles. The mass eigenstates of the CP-even neutral Higgs states are defined by

$$\left(\begin{array}{c}H^0\\ h^0\end{array}\right)=\left(\begin{array}{cc}\mathrm{cos}\alpha & \mathrm{sin}\alpha \\ \mathrm{sin}\alpha & \mathrm{cos}\alpha \end{array}\right)\left(\begin{array}{c}\varphi _1^0\\ \varphi _2^0\end{array}\right),$$
(2.10)
where the rotation angle $`\alpha `$ is determined in such a way that the Higgs mass matrix becomes diagonal. Since any of the known particles has its superpartner, the number of the elementary particles in the MSSM is more than double the number of particles appearing in the SM. In fact, for each known fermion $`f`$ (except for neutrinos), there are two kinds of sfermions, the left-handed sfermion $`\stackrel{~}{f}_L`$ and the right-handed sfermion $`\stackrel{~}{f}_R`$. The mass eigenstates of sfermion $`\stackrel{~}{f}_i`$ with mass $`m_{\stackrel{~}{f}_i}`$ ( $`m_{\stackrel{~}{f}_1}<m_{\stackrel{~}{f}_2}`$) is defined by the orthogonal matrix

$$\left(\begin{array}{c}\stackrel{~}{f}_1\\ \stackrel{~}{f}_2\end{array}\right)=\left(\begin{array}{cc}\mathrm{cos}\theta _f& \mathrm{sin}\theta _f\\ \mathrm{sin}\theta _f& \mathrm{cos}\theta _f\end{array}\right)\left(\begin{array}{c}\stackrel{~}{f}_L\\ \stackrel{~}{f}_R\end{array}\right),f=e,u,d,\mathrm{}$$
(2.11)
where the rotation angle is determined by the mass diagonalization condition.

For each gauge boson there is a gaugino, and to each Higgs particle corresponds a higgsino. Two neutral gauginos, $`\lambda `$ and $`\lambda ^0`$, and two neutral higgsinos, $`\stackrel{~}{H}_1^0`$ and $`\stackrel{~}{H}_2^0`$, are combined to form four Majorana fermions called neutralinos, $`\stackrel{~}{\chi }_i^0\left(\begin{array}{c}\overline{\lambda }_i\\ \lambda _i\end{array}\right)`$. The two-component spinors $`\lambda _i(i=14)`$ are linear combinations of $`\lambda `$, $`\lambda ^0`$, $`\stackrel{~}{H}_1^0`$ and $`\stackrel{~}{H}_2^0`$,

$$\stackrel{~}{\lambda }_i=\eta _i(𝒪_N)_{ij}\chi _j,i=14$$
(2.12)
where $`\chi _j=(\lambda ,\lambda ^0,\stackrel{~}{H}_1^0,\stackrel{~}{H}_2^0)`$, and masses are ordered as $`m_{\stackrel{~}{\chi }_1^0}m_{\stackrel{~}{\chi }_2^0}m_{\stackrel{~}{\chi }_3^0}m_{\stackrel{~}{\chi }_4^0}`$. The sign factor $`\eta _i`$ takes the value of $`+1`$ or $`+i`$, and $`𝒪_N`$ is the orthogonal matrix which diagonalizes the mass matrix of the $`(\lambda ,\lambda ^0,\stackrel{~}{H}_1^0,\stackrel{~}{H}_2^0)`$ system.

The other neutral gauginos of the strong interaction, gluinos, remain unmixed with other particles.

Corresponding to the charged Higgs particles $`H^\pm `$ and charged gauge bosons $`W^\pm `$, there are two charged higgsinos, $`\stackrel{~}{H}_1^{}`$ and $`\stackrel{~}{H}_2^+`$, and a pair of charged gauginos (wino), $`\lambda ^\pm `$. Two charged fermions called charginos, $`\stackrel{~}{\chi }_i^+=\left(\begin{array}{c}\overline{\lambda _{iR}^{}}\\ \lambda _{iL}^+\end{array}\right)`$, are linear combination of these particles,

$`\left(\begin{array}{c}\lambda _{1R}^{}\\ \lambda _{2R}^{}\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{cos}\varphi _R& \mathrm{sin}\varphi _R\\ \mathrm{sin}\varphi _R& \mathrm{cos}\varphi _R\end{array}\right)\left(\begin{array}{c}\lambda ^{}\\ \stackrel{~}{H}_1^{}\end{array}\right),`$ (2.13)
$`\left(\begin{array}{c}\lambda _{1L}^+\\ \lambda _{2L}^+\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}1& 0\\ 0& ϵ_L\end{array}\right)\left(\begin{array}{cc}\mathrm{cos}\varphi _L& \mathrm{sin}\varphi _L\\ \mathrm{sin}\varphi _L& \mathrm{cos}\varphi _L\end{array}\right)\left(\begin{array}{c}\lambda ^+\\ \stackrel{~}{H}_2^+\end{array}\right).`$ (2.14)
The mixing angle $`\varphi _R`$, $`\varphi _L`$ and the sign factor $`ϵ_L`$ are determined so that their masses satisfy the ordering $`0m_{\stackrel{~}{\chi }_1^\pm }m_{\stackrel{~}{\chi }_2^\pm }`$ 111Note that the value of the mixing angles $`\varphi _L`$ and $`\varphi _R`$ become unstable when $`\mathrm{tan}\beta =\mu /M_2`$ or $`\mathrm{cot}\beta =\mu /M_2`$. If this happens, $`\mathrm{tan}\varphi _L`$ or $`\mathrm{tan}\varphi _R`$ become $`\pm \mathrm{}`$..

Together with the particles appearing in the SM, these SUSY-particles constitute the ingredient of the MSSM. The list of the particles and their coding names is given in Appendix A.

### 2.3 Basic parameters of the MSSM

The basic parameters of the MSSM come in when the lagrangian is expressed in terms of superfields . They are

| coupling constants | $`g,g^{},g_s,\mu `$ |
| --- | --- |
| gaugino mass parameters | $`M_1,M_2,M_3`$ |
| fermion masses (for each generation) | $`m_e,m_u,m_d`$ |
| Higgs mass parameters | $`\stackrel{~}{m}_1,\stackrel{~}{m}_2,\stackrel{~}{m}_{12}`$ |
| trilinear couplings(for each generation) | $`A_e,A_u,A_d`$ |
| sfermion mass terms(for each generation) | $`\stackrel{~}{m}_{\stackrel{~}{e}_L},\stackrel{~}{m}_{\stackrel{~}{e}_R},\stackrel{~}{m}_{\stackrel{~}{u}_L},\stackrel{~}{m}_{\stackrel{~}{u}_R},\stackrel{~}{m}_{\stackrel{~}{d}_R}`$ |

The vacuum expectation values $`v_1`$ and $`v_2`$ are determined at the minimum of the Higgs potential and they are expressed in terms of these basic parameters. Assuming massless neutrinos, there are in all 4 coupling constants, $`3N_G`$ fermion masses, $`6+8N_G`$ soft supersymmetry breaking parameters, where $`N_G`$ stands for the number of fermion generations.

### 2.4 The MSSM parameters in GRACE/SUSY

As alternatives to ten MSSM parameters for the coupling constants and gaugino and Higgs masses shown above, we can use in GRACE/SUSY the following set of ten parameters as inputs,

$$\alpha _e,\alpha _s,M_W,M_Z,\mathrm{tan}\beta ,M_{A^0},\mu ,M_1,M_2,M_3.$$
(2.15)
Unless specified explicitly by the user (see Appendix B), $`M_1`$ and $`M_3`$ are computed from $`M_2`$ by the following unification conditions

$`M_1`$ $`=`$ $`{\displaystyle \frac{5}{3}}\mathrm{tan}^2\theta _WM_2,`$ (2.16)
$`M_3`$ $`=`$ $`{\displaystyle \frac{g_s^2}{e^2}}\mathrm{sin}^2\theta _WM_2.`$ (2.17)
Concerning the eight sfermion mass parameters in each generation, we use three mixing angles, and two up-type squark masses, the lighter down-type squark mass, and two down-type slepton masses as independent input parameters. Due to the $`SU(2)_L`$ invariance of the $`_{\mathrm{soft}}`$, the heavier down-type squark mass and the sneutrino mass are given as the solution of the following two equations (shown here only for the third generation sfermions),

$`\mathrm{sin}^2\theta _bm_{\stackrel{~}{b}_2}^2`$ $`=`$ $`\mathrm{cos}^2\theta _tm_{\stackrel{~}{t}_1}^2+\mathrm{sin}^2\theta _tm_{\stackrel{~}{t}_2}^2\mathrm{cos}^2\theta _bm_{\stackrel{~}{b}_1}^2+m_b^2m_t^2M_W^2\mathrm{cos}2\beta ,`$ (2.19)
$`m_{\stackrel{~}{\nu }_\tau }^2`$ $`=`$ $`\mathrm{cos}^2\theta _\tau m_{\stackrel{~}{\tau }_1}^2+\mathrm{sin}^2\theta _\tau m_{\stackrel{~}{\tau }_2}^2m_\tau ^2+M_W^2\mathrm{cos}2\beta .`$ (2.20)
Since we restrict our consideration only to the CP conserving processes, all the parameters are real. In particular, all the physical masses used as the input are real and positive. The allowed range of the SUSY parameters used in GRACE/SUSY whose values are not yet fixed by experiments are as follows222 In ref. the value of $`\mu `$ is restricted in the range $`0\mu `$, while $`\mathrm{}<\mathrm{tan}\beta <+\mathrm{}`$ is allowed. Due to the invariance of the lagrangian with respect to an interchange $`(\mu ,\mathrm{tan}\beta )(\mu ,\mathrm{tan}\beta )`$, the physics with $`\mu =\mu _0>0`$ and $`\mathrm{tan}\beta =\mathrm{tan}\beta _0<0`$ in ref. corresponds to the one with $`\mu =\mu _0<0`$ and $`\mathrm{tan}\beta =\mathrm{tan}\beta _0>0`$ in GRACE/SUSY.

$`\mathrm{}<\mu <+\mathrm{},0\mathrm{tan}\beta ,0M_1,M_2,M_3,M_{A^0},`$ (2.21)
$`0m_{\stackrel{~}{f}_1}m_{\stackrel{~}{f}_2},{\displaystyle \frac{\pi }{2}}<\theta _f+{\displaystyle \frac{\pi }{2}}.`$ (2.22)
The CP-even Higgs, charginos and neutralinos are numbered as

$$M_{h^0}<M_{H^0},m_{\stackrel{~}{\chi }_1^\pm }<m_{\stackrel{~}{\chi }_2^\pm },m_{\stackrel{~}{\chi }_1^0}<m_{\stackrel{~}{\chi }_2^0}<m_{\stackrel{~}{\chi }_3^0}<m_{\stackrel{~}{\chi }_4^0}.$$
(2.23)
The basic parameters of the MSSM are expressed in terms of these GRACE/SUSY variables given in (2.15) as follows,

$`g`$ $`=`$ $`e\sqrt{{\displaystyle \frac{M_Z^2}{M_Z^2M_W^2}}},`$
$`g^{}`$ $`=`$ $`e\sqrt{{\displaystyle \frac{M_Z^2}{M_W^2}}},`$
$`A_f`$ $`=`$ $`\mathrm{sin}\theta _f\mathrm{cos}\theta _f{\displaystyle \frac{m_{\stackrel{~}{f}_2}^2m_{\stackrel{~}{f}_1}^2}{m_f}}\{\begin{array}{cc}\mu \mathrm{cot}\beta ,\hfill & f=u,c,t\hfill \\ \mu \mathrm{tan}\beta ,\hfill & f=d,s,b,e,\mu ,\tau \hfill \end{array},`$ (2.24)
$`\stackrel{~}{m}_{\stackrel{~}{f}_L}^2`$ $`=`$ $`\mathrm{cos}^2\theta _fm_{\stackrel{~}{f}_1}^2+\mathrm{sin}^2\theta _fm_{\stackrel{~}{f}_2}^2m_f^2M_Z^2\mathrm{cos}2\beta (T_{3f}Q_fs_W^2),`$ (2.26)
$`f=\nu _e,\nu _\mu ,\nu _\tau ,e,\mu ,\tau ,u,d,c,s,t,b,`$
$`\stackrel{~}{m}_{\stackrel{~}{f}_R}^2`$ $`=`$ $`\mathrm{sin}^2\theta _fm_{\stackrel{~}{f}_1}^2+\mathrm{cos}^2\theta _fm_{\stackrel{~}{f}_2}^2m_f^2M_Z^2\mathrm{cos}2\beta Q_fs_W^2,`$ (2.28)
$`f=e,\mu ,\tau ,u,d,c,s,t,b.`$
In the last two equations of (2.28) the identities $`\stackrel{~}{m}_{\stackrel{~}{e}_L}^2=\stackrel{~}{m}_{\stackrel{~}{\nu }_e}^2`$, $`\stackrel{~}{m}_{\stackrel{~}{u}_L}^2=\stackrel{~}{m}_{\stackrel{~}{d}_L}^2`$, etc. are guaranteed by the relations (2.19) and (2.20). The vacuum expectation values $`v_1`$ and $`v_2`$ are given by

$`v_1`$ $`=`$ $`2{\displaystyle \frac{M_W}{M_Z}}\sqrt{M_Z^2M_W^2}{\displaystyle \frac{1}{e\sqrt{1+\mathrm{tan}^2\beta }}},`$ (2.29)
$`v_2`$ $`=`$ $`2{\displaystyle \frac{M_W}{M_Z}}\sqrt{M_Z^2M_W^2}{\displaystyle \frac{\mathrm{tan}\beta }{e\sqrt{1+\mathrm{tan}^2\beta }}}.`$ (2.30)
Note that the masses of Higgs particles, neutralinos and charginos are determined by the mass diagonalization matrices introduced in (2.10), (2.12) and (2.14), respectively. In particular, the Higgs masses are given at tree level by

$`M_{H^0,h^0}^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[M_{A^0}^2+M_Z^2\pm \sqrt{(M_{A^0}^2+M_Z^2)^24M_{A^0}^2M_Z^2\mathrm{cos}^22\beta }\right],`$ (2.31)
$`M_{H^\pm }^2`$ $`=`$ $`M_{A^0}^2+M_W^2.`$ (2.32)
According to (2.32), the MSSM at Born order predicts the existence of a CP-even Higgs particle which is lighter than $`Z`$. This possibility having already been experimentally excluded, we have to make a minimal modification of the Higgs potential so that $`M_{H^0}`$, $`M_{h^0}`$ etc. can be treated as free parameters. For this purpose, inspired by the two-Higgs-doublet model , we introduce the effective Higgs potential (see also ) in the following way:

$`V^{\mathrm{eff}}`$ $`=`$ $`C_1𝐇_1^{}𝐇_1+C_2𝐇_2^{}𝐇_2+2C_3\mathrm{Re}(𝐇_1𝐇_2)+2C_4[\mathrm{Im}(𝐇_1𝐇_2)]^2`$ (2.33)
$`+C_5(𝐇_1^{}𝐇_1𝐇_2^{}𝐇_2)^2+C_6|𝐇_1^{}𝐇_2|^2+C_7(𝐇_1^{}𝐇_1)^2+C_8(𝐇_2^{}𝐇_2)^2`$
$`+4C_9(𝐇_1^{}𝐇_1)\mathrm{Re}(𝐇_1𝐇_2)+4C_{10}(𝐇_2^{}𝐇_2)\mathrm{Re}(𝐇_1𝐇_2)`$
$`+C_{11}.`$
When the coefficients $`C_i`$ take special values

$`C_1`$ $`=`$ $`\stackrel{~}{m}_1^2+\mu ^2=M_{A^0}^2\mathrm{sin}^2\beta {\displaystyle \frac{M_Z^2}{2}}\mathrm{cos}2\beta ,`$ (2.34)
$`C_2`$ $`=`$ $`\stackrel{~}{m}_2^2+\mu ^2=M_{A^0}^2\mathrm{cos}^2\beta +{\displaystyle \frac{M_Z^2}{2}}\mathrm{cos}2\beta ,`$ (2.35)
$`C_3`$ $`=`$ $`\stackrel{~}{m}_{12}^2=\mathrm{sin}\beta \mathrm{cos}\beta M_{A^0}^2,`$ (2.36)
$`C_5`$ $`=`$ $`{\displaystyle \frac{1}{8}}(g^2+g^2),`$ (2.37)
$`C_6`$ $`=`$ $`{\displaystyle \frac{g^2}{2}},`$ (2.38)
$`C_4`$ $`=`$ $`C_7=C_8=C_9=C_{10}=C_{11}=0,`$ (2.39)
the above effective Higgs potential reduces to the standard one of the MSSM. Since $`C_{11}`$ does not contribute to the Higgs masses nor to the interactions, the effective Higgs potential depends on ten parameters $`C_i,(i=110)`$, allowing more freedom than in the case of the original MSSM Higgs potential which contains three parameters, $`\stackrel{~}{m}_1^2`$, $`\stackrel{~}{m}_2^2`$ and $`\stackrel{~}{m}_{12}^2`$. The conditions at the potential minimum relate $`C_i`$ to $`v_j,(j=1,2)`$. Therefore, the effective Higgs potential contains eight free input parameters 333 Compared with the MSSM, the effective Higgs potential (2.33) contains seven more parameters and consequently we need seven more new inputs., for which we take

$$M_{A^0},M_{H^0},M_{h^0},M_{H^\pm },\mathrm{tan}\alpha ,C_4,C_9,C_{10}.$$
(2.40)
with the convention adopted in GRACE/SUSY

$$\frac{\pi }{2}\alpha 0.$$
(2.41)
In terms of our input parameters, the coefficients $`C_i`$ can be expressed as

$`C_1`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left[M_{H^0}^2+M_{h^0}^2+D\mathrm{cos}2\alpha \right]{\displaystyle \frac{v_2}{4v_1}}D\mathrm{sin}2\alpha {\displaystyle \frac{v_2}{v_1}}C_3,`$ (2.42)
$`C_2`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left[M_{H^0}^2+M_{h^0}^2D\mathrm{cos}2\alpha \right]{\displaystyle \frac{v_1}{4v_2}}D\mathrm{sin}2\alpha {\displaystyle \frac{v_1}{v_2}}C_3,`$ (2.43)
$`C_3`$ $`=`$ $`\mathrm{sin}\beta \mathrm{cos}\beta M_{A^0}^2+v_1v_2C_4v_1^2C_9v_2^2C_{10},`$ (2.44)
$`C_5`$ $`=`$ $`{\displaystyle \frac{D\mathrm{sin}2\alpha }{4v_1v_2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{M_{A^0}^2}{v_1^2+v_2^2}}{\displaystyle \frac{C_3}{v_1v_2}}+{\displaystyle \frac{3}{2}}C_4,`$ (2.45)
$`C_6`$ $`=`$ $`2({\displaystyle \frac{M_{H^\pm }^2M_{A^0}^2}{v_1^2+v_2^2}}+C_4),`$ (2.46)
$`C_7`$ $`=`$ $`{\displaystyle \frac{1}{v_1^2}}\left[C_1+{\displaystyle \frac{v_2}{v_1}}C_3+(v_1^2v_2^2)C_5+3v_1v_2C_9+{\displaystyle \frac{v_2^3}{v_1}}C_{10}\right],`$ (2.47)
$`C_8`$ $`=`$ $`{\displaystyle \frac{1}{v_2^2}}\left[C_2+{\displaystyle \frac{v_1}{v_2}}C_3+(v_2^2v_1^2)C_5+{\displaystyle \frac{v_1^3}{v_2}}C_9+3v_1v_2C_{10}\right],`$ (2.48)
where the vacuum expectation values $`v_i`$ are given by (2.30) 444 (2.47) and (2.48) are the conditions for the vanishing of the linear terms in the effective Higgs potential, thus assuring that $`v_1`$ and $`v_2`$ correspond to the minimum of the potential (2.33).  and

$$D=M_{H^0}^2M_{h^0}^2.$$
(2.49)
Combining (2.15) and (2.40), we use the following quantities as the inputs of GRACE/SUSY in the gauge and Higgs sectors,

$`\alpha _e,\alpha _s,M_W,M_Z,\mathrm{tan}\alpha ,\mathrm{tan}\beta ,M_{A^0},M_{H^0},M_{h^0},M_{H^\pm },\mu ,M_1,M_2,M_3,`$
$`C_4,C_9,C_{10}.`$ (2.50)
In addition to these quantities, all the widths of the unstable particles must also be specified as input parameters. The set of the basic input parameters used in the program system are summarized in Table 1. The magnitudes of these inputs can be changed freely by the user. A list of the names of the variables created by the system in FORTRAN codes are given in Appendix B.

A possible minimal extension of the MSSM is obtained by setting

$$C_4=C_9=C_{10}=0,$$
(2.51)
and using for $`\mathrm{tan}\alpha `$ and $`M_{H^\pm }^2`$ the same form as given by the tree MSSM,

$`\mathrm{tan}2\alpha `$ $`=`$ $`\mathrm{tan}2\beta {\displaystyle \frac{M_{A^0}^2+M_Z^2}{M_{A^0}^2M_Z^2}},`$ (2.52)
$`M_{H^\pm }^2`$ $`=`$ $`M_W^2+M_{A^0}^2,`$ (2.53)
respectively. Concerning $`M_{H^0}`$ and $`M_{h^0}`$, one can use any desired value, such as those evaluated by the one-loop MSSM.

In GRACE/SUSY, two modes are prepared for the choice of Higgs parameters. In mode 1 the strict MSSM Higgs potential is used, while in mode 2 the extended Higgs potential, (2.33), is used with the input parameters given in (2.40). The mode and the input parameters of the Higgs potential must be specified in subroutine sethptal. Note that in editing the subroutine, one can also use a different set of input parameters for the Higgs potential. See also section 4.4.2.

### 2.5 Treatment of the fermion number clashing vertices

In models of SUSY the fermion number may not be conserved at the interaction vertices when they couple to neutralinos and charginos. Note that since the fermion number of chargino is not determined by their interactions, we define the positively charged charginos as Dirac fermions.

In GRACE/SUSY we manipulate the fermion number clashing interactions by introducing the concept of G-line (GRACE-line), F-direction (fermion-direction) and M-direction (momentum-direction). The detail of the treatment for the tree case being given in ref., we present here only a brief sketch of the procedure.

Since the flow of a spin one-half particle is not disconnected nor has any branch, all the external spinors are connected pairwise. When such a spinor line connecting two external spinors is stretched in a straight line and displayed horizontally, we call it a G-line. According to the G-line rule to be discussed below, we compute the G-line amplitude corresponding to each G-line. The Feynman amplitude for each graph is then obtained by multiplying all the so-constructed G-line amplitudes, the bosonic parts which connect vertices in different G-lines and the relative phase of the graph at issue.

G-line rule
Let us assume that the process contains $`n_i`$ incoming spinors and $`n_f`$ outgoing spinors. We assign the particle number 1 to $`n_i`$ to the $`n_i`$ incoming spinor particles, and $`n_i+1`$ to $`n_i+n_f`$ to the $`n_f`$ outgoing spinor particles. For a given fermion line connecting two external spinor particles with particle number $`m`$ and $`n`$ ( $`m>n`$), we define a G-line in such a way that the particle $`m`$ comes from the left and the particle $`n`$ comes from the right endpoint. This defines the G-line uniquely for each fermion line belonging to a Feynman graph and there are $`(n_i+n_f)/2`$ such G-lines at tree level.

On each segment of the G-line separated by the interaction vertices, we assign two directions, an M-direction and an F-direction . For the notational convenience, we also introduce the notion of an A-direction, which is the direction in which spinors are ordered in the Feynman amplitude when a fermion comes in and goes out after some interactions, namely, the direction indicated by the left-oriented arrow.

M-direction
An M-direction is assigned in a natural way. For external incoming particles, their M-direction is the direction towards the vertex, while for the external outgoing particles it is the direction from the vertex. For the internal fermions, the M-direction is always taken in the A-direction (irrespective of whether the internal propagator is that of fermion or anti-fermion when seen in the A-direction).

F-direction
For Dirac particles we assign the F-direction in a natural way; for an external fermion the F-direction is the same as its M-direction, while for an external anti-fermion it is opposite to its M-direction. For the internal Dirac particle it is taken to be equal to the A-direction for fermions (when seen in the A-direction) and to the opposite direction for anti-fermions. Concerning Majorana particles which are self-charge-conjugate, how to assign their F-direction is a matter of convention. In our convention, we define the F-direction of Majorana particles by their M-direction. This means that the F-direction of the internal Majorana particle coincides with its A-direction and only the following type of the Majorana propagator appears in the G-line amplitude;

$$<\mathrm{\Psi }\overline{\mathrm{\Psi }}>==\frac{i}{/pm},$$
(2.54)
where the momentum $`p`$ runs from right to left. Clashing of the F-direction occurs only at vertices, and the charge conjugation matrix appears only at the fermion number clashing vertices.

The above assignments of the F-direction and the M-direction uniquely identify the external spinors as shown in Table 2.

The general rule of flipping one of the F-direction at the vertex is summarized as follows

$``$ $`=`$ $`C(),`$
$``$ $`=`$ $`C^1(),`$
$``$ $`=`$ $`()C^1,`$
$``$ $`=`$ $`()C,`$
where $``$ is a fermion line whose F-direction can be arbitrary. The charge conjugation matrix $`C`$ satisfies the following properties

$$C^{}=C^1,C=C^T,C[1,\gamma _5,\gamma _\mu ,\gamma _\mu \gamma _5]^TC^1=[1,\gamma _5,\gamma _\mu ,\gamma _\mu \gamma _5].$$
(2.55)
It is also important to note that we can choose the relative phase of the spinor wave functions such that

$$C\overline{v}(k,\lambda )^T=u(k,\lambda ),C\overline{u}(k,\lambda )^T=v(k,\lambda ).$$
(2.56)
The rule of reversing the particle order at a vertex is given as follows;

$`BA`$ $`=`$ $`(AB)^T,`$
$`BA`$ $`=`$ $`(AB)^T,`$
$`BA`$ $`=`$ $`(AB)^T.`$
Applying the above rules to the lagrangian of fermion number conserving Dirac fermions interactions (including chargino),

$$_{DD}=\overline{\mathrm{\Psi }}(D_\alpha )\mathrm{\Gamma }_{D_\alpha ,D_\beta }\mathrm{\Psi }(D_\beta )+h.c.,$$
(2.57)
one obtains the G-line rule for the Dirac particles (neglecting the overall factor $`i`$)

$`D_\alpha D_\beta `$ $`\mathrm{\Gamma }_{D_\alpha ,D_\beta },`$ (2.58)
$`D_\beta D_\alpha `$ $`\mathrm{\Gamma }_{D_\alpha ,D_\beta }^T.`$ (2.59)
Concerning the vertices including Majorana particles, corresponding to the lagrangian

$`_{MD}`$ $`=`$ $`\overline{\mathrm{\Psi }}(M_\alpha )\mathrm{\Gamma }_{M_\alpha ,D_\beta }\mathrm{\Psi }(D_\beta )+\overline{\mathrm{\Psi }}(D_\beta )\mathrm{\Gamma }_{D_\beta ,M_\alpha }\mathrm{\Psi }(M_\alpha ),`$ (2.60)
$`_{MM}`$ $`=`$ $`{\displaystyle \underset{\alpha \beta }{}}\overline{\mathrm{\Psi }}(M_\alpha )\mathrm{\Gamma }_{M_\alpha ,M_\beta }\mathrm{\Psi }(M_\beta ),`$ (2.61)
where $`\mathrm{\Gamma }_{D_\alpha ,M_\beta }\gamma _0\mathrm{\Gamma }_{M_\beta ,D_\alpha }^{}\gamma _0`$ and $`D`$ stands for a Dirac particle which can be a chargino, we obtain the following G-line rule,

$`M_\alpha D_\beta `$ $`\mathrm{\Gamma }_{M_\alpha ,D_\beta },`$ (2.62)
$`D_\alpha M_\beta `$ $`\mathrm{\Gamma }_{D_\alpha ,M_\beta },`$ (2.63)
$`M_\alpha M_\beta `$ $`\mathrm{\Gamma }_{M_\alpha ,M_\beta },(\alpha \beta )`$ (2.64)
$`M_\beta D_\alpha `$ $`C\mathrm{\Gamma }_{D_\alpha ,M_\beta }^T,`$ (2.65)
$`D_\beta M_\alpha `$ $`\mathrm{\Gamma }_{D_\beta ,M_\alpha }^TC^1.`$ (2.66)
Because the F-direction of Majorana particles is fixed to be equal to their A-direction, other vertices with Majorana particle having F-directions opposite to the A-direction do not appear in the G-lines. We also need the G-line rule for the fermion number violating chargino interactions with a conventional Dirac particle, whose lagrangian is given by

$$_{CD}=\overline{\mathrm{\Psi }}^c(C)\mathrm{\Gamma }_{C,D}\mathrm{\Psi }(D)+\overline{\mathrm{\Psi }}(D)\mathrm{\Gamma }_{D,C}\mathrm{\Psi }^c(C),$$
(2.67)
where $`\mathrm{\Gamma }_{D,C}=\gamma _0\mathrm{\Gamma }_{C,D}^{}\gamma _0`$. The rules are

$`CD`$ $`C^1\mathrm{\Gamma }_{C,D},`$ (2.68)
$`DC`$ $`\mathrm{\Gamma }_{D,C}C.`$ (2.69)
Table 2 and eq.(2.57)-(2.69) are what we need to compute the G-line amplitudes.

As Denner et al.  have done, it is possible to formulate the fermion number clashing G-line amplitude without using the charge conjugation matrix $`C`$ if one observes that the propagator with its F-direction opposite to the A-direction is given by

$``$ $`=`$ $`S^T(k)=C^1S(k)C,`$ (2.70)
$`k`$
and that the right-hand side of (2.65) and (2.66) is expressed in terms of the charge conjugated vertex function

$$\mathrm{\Gamma }_{A,B}^cC\mathrm{\Gamma }_{B,A}C^1,$$
(2.71)
as

$`M_\beta D_\alpha `$ $`\mathrm{\Gamma }_{M_\beta ,D_\alpha }^cC,`$ (2.72)
$`D_\beta M_\alpha `$ $`C^1\mathrm{\Gamma }_{D_\beta ,M_\alpha }^c,`$ (2.73)
In ref. the equivalence of these two approaches has been demonstrated.

## 3 Comments

Several comments are in order.

We have checked and debugged the possible coding errors in GRACE/SUSY by several methods. First, the interaction vertices coded in the system, which are based on the expression presented in , are compared with the result computed by the algebraic manipulation program REDUCE. As another check, we have used the gauge invariance of the processes which include photon, $`W^\pm `$, $`Z`$ and gluon. In this check one has to make the widths of the heavy gauge bosons vanish. By changing the value of gauge parameters, we have confirmed at one phase space point that the amplitudes of all the possible 582,102 processes with up to 6 external legs remain unchanged within the error of 1 in $`10^{15}`$ or less in double-precision, and 1 in $`10^{30}`$ or less in quadruple-precision calculation. Also further checks have been made by Boudjema and Bélanger of LAPTH , who compared the GRACE/SUSY results of several SUSY processes with those of CompHEP , the automatic computation program for the high energy processes developed by a group of Moscow State University.

Using GRACE/SUSY we have computed a number of SUSY processes and some of the results have been published. As this system allows us to deal with any SUSY particle production process at various types of high energy colliders, there are analyses for $`e^+e^{}`$ , $`ep`$ , $`e\gamma `$, $`\gamma \gamma `$  and $`pp`$ colliders . Note, moreover, that some works have been done with a version of GRACE/SUSY extended to the $`R`$-parity violating SUSY models .

GRACE/SUSY can treat as many final-state particles as possible, provided that one prepares by oneself the kinematics routine necessary for the phase space integration and the computer has an enough capacity to store the gigantic output files. The system is, however, equipped with its own kinematics library which has all the kinematics routines necessary for the analysis of up to 4-body final-state (see Appendix C).

Users can modify the program for their own interest but at their own risk. In particular, if one would like to use the other set of input parameters for SUSY lagrangian, for example, the trilinear coupling constants $`A_f`$ in place of $`\theta _f`$, one should prepare a program in which the parameter values are converted from the user’s set to that of GRACE/SUSY given in Table 1, and call the program before the amplitude calculation starts. Note that in the system it is $`m_fA_f`$ but not $`A_f`$ itself that has the variable name (see Appendix B). Note also that in the sfermion sector, the unmixed physical state is obtained by setting either $`\theta _f=0`$ or $`\theta _f=\pi /2`$. In the former case, however, the identification corresponds to $`\stackrel{~}{f}_1=\stackrel{~}{f}_L`$ and $`\stackrel{~}{f}_2=\stackrel{~}{f}_R`$, while in the latter $`\stackrel{~}{f}_1=\stackrel{~}{f}_R`$ and $`\stackrel{~}{f}_2=\stackrel{~}{f}_L`$. Since $`m_{\stackrel{~}{f}_1}m_{\stackrel{~}{f}_2}`$ is assumed in the system, we have to set $`\theta _f=0`$ $`(\theta _f=\pi /2)`$ when one considers the case in which there is no mixing in the sfermion sector and the left-handed sfermion is lighter (heavier) than the right-handed one.

The program GRACE/SUSY is designed strictly within the framework of the tree level. To be more realistic, however, one can arrange the value of the input parameters so that they are subject to the renormalization group equations (RGE). A program MUSE developed by Lafage  can compute the SUSY particle masses as solutions of the coupled RGE. One can combine such a program with GRACE/SUSY to invoke improved values for the input parameters. See also ref..

It is also possible to extend GRACE/SUSY so as to include new particles or new interactions such as the $`R`$-parity violating interaction (see, for example, ).

The algorithm of G-line works without any change for loop diagrams. Hence the extension of GRACE/SUSY for one-loop process is possible once one had appropriate loop evaluation routines.

In its current version of GRACE/SUSY the quark-mixing matrix is taken to be the unit matrix. Neutrinos are massless and left-handed. One may introduce the quark mixing by modifying the model files in an appropriate way.

## 4 How to use GRACE/SUSY

In this section we show how to use the program GRACE/SUSY by adopting a SUSY process $`e^+e^{}\gamma \stackrel{~}{\chi }_1^+\stackrel{~}{\chi }_1^{}`$ as an example. This process has 28 Feynman graphs at tree level, among which 4 graphs contain even number of fermion number violating vertices.

### 4.1  Definition of physical process

First the physical process to be computed must be defined, which is given in the file ”in.prc”. The content of the file in.prc for the process $`e^+e^{}\gamma \stackrel{~}{\chi }_1^+\stackrel{~}{\chi }_1^{}`$ is as follows;

> ```
> Model="mssm.mdl";
> Process;
> ELWK=3;
> QCD=0;
>   Initial={electron,positron};
>   Final  ={photon,chargino1,anti-chargino1};
> Kinem="2302";
> Pend;
>
> ```

#### 4.1.1 Model file

In the first line the model file is defined. Two model files are prepared, one is mssm.mdl for the minimal SUSY standard model and the other is ms.mdl for the standard model. The latter will be useful when one wants to separate the SUSY contributions from those of non-SUSY in such processes where only the SM particles appear in the external legs.

#### 4.1.2 Order of perturbation

The orders of perturbation are to be defined; ELWK for the electroweak interaction and QCD for the QCD.

#### 4.1.3 Initial and final states

In the lines for the initial and final states, the particle names of the process are to be given for each state according to the list given in Appendix A. One should be careful about the ordering of the particles (in particular, of massless gauge bosons) in each state. The program reads this input list and numbers all the particles from 1 to $`n`$ for $`n`$ external legs according to the order in which they appear in these two lines. These assigned numbers are in turn used to specify the integral variables in the phase space integration by the kinematic library. For example, for the 3-body final state including a single photon, the photon must appear first in the final state specification as our example shows. Otherwise the best choice of the integral variables will not be taken. See the next subsection and Appendix C for detail. In our example of $`e^+e^{}\gamma \stackrel{~}{\chi }_1^+\stackrel{~}{\chi }_1^{}`$ the input file in.prc defines the particle numbers 1$``$5 as $`e^{}`$ ($`p_1`$), $`e^+`$ ($`p_2`$), $`\gamma `$ ($`p_3`$), $`\chi _1^+`$ ($`p_4`$) and $`\chi _1^{}`$ ($`p_5`$).

#### 4.1.4 Kinematics code number

GRACE/SUSY has a built-in kinematics library which contains a set of kinematics routines for the phase space integration. In order to specify the kinematics routine to be used, one has to choose an appropriate code number, 2302 in our example. A list of code numbers is given in Appendix C.

### 4.2  Graph generation

The graph generation procedure starts by typing the command grc:

> ```
> grace% grc
>
> ```

All Feynman graphs for the process are generated and their information is saved in the file out.grf.

### 4.3  Drawing graphs

A Feynman graph drawer gracefig is initiated by the command:

> ```
> grace% gracefig
>
> ```

The procedure gracefig reads the graph information from the file out.grf and analyzes the structure of each graph, and then creates two windows. The output of our example process is shown in Fig.1.

#### 4.3.1  Drawing window of gracefig

On the drawing window the first $`m\times n`$ Feynman graphs are drawn, in the unitary gauge by default, where the number of rows $`m`$ and the number of columns $`n`$ are determined automatically by the size of screen and the total number of graphs. When the size of the window is changed, the number of columns $`n`$ is conserved while the number of rows $`m`$ is adjusted so as to fit the height of the window.

#### 4.3.2  Control window of gracefig

In the control window several buttons are set with which one can select through the pull-down menus several options which are listed in Table 3.

Particularly one can draw all the Feynman graphs in the covariant gauge by choosing an appropriate button on the menu.

#### 4.3.3  To draw selected graphs

There are two modes, “Drawing Mode” and “Select Mode”. When gracefig is invoked, all Feynman graphs in the unitary gauge are displayed in the drawing mode. If one wants to see or print only some specific graphs, one can select them in the select mode by clicking their pictures in the drawing window. The background color of the selected graphs are darkened in this mode. After this selection only the selected graphs can be displayed in the drawing mode by choosing “Selected graphs” in the graph menu.

#### 4.3.4  EPS output and insertion of a picture in a text

When “EPS output” button is clicked, the EPS file for the current page in the drawing window is created. The name of EPS file is “grcfignn.eps”, where nn is the sequential number starting from 0, and one EPS file corresponds to one page of the graphs.

### 4.4  Source code generation

After the graph generation, by typing the command

> ```
> grace% grcfort
>
> ```

FORTRAN code and Makefile are generated, which are necessary to perform gauge independence check, the numerical integration and the event generation.

#### 4.4.1 Generated source code

Three kinds of program components are created by the system. They are used in the amplitude calculation, in the phase space integration by BASES and in the event generation by SPRING, respectively.

The interrelation among the subprograms generated by GRACE/SUSY is depicted in Fig.2, where those subprograms in the solid box are automatically generated by GRACE/SUSY, while those in the dashed box are already contained in other program packages BASES/SPRING, kinematics library, interface program library to CHANEL and program package CHANEL. This figure omits a few minor components. The components vertically located in the middle of the figure are commonly used in BASES and SPRING. The program specifications of BASES/SPRING, the interface and CHANEL are described in references . A brief description of generated subprograms is given below:

* Main programs.

| mainbs | ( main ) | is the main program for the integration. |
| --- | --- | --- |
| mainsp | ( main ) | is the main program for the event generation. |
| gauge | ( main ) | is the main program for the check of gauge independence. |
* A set of program components for the integration by BASES and for the event generation by SPRING.

| userin | ( subroutine ) | controls initialization. |
| --- | --- | --- |
| setmdl | ( subroutine ) | defines MSSM parameters, masses and decay widths of particles. |
| sethptal | ( subroutine ) | defines the parameters of the extended Higgs potential. |
| setcpl | ( subroutine ) | defines coupling constants. |
| gfinit | ( subroutine ) | initializes particle polarization and selects Feynman graphs (amplitudes). |
| kinit | ( subroutine ) | initializes integration parameters, histograms, kinematics and user’s parameters. |

| func | ( function ) | calculates the numerical values of differential cross section. |
| --- | --- | --- |
| kinem | ( subroutine ) | derives particle four momenta from the integration variables. |
| kfill | ( subroutine ) | fills values in histograms and scatter-plots. |

| amptbl | ( subroutine ) | calls amnnnn to calculate amplitudes. |
| --- | --- | --- |
| ampsum | ( subroutine ) | sums the square of the matrix elements |
|  |  | over the helicity states. |
| amps | ( subroutine ) | calculates amplitudes. |
| ampord | ( subroutine ) | arranges amplitudes. |

| incl1.h | ( include file ) | defines the common variables for masses, amplitude tables etc. |
| --- | --- | --- |
| incl2.h | ( include file ) | defines the work space for amptbl. |
| inclk.h | ( include file ) | transfers values of masses and a few constants. |

| userout | ( subroutine ) | prints the amplitude summary table. Used only by gauge. |
| --- | --- | --- |

The subprograms userin and func are used both in the numerical integration by BASES and in the event generation by SPRING. In userin the subroutine kinit is called to initialize the kinematics routine and other parameters. The subroutines setmdl, setcpl and gfinit are also called for the initialization of physical parameters and amplitude calculation.

The function subprogram func is used for calculating the numerical value of differential cross section, where the subroutine kinem is called for calculating four-momentum vectors of external particles, and amptbl and ampsum are called for the amplitude calculation. The detailed specifications of these program components are described in refs. .

#### 4.4.2  Setting physical parameters for numerical calculation

Prior to running the program, one has to set the values of the input physical parameters by editing appropriate FORTRAN source codes. The parameters and corresponding subroutines to be edited are categorized as follows.

kinit, kfill

setmdl, sethptal

gfinit

* Initialization routine kinit

The values of the following parameters can be changed by editing kinit.

+ Modification of physical parameters

    - Center of mass energy in GeV. $`\sqrt{s}=1000`$GeV is a default setting, which corresponds to a description w = 1000.d0 in kinit.
        - Various kinematical cuts etc.
        - Option flags in built-in kinematics routines ().
    + Integration parameters which control the action of BASES can be altered in kinit. The following parameters can be changed from their default values according to the required integration accuracy:

    | ncall | The number of sampling points per iteration |
    | --- | --- |
    | acc1 | The expected accuracy for the grid optimization step |
    | itmx1 | The maximum iteration number for the grid optimization |
    | acc2 | The expected accuracy for the integration step |
    | itmx2 | The maximum iteration number for the integration step |

    However, the following parameters are usually not to be changed:

    | ndim | The dimensions of integral |
    | --- | --- |
    | nwild | The number of wild variables |
    | xl(i),xu(i) | The lower and upper bounds of integration variable x(i) |
    | ig(i) | The grid optimization flag for i-th variable |
    + Initialization of histograms and scatter-plots
    If the user wants histograms for some physical variables, the initialization routines xhinit and dhinit should be called for each histogram and scatter-plot, respectively.

    The system automatically generates source code to call these routines so as to construct histograms with respect to each integration variable, x(i), i=1 $``$ ndim and the energies and angles of final particles. If the user does not need these default histograms, one may delete the corresponding FORTRAN statements in kinit. When the user needs more histograms for some distributions, one should add lines to call xhinit or dhinit. These routines are to be written in the following format

    ```
     call xhinit( id#,
    .             lower_limit, upper_limit, # of bins,
    .             ’Title of this histogram’ )

    ```

    and

    ```
     call dhinit( id#,
    .             x-lower_limit, x-upper_limit, # of x bins,
    .             y-lower_limit, y-upper_limit, # of y bins,
    .             ’Title of this scatter plot’ )

    ```

    where $`1\mathtt{id\#}99`$.
* Filling histograms and scatter-plots in kfill

In the generated kfill by GRACE/SUSY, histograms for all integration variables and scatter-plots for all combinations of them are filled by calling xhfill and dhfill, respectively. If one changes the initialization of histograms and scatter-plots in kinit, their filling parts should be changed accordingly. These routines are to be called as

```
         call xhfill( id#, v, func )
         call dhfill( id#, vx, vy, func )

```
* Subprogram setmdl

The MSSM parameters, masses and widths of particles are defined in setmdl subroutine. The variable names of all the parameters are listed in Appendix B. As an example we show how to set the MSSM parameters in the case ($`\mu `$, $`M_2`$, $`\mathrm{tan}\beta `$, $`m_{\stackrel{~}{\nu _e}}`$) = ($`300`$GeV, $`200`$GeV, $`10`$, $`150`$GeV).

```
          tanbe   =     10.d0
          xmu     =   -300.d0
          xm2     =    200.d0
          amsn(1) =    150.d0

```
* Subprogram sethptal

In this subprogram one has to specify the Higgs potential to be used in the calculation by assigning a value 1 or 2 to the variable ihiggs. See section 2.4 for detail.

+ ihiggs = 1: mode 1; MSSM Higgs potential
    + ihiggs = 2: mode 2; general extended Higgs potential (default)

The default value is ihiggs=2, in which case one has to give values of the additional seven input parameters, for example,

```
          amh1  = 115.d0
          amh2  = 300.d0
          amhc  = 350.d0
          al    = -0.37d0
          c4    = 0.d0
          c9    = 0.d0
          c10   = 0.d0

```
* Subprogram gfinit

+ Selection of graphs.

    The contribution from a graph is calculated when its flag is set to be jselg = 1  in gfinit. For example, if one wants to know the contributions only from the sneutrino exchange (graphs 19, 20, 21 and 28 in Fig.1), one should write

    ```
          do 10 n1 = 1, ngraph
            jselg(n1) = 0
       10 continue
            jselg(19) = 1
            jselg(20) = 1
            jselg(21) = 1
            jselg(28) = 1

    ```

    The selection of graphs should respect gauge invariance, otherwise the cross section is physically not sensible. It should be noted that this specification of graphs is independent from the graph selection in the Feynman graph drawer gracefig.
    + Set of polarization.

    As the helicity amplitude technique is used in GRACE/SUSY, the beam polarization can be specified. The polarization sum is performed by setting the flags jhs(i) and jhe(i) for particle $`i`$ having spin $`s_i`$, where $`0\mathrm{𝚓𝚑𝚜}(𝚒)\mathrm{𝚓𝚑𝚎}(𝚒)2s_i`$. In default the unpolarized cross section is assumed, which corresponds to setting the flags for the initial $`e^{}`$ (particle 1) and $`e^+`$ (particle 2) as

    ```
    *     1: initial electron mass=amlp(1)
             ...........
          jhs(1) = 0
          jhe(1) = 1
             ............
    *     2: initial positron mass=amlp(1)
             ...........
          jhs(2) = 0
          jhe(2) = 1
             ...........

    ```

    When we consider the initial $`e^{}`$ beam polarization, for example, they should be changed to

    ```
          jhs(1) = 0
          jhe(1) = 0

    ```

    for helicity=-1 and

    ```
          jhs(1) = 1
          jhe(1) = 1

    ```

    for helicity=+1.

Similarly, one can choose polarization vectors of external gauge bosons by putting jhs(i), jhe(i) = 0, 1(transverse), 2(longitudinal: only for massive gauge bosons). In the default the linearly polarized states are realized as combinations of the transverse polarizations. When one wants to take the circular polarization for a photon, an option switch jcpol should be changed from 0 to 1.
* Gauge parameters.

Without a specification the unitary gauge is used in the integration over the phase space. However, one can select covariant gauge ($`R_\xi `$-gauge) with the gauge parameters having arbitrary values.

The selection of gauge is controlled by integer variables igauwb, igauzb, igauab and igaugl which correspond to gauge bosons $`W^\pm ,Z,\gamma `$ and gluon, respectively. When these parameters are set to be 0, unitary gauge is selected. Otherwise they should be igauwb = 1, igauzb = 2, igauab = 3 and igaugl = 4. The value of each gauge parameter should be given by agauge(igauxx), where xx is one of wb, zb, ab or gl.

In the file gfinit.f one will find the following lines which select the unitary gauge.

```
* Gauge parametes (default is unitary gauge)
      igauab = 0
      igauwb = 0
      igauzb = 0
      igaugl = 0

      agauge(0) = 1.0d20

```

In a general covariant gauge these lines should be changed as follows;

```
* Gauge parametes (covariant gauge)
      igauwb = 1
      igauzb = 2
      igauab = 3
      igaugl = 4

      agauge(igauwb) = ...
      agauge(igauzb) = ...
      agauge(igauab) = ...
      agauge(igaugl) = ...

```

### 4.5  Compilation and gauge independence check

As described in a previous subsection, Makefile is also generated by the command grcfort. The command make creates the executable files gauge, integ and spring.

Before calculating the cross section a preliminary test of the program source should be done. In GRACE/SUSY system gauge independence test is prepared for this purpose. The numerical values of the differential cross section at one sample point in the phase space are calculated both in the covariant and unitary gauge, and the consistency between these two values is examined. By executing gauge, the contributions to the cross section at one phase space point in unitary gauge (ans1) and in covariant gauge (ans2), and their difference (ans1/ans2-1) are shown as the following example demonstrates.

> ```
>     . . . . .
>    ans1   =  0.1255918104665936E+01
>     . . . . .
>           ...............
>    ans2   =  0.1255918104665936E+01
>     . . . . .
>           ...............
>    ans1/ans2 - 1 =  2.22044605E-16
>     . . .. . .
>
> ```

The relative difference should be around $`10^{16}`$ in the double-precision computation and around $`10^{30}`$ in the quadruple-precision.

It should be noted that even though this gauge independence check shows consistency between the two gauges, it does not necessarily guarantee a complete gauge independence, since it tests only at a specific point in the phase space. It is recommended to repeat a similar check at several points in phase space by changing the variable x(i) in gauge.f from their default value of 0.45 to other value lying in the range 0 \< x(i) \< 1.

### 4.6 Integration

After a successful test of gauge independence, one can proceed to the numerical integration by BASES, whose details can be found in ref.. By the make command one has already created the executable files integ for the integration and spring for event generation.

For integration the command integ is used.

> ```
> grace% integ
>
> ```

The integration parameters and the result of the integration at each iteration are shown on the screen, and at the same time they are automatically stored in the file bases.result together with the histogram outputs. If one uses a built-in kinematics routine, the cross section is given in unit of pico-barn(pb).

The numerical integration is performed by two successive steps, the first is the grid optimization step and the second the integration step. In the former the gridding of the phase space in hypercubes is optimized for the integrand, while in the latter an accurate estimate of the integral is evaluated with frozen grids determined in the first step.

Before termination of the integration procedure, BASES creates a binary file bases.data and writes the probability information on it, which is used for the event generation.

The system will issue WARNING messages when the convergence is not well established. It is always recommended to look at the integration result carefully, particularly over the convergence behavior both for the grid optimization and integration step. When the convergence accuracy fluctuates largely from one iteration to another, and/or, when in some cases it suddenly jumps up to a value far from the preceding iterations, the resultant estimate of integral is not reliable. There are two possible origins of such a behavior; one is due to too small number of sampling points and the other, more serious, is due to an unsuitable choice of the integration variables for the integrand. In the latter case one should use a better set of the integral variables.

### 4.7 Event generation

The event generation starts with the command

> ```
> grace% spring
>
> ```

Then SPRING reads the probability information from the binary file bases.data and asks the number of events with the following prompt:

> ```
> Number of events ?
>
> ```

The user must type the number of events to be generated.

If the sample point does not meet the generation criteria, such a point is discarded as a ”failed generation”. The event generation will run until a given number of events are generated or the number of failed generations exceeds its given maximum mxtry, whose default value is 50.

In order to estimate the computing time for the event generation, one can refer to the expected generation time given in the information of BASES output, where the generation efficiency of 70% is assumed.

#### 4.7.1  Analysis of events

In the event generation four-momentum vectors of the final particles are automatically stored in vec(i,n) in mainsp, where n stands for the particle number and $`𝚒=\mathrm{𝟷}\mathrm{𝟺}`$ correspond to the three-momentum (1, 2, 3) and the energy (4) of the particle.

To create histograms of some variables of the generated events, one has to initialize them before the generation loop and fill them after each call of spring. This kind of histograms is called as an “additional histogram” in contrast to the “original histogram” which is filled both in the integration and event generation steps. In the original histograms one can see how the generated events reproduce those distributions produced by the integration, while in the additional histogram only the frequency distribution can be seen as in the usual histogram. An example for saving the four-momentum vectors and creating additional histograms is given below.

> ```
>
> *     Open the output file for the events
>       . . . . . .
>
> *     Initialize the additional histograms by
>         call xhinit( id, ... ).
>
>       do 100 nevnt = 1, mxevnt
>
>          call spring( func, mxtry )
>
> *    ------------------------------------------------
> *        analyze the event and
> *        fill the additional histograms by
>           call xhfill( id, .... ).
> *    ------------------------------------------------
>
>          do 90 k = 1 , nextrn
>             write(6,*) (vec(j,k),j=1,4)
>    90    continue
>   100 continue
>
> ```

#### 4.7.2  Output from SPRING

The output from SPRING is written on the file named spring.result, which consists of the general information, original and additional histograms, scatter-plots, and the number-of-trials distribution. In the last distribution one can see the generation efficiency. It is recommended to check if this distribution has a sharp peak in the first bin. If this is not the case, the integration should be carried out again with more sampling points or with more iterations.

We would like to thank F. Boudjema and G. Bélanger for checking parts of the program and reading the manuscript. We would like to thank also D. Perret-Gallix for his interest in this work and for encouragement. This work was partly supported by Japan Society for Promotion of Science under the Grant-in-Aid for Scientific Research B ( No.11440083 ).

## Appendix A List of particles in MSSM

The particles appearing in MSSM and their variable names are listed in two tables below. These variable names are to be used in specifying the process in in.prc.

| particle | variable name |
| --- | --- |
| photon | photon |
| $`W^+(W^{})`$ | W-plus (W-minus) |
| $`Z`$ | Z |
| gluon | gluon |
| $`\nu _e(\overline{\nu _e})`$ | nu-e (nu-e-bar) |
| $`e^{}(e^+)`$ | electron (positron) |
| $`\nu _\mu (\overline{\nu _\mu })`$ | nu-mu (nu-mu-bar) |
| $`\mu ^{}(\mu ^+)`$ | muon (anti-muon) |
| $`\nu _\tau (\overline{\nu _\tau })`$ | nu-tau (nu-tau-bar) |
| $`\tau ^{}(\tau ^+)`$ | tau (anti-tau) |
| $`u(\overline{u})`$ | u (u-bar) |
| $`d(\overline{d})`$ | d (d-bar) |
| $`c(\overline{c})`$ | c (c-bar) |
| $`s(\overline{s})`$ | s (s-bar) |
| $`t(\overline{t})`$ | t (t-bar) |
| $`b(\overline{b})`$ | b (b-bar) |
| $`h^0`$ | Higgs1 |
| $`H^0`$ | Higgs2 |
| $`A^0`$ | Higgs3 |
| $`H^+(H^{})`$ | Higgs-plus (Higgs-minus) |

| particle | variable name |
| --- | --- |
| $`\stackrel{~}{\chi _1^+}(\stackrel{~}{\chi _1^{}})`$ | chargino1 (anti-chargino1) |
| $`\stackrel{~}{\chi _2^+}(\stackrel{~}{\chi _2^{}})`$ | chargino2 (anti-chargino2) |
| $`\stackrel{~}{\chi }_1^0`$ | neutralino1 |
| $`\stackrel{~}{\chi }_2^0`$ | neutralino2 |
| $`\stackrel{~}{\chi }_3^0`$ | neutralino3 |
| $`\stackrel{~}{\chi }_4^0`$ | neutralino4 |
| $`\stackrel{~}{\nu _e}(\overline{\stackrel{~}{\nu _e}})`$ | snu-e (anti-snu-e) |
| $`\stackrel{~}{\nu _\mu }(\overline{\stackrel{~}{\nu _\mu }})`$ | snu-mu (anti-snu-mu) |
| $`\stackrel{~}{\nu _\tau }(\overline{\stackrel{~}{\nu _\tau }})`$ | snu-tau (anti-snu-tau) |
| $`\stackrel{~}{e_1^{}}(\stackrel{~}{e_1^+})`$ | selectron1 (anti-selectron1) |
| $`\stackrel{~}{e_2^{}}(\stackrel{~}{e_2^+})`$ | selectron2 (anti-selectron2) |
| $`\stackrel{~}{\mu _1^{}}(\stackrel{~}{\mu _1^+})`$ | smuon1 (anti-smuon1) |
| $`\stackrel{~}{\mu _2^{}}(\stackrel{~}{\mu _2^+})`$ | smuon2 (anti-smuon2) |
| $`\stackrel{~}{\tau _1^{}}(\stackrel{~}{\tau _1^+})`$ | stau1 (anti-stau1) |
| $`\stackrel{~}{\tau _2^{}}(\stackrel{~}{\tau _2^+})`$ | stau2 (anti-stau2) |
| $`\stackrel{~}{u_1}(\overline{\stackrel{~}{u_1}})`$ | su1 (anti-su1) |
| $`\stackrel{~}{u_2}(\overline{\stackrel{~}{u_2}})`$ | su2 (anti-su2) |
| $`\stackrel{~}{d_1}(\overline{\stackrel{~}{d_1}})`$ | sd1 (anti-sd1) |
| $`\stackrel{~}{d_2}(\overline{\stackrel{~}{d_2}})`$ | sd2 (anti-sd2) |
| $`\stackrel{~}{c_1}(\overline{\stackrel{~}{c_1}})`$ | sc1 (anti-sc1) |
| $`\stackrel{~}{c_2}(\overline{\stackrel{~}{c_2}})`$ | sc2 (anti-sc2) |
| $`\stackrel{~}{s_1}(\overline{\stackrel{~}{s_1}})`$ | ss1 (anti-ss1) |
| $`\stackrel{~}{s_2}(\overline{\stackrel{~}{s_2}})`$ | ss2 (anti-ss2) |
| $`\stackrel{~}{t_1}(\overline{\stackrel{~}{t_1}})`$ | st1 (anti-st1) |
| $`\stackrel{~}{t_2}(\overline{\stackrel{~}{t_2}})`$ | st2 (anti-st2) |
| $`\stackrel{~}{b_1}(\overline{\stackrel{~}{b_1}})`$ | sb1 (anti-sb1) |
| $`\stackrel{~}{b_2}(\overline{\stackrel{~}{b_2}})`$ | sb2 (anti-sb2) |

## Appendix B Name List of GRACE/SUSY Parameters

The variable names used in GRACE/SUSY are summarized. All mass variables. widths and the dimensional coupling constant $`\mu `$ must be given in unit of GeV. Those variables whose magnitude can be specified by the user in setmdl.f and sethptal.f are indicated by an asterisk ”\*” at the end of the row. The value of other variables should not be changed by the user. Such a modification may cause internal inconsistencies.

gauge bosons

| category | notation | variable name in GRACE/SUSY |  |
| --- | --- | --- | --- |
| gauge boson masses | $`M_Z`$, $`M_W`$ | amz, amw | * |
| gauge boson widths | $`\mathrm{\Gamma }_Z`$, $`\mathrm{\Gamma }_W`$ | agz, agw | * |
| coupling | $`\alpha _e`$, $`\alpha _s`$ | alpha, alphas | * |
|  | $`e`$, $`g_s`$ | ge, gs |  |
|  | $`g`$, $`g^{}`$ | gg, gz |  |
| mixing angle | $`\mathrm{cos}\theta _W`$, $`\mathrm{sin}\theta _W`$, $`\mathrm{tan}\theta _W`$ | cw, sw, tw |  |

Ordinary fermions

| fermion mass | $`m_u`$, $`m_c`$, $`m_t`$ | amuq(1), amuq(2), amuq(3) | * |
| --- | --- | --- | --- |
|  | $`m_d`$, $`m_s`$, $`m_b`$ | amdq(1), amdq(2), amdq(3) | * |
|  | $`m_e`$, $`m_\mu `$, $`m_\tau `$ | amlp(1), amlp(2), amlp(3) | * |
| fermion widths | $`\mathrm{\Gamma }_u`$, $`\mathrm{\Gamma }_c`$, $`\mathrm{\Gamma }_t`$ | aguq(1), aguq(2), aguq(3) | * |
|  | $`\mathrm{\Gamma }_d`$, $`\mathrm{\Gamma }_s`$, $`\mathrm{\Gamma }_b`$ | agdq(1), agdq(2), agdq(3) | * |
|  | $`\mathrm{\Gamma }_e`$, $`\mathrm{\Gamma }_\mu `$, $`\mathrm{\Gamma }_\tau `$ | aglp(1) , aglp(2), aglp(3) | * |

Of course one should set aguq(1) = agdq(1) = aglp(1) = 0.

Higgs bosons

| SUSY parameters | $`\mu `$ | xmu | * |
| --- | --- | --- | --- |
| Higgs masses | $`M_{A^0}`$ | amh3 | * |
|  | $`M_h`$, $`M_H`$, $`M_{H^\pm }`$ | amh1, amh2, amhc | * |
| Higgs widths | $`\mathrm{\Gamma }_h`$, $`\mathrm{\Gamma }_H`$, $`\mathrm{\Gamma }_A`$, $`\mathrm{\Gamma }_{H^\pm }`$ | agh1, agh2, agh3, aghc | * |
| mixing angle | $`\alpha `$, $`\mathrm{cos}\alpha `$, $`\mathrm{sin}\alpha `$ | al, ca, sa | * |
|  | $`\mathrm{cos}\beta `$, $`\mathrm{sin}\beta `$,$`\mathrm{tan}\beta `$, $`\mathrm{cot}\beta `$ | cb, sb, tanbe, cotbe | * |

Chargino

| SUSY parameters | $`\mu `$, $`M_2`$ | xmu, xm2 | * |
| --- | --- | --- | --- |
| chargino masses | $`m_{\stackrel{~}{\chi }_1^\pm }`$, $`m_{\stackrel{~}{\chi }_2^\pm }`$, | amsw(1), amsw(2) |  |
| chargino widths | $`\mathrm{\Gamma }_{\stackrel{~}{\chi }_1^\pm }`$, $`\mathrm{\Gamma }_{\stackrel{~}{\chi }_2^\pm }`$, | agsw(1), agsw(2) | * |
| chargino mixing angles | $`\mathrm{sin}\varphi _L`$, $`\mathrm{cos}\varphi _L`$, $`ϵ_L`$ | sphl, cphl, epsl |  |
|  | $`\mathrm{sin}\varphi _R`$, $`\mathrm{cos}\varphi _R`$ | sphr, cphr |  |

Neutralino

| SUSY parameters | $`\mu `$, $`M_1`$, $`M_2`$ | xmu, xm1, xm2 | * |
| --- | --- | --- | --- |
| neutralino masses | $`m_{\stackrel{~}{\chi }_i^0}`$ | amnl(i) |  |
| neutralino widths | $`\mathrm{\Gamma }_{\stackrel{~}{\chi }_i^0}`$ | agnl(i) | * |
| neutralino mixing matrix | $`(𝒪_N)_{ij}`$, $`\eta _i`$ | ogmn(i,j), etan(i)(complex) |  |

Gluino

| gluino mass | $`m_{\stackrel{~}{g}}=M_3`$ | amgl $`=`$ xm3 | * |
| --- | --- | --- | --- |
| gluino width | $`\mathrm{\Gamma }_{\stackrel{~}{g}}`$ | aggl | * |

Sfermion

| sfermion masses | $`m_{\stackrel{~}{u}_i}`$, $`m_{\stackrel{~}{c}_i}`$, $`m_{\stackrel{~}{t}_i}`$ | amsu(1,i), amsu(2,i), amsu(3,i) | * |
| --- | --- | --- | --- |
|  | $`m_{\stackrel{~}{d}_i}`$, $`m_{\stackrel{~}{s}_i}`$, $`m_{\stackrel{~}{b}_i}`$ | amsd(1,i), amsd(2,i), amsd(3,i) | * |
|  | $`m_{\stackrel{~}{\nu }_e}`$, $`m_{\stackrel{~}{\nu }_\mu }`$, $`m_{\stackrel{~}{\nu }_\tau }`$ | amsn(1), amsn(2), amsn(3) | * |
|  | $`m_{\stackrel{~}{e}_i}`$, $`m_{\stackrel{~}{\mu }_i}`$, $`m_{\stackrel{~}{\tau }_i}`$ | amsl(1,i), amsl(2,i), amsl(3,i) | * |
| sfermion widths | $`\mathrm{\Gamma }_{\stackrel{~}{u}_i}`$, $`\mathrm{\Gamma }_{\stackrel{~}{c}_i}`$, $`\mathrm{\Gamma }_{\stackrel{~}{t}_i}`$ | agsu(1,i), agsu(2,i), agsu(3,i) | * |
|  | $`\mathrm{\Gamma }_{\stackrel{~}{d}_i}`$, $`\mathrm{\Gamma }_{\stackrel{~}{s}_i}`$, $`\mathrm{\Gamma }_{\stackrel{~}{b}_i}`$ | agsd(1,i), agsd(2,i), agsd(3,i) | * |
|  | $`\mathrm{\Gamma }_{\stackrel{~}{\nu }_e}`$, $`\mathrm{\Gamma }_{\stackrel{~}{\nu }_\mu }`$, $`\mathrm{\Gamma }_{\stackrel{~}{\nu }_\tau }`$ | agsn(1), agsn(2), agsn(3) | * |
|  | $`\mathrm{\Gamma }_{\stackrel{~}{e}_i}`$, $`\mathrm{\Gamma }_{\stackrel{~}{\mu }_i}`$, $`\mathrm{\Gamma }_{\stackrel{~}{\tau }_i}`$ | agsl(1,i), agsl(2,i), agsl(3,i) | * |
| SUSY-breaking | $`m_uA_u`$, $`m_cA_c`$, $`m_tA_t`$ | xauq(1), xauq(2), xauq(3) |  |
| parameters | $`m_dA_d`$, $`m_sA_s`$, $`m_bA_b`$ | xadq(1), xadq(2), xadq(3) |  |
|  | $`m_eA_e`$, $`m_\mu A_\mu `$, $`m_\tau A_\tau `$ | xalp(1), xalp(2), xalp(3) |  |
| sfermion mixing | $`\mathrm{sin}\theta _{\stackrel{~}{u}}`$, $`\mathrm{sin}\theta _{\stackrel{~}{c}}`$, angles $`\mathrm{sin}\theta _{\stackrel{~}{t}}`$ | shuq(1), shuq(2), shuq(3) | * |
|  | $`\mathrm{cos}\theta _{\stackrel{~}{u}}`$, $`\mathrm{cos}\theta _{\stackrel{~}{c}}`$, $`\mathrm{cos}\theta _{\stackrel{~}{t}}`$ | chuq(1), chuq(2), chuq(3) | * |
|  | $`\mathrm{sin}\theta _{\stackrel{~}{d}}`$, $`\mathrm{sin}\theta _{\stackrel{~}{s}}`$, $`\mathrm{sin}\theta _{\stackrel{~}{b}}`$ | shdq(1), shdq(2), shdq(3) | * |
|  | $`\mathrm{cos}\theta _{\stackrel{~}{d}}`$, $`\mathrm{cos}\theta _{\stackrel{~}{s}}`$, $`\mathrm{cos}\theta _{\stackrel{~}{b}}`$ | chdq(1), chdq(2), chdq(3) | * |
|  | $`\mathrm{sin}\theta _{\stackrel{~}{e}}`$,$`\mathrm{sin}\theta _{\stackrel{~}{\mu }}`$, $`\mathrm{sin}\theta _{\stackrel{~}{\tau }}`$ | shlp(1), shlp(2), shlp(3) | * |
|  | $`\mathrm{cos}\theta _{\stackrel{~}{e}}`$,$`\mathrm{cos}\theta _{\stackrel{~}{\mu }}`$, $`\mathrm{cos}\theta _{\stackrel{~}{\tau }}`$ | chlp(1), chlp(2), chlp(3) | * |

## Appendix C Kinematics Library

### C.1 Cross section formula

The volume element for each of final particles is given by

$$\frac{\mathrm{d}^3p}{2E(2\pi )^3},$$
(3.74)
where $`𝐩`$ is the three momentum of the particle and $`E`$ is its energy.

In a collision of two particles labeled by 1 and 2, the flux is given by

$$\text{flux}=v_{rel}2E_12E_2,v_{rel}=\left|\frac{𝐩_1}{E_1}\frac{𝐩_2}{E_2}\right|,$$
(3.75)
with $`v_{rel}`$ being the relative velocity of particles 1 and 2. The total cross section is then written as

$$\sigma =\frac{1}{\text{flux}}\underset{j}{}\frac{\mathrm{d}^3p_j}{2E_j(2\pi )^3}(2\pi )^4\delta ^4(p_1+p_2\underset{j}{}p_j)\underset{h_f}{}\underset{h_i}{}|T_{if}|^2,$$
(3.76)
where $`T_{if}`$ is the helicity amplitude for the process $`if`$. Initial helicities are averaged and the final ones are summed. This formula is used for the cross section in the GRACE/SUSY system. The cross section is given in the unit of pico-barn(pb) in the kinematics library.

For the decay of a particle with momentum $`p`$, the total decay rate is given by

$$\mathrm{\Gamma }=\frac{1}{2E}\underset{j}{}\frac{\mathrm{d}^3p_j}{2E_j(2\pi )^3}(2\pi )^4\delta ^4(p\underset{j}{}p_j)\underset{h_f}{}\underset{h_i}{}|T_{if}|^2.$$
(3.77)
The independent integration variables expressed in terms of momenta and angles must be mapped onto those of the integration package BASES which are normalized as $`0𝚡(𝚒)1`$, where i runs from 1 to ndim, the dimension of the integration,

One should remember that the statistical factor for the identical particles is not generated automatically and should be supplied by hand. This is because this factor depends on the observables one is interested in. For the total cross section one should divide the result obtained by BASES by the factorial of the number of identical particles.

### C.2 Kinematics library

GRACE/SUSY assumes the following two types of processes:

$$p_1+p_2p_3+p_4+\mathrm{}(\text{2-body scattering}),$$

$$p_1p_2+p_3+p_4+\mathrm{}(\text{decay of a particle}).$$

It should be noted that the particle assignment follows the order in which the particles are originally defined in in.prc.

The Lorentz invariant phase space element for a final $`n`$-body ($`A1+2+3+\mathrm{}+n`$) is defined by

$`d\stackrel{~}{\mathrm{\Gamma }}_n`$ $`=`$ $`(2\pi )^4\delta ^4\left({\displaystyle \underset{in}{}}p{\displaystyle \underset{out}{\overset{n}{}}}p\right){\displaystyle \underset{out}{\overset{n}{}}}{\displaystyle \frac{d^3p_j}{(2\pi )^32E_j}}{\displaystyle \frac{1}{(2\pi )^{3n4}}}d\mathrm{\Gamma }_n,`$ (3.78)
$`d\mathrm{\Gamma }_n`$ $`=`$ $`\delta ^4\left({\displaystyle \underset{in}{}}p{\displaystyle \underset{out}{\overset{n}{}}}p\right){\displaystyle \underset{out}{\overset{n}{}}}{\displaystyle \frac{d^3p_j}{2E_j}}.`$ (3.79)
The following chain rule is useful($`0<k<n1`$):

$`d\mathrm{\Gamma }_n(A1+2+3+\mathrm{})`$ (3.80)
$`=d\mathrm{\Gamma }_{k+1}(A1+\mathrm{}+k+q){\displaystyle \frac{dQ^2}{2\pi }}`$
$`\times d\mathrm{\Gamma }_{nk}(q(k+1)+(k+2)+\mathrm{}),`$
where $`q^2=Q^2`$.

When $`p_a`$ is in the center-of-mass system, i.e., $`𝐩_a=\mathrm{𝟎}`$, and if $`p_a=p_b+p_c`$ we introduce the 2-body phase space

$`d\mathrm{\Gamma }_{\mathrm{CM}}(a;bc)d\mathrm{\Gamma }_2=\delta ^4(p_ap_bp_c){\displaystyle \frac{d^3p_b}{2E_b}}{\displaystyle \frac{d^3p_c}{2E_c}},`$ (3.81)
and write it by angular variables:

$`d\mathrm{\Gamma }_{\mathrm{CM}}(a;bc)={\displaystyle \frac{\beta (a;bc)}{8}}d\mathrm{\Omega }_{\mathrm{CM}}(a;bc)={\displaystyle \frac{\beta (a;bc)}{8}}d\mathrm{cos}\theta _{b,(bc)}d\varphi _{b,(bc)}.`$ (3.82)
Here $`\beta `$ is given by

$`\beta (a;bc)`$ $`=`$ $`{\displaystyle \frac{2P}{E_a}}`$
$`=`$ $`{\displaystyle \frac{1}{E_a}}\sqrt{\left(E_a+m_b+m_c\right)\left(E_am_bm_c\right)\left(E_a+m_bm_c\right)\left(E_am_b+m_c\right)}.`$
The subscript $`(bc)`$ of the angles indicates that $`\theta _b,\varphi _b`$ are defined in the center-of-mass frame. Angles in the laboratory frame have no subscript.

Kinematics routines prepared in the library are summarized in Table C.2.

| code number | contents |
| --- | --- |
| 1201 | 1-body $``$ 2 body decay |
| 1301 | 1-body $``$ 3 body decay |
|  | Sequential decay $`12+(3+4)2+3+4`$ can be treated. |
| 2201 | 2-body $``$ 2 body in CM frame |
|  | $`t`$\- and $`u`$-channel singularities can be treated. |
| 2301 | 2-body $``$ 3 body in CM frame , |
|  | Sequential decay type $`1+23+(4+5)3+4+5`$. |
|  | Resonance on particles 4 and 5 can be treated. |
| 2302 | 2-body $``$ 3 body in CM frame , |
|  | Radiative processes $`1+23(\gamma )+4+5`$, |
|  | both initial and final radiation can be treated. |
| 2303 | 2-body $``$ 3 body in CM frame , |
|  | Double-radiative processes $`1+23(\gamma )+4(\gamma )+5`$ |
| 2304 | 2-body $``$ 3 body in CM frame , |
|  | Three photon processes $`1+23(\gamma )+4(\gamma )+5(\gamma )`$ |
| 2401 | 2-body $``$ 4 body in CM frame, a pair of sequential |
|  | decay type $`1+2(3+4)+(5+6)3+4+5+6`$ |
|  | $`t`$-channel singularity can be treated. |
| 2402 | 2-body $``$ 4 body in CM frame, |
|  | ‘fusion’ type $`1+2(3+A)+(4+B);A+B5+6`$ |

Table C.2 Kinematics library prepared in GRACE/SUSY system.


format: text

text:
# 1 Feynman diagrams for 𝑒⁺𝑒⁻→ℎ⁰𝜏̃⁺₁𝜏̃⁻₁ production.

BIHEP-TH–2002-31

Associated production of a Higgs boson with tau sleptons in the *CP* violating MSSM at future $`e^+e^{}`$ colliders
Wei Min Yang***Email address: yangwm@mail.ihep.ac.cn and Dong Sheng Du
CCAST(World Laboratory), P.O.Box 8730, Beijing 100080, China and Institute of High Energy Physics, P.O.Box 918(4), Beijing 100039, China
ABSTRACT

We investigate the production of the lightest neutral Higgs boson in association with tau sleptons in the *CP* violating minimal supersymmetric standard model (MSSM) at future high-energy $`e^+e^{}`$ linear colliders. In parameter space of the constrained MSSM, the production cross section of $`e^+e^{}h^0\stackrel{~}{\tau }_1^+\stackrel{~}{\tau }_1^{}`$ can be very substantial at high energies. This process would provide a new production mechanism of neutral Higgs bosons and open a window to probe the Higgs-stau coupling and some soft supersymmetric breaking parameters at next linear colliders (NLC).

PACS numbers: 12.60.Jv, 14.80.Cp, 14.80.Ly

I. INTRODUCTION

Supersymmetric (SUSY) theories, in particular the MSSM , are currently considered as the most theoretically well motivated extensions of the standard model (SM). It provides a natural solution to the hierarchy problem , the generation of the electroweak symmetry breaking , the grand unification of the gauge coupling constants , as well as the elucidation of the baryon asymmetry and the cold dark matter in the universe . The MSSM contains five physical Higgs bosons : two neutral *CP*-even scalars ($`h^0`$ and $`H^0`$), one *CP*-odd neutral scalar $`A^0`$, and a pair of charged scalars ($`H^+`$, $`H^{}`$). Such particles break the electroweak symmetry and are responsible for the masses of the fermions and bosons. The search for the Higgs bosons and the study of their properties are one of the major goals of present and future colliders . Contrary to hadron machines, at $`e^+e^{}`$ colliders the background is relatively clean. A future high-energy and high-luminosity $`e^+e^{}`$ linear collider operating at energies $`\sqrt{s}=5001000`$ GeV would provide an ideal environment for searching the lightest Higgs boson and performing precision measurements of it .

In the context of the MSSM, some typical production mechanisms of Higgs bosons, such as Higgsstrahlung , vector bosons fusion , Higgs pair production , and associated production of Higgs boson with top or stop , have all been extensively studied at NLC energies. In addition, it was realized that couplings of Higgs bosons to sfermions can also provide additional sources for Higgs bosons production , in particular, when the mixing between the third generation scalar taus is very large, the associated production Higgs bosons with tau sleptons might provide a important information on the soft SUSY breaking scalar potential. In this paper, we will discuss in more detail the last type of the prospecting processes in the *CP* violating MSSM. Our motivation is as the following: First, for the large Yukawa coupling of the tau slepton, the mixing between tau sleptons can be very large. As a result, the lighter stau can not only be significantly lighter than the top quark and all scalar quarks, but also at the same time their couplings to Higgs bosons can be strongly enhanced. Therefore, the reaction $`e^+e^{}h^0\stackrel{~}{\tau }_1^+\stackrel{~}{\tau }_1^{}`$ can be more phase-space favored, and also it would be accessible to produce and measure them at the NLC. Second, it has received growing attention that *CP*-violating phases associated with the third generation trilinear soft term can induce *CP* violation in the Higgs and the third generation sfermion sectors through quantum corrections . Such phases may allow baryogenesis and do not necessarily violate the stringent bound from the nonobservation of electric dipole moments (EDM) . In fact, some of these phases can be $`𝒪(1)`$, so as to provide non-SM sources of *CP* violation required for dynamical generation of the baryon asymmetry of the universe . Since these phases directly affect the couplings of Higgs bosons to third generation sfermions, moreover, these couplings play an important role in the phenomenology of the third generation sfermions and Higgs bosons at colliders, we will consider these effects of the *CP*-violating phases on the associated Higgs production with the stau\*-stau pair. Last, for most of the SUSY parameter space allowed by present data constrains, the MSSM Higgs sector is in the decoupling regime . In this scenario, only the lightest Higgs boson would be accessible at the NLC, final states with the other Higgs bosons or with $`\stackrel{~}{\tau }_1^+\stackrel{~}{\tau }_2^{}`$ pair should be phase-space suppressed. In what follows, therefore, we will focus on the case of the lightest $`h^0`$ boson of MSSM in the decoupling regime, and discuss only the production in association with lighter tau sleptons.

The remainder of this paper is organized as follows. In Section II we outline the masses and mixing of the tau sleptons, and then we present the relevant couplings and analytical expressions of the production cross section for the process $`e^+e^{}h^0\stackrel{~}{\tau }_1^+\stackrel{~}{\tau }_1^{}`$. In Sec.III, a detailed numerical analysis of the cross section is given in an appropriate MSSM scenario satisfying the experimental and theoretical constraints. Sec.IV is devoted to the conclusions.

II. MASSES, COUPLINGS AND CROSS SECTION

In this section, first, to fix notation we will simply summarize the masses and mixing of the tau slepton sector in the MSSM, and then list the relevant couplings involved in the associated production of the lightest neutral Higgs boson with tau sleptons. Finally, we will give analytical results of the cross section of $`e^+e^{}h^0\stackrel{~}{\tau }_1^+\stackrel{~}{\tau }_1^{}`$ production.

A. Tau slepton mass and mixing

In the tau slepton sector , the left- and right-handed tau sleptons can strongly mix, the mass-squared matrix is given by

$`M_{\stackrel{~}{\tau }}^2=\left(\begin{array}{cc}M_L^2& m_\tau (A_\tau +\mu \mathrm{tan}\beta )\\ m_\tau (A_\tau ^{}+\mu ^{}\mathrm{tan}\beta )& M_R^2\end{array}\right)`$ (3)
with
$`M_L^2`$ $`=m_{\stackrel{~}{\tau }_L}^2+m_\tau ^2+({\displaystyle \frac{1}{2}}M_Z^2M_W^2)\mathrm{cos}2\beta ,`$    $`M_R^2`$ $`=m_{\stackrel{~}{\tau }_R}^2+m_\tau ^2+(M_Z^2M_W^2)\mathrm{cos}2\beta ,`$
where $`m_{\stackrel{~}{\tau }_{L,R}}`$ are the left- and right-handed soft SUSY-breaking tau slepton masses, respectively. The soft breaking trilinear stau coupling $`A_\tau `$ and Higgs mass mixing parameter $`\mu `$ are complex,

$`A_\tau =|A_\tau |e^{i\phi },\mu =|\mu |e^{i\eta }.`$ (5)
The phases $`\phi `$ and $`\eta `$ are the sources of *CP* violation. In fact, the phase $`\eta `$ can be transferred into the phase $`\phi `$ in the form of $`\phi \eta `$ by the field redefinition. So we will take $`\mu `$ to be real, and vary $`\phi `$ from $`0`$ to $`2\pi `$ in the following. The tau slepton mass eigenstates can be realized by a unitary transformation $`U`$ which diagonalizes the mass-squared matrix $`M_{\stackrel{~}{\tau }}^2`$ ,

$`UM_{\stackrel{~}{\tau }}^2U^{}=\text{diag}(m_{\stackrel{~}{\tau }_1}^2,m_{\stackrel{~}{\tau }_2}^2),`$ (6)
where the diagonalization matrix can be parameterized as

$`U=\left(\begin{array}{cc}\mathrm{cos}\theta _\tau & \mathrm{sin}\theta _\tau e^{i\delta }\\ \mathrm{sin}\theta _\tau e^{i\delta }& \mathrm{cos}\theta _\tau \end{array}\right)`$ (9)
with

$`\delta =\mathrm{arg}(A_\tau +\mu \mathrm{tan}\beta ).`$ (10)
The tau slepton mixing angle and mass eigenvalues are then given as
$`\mathrm{tan}2\theta `$ $`={\displaystyle \frac{2m_\tau |A_\tau +\mu \mathrm{tan}\beta |}{M_L^2M_R^2}},`$    $`m_{\stackrel{~}{\tau }_{1,2}}^2`$ $`={\displaystyle \frac{1}{2}}\left[M_L^2+M_R^2\sqrt{(M_L^2M_R^2)^2+4m_\tau ^2|A_\tau +\mu \mathrm{tan}\beta |^2}\right]`$
with convention $`0\theta \pi /2,m_{\stackrel{~}{\tau }_1}^2m_{\stackrel{~}{\tau }_2}^2`$. For large values of $`|A_\tau |`$, $`|\mu |`$ and $`\mathrm{tan}\beta `$, the mixing in the tau slepton sector can be very strong. In the constrained MSSM or minimal supergravity (mSUGRA) scenario , the soft SUSY breaking scalar fermion masses and gaugino masses are respectively unified as $`m_0`$ and $`m_{1/2}`$ at the grand unified theory (GUT) scale $`M_{GUT}`$. The left- and right-handed tau slepton masses at the weak scale $`𝒪(M_Z)`$ are then given in terms of the universal masses $`m_0`$ and $`m_{1/2}`$ through running renormalization group equations (RGE) as
$`m_{\stackrel{~}{\tau }_L}^2`$ $`m_0^2+0.52m_{1/2}^20.27M_Z^2\mathrm{cos}2\beta ,`$    $`m_{\stackrel{~}{\tau }_R}^2`$ $`m_0^2+0.15m_{1/2}^20.23M_Z^2\mathrm{cos}2\beta .`$
Inputting values of the parameters $`|A_\tau |,|\mu |,\phi ,\mathrm{tan}\beta ,m_0,m_{1/2}`$, the mass spectra and mixing angle of the tau sleptons can then be worked out. For smaller values of $`m_0`$ and $`m_{1/2}`$ , the stau masses are lower.

B. $`𝐞^+𝐞^{}𝐡^\mathrm{𝟎}\stackrel{~}{\tau }_\mathrm{𝟏}^+\stackrel{~}{\tau }_\mathrm{𝟏}^{}`$ production

$`e^+e^{}h^0\stackrel{~}{\tau }_1^+\stackrel{~}{\tau }_1^{}`$ production proceeds via three kinds of diagrams, as shown in Fig. 1. When the $`h^0\stackrel{~}{\tau }_1^{}\stackrel{~}{\tau }_1`$ couplings are large, the most important is Higgs boson emission from the $`\stackrel{~}{\tau _1}`$ states which are produced through s-channel photon and Z-boson exchange (Diag. a and b in Fig. 1). One also has the conversion of $`\stackrel{~}{\tau _2}`$ to $`\stackrel{~}{\tau _1}`$, which involves the $`h^0\stackrel{~}{\tau }_1^{}\stackrel{~}{\tau }_2`$ vertex (Diag. c and d in Fig. 1). The last type is $`h^0Z^{}`$ production with $`Z^{}\stackrel{~}{\tau }_1^+\stackrel{~}{\tau }_1^{}`$ (Diag. e in Fig. 1). Though the two last types of Feynman diagrams turn out to give small contributions, there are some regions in the parameter space where they cannot be neglected. The relevant couplings for the process are described by the following set of interaction Lagrangian in the MSSM :
(a). The couplings of a neutral Higgs boson $`h^0`$ to a pair of tau sleptons $`\stackrel{~}{\tau }_i\stackrel{~}{\tau }_j(i,j=1,2)`$ is given by

$`_{h^0\stackrel{~}{\tau }_i^{}\stackrel{~}{\tau }_j}=g_{ij}h^0\stackrel{~}{\tau }_i^{}\stackrel{~}{\tau }_j`$ (13)
with

$`g_{ij}=`$ $`{\displaystyle \frac{gm_Z}{c_W}}\left[(s_W^2{\displaystyle \frac{1}{2}})\mathrm{sin}(\alpha +\beta )+\left({\displaystyle \frac{m_\tau }{m_Z}}\right)^2{\displaystyle \frac{\mathrm{sin}\alpha }{\mathrm{cos}\beta }}\right](U_{1i}^{})^{}U_{1j}^{}`$
$`+{\displaystyle \frac{gm_Z}{c_W}}\left[s_W^2\mathrm{sin}(\alpha +\beta )+\left({\displaystyle \frac{m_\tau }{m_Z}}\right)^2{\displaystyle \frac{\mathrm{sin}\alpha }{\mathrm{cos}\beta }}\right](U_{2i}^{})^{}U_{2j}^{}`$
$`+{\displaystyle \frac{gm_\tau }{2c_Wm_Z}}\left[A_\tau ^{}{\displaystyle \frac{\mathrm{sin}\alpha }{\mathrm{cos}\beta }}\mu ^{}{\displaystyle \frac{\mathrm{cos}\alpha }{\mathrm{cos}\beta }}\right](U_{2i}^{})^{}U_{1j}^{}`$
$`+{\displaystyle \frac{gm_\tau }{2c_Wm_Z}}\left[A_\tau {\displaystyle \frac{\mathrm{sin}\alpha }{\mathrm{cos}\beta }}\mu {\displaystyle \frac{\mathrm{cos}\alpha }{\mathrm{cos}\beta }}\right](U_{1i}^{})^{}U_{2j}^{},`$ (14)
where $`s_W=\mathrm{sin}\theta _W,c_W=\mathrm{cos}\theta _W`$, $`\alpha `$ is the *CP*-even Higgs mixing angle at the tree level. When the loop corrections are included, the tree level parameters are deformed, even inducing the mixing among Higgs with different *CP* properties . In the decoupling limit and up to very small radiative corrections, however, the mixing angle $`\alpha `$ reaches the limit $`\beta \pi /2`$.
(b). The couplings of neutral vector gauge bosons ($`\gamma `$ and $`Z^0`$) to tau slepton pairs are written as

$`_{V\stackrel{~}{\tau }_i^{}\stackrel{~}{\tau }_j}={\displaystyle \frac{ig}{c_W}}\left[s_Wc_W\delta _{ij}A^\mu +\left(s_W^2\delta _{ij}+{\displaystyle \frac{1}{2}}(U_{1i}^{})^{}U_{1j}^{}\right)Z^\mu \right]\stackrel{~}{\tau }_i^{}\underset{\mu }{\overset{}{}}\stackrel{~}{\tau }_j.`$ (15)
(c). The couplings of neutral vector gauge bosons to charged lepton pairs are described by the same one as in the SM

$`_{V\overline{f}f}={\displaystyle \frac{g}{c_W}}\left\{s_Wc_WA^\mu \overline{f}\gamma _\mu f+Z^\mu \overline{f}\gamma _\mu \left[\left({\displaystyle \frac{1}{2}}s_W^2\right){\displaystyle \frac{1\gamma ^5}{2}}s_W^2{\displaystyle \frac{1+\gamma ^5}{2}}\right]f\right\}.`$ (16)
(d). The couplings of a neutral Higgs boson to a pair of massive neutral gauge bosons is given by

$`_{h^0ZZ}={\displaystyle \frac{gm_Z}{2c_W}}\mathrm{sin}(\beta \alpha )h^0Z^\mu Z_\mu .`$ (17)
Using Feynman rules, the Dalitz plot density of the production cross section can be obtained. In terms of the energies $`E_1`$, $`E_2`$, $`E_h`$ of the three final state particles $`\stackrel{~}{\tau }_1^{}`$, $`\stackrel{~}{\tau }_1^+`$ and $`h^0`$, respectively, it is written as

$`{\displaystyle \frac{\text{d}\sigma }{\text{d}E_1\text{d}E_2}}=`$ $`{\displaystyle \frac{m_e^2}{8(2\pi )^3\sqrt{s^24sm_e^2}}}\{{\displaystyle \frac{s}{8}}|A|^2+{\displaystyle \frac{\sqrt{s}E_h}{4}}\text{Re}(AB^{})`$
$`+[{\displaystyle \frac{m_h^2}{8}}+({\displaystyle \frac{1}{8}}{\displaystyle \frac{1}{2}}s_W^2+s_W^4)(E_h^2m_h^2)({\displaystyle \frac{s}{3m_e^2}}+{\displaystyle \frac{2}{3}})]|B|^2\},`$ (18)
where as usual $`s`$ is the total center-of-mass energy squared, the coefficients $`A`$ and $`B`$ are defined as follows:

$`A=`$ $`\left({\displaystyle \frac{g}{c_Wm_Z}}\right)^2[g_{11}({\displaystyle \frac{|U_{11}^{}|^2}{2}}s_W^2)({\displaystyle \frac{1}{s2\sqrt{s}E_1}}+{\displaystyle \frac{1}{s2\sqrt{s}E_2}})`$
$`+g_{12}{\displaystyle \frac{(U_{12}^{})^{}U_{11}^{}}{2}}{\displaystyle \frac{1}{s2\sqrt{s}E_2+m_1^2m_2^2}}`$
$`+g_{21}{\displaystyle \frac{(U_{11}^{})^{}U_{12}^{}}{2}}{\displaystyle \frac{1}{s2\sqrt{s}E_1+m_1^2m_2^2}}`$
$`+i({\displaystyle \frac{|U_{11}^{}|^2}{2}}s_W^2){\displaystyle \frac{g\mathrm{sin}(\beta \alpha )(s\sqrt{s}E_hm_Z^2)}{c_Wm_Z(sm_Z^2)}}],`$
$`B=`$ $`i\left({\displaystyle \frac{|U_{11}^{}|^2}{2}}s_W^2\right){\displaystyle \frac{g^3\mathrm{sin}(\beta \alpha )}{c_W^3m_Z(sm_Z^2)}}.`$ (19)
To obtain the total production cross section $`\sigma (e^+e^{}h^0\stackrel{~}{\tau }_1^+\stackrel{~}{\tau }_1^{})`$, one can rewrite equation (14) in the standard form for the Dalitz plot density , and then integrate over the variables $`s_{12}`$, $`s_{23}`$ whose integration bounds are given from kinematics. For a given center of mass energy the cross section depends on the masses of created particles and vanishes at the kinematic boundary.

III. NUMERICAL RESULTS

In this section, we will illustrate our numerical results of the total production cross section for $`e^+e^{}h^0\stackrel{~}{\tau }_1^+\stackrel{~}{\tau }_1^{}`$ based on a specific constrained MSSM scenario for the relevant SUSY parameters. The whole parameters appearing in our analyses contain $`|A_\tau |,|\mu |,\phi ,\mathrm{tan}\beta ,m_0,m_{1/2},m_h`$. Including the two-loop radiative corrections, the mass of the lightest Higgs boson, $`m_h`$, depends mainly on the following variables : the top quark mass, the top and bottom squark masses, the mixing in the stop sector, the pseudoscalar Higgs mass $`m_A`$ and $`\mathrm{tan}\beta `$. In particular scenarios, like in SUGRA or gauge and anomaly mediation , some of the soft parameters may be related to each other. However, since the mechanism of SUSY breaking is unknown, the soft SUSY breaking parameters associated with the tau slepton sector can be considered to be independent of the parameters in the Higgs and squark sector . In what follows, therefore, we consider the former seven parameters as free phenomenological parameters to be constrained by experimental and theoretical considerations. Among the constraints that we are going to impose on the MSSM are those which follow from the comparison of the SM with experimental data, from the experimental limits on the masses of yet unobserved particles , etc, and also those that follow from the ideas of unification and from SUSY GUT models .

In the following numerical analysis, we fix the center of mass energy to $`\sqrt{s}=1`$ TeV, and then demonstrate in turn the dependence of the cross section on the choice of parameters. In Fig. 2 and Fig. 3, the cross section is shown as a function of the unified scalar fermion mass $`m_0`$ for the case of low and high $`\mathrm{tan}\beta `$ value: $`\mathrm{tan}\beta =5`$ and $`\mathrm{tan}\beta =40`$, respectively. Three curves in each figure correspond respectively to three chosen representative values of the unified gaugino mass $`m_{1/2}`$ . The other parameters are typically taken, for example, as the following: $`|A_\tau |=1`$ TeV, $`\phi =\pi /4`$, $`m_h=120`$ GeV, $`|\mu |=500`$ GeV (for $`\mathrm{tan}\beta =5`$), $`|\mu |=300`$ GeV (for $`\mathrm{tan}\beta =40`$). As can be seen, a typical magnitude of the cross section is $`0.0020.03`$ fb. The variation of magnitude depends mainly on $`m_0`$ . As the value of $`m_0`$ is increased, the cross section drops rather dramatically. The values of $`m_{1/2}`$ do not change the whole trends of the plots, but can slightly shift the value of the cross section in the same order of magnitude. This is due to the fact that the tau slepton masses basically depend on $`m_0`$ . In addition, the value of $`\mathrm{tan}\beta `$ can also influence the value of the cross section obviously, which can be seen from the position interchange of the solid line and the long-dashed line in the two graphs, in particular, for the large $`\mathrm{tan}\beta `$, smaller $`m_0`$ and $`m_{1/2}`$ are excluded by experimental limits. As the value of $`\mathrm{tan}\beta `$ increases from low to large, the cross section decreases significantly. However, for no more than about 200 GeV value of $`m_0`$ in the case of the low $`\mathrm{tan}\beta `$ (about 150 GeV in the case of the large $`\mathrm{tan}\beta `$), the production cross section can exceed the value $`\sigma (e^+e^{}h^0\stackrel{~}{\tau }_1^+\stackrel{~}{\tau }_1^{})0.01`$ fb. With a integrated luminosity 500 $`\text{fb}^1`$, this provides several dozen events in a few years, which should be sufficient as a sample to isolate the final state and measure the $`g_{h^0\stackrel{~}{\tau }_1^{}\stackrel{~}{\tau }_1}`$ coupling with some accuracy.

In the Fig. 4, the cross section is displayed as a function of the Higgs mass mixing parameter $`|\mu |`$ for three representative values of the lightest Higgs boson mass $`m_h`$ , and the other parameters are fixed to $`\mathrm{tan}\beta =5`$, $`|A_\tau |=1`$ TeV, $`\phi =\pi /4`$, $`m_0=200`$ GeV, $`m_{1/2}=250`$ GeV. The value of $`|\mu |`$ is constrained by the requirement of radiative electroweak symmetry breaking , while the values of $`m_h`$ are chosen from the current experimental and theoretical considerations . The plots show that the cross section changes somewhat large with values of $`|\mu |`$ but is not too sensitive to the values of $`m_h`$. For a larger value of $`|\mu |`$, the cross section is smaller. When the value of $`m_h`$ is increased by 20 GeV, the cross section drops approximately 0.001 fb.

The *CP* violating phases in $`A_\tau `$, arg($`A_\tau `$), can modify the masses and the mixing parameters of stau as well as the Higgs-stau\*-stau couplings. This would also have large impact on the production cross section of the $`h^0\stackrel{~}{\tau }_1^+\stackrel{~}{\tau }_1^{}`$ final state since it is directly proportional to the square of the $`h^0\stackrel{~}{\tau }_1^+\stackrel{~}{\tau }_1^{}`$ couplings. In Fig. 5, the cross section is shown as a function of the *CP*-violating phase $`\phi `$ (=arg($`A_\tau `$)) for three representative values of the soft breaking trilinear stau coupling $`|A_\tau |`$ , and the other parameters are fixed to $`\mathrm{tan}\beta =5`$, $`|\mu |=500`$ GeV, $`m_h=120`$ GeV, $`m_0=200`$ GeV, $`m_{1/2}=250`$ GeV. The graph shows that the cross section has a symmetry about $`\phi =\pi `$. This results from that the stau mass $`m_{\stackrel{~}{\tau }_1}`$ are symmetric under $`\phi \phi 2\pi `$. Furthermore, the cross section reaches a maximum near $`\phi =\pi `$. This is because the coupling $`g_{h^0\stackrel{~}{\tau }_1^{}\stackrel{~}{\tau }_1}`$ has a maximum at $`\phi =\pi `$ and its phase dependence is stronger than that of the masses of stau. However, the trends of plots and the varying magnitude of the cross section are obviously subject to the values of $`|A_\tau |`$, in particular, it should be noted that the cross section might increase or also decrease with the values of $`|A_\tau |`$ in deferent interval of the phase $`\phi `$.

IV. CONCLUSIONS

In summary, in the framework of the *CP* violating MSSM, we have discussed the production of the lightest neutral Higgs particle in association with the lighter tau slepton pair, $`e^+e^{}h^0\stackrel{~}{\tau }_1^+\stackrel{~}{\tau }_1^{}`$, at future high-energy $`e^+e^{}`$ linear colliders. Simple analytical formulae for the cross section of the three-body process have been given. In a specific constrained MSSM scenario, where the lightest Higgs boson is in the decoupling regime, we have analyzed in detail the dependence of the cross section on the relevant SUSY parameters. On the one hand, for large values of $`|A_\tau |`$, $`|\mu |`$ and $`\mathrm{tan}\beta `$, the mixing between tau sleptons can be very large, this results in the lighter stau can not only be rather light but also at the same time their couplings to Higgs bosons can be strongly enhanced. As a result, the cross section can be very substantial. On the other hand, the effect of the *CP* violating phase in $`A_\tau `$ on the cross section is also demonstrated. It can induce a sizable change of the cross section through the significant dependence of the stau masses and maxing as well as the $`h^0\stackrel{~}{\tau }_1^+\stackrel{~}{\tau }_1^{}`$ couplings on the phase. At $`e^+e^{}`$ colliders with c.m.energy around $`\sqrt{s}=1`$ TeV and with very high luminosity $`\text{d}t500`$ $`\text{fb}^1`$, the typical cross section is of the order of $`𝒪(0.0020.03)`$ fb depending mainly on the universal scalar fermion mass $`m_0`$, in particular for light $`m_0`$ ($`m_0150`$ GeV), the cross section can exceed the level of a 0.01 fb, leading to several dozen events in a few years. In the case where the tau slepton decays into a neutrino and a chargino, the final state topology should be easy to be seen experimentally owing to the clean environment of these colliders. Therefore, all of these results can be expected to be detectable at future $`e^+e^{}`$ linear colliders. Analyzing these final states will allow to measure this important coupling and determine the value of the *CP* violating phase, thus opening a window to probe directly some of the soft-SUSY breaking scalar potential and also to shed some light on the SUSY breaking mechanism.

ACKNOWLEDGMENTS

One of the authors, W. M. Yang, thanks A. Djouadi for useful comments very much, and thanks M. Z. Yang for helpful discussions. This work is in part supported by National Natural Science Foundation of China.