# Triple and Quartic Interactions of Higgs Bosons in the General Two-Higgs-Doublet Model

## 1 Introduction

A particularly simple extension of the Standard Model containing two scalar doublets  has been very extensively investigated in the framework of minimal supersymmetry. In order to cancel gauge anomalies introduced by the fermionic superpartners of gauge bosons and to generate masses of up- and down- quarks in a consistent manner two doublets of Higgs fields are necessary.

Soft supersymmetry-breaking terms  introduce large radiative corrections to the tree-level Higgs boson masses and couplings  and the effective lagrangian of the Higgs sector at the electroweak scale does not satisfy the supersymmetry constraints valid at the SUSY scale. In the most general case when the supersymmetry scale and the scale of heavy Higgs boson mass (usually defined by the mass of the pseudoscalar) are different, the effective theory at the electroweak scale is a two-Higgs-doublet model, where the self-interaction couplings are defined by the renormalization group evolution of the supersymmetric potential couplings from the SUSY scale down to the electroweak scale .

The investigation of direct phenomenological consequences of a two-doublet Higgs sector at a future high luminosity colliders, such as LHC and TESLA, could provide a possibility to study in detail the structure of effective Higgs potential, mass spectrum and couplings of the scalar particles. As usual, the variety of channels where scalars could be produced individually or in association with vector bosons requires a systematical calculation in order to find out what particular channels could have a sufficient counting rates for experimental detection at a given collider luminosity.

We propose a convenient compact form of Feynman rules for a general two-Higgs-doublet model that can be used in the following systematical study of the Higgs boson production channels and use these rules for the calculation of two and three Higgs boson production at a high energy $`e^+e^{}`$ collider.

## 2 Diagonalisation of the mass matrix in the general two-Higgs-doublet model

General form of the (nonsupersymmetric) $`SU(2)\times U(1)`$ invariant potential in the case of two doublets of complex scalar fields $`\phi _1`$, $`\phi _2`$ can be found in 

$`V(\phi _1,\phi _2)=`$ $`\lambda _1(\phi _1^+\phi _1{\displaystyle \frac{v_1^2}{2}})^2+\lambda _2(\phi _2^+\phi _2{\displaystyle \frac{v_2^2}{2}})^2`$   
$`+\lambda _3[(\phi _1^+\phi _1{\displaystyle \frac{v_1^2}{2}})+(\phi _2^+\phi _2{\displaystyle \frac{v_2^2}{2}})]^2`$   
$`+\lambda _4[(\phi _1^+\phi _1)(\phi _2^+\phi _2)(\phi _1^+\phi _2)(\phi _2^+\phi _1)]`$   
$`+\lambda _5[Re(\phi _1^+\phi _2){\displaystyle \frac{v_1v_2}{2}}Re(e^{i\xi })]^2+\lambda _6[Im(\phi _1^+\phi _2){\displaystyle \frac{v_1v_2}{2}}Im(e^{i\xi })]^2`$   
where $`\lambda _i`$ are real constants. Components of scalar doublets $`\phi _{1,2}`$ are

$$\phi _1=\{iw_1^+,\frac{1}{\sqrt{2}}(v_1+h_1+iz_1)\},\phi _1=\{iw_2^+,\frac{1}{\sqrt{2}}(v_2+h_2+iz_2)\}.$$
(2)   
where $`w`$ is a complex field and $`z`$, $`h_{1,2}`$ are real scalar fields. Vacuum expectation values $`v_1`$, $`v_2`$ correspond to the minimum of the potential

$$\phi _1=\frac{1}{\sqrt{2}}\{0,v_1\},\phi _2=\frac{1}{\sqrt{2}}\{0,v_2e^{i\xi }\}$$
(3)   
where the phase $`\xi `$ can be removed by the rotation of $`\phi _1^+\phi _2`$ not affecting the $`\lambda _4`$ term in (1). Substiution of (2) to (1) gives a bilinear form of the mass term with mixed components $`w,h_{1,2},z`$, which can be diagonalized by an orthogonal transformation of the fields in order to define the tree level masses of physical bosons. The resulting spectrum of scalars consists of two charged $`H^\pm `$, three neutral $`h`$, $`H`$, $`A^0`$ scalar fields, and three Goldstone bosons $`G`$. This procedure is described in many papers (for instance, ). The $`w_{1,2}`$ sector is diagonalized by the rotation of $`w_1,w_2H,G`$

$$w_1^\pm =H^\pm s_\beta +G^\pm c_\beta ,w_2^\pm =H^\pm c_\beta +G^\pm s_\beta $$
(4)   
defined by the angle

$$\mathrm{𝚝𝚐}\beta =\frac{v_2}{v_1}$$
(5)   
and leading to the massless $`G`$ field and the field of massive charged Higgs boson $`H^\pm `$, $`m_{H\pm }^2=\lambda _4(v_1^2+v_2^2)/2`$. The $`z_{1,2}`$ sector is diagonalized by the rotation $`z_1,z_2A^0,G^{^{}}`$ defined by the angle $`\beta `$ (4) and giving again one massless field $`G^{^{}}`$ and the field of CP-odd Higgs boson $`A^0`$ with the mass $`m_A^2=\lambda _5(v_1^2+v_2^2)/2`$. Finally, the $`h_1,h_2`$ sector is diagonalised by the rotation $`h_1,h_2h,H`$ defined by the angle $`\alpha `$

$$\mathrm{𝚜𝚒𝚗}2\alpha =\frac{2m_{12}}{\sqrt{(m_{11}m_{22})^2+4m_{12}^2}},\mathrm{𝚌𝚘𝚜}2\alpha =\frac{m_{11}m_{22}}{\sqrt{(m_{11}m_{22})^2+4m_{12}^2}}$$
(6)   
where

$`m_{11}=`$ $`{\displaystyle \frac{1}{4}}[4v_1^2(\lambda _1+\lambda _3)+v_2^2\lambda _5]`$   
$`m_{22}=`$ $`{\displaystyle \frac{1}{4}}[4v_2^2(\lambda _2+\lambda _3)+v_1^2\lambda _5]`$   
$`m_{12}=`$ $`{\displaystyle \frac{1}{4}}(4\lambda _3+\lambda _5)v_1v_2`$   
giving two massive fields of CP-even Higgs bosons $`H,h`$ with the mass values

$$m_{H,h}^2=m_{11}+m_{22}\pm \sqrt{(m_{11}m_{22})^2+4m_{12}^2}$$
(7)   
In the explicit form the diagonal mass matrix of scalar fields and the physical boson interaction vertices are obtained after the following substitution of $`\lambda _i`$ to the potential $`V(\phi _1,\phi _2)`$ (1):

$`\lambda _1=`$ $`{\displaystyle \frac{1}{2v^2}}{\displaystyle \frac{1}{c_\beta ^2}}[{\displaystyle \frac{s_\alpha }{s_\beta }}c_{\alpha \beta }m_h^2{\displaystyle \frac{c_\alpha }{s_\beta }}s_{\alpha \beta }m_H^2]+{\displaystyle \frac{c_{2\beta }}{4c_\beta ^2}}\lambda _5`$ (8)   
$`\lambda _2=`$ $`{\displaystyle \frac{1}{2v^2}}{\displaystyle \frac{1}{s_\beta ^2}}[{\displaystyle \frac{c_\alpha }{c_\beta }}c_{\alpha \beta }m_h^2+{\displaystyle \frac{s_\alpha }{c_\beta }}s_{\alpha \beta }m_H^2]{\displaystyle \frac{c_{2\beta }}{4s_\beta ^2}}\lambda _5`$   
$`\lambda _3=`$ $`{\displaystyle \frac{1}{2v^2}}[{\displaystyle \frac{s_{2\alpha }}{s_{2\beta }}}m_h^2+{\displaystyle \frac{s_{2\alpha }}{s_{2\beta }}}m_H^2]{\displaystyle \frac{1}{4}}\lambda _5`$   
$`\lambda _4=`$ $`{\displaystyle \frac{2}{v^2}}m_{H^\pm }^2`$   
$`\lambda _6=`$ $`{\displaystyle \frac{2}{v^2}}m_{A^0}^2`$   
where we used the notation $`v^2=v_1^2+v_2^2`$, $`s_\alpha =\mathrm{𝚜𝚒𝚗}\alpha `$, $`c_\alpha =\mathrm{𝚌𝚘𝚜}\alpha `$. Diagonalization of the mass term takes place if $`\lambda _5`$ is arbitrary, but the necessary condition for the $`CP`$-invariance of potential (1) is $`\lambda _5=\lambda _6`$ . Unfortunately, after the substitution of (8) to the potential (1) the intermediate expressions for the four scalar boson interaction vertices turn out to be extremely cumberous and it is very difficult to reduce them to some compact convenient form, where the dependence of the coupling from the parameters could be clearly seen. This is a technical problem of the symbolic manipulaion program  that we used. However, symbolic transformations of the intermediate expressions are simpler, if we rewrite the potential $`V(\phi _1,\phi _2)`$ in the oftenly used representation

$`U(\phi _1,\phi _2)=`$ $`\mu _1^2(\phi _1^+\phi _1)\mu _2^2(\phi _2^+\phi _2)\mu _{12}^2(\phi _1^+\phi _2+\phi _2^+\phi _1)`$   
$`+\overline{\lambda }_1(\phi _1^+\phi _1)^2+\overline{\lambda }_2(\phi _2^+\phi _2)^2+\overline{\lambda }_3(\phi _1^+\phi _1)(\phi _2^+\phi _2)`$   
$`+\overline{\lambda }_4(\phi _1^+\phi _2)(\phi _2^+\phi _1)+{\displaystyle \frac{\overline{\lambda }_5}{2}}[(\phi _1^+\phi _2)(\phi _1^+\phi _2)+(\phi _2^+\phi _1)(\phi _2^+\phi _1)]`$   
It is easy to check that in the case of zero $`\phi _1^+\phi _1`$ phase the potentials (1) and (9) are equivalent if the constants $`\overline{\lambda }_i`$, $`\mu `$ and $`\lambda _i`$ are related by the formulas

$`\overline{\lambda }_1=\lambda _1+\lambda _3,\overline{\lambda }_2=\lambda _2+\lambda _3,\overline{\lambda }_3=2\lambda _3+\lambda _4,`$ (10)   
$`\overline{\lambda }_4=\lambda _4+{\displaystyle \frac{\lambda _5}{2}}+{\displaystyle \frac{\lambda _6}{2}},\overline{\lambda }_5={\displaystyle \frac{\lambda _5}{2}}{\displaystyle \frac{\lambda _6}{2}}`$   
and

$$\mu _{12}^2=\lambda _5\frac{v_1v_2}{2},\mu _1^2=\lambda _1v_1^2+\lambda _3v_1^2+\lambda _3v_2^2,\mu _2^2=\lambda _2v_2^2+\lambda _3v_1^2+\lambda _3v_2^2$$
(11)   
The expressions (11) are sometimes called ’minimization conditions’, if one starts from the potential $`U(\phi _1,\phi _2)`$ (9), where the symbolic structure does not show clearly a possible minimum. In the MSSM $`\overline{\lambda }_5=`$0 and it follows that $`\mu _{12}^2`$ is fixed and equal to $`m_A^2s_\beta c_\beta `$. If this equality is not satisfied (or, equivalently, $`\lambda _5\lambda _6`$ in (1)), $`CP`$-violation in the Higgs sector can be introduced. The diagonal form of $`U(\phi _1,\phi _2)`$ and the physical scalar boson interaction vertices are obtained by the substitution of the following expressions for $`\overline{\lambda }_i`$ and $`\mu _i`$ in the potential (9):

$`\overline{\lambda }_1=`$ $`{\displaystyle \frac{1}{2v^2}}[({\displaystyle \frac{s_\alpha }{c_\beta }})^2m_h^2+({\displaystyle \frac{c_\alpha }{c_\beta }})^2m_H^2{\displaystyle \frac{s_\beta }{c_\beta ^3}}\mu _{12}^2]`$ (12)   
$`\overline{\lambda }_2=`$ $`{\displaystyle \frac{1}{2v^2}}[({\displaystyle \frac{c_\alpha }{s_\beta }})^2m_h^2+({\displaystyle \frac{s_\alpha }{s_\beta }})^2m_H^2{\displaystyle \frac{c_\beta }{s_\beta ^3}}\mu _{12}^2]`$   
$`\overline{\lambda }_3=`$ $`{\displaystyle \frac{1}{v^2}}[2m_{H^\pm }^2{\displaystyle \frac{\mu _{12}^2}{s_\beta c_\beta }}{\displaystyle \frac{s_{2\alpha }}{s_{2\beta }}}(m_H^2m_h^2)]`$   
$`\overline{\lambda }_4=`$ $`{\displaystyle \frac{1}{v^2}}({\displaystyle \frac{\mu _{12}^2}{s_\beta c_\beta }}+m_A^22m_{H^\pm }^2)`$   
$`\overline{\lambda }_5=`$ $`{\displaystyle \frac{1}{v^2}}({\displaystyle \frac{\mu _{12}^2}{s_\beta c_\beta }}m_A^2)`$   
$`\mu _1^2=`$ $`{\displaystyle \frac{1}{2}}[{\displaystyle \frac{s_\alpha }{c_\beta }}s_{\alpha \beta }m_h^2+{\displaystyle \frac{c_\alpha }{c_\beta }}c_{\alpha \beta }m_H^22tg\beta \mu _{12}^2]`$   
$`\mu _2^2=`$ $`{\displaystyle \frac{1}{2}}[{\displaystyle \frac{c_\alpha }{s_\beta }}s_{\alpha \beta }m_h^2+{\displaystyle \frac{s_\alpha }{s_\beta }}c_{\alpha \beta }m_H^22ctg\beta \mu _{12}^2]`$   
Our expressions for $`\overline{\lambda }_4`$ and $`\overline{\lambda }_5`$ are the same as given in  for the case of zero $`\lambda _6`$ and $`\lambda _7`$. Complete sets of Feynman rules (unitary gauge) for the triple and quartic Higgs boson interactions in the general two-Higgs-doublet model with a possibility of $`CP`$-violation in the Higgs sector (defined by $`\mu _{12}`$ parameter), are shown in Tables 1-2. These sets were obtained by means of LanHEP package . LanHEP package is a specialized symbolic manipulation system capable to generate Feynman rules for the $`SU(2)`$, $`SU(3)`$ gauge invariant lagrangians with arbitrary sets of particle multiplets, in the standard input lagrangian format of CompHEP package . We do not show here a rather long set of Feynman rules in the ’tHooft-Veltman gauge, that can be also generated after the introduction of ghost and ghost-goldstone lagrangian terms to LanHEP program. <sup>1</sup><sup>1</sup>1The generation process takes 15 sec. of CPU time (i686). Complete lagrangian tables in CompHEP format and LanHEP package are available at http://theory.npi.msu.su/~semenov/lanhep.html

We assume that in the Yukawa sector $`<\phi _1>`$ couples only to down fermions

$$V_{ud}\frac{em_d}{2\sqrt{2}m_Ws_Wc_\beta }[\overline{\psi }_1(1+\gamma _5)\psi _2\phi _1+\overline{\psi }_2(1\gamma _5)\psi _1\phi _1^+]$$
(13)   
(here for the $`u`$, $`d`$ quarks $`\overline{\psi }_1=\{\overline{u},V_{ud}\overline{d}+V_{us}\overline{s}+V_{ub}\overline{b}\},\psi _2=d`$ and analogous structures for $`s`$,$`b`$ quarks and leptons, in the case of quarks $`V_{ab}`$ denotes the CKM matrix elements), and $`<\phi _2>`$ couples only to up fermions (so-called model of type II ):

$$\frac{em_u}{2\sqrt{2}m_Ws_Ws_\beta }[\overline{\psi }_1(1+\gamma _5)i\tau _2\psi _2\phi _2^++\overline{\psi }_2(1\gamma _5)i\tau _2\psi _1\phi _2]$$
(14)   
(here $`\overline{\psi }_1=\{\overline{u},V_{ud}\overline{d}+V_{us}\overline{s}+V_{ub}\overline{b}\},\psi _2=u`$ and analogous structures for $`c`$ and $`t`$ quarks). Higgs-gauge boson interaction is defined by the straightforward extension of the covariant derivative in the case of two scalar doublets. It is easy to find the relation between the vacuum expectation values of the potential and the $`W`$-boson mass and coupling $`g=e/\mathrm{𝚜𝚒𝚗}\vartheta _W`$

$$v^2=v_1^2+v_2^2=\frac{4m_W^2}{e^2}s_W^2$$
(15)   
following from the structure of scalar fields kinetic term $`D_\mu \phi D^\mu \phi `$.

From the phenomenological point of view the general multiparametric two-Higgs-doublet model is too flexible to be systematically used for data analysis. Practically no limits on the masses of individual scalars can be set if their couplings to gauge bosons and fermions depend on some free parameters, and can be very small in a rather large regions of parameter space. Recent discussion of the possible limits can be found in . However, the parameter space can be strongly restricted by the constraints imposed by the supersymmetry.

Let us consider the reduction of the general two-doublet model Feynman rules shown in Tables 1,2 to the case of minimal supersymmetry model (MSSM). The potential $`V(\phi _1,\phi _2)`$ (1) contains eight parameters: two VEV’s $`v_1`$, $`v_2`$ and six $`\lambda _i`$ ($`i`$=1,…6). Eight parameters of the potential $`U(\phi _1,\phi _2)`$ (9) $`\mu _1`$, $`\mu _2`$, $`\mu _{12}`$ and $`\overline{\lambda }_i`$ ($`i`$=1,…5) can be found using (10),(11). From the other side, in order to define the Higgs sector we need eight physical parameters: the mixing angle $`\beta `$ and $`W`$-boson mass $`m_W`$, mixing angle $`\alpha `$, the parameter $`\mu _{12}`$ and four masses of scalars $`m_h`$, $`m_H`$, $`m_A`$, $`m^\pm `$. Two VEV’s can be expressed through $`m_W`$, $`\mathrm{𝚝𝚐}\beta `$ by (5) and (15) and only one degree of freedom remains here. In the case of superpotential five additional constraints are imposed, relating all Higgs boson self couplings $`\overline{\lambda }_i`$, ($`i`$=1,…5) to the gauge coupling constants at the energy scale $`M_{SUSY}`$ :

$$\overline{\lambda }_1=\overline{\lambda }_2=\frac{g^2+g_1^2}{8},\overline{\lambda }_3=\frac{g^2g_1^2}{4},\overline{\lambda }_4=\frac{g^2}{2},\overline{\lambda }_5=0.$$
(16)   
As we already noticed, if $`\overline{\lambda }_5=`$0, $`\mu _{12}`$ is fixed and $`CP`$-parity is conserved. The remaining two independent parameters may be used to define all Higgs boson masses and mixing angles. One can choose, for instance, $`r_1,r_2`$ parametrization  ($`r_{1,2}=m_{h,H}^2/m_Z^2`$) or the well-known $`m_A`$, $`\mathrm{𝚝𝚐}\beta `$ parametrization. In order to reduce the general two-Higgs-doublet model vertices to the case of MSSM it is convenient to use the $`\alpha `$, $`\beta `$ parametrization:

$`m_h^2=m_Z^2\mathrm{𝚌𝚘𝚜}2\beta {\displaystyle \frac{\mathrm{𝚜𝚒𝚗}(\alpha +\beta )}{\mathrm{𝚜𝚒𝚗}(\alpha \beta )}},m_H^2=m_Z^2\mathrm{𝚌𝚘𝚜}2\beta {\displaystyle \frac{\mathrm{𝚌𝚘𝚜}(\alpha +\beta )}{\mathrm{𝚌𝚘𝚜}(\alpha \beta )}},`$ (17)   
$`m_A^2=m_Z^2{\displaystyle \frac{\mathrm{𝚜𝚒𝚗}2(\alpha +\beta )}{\mathrm{𝚜𝚒𝚗}2(\alpha \beta )}},\mu _{12}^2=m_A^2\mathrm{𝚜𝚒𝚗}\beta \mathrm{𝚌𝚘𝚜}\beta .`$   
Substitution of these expressions to the vertex factors in Tables 1,2 after trivial trigonometric transformations reduces them to simpler MSSM factors . Complete list of Feynman rules at the MSSM scale is shown in Table 3.

Renormalization group (RG) evolution of the coupling constants $`\lambda _i`$ from the energy scale $`M_{SUSY}`$ to the electroweak scale $`M_{EW}`$ violates the constraints (16)  and the effective low energy potential at the scale $`M_{EW}`$ is the potential of a general two-Higgs doublet model with RG evolved couplings $`\overline{\lambda }_i`$. At a given values of $`m_A`$, $`\mathrm{𝚝𝚐}\beta `$ (or $`\alpha `$, $`\beta `$), masses of Higgs bosons and the mixing angle $`\alpha `$ (or $`m_A`$) at the energy scale $`M_{SUSY}`$ can be obtained using (17). Detailed analysis of the following RG evolution and the calculation of leading-logarithmic radiative corrections to the mixing angles, masses and couplings of Higgs bosons can be found, for instance, in . We briefly point out that the additional input parameters to be defined in order to fix the scheme are the scale of SUSY breaking $`M_{SUSY}`$, the mass parameter in higgsino-gaugino sector $`\mu `$, and the squark mixing parameters $`A`$.

## 3 Multiple production of neutral Higgs bosons

The processes of multiple neutral Higgs boson production in the MSSM were considered in  in the framework of effective potential approach to the calculation of radiatively corrected scalar masses and couplings  of the SUSY Higgs sector. The reactions

$$e^+e^{}hhZ,e^+e^{}hhA,e^+e^{}\nu _e\overline{\nu }_ehh$$
(18)   
were considered and it was shown that the cross sections of double and triple Higgs boson production are not small and the experimental measurements of triple Higgs boson couplings are realistic.

We used the results of  to calculate the radiatively corrected masses of Higgs bosons and mixing angle $`\alpha `$ in the renormalization group approach to the Higgs potential couplings evolution from the SUSY scale $`M_{SUSY}`$ down to the electroweak scale (see also ). We set the $`M_{SUSY}=`$1 TeV and have not included the effects of squark mixing by setting the parameters $`A`$ and $`\mu `$ equal to zero. In the case of not too large $`\mathrm{𝚝𝚐}\beta `$ (we used $`\mathrm{𝚝𝚐}\beta =3`$) and the pseudoscalar mass $`m_A`$ of order 150–250 GeV, masses of heavy CP-even Higgs boson $`H`$ and charged Higgs boson $`H^\pm `$ are also at the scale 150-250 GeV. The lightest Higgs boson mass is approaching 100 GeV when the pseudoscalar mass surpasses 200 GeV. (Changes of the SUSY scale and mixing parameters can in principle shift $`m_h`$ by about 50 GeV, see the details in ). These radiatively corrected parameters were used in our set of Feynman rules. The following calculation of the complete tree level amplitude for the multiple Higgs boson production processes (18) was performed by means of CompHEP package , when the exact symbolic result for the matrix element squared is converted to FORTRAN code and integrated by multichannel Monte-Carlo method. The s-channel resonant peaks of the amplitude (see Fig.1) are regularized by phase space mappings  to ensure an efficient application of VEGAS integrator .

While in the Standard Model the cross section of $`hhZ`$ production  is of order $`210^1`$ fb at the Higgs boson mass 100 GeV and slowly decreasing when the mass of Higgs boson increases, the picture in the two-doublet MSSM sector is strongly changed by the availability of resonant production mechanisms, when the decays of on-shell $`Hhh`$ and $`A^0Zh`$ become possible. We show the dependence of total cross sections in the channels (18) from the masses of CP-even states $`m_h`$ and $`m_H`$ in Fig.2,3. In order to understand qualitatively the cross section behaviour we show also the $`h`$, $`H`$ branching ratio dependence (in the two-body decay channels with the contribution greater than 1%) from their masses in Fig.4,5. Rapid decrease of the total rate at $`m_h=`$ 60 GeV ($`m_H=`$ 120 GeV) and rapid increase at $`m_h=`$ 95 GeV ($`m_H=`$ 190 GeV) are directly connected with the resonant threshold of the heavy scalar decay $`Hhh`$ (see diagrams in Fig.1). The channel $`e^+e^{}hhZ`$ receives some enhancement at $`m_h=`$ 95 GeV ($`m_H=`$ 210 GeV) when the resonance threshold $`A^0Zh`$ is opened.

Our results are qualitatively consistent with the results of , where somewhat different regions of the two-Higgs-doublet model parameter space were explored. Radiatively corrected masses are rather sensitive to the input parameter values. At smaller value of $`\mathrm{𝚝𝚐}\beta `$ a mass interval between the closing and opening $`hh`$ thresholds decreases to a few GeV (, $`\mathrm{𝚝𝚐}\beta =`$ 1.5).

The reactions (18) do not include quartic Higgs boson interaction vertices. We calculated the cross-section of the simplest process

$`e^+e^{}hhhZ`$   
(see Fig.2,3), where quartic vertices $`hhhh`$ and $`hhhH`$ participate (21 diagrams in the unitary gauge). In a very limited region of parameter space the reaction has an observable cross-section if the luminosity is high, and the experimental reconstruction of multijet events is very efficient.

## 4 Conclusions

Large increase of the estimate of the possibly achievable integrated luminosity in the next linear colliders (especially L$`=`$ 500 fb<sup>-1</sup>/year for the TESLA project) makes quite realistic the experimental study of Higgs boson self-interaction. Such an investigation is especially interesting if the Higgs sector of the model includes more than one $`SU(2)`$ multiplet. Untrivial spectrum of scalars leads to the resonant multiple Higgs boson production mechanisms, when the final states with 4 or 6 $`b`$-jets from their decays will appear with the cross sections of one-two orders of magnitude greater than in the SM case of only one scalar boson in the Higgs sector.

For a systematical study of various production channels we derive in a compact form a complete set of Feynman rules for the general case of two-Higgs doublet model. We demonstrate that in the case of minimal supersymmetry, when additional constraints are imposed on the general parameter space, the interaction vertices are reduced to the well-known vertices of the MSSM at the scale $`M_{SUSY}`$. Useful connection of LanHEP output in the standard lagrangian format of CompHEP input makes possible the following efficient calculation of various reactions.

Acknowledgements

M.D. is grateful to S.Y.Choi and P.M.Zerwas for useful discussions. The authors would like to express their gratitude to H.S.Song and the Center for Theoretical Physics, Seoul National University, where this work was completed, for hospitality. The authors acknowledge thankfully discussions with the members of Minami-Tateya group (KEK, Tsukuba) and N.Skachkov. The research was partially supported by RFBR grant 96-02-19773a, INTAS grant INTAS-YSF 98-50 and St.-Petersburg University KCFE grant.

Table 1. Triple Higgs boson interaction vertices in the general two Higgs doublet model

.

Table 2. Quartic Higgs boson interaction vertices in the general two Higgs doublet model

Table 3. Triple and quartic Higgs boson interaction vertices at the scale $`M_{SUSY}`$

Figure captions

Fig. 1 Feynman diagrams for the process $`e^+e^{}hhZ`$

Fig. 2 Total cross sections for the reactions $`e^+e^{}hhZ`$, $`e^+e^{}hhA`$, $`e^+e^{}\nu _e\overline{\nu }_ehh`$ and $`e^+e^{}hhhZ`$ versus the mass of light $`CP`$-even Higgs boson at $`\sqrt{s}=`$500 GeV

Fig. 3 Total cross sections for the reactions $`e^+e^{}hhZ`$, $`e^+e^{}hhA`$, $`e^+e^{}\nu _e\overline{\nu }_ehh`$ and $`e^+e^{}hhhZ`$ versus the mass of heavy $`CP`$-even Higgs boson at $`\sqrt{s}=`$500 GeV

Fig. 4 Two-body branching ratios of heavy $`CP`$-even Higgs boson

Fig. 5 Two-body branching ratios of $`CP`$-odd Higgs boson

.