text:
# 1 Introduction

## 1 Introduction

Uncovering the mechanism of symmetry breaking is one of the major tasks of the high energy colliders. Most prominent is the search for the Higgs particle. Although the LHC should not miss this particle even if it weighed up to 1TeV, precision measurements on the Higgs properties will only be conducted in an $`e^+e^{}`$collider. There are two important mechanisms for Higgs production in $`e^+e^{}`$. The Higgs-strahlung process, $`e^+e^{}ZH`$and the $`W`$-fusion process, $`e^+e^{}\nu _e\overline{\nu }_eH`$. The former is the dominant one at small (in the LEP2 range say) to moderate energies but decreases rather fast with energy. At TeV energies the $`W`$-fusion process dominates by far for Higgs masses up to $`1`$TeV. Even at $`500`$GeV this $`t`$-channel process is dominant for Higgs masses in the range preferred by the indirect electroweak precision data and remains an important component for all Higgs masses at energies of the linear collider. Tree-level computations of Higgs production are rather well under control, including interference of the $`W`$ fusion process with the Higgs-strahlung process. Note however that the complexity of the process $`e^+e^{}\nu \overline{\nu }H`$precludes a full analytic result for the total cross section even at tree-level, although the differential cross section can be cast in a very compact form. Full radiative corrections for Higgs-strahlung have been considered by a number of groups, while a proper one-loop treatment of the fusion process is still lacking despite the importance of the process for the linear collider physics program. Some recipes have been suggested to include parts of the radiative corrections to the fusion process but considering the domain of validity of these approximations $`\sqrt{s},M_H2m_t`$, they are expected not to be precise for the interesting range of Higgs masses(preferred by the latest precision measurements) and next collider energies. One-loop contribution to the $`HWW`$ vertex has been considered on the basis that it might constitute a good approximation for the fusion process, but it rests to see how well this approximation fares in comparison of the full calculation. Very recently one-loop radiative corrections to this process have been investigated within the minimal supersymmetric model but again by only taking into account the contribution of the fermions and sfermions to the $`H/hWW`$ vertex. It is the aim of this letter to summarise the results of the full radiative corrections to single Higgs production in $`e^+e^{}`$, including both the fusion and Higgs-strahlung processes in the $`𝒮`$(Standard Model). We include both the virtual and soft corrections as well as the hard photon radiation. A longer paper will detail our computation and results and will look into the issue of finding approximations to the full result111Preliminary results have been presented at the Workshop RADCOR2002. At this meeting the FIRCLA group exposed their plans and techniques, different from ours concerning Feynman integration, for tackling the calculation of this process. While finalising this letter we also learnt of a calculation by A. Denner, S. Dittmaier, M. Roth and M. Weber, in preparation..

A standard hand calculation using the usual techniques could hardly be attempted for such $`23`$ processes at one-loop. Considering the ever increasing power of computers, the possibility of parallelisation and the fact that the whole procedure of perturbation theory consists of algorithms that can be directly translated on a computer it seems that most, if not all, complex calculations in high-energy physics can be automated. This is especially true for electroweak processes where various scales and masses enter the calculations. GRACE-loop  from which our results are derived is such a program. GRACE , the tree-level component of the system, has been tested and heavily used for tree-level cross sections up to 6-fermions in the final state. GRACE-loop has been exploited and checked thoroughly for a variety of $`22`$ processes in the electroweak theory . The system which requires as input, a model file that describes all the interaction vertices derived from a particular Lagrangian can generate all the necessary Feynman graphs together with their codes so that matrix elements can be generated before being processed for the calculation of the cross section and event generation. For loop processes, there is a symbolic manipulation stage (either FORM or REDUCE ) that handles all the Dirac and tensor algebra in $`n`$-dimension for all the interference terms between tree-level and 1-loop diagrams and automatically applies the Feynman trick for the propagator. This is then passed to a module that contains two libraries for the loop integration containing the FF package as well as an in-house numerical code. The system together with the one-loop renormalisation program is described in detail in . As far as the calculation of one-loop processes is concerned, a series of powerful tests are implemented in the code as described in  and as will be presented below for $`e^+e^{}\nu \overline{\nu }H`$.

## 2 Tree-level results, setting-up the loop calculation

Our input parameters for the calculation of $`e^+e^{}\nu \overline{\nu }H`$are the following. Throughout we expressed our results in terms of the fine structure constant in the Thomson limit $`\alpha ^1=137.0359895`$ and the $`Z`$ mass $`M_Z=91.1876`$GeV. Our on-shell renormalisation program uses $`M_W`$ as input parameter, nonetheless our numerical value of $`M_W`$ is derived through $`\mathrm{\Delta }r`$222We include NLO QCD corrections and two-loop Higgs effects. We take $`\alpha _s(M_Z^2)=.118`$ together with $`G_\mu =1.16639\times 10^5\mathrm{GeV}^2`$. $`M_W`$ thus changes as a function of $`M_H`$. For the the lepton masses we take $`m_e=0.510999`$MeV, $`m_\mu =105.6584`$MeV and $`m_\tau =1.777`$GeV. For the quark masses beside the top mass $`M_t=174`$GeV, we take the set $`M_u=M_d=58`$MeV, $`M_s=92`$MeV, $`M_c=1.5`$GeV and $`M_b=4.7`$GeV. With these values we calculate $`\mathrm{\Delta }\alpha (M_Z)=0.059258`$. With this we find for example that $`M_W=80.3767`$GeV for $`M_H=150`$GeV and $`M_W=80.3158`$GeV for $`M_H=350`$GeV. Especially for the Higgs-strahlung subprocess we require a $`Z`$-width. We have taken a constant fixed $`Z`$-width, $`\mathrm{\Gamma }_Z=2.4952`$GeV. Unless when otherwise stated our results refer to the full $`e^+e^{}\nu \overline{\nu }H`$, summing over all three types of neutrinos with, for electron neutrinos, the effect of interference between fusion and Higgs-strahlung.

We have checked that our tree-level results are in very good agreement with those in  after expressing them in terms of $`G_\mu `$. However since we are considering the effect of radiative corrections, within our scheme we prefer showing all our results using $`\alpha `$. We will only comment on the $`G_\mu `$ scheme at the end of this letter. We find for example for $`\sqrt{s}=500`$GeV and $`M_H=350`$GeV that $`\sigma _{\mathrm{tree},\mathrm{total}}=4.606`$fb at tree-level for the contribution of all three neutrinos. In an attempt to separate the different contributions to single Higgs production, we will refer to the $`s`$-channel as given by $`\sigma _s=3\times \sigma (e^+e^{}\nu _\mu \overline{\nu }_\mu H)`$. The bulk of this contribution is given by $`\sigma (e^+e^{}ZH)\times B_{Z\mathrm{inv}.},B_{Z\mathrm{inv}.}20\%`$. We will define the $`t`$-channel as $`\sigma _t=\sigma _{\mathrm{total}}\sigma _s`$, this implicitly means that the interference term is included in this contribution. These definitions will be carried over to the one-loop case as well.

In Fig. 2 we have also included the tree-level cross section. They are shown as a function of the centre-of-mass energy for a light Higgs of mass $`M_H=150`$GeV as well as a function of the Higgs mass at a centre-of-mass of $`500`$GeV. All integration over phase space are done with the help of BASES, see. These figures clearly show the importance of the $`t`$-channel contribution pointed out in the introduction. For a low Higgs mass of $`150`$GeV, although the $`s`$-channel still dominates at $`\sqrt{s}=300`$GeV, very quickly at $`\sqrt{s}=500`$GeV it is the $`t`$-channel that dominates. For the latter energy as the Higgs mass increases, the $`t`$-channel contribution drops much quickly than the $`s`$-channel, where both merge around $`M_H=380`$GeV to $`M_H=390`$GeV, but very quickly around the $`ZH`$ threshold for $`M_H408`$GeV, the $`s`$-channel drops precipitously leaving the $`t`$-channel as the sole contribution for the whole process.

Neglecting all Goldstone-electron coupling (proportional to the electron mass), one has at one-loop $`249`$ diagrams (for $`\nu _e\overline{\nu }_eH`$, and $`146`$ for $`\nu _\mu \overline{\nu }_\mu H`$) including $`15`$ pentagons (5-point functions) compared to only 2 diagrams at tree-level, one for the fusion process and one for the Higgs-strahlung. Keeping the electron Yukawa coupling one has a total of $`1350`$ diagrams (98 pentagons corresponding to 5-point functions) at one-loop. Although in running our program to derive cross sections we only use the set of 249 diagrams, to perform our extensive checks especially those of gauge-parameter independence, at the level of the differential cross section, we keep the full set of 1350 diagrams. It is impossible to show all the contributing diagrams here. They may be downloaded or visualised at this location. All these diagrams are generated and drawn by gracefig the Feynman diagrams generator of GRACE . A representative selection of diagrams is shown in Fig. 1.

The results of the calculation are checked by performing three kinds of tests at some random points in phase space. For these tests to be passed one works in quadruple precision. We first check the ultraviolet finiteness of the results. This test applies to the whole set of the virtual one-loop diagrams. In order to conduct this test we regularise any infrared divergence by giving the photon a fictitious mass (we set this at $`\lambda =10^{15}`$GeV). In the intermediate step of the symbolic calculation dealing with loop integrals (in $`n`$-dimension), we extract the regulator constant $`C_{UV}=1/\epsilon \gamma _E+\mathrm{log}4\pi `$, $`n=42\epsilon `$ and treat this as a parameter. The ultraviolet finiteness test gives a result that is stable over $`30`$ digits when one varies the dimensional regularisation parameter $`C_{UV}`$. This parameter could then be set to $`0`$ in further computation. The test on the infrared finiteness is performed by including both loop and bremsstrahlung contributions and checking that there is no dependence on the fictitious photon mass $`\lambda `$. We find results that are stable over $`23`$ digits when varying $`\lambda `$. An additional test concerns the bremsstrahlung part. It relates to the independence in the parameter $`k_c`$ which is a soft photon cut parameter that separates soft photon radiation and the hard photon performed by the Monte-Carlo integration.

Gauge parameter independence of the result is performed through a set of five gauge fixing parameters. For the latter a generalised non-linear gauge fixing condition has been chosen.

$`_{GF}`$ $`=`$ $`{\displaystyle \frac{1}{\xi _W}}|(_\mu ie\stackrel{~}{\alpha }A_\mu igc_W\stackrel{~}{\beta }Z_\mu )W^{\mu +}+\xi _W{\displaystyle \frac{g}{2}}(v+\stackrel{~}{\delta }H+i\stackrel{~}{\kappa }\chi _3)\chi ^+|^2`$ (2.1)
$`{\displaystyle \frac{1}{2\xi _Z}}(.Z+\xi _Z{\displaystyle \frac{g}{2c_W}}(v+\stackrel{~}{\epsilon }H)\chi _3)^2{\displaystyle \frac{1}{2\xi _A}}(.A)^2.`$
The $`\chi `$ represent the Goldstone. We take the ’t Hooft-Feynman gauge with $`\xi _W=\xi _Z=\xi _A=1`$ so that no “longitudinal” term in the gauge propagators contributes. Not only this makes the expressions much simpler and avoids unnecessary large cancelations, but it also avoids the need for high tensor structures in the loop integrals. The use of five parameters is not redundant as often these parameters check complementary sets of diagrams. For example the parameter $`\stackrel{~}{\beta }`$ is involved in all diagrams containing the gauge $`WWZ`$ and their Goldstone counterpart, whereas $`\stackrel{~}{\alpha }`$ checks $`WW\gamma `$ and $`\stackrel{~}{\delta }`$ is implicitly present in $`WWH`$. For each parameter of the set $`\zeta =(\stackrel{~}{\alpha },\stackrel{~}{\beta },\stackrel{~}{\delta },\stackrel{~}{\kappa },\stackrel{~}{ϵ})`$ the first check is made while freezing all other four parameters to $`0`$. We have also made checks with two parameters non-zero. In principle checking for $`2`$ or $`3`$ values of the gauge parameter should be convincing enough. We in fact go one step further and perform a complete gauge parameter independence. To achieve this we generate for each non-linear gauge parameter $`\zeta `$, the values of the loop correction to the total differential cross section as well as the contribution of each one-loop diagram contribution for the five values $`\zeta =0,\pm 1,\pm 2`$. The one-loop diagram contribution from each loop graph $`g`$, is defined as

$`\mathrm{d}\sigma _g=\mathrm{d}\sigma _g(\zeta )=\mathrm{}\left(𝒯_g^{loop}𝒯^{tree}\right).`$ (2.2)
$`𝒯^{tree}`$ is the tree-level amplitude summed over all tree-diagrams333Therefore the tree-level amplitude does not depend on any gauge parameter. For the process at hand nonetheless, some individual tree diagrams depend on the gauge parameter $`\stackrel{~}{\delta }`$ and $`\stackrel{~}{ϵ}`$ giving extremely small contributions proportional to the electron mass. Our numerical procedure to isolate the gauge parameter dependence detects these tiny variations, replacing $`𝒯^{tree}`$ by any tree-level diagram instead of the sum, one is able to differentiate between a variation in $`\mathrm{d}\sigma _g`$ due to the loop diagram or the residual one from $`𝒯^{tree}`$.. $`𝒯_g^{loop}`$ is the one-loop amplitude contribution of the one-loop diagram $`g`$. A rapid look at the structure of the Feynman rules of the non-linear gauge leads one to conclude that for $`e^+e^{}\nu \overline{\nu }H`$each contribution is a polynomial of (at most) third degree in the gauge parameter and thus, that each contribution, $`\mathrm{d}\sigma _g`$ may be written as

$`\mathrm{d}\sigma _g=\mathrm{d}\sigma _g^{(0)}+\zeta \mathrm{d}\sigma _g^{(1)}+\zeta ^2\mathrm{d}\sigma _g^{(2)}+\zeta ^3\mathrm{d}\sigma _g^{(3)},`$ (2.3)
For each contribution $`\mathrm{d}\sigma _g`$, it is a straightforward matter, given the values of $`\mathrm{d}\sigma _g`$ for the five input $`\zeta =0,\pm 1,\pm 2`$, to reconstruct $`\mathrm{d}\sigma _g^{(0,1,2,3)}`$. This is what we do. In fact for each set of parameters we automatically pick up all those diagrams that involve a dependence on the gauge parameter. The number of diagrams in this set depends on the parameter chosen. In some cases a very large number of diagrams is involved. For the process at hand this occurs with the parameter $`\stackrel{~}{\delta }`$ where about $`500`$ diagrams are involved in the check.

We then verify that the differential cross section is independent of $`\zeta `$

$`\mathrm{d}\sigma ={\displaystyle \underset{g}{}}\mathrm{d}\sigma _g={\displaystyle \underset{g}{}}\mathrm{d}\sigma _g^{(0)},`$ (2.4)
and therefore that

$`\mathrm{sum}_i={\displaystyle \frac{_g\mathrm{d}\sigma _g^{(i)}}{\mathrm{Max}_g\left(|\mathrm{d}\sigma _g^{(i)}|\right)}},i=1,2,3,`$ (2.5)
vanishes.

As seen from Table 1 agreement within $`20`$ to $`30`$ digits is observed. This agreement gets better if one gives the electron mass a higher value, say $`1`$GeV. The gauge parameter dependence check not only tests the various components of the input file (correct Feynman diagrams for example) but also the symbolic manipulation part and most important of all the correctness of all the reduction formulae and the proper implementation of all the $`N`$-point functions. This is quite useful when one deals with 5-point functions as is the case at hand. Talking of parametric integrals all tensor reductions are done following the standard procedure and then passing the scalar integrals to the FF package or to our own specially optimised routines when photon exchange is involved. The pentagon integrals are expressed in terms of boxes as is now standard, our procedure is outlined in the Appendix. We work in the on-shell renormalisation scheme closely following . Apart from masses and couplings, renormalisation is also carried for the fields. In particular we also require the residues of the renormalised propagators of all physical particles to be unity. As known, this procedure leads to a (very sharp) threshold singularity in the wave function of the Higgs at the thresholds corresponding to $`M_H=2M_W,2M_Z`$. Solutions to smooth this behaviour, like the inclusion of the finite width of the $`W`$ and $`Z`$, do exist but we have not implemented them yet in the present version of GRACE-loop . Therefore when scanning over $`M_H`$ it is sufficient to avoid these regions within $`1`$GeV around the thresholds.

As a separate check on our implementation we have also computed $`HWW`$ and $`e^+e^{}ZH`$, after tuning our parameters we find excellent agreement with the literature.

## 3 Results

### 3.1 Full $`𝒪(\alpha )`$ results

The results we show here include all 3 neutrino species. We first discuss the full $`𝒪(\alpha )`$ which includes the hard bremsstrahlung part. The separation between the $`s`$-channel and the $`t`$-channel is done in the same way as with the separation done at tree-level. It should be noted that in the one-loop diagrams that contribute to $`e^+e^{}\nu _\mu \overline{\nu }_\mu H`$, and which we classify as $`s`$-channel, there are diagrams which can not be deduced from the one-loop corrections to the $`s`$-channel $`22`$ process $`e^+e^{}ZH`$. Graph 249 in Fig.1 is one example. The relative correction for all three contributions is defined as $`\delta _{𝒪(\alpha )}\frac{\sigma _{O(\alpha )}}{\sigma _{tree}}1`$ and will be referred to as the full one-loop $`𝒪(\alpha )`$ correction.

One important remark is that the overall correction in the $`s`$-channel and $`t`$-channel are quite different. For a light Higgs of mass $`150`$GeV the correction to the $`s`$-channel Higgs-strahlung contribution is positive for practically all centre of mass energies of the next linear collider, see Fig. 2. It rises rather sharply to reach about $`+20\%`$ for a centre-of mass energy of 1TeV, however as will be argued below, the bulk of these large corrections are due to virtual and real QED corrections. Moreover in regions where these corrections are large, the Higgs-strahlung contribution is rather small. On the other hand the total correction in the $`t`$-channel for a small Higgs mass of 150GeV is negative throughout the range $`\sqrt{s}=350`$GeV to 1TeV, and is almost constant past $`\sqrt{s}=500`$GeV, reaching about $`2\%`$. Combining these two contributions, we see that the full correction to the whole process also remains small for a small Higgs mass. In fact at $`500`$GeV it is almost at its lowest of about $`0.2\%`$. This is an accidental cancellation between the contributions of both the $`t`$-channel and the $`s`$-channel at this energy. We have also studied the Higgs mass dependence of the corrections. First of all note that all our results capture the sharp spikes for $`M_H=2M_W,2M_Z`$, the top threshold is also visible when we plot the relative corrections. For $`500`$GeV the $`s`$-channel $`𝒪(\alpha )`$ correction decreases for higher Higgs masses, eventually turning negative with a value $`6\%`$ for $`M_H=350`$GeV. It then drops rather sharply to reach as much as $`30\%`$ at the $`ZH`$ threshold, past which this cross section is completely negligible. This behaviour is largely due to QED corrections and is driven by the kinematics of the two-body $`e^+e^{}ZH`$. At $`500`$GeV the correction in the $`t`$-channel contribution remains negative for all Higgs masses that we considered, i.e. in the range $`115450`$GeV. It drops steadily from about $`2\%`$ for $`M_H=115`$GeV to about $`10\%`$ at $`350`$GeV close to the top pair threshold. It then increases up to the $`ZH`$ threshold before dropping sharply around the $`ZH`$ threshold. Most of the large corrections are due to QED corrections.

### 3.2 Extraction the QED corrections

As is known large QED corrections require a higher order treatment. In order to quantify the effects of the genuine weak corrections, one could try subtract these QED corrections. This can be done rather easily for the $`s`$-channel contribution, where the correction can be readily extracted from the electromagnetic correction to the $`eeZ`$ vertex and the soft-photon bremsstrahlung part. Indeed our computation produces at an intermediate stage the result including the soft bremsstrahlung correction, that is before the inclusion of hard photons. The cut on the photon energy, $`k_c`$, has been taken sufficiently small, $`k_c=0.1`$GeV. These corrections without hard bremsstrahlung include thus the QED virtual and soft bremsstrahlung (which depend on $`k_c`$) as well as the genuine weak correction to the process. For this $`s`$-channel process the latter QED corrections are given by the universal soft photon factor that leads to a relative correction

$`\delta _{V+S}^{QED}={\displaystyle \frac{2\alpha }{\pi }}\left((L_e1)\mathrm{ln}{\displaystyle \frac{k_c}{E_b}}+{\displaystyle \frac{3}{4}}L_e+{\displaystyle \frac{\pi ^2}{6}}1\right),L_e=\mathrm{ln}(s/m_e^2).`$ (3.6)
where $`m_e`$ is the electron mass and $`E_b`$ the beam energy $`s=4E_b^2`$. Subtracting this contribution from our $`k_c`$ dependent (numerical) result reproduces the genuine weak correction, $`\delta _{W,schannel}`$.

To quantify in an unambiguous way the effect of the weak correction in the $`t`$-channel we have also subtracted this universal factor $`\delta _{V+S}^{QED}`$ from the full $`𝒪(\alpha )`$ correction. It can be shown that the leading (infrared and collinear) contributions are given by the universal factor, Eq. 3.6. This procedure also paves the way to a resummation of these large QED factors for the full process which is conducive to a Monte-Carlo implementation as could be done for instance through a QED parton shower. This will be treated elsewhere. Coming back to the weak corrections, we will denote by $`\delta _{W,tchannel}`$ and $`\delta _{W,total}`$ the weak correction for the $`t`$-channel and the full process based on the subtraction of the universal QED factor in Eq. 3.6.

For a light Higgs mass ($`M_H=150`$GeV), the weak correction for the $`s`$ and $`t`$-channels have a different behaviour as the energy increases, see Fig. 3. The former varies from about $`6\%`$ at $`\sqrt{s}=300`$GeV to $`2.5\%`$ at $`\sqrt{s}=1`$TeV. Past $`400`$GeV where it dominates, the weak correction to the $`t`$-channel varies rather slowly from $`7\%`$ at $`\sqrt{s}=400`$GeV to about $`5\%`$ at $`\sqrt{s}=1`$TeV. The dependence of the weak corrections on the Higgs mass for a moderate centre-of-mass energy, $`\sqrt{s}=500`$GeV, reveals that up to the $`ZH`$ threshold these corrections increase with the Higgs mass (most probably due to $`M_H^2`$ terms from the Higgs self-coupling as in $`e^+e^{}ZH`$), apart from the clearly visible spikes at the $`W,Z`$ and the top thresholds. Apart from the drop in the $`t`$-channel contribution around the $`ZH`$ threshold, the weak correction in the $`t`$-channel picks up again and as expected merges with the correction to the full process.

### 3.3 Expressing the weak corrections in terms of $`G_\mu `$

Expressing the corrections in the $`G_\mu `$ scheme or in other words had we expressed our tree-level results in terms of $`G_\mu ,`$ thus subtracting some universal weak corrections (essentially fermionic contributions) affecting two-point functions, we can have a quantitative measure of the non-universal weak radiative correction specific to this process. We thus define for the $`s`$-channel and $`t`$-channel contributions, these weak corrections as $`\delta _W^G=\delta _W3\mathrm{\Delta }r`$. Let us briefly summarise our findings for $`M_H=150`$GeV where with our input $`\mathrm{\Delta }r`$ contributes about $`3\%`$ (the leading Higgs mass dependence in $`\mathrm{\Delta }r`$ is logarithmic). For the full contribution with all three neutrinos we find $`\delta _W^G`$ to be slowly varying (with exactly the same “slope” as $`\delta _{W,total}`$ in Fig. 3), from about $`2\%`$ to about $`4\%`$ in the energy range from $`\sqrt{s}=300`$GeV to $`\sqrt{s}=1`$TeV. These genuine weak corrections remain therefore well contained in the full process, but in view of the precision of the $`e^+e^{}`$machine they must be taken into account. Applied to the $`s`$-channel with $`M_H=150`$GeV, the corrections with $`G_\mu `$ as an input, are moderate for energies up to $`400`$GeV but they quickly decrease below about $`12\%`$ at $`1`$TeV. Such behaviour had been observed in $`e^+e^{}ZH`$. This is another manifestation of the failure of the $`G_\mu `$ scheme to properly describe the weak corrections for such processes at high energies. For example it is known that in $`e^+e^{}ZH`$ the contribution of boxes is important. We do not attempt in this letter to make a thorough investigation of the different loop contributions to $`e^+e^{}\nu \overline{\nu }H`$, for example the fermionic and bosonic contributions. We leave this to a further study. This could be interesting in order to devise reliable approximations based on a small subset compared to the large number of contributions for such a complex process. For example, very recently, the fermionic contributions and especially the effect of the third generation have been investigated in  and  with differing results. It could be interesting to see how well these contributions can reproduce the full result. We also do not report here on how the distributions in the Higgs variables are affected by the radiative corrections. We have briefly discussed this in a previous note and leave the full discussion for a forthcoming paper.

## 4 Conclusions

We have calculated with the help of GRACE-loop the full radiative corrections including hard photon radiation to the important Higgs discovery channel at a future high energy $`e^+e^{}`$machine, $`e^+e^{}\nu \overline{\nu }H`$. Apart from the usual checks on the ultraviolet and infra-red finiteness of the result, we have performed tests on the gauge parameter independence of the results. To this end we have relied on a generalised non-linear gauge fixing condition where one has control over five independent gauge parameters. For a light Higgs of mass $`150`$GeV for energies ranging from $`300`$GeV to $`1`$TeV we find a modest total $`𝒪(\alpha )`$ correction which is within $`\pm 2\%`$, being negligible at $`500`$GeV ($`2`$ per-mil). We have also studied the Higgs mass dependence at $`\sqrt{s}=500`$GeV. For example with $`M_H=350`$GeV we find a larger $`𝒪(\alpha )`$ negative correction of about $`10\%`$. In order to quantify the weak correction we have subtracted the universal QED virtual and bremsstrahlung corrections. In the energy range $`\sqrt{s}=300`$GeV to $`\sqrt{s}=1`$TeV we find, for $`M_H=150`$GeV, that for the full process the correction ranges from $`+7\%`$ to $`+5\%`$ when the tree-level is expressed in terms of $`\alpha `$. Further investigations and details on this important process are left to a forthcoming publication.

Acknowledgment This work is part of a collaboration between the GRACE project in the Minami-Tateya group and LAPTH. D. Perret-Gallix and Y. Kurihara deserve special thanks for their contribution. This work was supported in part by Japan Society for Promotion of Science under the Grant-in-Aid for scientific Research B(No. 14340081) and PICS 397 of the French National Centre for Scientific Research.

## A Appendix

Five point functions are calculated as linear combinations of four point functions. Our method is based on an identity suitable for the Feynman parameter integration, which is similar to the one described in.

A five point function is expressed as

$`I_5`$ $`=`$ $`{\displaystyle \frac{d^4l}{(2\pi )^4}\frac{N(l)}{D_0D_1D_2D_3D_4}},`$ (A. 1)
where $`l`$ is the loop momentum and $`N(l)`$ is a polynomial of $`l^2`$ and inner products of $`l`$ with other four-vectors. The denominators of propagators are defined as

$`D_0`$ $`=`$ $`l^2m_0^2=l^2+X_0,`$
$`D_i`$ $`=`$ $`(l+r_i)^2m_i^2=l^2+2l.r_i+X_i,i=1,\mathrm{},4.`$ (A. 2)
We take a set $`r_i`$ ($`i=1,\mathrm{},4`$) of linearly independent momenta. The latter form a basis for vectors in 4-dimensional space. Therefore with the Gram matrix $`A_{ij}=r_i.r_j`$ one has the following identity

$`g^{\mu \nu }={\displaystyle \underset{i,j=1}{\overset{4}{}}}r_i^\mu A_{ij}^1r_j^\nu l^2={\displaystyle \underset{i,j=1}{\overset{4}{}}}l.r_iA_{ij}^1l.r_j.`$ (A. 3)
Combining this identity with Eq.(A. 2) we obtain

$`1`$ $`=`$ $`{\displaystyle \underset{\alpha =0}{\overset{4}{}}}[a_\alpha +{\displaystyle \underset{i=1}{\overset{4}{}}}l.r_ib_{\alpha ,i}]D_\alpha ,`$ (A. 4)
where

$`a_i`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }}}{\displaystyle \underset{j=1}{\overset{4}{}}}A_{ij}^1(X_jX_0),`$ (A. 5)
$`a_0`$ $`=`$ $`{\displaystyle \frac{4}{\mathrm{\Delta }}}{\displaystyle \underset{i=1}{\overset{4}{}}}a_i,`$ (A. 6)
$`b_{i,k}`$ $`=`$ $`{\displaystyle \frac{2}{\mathrm{\Delta }}}A_{ik}^1,`$ (A. 7)
$`b_{0,k}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{4}{}}}b_{i,k},`$ (A. 8)
$`\mathrm{\Delta }`$ $`=`$ $`4X_0+{\displaystyle \underset{i,j}{}}(X_iX_0)A_{ij}^1(X_jX_0).`$ (A. 9)
This immediately shows that the five point tensor integral can be reduced to $`5`$ box integrals.

Now we introduce the Feynman parameters. It is easy to see that

$`I_5={\displaystyle }{\displaystyle \frac{d^4l}{(2\pi )^4}}{\displaystyle }{\displaystyle \underset{\lambda =0}{\overset{4}{}}}dx_\lambda \delta (1{\displaystyle \underset{\beta =0}{\overset{4}{}}}x_\beta ){\displaystyle \frac{N(l)}{D^4}}{\displaystyle \underset{\alpha }{}}(a_\alpha +{\displaystyle \underset{i=1}{\overset{4}{}}}l.r_ib_{\alpha ,i})\delta (x_\alpha ).`$ (A. 10)
Making a shift in the loop momentum , $`lltt=_{i=1}^4x_ir_i`$, so as to eliminate linear terms in the loop momentum in $`D`$, we obtain our reduction formula

$`I_5`$ $`=`$ $`{\displaystyle }{\displaystyle \frac{d^4l}{(2\pi )^4}}{\displaystyle }{\displaystyle \underset{\lambda =0}{\overset{4}{}}}dx_\lambda \delta (1{\displaystyle \underset{\beta =0}{\overset{4}{}}}x_\beta ){\displaystyle \frac{N(lt)}{D^4}}{\displaystyle \underset{\alpha =0}{\overset{4}{}}}(\overline{a}_\alpha +{\displaystyle \underset{i=1}{\overset{4}{}}}l.r_ib_{\alpha ,i})\delta (x_\alpha )`$
$`\overline{a}_0`$ $`=`$ $`a_02/\mathrm{\Delta },\overline{a}_i=a_i,D=l^2+{\displaystyle \underset{\alpha _=0}{\overset{4}{}}}X_\alpha x_\alpha {\displaystyle \underset{i,j=1}{\overset{4}{}}}A_{ij}x_ix_j.`$ (A. 11)
For the scalar pentagon, $`N(l)=1`$, only $`\overline{a}_\alpha `$ in the previous equation contributes.


format: text

text:
# Contact Pair Dynamics During Folding of a Model Globular Protein, Hp-36

## Abstract

The dynamics of contact pair formation between various hydrophobic residues during folding of a model protein Hp-36 is investigated by Brownian dynamics simulation. Hydropathy scale and non-local helix propensity of amino acids are used to model the complex interaction potential. The resulting structure of the model protein mimics the native state of the real protein with a $`RMSD`$ of 4.5 Å. A contact pair distance time correlation function (CPCF), $`C_P^{ij}(t)`$, is introduced which shows multistage decay, including a slow late stage dynamics for a few specific pairs. These pairs determine the long time folding rate. Dynamics can be correlated with the landscape, relative contact order and topological contact.

The dynamics of folding of an extended protein chain at high temperature (or high urea concentration) to its unique folded state at low temperature (or low urea concentration) is a highly complex problem with many interesting aspects. Recent experimental, theoretical and computer simulations studies $`[17]`$ have unearthed and explained many fascinating aspects of folding, although still many others remain to be explored. The paradigm of landscape (with the idea of folding funnel) has provided new insight into the problem pgscience ; bw . Experimental data on the rate of folding of a large number of small proteins have suggested a close relation between the the relative contact order and the rate of folding baker . The relative contact order denotes the average sequence distance between the hydrophobic pair contacts and is defined as plaxo ,

$$RCO=\frac{\underset{i,j}{}(s_js_i)}{LN_c}$$
(1)
where ($`i,j)`$ are the specific hydrophobic pair contacts, $`N_c`$ is the number of contacts while $`L`$ is total number of hydrophobic amino acids present in the protein. $`s_j`$ and $`s_i`$ are the sequence number along the contour of the chain. The rate of folding was found to decrease nearly exponentially with $`RCO`$. The dynamics of such non-local contact formation holds the key to the understanding of dynamics of folding. However, this aspect has remained largely unexplored.

An attractive way to explore pair dynamics of hydrophobic contact is via the technique of fluorescence resonance energy transfer (FRET). In FRET, one measures the time dependence of energy transfer from a chosen donor fluorophore to a chosen acceptor. The rate of transfer may be due to dipolar interactions and the rate of transfer is given by the well known Förster expression forster ,

$$k_f=k_{rad}\left(\frac{R_F}{R}\right)$$
(2)
where $`k_{rad}`$ is the radiative rate and $`R_F`$ is the Förster radius. By suitably choosing donor-acceptor pair, $`R_F`$ can be varied over a wide range. This allows the study of the dynamics of pair separation, essential to understand protein folding gray . $`k_{rad}`$ is typically less than (but of the order of) $`10^9`$ $`sec^1`$. Thus Förster transfer provides us with a sufficiently fast camera to take snapshots of the dynamics of contact pair formation.

In this study, we have studied contact formation by Brownian dynamics simulations of a model protein Hp-36 which is one of the smallest protein that folds autonomously to a stable compact ordered structure, with a large helix content mcknight . Hp-36 is a subdomain of chicken villin which is implicated in the formation of microvilli in the absorbtive epithelium of the gut and the proximal tube of the kidney bretscher . All atom simulation study on this Hp-36 have revealed at least two pathways of folding kollman . Earlier, several studies of Hp-36 were presented using Monte Carlo technique hansmannprl  and Brownian dynamics srini .

The model studied here is constructed by taking two atoms for a particular amino acid. The smaller atom represents the backbone $`C_\alpha `$ atom of real protein while the bigger atom mimics the whole side chain residue. Construction of the model protein has been described in detail elsewherearnjcp . Similar types of model (with more rigorous force field) have been introduced by Scheraga et al. recently scheraga . The total potential energy function of the model protein $`V_{Total}`$ is written as,

$$V_{Total}=V_B+V_\theta +V_T+V_{LJ}+V_{helix}$$
(3)
where $`V_B`$ and $`V_\theta `$ are the potential contributions due to vibration of bonds and bending motions of the bond angles. Standard harmonic potential is assumed for the above two potentials with spring constants 43.0 kJ mole-1Å-2 and 8.6 kJ mole-1Å-2 for the bonds between backbone atoms and bonds joining side residues with the backbone atoms, respectively. In case of the bending potential, spring constant is taken to be 10.0 kJ mol-1 rad-2. $`V_T`$ ($`=ϵ_T_\varphi (1/2)[1+cos(3\varphi )]`$) is taken as torsional potentials for the rotations of the bonds. $`ϵ_T=`$ 1 kJ mol-1.

The nonbonding potential $`V_{LJ}`$ is the sum of the pair interactions between the atoms and is given by,

$`V_{LJ}=4{\displaystyle \underset{i,j}{}}ϵ_{ij}\left[\left({\displaystyle \frac{\sigma _{ij}}{r_{ij}}}\right)^{12}\left({\displaystyle \frac{\sigma _{ij}}{r_{ij}}}\right)^6\right]`$ (4)
where $`r_{ij}`$ and $`ϵ_{ij}`$ are the distance and interaction between the $`i`$-th and $`j`$-th atom. $`\sigma _{ij}=\frac{1}{2}(\sigma _{ii}+\sigma _{jj})`$ and $`ϵ_{ij}=\sqrt{ϵ_{ii}ϵ_{jj}}`$. Sizes and interactions are taken to be the same (1.8 Å and 0.05 kJ mol-1, respectively) for all the backbone atoms as they represent the $`C_\alpha `$ atoms in case of real proteins. Side residues, on the other hand, carry the characteristics of a particular amino acid. Different sizes of the side residues are taken from the values given by Levitt levwar . Interactions of the side residues are obtained from the hydrophobicities of the amino acids. We construct effective potential guided by the well-known statistical mechanical relation between potential of mean force and radial distribution function as $`V_{ij}^{eff}=k_BTlng_{ij}(r)`$ hansen . Strong correlation among the hydrophobic groups (absent among the hydrophilic amino acids) implies that the hydrophobic amino acids should have stronger effective interaction than the hydrophilic groups. So the interaction parameters of the side residues can be mapped from the hydropathy scale kyte  by using a linear equation as given below,

$$ϵ_{ii}=ϵ_{min}+(ϵ_{max}ϵ_{min})\left(\frac{H_{ii}H_{min}}{H_{max}H_{min}}\right)$$
(5)
where, $`ϵ_{ii}`$ is the interaction parameter of the $`i`$th amino acid with itself. $`ϵ_{min}`$(=0.2 kJ mol-1) and $`ϵ_{max}`$(=11.0 kJ mol-1) are the minimum and maximum value of the interaction strength chosen for the most hydrophilic(arginine) and most hydrophobic(isoleucine) amino acids, respectively. $`H_{ii}`$ is the hydropathy index of $`i`$th amino acid given by Kyte and Doolittle kyte  and $`H_{min}`$(=-4.5) and $`H_{max}`$(=4.5) are the minimum and maximum hydropathy index among all the amino acids. Further details are available in Ref. $`16`$.

An important part of secondary structure of the real protein is the formation of $`\alpha `$ helix. In the absence of hydrogen bonding, we introduce the following effective potential among the backbone atoms to mimic the helix formation along the chain of residues that show high helix propensity,

$$V_{helix}=\underset{i=3}{\overset{N3}{}}\left[\frac{1}{2}K_i^{13}(r_{i,i+2}r_h)^2+\frac{1}{2}K_i^{14}(r_{i,i+3}r_h)^2\right]$$
(6)
where $`r_{i,i+2}`$ and $`r_{i,i+2}`$ are the distances of $`i`$th atom with $`i+2`$ and $`i+3`$ th atoms, respectively. $`r_h`$ is the equilibrium distance and is taken as 5.5 Å, motivated by the observation that the distance of $`r_i`$ with $`r_{i+2}`$ and $`r_{i+3}`$ is nearly constant at 5.5 Å in an $`\alpha `$ helix. The summation excludes the first and last three amino acids as there is less helix formation observed in the ends of the protein chain zimm . The force constant for the above harmonic potential is mapped from the helix propensities $`Hp_i`$ taken from Scholtz et al. pace , $`𝒦_i=𝒦_{alanine}Hp_i\times (𝒦_{alanine}𝒦_{glycine})`$. $`𝒦_{alanine}`$ and $`𝒦_{glycine}`$ are the force constants for alanine and glycine, 17.2 and 0.0 kJ mol-1, respectively. Next, the influence of the neighboring amino acids for the formation of helix has been considered by taking an average of the spring constants as $`K_i^{13}=\frac{1}{3}[𝒦_i+𝒦_{i+1}+𝒦_{i+2}]`$ and $`K_i^{14}=\frac{1}{4}[𝒦_i+𝒦_{i+1}+𝒦_{i+2}+𝒦_{i+3}]`$, with the condition that $`K_i^{13},K_i^{14}0`$ as the force constant must remain positive. The above formulation of helix potential is motivated by the work of Chou and Fasman about the prediction of helix formation that the neighbors of a particular amino acid should be considered rather than its own helix propensity chou .

Figure $`1(a)`$ shows the best folded structure obtained in our simulations while figure $`1(b)`$ shows the real one. The $`RMSD`$ calculated over backbone with the real native Hp-36 is 4.5 Å which is reasonable for such a simplified model. The initial configuration of the model protein was generated by configurational bias Monte Carlo technique frenkel . Atoms attached to a single branch point were generated simultaneously. Then the initial configuration was subjected to Brownian dynamics simulation for the study the folding. Time evolution of the model protein was carried out according to the motion of each bead as below,

$$𝐫_i(t+\mathrm{\Delta }t)=𝐫_i(t)+\frac{D_i}{k_BT}𝐅_i(t)\mathrm{\Delta }t+\mathrm{\Delta }𝐫_i^G$$
(7)
where each component of $`\mathrm{\Delta }𝐫_{i\alpha }^G`$ is taken from a Gaussian distribution with mean zero and variance $`(\mathrm{\Delta }r_{i\alpha }^G)^2=2D\mathrm{\Delta }t`$ mccammon ; hansen . $`𝐫_i(t)`$ is the position of the $`i`$th atom at time t and the systematic force on $`i`$th atom at time t is $`𝐅_i(t)`$. The time step $`\mathrm{\Delta }t`$ is taken as 0.001. $`D_i`$ is the diffusion coefficient of the $`i`$-th particle calculated from the Stokes-Einstein relation $`D_i=\frac{k_BT}{6\pi \eta R_i}`$. $`R_i`$ is the radius of the $`i`$-th atom and $`\eta `$ is the viscosity of the solvent. $`k_B`$ and $`T`$ are the Boltzmann constant and temperature, respectively. Simulations have been carried out for $`𝒩`$ number of different initial configurations, where $`𝒩=584`$.

The potential energy $`E_N`$, the radius of gyration $`R_\gamma `$ and the $`RMSD`$ (calculated over backbone from the real native structure of Hp-36 obtained from protein data bank pdb ; mcknight ) all exhibit the multistage dynamics as shown in figure $`2(a)`$ in the same time axis. There is an initial sharp hydrophobic collapse followed by a slower decay is observed till 500$`\tau `$ ($``$ 200 ns). A long plateau follows in the final stage which exists for a very long time ( 2000$`\tau `$ 1$`\mu `$s) before the protein reaches its final lowest energy state at around 2400$`\tau `$. Energy values of the corresponding inherent structures obtained by conjugate gradient is shown by the symbols which depicts the decreasing local minima attained by the system until it reaches the final folded state. Total number of hydrophobic topological contact $`N_{topo}`$ and relative contact order $`RCO`$ are plotted in figure $`2(b)`$. $`N_{topo}`$ is defined to be formed if two hydrophobic side chain residues come within a distance of 8.5 Å. $`RCO`$ is calculated from Eq. $`1`$. Both $`N_{topo}`$ and $`RCO`$ show a multistage increase signifying the participation of nonlocal contacts.

Folding can be probed microscopically by monitoring the dynamics of separation between different amino acid pairs. The widely different time scales of movement of all the different pairs together give rise to an overall picture of the dynamics of folding which is reflected in the macroscopic quantities. The effective dynamics of pair separation can be described by introducing a new pair correlation function defined below srini ,

$$C_P^{ij}(t)=\frac{d^{ij}(t)d^{ij}(\mathrm{})}{d^{ij}(0)d^{ij}(\mathrm{})}$$
(8)
where, $`d^{ij}(t)=|𝐫_i(t)𝐫_j(t)|`$. $`𝐫_i`$ and $`𝐫_j`$ are the positions of the $`i`$-th and $`j`$-th atom, respectively. Figure $`3`$ shows the $`C_P^{ij}(t)`$ of the $`9`$th side residue with many other hydrophobic side residues. The three stages of the folding process are reflected by mainly three different dynamical behavior seen amongst the side residues. Side residues closed to the tagged one collapse very fast. Some show an initial shoulder and only a few show the plateau in the long time decay that correlates with the similar plateau observed in case of other dynamical quantities. The final decrease in energy is observed when all the different contact pair correlation functions decay at around 2400$`\tau `$.

We have calculated survival probability $`S_P(t)`$ in FRET using Förster energy transfer rate from Eq. $`2`$ for the $`935`$ pair along the Brownian dynamics trajectory leading to the most stable structure. The Förster radius $`R_F`$ is taken as 10 Å. Figure $`4`$ shows initial very slow decrease of $`S_P(t)`$ to be followed by a sudden drop at around 2400$`\tau `$ as seen in case of the different dynamical properties discussed above. $`S_P(t)`$ is found to be relatively insensitive to $`k_{rad}`$.

The resemblance with the real native state of the protein and the dynamics of $`RMSD`$ show the validity of the model used for this protein. Contact pair dynamics and the time evolution of energy, radius of gyration, relative contact order formation etc. brings out the rich and diverse dynamics of protein folding. The initial ultrafast hydrophobic collapse signify that the upper part of the funnel is steep – followed by a change in slope. Rate determining step, however, arises from the final stage of folding on a very flat and rugged underlying landscape marked by the large conformational entropy barrier with little energy change bob . This entropic bottleneck arises from the necessity to form long range hydrophobic contacts, as envisaged by Dill and Wolynes. The atoms mimicking the whole side chain of the real protein play a very important role for structural and dynamical aspects in this study. Moreover, the new contact pair correlation function and FRET probes the folding events in minute detail.

A. Mukherjee thanks Ashwin, Kausik and Prasanth for technical discussions. Authors thank CSIR, New Delhi, India and DST, India for financial support.


text:
# Particle acceleration in ultra-relativistic oblique shock waves

## 1 Introduction

Detection of $`10^610^{20}`$eV nuclei within the solar system, implying their presence through much of the observable universe, has stimulated the development of the diffusive shock acceleration model, whereby particles are repeatedly accelerated in multiple crossings of a shock interface, due to collisions with upstream and downstream magnetic scattering centres. The pioneering work of the late 70s (Krymsky, 1977; Bell, 1978a,b; Axford et al., 1978; Blandford and Ostriker, 1978) established the basic mechanism of particle diffusive acceleration in non-relativistic flows. However, since then, extension of the model to oblique shocks, relativistic plasma flows and the incorporation of non-linear effects, have rendered the emergence of a consensus physical picture more difficult. At oblique shock fronts the motion of the particles is more complicated, compared to a parallel shock configuration, as they may either be transmitted or reflected by the shock surface. Apart from the diffusive acceleration mechanism, an oblique shock could also accelerate particles in the absence of scattering (Armstrong and Decker, 1979), and this mechanism is called ’shock drift’ acceleration. Provided that the motion of the particles is diffusive in the oblique shock configurations it has been shown, for an analytical non-relativistic flow approach, that the same result as in parallel shocks, connects the power-law index of accelerated particles with the compression ratio and the acceleration rate with the diffusion coefficients in the shock normal direction (e.g. Axford et al., 1978; Bell, 1978a; Drury, 1983).
Also, following energy gain in the shock frame on shock crossing via the drift approach, is equivalent to following the gain via Lorentz transformations from the shock frame into an E=0 frame (defined by de Hoffmann and Teller, 1950) at the shock, and then back into the shock frame. On the other hand, the theory of oblique shocks considers features that are not involved in the parallel shock case, in particular the ability of the shock to reflect particles meeting a weaker to stronger field transition and the dependence of the acceleration time on the field-shock normal angle (e.g. Jokipii, 1982). Oblique shocks can be classified into sub-luminal and super-luminal. The difference between these two categories is that, in the sub-luminal case, it is possible to find a relativistic transformation to the frame of reference (de Hoffmann-Teller frame), in which the shock front is stationary and the electric field is zero (E=0) in both upstream and downstream regions. On the other hand, super-luminal shock fronts do not admit a transformation to such a frame, as they correspond to shock fronts in which the point of intersection of the front with a magnetic field line moves with a speed greater than $`c`$. Kirk and Heavens (1989) employed a numerical solution to the transport equation (in the test particle limit) and found spectral flattening at relativistic speeds in inclined shocks as well as an increased anisotropy both in fluid frames and the E=0 frame. These results were confirmed by Monte Carlo simulations (Ostrowski, 1991 and Lieu et al., 1994). The former author investigated the gyromotion across the shock to check the region of applicability of magnetic moment conservation in the E=0 frame at the shock while, the later authors simply relied on this assumption. Baring et al. (1994) noticed ’bumps’ in the accelerated spectra just upstream of the shock in the non-relativistic oblique case. The reason a Monte Carlo approach became favoured, lies in the large anisotropies and gradients which develop around the shock at relativistic flow speeds. These render direct solution of the transport equation difficult, because the Boltzmann equation must be applied to a fine-grained bundle of particle trajectories in phase space and the diffusion approximation, employing a simple spatial diffusion tensor is inadequate (e.g. Axford, 1981).
All the work mentioned so far is sub-luminal, but Bednarz and Ostrowski (1998) used a small-angle scattering model to include super-luminal situations with field-shock normal angles, $`\psi \pi /2`$ and flow $`\mathrm{\Gamma }243`$ and also a simulation of varying amounts of cross-field diffusion. These authors found spectral indices which could be very large at intermediate $`\mathrm{\Gamma }`$ values with small $`K_{}`$, but at high $`\mathrm{\Gamma }`$, all results appeared to tend towards a differential exponent of -2.2.
Earlier discussion of an extreme example of a super-luminal case, the pulsar wind shock acceleration scheme for the termination shock of an isolated pulsar magnetosphere where the $`\mathrm{\Gamma }10^4`$ polar wind meets the slowly expanding nebula envelope, has emphasized the apparent difficulty in achieving multiple shock crossings of accelerated particles. This can be seen if we remember that for super-luminal shocks, a transformation to a frame with B parallel to the shock plane is possible (Hudson, 1965) so in the absence of scattering close to the shock, particles are swept at an angle to the shock normal into the shock, shock drift and then exit downstream with no return possible. Begelman and Kirk (1990) note that only shock compression of the upstream population results in the low turbulence case. However, Lucek and Bell (1994) find special trajectories if the field is ordered but not uniform, allowing energy increases by factors of the order of hundreds per crossing, which coupled with particle motion in the electromagnetic field of the nebula can render isolated pulsars accelerators of high energy cosmic rays. Gallant and Arons (1994) provide a rare example of a hybrid model of relativistic particle acceleration, applied to the Crab Nebula, using a fluid approach for an $`e^\pm `$ pair plasma together with a kinematic treatment of ions to compute the synchrotron spectrum.
Acceleration time under non-relativistic theory clearly depends on shock obliquity through variation of $`K_n`$, the diffusion coefficient along the shock normal, with $`\psi `$, where $`K_n=K_{||}cos^2\psi +K_{}sin^2\psi `$. $`K_{||},K_{}`$ are respectively the diffusion coefficients parallel and perpendicular to the field. The acceleration time is given by,

$$\tau =\frac{3}{V_1V_2}\left(\frac{K_{n,1}}{V_1}+\frac{K_{n,2}}{V_2}\right)$$
(1)
Here the suffixes ’1’ and ’2’ refer to upstream and downstream, and $`V`$ is plasma flow velocity. Because this equation (1) is clearly an inadequate measure of acceleration time in the highly relativistic case, we note it may also be written as $`\tau `$ lying in the range $`(48/3)(nr_g/c)cos^2\psi `$ for a parallel mean free path $`\lambda _{}=nr_g`$, a shock frame compression ratio of 4 and a downstream reduction in $`\lambda _{}`$ by a factor 4. Hence while it will be convenient to measure any ’speed-up’ of acceleration relative to equation (1), we may equally regard the time to move over the relativistic Larmor radius as the unit of acceleration time for comparison purposes. Quenby and Lieu (1989) and Ellison et al. (1990) employed Monte Carlo guiding centre approximation calculations to show that relativistic parallel shock flows with $`\mathrm{\Gamma }10`$ produced a ’speed-up’ of acceleration by factors $`310`$ when compared with equation (1) for $`\tau `$. The former authors argued for large angle scattering while the latter found essentially the same result with either large or small angle scattering. Lieu et al. (1994) showed the speed-up as a function of flow velocity was similar for inclined, sub-luminal shocks, provided the results are plotted against flow seen in the E=0 frame. Quenby and Drolias (1995) using an inclined shock model which included the expected shock structure due to cosmic ray pressure back reaction, found some speed-up still occurred. Bednarz and Ostrowski (1996) employ a Monte Carlo computation scheme to investigate acceleration times for shock speeds up to 0.9$`c`$ for the inclined case, with both small angle and large angle scattering and allowing for cross field motion of guiding centres. Speed-up is noticed especially at high shock speed and turbulence although, cross-field diffusion eventually reverses the correlation with turbulence level. Power law spectra only occur in the super-luminal case for isotropic scattering. Acceleration time scales, shorter than the upstream gyroperiod have been obtained.
In a companion work, Meli and Quenby (2002a), have studied parallel shock acceleration for flow $`\mathrm{\Gamma }10^3`$, for both small and large angle scattering models using the guiding centre Monte Carlo approach. They confirm a $`\mathrm{\Gamma }^2`$ energy enhancement for the first shock crossing cycle, first emphasized by Quenby and Lieu (1989,) and they find significant energy gains subsequently. Also, the structure in the accelerated particle spectrum at the shock, indicated in previous lower $`\mathrm{\Gamma }`$ simulations, became enhanced for large angle scattering, although not in the small angle case, rendering the idea of a simple power law output spectrum for photon production uncertain and they confirm the existence of the speed-up effect. It is the purpose of this work to extend the parallel shock studies of Meli and Quenby (2002a) to provide answers concerning the acceleration time scales, particle spectral shape at the source, angular distribution and energy gain per cycle of the accelerated particles for different shock inclinations using very high gamma flows. We consider both sub-luminal shocks (where an E=0 frame can be found) and super-luminal ones (where E=0 is impossible) and the helical particle trajectories are followed in the appropriate frames to determine the characteristics of the shock crossing. Streaming instabilities, drift waves and Rayleigh-Taylor like instabilities all can increase scattering around the shock and contribute to the back reaction of the relativistic gas, but these are not enough to destroy the basic relativistic beaming in flows, where $`Vc`$. Hence a Monte Carlo approach is employed. The speed-up for a true discontinuity in the test-particle strong scattering limit, is due to a $`\mathrm{\Gamma }^2`$ change in energy in one complete shock crossing cycle for shock $`\mathrm{\Gamma }`$, and this effect is the basic cause of the deviation of the spectral shape from the non-relativistic result. Our current interest, regarding a real astrophysical scenario, lies in relativistic flows in AGN jets where $`\mathrm{\Gamma }510`$ and in GRB fireballs where a $`\mathrm{\Gamma }10^210^3`$ is estimated.

## 2 Numerical Method

In these Monte Carlo codes we use both large angle and pitch angle scattering, assuming elastic scattering in the fluid frames and the important aim is to investigate the role of different scattering models in reference to spectral shape, acceleration time scales and angular distribution in different frames of reference with varying shock obliquity at very high gamma plasma flows, ranging from $`\mathrm{\Gamma }=1010^3`$. We consider oblique, sub-luminal and super-luminal shocks with $`15^o\psi 75^o`$ where $`\psi `$ is the angle between the shock normal and the magnetic field, seen in the shock frame. The relativistic frames of references we use during the calculations are, the local fluid frame (1,2), the normal shock frame (s) and the de Hoffman-Teller frame (HT). In general, flow into and out of the shock discontinuity is not along the shock normal, but a transformation is possible into the normal shock frame to render the flows along the normal (e.g. Begelman and Kirk, 1990) and for simplicity we assume such a transformation has already been made.

### 2.1 Sub-luminal case

The Monte Carlo codes for the simulations in the sub-luminal case are similar to the treatment we followed in Meli and Quenby (2002a) for the parallel shock simulations, apart from the fact that more complicated relativistic formulas and different frames of references are used (which make the runs of the codes considerably more CPU-time consuming). First, relativistic particles are injected upstream and are allowed to scatter in either fluid frame and across the shock discontinuity. We follow each particle trajectory across the shock according to the jump conditions and each particle leaves the system when it ’escapes’ far downstream at the spatial boundary or if it reaches a well defined maximum energy $`E_{max}`$, equal to $`10^{10}E_o`$. The spatial boundary is placed $`100\lambda `$ downstream in the large angle scatter case which proved an optimum number, but for small angle scattering, the possible dependence of the results on this parameter needed investigation (see the ’Results’ section). The compression ratio $`r=\frac{V_1}{V_2}`$ takes the values of 4 and 3 for comparison with previous work with varying relativistic Rankine-Hugoniot conditions and a similar version of Bednarz and Ostrowski (1998) splitting technique is used (see Meli and Quenby, 2002a) so that when an energy is reached such that only a few accelerated particles remain, each particle is replaced by $`N`$ of statistical weight $`1/N`$ so as to keep roughly the same number of particles being followed. A guiding centre approximation is used, where the particle trajectory is followed in two-dimensional, pitch angle and distance along B space. The mean free path is calculated in the respective fluid frames by the formula: $`\lambda =\lambda _{}p_{1,2}`$, assuming a momentum dependence to this mean free path for scattering along the field and related to the spatial diffusion coefficient, $`K`$, in the shock normal or $`x`$ direction by $`K_{}=\lambda v/3`$ and $`K=K_{}cos^2\psi `$. Thus we specifically neglect cross-field diffusion, under conditions to be discussed in the next paragraphs.
There are three conditions for the guiding centre approximation to be fulfilled relating to the ratio $`K_{}/K_{}`$, the probability of hitting the shock at the last scattering including the effect of field line wandering near the shock and the across field jumps each scatter and we mention these, employing entirely plasma frame variables. If $`K_{}/K_{}0.1`$, diffusion away from the shock is clearly predominantly along the mean field. Adapting quasi-linear theory for gyro-resonance scattering by a wave number, $`k`$, spectrum of fluctuations perpendicular to the mean field of power $`P(k)=Ak^n`$ we get (Jokipii, 1966),

$$K_{}=\frac{vB^2r_g^{2n}}{\pi (2n)(4n)A}$$
(2)
To match the usual assumption, $`\lambda _{}r_g`$, we need $`n`$=1. Then if $`k_{},k_1`$ are the lower and upper limits to the wave number spectrum, chosen to avoid divergence in fluctuation power,

$$\lambda _{}=\frac{r_g}{\pi }\frac{ln(k_1/k_{})}{\delta b^2}$$
(3)
where $`\delta b^2`$ is the power in transverse fluctuations normalized to the mean field. Perpendicular diffusion is likely to be predominantly due to field line wandering with (Jokipii and Parker, 1969),

$$K_{}=\frac{v}{4B^2}P_{}(k0)=\frac{v}{4k_{}}\frac{\delta b_{}^2}{ln(k_1/k_{})}$$
(4)
where $`P_{}(k0)`$ is the power in the parallel field at wave numbers near zero. Even if the fluctuations are isotropic parallel and perpendicular to the field and if $`k_{}^1=L_s/2\pi `$ where $`L_s`$ is the scale size of the acceleration region, with $`\delta b0.3`$ and $`k_1/k_{}10`$, since at the highest reasonable $`r_g`$, $`r_gL_s/10`$, we expect $`K_{}/K_{}0.05`$. Hence at high $`r_g`$ the guiding centre approximation for the mean field is reasonable.
At low $`r_g`$, perpendicular motion is entirely determined by the wandering of the field line bundle to which the particle is ’attached’ and it is the local value of $`\psi `$, within a few $`\lambda _{}`$ of the shock, which determines the acceleration details. However, a second criterion that the guiding centre approximation remains reasonable is that $`tan\psi \lambda _{}/r_g`$. This means physically that the gyrating particle does not on average scatter far enough off from the original ’bundle’ to intersect the shock front, as it scatters at a distance of one parallel mean free path from the front along the actual field line. The reason for the scattering, large angle or a reflection due to multiple small-angle scatterings, is not crucial. For the above parameters this criterion becomes $`\psi 84^{}`$. However, since the field itself ’wanders’ by an rms value $`18^{}`$, the guiding centre approximation is probably limited to field-shock normal angles $`66^{}`$, at all $`r_g`$.
As explained by Baring et al. (1995), each scatter induces a jump across the field which changes the distance to the shock by $`r_g`$, but this is unimportant in the numerical regime we consider where, $`\lambda _{}10r_g`$.
The probability that a particle will move a distance $`\mathrm{\Delta }s`$ along the field lines at pitch angle $`\theta `$ before it scatters, is given by the expression $`Prob(\mathrm{\Delta }s)exp(\mathrm{\Delta }s/\lambda )`$. The probability of finding a particle in a pitch angle interval $`d\mu `$ where $`\mu =cos\theta `$ is $`\mu `$, and the time between collisions is $`\upsilon |\mu |`$ for particle velocity $`\upsilon `$ and $`\mu `$ chosen randomly between $`1`$ and $`+1`$. Since cross field diffusion is neglected, $`\mathrm{\Delta }x=\mathrm{\Delta }scos\psi `$. To take into account increased downstream scattering, we take $`\lambda _{,1}=4\lambda _{,2}`$.
In the shock rest frame, the flow velocity ($`V_1`$,$`V_2`$) for upstream and downstream respectively, is parallel to the shock normal and the magnetic fields $`𝐁_1`$ and $`𝐁_2`$ are at an angle $`\psi _1`$ and $`\psi _2`$ to the shock normal respectively. We adopt a geometry with $`x`$ in the flow direction, positive downstream, $`𝐁_1,𝐁_2`$ in the $`xy`$ plane and directed in the negative $`x`$ and $`y`$ directions and only $`𝐄_z`$ finite and in the positive $`z`$ direction. A relativistic transformation is performed to the local fluid frame each time the particle scatters across the shock.
Now, for the transformations between the shock frame and the de Hoffmann-Teller frame, we need to boost by a $`V_{HT}`$ speed along the shock frame where, $`V_{HT}`$ is equal to $`V_1tan\psi _1`$. The transformation to the de Hoffmann-Teller frame though, is only possible if $`V_{HT}`$ is less or equal to the speed of light. This means that $`tan\psi _11`$ (sub-luminal case). The particles ($`10^5`$) of a weight equal to 1.0, are injected far upstream at a constant energy of high gamma, which supposes that a pre-acceleration of the particles has already taken place, following for example the astrophysical scenario of Gallant and Achterberg (1999) and as discussed in more detail in Meli and Quenby (2002a). The particles are left to move towards the shock, along the way they collide with the scattering centers and consequently as they keep scattering between the upstream and downstream regions of the shock (its width is much smaller than the particle’s gyroradius) they gain a considerable amount of energy (see figures of the next section). For the case of highly relativistic flows, the definition of particle pitch angle diffusion is that the scattering is limited within an upstream frame angle $`0.1/\mathrm{\Gamma }`$ where, $`\mathrm{\Gamma }`$ is the upstream gamma, measured in the shock frame (Gallant and Achterberg, 1999) and investigated by Protheroe et al. (2002). This condition is also included in the current Monte Carlo for appropriate comparisons. The energy of the particle is kept constant at the scattering centers and only the direction of the velocity vector is randomized. We put forward a condition for large angle scattering to occur in Meli and Quenby (2002a), $`r_{g,1}\lambda _1/2\mathrm{\Gamma }_1^2`$ where the mean free path and gyroradius within the scatter centre, $`r_{g,1}`$, are measured in the upstream frame and will also consider this case as possible, but the reasoning, both for the sub-luminal and the super-luminal situations is modified. However, we note Achterberg et al. (2001) exclude large angle scatterings in an extensive investigation of a variety of models.
Provided the field directions encountered are reasonably isotropic in the shock frame, we are about to show $`tan\psi _1=\mathrm{\Gamma }_1^1tan\psi _s\mathrm{\Gamma }_1^1\psi _1`$ where ’1’ and ’s’ refer to the upstream and shock frames. Now, $`V_1=c(10.5\mathrm{\Gamma }_1^2)`$ while, a particle at velocity $`c`$ in the upstream frame starting along the shock normal, needs a deflection of $`\mathrm{\Gamma }_1^1`$ radians to be reduced to $`V_1`$ resolved in the shock normal direction and hence, be overtaken by the shock. This angle is the limit of gyration allowed in the upstream field perpendicular to the shock normal before some large angle scatter must take place. If $`t_u`$ is the time to the scatter and $`\omega _{}`$ is the perpendicular gyrofrequency, the angular deflection limit is $`\omega _{}t_u=\mathrm{\Gamma }_1^1`$ where, $`\omega _{}=u_{}/r_{a,,1}`$ and $`u_{}=c\psi _1`$. Then the allowed distance before large angle scatter is $`\lambda _{}=ct_ur_{a,,1}`$. This means the Larmor radius of the high field ’bullet’ scatter centre is, $`r_{g,1}r_{a,1}/\mathrm{\Gamma }_1^2`$ using the previous condition, which requires field strength regimes as discussed in Meli and Quenby (2002a) for possible GRB situations.
Between the ’i th’ and ’i+1 th’ scattering, we check for a shock crossing up to down via ($`c=1`$ units)

$$\mathrm{\Delta }x_s^{i+1}=\mathrm{\Gamma }_1(\mathrm{\Delta }x_i^i+V_1|\frac{\mathrm{\Delta }x_1^isec\psi _1}{\upsilon _1cos\theta _1}|)$$
(5)
$$\mathrm{\Delta }t_s^{i+1}=\mathrm{\Gamma }_1(|\frac{\mathrm{\Delta }x_1^isec\psi _1}{\upsilon _1cos\theta _1}|+V_1\mathrm{\Delta }x_1)$$
(6)
where $`\psi _1`$ is measured in the upstream frame. Defining $`𝐁_s`$ and $`𝐄_s`$ in the shock normal frame with compression ratio $`r`$, we arrive at upstream frame quantities; $`B_{1,x}=B_{s,x}`$, $`B_{1,y}=\mathrm{\Gamma }_1(B_{s,y}+V_1E_{s,z})`$, $`E_{1,z}^{^{}}=0`$ where $`E_{s,z}=V_1B_ssin\psi _{1,s}`$. Also $`tan\psi _1=\mathrm{\Gamma }_1^1tan\psi _s`$. Note the high relative velocity of the frames in the $`x`$ direction induces an electric field which via a back induction suppresses the $`y`$ component of B and swings the field lines towards the $`x`$ axis in the plasma frame.
Physically, for the occurrence of a shock crossing the conservation of the adiabatic invariant in the de Hoffman-Teller frame is used. However, to get into this frame we need all components of momentum and so need to assign a particle phase angle $`\varphi `$ at random. Then, in $`m=1`$ units, the momentum components for particle mass $`\gamma `$ are,

$$p_{1,x}=\gamma _1\upsilon _1cos\theta _1cos\psi _1\gamma _1\upsilon _1sin\theta _1cos\varphi _1sin\psi _1$$
(7)
$$p_{1,y}=\gamma _1\upsilon _1cos\theta _1sin\psi _1+\gamma _1\upsilon _1sin\theta _1cos\varphi _1cos\psi _1$$
(8)
$$p_{1,z}=\gamma _1\upsilon _1sin\theta _1sin\varphi _1$$
(9)
A two stage transform is the initiated, first to the shock frame,

$$\gamma _s=\mathrm{\Gamma }_1(\gamma _1+V_1p_{x,1})$$
(10)
$$p_{s,x}=\mathrm{\Gamma }_1(p_{x,1}+V_1\gamma _1)$$
(11)
$$p_{s,y}=p_{y,1}$$
(12)
$$p_{s,z}=p_{z,1}$$
(13)
A transformation is now made into the de Hoffmann-Teller frame with a boost along the negative $`y`$ axis with $`V_{HT}=V_1tan\psi _1`$. Then in this new frame, the $`B`$ components become $`B_{HT,x}=\mathrm{\Gamma }_{HT}(B_{s,x}V_{HT}E_{s,z})`$, $`B_{HT,y}=B_{s,y}`$ and hence $`tan\psi _{HT,1}=\mathrm{\Gamma }_{HT,1}tan\psi _1`$. Here the field lines are swept towards the transformation velocity in the $`y`$ direction under frame transform. Energy and the $`y`$ component of momentum are then transformed into the de Hoffmann-Teller frames according to, $`\gamma _{HT}=\mathrm{\Gamma }(\gamma _s+V_{HT}p_{s,y})`$ and $`p_{HT,y}=\mathrm{\Gamma }_{HT}(p_{s,y}+V_{HT}\gamma _s)`$ allowing the particle pitch angle to be obtained from $`cos\theta _{HT}=𝐩𝐁/pB`$. Then the particle that crosses the shock from upstream, is transmitted only if its pitch angle is less than the critical pitch angle:

$$\theta _c=arcsin\sqrt{\frac{𝐁_{HT,1}}{𝐁_{HT,2}}}$$
(14)
From the conservation of the first adiabatic invariant we can find the new pitch angle in the downstream frame and similar transformations allow the particle scattering to be followed in this frame.

### 2.2 Super-luminal case

For the super-luminal case a variation of the computational model is applied, compared to the one depicted above. For this investigation we employ once again a Monte Carlo scheme and consider the motion of a particle of momentum $`p`$ in the magnetic field B. As we have mentioned, in the super-luminal situation it is not possible to transform to a single frame where E=0 (de Hoffmann-Teller frame). Hence, for our investigation the most convenient frames of reference to use are the fluid frames (where still the electric field is zero) and the shock frame which we will use only as a frame to check whether upstream or downstream conditions apply. We inject $`8\times 10^5`$ particles in order to keep reasonable statistics throughout the simulations. We consider large angle scattering, which is calculated in the respective fluid rest frames, based on the possible situation previously discussed. Also Bednarz and Ostrowski (1996) find that only in the large angle scattering case does a power law develop and thus this is the most favourable case to further investigate. It will at least yield a limit to the acceleration a super-luminal shock may provide. Initially we inject the particles $`100\lambda `$ from the shock and we follow their guiding center in the upstream frame as in the sub-luminal case until after the appropriate transformation to the shock frame the particle reaches the shock at $`x_{sh}=0`$. At this juncture, there is no easy approximation to determine the probability of shock crossing or reflection. We change the model to following a helical trajectory, in the fluid frames upstream (’1’) or downstream (’2’) respectively, where the velocity coordinates of the particle are calculated in a three dimensional space as follows:

$$\upsilon _{x_1}=\upsilon _1cos\theta _1cos\psi _1\upsilon _1sin\theta _1cos\varphi _1sin\psi _1,$$
(15)
$$\upsilon _{y_1}=\upsilon _1cos\theta _1sin\psi _1+\upsilon _1sin\theta _1cos\varphi _1cos\psi _1$$
(16)
and,

$$\upsilon _{z_1}=\upsilon _1sin\theta _1sin\varphi _1$$
(17)
where $`\theta _1`$ is the pitch angle and $`\psi _1`$ is the angle between the magnetic field and the shock normal.

We follow the trajectory in time, using $`\varphi _1=\varphi _{}+\omega t`$, where $`\varphi _1(0,2\pi )`$ and $`t`$ is the time from detecting shock presence at $`x_{sh}`$, $`y_{sh}`$, $`z_{sh}`$ by using,

$$dx=x_{sh}+\upsilon _{x_1}\delta t$$
(18)
$$dy=y_{sh}+\upsilon _{y_1}\delta t$$
(19)
and,

$$dz=z_{sh}+\upsilon _{z_1}\delta t$$
(20)
assuming that $`\delta t=r_g/Hc`$, where $`r_g`$ is the Larmor radius, $`H100`$. The particle’s gyrofrequency $`\omega `$ is given by the relation:

$$\omega _1=\frac{e|𝐁_1|}{\gamma _1}$$
(21)
$`𝐁_1`$ is the magnetic field, $`\gamma _1`$ is the particle’s gamma and $`e`$ is its charge in gaussian units. For a matter of convenience though, the last relation is transformed in units of $`c/sec`$. All the above relations apply to the downstream case by only changing respectively the signs from 1 to 2, and all the calculations are performed in the upstream or downstream frames. Because of the peculiar properties of the helix we need to establish where a particle -in the upstream frame- of a particular $`\theta ,\varphi `$ first encounters the shock. To establish when this happens, we choose to go back a whole period,

$$T_1=\frac{2\pi }{\omega _1}$$
(22)
by reversing signs of the helix velocity coordinates and keep checking throughout the simulation to see if the particle’s trajectory encounters the shock again. The starting point for transforming to the downstream frame is the furthest upstream shock crossing. After making the suitable relativistic transformations to the downstream fluid frame by calculating the ($`x_2,y_2,z_2,t_2`$) and ($`\upsilon _{x_2},\upsilon _{y_2},\upsilon _{z_2}`$) coordinates, the momentum $`P_2`$ and the gamma $`\gamma _2`$ of the particle, we follow the trajectory of the particle for a whole downstream period,

$$T_2=\frac{2\pi }{\omega _2}$$
(23)
checking to see whether the particle meets the shock again, by transforming to the shock frame. If the particle meets the shock then the suitable transformations to the upstream frame are made again and we follow the particle’s trajectory as described above. If the particle never meets the shock, than its guiding center is followed, the same way as mentioned earlier for the upstream side after the injection and it is left to leave the system if it reaches a well defined $`E_{max}`$ momentum boundary or, a spatial boundary of $`100\lambda `$.

## 3 Results

In this section we present results of the simulations described above. We show similar figures regarding the acceleration time decrease, mean energy gain, spectral shapes and angular distributions for both sub- and super-luminal situations. First, the code runs have been tested in the non-relativistic and mild relativistic regime using $`V_10.1c0.6c`$ and are found to be in accordance with non-relativistic theory and mild relativistic theoretical predictions (e.g. Kirk and Schneider, 1987a,b). In the highly relativistic regime for the sub-luminal case, the Lorentz factors of the flow used are equal to 10, 50, 200, 500 and 990. It is presumed that the injected particles are already relativistic with gammas just higher than the flow gammas. The compression ratio, $`r`$, has values of 3 and 4. The value of $`r`$=3 is expected to mimic the physical condition of an ultra relativistic gas. The value $`r`$=4 corresponds to the extreme hydrodynamic limit, where the magnetic field pressure is unimportant, but the shock speed in the upstream frame tends to $`c`$. However we will see from the figures that the value 3 or 4 does not alter the results (possibly due to the mildly-relativistic nature of the downstream plasma).

Two sub-luminal inclination cases are studied, equal to $`15^{}`$ and $`35^{}`$. Particles are injected isotropically at 100$`\lambda `$ upstream in the shock frame and when they reach the spatial boundary placed at 100$`\lambda _{}`$ downstream they exit the simulation for the large angle case. We apply large angle and pitch angle scatters. For pitch angle scatter, the stability of the results against boundary position enabled a suitable downstream exit to be chosen (1000$`\lambda `$) where here, $`\lambda `$ is the distance between scatterings. The mean free path ($`\lambda `$), is chosen to have the value of $`\lambda 20r_g`$, close to the value that Quenby and Lieu (1989) used, based on a summary of interplanetary transport simulations referring to the only astrophysical plasma where we have any detailed knowledge of the magnetic field properties and performed by Moussas et al., 1992 (and references therein), but in most investigations it is only the ratio of up/down stream mean free paths and their ratio to the spatial boundary which matters. The momentum boundary is equal to $`E_o10^{14}`$, where $`E_o`$ is equal to the initial energy of the particle at the injection.

For oblique sub-luminal shocks, we show in figure 1 the ratio of computational acceleration time to the analytic non-relativistic time constant, versus the logarithm of the upstream flow gammas used. The measurements are all made in the de Hoffmann-Teller frame as in Lieu et al. (1994) where, it was found that this was a relevant frame to observe acceleration ’speed-up’. Then the non-relativistic acceleration time is,

$$\tau =\frac{3}{V_1cos\psi _1V_2cos\psi _2}\left(\frac{K_{n,1}}{V_1cos\psi _1}\frac{K_{n,2}}{V_2cos\psi _2}\right)$$
(24)
When considering pitch or small angle scattering, the value of $`\lambda _{||}`$ employed is that for $`\pi /2`$ scatter after random pitch change so that, $`\lambda _{}=(\pi ^2/4)(\lambda /\delta \theta ^2)`$ where $`\lambda `$ is the distance between $`\delta \theta `$ collisions. We observe a significant acceleration time decrease of a factor $``$ 5-10, which is comparable to that of the parallel shock ’speed up’ found in Ellison et al. (1990) and Meli and Quenby (2002a). Basically, this ’speed up’ effect observed is due to the $`\mathrm{\Gamma }_1`$ energy increase each shock crossing, $`\mathrm{\Gamma }_1^2`$ per cycle although it can be reduced by the very small pitch angles the particle may have in the small angle case as it crosses the shock upstream to downstream. In figures 2-4 the mean energy gain per crossing is presented for large angle scattering and pitch angle diffusion, with $`\psi _1`$ equal to $`15^{}`$, $`35^{}`$, $`r`$=4 and $`r`$=3 and $`\mathrm{\Gamma }=`$10, 200 and 990. Particle energy is measured in $`\gamma `$ units and so the plots are independent of mass and sign of charge. The relative independence on $`r`$ but the similarity in the trend from $`\mathrm{\Gamma }^2`$ to $`\mathrm{\Gamma }`$ in energy dependence and finally ’saturation’ found in the parallel case (Meli and Quenby, 2002a) are seen. The gain per crossing falls off noticeably faster as $`\mathrm{\Gamma }_1`$ increases.

It must be emphasized that a full model would include loss processes which, would limit the maximum $`\mathrm{\Gamma }`$ obtainable. Note that proton energies $`10^{20}`$ eV can be reached in all models computed before the saturation sets in, provided faster loss mechanisms are not present. This is despite the relatively small amount of upstream deflection allowed before particles are swept downstream on all cycles except the first which, limits the first order ’Fermi’ gain available (see Gallant and Achterberg, 1999, where the full $`\mathrm{\Gamma }_1^2`$ factor requires isotropization in each frame). Rather different spectra shapes are seen in figures 5-8 for $`\mathrm{\Gamma }`$=10, 50, 500 and 990, depending on whether it is pitch angle or large angle scattering, at $`\psi _1=15^{}`$ and $`\psi _1=35^{}`$. The sensitivity to downstream boundary position has been checked up to a distance $`r_b=210^5\lambda `$ where, $`\lambda `$ is distance between small angle scattering for this model.

The plateau-like shapes are again prominent and more ’disrupted’, compared to parallel shock cases found in the similar computations of Meli and Quenby (2002a), for large angle scatter. In the case of pitch angle diffusion the shapes are much smoother, particularly at lower $`\mathrm{\Gamma }_1`$. The disrupted shapes are probably due to high anisotropy allowed in some situations. We qualitatively understand the spectral shapes in terms of each plateau exhibiting the acceleration gain in one shock cycle, down to up to down, with a $`\mathrm{\Gamma }_1^2`$ factor boosting energy, followed by a high probability of loss downstream due to the beaming of the particles away from the shock. Examples of the beaming are presented in figure 9 where, we see the angular distributions of the transmitted particles at the shock front in the de Hoffmann-Teller frame. We only show particles which have just crossed the shock, not the complete, time averaged, distribution function. The observed effect is due to the anisotropies in pitch angle space expected for highly relativistic flows, which are in accordance with lower $`\mathrm{\Gamma }_1`$ results of Ellison et al., (1990) and Lieu et al., (1994).

These results show that, due to very relativistic velocities the distribution functions are very anisotropic and so the spatial diffusion approximation cannot apply. In figures 10-11 we present results for the super-luminal case, where $`\psi _1=48^{}`$ and $`75^{}`$ and Lorentz factors range from 10 to 40, expecting to apply in the relativistic flows of AGN environments.

The compression ratio is equal to 4 and large angle scattering is used.

To obtain the spectral shapes for the accelerated particles as the trajectories are followed during the simulations, we add the momentum of each particle in the corresponding bin every time the particle crosses the shock as before to obtain the mean energy gain and distribution function. We observe that the spectral shape falls away from a power law and descends rapidly to an upper cut off which is expected, because the particles are swept away by the flow after one cycle and have a very limited chance of return upstream, especially due to the high field inclination to the shock normal.

In the lower plots of figure 11, the energy gain of the particle is shown, versus the number of the helix-trajectory shock crossings for Lorentz factor 10 and $`\psi _1=48^{},75^{}`$ between the magnetic field and the shock normal. It is seen that after the 3<sup>rd</sup> crossing (after one cycle) the particle gains a considerable amount of energy ($`\mathrm{\Gamma }^2`$), but further gain is limited due to the fact that the particles are swept away by the flow very rapidly with no chance of ’re-cycling’ the shock again. Our calculations show that $`60\%`$ of the particles cross the shock three times (upstream$``$downstream $``$ upstream $``$downstream). It is possible to conclude that ’one-shot’ drift acceleration has taken place, that is on one up to down crossing and there is energy increase as expected from drift in the shock frame electric field under the magnetic field gradient. Some particle return upstream has been possible under the action of large angle, downstream scattering. However this return flux is so limited that a power law spectrum cannot develop and a cutoff is soon reached. Computational time limitations have currently prevented an extension of the modeling to higher $`\mathrm{\Gamma }`$ flows and pitch angle scattering models for the super-luminal case.

## 4 Discussion and Conclusions

Oblique particle shock acceleration simulations for sub-luminal and super-luminal cases have been presented. High gamma upstream flows have been used ranging from $`\mathrm{\Gamma }`$=10 up to $`10^3`$. The results found are important in answering questions concerning the efficiency and relativistic effects of first order Fermi acceleration. A major finding of both this work and the companion paper (Meli and Quenby, 2002a), is that in contrast to parallel shocks with small angle scattering, all other types of shocks and scattering models do not yield smooth, power law spectra at the shock in the $`\mathrm{\Gamma }>>1`$ regime, even though subsequent smoothing during particle escape may arise. Reasons for this will be re-iterated later. However, sub-luminal, inclined relativistic shocks provide spectra more extended in energy and therefore more approximating to a power law distribution than super-luminal shocks, the latter exhibiting a pronounced high energy cutoff feature. The result that up to 3 crossings of the shock can often occur in the adopted large angle scattering model implies that the super-luminal case is not simply ’one-shot’ shock drift. However the large gain crossing $`23`$ followed by loss to the system can well be described as the energy gain experienced on the last shock drift pass through the discontinuity standing out from the smaller, average diffusive gain.
Concerning the important ’speed up’ effect in the sub-luminal case, this appears to be due to a combination of $`\mathrm{\Gamma }^2`$ or $`\mathrm{\Gamma }`$ energy gain and the anisotropic angular distributions at the shock front (Ellison et al., 1990) which as we showed, is confined to a very small angle transmission cone. Thus, only a narrow range of pitch angle particles return to the shock and are ’counted’ for the acceleration time calculation, as the ones with larger pitch angles are lost downstream. If the particle distribution is made isotropic each side of the shock before re-crossing, the energy gain per cycle is $`\mathrm{\Gamma }_1^2`$ Quenby and Lieu (1989), but a distribution remaining anisotropic will experience less energy gain (Gallant and Achterberg, 1999). Hence some part of the ’speed up’ compared with the non-relativistic situation is expected to be connected with the un-allowed for energy enhancement. The anisotropy is of course linked to the ’beaming effect’ found also in mild relativistic cases, in accordance with Peacock (1981) and Kirk and Schneider (1988). Lieu and Quenby (1990) and Lieu et al. (1994) note that the distribution of the particle’s pitch angle is peaked towards the direction of the relativistic transformation to the de Hoffmann-Teller frame. In other words this beaming which enhances the number of particles close to the field line, may provide an additional reason for the decrease in acceleration time as the flow becomes more relativistic for both large angle and pitch angle diffusion. This ’speed-up’ resembles a trend found in the parallel shock particle acceleration case. Furthermore, the last arguments could explain the prominent spectral ’plateau’ shapes found in many oblique, sub-luminal and parallel shock cases where each shock cycle crossing seems to stand out separately. Enhanced single cycle acceleration is followed by large flux loss downstream, preventing statistical smoothing into a single power law accelerated spectrum.
The ’speed up’ effect -which has also been observed in the mildly relativistic regime Quenby and Lieu (1989) and Ellison et al. (1990), is important for a complete understanding of particle acceleration to the highest energies. An understanding of this crucial effect is necessary in any explanation for the maximum energy that particles (protons, electrons, iron nuclei) can achieve in GRB fireballs, AGN hotspots, and pulsar ultra relativistic polar winds. For example results found in the previous sections dealing with the sub-luminal case for $`\mathrm{\Gamma }10^3`$, show only few crossings are needed for the particle to achieve energies up to $`10^{19}10^{20}`$eV, although for $`\mathrm{\Gamma }`$ 10 more than ten crossings may be needed to gain such high energies. It is the competition between loss time and energy gain time which determines the maximum energy in many situations. Generally speaking the percentage increase of the maximum energy would depend on many factors. One of those is the velocity of the upstream flow. Also the dimensions of the shock could play a vital role in order to define the escape losses in a specific shock acceleration site. In addition another factor in the physical conditions which could alter the cosmic ray spectrum and decrease the spectral index is due to different energy losses occurring in the regime; for example the presence of high or low energy photons (Greisen-Zatsepin-Kuz’min cutoff effect<sup>1</sup><sup>1</sup>1This is called also the microwave background effect, where the interaction between the cosmic rays and the microwave background results in a cutoff in the energy spectrum of cosmic rays (Greisen 1966; Zatsepin-Kuz’min 1966)) or $`\gamma \gamma `$ flux interactions<sup>2</sup><sup>2</sup>2If GRB high energy photons come from electron acceleration in the $`\mathrm{\Gamma }10^3`$ regime, we might expect some ’structure’ in the spectrum at the higher energies on any scattering model. in GRB will decrease the maximum energy of the particles. In the context of AGN it is well known from a number of interesting papers such as Protheroe and Kazanas (1983) and Chakrabarti and Moltoni (1993) that a shock in the central engine of the accreting flow ($`\mathrm{\Gamma }510`$) is likely to appear and particle acceleration can consequently take place. Also there is observational evidence that a common feature of many AGN is the formation of one or two jets where shock formations should appear in conjunction with $`\mathrm{\Gamma }510`$ flows. Thus, the computed flattening of the spectra -they are all flatter than a -2 power law- as well as the acceleration ’speed up’, should have important consequences when compared with predictions and observational data, regarding for example, the photons produced by emission from electrons which can be accelerated at the shocks or produced from accelerated protons. GRB seem to be another potential candidate for the acceleration of Ultra High Energy Cosmic Rays (UHECRs), but the issue is debated by many workers. As we have mentioned Vietri (1995) and Waxman (1995), who have exploited the ’gamma-squared’ factor energy gain of the particle, were first to predict theoretically that indeed UHECRs can be produced in GRB where a pre-acceleration of the particles has been implied. Briefly, we note that observationally there is very strong evidence that pulsars (e.g. PSR 1913+16, PSR 2127+11C) are capable of injecting continuously relativistic particles into the surrounding medium over their lifetime. While this plasma will probably be predominantly in the form of $`e^+e^{}`$ pairs created in the pulsar magnetosphere, it has been argued that pulsar winds must also contain ions in order to account for the electrical current required in the Crab Nebula (Hoshino et al., 1992; Gallant and Arons, 1994). The termination shock as the $`\mathrm{\Gamma }10^5`$ flow meets the nebula is likely to be quasi-perpendicular and hence the super-luminal situation applies for the additional $`e^+e^{}`$ acceleration. The enhanced production of neutrinos in astrophysical sites such as in AGN and GRB is another consequence of the simulation findings of the ’speed up’ effect and ’gamma squared’ factor. Mastichiadis and Kirk (1992), Protheroe and Szabo (1992) and Vietri (1998) showed that the neutrino flux depends on the maximum energy attained from the primary protons, thus a decrease of the acceleration time constant could produce higher maximum energy limits for the accelerated particles and consequently an increase in the cutoff energy of the neutrino flux will be noted.

Acknowledgements

We would like to express our thanks to Prof. Drury and the unknown referee for valuable comments and suggestions.