c PROGRAMS FOR ONE-LOOP ON-SHELL CALCULATIONS IN THE MSSM c Masses/vertex diagrams calculations: P.H.Chankowski, S.Pokorski, J.Rosiek c e-mail: rosiek@fuw.edu.pl c chank@fuw.edu.pl c Box contribution calculations: V. Driesen, W. Hollik c e-mail: vd@itpaxp1.physik.uni-karlsruhe.de c hollik@cernvm.cern.ch c Example of physical masses, cross sections and decay widths calculations program mssm_example implicit double precision (a-h,o-z) logical ino_check dimension hgam_arr(100),hbr_arr(100) dimension agam_arr(100),abr_arr(100) dimension tb_arr(10) dimension crf(12),errout(12) common/hm_phys/frm(2),frm2(2) open(1,file='example.new',status='unknown') c Center of mass energy (LEP 200 here, can be anything) s = 210.d0**2 c Accuracy of physical Higgs mass calculations [GeV^2] errin = 1.d-2 c Additional parameter (step in numerical differentiation) zeps = 1.d-6 c Set parameters now. Most of the parameters is initialized in block c datas in file MH_INIT.FOR - look into it if you want c change something like b quark mass, Z mass etc. Below free parameters c are initialized. Convention used for the MSSM parameters and Feynman c rules are described in the paper Phys. Rev. D41, 3464 (1990) c Pseudoscalar mass [GeV] pm = 50 c tan(beta) (don't put exactly 1 - m_h = 0 can result in program crash) tanb = 5.d0 c Top mass [GeV] tm = 175 c Squark and slepton mass parameter (for definition see SET_PAR subroutine c and Phys.Rev D41,3464 (1990) - the same for other parameters below) amsq = 1000 amsl = 300 c Gaugino mass parameter c CAUTION: amg = 250 below is equivalent to gaugino mass 500 in c more frequently used convention (we use 2*amg as the soft gaugino c mass term in the lagrangian) amg = 250 c Mu parameter in Higgs potential am = 200 c Mixing parameter in sfermion sector (A parameter) ys = 0.5d0 c Subroutine set_par has been designed for simplified initialization c of SUSY parameters. Program is however very general and can c incorporate much more general set of input parameters, for example c independent mass parameters for binos, winos and gluinos, c independent non-diagonal entries in soft sfermion mass matrices, c complex CP breaking miu parameter in Higgs potential etc. If you c want to perform such calculations, don't use set_par but initialize c yourself necessary commons. To understand how to do it, look into c SET_PAR.FOR and MH_DIAG.FOR files and the paper in PRD41 listed above. call set_par(tanb,pm,tm,am,amsq,amsl,amg,ys,ierr) c If ierr<>0 something is wrong - negative sfermion masses or very c light -ino masses (below 1 MeV). Possible ierr values are described c in set_par header. Choose more reasonable parameter set and c try once more (actually very low gaugino masses are unreasonable, c but should not cause the program crash) if (ierr.ne.0) stop c Now initialize renormalization constants in gauge and Higgs sectors call gr_const call hr_const c Calculate physical one-loop Higgs boson masses. They are stored c after this procedure in common/hm_phys/. Array frm gives masses, c frm2 mass squares (m_H = frm(1), m_h = frm(2)). Lighter scalar mass c square for some parameters sets may be negative - this is not a c numerical error, but the sign that the initial choice of vacuum was bad. call hm_solve(ierr,errin,errout) c If ierr>0, error in physical Higgs masses calculations is greater c then required (and equal to errout). if (ierr.gt.0) stop c For some cases (very rare - happens sometimes for 80 < M_A < 150, c tan(beta)>25) physical masses calculated in hm_phys are reversed: c m_h > m_H. Procedure switches them automatically m_H <-> m_h, c but cannot decide if to switch also the couplings associated to c the given scalar, so then there is a danger of using object with c mass of h and couplings of H etc. In such cases ierr flag is set c to ierr = - 1. One can treat it as a warning or fatal error, c depending how careful one wants to be. if (ierr.eq.-1) write(1,'(a)') 1 ' Possible error in mass ordering: m_h > m_H?' c Stop here if you need masses only. The rest of the code calculates c the cross section and branching ratios write(1,99)'tan(beta)=',tanb,'M_A=',pm,'M_h=',frm(2),'M_H=',frm(1) c Also h (or H) mass can be used directly as input. c Before choosing some mh (mH) value one can wish to know allowed mass c ranges for the given set of SUSY parameters. This is done by the c following subroutine: c mh_limits(hmin,hmax,pm,tm,am,amsq,amsl,amg,ys,scalar_no,iout) c hmin and hmax return the minimal and maximal scalar masses achievable c for 0.5 m_H?' c And now the examples of cross sections/branching ratios calculations. c Initialize external Z factors on Z, W and Higgs lines. c They play a role of effective (1-loop) alfa angle call z_ext(zeps,ierr) c For very large tan(beta) scalar propagator can have negative residuum. c One need to include second order corrections to cure it. if (ierr.ne.0) stop c Initialize renormalization constants in the fermion sector call fr_const c Check if contribution from gauginos/higgsinos to the invisible Z width c does not exceed LEP limits (warning only) if (ino_check()) write(1,'(a)')' Too light -inos!' c Initialization of box contribution calculation. If you DO NOT CALL c set_par_box procedure (comment it out for example), box diagrams WILL NOT c BE CALCULATED. This leads to some error in final result for cross c sections (typically 3-4% for LEP200 energies, 15% for NLC energies) c but significantly speed up the program: box contribution calculations c are at least 10 times slower all other calculations together c CAUTION: c 1. set_par_box has to be called for each new set of parameters c (set_par reset box calculations flag). c 2. set_par_box can be called only AFTER mass calculations (hm_solve) c 3. Box diagrams are included only in the on-shell cross section c calculations, they are always neglected for LEP100 process ee->Z*h->ffh. call set_par_box(s,tanb,pm,am,amsl,amg) c If box diagrams are switched off, calculations of vertex corrections c take more than 90% of CPU time. Neglecting them can lead typically c to 10-20% error in the cross section values (for LEP200 energies; more c for LEP100 and NLC!). If you do not need full accuracy, you can switch c off also the vertex corrections by uncommenting the line below: c call vert_stat(.false.,.false.) c First argument switches on/off more important ZZS and ZPS vertex c corrections, second argument switches off Zee vertex corrections, c usually small. After calling vert_stat(.false.,.false.) c vertex corrections are calculated only in places where they are c always important (some decay branching ratios). c CAUTION: vertex corrections flags are not reset by set_par, one must c explicitly call vert_stat(.true.,.true.) to switch them on back. c Of course, one can keep only the ZZS/ZPS vertex diagrams: c call vert_stat(.true.,.false.) c Now on-shell cross sections. Higgs bosons are numbered: H has number 1, c h number 2 (CP-even scalars), A number 1, G (Goldstone) number 2 c (CP-odd scalars). Cross sections have arguments: c zzs_cross(scalar_no,center_of_mass_energy) c zps_cross(scalar_no,pseudoscalar_no,center_of_mass_energy) c For example: c e+e- -> Zh crzh = zzs_cross(2,s) c e+e- -> Ah crah = zps_cross(2,1,s) write(1,100)crzh,crah c Bjorken process cross section for Z produced of shell. Scalar c indices as above, arguments of the procedure are: c zzs_cr_virt(scalar_no, c center_of_mass_energy, c array[1..12] with cross section for the 12 possible c fermions in the final state S_i+ff, c assumed relative error of integration, c array[1..12] with achieved relative errors of c integration each final state, c maximal number of iterations during integration) c Example: c LEP center of mass energy: s1 = 91.175**2 c Assumed relative error of integration over Z* energy errin = 1.d-2 c Maximal number of iteration niter = 10 c e+e- -> Z*h crzhv = zzs_cr_virt(2,s1,crf,errin,errout,niter) c Estimation of relative error of integration for total cross c section err = 0 do 10 i=1,12 10 err = err + (crf(i)*errout(i))**2 if (crzhv.ne.0) err = sqrt(err)/crzhv write(1,'(a)')' Total cross section Z*h for E_CM = M_Z^2:' write(1,101)'ee->Z*h=',crzhv,'relative error:',err write(1,'(a)')' Cross sections in [pb] for different channels:' write(1,102)'nunu(e): ',crf(1),errout(1) write(1,102)'nunu(miu):',crf(2),errout(2) write(1,102)'nunu(tau):',crf(3),errout(3) write(1,102)'ee: ',crf(4),errout(4) write(1,102)'mumu: ',crf(5),errout(5) write(1,102)'tautau: ',crf(6),errout(6) write(1,102)'uu: ',crf(7),errout(7) write(1,102)'cc: ',crf(8),errout(8) write(1,102)'tt: ',crf(9),errout(9) write(1,102)'dd: ',crf(10),errout(10) write(1,102)'ss: ',crf(11),errout(11) write(1,102)'bb: ',crf(12),errout(12) write(1,*) c Decay branching ratios for scalars. c Function gam_s(scalar_no,gam_arr,br_arr) returns c total width of h or H decay, depending on scalar_no argument value. c Partial decay widths and branching ratios are stored in arrays c gam_arr and br_arr respectively. c List of entries in gam_arr and br_arr arrays for scalar decays: c Entry Decay mode c 1: S -> ZZ c 2: S -> Zgamma c 3: S -> gamma + gamma c 4-6: S -> ll c 7-9: S -> uu c 10-12: S -> dd c 13: S -> gg c 14-17: S -> 2 neutralino (identical) c 18-23: S -> neutralino + neutralino (not identical) c 24-25: S -> 2 chargino (identical with opposite charge) c 26: S -> chargino + chargino (not identical) c attention: decays into c1+c2- ad c1-c2+ summed c 27: S -> AA c 28: S -> SS c 29: S -> H+H- c 30: S -> W+H-,W-H+ c 31: S -> W+W- c 32: S -> ZP c 33-35: S -> 2 sneutrino (identical) c 36-41: S -> 2 sleptons (identical) c 42-56: S -> slepton + slepton (not identical) c attention: decays into Li+Lj- and Li-Lj+ summed c 57-62: S -> 2 up-squarks (identical) c 63-77: S -> up-squark + up-squark (not identical) c attention: decays into Ui+Uj- and Ui-Uj+ summed c 78-83: S -> 2 down - squarks (identical) c 84-98: S -> down-squark + down-squark (not identical) c attention: decays into Di+Dj- and Di-Dj+ summed c h decays h_total_width = gam_s(2,hgam_arr,hbr_arr) c Finally pseudoscalar decays. c Function gam_p(gam_arr,br_arr) returns total width of A decay. c Partial decay widths and branching ratios are as for scalar stored c in arrays gam_arr and br_arr respectively. CAUTION: decays into c squarks not included (even if kinematically allowed). c List of entries in gam_arr and br_arr arrays for scalar decays: c Entry Decay mode c 1: P -> ZZ c 2: P -> Zgamma c 3: P -> gamma + gamma c 4-6: P -> ll c 7-9: P -> uu c 10-12: P -> dd c 13: P -> gg c 14-17: P -> 2 neutralino (identical) c 18-23: P -> neutralino + neutralino (not identical) c 24-25: P -> 2 chargino (identical with opposite charge) c 26: P -> chargino + chargino (not identical) c attention: decays into c1+c2- ad c1-c2+ summed c 27: P -> W+H-,W-H+ c 28-29: P -> ZS c 30-35: P -> 2 sleptons (identical) c 36-50: P -> slepton + slepton (not identical) c attention: decays into Li+Lj- and Li-Lj+ summed c 51-56: P -> 2 up-squarks (identical) c 57-71: P -> up-squark + up-squark (not identical) c attention: decays into Ui+Uj- and Ui-Uj+ summed c 72-78: P -> 2 down - squarks (identical) c 79-92: P -> down-squark + down-squark (not identical) c attention: decays into Di+Dj- and Di-Dj+ summed c A decays p_total_width = gam_p(agam_arr,abr_arr) write(1,103)hbr_arr(6),hbr_arr(12),abr_arr(6),abr_arr(12), 1 hbr_arr(3),hbr_arr(13),h_total_width,p_total_width c That's all, we hope not very complicated and easy in use. c Basically the library procedures can be used as a black box, c at least if you don't need modify anything yourself. c Be careful - whole program is very complicated and calculations are c time-consuming, so using it call only the procedures which are really c useful for your purposes. c There are some dangers. Due to accidental cancellations, for very c special choices of parameters program crash can occur. This effect c is also machine/compiler dependent: some compilers manage to handle with c such situations, others are not so good. We have some c experience from the long runs scanning over SUSY parameters - it c happened several times per more than 100 000 parameter choices and c it is rather difficult to remove entirely this possibility. c The second problem is theoretical rather than the numerical one. For c pseudoscalar mass in the range (80-140) GeV, large tan(beta) (>20) c and some choices of SUSY parameters 1-loop h and H masses can be c very close each other (specially for large mtop). In such case c mixing on the external Higgs line is very strong and accuracy c of the "external Z factors" calculations (z_ext subroutine) fails. c The same problem can appear for large tan(beta) and M_h (physical) c close to m_H (tree level value). c Numerically it is connected with the fact that the scalar propagator c has the pole of the (almost) second order (normally it has two c separated poles of the first order for m_h and m_H). Program works c but numbers produced for cross sections are meaningless. c Fixing this problem is complicated and requires extremely c complicated two-loop calculations. c Of course reports on any bugs, errors, inaccuracies are appreciated c very much. Please report them to Janusz.Rosiek@fuw.edu.pl close(1) 97 format(1x,'minimal and maximal h mass (GeV):',2(1pe9.2,3x),/) 98 format(1x,'tan(beta) solutions for fixed h mass (mh=', 1 1pe9.2,' GeV):',5(1pe9.2,3x),/) 99 format(1x,2(a,1pe9.2,3x),2(a,1pe11.4,3x),/) 100 format(1x,'Cross sections in [pb] for E_CM = (210 GeV)^2:',/, 1 ' ee->Zh=',1pe10.3,/,' ee->Ah=',1pe10.3,/) 101 format(1x,a,1pe11.4,4x,a,1pe9.2,/) 102 format(1x,'ee->h Z*-> ',a,1pe11.4,4x,'rel. error:',1pe9.2) 103 format(1x,'BR(h->2tau)=',1pe10.3,/,' BR(h->bb)= ',1pe10.3,/, 1 ' BR(A->2tau)=',1pe10.3,/,' BR(A->bb)= ',1pe10.3,/, 2 ' BR(h->2photon)= ',1pe10.3,/, 3 ' BR(h->2gluon)= ',1pe10.3,/, 4 ' h_width= ',1pe12.5,/,' A_width= ',1pe12.5) end