Neutral Higgs Boson Production at e+e- Colliders in the Complex MSSM: A Full One-Loop Analysis

For the search for additional Higgs bosons in the Minimal Supersymmetric Standard Model (MSSM) as well as for future precision analyses in the Higgs sector a precise knowledge of their production properties is mandatory. We evaluate the cross sections for the neutral Higgs boson production at e^+e^- colliders in the MSSM with complex parameters (cMSSM). The evaluation is based on a full one-loop calculation of the production channels e^+e^- ->h_i Z, h_i gamma, h_i h_j (i,j = 1,2,3), including soft and hard QED radiation. The dependence of the Higgs boson production cross sections on the relevant cMSSM parameters is analyzed numerically. We find sizable contributions to many cross sections. They are, depending on the production channel, roughly of 10-20% of the tree-level results, but can go up to 50% or higher. The full one-loop contributions are important for a future linear e^+e^- collider such as the ILC or CLIC. There are plans to implement the evaluation of the Higgs boson production cross sections into the code FeynHiggs.


Introduction
combined with the appropriate effective couplings [32]. The full one-loop corrections in the cMSSM listed here together with resummed SUSY corrections have been implemented into the code FeynHiggs [32][33][34][35][36]. Corrections at and beyond the one-loop level in the MSSM with real parameters (rMSSM) are implemented into the code HDECAY [37,38]. Both codes were combined by the LHC Higgs Cross Section Working Group to obtain the most precise evaluation for rMSSM Higgs boson decays to SM particles and decays to lighter Higgs bosons [39].
The most advanced SUSY Higgs boson production calculations at the LHC are available via the code SusHi [40], which are, however, so far restricted to the rMSSM. However, particularly relevant are higher-order corrections also for the Higgs boson production at e + e − colliders, where a very high accuracy in the Higgs property determination is anticipated [19]. In this paper we concentrate on the neutral Higgs boson production at e + e − colliders in association with a SM gauge boson or another cMSSM Higgs boson, i.e. we calculate (i, j = 1, 2, 3), (2) σ(e + e − → h i γ) . ( The processes e + e − → h i h i and e + e − → h i γ are purely loop-induced. The evaluation of the channels (1) -(3) is based on a full one-loop calculation, i.e. including electroweak (EW) corrections, as well as soft and hard QED radiation. Results for the cross sections (1) - (3) have been obtained over the last two decades. A first (nearly) full calculation of the production channels (1) and (2) in the rMSSM was presented in Ref. [41] (leaving out only a detailed evaluation of the initial state radiation). 1 A tree-level evaluation of the channels (1) and (2) in the cMSSM was presented in Ref. [42], where higher-order corrections were included via effective couplings. Higher-order corrections to the channels (1) and (2) in the cMSSM were given in Ref. [22], where the third generation (s)fermion contributions to the production vertex as well as Higgs boson propagator corrections were taken into account. Another full one-loop calculation of e + e − → hZ was given in Ref. [43]. The production of two equal Higgs bosons in e + e − collisions in the rMSSM, where only box-diagrams contribute, were presented in Refs. [44] and further discussed in Ref. [45]. Finally, the channel (3) in the rMSSM was evaluated in Ref. [46]. A short numerical comparison with the literature will be given in Sect. 4.
In this paper we present a full one-loop calculation for neutral Higgs boson production at e + e − colliders in association with a SM gauge boson or another cMSSM Higgs boson, taking into account soft and hard QED radiation. In Sect. 2 we very briefly review the renormalization of the relevant sectors of the cMSSM. Details about the calculation can be found in Sect. 3. In Sect. 4 various comparisons with results from other groups are given. The numerical results for all production channels (1) -(3) are presented in Sect. 5. The conclusions can be found in Sect. 6. There are plans to implement the evaluation of the production cross sections into the Fortran code FeynHiggs [32][33][34][35][36].
They will be further explained in the text below.

The complex MSSM
The cross sections (1) - (3) are calculated at the one-loop level, including soft and hard QED radiation, see the next section. This requires the simultaneous renormalization of the Higgs and gauge boson sector as well as the fermion sector of the cMSSM. We give a few relevant details about these sectors and their renormalization. More information can be found in Refs. [23,24,[47][48][49][50][51][52][53][54].
The renormalization of the Higgs and gauge-boson sector follow strictly Ref. [47] and references therein (see especially Ref. [32]). This defines in particular the counterterm δt β , as well as the counterterms for the Z boson mass, δM 2 Z , and for the sine of the weak mixing angle, δs w (with s w = 1 − c 2 w = 1 − M 2 W /M 2 Z , where M W and M Z denote the W and Z boson masses, respectively).
The renormalization of the fermion sector is described in detail in Ref. [47] and references therein. For simplification we use the DR renormalization for all three generations of downtype quarks and leptons, in the notation of Ref. [47]:

Calculation of diagrams
In this section we give some details about the calculation of the tree-level and higher-order corrections to the production of Higgs bosons in e + e − collisions. The diagrams and corresponding amplitudes have been obtained with FeynArts (version 3.9) [55], using the MSSM model file (including the MSSM counterterms) of Ref. [47]. The further evaluation has been performed with FormCalc (version 8.4) and LoopTools (version 2.12) [56]. The Higgs sector quantities (masses, mixings,Ẑ factors, etc.) have been evaluated using FeynHiggs [32][33][34][35][36] (version 2.11.0).

Contributing diagrams
Sample diagrams for the process e + e − → h i h j (i, j = 1, 2, 3) are shown in Fig. 1, for the process e + e − → h i Z (i = 1, 2, 3) are shown in Fig. 2 and for the process e + e − → h i γ (i = 1, 2, 3) in Fig. 3. Not shown are the diagrams for real (hard and soft) photon radiation.  Figure 1: Generic tree, self-energy, vertex, box and counterterm diagrams for the process e + e − → h i h j (i, j = 1, 2, 3). F can be a SM fermion, chargino or neutralino; S can be a sfermion or a Higgs/Goldstone boson; V can be a γ, Z or W ± . It should be noted that electron-Higgs couplings are neglected.  Figure 2: Generic tree, self-energy, vertex, box and counterterm diagrams for the process e + e − → h i Z (i = 1, 2, 3). F can be a SM fermion, chargino or neutralino; S can be a sfermion or a Higgs/Goldstone boson; V can be a γ, Z or W ± . It should be noted that electron-Higgs couplings are neglected. They are obtained from the corresponding tree-level diagrams by attaching a photon to the electrons/positrons. The internal particles in the generically depicted diagrams in Figs. 1 -3 are labeled as follows: F can be a SM fermion f , charginoχ ± c or neutralinoχ 0 n ; S can be a sfermionf s or a Higgs (Goldstone) boson h i , H ± (G, G ± ); U denotes the ghosts u V ; V can be a photon γ or a massive SM gauge boson, Z or W ± . We have neglected all electron-Higgs couplings via the FeynArts command [55] Restrictions -> NoElectronHCoupling and terms proportional to the electron mass ME (and the squared electron mass ME2) via the FormCalc command [56] Neglect[ME] = Neglect[ME2] = 0 which allows FormCalc to replace ME by zero whenever this is safe, i.e. except when it appears in negative powers or in loop integrals. We have verified numerically that these contributions are indeed totally negligible. For internally appearing Higgs bosons no higherorder corrections to their masses or couplings are taken into account; these corrections would correspond to effects beyond one-loop order. 2 For external Higgs bosons, as discussed in Ref. [32], the appropriateẐ factors are applied and on-shell (OS) masses (including higherorder corrections) are used [32], obtained with FeynHiggs [32][33][34][35][36].
Also not shown are the diagrams with a Z/Goldstone-Higgs boson self-energy contribution on the external Higgs boson leg. They appear in e + e − → h i h j , with a Z/G-h i,j transition and have been calculated explicitly as far as they are not proportional to the electron mass. It should be noted that for the process e + e − → h i Z all these contributions are proportional to the electron mass and have consistently be neglected.
Furthermore, in general, in Figs. 1 -3 we have omitted diagrams with self-energy type corrections of external (on-shell) particles. While the contributions from the real parts of the loop functions are taken into account via the renormalization constants defined by OS renormalization conditions, the contributions coming from the imaginary part of the loop functions can result in an additional (real) correction if multiplied by complex parameters. In the analytical and numerical evaluation, these diagrams have been taken into account via the prescription described in Ref. [47].
Within our one-loop calculation we neglect finite width effects that can help to cure threshold singularities. Consequently, in the close vicinity of those thresholds our calculation does not give a reliable result. Switching to a complex mass scheme [57] would be another possibility to cure this problem, but its application is beyond the scope of our paper.

Ultraviolet divergences
As regularization scheme for the UV divergences we have used constrained differential renormalization [58], which has been shown to be equivalent to dimensional reduction [59] at the one-loop level [56]. Thus the employed regularization scheme preserves SUSY [60,61] and guarantees that the SUSY relations are kept intact, e.g. that the gauge couplings of the SM vertices and the Yukawa couplings of the corresponding SUSY vertices also coincide to one-loop order in the SUSY limit. Therefore no additional shifts, which might occur when using a different regularization scheme, arise. All UV divergences cancel in the final result. 3

Infrared divergences
Soft photon emission implies numerical problems in the phase space integration of radiative processes. The phase space integral diverges in the soft energy region where the photon momentum becomes very small, leading to infrared (IR) singularities. Therefore the IR divergences from diagrams with an internal photon have to cancel with the ones from the corresponding real soft radiation. We have included the soft photon contribution via the code already implemented in FormCalc following the description given in Ref. [62]. The IR divergences arising from the diagrams involving a photon are regularized by introducing a photon mass parameter, λ. All IR divergences, i.e. all divergences in the limit λ → 0, cancel once virtual and real diagrams for one process are added. We have (numerically) checked that our results do not depend on λ.
We have also numerically checked that our results do not depend on ∆E = δ s E = δ s √ s/2 defining the energy cut that separates the soft from the hard radiation. As one can see from the example in the upper plot of Fig. 4 this holds for several orders of magnitude. Our numerical results below have been obtained for fixed δ s = 10 −3 .

Collinear divergences
Numerical problems in the phase space integration of the radiative process arise also through collinear photon emission. Mass singularities emerge as a consequence of the collinear photon emission off massless particles. But already very light particles (such as e.g. electrons) can produce numerical instabilities. There are several methods for the treatment of collinear singularities. In the following, we give a very brief description of the so-called phase space slicing (PSS) method [63], which we adopted. The treatment of collinear divergences is not (yet) implemented in FormCalc, and therefore we have developed and implemented the code necessary for the evaluation of collinear contributions.
In the PSS method, the phase space is divided into regions where the integrand is finite (numerically stable) and regions where it is divergent (or numerically unstable). In the stable regions the integration is performed numerically, whereas in the unstable regions it is carried out (semi-) analytically using approximations for the collinear photon emission.
The collinear part is constrained by the angular cut-off parameter ∆θ, imposed on the angle between the photon and the (in our case initial state) electron/positron.
The differential cross section for the collinear photon radiation off the initial state e + e − pair corresponds to a convolution with P ee (z) = (1+z 2 )/(1−z) denoting the splitting function of a photon from the initial e + e − pair. The electron momentum is reduced (because of the radiated photon) by the fraction z such that the center-of-mass frame of the hard process receives a boost. The integration over the divergence, the numerical difference between the two calculations was found to be negligible. Therefore we used the (faster) simplified code with neglected electron mass for our numerical analyses below. all possible factors z is constrained by the soft cut-off δ s = ∆E/E, to prevent over-counting in the soft energy region. We have (numerically) checked that our results do not depend on ∆θ over several orders of magnitude; see the example in the lower plot of Fig. 4. Our numerical results below have been obtained for fixed ∆θ/rad = 10 −2 .
The one-loop corrections of the differential cross section are decomposed into the virtual, soft, hard and collinear parts as follows: The hard and collinear parts have been calculated via the Monte Carlo integration algorithm Vegas as implemented in the CUBA library [64] as part of FormCalc.

Comparisons
In this section we present the comparisons with results from other groups in the literature for neutral Higgs boson production in e + e − collisions. Most of these comparisons were restricted to the MSSM with real parameters. The level of agreement of such comparisons (at one-loop order) depends on the correct transformation of the input parameters from our renormalization scheme into the schemes used in the respective literature, as well as on the differences in the employed renormalization schemes as such. In view of the non-trivial conversions and the large number of comparisons such transformations and/or change of our renormalization prescription is beyond the scope of our paper.
• A numerical comparison with the program FeynHiggsXS [41] can be found in Tab. 1.
We have neglected the initial state radiation and diagrams with photon exchange, as done in Ref. [41]. In Tab. 1 "self", "self+vert" and "full" denotes the inclusion of only self-energy corrections, self-energy plus vertex corrections or the full calculation including box diagrams. The comparison for the production of the light Higgs boson is rather difficult, due to the different FeynHiggs versions. As input parameters we used our scenario S; see Tab. 2 below. (We had to change only A t,b,τ to A t,b = 1500 and A τ = 0 to be in accordance with the input options of FeynHiggsXS). It can be observed that the level of agreement for the "self+vert" calculation is mostly at the level of 5% or better. However, the box contributions appear to go in the opposite direction for the first three cross sections in the two calculations. This hints towards a problem in the box contributions in Ref. [41], where the box contributions were obtained independently from the rest of the loop corrections, whereas using FeynTools all corrections are evaluated together in an automated way. It should be noted that a self-consistent check with the program FeynHiggsXS gave good agreement with Ref. [41] as expected (with tiny differences due to slightly different SM input parameters).
• In Ref. [65] the processes e + e − → HA, hA (and e + e − → H + H − ) have been calculated in the rMSSM. Unfortunately, in Ref. [65] the numerical evaluation (shown in their • In Ref. [66] the processes e + e − → AH, HZ, Ah have been calculated in the rMSSM at tree-level. As input parameters we used their parameters as far as possible. For the comparison with Ref. [66] we successfully reproduced their upper Fig. 2. • In Ref. [42] a tree-level evaluation of the channels (1) and (2) in the cMSSM was presented, where higher-order corrections were included via (CP violating) effective couplings. Unfortunately, no numbers are given in Ref. [42], but only two-dimensional parameter scan plots, which we could not reasonably compare to our results. Consequently we omitted a comparison with Ref. [42].
• We performed a comparison with Ref. [22] for in the cMSSM. In Ref. [22] only self-energy and vertex corrections involving t,t, b,b were included, and the numerical evaluation was performed in the CPX scenario [67] (with M H ± chosen to yield m h 1 = 40 GeV) which is extremely sensitive to the chosen input parameters. Nevertheless, using their input parameters as far as possible, we found qualitative agreement for t β < 15 with their Fig. 20.
• In Refs. [68,69] "supersimple" expressions have been derived for the processes e + e − → hZ, hγ in the rMSSM. We successfully reproduced Fig. 4 (right panels) of Ref. [68] in the upper plots of our Fig. 5 and Fig. 5 (right panels) of Ref. [69] in the lower plots of our Fig. 5. As input parameters we used their (SUSY) parameter set S1. The small differences in the differential cross sections are caused by the SM input parameters (where we have used our parameters; see Sect. 5.1 below) and the slightly different renormalization schemes and treatment of the Higgs boson masses.
• e + e − → hZ at the full one-loop level (including hard and soft photon bremsstrahlung) in the rMSSM has been analyzed in Ref. [43]. They also used FeynArts, FormCalc and LoopTools to generate and simplify their code. Unfortunately no numbers are given in Ref. [43], but only two-dimensional parameter scan plots, which we could not reasonably compare to our results. Consequently, we omitted a comparison with Ref. [43].
• The Higgsstrahlung process e + e − → hZ with the expected leading corrections in "Natural SUSY" models (i.e. a one-loop calculation with third-generation (s)quarks) has been computed in Ref. [70]. They also used FeynArts, FormCalc and LoopTools to generate and simplify their code. Unfortunately, again no numbers are given in Ref. [70], but mostly two-dimensional parameter scan plots, which we could not reasonably compare to our results. Only in the left plot of their Fig. 4 they show (fractional) corrections to the Higgsstrahlung cross section. However, the MSSM input parameters are not given in detail, rendering a comparison again impossible.
• In Ref. [71] the processes e + e − → hZ, hA are computed within a complete one-loop calculation. Only the QED (including photon bremsstrahlung) has been neglected. We  Figure 5: σ(e + e − → hZ, hγ). The left (right) plots show the differential cross section with √ s = 1 TeV (5 TeV) and cos ϑ varied. Upper row: Tree-level and one-loop corrected differential cross sections (in fb) are shown with parameters chosen according to S1 of Ref. [68]. Lower row: Loop induced differential cross sections (in ab) are shown with parameters chosen according to S1 of Ref. [69].
used their input parameters as far as possible and (more or less successfully) reproduced Fig. 5 and Fig. 6 (upper rows, solid lines) of Ref. [71] qualitatively in our Fig. 6. The differences are mainly due to different Higgs boson masses and the use of Higgs boson wave function corrections in Ref. [71], while we used an effective mixing angle α eff . In order to facilitate the comparison we used the same simple formulas for our Higgs boson masses and α eff as in their Eqs. (4)- (7). Therefore our σ tree correspond rather to their σ ǫ and our σ full rather to their σ FDC . It should be noted that the code of Ref. [71] is also part of the code from Ref. [41]. Using FeynHiggsXS with the input parameters of Ref. [71] (as far as possible) gave also only qualitative agreement with the Figs. of Ref. [71].
• In Ref. [72] the loop induced processes e + e − → hγ, Hγ, Aγ have been computed. We used the same simple formulas for our Higgs boson masses and α eff as in their Eqs.  (3.48)-(3.50). We also used their input parameters as far as possible, but unfortunately they forgot to specify the trilinear parameters A f . Therefore we chose arbitrarily A f = 0 for our comparison. In view of this problem the comparison is acceptable; see our Fig. 7 vs. Fig. 4, Fig. 5 and Fig. 7 of Ref. [72]. It should be noted that the code of Ref. [72] is also part of the code from Ref. [41] and Ref. [71].
A final comment is in order. We argue that the problems in the comparison with Ref. [41] (i.e. FeynHiggsXS), Ref. [71] and Ref. [72] are due to the fact that all three papers are based (effectively) on the same calculation/source. Therefore, these three papers should be considered as one rather than three independent comparisons, and thus do not disprove the reliability of our calculation. It should also be kept in mind that our calculational method/code has already been successfully tested and compared with quite a few other programs, see Refs. [23,24,[47][48][49][50][51][52][53].

Numerical analysis
In this section we present our numerical analysis of neutral Higgs boson production at e + e − colliders in the cMSSM. In the various figures below we show the cross sections at the treelevel ("tree") and at the full one-loop level ("full"). In case of extremely small tree-level cross sections we also show results including the corresponding purely loop induced contributions ("loop"). These leading two-loop contributions are ∝ |M 1-loop | 2 , where M 1-loop denotes the one-loop matrix element of the appropriate process.
The SUSY parameters are chosen according to the scenario S, shown in Tab. 2, unless otherwise noted. This scenario constitutes a viable scenario for the various cMSSM Higgs production modes, i.e. not picking specific parameters for each cross section. The only variation will be the choice of √ s = 500 GeV for production cross sections involving the  4 This will be clearly indicated below. We do not strictly demand that the lightest Higgs boson has a mass around ∼ 125 GeV, although for most of the parameter space this is given. We will show the variation with √ s, M H ± , t β and ϕ At , the phase of A t .
Concerning the complex parameters, some more comments are in order. No complex parameter enters into the tree-level production cross sections. Therefore, the largest effects are expected from the complex phases entering via the t/t sector, i.e. from ϕ At , motivating our choice of ϕ At as parameter to be varied. Here the following should be kept in mind. When performing an analysis involving complex parameters it should be noted that the results for physical observables are affected only by certain combinations of the complex phases of the parameters µ, the trilinear couplings A f and the gaugino mass parameters M 1,2,3 [76,77]. It is possible, for instance, to rotate the phase ϕ M 2 away. Experimental constraints on the (combinations of) complex phases arise, in particular, from their contributions to electric dipole moments of the electron and the neutron (see Refs. [78,79] and references therein), of the deuteron [80] and of heavy quarks [81]. While SM contributions enter only at the three-loop level, due to its complex phases the MSSM can contribute already at one-loop order. Large phases in the first two generations of sfermions can only be accommodated if these generations are assumed to be very heavy [82] or large cancellations occur [83]; see, however, the discussion in Ref. [84]. A review can be found in Ref. [85]. Accordingly (using the convention that ϕ M 2 = 0, as done in this paper), in particular, the phase ϕ µ is tightly constrained [86], while the bounds on the phases of the third generation trilinear couplings are much weaker. Setting ϕ µ = 0 and ϕ A f =t = 0 leaves us with ϕ At as the only complex valued parameter.
Since now the complex trilinear coupling A t can appear in the couplings, contributions from absorptive parts of self-energy type corrections on external legs can arise. The corresponding formulas for an inclusion of these absorptive contributions via finite wave function correction factors can be found in Refs. [47,49].
The numerical results shown in the next subsections are of course dependent on the choice of the SUSY parameters. Nevertheless, they give an idea of the relevance of the full one-loop corrections.

√ s, M H ± , t β , and ϕ A t
The results shown in this and the following subsections consist of "tree", which denotes the tree-level value and of "full", which is the cross section including all one-loop corrections as described in Sect. 3. We begin the numerical analysis with the cross sections of e + e − → h i h j (i, j = 1, 2, 3) evaluated as a function of We begin with the process e + e − → h 1 h 2 as shown in Fig. 8. As a general comment it should be noted that in S one finds that h 1 ∼ h, h 2 ∼ A and h 3 ∼ H. The hAZ coupling is ∝ c β−α which goes to zero in the decoupling limit [88], and consequently relatively small cross sections are found. In the analysis of the production cross section as a function of √ s (upper left plot) we find the expected behavior: a strong rise close to the production threshold, followed by a decrease with increasing √ s. We find a relative correction of ∼ −15% around the production threshold. Away from the production threshold, loop corrections of ∼ +27% at √ s = 1000 GeV are found in S (see Tab. 2). The relative size of loop corrections decrease with increasing √ s and reach ∼ +61% at √ s = 3000 GeV where the tree-level becomes very small. With increasing M H ± in S (upper right plot) we find a strong decrease of the production cross section, as can be expected from kinematics, but in particular from the decoupling limit discussed above. The loop corrections reach ∼ +27% at M H ± = 300 GeV and ∼ +62% at M H ± = 500 GeV. These large loop corrections are again due to the (relative) smallness of the tree-level results. It should be noted that at M H ± ≈ 350 GeV the limit of 0.01 fb is reached, corresponding to 10 events at an integrated luminosity of L = 1 ab −1 . The cross sections decrease with increasing t β (lower left plot), and the loop corrections reach the maximum of ∼ +38% at t β = 36 while the minimum of ∼ +26% is at t β = 5. The phase dependence ϕ At of the cross section in S (lower right plot) is at the 10% level at tree-level. The loop corrections are nearly constant, ∼ +28% for all ϕ At values and do not change the overall dependence of the cross section on the complex phase.
Not shown is the process e + e − → h 1 h 3 . In this case, for our parameter set S (see Tab. 2), one finds h 3 ∼ H. Due to the absence of the hHZ coupling in the MSSM (see Ref. [87]) this leads to vanishing tree-level cross sections in the case of real parameters. For complex  parameters (i.e. ϕ At ) the tree-level results stay below 10 −5 fb. Also the loop induced cross sections ∝ |M 1-loop | 2 (where only the vertex and box diagrams contribute in the case of real parameters) stay below 10 −5 fb for our parameter set S. Consequently, in this case we omit showing plots to the process e + e − → h 1 h 3 .
In Fig. 9 we present the cross section e + e − → h 2 h 3 with h 2 ∼ A and h 3 ∼ H in S. The HAZ coupling is ∝ s β−α , which goes to one in the decoupling limit, and consequently relatively large cross sections are found. In the analysis as a function of √ s (upper left plot) we find relative corrections of ∼ −37% around the production threshold, ∼ −5% at √ s = 1000 GeV (i.e. S), and ∼ +6% at √ s = 3000 GeV. The dependence on M H ± (upper right plot) is nearly linear above M H ± > ∼ 250 GeV, and mostly due to kinematics. The loop corrections are ∼ −8% at M H ± = 160 GeV, ∼ −5% at M H ± = 300 GeV (i.e. S), and ∼ −52% at M H ± = 500 GeV where the tree-level goes to zero. As a function of t β (lower left plot) the tree-level cross section is rather flat, apart from a dip at t β ≈ 10, corresponding to the threshold mχ± corrections only via theẐ matrix contribution (calculated by FeynHiggs). The relative corrections increase from ∼ −5% at t β = 7 to ∼ +7% at t β = 50. The dependence on ϕ At (lower right plot) is very small, below the percent level. The loop corrections are found to be nearly independent of ϕ At at the level of ∼ −4.6%.
We now turn to the processes with equal indices. The tree couplings h i h i Z (i = 1, 2, 3) are exactly zero; see Ref. [87]. Therefore, in this case we show the pure loop induced cross sections ∝ |M 1-loop | 2 (labeled as "loop") where only the box diagrams contribute. These box diagrams are UV and IR finite.
In Fig. 10 we show the results for e + e − → h 1 h 1 . This process might have some special interest, since it is the lowest energy process in which triple Higgs boson couplings play a role, which could be relevant at a high-luminosity collider operating above the two Higgs boson production threshold. In our numerical analysis, as a function of √ s we find a maximum of ∼ 0.014 fb, at √ s = 500 GeV, decreasing to ∼ 0.002 fb at √ s = 3 TeV. The dependence on M H ± is rather small, as is the dependence on t β and ϕ At in S. However, with cross sections  found at the level of up to 0.015 fb this process could potentially be observable at the ILC running at √ s = 500 GeV or below (depending on the integrated luminosity).
We finish the e + e − → h i h i analysis in Fig. 11 in which the results for i = 2, 3 are displayed. Both production processes have rather similar (purely loop-induced) production cross sections. As a function of √ s we find a maximum of ∼ 0.0035 fb at √ s = 1.4 TeV. In S, but with M H ± varied we find the highest values of ∼ 0.007 fb at the lowest mass scales, going down below 0.001 fb at around M H ± ∼ 380 GeV. The production cross sections depend only very weakly on t β and ϕ At , where in S values of ∼ 0.0026 fb are found, leading only to about 5 events for an integrated luminosity of L = 2 ab −1 . Furthermore, due to the similar decay patterns of h 2 ∼ A and h 3 ∼ H and the similar masses of the two states it might be difficult to disentangle it from e + e − → h 2 h 3 , and a more dedicated analysis (beyond the scope of our paper) will be necessary to determine its observability. The large dip at t β ≈ 10 (red solid line) is the threshold mχ± 1 + mχ± 1 = m h 2 in e + e − → h 2 h 2 . The dip at t β ≈ 5 (blue dotted  Overall, for the neutral Higgs boson pair production we observed an increasing cross section ∝ 1/s for s → ∞; see Eq. (4). The full one-loop corrections reach a level of 10% -20% or higher for cross sections of 0.01 -10 fb. The variation with ϕ At is found rather small, except for e + e − → h 1 h 2 , where it is at the level of 10%. The results for h i h i production turn out to be small (but not necessarily hopeless) for i = 1, and negligible for i = 2, 3 for Higgs boson masses above ∼ 200 GeV. (see [87]). In the case of real parameters this leads to vanishing tree-level cross sections if h i ∼ A.
We start with the process e + e − → h 1 Z shown in Fig. 12. In S one finds h 1 ∼ h, and since the ZZh coupling is ∝ s β−α → 1 in the decoupling limit, relative large cross sections are found. As a function of √ s (upper left plot) a maxium of more than 200 fb is found at √ s ∼ 250 GeV with a decrease for increasing √ s. The size of the corrections of the cross section can be especially large very close to the production threshold 5 from which on the considered process is kinematically possible. At the production threshold we found relative corrections of ∼ −60%. Away from the production threshold, loop corrections of ∼ +20% at √ s = 500 GeV are found, increasing to ∼ +30% at √ s = 3000 GeV. In the following plots we assume, deviating from the definition of S, √ s = 500 GeV. As a function of M H ± (upper right plot) the cross sections strongly increases up to M H ± < ∼ 250 GeV, corresponding to s β−α → 1 in the decoupling limit discussed above. For higher M H ± values it is nearly constant, and the loop corrections are ∼ +20% for 160 GeV < M H ± < 500 GeV. Hardly any variation is found for the production cross section as a function of t β or ϕ At . In both cases the one-loop corrections are found at the level of ∼ +20%.
Not shown is the process e + e − → h 2 Z. In this case, for our parameter set S (see Tab. 2), one finds h 2 ∼ A. Because there are no AZZ couplings in the MSSM (see [87]) this leads to vanishing tree-level cross sections in the case of real parameters. For complex parameters (i.e. ϕ At ) the tree-level results stay below 10 −5 fb. Also the loop induced cross sections ∝ |M 1-loop | 2 (where only the vertex and box diagrams contribute in the case of real parameters) are below 10 −3 fb for our parameter set S. Consequently, in this case we omit showing plots to the process e + e − → h 2 Z.
We finish the e + e − → h i Z analysis in Fig. 13 in which the results for e + e − → h 3 Z are shown. In S one has h 3 ∼ H, and with the ZZH coupling being proportional to c β−α → 0 in the decoupling limit relatively small production cross sections are found for M H ± not too small. As a function of √ s (upper left plot) a dip can be seen at √ s ≈ 540 GeV, due to the threshold mχ± 2 + mχ± 2 = √ s. Around the production threshold we found relative corrections of ∼ 3%. The maxium production cross section is found at √ s ∼ 500 GeV of about 0.065 fb including loop corrections, rendering this process observable with an accumulated luminosity L < ∼ 1 ab −1 . Away from the production threshold, one-loop corrections of ∼ 47% at √ s = 1000 GeV are found in S (see Tab. 2), with a cross section of about 0.03 fb. The cross section further decreases with increasing √ s and the loop corrections reach ∼ 45% at √ s = 3000 GeV, where it drops below the level of 0.0025 fb. As a function of M H ± we find the afore mentioned decoupling behavior with increasing M H ± . The the loop corrections reach ∼ 26% at M H ± = 160 GeV, ∼ 47% at M H ± = 300 GeV and ∼ +56% at M H ± = 500 GeV. These large loop corrections (> 50%) are again due to the (relative) smallness of the treelevel results. It should be noted that at M H ± ≈ 360 GeV the limit of 0.01 fb is reached; see the line in the upper right plot. The production cross section decreases strongly with t β (lower right plot). The loop corrections reach the maximum of ∼ +95% at t β = 50 do to the very small tree-level result, while the minimum of ∼ +47% is found at t β = 7. The phase dependence ϕ At of the cross section (lower right plot) is at the level of 5% at tree-level, but increases to about 10% including loop corrections. Those are found to vary from ∼ +47% at ϕ At = 0 • , 360 • to ∼ +39% at ϕ At = 180 • .
Overall, for the Z Higgs boson production we observed an increasing cross section ∝ 1/s for s → ∞; see Eq. (5). The full one-loop corrections reach a level of 20% (50%) for cross sections of 60 fb (0.03 fb). The variation with ϕ At is found to be small, reaching up to 10% for e + e − → h 3 Z, after including the loop corrections.  couplings in the MSSM; see Ref. [87]. In the following analysis e + e − → h i γ (i = 1, 2, 3) are purely loop induced processes (via vertex and box diagrams) and therefore ∝ |M 1-loop | 2 .
We start with the process e + e − → h 1 γ shown in Fig. 14  the production cross sections from 0.023 fb at M H ± ≈ 160 GeV to about 0.03 fb in the decoupling regime. This dependence shows the relevance of the SM gauge boson loops in the production cross section, indicating that the top quark loops dominate this production cross section. The variation with t β and ϕ At (lower row) is rather small, and values of 0.03 fb are found in S.
We finish the e + e − → h i γ analysis in Fig. 15   is given roughly by the sum of the two. This renders these loop-induced processes at the border of observability. The peaks oberved are found at √ s ≈ 540 GeV due to the threshold mχ± 2 + mχ± 2 = √ s for both production cross sections. They drop to the unobservable level for √ s > ∼ 1 TeV. As a function of M H ± (upper right plot) one can observe the decoupling of h 3 ∼ H of the SM gauge bosons with increasing M H ± , lowering the cross section for larger values. The "knee" at M H ± ≈ 294 GeV is the threshold mχ± 1 + mχ± 1 = m h 2 . This threshold enter into the loop corrections only via theẐ matrix contribution (calculated by FeynHiggs). The loop corrections vary between 0.008 fb at M H ± ≈ 160 GeV and far below 0.001 fb at M H ± ≈ 500 GeV. The dependence on t β (lower left plot) is rather strong for the h 2 γ production going from 0.007 fb at t β = 4 down to 0.0035 fb at t β = 50. The dip at t β ≈ 10 is the threshold mχ± 1 + mχ± 1 = m h 2 . This threshold enter into the loop corrections again only via theẐ matrix contribution (calculated by FeynHiggs). For the h 3 γ production the cross section stays at the very low level of 0.001 fb for all t β values. The dependence on the phase ϕ At of the cross sections (lower right plot) is very small in S, with no visible variation in the plot.
Overall, for the γ Higgs boson production the leading order corrections can reach a level of 0.1 fb, depending on the SUSY parameters. This renders these loop-induced processes in principle observable at an e + e − collider. The variation with ϕ At is found to be extremely small.

Conclusions
We evaluated all neutral MSSM Higgs boson production modes at e + e − colliders with a twoparticle final state, i.e. e + e − → h i h j , h i Z, h i γ (i, j = 1, 2, 3), allowing for complex parameters. In the case of a discovery of additional Higgs bosons a subsequent precision measurement of their properties will be crucial determine their nature and the underlying (SUSY) parameters. In order to yield a sufficient accuracy, one-loop corrections to the various Higgs boson production modes have to be considered. This is particularly the case for the high anticipated accuracy of the Higgs boson property determination at e + e − colliders [19].
The evaluation of the processes (1) -(3) is based on a full one-loop calculation, also including hard and soft QED radiation. The renormalization is chosen to be identical as for the various Higgs boson decay calculations; see, e.g., Refs. [23,24].
We first very briefly reviewed the relevant sectors including some details on the one-loop renormalization procedure of the cMSSM, which are relevant for our calculation. In most cases we follow Ref. [47]. We have discussed the calculation of the one-loop diagrams, the treatment of UV, IR and collinear divergences that are canceled by the inclusion of (hard, soft and collinear) QED radiation. We have checked our result against the literature as far as possible, and in most cases we found accaptable or qualitative agreement, where parts of the differences can be attributed to problems with input parameters (conversions) and/or special szenarios. Once our set-up was changed successfully to the one used in the existing analyses we found good agreement.
For the analysis we have chosen a parameter set that allows simultaneously a maximum number of production processes. In this scenario (see Tab. 2) we have h 1 ∼ h, h 2 ∼ A and h 3 ∼ H. In the analysis we investigated the variation of the various production cross sections with the center-of-mass energy √ s, the charged Higgs boson mass M H ± , the ratio of the vacuum expectation values t β and the phase of the trilinear Higgs-top squark coupling, ϕ At . For light (heavy) Higgs production cross sections we have chosen √ s = 500 (1000) GeV.
In our numerical scenarios we compared the tree-level production cross sections with the full one-loop corrected cross sections. In certain cases the tree-level cross sections are identical zero (due to the symmetries of the model), and in those cases we have evaluated the one-loop squared amplitude, σ loop ∝ |M 1-loop | 2 .
We found sizable corrections of ∼ 10 − 20% in the h i h j production cross sections. Substantially larger corrections are found in cases where the tree-level result is (accidentally) small and thus the production mode likely is not observable. The purely loop-induced processes of e + e − → h i h i could be observable, in particular in the case of h 1 h 1 production. For the h i Z modes corrections around 10 − 20%, but going up to ∼ 50%, are found. The purely loop-induced processes of h i γ production appear observable for h 1 γ, but very challenging for h 2,3 γ.
Only in very few cases a relevant dependence on ϕ At was found. Examples are e + e − → h 1 h 2 and e + e − → h 3 Z, where a variation, after the inclusion of the loop corrections, of up to 10% with ϕ At was found. In those cases neglecting the phase dependence could lead to a wrong impression of the relative size of the various cross sections.
The numerical results we have shown are, of course, dependent on the choice of the SUSY parameters. Nevertheless, they give an idea of the relevance of the full one-loop corrections. Following our analysis it is evident that the full one-loop corrections are mandatory for a precise prediction of the various cMSSM Higgs boson production processes. The full one-loop corrections must be taken into account in any precise determination of (SUSY) parameters from the production of cMSSM Higgs bosons at e + e − linear colliders. There are plans to implement the evaluation of the Higgs boson production into the public code FeynHiggs.