text:
# Electron energy losses near pulsar polar caps: a Monte Carlo approach

## 1 Introduction

Curvature radiation (CR) and magnetic inverse Compton scattering (ICS) are usually considered to be the most natural ways of hard gamma-rays production operating at the expense of pulsar rotational energy (eg. Zhang & Harding zh (2000), and references therein). These two radiation mechanisms dominate within two different ranges of Lorentz factors $`\gamma `$ of beam particles (ie. those leaving the polar cap). When $`\gamma <10^6`$, magnetic inverse Compton scattering plays a dominant role in braking beam particles (Xia et al. xqwh (1985); Chang chang (1995); Sturner sturner (1995)) and is the main source of hard gamma-ray photons (Sturner & Dermer sd (1994); Sturner et al. sdm (1995)).

When $`\gamma >10^6`$, the curvature radiation becomes responsible for cooling beam particles (eg. Daugherty & Harding 1982), yet, the inverse Compton scattering may still be important for secondary $`e^\pm `$-pair plasma. Zhang & Harding (zh (2000)) show that different dependencies of X-ray and gamma-ray luminosities on the pulsar spin-down luminosity as inferred from observations by CGRO, ROSAT and ASCA (Arons arons (1996); Becker & Trümper bt (1997); Saito et al. skks (1998), respectively) can be well reproduced when this effect is included.

So far, the influence of ICS on electron energy have been investigated most thoroughly by Chang (chang (1995)), Sturner (sturner (1995)), Harding & Muslimov (hm (1998)), and Supper & Trümper (supper (2000)). One of several interesting results they found was that energy losses due to resonant ICS can limit the Lorentz factors $`\gamma `$ of electrons to a value which depends on both electric field strength $`E_{}`$, temperature $`T`$ and radius $`R_{\mathrm{th}}`$ of thermal polar cap, and on the magnetic field strength $`B_{\mathrm{pc}}`$ at a polar cap. For example Sturner (sturner (1995)) shows that assuming the acceleration model of Michel (michel (1974)) and a thermal polar cap size comparable to that determined by the open field lines, the Lorentz factors $`\gamma `$ are limited to $`10^3`$ if $`B_{\mathrm{pc}}>10^{13}`$ G, and $`T>3\times 10^6`$ K. This acceleration stopping effect becomes more efficient for stronger magnetic fields, and it led Sturner (sturner (1995)) to propose it as a possible explanation for an apparent cutoff around 10 MeV in gamma radiation of B1509$-$58 (Kuiper et al. kuiper (1999)).

Those results were obtained with a numerical method which assumes smooth changes in electron energy and makes use of expressions for *averaged* energy losses due to the magnetic ICS (Dermer dermer (1990)). The continuous changes of the Lorentz factor with height as determined with this fully deterministic treatment are considered as representative for a behaviour of a large number of electrons. Actually however, electrons lose their energy in discontinuous scattering events. Between some of these events (which *occasionally* can be very distant even for a short mean free path) the energy of an accelerated electron may increase considerably (even a few times) and therefore, the electron can find itself in quite different conditions than if it were losing its energy continuously (the mean free path depends very strongly on $`\gamma `$). Since the energy loss rate due to the resonant ICS is not a monotonic function of $`\gamma `$ this effect may be essential for future behaviour of the particle.

The aim of our paper is to investigate the energetics of electrons on individual basis with a Monte Carlo treatment. In particular, a question of the ICS-limited acceleration is adressed. In Section 2 we calculate electron energy losses near a neutron star using both, the Monte Carlo approach and then the smooth method (after Dermer 1990; Chang 1995; Sturner 1995) for comparison. We consider two distinct cases of electron acceleration: 1) in the constant electric field, proposed by Michel (michel (1974)) (hereafter M74); 2) in the electric field elaborated by Harding & Muslimov (hm (1998)) (hereafter HM98). Section 3 presents the results for electron energetics, and the finding that within some range of pulsars parameters, the averaged Monte Carlo behaviour of electrons does not coincide with a solution found with the smooth method. It also contains the explanation of these differences as well as consequences for predicted efficiencies of gamma-ray emission. Conclusions are given in Sect. 4.

## 2 The model

We consider electrons accelerated in a longitudinal electric field induced rotationally within the region adjacent to the surface of a neutron star. The particles lose their energy due to scattering off soft thermal X-ray photons through magnetic inverse Compton mechanism and (marginally) due to emission of curvature radiation. The thermal photons are assumed to originate from a flat ‘thermal polar cap’ with a temperature $`T\mathrm{a}\mathrm{few}\times 10^6`$ K and with a radius $`R_{\mathrm{th}}\mathrm{a}\mathrm{few}\times r_{\mathrm{pc}}`$ where $`r_{\mathrm{pc}}=(2\pi R_{\mathrm{ns}}^3/cP)^{1/2}`$ is the standard polar cap radius as determined for an aligned rotator by the open lines of purely dipolar magnetic field; ($`R_{\mathrm{ns}}`$ denotes the neutron star radius and $`P`$ is the pulsar period).

Since the ICS losses are to dominate over the curvature losses, one has to ensure that electrons do not attain extremely high Lorentz factors ($`\gamma <10^6`$). Therefore, either the electric field should be weak or at least the size of the accelerating region should be small enough to prevent $`\gamma >10^6`$. For this reason most papers focusing on the resonant ICS made use of a relatively weak electric field after Michel michel (1974) (for example: Sturner sturner (1995) and recently Supper & Trümper supper (2000)). Michel’s model is considered as unrealistic since it ignores magnetic field line curvature (detailed treatment of this problem was first presented by Arons & Scharlemann as (1979)) as well as inertial frame dragging (its importance was acknowledged and the effect was worked out for the first time by Muslimov & Tsygan mt (1992)). Both effects lead to significantly stronger electric fields operating in polar-cap accelerators.

However, for the sake of comparison with previous papers on the resonant ICS we first consider the weak electric field after Michel michel (1974). Accordingly, we assume that electrons are accelerated in a constant $`E_{}`$ extending between stellar surface ($`h=0`$) and altitude $`h_{\mathrm{acc}}`$ :

$$E_{}=\{\begin{array}{c}\left(\frac{8\pi mcB_{\mathrm{pc}}}{eP}\right)^{1/2}\hfill \\ \text{ }3.5810^5\left(\frac{B_{12}}{P}\right)^{1/2}\mathrm{V}\mathrm{cm}^1,\hfill \\ \text{ }\mathrm{for}0hh_{\mathrm{acc}};\hfill \\ \text{ }0,\text{ }\mathrm{for}h>h_{\mathrm{acc}}.\hfill \end{array}$$
(1)
with $`h_{\mathrm{acc}}=r_{\mathrm{pc}}`$, and where $`P`$ is in seconds, $`B_{12}`$ is $`B_{\mathrm{pc}}`$ in Teragauss. Fixed location of $`h_{\mathrm{acc}}`$ means that the model is not self-consistent in the sense that the accelerating field will not be shorted out by a pair formation front (should it occur at a height below $`h_{\mathrm{acc}}=r_{\mathrm{pc}}`$).

As the second case for $`E_{}`$ we take the advanced model of an electric field in the form elaborated by Harding & Muslimov (hm (1998)). HM98 were the first to incorporate electron energy losses due to the ICS in a strong-$`E_{}`$ acceleration model. However, this was done to introduce a high altitude accelerator with negligible resonant ICS losses (an idea of an unstable ICS-induced pair formation fronts). In that work Harding and Muslimov also calculated self-consistent values of the acceleration height $`h_{\mathrm{acc}}`$ using an accelerating electric field which takes into account the upper pair formation front (eqs. (18) and (23) in HM98).

For $`h_{\mathrm{acc}}<0.1r_{\mathrm{pc}}`$ the field is well approximated with

$`E_{}(h)3{\displaystyle \frac{\mathrm{\Omega }R_{\mathrm{ns}}}{c}}{\displaystyle \frac{B_{\mathrm{pc}}}{1ϵ_{\mathrm{GR}}}}{\displaystyle \frac{h}{R_{\mathrm{ns}}}}{\displaystyle \frac{h_{\mathrm{acc}}}{R_{\mathrm{ns}}}}(1{\displaystyle \frac{h}{h_{\mathrm{acc}}}})\times `$ (2)
$`\text{ }\times \left[\kappa \mathrm{cos}\chi +A\mathrm{sin}\chi \right],`$
where $`\mathrm{\Omega }=2\pi /P`$, $`(1ϵ_{\mathrm{GR}})0.6`$, $`\kappa 0.15`$, $`\chi `$ is the angle between the magnetic dipole and rotation axes, $`h`$ is the altitude, and the factor $`A\mathrm{sin}\chi `$ is negligible for nonorthogonal rotators (see HM98 for details). For $`h_{\mathrm{acc}}`$ approaching $`r_{\mathrm{pc}}`$ the linear dependence on $`h_{\mathrm{acc}}`$ in Eq.(2) disappears. We have found that for $`0.5r_{\mathrm{pc}}<h_{\mathrm{acc}}<3r_{\mathrm{pc}}`$ and $`h<r_{\mathrm{pc}}/3`$ the electric field given by eq.(18) in HM98 is well aproximated with a formula

$`E_{}(h)3{\displaystyle \frac{\mathrm{\Omega }R_{\mathrm{ns}}}{c}}{\displaystyle \frac{B_{\mathrm{pc}}}{(1ϵ_{\mathrm{GR}})^{1/2}}}{\displaystyle \frac{h}{R_{\mathrm{ns}}}}{\displaystyle \frac{r_{\mathrm{pc}}}{R_{\mathrm{ns}}}}(1{\displaystyle \frac{h}{h_{\mathrm{acc}}}})\times `$ (3)
$`\text{ }\times \left[\kappa \mathrm{cos}\chi +A\mathrm{sin}\chi \right].`$
For purposes of comparison as well as to simplify our calculations, we inject electrons from the center of polar cap and propagate them along the straight magnetic axis field line. The propagation proceeds in steps of a size $`dh`$, which is one hundred times smaller than any of characteristic length scales involved in the problem, like the acceleration length scale or the mean free path for the magnetic ISC. In each of the distance steps the energy of an electron is increased by a value $`eE_{}dh`$. For the sake of completeness, we also include the energy loss due to the curvature radiation $`|\dot{\gamma }_{\mathrm{cr}}|mc^2dh/v`$ where $`\dot{\gamma }_{\mathrm{cr}}`$ is the CR cooling rate given by

$$\dot{\gamma }_{\mathrm{cr}}=\frac{2e^2\beta ^4\gamma ^4}{3mc\rho _{\mathrm{curv}}^2}$$
(4)
and $`v=c\beta `$ is the electron velocity. Nevertheless, the CR is a very inefficient cooling mechanism in the presence of acceleration models we consider below and it could be neglected equally well. When computing energy losses due to the CR we assume artificially (after Sturner sturner (1995)) that the magnetic axis has a fixed radius of curvature $`\rho _{\mathrm{curv}}=10^7`$ cm, though we keep the electrons all the time directly over the polar cap center.

Magnetic inverse Compton scattering has been treated with a Monte Carlo simulations. The framework of our numerical code is based on an approach proposed by Daugherty & Harding (dh89 (1989)). We improve their method of sampling the parameters of incoming photon and account for the Klein-Nishina regime in an approximate way.

At each step in the electron trajectory, the optical depth for magnetic ICS is calculated as $`d\tau =dh/c`$, where $``$ is a scattering rate, to decide whether a scattering event is to occur. The scattering rate $``$ in the observer frame OF is calculated as

$$=c𝑑\mathrm{\Omega }𝑑\sigma \left(\frac{dn_{\mathrm{ph}}}{dd\mathrm{\Omega }}\right)(1\beta \mu )$$
(5)
(eg. Ho & Epstein he (1989)) where $`\mathrm{\Omega }=d\mu d\varphi `$ is the solid angle subtended by the source of soft photons, ($`\mu =\mathrm{cos}\theta `$), $`\sigma `$ is a total cross section (see below), and $`dn_{\mathrm{ph}}/d/d\mathrm{\Omega }`$ is the density of the soft photons per unit energy and per unit solid angle. The symbol $``$ denotes photon energy in dimensional energy units. Hereafter we will use its dimensionless counterpart $`ϵ=/(mc^2)`$ to denote the photon energy in the observer frame OF and the primed symbol $`ϵ^{}=ϵ\gamma (1\beta \mu )`$ in the electron rest frame ERF.

For temperatures of the thermal polar cap and Lorentz factors of electrons considered below, photon energies $`ϵ^{}`$ may fall well above $`ϵ_B=B/B_{\mathrm{cr}}`$, a local magnetic field strength in units of the critical magnetic field $`B_{\mathrm{cr}}=m^2c^3(e\mathrm{})^1`$. This suggests a full form of the relativistic cross section for Compton scattering in strong magnetic fields to be used (Daugherty & Harding dh86 (1986)). However, incoming photons that propagate in the OF at an angle $`\theta =\mathrm{cos}^1\mu `$ are strongly collimated in the ERF with a cosine of polar angle $`\mu ^{}=(\mu \beta )/(1\beta \mu )`$ close to $`1`$. As was shown by Daugherty & Harding (dh86 (1986)), when $`|\mu ^{}|`$ approaches $`1`$, resonances at higher harmonics become narrower and weaker, and scattering into higher Landau states becomes less important. In such conditions, the polarization-averaged relativistic magnetic cross section in the Thomson regime is reasonably well approximated with a nonrelativistic, classical limit:

$$\sigma =\frac{\sigma _T}{2}\left(1\mu _{}^{}{}_{}{}^{2}+(1+\mu _{}^{}{}_{}{}^{2})\left[g_1+\frac{g_2g_1}{2}\right]\right)$$
(6)
where $`\sigma _T`$ is the Thomson cross section, and $`g_1`$ and $`g_2`$ are given by

$$g_1(u)=\frac{u^2}{(u+1)^2},\text{ }g_2(u)=\frac{u^2}{(u1)^2+a^2}$$
(7)
with $`uϵ^{}/ϵ_B`$, $`a2\alpha _fϵ_B/3`$, where $`\alpha _f`$ is a fine-structure constant (eg. Herold herold (1979); Dermer dermer (1990)).

In the Klein-Nishina regime ($`ϵ^{}>1`$) the relativistic magnetic cross section for the $`|\mu ^{}|1`$ case becomes better approximated with the well known Klein-Nishina relativistic nonmagnetic total cross section $`\sigma _{KN}`$ (Daugherty & Harding dh86 (1986); Dermer dermer (1990)). When $`ϵ^{}ϵ_B`$, (the condition fulfilled in the K-N regime since $`ϵ_B<1`$ holds throughout this paper), the resonant term $`(g_2g_1)/2`$ in Eq. (6) becomes negligible and the nonresonant term $`g_1`$ approaches unity which results in $`\sigma \sigma _T`$. Therefore, to approximate the Klein-Nishina decline we replace the single nonresonant term $`g_1`$ in the square bracket of Eq. (6) with $`\sigma _{KN}/\sigma _T`$ for $`ϵ^{}>2ϵ_B`$.

To make calculations less time-consuming, we have used the delta-function approximation for the resonant part of the cross section (Dermer dermer (1990)), when calculating the optical depth for a scattering and the averaged energy loss of electron (Eq. 11, see below).

As the density of soft photons $`dn_{\mathrm{ph}}/d/d\mathrm{\Omega }`$ we take the spectral density of blackbody radiation given by

$$\left(\frac{dn_{\mathrm{ph}}}{dd\mathrm{\Omega }}\right)dd\mathrm{\Omega }=\frac{2}{\lambda _C^3}\frac{ϵ^2dϵd\mathrm{\Omega }}{[\mathrm{exp}(ϵ/𝒯)1]}$$
(8)
where $`𝒯`$ is a dimensionless temperature ($`𝒯\frac{kT}{mc^2}`$) and $`\lambda _C=h/(mc)`$ is the electron Compton wavelength. Hereby we neglect anisotropy of the thermal emission expected in the strong magnetic field (eg Pavlov et al. psvz (1994)). For $`ϵϵ_B`$ a preferred photon propagation direction is that of the magnetic field because the opacity is reduced for radiation polarized across $`B`$. The probability of collision with an electron is greatly reduced for photons propagating at $`\mu 1`$, thus, the efficiency of ICS as calculated below should be treated as an upper limit.

If a scattering event occurs, we draw the energy $`ϵ`$ and the cosine of polar angle $`\mu `$ for an incoming photon from a distribution determined with the integrand of Eq. (5). We do this with the simple two-dimensional rejection method which is reliable though more time-consuming (cf. an approach by Daugherty & Harding dh89 (1989)). Next, we transform $`ϵ`$ and $`\mu `$ to the ERF values $`ϵ^{}`$ and $`\mu ^{}`$, and sample the cosine of polar angle of an outgoing photon $`\mu _s^{}`$ using the differential form of the cross section (6) in the Thomson limit:

$`{\displaystyle \frac{d\sigma }{dϵ_s^{}d\mathrm{\Omega }_s^{}}}={\displaystyle \frac{3\sigma _T}{16\pi }}\delta (ϵ_s^{}ϵ_{\mathrm{scat}}^{})[(1\mu _{}^{}{}_{}{}^{2})(1\mu _{s}^{}{}_{}{}^{2})+`$ (9)
$`\text{ }+{\displaystyle \frac{1}{4}}(1+\mu _{}^{}{}_{}{}^{2})(1+\mu _{s}^{}{}_{}{}^{2})(g_1+g_2)]`$
(eg. Herold herold (1979)), where $`d\mathrm{\Omega }_s^{}=d\varphi _s^{}d\mu _s^{}`$ is an increment of solid angle into which outgoing photons with energy $`ϵ_s^{}`$ in the ERF are directed. As in the case of total cross section (6), when $`ϵ^{}>2ϵ_B`$ we replace the factor $`(g_1+g_2)`$ in (9) with $`(2\sigma _{KN}/\sigma _T+g_2g_1)`$. The sampled value of $`\mu _s^{}`$ determines the energy $`ϵ_{\mathrm{scat}}^{}`$ for the particular scattered photon with the relativistic formula:

$`ϵ_{\mathrm{scat}}^{}=(1\mu _{s}^{}{}_{}{}^{2})^1\{\text{}1+ϵ^{}(1\mu ^{}\mu _s^{})+`$ (10)
$`\text{ }[1+2ϵ^{}\mu _s^{}(\mu _s^{}\mu ^{})+ϵ_{}^{}{}_{}{}^{2}(\mu _s^{}\mu ^{})^2]^{1/2}\}`$
appropriate for collisions with a recoiled electron remaining at the ground Landau level (Herold herold (1979)). Finally, a value of the electron’s longitudinal momentum in the ERF is changed due to recoil from zero to $`(ϵ^{}\mu ^{}ϵ_{\mathrm{scat}}^{}\mu _s^{})mc`$, and transformed back to the OF.

The Monte Carlo method will be compared with a smooth integration treatment of the magnetic ICS (Chang chang (1995); Sturner sturner (1995)). In the integration method we assume the same procedure to accelerate an electron and to subtract its energy losses due to the curvature radiation as described above (Eq. 4). The only difference is in accounting for electron energy losses due to the magnetic inverse Compton scattering. These are estimated in *each* step (regardless the value of the optical depth for the scattering process) from the following formula for the mean electron energy loss rate:

$`\dot{\gamma }_{\mathrm{ICS}}=c{\displaystyle }dϵ{\displaystyle }d\mathrm{\Omega }\left({\displaystyle \frac{dn_{\mathrm{ph}}}{dϵd\mathrm{\Omega }}}\right)(1\beta \mu )\times `$ (11)
$`\text{ }\times {\displaystyle 𝑑ϵ_s^{}𝑑\mathrm{\Omega }_s^{}\left(\frac{d\sigma }{dϵ_s^{}d\mathrm{\Omega }_s^{}}\right)(ϵ_sϵ)}`$
where $`ϵ_s=ϵ_s^{}\gamma (1+\beta \mu _s^{})`$ is the scattered photon energy in the OF (eg. Dermer dermer (1990)). In other words, the electron Lorentz factor is determined as a solution of the differential equation:

$$\frac{d\gamma }{dh}=v^1(\dot{\gamma }_{\mathrm{acc}}+\dot{\gamma }_{\mathrm{ICS}}+\dot{\gamma }_{\mathrm{cr}}).$$
(12)
At each step, as an electron moves upwards, a decrease in both the dipolar magnetic field strength and the solid angle subtended by the thermal polar cap is taken into account in both methods.

## 3 Results

We compare the Monte Carlo method described in Sect. 2 with the integration method (numerical integration of Eq. (12)), for a pulsar with the polar magnetic field strength and rotation period as for B1509$-$58 ($`B_{\mathrm{pc}}=15.8\times 10^{12}`$ G, $`P=0.15`$ s).

### 3.1 The case of the M74 electric field

In this subsection the results are presented for the electric field of Michel (michel (1974)) (see Eq.1). The assumed values of $`B_{\mathrm{pc}}`$ and $`P`$ give $`E_{}=3.67\times 10^6`$ V cm<sup>-1</sup> and the standard polar cap radius $`r_{\mathrm{pc}}=3.74\times 10^4`$ cm. Following Sturner (sturner (1995)), we assumed the thermal polar cap radius $`R_{\mathrm{th}}=10^5`$ cm ($`2.7\times r_{\mathrm{pc}}`$) and $`h_{\mathrm{acc}}=r_{\mathrm{pc}}`$. Calculations have been performed for four different temperatures of the thermal polar cap: $`T_6=3.0`$, $`3.5`$, $`4.0`$ and $`4.5`$, where $`T_6=T/(10^6\mathrm{K})`$.

Each panel in Fig. 1 presents in grey the curves $`\gamma (h)`$ as calculated in the Monte Carlo way for one hundred individual electrons. The behaviour determined with Eq. (12) is overplotted as the thick solid line. The inclined dotted line presents changes of electron energy if there were no losses. It becomes horizontal at altitude $`h_{\mathrm{acc}}=r_{\mathrm{pc}}=3.74\times 10^4`$ cm above which the longitudinal component of electric field $`E_{}`$ is assumed to be screened by a charge-separated plasma distribution (see Eq. 1). This height is denoted by the vertical long dashed line.

The thick solid trajectories $`\gamma (h)`$ in four panels of Fig. 1 represent the ”averaged treatment” solutions and are in good agreement with the results of previous calculations by Sturner (sturner (1995)) (cf his Fig. 4). Below $`h10`$ cm they overlap with the case with no energy losses (dotted), since the acceleration rate of electron (as given by $`\dot{\gamma }_{\mathrm{acc}}=\frac{eE_{}}{mc^2}\beta c`$) significantly exceeds the energy loss rate due to the ICS (see Fig. 2a). Energy losses due to the CR are negligible for any Lorentz factor accessible for an electron, given the assumed acceleration model of Michel (Eq. 1). At larger altitudes the thick trajectories illustrate the acceleration stopping effect noted by Chang (chang (1995)) and Kardashëv et al. (kmn (1984)): the energy that the electrons would have gained between altitudes $`h10^2`$ and $`10^4`$ cm due to the electric field $`E_{}`$ is transferred to the thermal photons through the resonant inverse Compton scatterings. Depending on which process – either the resonant ICS cooling or the acceleration – ceases first, the electrons end up with energy which is either closer to $`eE_{}h_{\mathrm{acc}}`$ (ie. the energy acquired due to the full voltage drop at no radiative losses, the case $`T_6=3.0`$) or decreases towards a few tens $`\times mc^2`$ (the cases with $`T_63.5`$).

As can be seen in Fig. 1, the Monte Carlo tracks generally behave qualitatively in the same way: the increasing temperature of the thermal cap increases the altitude at which the acceleration takes over, and eventually, the Lorentz factors of most of electrons become limited below $`10^2`$. However, there are strong differences between the average energy of the Monte Carlo electrons and the value obtained with the integration method. They are especially pronounced for $`T_6=3.5`$, which is the case where a maximum altitude at which the acceleration can be counterbalanced by the electron energy losses due to the resonant ICS is equal to $`h_{\mathrm{acc}}`$. The bulk of the Monte Carlo tracks ends up with Lorentz factors $`\gamma >10^5`$, whereas the integration-method solution gives $`\gamma 64`$. To understand these differences it is worth investigating closely the behaviour of one exemplary electron as determined by the two methods.

In the integration method, the electron’s energy increases monotonically until an equilibrium settles between the energy loss process and the acceleration. Fig. 2b shows $`|\dot{\gamma }_{\mathrm{ICS}}|`$ and $`\dot{\gamma }_{\mathrm{acc}}`$ as a function of the Lorentz factor $`\gamma `$ for three different altitudes in the case $`T_6=3.0`$. One can see that $`|\dot{\gamma }_{\mathrm{ICS}}|`$ exceeds $`\dot{\gamma }_{\mathrm{acc}}`$ between the two equilibrium Lorentz factors $`\gamma _{\mathrm{eq1}}`$ and $`\gamma _{\mathrm{eq2}}`$ the values of which are determined by the condition $`|\dot{\gamma }_{\mathrm{ICS}}|=\dot{\gamma }_{\mathrm{acc}}`$. We have calculated them using the approximation $`\dot{\gamma }_{\mathrm{ICS}}\dot{\gamma }_{\mathrm{res}}`$ where $`\dot{\gamma }_{\mathrm{res}}`$ is the energy loss rate due to the resonant part of the cross section alone (cf. Eq. (22) in Sturner sturner (1995)). The two equilibrium Lorentz factors $`\gamma _{\mathrm{eq1}}`$ and $`\gamma _{\mathrm{eq2}}`$ as calculated for different altitudes $`h`$ are presented in Fig. 1 by the thick dashed line. For increasing altitude their values approach each other because a decrease in the density of black-body photons moves down the whole curve $`|\dot{\gamma }_{\mathrm{ICS}}(\gamma )|`$ in Fig. 2b. Thus, the energy-altitude space is divided into two regimes: the region where $`|\dot{\gamma }_{\mathrm{ICS}}|>\dot{\gamma }_{\mathrm{acc}}`$, hereafter called the ”energy loss dominated” region, (it is surrounded by the thick dashed line in Fig. 1), and the ”acceleration dominated” region where $`|\dot{\gamma }_{\mathrm{ICS}}|<\dot{\gamma }_{\mathrm{acc}}`$. As can be seen in Fig. 1, according to the integration method the trajectory of an electron in the energy-altitude space does not penetrate the energy loss dominated region. As the electron moves upwards, its Lorentz factor settles at the lower equilibrium value $`\gamma _{\mathrm{eq1}}`$, and stays there as long as the equilibrium can exist ie. until the curves $`|\dot{\gamma }(\gamma )|`$ and $`\dot{\gamma }_{\mathrm{acc}}(\gamma )`$ in Fig. 2b disconnect. For $`T_6=3.0`$ this happens near the altitude $`h_{\mathrm{eq}}=2\times 10^4`$. Above the height $`h_{\mathrm{eq}}`$ the acceleration cannot be counterballanced by the ICS energy losses at any Lorentz factor accessible for the electron and the increase in the electron energy resumes.

In the Monte Carlo method, the electron’s energy losses due to the scattering occur in a discontinuous way, thus, $`\gamma `$ oscillates around the value $`\gamma _{\mathrm{eq1}}`$ (see Fig. 3). When the energy loss rate $`|\dot{\gamma }_{\mathrm{ICS}}|`$ hardly exceeds the acceleration rate (see Fig. 2) the difference between the equilibrium Lorentz factors $`(\gamma _{\mathrm{eq1}}\gamma _{\mathrm{eq2}})`$ is small and the electron is able to pass *through* the energy-loss dominated region. Eg. for the case $`T_6=3.0`$, $`\gamma _{\mathrm{eq2}}22004\times \gamma _{\mathrm{eq1}}`$ at $`h=10^3`$ cm. To increase its energy from $`\gamma _{\mathrm{eq1}}`$ to $`\gamma _{\mathrm{eq2}}`$ it is enough for the electron to avoid a scattering over a ‘break-through’-distance $`\mathrm{\Delta }_{\mathrm{bt}}=\frac{\gamma _{\mathrm{eq2}}\gamma _{\mathrm{eq1}}}{eE_{}/mc^2}=228\mathrm{cm}`$ which is only ‘a few’ times larger than the local mean free path for scattering $`\lambda _{\mathrm{ICS}}`$. The local value of $`\lambda _{\mathrm{ICS}}`$ depends strongly on the electron Lorentz factor and for $`\gamma `$ between $`\gamma _{\mathrm{eq1}}`$ and $`\gamma _{\mathrm{eq2}}`$ it ranges from $`28`$ to $`108`$ cm (the case $`T_6=3.0`$, $`h=10^3`$ cm). This gives $`\mathrm{\Delta }_{\mathrm{bt}}(82)\times \lambda _{\mathrm{ICS}}`$. Thus, there is a very large probability for the electron to gain $`\gamma >\gamma _{\mathrm{eq2}}`$ at an altitude $`h_{\mathrm{MC}}`$ lower than $`h_{\mathrm{eq}}`$. Once the electron enters the acceleration dominated region at $`h_{\mathrm{MC}}h_{\mathrm{eq}}`$, its energy starts to increase up to a value much larger than obtained in the integration method.

For increasing temperatures $`T`$, the width $`(\gamma _{\mathrm{eq2}}\gamma _{\mathrm{eq1}})`$ of the energy-loss dominated region increases (Fig. 1) whereas the mean free path for $`\gamma =\gamma _{\mathrm{eq1}}`$ decreases. This makes the diffusion of electrons through the energy loss dominated region more difficult: to increase its energy from $`\gamma _{\mathrm{eq1}}`$ to $`\gamma _{\mathrm{eq2}}`$ an electron must avoid a scattering over a distance which is increasing multiplicity of the local mean free path. As a result, the energy distribution for outgoing electrons becomes softer (Fig. 4).

It should be emphasized that the strong disagreement between the final electron energies as determined with the two methods only appears if conditions similar to those for the case $`T_6=3.5`$ in Fig. 1 are fulfilled. These include the equilibrium between the maximum rate of resonant energy losses and the rate of acceleration at $`h=h_{\mathrm{acc}}`$. Making use of Dermer’s approximation for $`|\dot{\gamma }_{\mathrm{res}}|`$ (Dermer dermer (1990)) one can easily find that it has the maximum at $`\gamma _{\mathrm{res}}=ϵ_B/(w(1\beta \mu _{\mathrm{min}})𝒯)`$ where $`\mu _{\mathrm{min}}=h/(h^2+R_{\mathrm{th}}^2)^{1/2}`$ is a cosine of an angle at which the thermal cap radius is seen from the position of electron and $`w=\mathrm{ln}0.5`$. The equilibrium $`|\dot{\gamma }_{\mathrm{res}}|=\dot{\gamma }_{\mathrm{acc}}`$ holds for

$$T_6^{\mathrm{eq}}=4.14\left(\frac{\beta ^2\frac{eE_{}}{mc^2}}{(1\beta \mu _{\mathrm{min}})B_{12}(h)}\right)^{1/2}$$
(13)
where $`B_{12}(h)=B(h)/(10^{12}\mathrm{G})`$ and $`\frac{eE_{}}{mc^2}`$ is in cm<sup>-1</sup> (cf Eq. (29) in Chang chang (1995)). For the considered accelerating field (Eq. 1) this condition becomes

$$T_6^{\mathrm{eq}}=3.46\beta (1\beta \mu _{\mathrm{min}})^{1/2}\left(B_{12}P\right)^{1/4}$$
(14)
and is shown in Fig. 5 for $`h=0`$ and $`P=0.15`$ s as the solid line.

Other conditions required for the diffusion are low thermal cap temperatures and high magnetic field strengths. Both requirements partially stem from the dependence $`\lambda _{\mathrm{ICS}}(\gamma _{\mathrm{res}})BT^3`$ which results in $`\lambda _{\mathrm{ICS}}\mathrm{\Delta }_{\mathrm{bt}}`$ for high $`T`$ and low $`B_{\mathrm{pc}}`$. \[Note that at the resonance $`\lambda _{\mathrm{ICS}}`$ is *proportional to* $`B`$ just as $`|\dot{\gamma }_{\mathrm{res}}|`$\].

Additionally, high temperatures and low $`B`$-field strengths preclude the diffusion because of relatively large energy losses due to scatterings in the Klein-Nishina regime (see Fig 2a). First, for $`\gamma >\gamma _{\mathrm{res}}`$, $`|\dot{\gamma }_{\mathrm{KN}}|\dot{\gamma }_{\mathrm{acc}}`$ must hold, which requires low $`T`$. Second, the resonant bump in the $`|\dot{\gamma }_{\mathrm{ICS}}(\gamma )|`$ curve must be present, which is possible when $`|\dot{\gamma }_{\mathrm{res}}(\gamma _{\mathrm{res}})||\dot{\gamma }_{\mathrm{KN}}|`$. Since $`\dot{\gamma }_{\mathrm{res}}(\gamma _{\mathrm{res}})/\dot{\gamma }_{\mathrm{KN}}B`$, strong magnetic fields are preferred.

For $`P0.15`$ s and $`R_{\mathrm{th}}>h_{\mathrm{acc}}`$ all these constraints limit the diffusion regime to $`T<5\times 10^6`$ K and $`B_{\mathrm{pc}}>3\times 10^{12}`$ G with the $`T`$ and $`B_{\mathrm{pc}}`$ roughly fulfilling Eq. (14) taken for $`h=h_{\mathrm{acc}}`$. The region of strong discrepancies between the two methods is shown in Fig. 5 (hatched).

Among a few pulsars exhibiting a two-component thermal X-ray spectra only the Vela pulsar has an inferred temperature of the thermal polar cap ($`T15.8\times 10^6`$ K, Ögelman ogelman (1995)) which is sufficient for the acceleration to be initially halted by the resonant ICS (ie. it exceeds $`T_6^{\mathrm{eq}}`$ for $`hR_{\mathrm{th}}`$, see Fig. 5). The equilibrium $`|\dot{\gamma }_{\mathrm{res}}|=\dot{\gamma }_{\mathrm{acc}}`$ holds up to the altitude $`h_{\mathrm{eq}}2.2\times 10^4`$ cm which is two times larger than the inferred thermal cap radius ($`10^4`$ cm, Sturner sturner (1995)) and two times smaller than $`r_{\mathrm{pc}}=4.8\times 10^4`$ cm. Nevertheless, because of the high $`T`$, at $`hh_{\mathrm{eq}}`$ the diffusion is precluded by the extremely small $`\lambda _{\mathrm{ICS}}0.1`$ cm. It becomes efficient only at $`h>0.9h_{\mathrm{eq}}`$, thus, most Monte-Carlo electrons follow closely the behaviour predicted by the continuous approach by Chang (chang (1995)) or Sturner (sturner (1995)).

For other pulsars with a likely thermal X-ray emission from heated polar caps (eg. Geminga, B1055$-$52, and B0656$+$14) inferred temperatures are much lower (between 2 and 4 $`\times 10^6`$ K). Therefore, $`|\dot{\gamma }_{\mathrm{ICS}}|\dot{\gamma }_{\mathrm{acc}}`$ for any Lorentz factor accessible for an electron (Fig. 5) and the Compton scattering cannot considerably redistribute an initial electron energy spectrum.

Both, in the case of negligible ICS losses, ie. for $`TT^{\mathrm{eq}}(h=0)`$, and in the cases when $`|\dot{\gamma }_{\mathrm{res}}(\gamma _{\mathrm{res}})|\dot{\gamma }_{\mathrm{acc}}`$ holds up to $`h=h_{\mathrm{acc}}`$, ie. for $`TT^{\mathrm{eq}}(h=h_{\mathrm{acc}})`$, we find the energy distribution of electrons to be quasi-monoenergetic with energy well approximated by Eq.12. In the case of negligible ICS cooling the energy distribution of electrons is a narrow peak at the energy $`\gamma _{\mathrm{max}}mc^2=eE_{}h_{\mathrm{acc}}`$, whereas for very large efficiency of ICS losses (close to 100% of $`\gamma _{\mathrm{max}}mc^2`$) nearly all particles are cooled down to the energy $`\gamma _{\mathrm{min}}`$ for which the ICS loss rate decreases sharply (see Fig. 2a and Fig. 1). Between these limiting cases, ie. when the ICS cooling is comparable to the acceleration, broad energy distributions of electrons (ranging from $`\gamma _{\mathrm{min}}`$ to $`\gamma _{\mathrm{max}}`$) emerge (Fig. 4).

### 3.2 The case of the HM98 electric field

To enable electron energy losses due to the resonant ICS to compete with acceleration (the resonant ICS damping becomes then important and the diffusion effect does occur) the strength of the electric field $`E_{}`$ should not exceed a critical value which may be roughly estimated as $`E_{}^{\mathrm{eq}}310^4B_{12}T_6^2\mathrm{V}\mathrm{cm}^1`$ (cf. Chang chang (1995)). For $`E_{}(h)`$ given by eq. (18) of HM98 this may only occur in the case of small acceleration length $`h_{\mathrm{acc}}r_{\mathrm{pc}}`$ (see Eq. (2)), when the acceleration height is limited by the ICS-induced pair formation front. Self consistent values of $`h_{\mathrm{acc}}`$ are then of the order of a few$`\times 10^3`$ cm, i.e. much smaller than for a CR-induced case (cf. figs. 5 and 6 in HM98). A value of $`h_{\mathrm{acc}}`$ is considered self consistent if the acceleration of primary electrons by the electric field induced with this $`h_{\mathrm{acc}}`$ generates pair formation front at $`h_{\mathrm{acc}}`$.

We calculated self consistent $`h_{\mathrm{acc}}`$ in an approximate way, following HM98, by finding a minimum of a sum of acceleration length and a mean free path for one photon absorption:

$$h_{\mathrm{acc}}=\mathrm{min}\left\{\left(\frac{2\gamma }{N}\right)^{1/2}+\left(\frac{0.2\rho _{\mathrm{curv}}}{ϵ_Bϵ}\right)\right\}$$
(15)
where to get the first term we used $`E_{}`$ in the limit $`hh_{\mathrm{acc}}`$, given by Eq. (2): $`E_{}(h)=Nh`$ with $`N=3\mathrm{\Omega }R_{\mathrm{ns}}B_{\mathrm{pc}}\kappa \mathrm{cos}\chi (c(1ϵ_{\mathrm{GR}})R_{\mathrm{ns}}^2)^1h_{\mathrm{acc}}`$. As an energy of pair producing photons in the second term we assumed that for the resonant scattering: $`ϵ\gamma ϵ_B`$. In Eq. (15) we neglected a contribution from the mean free path for the resonant ICS which is of the order of $`150(1\beta \mu _{\mathrm{min}})^1B_{12}T_6^3`$ cm near the resonance. Since in the regime $`h_{\mathrm{acc}}r_{\mathrm{pc}}`$ the electric field $`E_{}`$ depends itself on $`h_{\mathrm{acc}}`$ we had to repeat the calculation of $`h_{\mathrm{acc}}`$ starting from some guess value until convergence. For $`B_{\mathrm{pc}}=15.8`$ TG, $`P=0.15`$ s, $`\rho _{\mathrm{curv}}=10^7`$ cm and $`\chi =0.1`$ rad we obtained $`h_{\mathrm{acc}}2\times 10^3`$ cm.

Then we calculated the changes of $`\gamma `$ with $`h`$ for $`R_{\mathrm{th}}=10^5`$ cm and for different thermal cap temperatures (as listed in the preceding subsection). We found that the most notable divergence between the Monte Carlo results and the integration method results occurs for $`T_6=3.5`$. As can be seen in Fig. 6a, at $`h10^3`$ cm most Monte Carlo electrons reach Lorentz factor values between $`10^3`$ and $`10^4`$, contrary to $`470`$ as predicted by the integration method. The difference disappears only above the accelerator (i.e. for $`h>h_{\mathrm{acc}}`$) because such altitude is still low enough (i.e. $`hR_{\mathrm{th}}`$) for the resonant ICS to damp $`\gamma `$ to its final value of about $`50`$ (at which it becomes eventually negligible).

An effect similar to that shown in Fig. 1 appears, however, when the resonant ICS losses are no longer effective above $`h=h_{\mathrm{acc}}`$. Fig. 6b shows the case with $`R_{\mathrm{th}}=4.0\times 10^3`$ cm $`=2\times h_{\mathrm{acc}}`$ and $`T_6=3.0`$. Broad electron energy distribution (with average energy much exceeding a value from the integration method) emerges.

Similarly as in the weak electric field case, the position of the diffusion region in $`T`$-$`B_{\mathrm{pc}}`$ plane (for a fixed value of $`R_{\mathrm{th}}`$) is here determined with a condition analogous to Eq. (13) (but $`E_{}`$ of Eq.(2) requires a more accurate determination of $`h_{\mathrm{acc}}`$) along with a requirement for strong $`B_{\mathrm{pc}}`$ and low $`T`$.

For stronger electric fields ($`E_{}>E_{}^{\mathrm{eq}}`$), or for larger acceleration lengths the resonant ICS has negligible influence on final electron energies. In the first case acceleration greatly dominates the ICS damping; in the second case acceleration is stopped only over a negligible, initial part of accelerator height.

## 4 Conclusions

We have calculated the efficiency of energy losses for electrons accelerated over a hot polar cap of a neutron star. As cooling mechanisms we have considered the inverse Compton scattering and the curvature radiation, (though the latter has been negligible for the considered acceleration models).

The electron energy losses due to the ICS have been calculated with the Monte Carlo method and compared with the integration method (after Chang chang (1995); Sturner sturner (1995)) based on the prescription for the averaged ICS cooling (Dermer dermer (1990)). We confirm general predictions of the integration approach (eg. the limitation of electron Lorentz factors below $`10^2`$ for high temperatures), nevertheless, the “stopping acceleration effect” occurs at a slightly higher temperature than that one resulting from integrating the differential equation (12).

For the considered pulsar parameters ($`B_{\mathrm{pc}}=15.8\times 10^{12}`$ G, $`R_{\mathrm{th}}=10^5`$ cm) and the acceleration potential by Michel (michel (1974)), the Monte Carlo-based value is approximately equal to $`T_{\mathrm{MC}}=4.5\times 10^6`$ K (at this $`T`$ only a few percent of electrons reach $`\gamma >10^2`$ at $`h=10^6`$ cm) in comparison with $`T_{\mathrm{int}}=3.3\times 10^6`$ as given by the integration approach. Although the temperature difference is small, the final electron Lorentz factors as determined by the two methods may easily differ by a few orders of magnitude (Fig. 1, the case $`T_6=3.5`$) if the observed temperature matches this range. Moreover, it should be kept in mind that for most objects with hard X-ray emission identified putatively as originating from the hot polar cap, temperatures lie within the similar range $`(24)\times 10^6`$ K (Ögelman ogelman (1995)).

We find that for high-$`B`$ and low-$`T`$ pulsars the integration method gives a poor estimation of electron Lorentz factors if the energy loss rate due to the resonant ICS hardly exceeds the rate of acceleration. In the case of Michel’s model and for $`R_{\mathrm{th}}>h_{\mathrm{acc}}`$, this occurs if $`T4(B_{12}P)^{1/4}`$ which may be the case for PSR B1509$-$58. In the cases when the ICS loss rate dominate firmly over the acceleration or is negligible the energy distribution of outgoing electrons is quasi-monoenergetic around the value which is well approximated with the continuous approach.

We find that preventing electrons from achieving large Lorentz factors occurs also for accelerating potential by Harding & Muslimov (hm (1998)) in its nonsaturated version with $`h_{\mathrm{acc}}r_{\mathrm{pc}}`$. The discrepancy between the Monte Carlo and the integration method also occurs within the high-$`B`$ and low-$`T`$ regime and is especially pronounced for $`R_{\mathrm{th}}h_{\mathrm{acc}}`$.

###### Acknowledgements.

This work was supported by KBN grants 2P03D 01016 and 2P03D 02117. JD thanks T.Bulik for useful advice in programming. We appreciate sugestions by the anonymous referee and comments by Bing Zhang on the manuscript which helped to improve the paper.


format: text

text:
# 1 Introduction

## 1 Introduction

One of the basic open questions of particle physics is the nature of electroweak symmetry breaking. Beyond the discovery of the Higgs boson(s), however, reconstructing the Higgs potential will be necessary, and that requires the experimental study of the self–couplings of the Higgs bosons.

Pair production of neutral Higgs particles in gluon fusion sensitive to the trilinear couplings of Higgs bosons was studied in the SM  and the Minimal Supersymmetric Standard Model (MSSM) \[2–5\]. QCD corrections were included in the limiting case of a heavy top mass .

In the SM the trilinear and the quartic couplings of the physical Higgs particle are fixed by the Higgs mass. There is no resonance in the pair production of SM Higgses, and the cross section is only $``$20–50 fb in the intermediate mass range . No resonance effect is present in the production of three Standard Model Higgs bosons via one Higgs in gluon fusion and it is estimated to be under the discovery limit at LHC.

In the MSSM there are five Higgs bosons ($`h,H,A,H^\pm `$) and many trilinear and quartic couplings among them. Two parameters describe the Higgs sector of the MSSM: $`M_A`$, the mass of the CP-odd Higgs boson A, and tan$`\beta `$, the ratio of the two vacuum expectation values. For a wide range of $`M_A`$, in $`pp`$ collision only the cross section of the hh production is large . Possible other processes, such as the gauge boson fusion and Higgs–strahlung off W bosons or heavy quarks (top) provide less event number.

In order to provide further tests of trilinear and quartic Higgs self–couplings, in this paper we calculate the resonance enhanced production of three lightest supersymmetric Higgs bosons (h) in gluon fusion

$$ppgghhh.$$
(1)
Generally four different graphs contribute to the production of three Higgses. Gluons couple to triangle, box or pentagon quark loops emitting 1,2 or 3 Higgs bosons, as is seen in Fig.1. Squark loops are neglected. We have omitted the pentagon graph producing a flat continuum background for the resonance production and we will give the complete cross section valid in a wider range in a subsequent publication . We learn that the quartic coupling is not accessible by LHC experiment. The resonance contribution to (1), however, sizeable at LHC For instance at $`\mathrm{tan}\beta =3`$ it yields a cross section of about 300 fb. This is about $`1/10`$ times the corresponding $`hh`$–case .

The masses, widths and the couplings were calculated using the complete one–loop and the leading two–loop radiative corrections from . The relevant Higgs self–couplings are the following

$`\lambda _{hhh}`$ $`=3\mathrm{cos}(2\alpha )\mathrm{sin}(\beta +\alpha )+{\displaystyle \frac{3ϵ}{M_Z^2}}{\displaystyle \frac{\mathrm{cos}^3\alpha }{\mathrm{sin}\beta }}`$ (2)
$`\lambda _{Hhh}`$ $`=2\mathrm{sin}(2\alpha )\mathrm{sin}(\beta +\alpha )\mathrm{cos}(2\alpha )\mathrm{cos}(\beta +\alpha )+{\displaystyle \frac{3ϵ}{M_Z^2}}{\displaystyle \frac{\mathrm{sin}\alpha \mathrm{cos}^2\alpha }{\mathrm{sin}\beta }}`$
$`\lambda _{HHh}`$ $`=2\mathrm{sin}(2\alpha )\mathrm{cos}(\beta +\alpha )\mathrm{cos}(2\alpha )\mathrm{sin}(\beta +\alpha )+{\displaystyle \frac{3ϵ}{M_Z^2}}{\displaystyle \frac{\mathrm{sin}^2\alpha \mathrm{cos}\alpha }{\mathrm{sin}\beta }}`$
$`\lambda _{Hhhh}`$ $`=3\mathrm{sin}(2\alpha )\mathrm{cos}(2\alpha )+{\displaystyle \frac{ϵ}{8M_Z^2}}{\displaystyle \frac{\mathrm{cos}\alpha \mathrm{sin}^2\alpha }{\mathrm{sin}^2\beta }}`$
$`\lambda _{hhhh}`$ $`=3\mathrm{cos}^2(2\alpha )+{\displaystyle \frac{ϵ}{4M_Z^2}}{\displaystyle \frac{\mathrm{sin}^4\alpha }{\mathrm{sin}^2\beta }}`$
The trilinear and quartic couplings are normalized to $`\lambda _3=[\sqrt{2}G_F]^{1/2}M_Z^2`$ and $`\lambda _4=\sqrt{2}G_FM_Z^2`$, respectively. The couplings depend on $`\beta `$ and the mixing angle $`\alpha `$ of the CP–even Higgs sector

$$\mathrm{tan}2\alpha =\frac{M_A^2+M_Z^2}{M_A^2M_Z^2+ϵ/\mathrm{cos}2\beta }\mathrm{tan}2\beta .$$
(3)
The leading $`m_t^4`$ one–loop corrections  are parametrized by

$$ϵ=\frac{3G_F}{\sqrt{2}\pi ^2}\frac{m_t^4}{\mathrm{sin}^2\beta }\mathrm{log}\left[1+\frac{M_S^2}{m_t^2}\right].$$
(4)
The common squark mass is fixed to $`M_S`$= 1 TeV. The quartic couplings depend on $`M_A`$ and are shown for two values of $`\text{tan}\beta =3`$ and 30 in Fig.2.

The paper is organized as follows. In the next section we give the cross section of the resonance production of three h’s in gluon fusion and in section 3 we show the numerical results.

## 2 The cross section

Three different channels contribute to the resonance production of three lightest supersymmetric Higgs particles (h) in gluon fusion. (a) One virtual Higgs boson (h,H) is produced by the heavy quark triangle and it decays into 3 h’s via the quartic coupling $`\lambda _{Hhhh}`$ or $`\lambda _{hhhh}`$, Fig. 1a. (b) One virtual Higgs boson decays into 3 h’s in two steps testing trilinear couplings, $`ggH,h(H,h)hhhh`$, Fig. 2b. (c) Heavy quark box diagram couples to two Higgses, one of them decays subsequently into two h’s, $`gg(H,h)hhhh`$, Fig. 2c. We have omitted the pentagon graphs contributing only to the continuum production. We believe it increases the production at large tan$`\beta `$ due to the large hbb–coupling, but does not change our conclusion about the measureability of the quartic Higgs–couplings.

There are two different helicity amplitudes contributing to the total cross section: the total spin of the two gluons along the collision axis can be $`S_z=0`$ (contribution F) or $`S_z=2`$ (contribution G). The triangle graphs give only contributions $`F_3`$ ($`S_z=0`$) as they involve a single spin–0 Higgs intermediate state. Box graphs give both contributions $`F_4`$ ($`S_z=0`$) and $`G_4`$ ($`S_z=2`$). $`F_4`$ and $`G_4`$ are Lorentz and gauge invariant decompositions of the box amplitude. The tensor basis is the following

$`S_z=0A^{\mu \nu }=g^{\mu \nu }{\displaystyle \frac{p_1^\nu p_2^\mu }{(p_1p_2)}}`$ (5)
$`S_z=2B^{\mu \nu }=g^{\mu \nu }+{\displaystyle \frac{1}{p_T^2(p_1p_2)}}\left(p_3^2p_1^\nu p_2^\mu 2(p_2p_3)p_1^\nu p_3^\mu 2(p_1p_3)p_2^\mu p_3^\nu +2p_3^\mu p_3^\nu (p_1p2)\right)`$
where $`p_1,p_2`$ are the momenta of the incoming gluons and $`p_3`$ is one of the momenta of the outgoing Higgs bosons ($`p_3,p_4,p_5`$). Here $`p_T^2=2\frac{(p_1p_3)(p_2p_3)}{(p_1p_2)}p_3^2`$ is the transverse momentum of the third particle. Tensors $`A^{\mu \nu }`$ and $`B^{\mu \nu }`$ are orthogonal and normalized to $`2`$ .

The M–matrix of the process is

$$=\frac{(\sqrt{2}G_F)^{3/2}\alpha _s\widehat{s}}{4\pi }ϵ_a^\mu ϵ_b^\nu \delta _{ab}(FA_{\mu \nu }+GB_{\mu \nu }).$$
(6)
The spin and color averaged parton level cross section of the process $`gghhh`$ is

$$d\widehat{\sigma }(gghhh)=\frac{\sqrt{2}G_F^3\alpha _s^2}{1024(2\pi )^6}\left[|F|^2+|G|^2\right]\frac{\lambda ^{1/2}(\widehat{s}_{45},p_4^2,p_5^2)}{\widehat{s}_{45}}d\widehat{s}_{45}d\widehat{s}_{13}d\mathrm{\Omega }_{45}^{CM}.$$
(7)
Here we used the Chew-Low parametrization of the three particle phase space , a trivial angle integration was carried out, $`\widehat{s}_{ik}=(p_i+p_k)^2`$, and the usual lambda function is $`\lambda (x,y,z)=(x(\sqrt{y}+\sqrt{z})^2)(x(\sqrt{y}\sqrt{z})^2)`$. $`d\mathrm{\Omega }_{45}^{CM}`$ is the differential solid angle in the center of mass system of particles 4 and 5.

$`F`$ receives contributions from the triangle and box graphs, while only box diagrams give contributions $`S_z=2`$ to $`G`$

$`F`$ $`=`$ $`{\displaystyle \underset{t,b}{}}(C_3F_3+{\displaystyle \underset{i=3,4,5}{}}C_4^{(i)}F_4^{(i)})`$
$`G`$ $`=`$ $`{\displaystyle \underset{t,b}{}}{\displaystyle \underset{i=3,4,5}{}}C_4^{(i)}G_4^{(i)}.`$ (8)
Here we separated the functions $`F_3,F_4,G_4`$ responsible for the heavy quark triangle and boxes, $`C_3,C_4`$ contain the coupling constants and the propagators. The amplitude has to be summed over all quark flavours. However, the four light quarks can be neglected having very small Higgs couplings. Supersymmetric particles are assumed to be too heavy to participate in the triangle or box loop. The second sum in the box contribution corresponds to the outgoing Higgs particle that couples directly to the quark loop.

We calculated $`F_3,F_4,G_4`$ and found them in agreement with the results of , . For instance,

$$F_3=2\frac{m_Q^2}{\widehat{s}}\left(2+(4m_Q^2\widehat{s})C_{12}\right),\widehat{s}=\widehat{s}_{12}.$$
(9)
where

$$C_{ij}=\frac{d^4q}{i\pi ^2}\frac{1}{(q^2m_Q^2)[(q+p_i)^2m_Q^2][(q+p_i+p_j)^2m_Q^2]}.$$
(10)
$`F_4,G_4`$ have long analytic expressions which can be found in ref . The box functions $`F_4^{(i)},G_4^{(i)}`$ in (8) depend on three momenta $`(p_1,p_2,p_i),i=3,4,5`$ and $`F_4^{(i)}=F_4(p_1,p_2,p_i)`$ etc.

In case of a large quark mass ($`m_Q^2\widehat{s}M_{h,H}^2`$) we get as in

$$F_3=\frac{2}{3}+𝒪\left(\widehat{s}/m_Q^2\right),F_4^{(i)}=\frac{2}{3}𝒪\left(\widehat{s}/m_Q^2\right),G_4^{(i)}=𝒪\left(\widehat{s}/m_Q^2\right)$$
(11)
In the limit of light quark masses ($`m_Q^2\widehat{s}M_{h,H}^2`$) all the form factors vanish to $`𝒪\left(m_Q^2/\widehat{s}\right)`$.

The generalized triangle couplings are the following

$`C_3=C_3^h+C_3^H+C_3^{hh}+C_3^{hH}+C_3^{Hh}+C_3^{HH}`$
$`C_3^{H_a}=\lambda _{hhhH_a}{\displaystyle \frac{M_Z^2}{\widehat{s}M_{H_a}^2+iM_{H_a}\mathrm{\Gamma }_{H_a}}}g_Q^{H_a}`$ (12)
$`C^{H_aH_b}=\lambda _{H_aH_bh}{\displaystyle \frac{M_Z^2}{\widehat{s}M_{H_a}^2+iM_{H_a}\mathrm{\Gamma }_{H_a}}}\lambda _{H_bhh}{\displaystyle \frac{M_Z^2}{\widehat{s}M_{H_b}^2+iM_{H_b}\mathrm{\Gamma }_{H_b}}}g_Q^{H_a},H_{a,b}=h/H.`$
The box couplings are

$`C_4=C_4^h+C_4^H`$
$`C_4^{(h,H)(i)}=\lambda _{hh(h/H)}{\displaystyle \frac{M_Z^2}{\widehat{s}_{kl}M_{h/H}^2+iM_{h/H}\mathrm{\Gamma }_{h/H}}}g_Q^hg_Q^{h/H},`$ (13)
$`i,k,l`$ are cyclic. $`g_Q^{h/H}`$ denotes the Higgs–quark couplings normalized to the SM Yukawa coupling $`[\sqrt{2}G_F]^{1/2}m_Q`$,

$$g_t^h=\mathrm{cos}\alpha /\mathrm{sin}\beta ,g_b^h=\mathrm{sin}\alpha /\mathrm{cos}\beta ,g_t^H=\mathrm{sin}\alpha /\mathrm{sin}\beta ,g_b^H=\mathrm{cos}\alpha /\mathrm{cos}\beta .$$
(14)
## 3 Results

Now, we calculate the production cross section of three lightest neutral supersymmetric Higgs bosons $`pphhh+X`$ via gluon fusion at the LHC energy of $`\sqrt{s}=14`$ TeV in the above picture by folding the parton level fusion cross section with the gluon luminosity. We used the GRV structure functions  for the gluon luminosity at $`Q^2=\widehat{s}`$. The numerical integration was made by the VEGAS package . The total cross section vs. $`M_A`$ is shown in Fig.3 for two represantative values of tan$`\beta =3`$ and $`\text{tan}\beta =30`$. In the plotted range of $`M_A`$ the heavy CP–even Higgs boson (H) is nearly degenerate in mass with A while the mass of the lightest neutral Higgs h is quickly approaches approximately 104 GeV from below.

For $`\text{tan}\beta =3`$ in the range of $`M_A=`$ 200–350 GeV the cross section is enhanced by resonance effect and large. The cross section goes up to nearly 300 fb and slowly decreases to about 150 fb giving a large number of events at the expected luminosity $`L𝑑t100fb^1/year`$. The main contribution comes from the box diagram where the heavy CP–even Higgs boson H couples to the quark loop and decays into two bosons h. $`M_h`$ slowly increases around 100 GeV when $`M_HM_A`$ reaches 200 GeV the H propagator enters the resonance region. At $`M_H2m_t`$ there is a threshold effect due to the fall–off of the branching ratio $`BR(Hhh)`$ and partly because on–shell top quark pairs can be produced in the quark loop. For large values of $`M_A`$ the cross section reaches a continuum value of a few fb. The dip at small $`M_A`$, similarly to the case of $`\text{tan}\beta =30`$, is the result of the zero in the trilinear couplings $`\lambda _{(H/h)hh}`$.

For $`\text{tan}\beta =30`$ the main contribution comes from the box diagram where an h couples to the quark loop, propagates and decays into two h’s via $`\lambda _{hhh}`$. The Hhh coupling is small compared to the case of $`\text{tan}\beta =3`$ while the hhh coupling is larger giving a sizeable continuum. There is no observable resonance effect and the cross section is nearly constant, 6 fb, versus $`M_A`$ after the coupling constant changed sign. The K factor is expected to increase the cross section as in other cases .

The contributions of the diagrams containing a heavy quark triangle are also shown in Fig.3 at $`\text{tan}\beta =3`$. The lowest curve corresponds to the diagram with quartic couplings, the one in the middle sums up all the contributions of the graphs involving a quark triangle. There is a clear resonance effect in both curves when $`M_H`$ reaches $`3M_h`$ at $``$310 GeV. Here the first propagator of $`ggHHhhhh`$ becomes resonant too. The first rise in the middle curve is the result of the resonance in the second propagator of Fig. 2b, similarly to the box. The triangle loops yield, however, only a small fraction of the cross section even at their peak at 300-350 GeV.

The suppressed contribution of the quartic coupling implies that in $`pphhh+X`$ the measurement of the quartic coupling is not possible. Not only the box contribution is larger by two orders of magnitude but also the other triangle diagrams are 10 times larger. For $`\text{tan}\beta =30`$ the case is even worse.

In conclusion, in this paper we have calculated the resonance contribution to the cross section of 3h–production at LHC via gluon fusion. The main contribution comes from the trilinear Higgs coupling while the quartic ones turn out to be negligible. In the resonance region at small tan$`\beta `$ the ratio of 2h and 3h production is about 10.

This work was partially supported by Hungarian Science Foundation Grants under Contract Nos. OTKA T029803 and F022998.

FIGURE CAPTIONS

Fig.1: Generic diagrams contributing to the production of three CP–even MSSM Higgs bosons in gluon–gluon collisions, $`gghhh`$: a) triangle with quartic coupling, b)triangle with trilinear couplings c) box , d) pentagon contributions.

Fig.2: Quartic Higgs couplings in the MSSM as functions of the pseudoscalar Higgs mass $`M_A`$ for two representative values of $`\text{tan}\beta =3`$ and 30.

Fig.3: Total cross sections for production of three lightest CP–even MSSM Higgs bosons in gluon–gluon collisions at the LHC for $`\text{tan}\beta =3`$ and 30 (upper two curves). The upper axis presents the scalar Higgs masses $`M_h`$ for $`\text{tan}\beta =3`$ corresponding to the pseudoscalar masses $`M_A`$. The lower two dashed curves represent the resonance contributions of the triangle graphs of Fig. 2a and 2b for $`\text{tan}\beta =3`$.