text:
# References

Mass and Scalar Cross-sections for Neutralino Dark Matter in Anomaly Mediated Supersymmetry Breaking Model

Debasish Majumdar $`^{a,}`$ ***E-mail address: debasish@theory.saha.ernet.in aTheory Division, Saha Institute of Nuclear Physics,
1/AF Bidhannagar, Kolkata 700064, India.

We have considered neutralino to be the lightest supersymmetric particle (LSP) in the framework of minimal Anomaly Mediated Supersymmetric (mAMSB) model. We have studied variation of neutralino mass with the supersymmetric parameters. Considering these neutralinos to be the candidates for weakly interacting massive particle (WIMP) or cold dark matter (CDM), we have calculated the neutralino nucleon scalar cross-sections and compared them with DAMA-NaI neutralino direct detection search results. From this study we observe that the mAMSB model results cannot explain the allowed region in WIMP mass and WIMP-nucleon scalar cross-section space obtained from annual modulation signature in DAMA-NaI experiment.

Observational evidence like velocity curves of spiral galaxies, measurement of X-ray emissions of the cluster of galaxies, gravitational lensing and other theoretical calculations strongly suggest the existence of enormous amount of matter in the universe which is invisible to us. This invisible, nonluminous matter or ‘Dark Matter’ constitute major fraction of the mass of the universe. While there may well be more than one type of dark matter, like baryonic or non-baryonic, hot (relativistic) or cold (non-relativisitic), large scale structure calculations suggest non-relativistic massive particles (cold) that were decoupled from the plasma when the universe was entering into the epoch of matter dominance, constitute a dominant fraction of dark matter. These cold dark matters should be made up of weakly interacting massive ( $`>`$GeV) particles or WIMPs.

The favourite candidate for the WIMPs for cold dark matter, are lightest supersymmetric particles or LSP. As generally WIMPs are to be electrically neutral, the widespread choice for LSP is neutralino. Neutralino is a neutral Majorana particle and a superposition of fermionic super-partners of gauge and Higgs bosons namely bino, wino and two neutral higgsinos.

The direct detection of dark matter WIMPs can be realised by investigating their elastic scattering on the target nuclei of the detector. Hence for theoretical prediction of the WIMP detection rate it is important to calculate this cross-section for a choice of dark matter candidate.

The WIMPs in the galactic halo are assumed to have an approximately Maxwellian velocity distribution in the galactic rest frame. But for an observer on the earth (for one who is moving with the velocity of the earth) one has to include earth’s motion with respect to the halo. This has two components. One is sun’s motion in galaxy and the other is earth’s motion with respect to the sun. The directionality of this latter motion changes over the year as earth rotates around the sun. This in turn induces an annual variation of the WIMP speed relative to earth (maximum when earth’s rotational velocity adds up to the velocity of the sun and minimum when these velocities are in opposite directions). This imparts an annual modulation of the direct detection rates of the WIMPs constatnly encountered by earth . DAMA- NaI detector at Gran Sasso in Italy  looks for this annual modulation signature in WIMPs direct search. They also have published allowed 3$`\sigma `$ confidence level contours for their dark matter search in $`M_W`$ \- $`\xi \sigma _p`$ plane where $`M_W`$ is the WIMP mass, $`\sigma _p`$ is the WIMP -nucleon scattering cross section and $`\xi =\rho _\chi /\rho _{\mathrm{}}`$, $`\rho _\chi `$ being the local dark matter density and $`\rho _{\mathrm{}}`$ is the halo density (taken to be 0.3 GeV/cm3).

In this letter minimal Anomaly Mediated Supersymmetry Breaking (mAMSB) model is cosidered for generating neutralino ($`\chi `$) LSP. It is our endeavour here to investigate the mass and scalar cross-sections for these neutralinos within the range of parameters given by mAMSB model. We also investigate whether these mass and scalar cross-sections can explain DAMA-NaI dark matter search results.

In this letter we consider the supersymmetric particle neutralino as a WIMP particle which is a possible candidate for dark matter. Neutralino is the lowest mass eigenstate of linear superposition of photino ($`\stackrel{~}{\gamma }`$), zino ($`\stackrel{~}{Z}`$), and the two Higgsino states ($`\stackrel{~}{H_1^0},\stackrel{~}{H_2^0}`$)  and is written as

$$\chi =a_1\stackrel{~}{\gamma }+a_2\stackrel{~}{Z}+a_3\stackrel{~}{H_1^0}+a_4\stackrel{~}{H_2^0}.$$
(1)
The supersymmetric model we have chosen here is minimal Anomaly Mediated Supersymmetry Breaking (mAMSB) model . In this model, the observable sector (OS) and the hidden sector (HS) are in two distinct 3-brane separated by a finite bulk distance which is of the order of the compactification radius of a fifth compactified dimension. Unlike in ordinary gravity-mediated supersymmetry breaking where the breaking is transmitted from HS to OS by tree level exchanges, in anomaly mediated supersymmetry breaking model or AMSB model, the supersymmetry breaking is transmitted through loop generated superconformal anomaly. An sparticle spectrum in this model is fixed by three parameters namely $`m_{3/2}`$ which is equal to the gravitino mass, the ratio $`\mathrm{tan}\beta `$ of the vacuum expectation values of two Higgs fields $`H_1^0`$ and $`H_2^0`$ and sign($`\mu `$) ($`\mu `$ being the Higgsino mass). The problem of tachyonic sfermions in this model is remedied by cosidering a universal mass squared term $`m_0^2`$ that will make all sfermion masses positive. This is minimal anomaly mediated supersymmetry breaking or mAMSB model. Therefore, with $`m_0`$, we have four parameters in mAMSB model.

In Minimal Supersymmetric Standard Model (MSSM), the $`\chi `$ \- nucleon scattering cross-sections involve two independent parameters namely neutralino mass and squark mass. In AMSB model both the squark mass and the neutralino mass are determined by $`m_{3/2}`$ parameter. For phenomenological reasons, a universal mass term $`m_0`$ has to be added to the squark mass. However, the stability of SUSY potential  constrains the $`m_0m_{3/2}`$ parameter space severly. Similar constraints cannot be obtained in MSSM. Hence our calculation has more predictive power.

Recently several authors had used this model for various calculations . In the present calculations, with mAMSB model, we need to fix the range of parameters namely $`m_{3/2}`$, $`\mathrm{tan}\beta `$, $`m_0`$ and the sign of $`\mu `$. We have chosen Ref. for this purpose. Datta et al in Ref.  give the allowed parameter space for the mAMSB model for $`\mathrm{tan}\beta =5`$ in Fig 1. In Tables 1 and 2 of the same reference, Datta et al furnished some indicative values of the parameters (including lower bounds on $`m_0`$ and upper bounds on $`m_{3/2}`$) from their calculations. Using them we have chosen the following ranges for $`\mathrm{tan}\beta `$ and $`m_0`$.

$`5`$ $`\mathrm{tan}\beta `$ $`15`$
$`405\mathrm{GeV}`$ $`m_0`$ $`905\mathrm{GeV}`$ (2)
The range of the parameter $`m_{3/2}`$ (in TeV) for each choice of $`m_0`$ is fixed by the allowed region in Fig 1 of Ref. . The sign of $`\mu `$ is taken to be negative.

We attempt here to calculate the scalar cross-section for different WIMP mass. The WIMP here is neutralino dark matter. The neutralino mass and neutralino nucleon scalar scattering cross-sections are calculated in the framework of mAMSB model.

Bottino et al, who have done many pioneering work in dark matter calculations gave an expression for neutralino-nucleon scalar scattering cross-section . This expression (Eq. 2 of  and Eq. 1 of ) reads as

$$\sigma _{\mathrm{scalar}}=\frac{8G_F^2}{\pi }M_Z^2m_{red}^2\left[\frac{F_hI_h}{m_h^2}+\frac{F_HI_H}{m_H^2}+\mathrm{squark}\mathrm{exchange}\mathrm{term}\right]^2$$
(3)
The two terms in bracket in the above equations are $`h`$ and $`H`$ exchange terms where $`h`$ and $`H`$ are CP even neutral Higgs bosons.

$`F_h`$ $`=`$ $`(a_1\mathrm{sin}\theta _W+a_2\mathrm{cos}\theta _W)(a_3\mathrm{sin}\alpha +a_4\mathrm{cos}\alpha )`$
$`F_H`$ $`=`$ $`(a_1\mathrm{sin}\theta _W+a_2\mathrm{cos}\theta _W)(a_3\mathrm{cos}\alpha +a_4\mathrm{sin}\alpha )`$
$`I_{h,H}`$ $`=`$ $`{\displaystyle \underset{q}{}}k_q^{h,H}m_qN|\overline{q}q|N`$ (4)
In the above $`a_1`$, $`a_2`$, $`a_3`$, $`a_4`$ are the coefficients as shown in Eq.(1), the angle $`\alpha `$ is the Higgs mixing angle that rotates $`H_1^0`$ and $`H_2^0`$ into $`h`$ and $`H`$, $`m_h`$ and $`m_H`$ are the masses of $`h`$ and $`H`$ respectively, $`m_{\mathrm{red}}^2`$ is the neutralino-nucleon reduced mass. The last term in Eq. (3) can be written as

$$I_{h,H}=k_{u\mathrm{type}}^{h,H}g_u+k_{d\mathrm{type}}^{h,H}g_d$$
(5)
where $`k_{u\mathrm{type}}^h=\mathrm{cos}\alpha /\mathrm{sin}\beta `$, $`k_{d\mathrm{type}}^h=\mathrm{sin}\alpha /\mathrm{cos}\beta `$, $`k_{u\mathrm{type}}^H=\mathrm{sin}\alpha /\mathrm{sin}\beta `$, $`k_{d\mathrm{type}}^H=\mathrm{cos}\alpha /\mathrm{cos}\beta `$ . The values of $`g_u`$ and $`g_d`$, in terms of pion-nucleon sigma term, $`\sigma _{\pi N}`$, are given in Ref. . In the present calculation, for simplicity, we deal only with Higgs mediated terms. Also it is pointed out in Ref.  that in general these terms are largely dominant over the squark exchange term.

For the calculation of scalar cross-sections (Eq. 3-5) one has to compute $`a_1`$, $`a_2`$, $`a_3`$, $`a_4`$, the Higgs mixing angle $`\alpha `$, the Higgs masses $`m_h`$ and $`m_H`$, the neutralino mass $`m_\chi `$. These are calculated in the framework of mAMSB model for the range of supersymmetric parameters discussed above, using the code ISAJET 7.48 . We have checked that the neutralino is the LSP. In Fig. 1 we show the variation of neutralino mass $`m_\chi `$ with $`\mathrm{tan}\beta `$ for three sets of values of $`m_0`$ and $`m_{3/2}`$ with $`\mu <0`$. Fig. 1 shows that there are no significant variations of $`m_\chi `$ with $`\mathrm{tan}\beta `$. In Fig. 2 we plot $`m_0`$ (in GeV) vs $`m_\chi `$ for three fixed sets of values of $`m_{3/2}`$ and $`\mathrm{tan}\beta `$. These three sets are chosen by fixing $`m_{3/2}`$ at a lower value 31.5 TeV with three values of $`\mathrm{tan}\beta `$ namely 5.0, 10.0 and 15.0. It is observed from Fig. 2 that although the plot corresponding to $`\mathrm{tan}\beta =5`$ is rather flat the plots for $`\mathrm{tan}\beta =10`$ and 15.0 show some variations of neutralino mass with $`m_0`$. More so the mean difference of $`m_\chi `$’s corresponding to $`\mathrm{tan}\beta =5`$ and $`\mathrm{tan}\beta =10`$ is almost double than that of $`\mathrm{tan}\beta =10`$ and $`\mathrm{tan}\beta =15`$. The variation of $`m_\chi `$ is therefore a collective effect of variation of the three mAMSB model parameters.

For calculation of neutralino-nucleon scalar cross-sections using Eq. (3-5), we express all the quantities in the units of GeV and the cross-sections are expressed in pico barn (pb). The values of $`g_u`$ and $`g_d`$ (in Eq. 5) are obtained from  and they are $`g_u=123`$ GeV, $`g_d=288`$ GeV. The results of calculated scalar cross-sections, $`\sigma _{\mathrm{scalar}}`$ (in pb) for various netralino masses ($`m_\chi `$ in GeV) are shown as scatter plot in Fig. 3. The particular pattern of this scatter plot for cross-section calculations using mAMSB model is reflective of the shape of the allowed region for $`m_0m_{3/2}`$ parameter space for mAMSB model as shown in Fig 1. of Datta et al .

In order to study whether these results can represent the DAMA-NaI results stated above, we also show in Fig. 3, the results obtained from DAMA-NaI experiment. They are shown as contour plots in Fig. 3. The contours in Fig. 3 (continuous plots) are obtained from Fig. 4b of Ref.  where they have given 3$`\sigma `$ C.L. allowed regions for experiments in four different yearly cycles (DAMA-NaI-1 to DAMA-NaI-4, with total statistics of 57986 kg.day) alongwith DAMA-NaI-0 results . The plots suggest that although the neutralino mass range falls within that obtained from DAMA-NaI experiments but the cross-sections are about an order of magnitude below than the DAMA contour. An attempt to find out $`\xi =\rho _\chi /\rho _{\mathrm{}}`$ where $`\rho _{\mathrm{}}=0.3`$ GeV/cm3  using the relation $`\xi =(\xi \sigma )_{DAMA}/\sigma _{calc.}`$, leads to values of $`\rho _\chi 3`$ GeV/cm3 which is much more than the halo density $`\rho _{\mathrm{}}`$ and thus is unphysical.

To conclude, we have considered minimal Anomaly Mediated Supersymmetry Breaking model for studying the variation of neutralino mass with the supersymmetric parameters, namely $`\mathrm{tan}\beta `$, $`m_0`$ and $`m_{3/2}`$. We have calculated the netralino nucleon scalar cross-sections and studied their variations with neutralino mass. Considering neutralino to be cold dark matter candidate, we went on to compare the results with those obtained from DAMA-NaI detection of annual modulation of dark matter signal. It appears that the mAMSB model calculations cannot reproduce DAMA-NaI experimental findings for WIMPs. Although earlier results from calculations within MSSM frame-work can represent the DAMA results to some extent, for mAMSB model it is not so and it is due to the fact that the parameter space of mAMSB model is somewhat different from MSSM. As discussed earlier, mAMSB model has only four parameters and $`m_0m_{3/2}`$ parameter space is severly constrained. This affects the calculation of the coefficients $`a_1`$, $`a_2`$ etc. and the value of $`\alpha `$. The scalar cross-sections thus obtained using Eqs. (3 - 5) for different $`m_\chi `$’s are therefore different from those obtained in MSSM. We repeat our calculations within the framework of mAMSB model with different $`g_u`$, $`g_d`$ values given in Table 2 of a more recent work by Bottino et al . We find these cross-sections fall by orders of magnitudes from those calculated using . These results are further off from the results of DAMA-NaI experiments. We have plotted these results in Fig. 4 for demonstration. The calculation of neutralino relic density in this model requires the calculations of neutralino-neutralino annihilation cross-sections and these will constitute our future work.

The author is grateful to Anindya Dutta for introducing him to the computer code ISAJET and for many valuable discussions. The author also thanks Amitava Datta for some useful discussions.

Figure Captions

Fig. 1 Variation of neutralino mass $`m_\chi `$ (in GeV) with $`\mathrm{tan}\beta `$ in mAMSB model for three sets of values of $`m_0`$ (in GeV) and $`m_{3/2}`$ in TeV. $`\mu <0`$.

Fig. 2 Variation neutralino mass $`m_\chi `$ (in GeV) with $`m_0`$ in GeV for three sets of values of $`m_{3/2}`$ and $`\mathrm{tan}\beta `$. See text for detail.

Fig. 3 Variation of neutralino nucleon scalar cross-section $`\sigma _{\mathrm{scalar}}`$ in pb (pico barn) with neutralino mass $`m_\chi `$ in (GeV) (scatter plot) and their comparison with DAMA-NaI results (see text for details). For DAMA-NaI results (closed contours) $`\xi \sigma `$ in pb is plotted along the Y-axis while for the results of present calculations $`\sigma _{\mathrm{scalar}}`$ in pb are plotted along Y-axis.

Fig. 4 Comparison of cross-section ($`\sigma _{\mathrm{scalar}}`$ in pb) calculations for different sets of values for $`g_u`$ and $`g_d`$. The symbols ’X’ and circle corresponds to data sets (a) and (b) respectively in Table 2 of , the symbol ’.’ (dots) are for the same sets as in Fig. 3.


format: text

text:
# Measurements of 𝛼_𝑠 from event shapes and the four–jet rate Talk presented at the 31st International Conference on High Energy Physics, Amsterdam, 24-31 July 2002.

## 1 Introduction

The strong coupling constant $`\alpha _s`$ is the only free parameter of the QCD Lagrangian and many measurement have been performed at various energy scales in different processes in order to determine it precisely . The measurements at LEP have contributed significantly to this effort.

Two years after the closure of LEP still new results from QCD analyses become available . These results are based on data from LEP1, where the centre-of-mass energy was around 91 GeV and the data sample consisted of several million events, practically background free, as well as from LEP2, with energies from 133 GeV up to 206 GeV. Here the data samples are much smaller, of the order of several thousand events, and the backgrounds (mainly fully hadronic W-pair events) can amount up to about 15%. In addition, at LEP2 a special selection has to be applied in order to eliminate the so-called radiative return events, where a hard photon is radiated in the initial state, reducing the effective centre-of-mass energy to the Z mass. This leads to a significant loss in statistics. Therefore the $`\alpha _s`$ measurements at LEP2 have a much larger statistical uncertainty. However, systematic uncertainties related to non-perturbative corrections or unknown higher order terms are reduced because of the larger energy scale. This motivates the combination of all measurements performed by the four LEP experiments ALEPH, DELPHI, OPAL and L3, in order to obtain a precise test of the running of the coupling, as well as a competitive measurement at LEP2 compared to LEP1, and finally the best possible overall combination. Such a combination has now been performed by the LEP QCD working group (LEPQCDWG) for measurements based on event–shape variables.

## 2 Event–shape variables

Event–shape variables are a classical tool to measure the strong coupling constant. They are sensitive to the topology of an event, therefore to gluon radiation and thus to the coupling strength. The most thoroughly tested variables are thrust, heavy jet mass, C–parameter, total and wide jet broadening and the $`y_3`$–distribution, where $`y_3`$ is the resolution parameter where a three–jet event is clustered into a two–jet event (Durham clustering algorithm). For these variables the perturbative QCD predictions exist at next–to–leading order (NLO), $`𝒪(\alpha _s^2)`$. In addition, large logarithms $`L=\mathrm{ln}(1/y)`$ have been resummed at next–to–leading logarithmic (NLL) accuracy to all orders in $`\alpha _s`$, where $`y`$ is an event–shape variable. Only recently some missing NLL terms have been calculated for the $`y_3`$–distribution . Several matching schemes have been proposed over the years in order to combine the NLO and NLL calculations, such as the logR and R matching schemes . The LEPQCDWG adopts so–called modified matching schemes, where the logarithmic term is replaced by $`L\stackrel{~}{L}=\frac{1}{p}\mathrm{ln}\left[(1/y)^p(1/y_{\mathrm{max}})^p+1\right].`$ This ensures a physical behaviour of the prediction at the phase space boundary $`y_{\mathrm{max}}`$. The choice of the power $`p`$ is arbitrary, in the sense that different choices only give different subleading contributions, whereas the NLL terms remain unchanged. The modified logR matching scheme corresponds to $`p=1`$, whereas $`p=2`$ is called the second degree modification.

As mentioned earlier, the measurements at LEP1 are basically background free, whereas at LEP2 the W background is very difficult to suppress, in particular in the multi–jet region of the event–shape distribution which is also sensitive to hard gluon radiation. In order to avoid large systematic uncertainties related to the background subtraction, a careful choice of the fit range in the event–shape distribution is needed, over which the theoretical prediction is adjusted to the data. At the other end of the spectrum (the two–jet limit) the hadronization corrections become sizeable, and correspondingly the fit range has to be adapted. These corrections are obtained from phenomenological models implemented in Monte Carlo generators (JETSET, HERWIG, ARIADNE). The difference in the results obtained with different generators are quoted as systematic uncertainties due to hadronization, of the order of 1-2%.

The dominant systematic uncertainty is related to missing higher order terms in the perturbative calculations. In order to avoid a very bad fit quality when changing the renormalization scale, and in particular to reduce the multitude and arbitrariness of prescriptions used by the LEP experiments, the LEPQCDWG has proposed a new method, the so–called uncertainty band method. First a reference prediction is chosen, such as the modified logR matching scheme, together with a certain $`\alpha _s`$ value. Then various changes are applied to the prediction, keeping the same $`\alpha _s`$ value: The renormalization scale $`\mu `$ is varied in the range $`0.5<\mu /\sqrt{s}<2`$, the modified R matching scheme is applied, the second degree modification is tested, the logarithmic terms are changed to $`L^{}=\mathrm{ln}(1/(x_Ly))`$ for $`2/3x_L3/2`$, which introduces different subleading terms, and the kinematic boundary $`y_{\mathrm{max}}`$ is varied over a range given by the difference between NLO and parton shower predictions. The difference of the distributions thus obtained to the reference distribution defines an uncertainty band (Fig. 1). In the next step again the reference prediction is taken, but now $`\alpha _s(M_\mathrm{Z})`$ is varied such that the changes in the prediction reproduce the spread given by the uncertainty band over the fit range. The necessary variations in $`\alpha _s(M_\mathrm{Z})`$ define the theoretical uncertainty.

The LEPQCDWG has combined all measurements based on event–shape distributions at LEP1 and LEP2 for the six variables mentioned above. Correlations between variables, energies and experiments have been taken into account. The results are $`\alpha _s(M_\mathrm{Z})=0.1197\pm 0.0002_{\mathrm{stat}}\pm 0.0008_{\mathrm{exp}}\pm 0.0010_{\mathrm{had}}\pm _{0.0047}^{0.0048}{}_{\mathrm{theo}}{}^{}`$ for LEP1 alone, $`\alpha _s(M_\mathrm{Z})=0.1196\pm 0.0005_{\mathrm{stat}}\pm 0.0010_{\mathrm{exp}}\pm 0.0007_{\mathrm{had}}\pm _{0.0044}^{0.0043}{}_{\mathrm{theo}}{}^{}`$ for LEP2 alone, and $`\alpha _s(M_\mathrm{Z})=0.1198\pm 0.0003_{\mathrm{stat}}\pm 0.0009_{\mathrm{exp}}\pm 0.0008_{\mathrm{had}}\pm 0.0046_{\mathrm{theo}}`$ for all LEP energies combined (including measurements at an effective centre-of-mass energy below $`M_\mathrm{Z}`$, obtained by selecting events with hard initial state radiation). Excellent consistency between the LEP1 and LEP2 measurements is found. Taking the individual measurements at the various energy points, perfect agreement with the expected running of the strong coupling is observed (Fig. 2).

In an analysis by DELPHI the energy dependence of the mean values of event–shape variables has been analysed , but now the hadronization corrections are replaced by a power law as proposed in . It results in a rather precise determination of $`\alpha _s(M_\mathrm{Z})=0.1184\pm 0.0033_{\mathrm{tot}}`$. A more detailed discussion of power law studies can be found in .

## 3 Four–jet rate

Recently, the ALEPH collaboration has published a new determination of $`\alpha _s`$ , based on a measurement of the four-jet rate (Durham algorithm) at LEP1. The difference to the event–shape variables is that here the perturbative expansion starts at $`𝒪(\alpha _s^2)`$. The full NLO corrections are known, as well as the resummation of NLL terms . The quoted result is extremely precise, namely $`\alpha _s(M_\mathrm{Z})=0.1170\pm 0.0001_{\mathrm{stat}}\pm 0.0009_{\mathrm{exp}}\pm 0.0003_{\mathrm{had}}\pm 0.0008_{\mathrm{theo}}`$. This is mainly due to the larger sensitivity, since the leading term is of $`𝒪(\alpha _s^2)`$. It is observed that although the $`\chi ^2`$ of the fit changes substantially when varying the renormalization scale over a wide range, and also the location of the sharp minimum in this $`\chi ^2`$ varies after changes in the fitted prediction, the actually fitted $`\alpha _s(M_\mathrm{Z})`$ value remains remarkably stable.

## 4 Summary

New combined results for measurements of $`\alpha _s(M_\mathrm{Z})`$ at LEP have been presented, based on event–shape distributions and a new prescription for the theoretical uncertainty. Furthermore, it has turned out that the measurement of the four-jet rate at LEP allows for a very precise determination of $`\alpha _s(M_\mathrm{Z})`$.

## 5 Acknowledgements

I would like to thank the members of the LEPQCDWG for providing me with the relevant numbers and plots.


text:
# The Υ Spectrum from Lattice QCD with 2+1 Flavors of Dynamical Quarks

## 1 INTRODUCTION

Measurement of the $`\mathrm{{\rm Y}}`$ spectrum gives an excellent test for studying the dynamical content of the configurations being used. NRQCD can be used to give precise results which indicate how the light dynamical quarks feed in to the heavy meson properties. The program at CLEO-c will enable lattice results on $`\mathrm{{\rm Y}}`$ spectroscopy to be experimentally verified .

Here we present new results using dynamical configurations from the MILC collaboration . We have used the coarser set of available ensembles. These were generated using the improved staggered Asqtad action and have the advantage over previous configurations of the incorporation of 2+1 flavors of dynamical quarks. These are therefore the most realistic simulations of the physical vacuum to date. Several 2+1 flavor ensembles are available with one dynamical quark being held at approximately the value of the strange quark mass $`m_s`$ and two degenerate flavors varying from $`m_s`$ to $`m_s/5`$. Also available are a quenched ensemble and a two flavor ensemble with degenerate dynamical quark masses of $`m_s/2.5`$. We have used both the most chiral 2+1 flavor ensemble and the degenerate 3 flavor ensemble plus the quenched and 2 flavor ensembles to best demonstrate the effects of unquenching on the bottomonium spectrum.

We use the standard $`𝒪(v^4)`$ NRQCD Lagrangian

$``$ $`=`$ $`\psi ^{}(D_t+H_0)\psi +\psi ^{}\delta H\psi ;`$ (1)
$`H_0`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }^{(2)}}{2M}}`$ (2)
$`\delta H`$ $`=`$ $`{\displaystyle \frac{(𝚫^{(2)})^2}{8M^3}}+{\displaystyle \frac{ig}{8M^2}}(𝚫^{(\pm )}𝐄𝐄𝚫^{(\pm )})`$ (3)
$``$ $`{\displaystyle \frac{g}{8M^2}}\sigma (𝚫^{(\pm )}\times 𝐄𝐄\times 𝚫^{(\pm )})`$
$``$ $`{\displaystyle \frac{g}{2M}}\sigma 𝐁+{\displaystyle \frac{a^2\mathrm{\Delta }^{(4)}}{24M}}{\displaystyle \frac{a(\mathrm{\Delta }^{(2)})^2}{16nM^2}}`$
but with corrections to the $`E`$ and $`B`$ fields for $`𝒪(a^2)`$ errors. All terms are tadpole improved using an estimate of the Landau gauge link.

## 2 FITTING

For all ensembles, 210 configurations were used. For all but the degenerate 3 flavor ensemble, 4 different quark (antiquark) origins were used on each of 4 different timeslices. For the degenerate 3 flavor ensemble, 4 different origins were used on each of 2 different timeslices. In fitting the meson correlators, the Bayes technique  was used. $`3\times 3`$ matrix fits were done using local, ground and excited state smearings for the S states and $`2\times 2`$ matrix fits were done using ground and excited state smearings for P states. The Bayes parameters used were quite relaxed with $`ln(E_0)=1.1(7)`$, each of the amplitudes = 0.1(7), and splittings between excited states $`500\text{MeV}\pm \text{ factor 2}`$.

## 3 RESULTS FOR THE LATTICE SPACING

In Figure 1 results are plotted for the ratio of $`a^1`$ obtained from both the 1P-1S and 2S-1S splittings in the $`\mathrm{{\rm Y}}`$ system to that obtained from using $`r_0`$, assuming that $`r_0=0.5\text{fm}`$. The results from the dynamical configurations are plotted against the ratio of the pseudoscalar to vector light-light meson mass where the valence quark mass is equal to the $`u/d`$ sea quark mass. We will also be able to extract a reliable scale from the quantity $`2m_{B_s}m_\mathrm{{\rm Y}}`$ once $`m_s`$ is well determined . Dependence on the quantity used to obtain $`a^1`$ with quenching gives ambiguities in calculating other physical observables. The plot clearly illustrates these ambiguities disappearing as $`m_\pi /m_\rho `$ tends to its physical value. For the most chiral 2+1 flavor dynamical ensemble the scales from $`\mathrm{{\rm Y}}`$ have converged. The $`\mathrm{{\rm Y}}`$ results imply that $`r_0<0.5\text{fm}`$. These scales can be used to determine a value for $`\alpha _s`$ .

## 4 RESULTS FOR THE LEPTONIC WIDTH AND HYPERFINE SPLITTING

Figure 2 shows the ratio of the leptonic width, i.e. the width of decay to $`e^+e^{}`$, of the 1st radially excited $`\mathrm{{\rm Y}}`$ to that of the ground state $`\mathrm{{\rm Y}}`$. The experimental points for $`\psi `$ and $`\mathrm{{\rm Y}}`$ are plotted to demonstrate that this ratio is to a large extent independent of the heavy quark mass. At leading order $`\mathrm{\Gamma }_{ee}(nS)|\psi _n(0)|^2`$, where the ‘wavefunction at the origin’ $`\psi _n(0)=n|l`$, the overlap of the local vector current with a given state. Therefore we can calculate this ratio using appropriate amplitudes from our $`3\times 3`$ fits. The $`Z`$ factor for this current, which must be calculated in perturbation theory, cancels in the ratio. It can be seen that although even the most chiral result is still some way from experiment, the ratio is clearly falling with dynamical quark mass. Next, the currents must be improved by adding current correction operators as has been done for $`f_B`$ .

The hyperfine splitting between the $`\mathrm{{\rm Y}}`$ and $`\eta _b`$ is a spin dependent quantity which is proportional to the square of the coefficient of the $`\sigma .B`$ term in the action and also, in leading order in a potential model, to $`|\psi _n(0)|^2`$. In Figure 3 we show the ratio of hyperfine splitting of the 1st excited to that of the ground $`\mathrm{{\rm Y}}`$ state. The uncertainty from the coefficient of the $`\sigma .B`$ term, and, similarly to the leptonic width, a lot of the heavy quark mass dependence cancels in the ratio. The $`\eta _b`$ and $`\eta _b^{}`$ are experimentally unseen as yet so the experimental result given is for the $`\psi `$ system from the BELLE collaboration . The plot shows a tendency towards experiment with unquenching.

## 5 CONCLUSION

The dependence of the scale on the quantity used to obtain it has been a long standing problem of the quenched approximation. Full convergence of $`a^1`$ obtained from different splittings in the $`\mathrm{{\rm Y}}`$ spectrum is now clearly seen on the most chiral 2+1 flavor MILC ensemble. We conclude that this ensemble has good dynamical content. Ratios of leptonic widths and hyperfine splitting have been measured and are seen to tend toward experiment wih unquenching and decreasing dynamical quark mass.

## 6 ACKNOWLEDGEMENTS

We are grateful to MILC for use of their configurations. This work was supported by NSF, DoE, PPARC, the EU, and NSERC Canada.


text:
# Factors Determining Variability Time in Active Galactic Nucleus Jets

## 1 Introduction

In the standard picture of active galactic nuclei (AGN), accretion onto a super-massive black hole is via an accretion disk, and a significant fraction of the accretion power (possibly supplemented by tapping into the rotational energy of the black hole) produces twin opposing relativistic jets moving outward along the disk axis, with typical Lorentz factors $`\mathrm{\Gamma }2`$–10 as inferred from very long baseline interferometry (VLBI) observations. The objects observed in high energy $`\gamma `$ rays are “blazars”, AGN in which one of the jets is closely aligned toward the observer. It is natural that in $`\gamma `$-rays we should see preferentially those AGN with aligned jets because the emission from the jet is Doppler boosted in energy and relativistically beamed along the jet direction (for a discussion of relativistic effects see Urry and Padovani 1995). The $`\gamma `$ ray emission from blazars is variable (as it is also at optical, UV and X-ray energies). Relativistic effects also cause the observed variability time to be shorter than the time scale over which the emission changes in the jet frame.

The spectral energy distribution (SED) of blazars shows two broad peaks, the low energy peak extending from the infrared to the UV or X-ray region of the spectrum, and the high energy peak starting in the X-ray or $`\gamma `$ ray range. The usual interpretation is that relativistic electrons produce the low energy part by synchrotron emission, and that the same electrons produce the high energy part by Compton scattering the low energy part and/or external photons to higher energies. The 3rd EGRET catalog of high-energy $`\gamma `$-ray sources (Hartman et al. 1999) contains around 70 high confidence identifications of AGN, and all appear to be blazars (Montigny et al. 1995, Mukherjee et al. 1997). Clearly, the $`\gamma `$-ray emission is associated with AGN jets.

Four BL Lac objects have been detected in the TeV energy range: Mrk 421 (Punch et al. 1992), Mrk 501 (Quinn et al. 1996), 1E S2344+514 (Catanese et al. 1998) and PKS 2155-304 (Chadwick et al. 1999). Recently, the spectrum of Mrk 501 has been measured up to 24 TeV by the HEGRA telescopes (Konopelko et al. 1999). Several of the EGRET AGN show $`\gamma `$-ray variability with time scales of $`1`$ day (Kniffen et al. 1993) at GeV energies. The TeV $`\gamma `$-ray emission of two BL Lacs shows very rapid variability. For Mrk 421, variability on a time scale as short as $`15`$ minutes has been reported (Gaidos et al. 1996). In the case of Mrk 501, variability on a time scale of a few hours was observed during the 1997 high level of activity, and there is evidence of a 23 day periodicity (Protheroe et al. 1998, Hayashida et al. 1998) interpreted in terms of a binary black hole model for the central engine by Rieger and Mannheim (2000). These variability timescales place important constraints on the models. For example, the synchrotron self-Compton (SSC) model appears to be just consistent with recent multi-wavelength observations of Mrk 421 and Mrk 501 during flaring activity (Bednarek and Protheroe 1997, 1999). However, the allowed range of physical parameters (Doppler factor and magnetic field) is rather small, and this mechanism may well be excluded by future observations. For a recent review of TeV $`\gamma `$-ray astronomy see Kifune (2002).

Rapid variability in intra-day variability (IDV) sources is a long-standing problem as it implies apparent brightness temperatures in the radio regime which may exceed 1017 K or relativistic beaming with extremely high Doppler factors, coherent radiation mechanisms, or special geometric effects (Wagner and Witzel 1995). The very rapid flaring observed at TeV and X–ray energies during flaring activity in blazars also presents a challenge for any model and suggests a re-examination of mechanisms which may cause very rapid variability would be worthwhile. In this paper I concentrate on how the observed variability time is related to the geometry and motion of the emission region, and thus to the photon energy density in the emission region. The blazar emission mechanisms to be discussed include: a shock excited emission region, bent jets, a shock propagating along a jet containing a helical structure and illuminating parts of the helix by enhanced interactions/emission of radiation such that the emission regions move along helical paths, and highly oblique conical shocks in the jet. Together with geometry-specific time delays and variable Doppler boosting associated with relativistic motion of the emission region along a curved trajectory, it may well be possible to explain the observed flaring activity and high brightness temperatures. Another possibility briefly discussed in the context of bent jets and conical shocks is that jets may be fueled on an irregular time scale.

The observed variability time $`\mathrm{\Delta }t_{\mathrm{obs}}`$, and some assumed or estimated Doppler factor $`D`$ is often used to estimate the jet-frame source radius, $`R^{}Dc\mathrm{\Delta }t_{\mathrm{obs}}`$. The jet-frame photon energy density is then usually assumed to be $`U_{\mathrm{phot}}^{}L^{}/4\pi R_{}^{}{}_{}{}^{2}c`$, where $`L^{}`$ is the jet-frame bolometric luminosity given by $`L^{}=D^44\pi d_L^2F`$, $`F`$ being the observed bolometric flux, and $`d_L`$ being the luminosity distance. However, this approach can lead to large systematic errors in the jet-frame photon energy density. This is important because the energy density of photons of the low energy part of the SED may determine the energy losses of electrons, and the rate of up-scattering to $`\gamma `$-ray energies in SSC models, and the rate of proton-photon collisions in hadronic models (see, e.g., Mannheim and Biermann 1989, Protheroe 1997, Mannheim et al. 2001, Mücke and Protheroe 2001, Mücke et al. 2002). If the emission region is optically thin, it may also have consequences for IDV sources as I will show that it is quite possible that the photon energy density responsible for the so-called Compton catastrophe may actually be lower than usually estimated. In the following sections I shall discuss how the variability time is related to the emission region geometry and intrinsic jet-frame variability, show how this can lead to large systematic errors in the jet-frame photon energy density, and discuss other scenarios which can lead to rapid variability and high fluxes.

## 2 Relationship between variability time and emission region geometry

I shall discuss first a relativistic jet pointing at angle $`\theta `$ with respect to the line of sight to the observer, in which the emission region is a cylinder of radius $`r`$ and jet-frame length $`\mathrm{}^{}`$ moving along the jet with the jet’s Lorentz factor $`\mathrm{\Gamma }=1/\sqrt{1\beta ^2}`$. To work out how the observed lightcurve depends on the emission region geometry and the duration of the emission in the jet frame (primed coordinates) we first consider in the observer frame (unprimed coordinates) the events corresponding to: (i) the emission from centre of cylinder at $`(tt_\mathrm{i}=0,xx_\mathrm{i}=0,yy_\mathrm{i}=0,zz_\mathrm{i}=0)`$ with this event defining the origin of coordinates in both frames, i.e. $`(t^{}t_\mathrm{i}^{}=0,x^{}x_\mathrm{i}^{}=0,y^{}y_\mathrm{i}^{}=0,z^{}z_\mathrm{i}^{}=0)`$ ; (ii) the emission from some arbitrary point in the cylinder at a later time $`(t_{\mathrm{ii}},x_{\mathrm{ii}},y_{\mathrm{ii}},z_{\mathrm{ii}})`$ ; (iii) the arrival at the telescope, located in the $`x`$$`y`$ plane, of the photon emitted in event (i); (iv) the arrival at the telescope of the photon emitted in event (ii) (see Fig. 1).

In this paper, I shall assume that the distance between the AGN and the observer is very much larger than the dimensions of the emission region ($`R`$ and $`\mathrm{}`$ in Fig. 1), and that we are not concerned with investigating variability on time scales very much shorter than $`R/c`$ or $`\mathrm{}/c`$. Under this assumption, the time interval between arrival of the two photons at the telescope is independent of distance to the AGN and is simply given by

$`t_{\mathrm{obs}}(t_{\mathrm{iv}}t_{\mathrm{iii}})=\left(t_{\mathrm{ii}}{\displaystyle \frac{x_{\mathrm{ii}}}{c}}\mathrm{cos}\theta {\displaystyle \frac{y_{\mathrm{ii}}}{c}}\mathrm{sin}\theta \right).`$ (1)
The error made by using this approximation is $`(R/d_L)(R/c)`$ and is negligible compared with the variability time scale associated with the emission region geometry investigated in this paper which is $`(r/c)`$.

Lorentz transformation to the observer frame, gives $`t_{\mathrm{ii}}=\mathrm{\Gamma }(t_{\mathrm{ii}}^{}+\beta x_{\mathrm{ii}}^{})`$, $`x_{\mathrm{ii}}=\mathrm{\Gamma }(x_{\mathrm{ii}}^{}+\beta t_{\mathrm{ii}}^{})`$ and $`y_{\mathrm{ii}}=y_{\mathrm{ii}}^{}`$. Noting that the Doppler factor is $`D=[\mathrm{\Gamma }(1\beta \mathrm{cos}\theta )]^1=\mathrm{\Gamma }(1+\beta \mathrm{cos}\theta ^{})`$, and using the Aberration formulae,

$`\mathrm{cos}\theta `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }}{D}}(\mathrm{cos}\theta ^{}+\beta ),\mathrm{sin}\theta ={\displaystyle \frac{1}{D}}\mathrm{sin}\theta ^{}`$ (2)
$`\mathrm{cos}\theta ^{}`$ $`=`$ $`\mathrm{\Gamma }D(\mathrm{cos}\theta \beta ),\mathrm{sin}\theta ^{}=D\mathrm{sin}\theta ,`$ (3)
Eq. 1 becomes

$`t_{\mathrm{obs}}`$ $`=`$ $`\mathrm{\Gamma }(t_{\mathrm{ii}}^{}+\beta x_{\mathrm{ii}}^{})\mathrm{\Gamma }(x_{\mathrm{ii}}^{}+\beta t_{\mathrm{ii}}^{}){\displaystyle \frac{\mathrm{\Gamma }}{Dc}}(\mathrm{cos}\theta ^{}+\beta )y_{\mathrm{ii}}^{}{\displaystyle \frac{1}{Dc}}\mathrm{sin}\theta ^{}`$ (4)
$`Dt_{\mathrm{obs}}`$ $`=`$ $`t_{\mathrm{ii}}^{}{\displaystyle \frac{x_{\mathrm{ii}}^{}}{c}}\mathrm{cos}\theta ^{}{\displaystyle \frac{y_{\mathrm{ii}}^{}}{c}}\mathrm{sin}\theta ^{}.`$ (5)
Let us suppose that the cylindrical emission region emits radiation simultaneously and uniformly throughout its volume between times $`t^{}=\mathrm{\Delta }t^{}`$ and $`t^{}=\mathrm{\Delta }t^{}`$, as measured in the jet frame. While a simultaneous burst violates causality, it is nevertheless a useful case to consider because it enables us to determine clearly the contribution of emission region geometry to the observed variability time. For $`\theta ^{}\pi /2`$ the first photon to arrive would have been emitted at $`(t^{}=\mathrm{\Delta }t^{},x^{}=\mathrm{}^{}/2,y^{}=r,z^{}=0)`$ giving

$`t_{\mathrm{obs}}^{\mathrm{first}}=\mathrm{\Delta }t^{}(\mathrm{}^{}/2c)\mathrm{cos}\theta ^{}(r/c)\mathrm{sin}\theta ^{}.`$ (6)
Similarly, the last photon to arrive would have been emitted at $`(t^{}=\mathrm{\Delta }t^{},x^{}=\mathrm{}^{}/2,y^{}=r,z^{}=0)`$ giving

$`t_{\mathrm{obs}}^{\mathrm{last}}=\mathrm{\Delta }t^{}+(\mathrm{}^{}/2c)\mathrm{cos}\theta ^{}+(r/c)\mathrm{sin}\theta ^{}.`$ (7)
For $`\theta ^{}>\pi /2`$ the first photon to arrive would have been emitted at $`(t^{}=\mathrm{\Delta }t^{},x^{}=\mathrm{}^{}/2,y^{}=r,z^{}=0)`$, and the last photon to arrive would have been emitted at $`(t^{}=\mathrm{\Delta }t^{},x^{}=\mathrm{}^{}/2,y^{}=r,z^{}=0)`$. Hence, if we define the observer frame duration of the burst as $`2\mathrm{\Delta }t_{\mathrm{obs}}`$ then

$`D\mathrm{\Delta }t_{\mathrm{obs}}=\mathrm{\Delta }t^{}+(\mathrm{}^{}/2c)|\mathrm{cos}\theta ^{}|+(r/c)\mathrm{sin}\theta ^{}`$ (8)
and this is valid for all $`\theta ^{}`$. Note that Eq. 8 gives the usual formula, $`\mathrm{\Delta }t_{\mathrm{obs}}=\mathrm{\Delta }t^{}/D`$, if the emission region is point-like (i.e. $`\mathrm{}^{}=r=0`$). If one term in Eq. 8 dominates, $`D\mathrm{\Delta }t_{\mathrm{obs}}`$ gives one of: $`\mathrm{\Delta }t^{}`$, $`(\mathrm{}^{}/2c)|\mathrm{cos}\theta ^{}|`$, or $`(r/c)\mathrm{sin}\theta ^{}`$.

I shall next consider the case of a cylindrical emission region in the jet being rapidly energized by a plane shock with jet-frame speed $`\beta _{\mathrm{shock}}^{}c`$ travelling along the jet, such that photons are emitted immediately after shock passage from a thin disk-like region immediately downstream of the shock. In this case, the location of the emitting disk is defined by $`x^{}=\beta _{\mathrm{shock}}^{}ct^{}`$, and so the arrival times of photons at the telescope may be obtained from Eq. 5

$`cDt_{\mathrm{obs}}`$ $`=`$ $`x^{}\left({\displaystyle \frac{1}{\beta _{\mathrm{shock}}^{}}}\mathrm{cos}\theta ^{}\right)y^{}\mathrm{sin}\theta ^{}.`$ (9)
For $`(1/\beta _{\mathrm{shock}}^{}\mathrm{cos}\theta ^{})<1`$ the first and last photons to be received would have been emitted at $`(x^{}=\mathrm{}^{}/2,y^{}=r,z^{}=0)`$ and $`(x^{}=\mathrm{}^{}/2,y^{}=r,z^{}=0)`$, respectively. However, for $`(1/\beta _{\mathrm{shock}}^{}\mathrm{cos}\theta ^{})>1`$ the first and last photons to be received would have been emitted at $`(x^{}=\mathrm{}^{}/2,y^{}=r,z^{}=0)`$ and $`(x^{}=\mathrm{}^{}/2,y^{}=r,z^{}=0)`$, respectively. Hence, the observer frame duration of the burst is

$`D\mathrm{\Delta }t_{\mathrm{obs}}=(\mathrm{}^{}/2c)\left|{\displaystyle \frac{1}{\beta _{\mathrm{shock}}^{}}}\mathrm{cos}\theta ^{}\right|+(r/c)\mathrm{sin}\theta ^{}.`$ (10)
If $`\theta ^{}`$ is small ($`\theta `$ is very small), one finds

$`\mathrm{\Delta }t_{\mathrm{obs}}`$ $``$ $`{\displaystyle \frac{\mathrm{}^{}}{2Dc}}\left|{\displaystyle \frac{1\beta _{\mathrm{shock}}^{}}{\beta _{\mathrm{shock}}^{}}}\right|+{\displaystyle \frac{r\theta ^{}}{Dc}}`$ (11)
and for a reasonable shock speed, e.g. $`\beta _{\mathrm{shock}}^{}\pm 0.5`$,

$`\mathrm{\Delta }t_{\mathrm{obs}}`$ $``$ $`(21){\displaystyle \frac{\mathrm{}^{}}{2Dc}}+{\displaystyle \frac{r\theta ^{}}{Dc}}`$ (12)
where $`\beta _{\mathrm{shock}}^{}0.5`$ corresponds to a reverse shock.

What we have learned from this discussion is that multiplying $`c\mathrm{\Delta }t_{\mathrm{obs}}`$ by the Doppler factor might give the jet-frame intrinsic variability time or one of the dimensions of the emission region (possibly multiplied by some unknown factor). One is tempted to ask if $`R^{}Dc\mathrm{\Delta }t_{\mathrm{obs}}`$ is the right dimension to put in $`U_{\mathrm{phot}}^{}=L^{}/4\pi R_{}^{}{}_{}{}^{2}c`$ in order to estimate the jet-frame photon energy density. We shall discuss this point further in Section 4.

## 3 Monte Carlo investigation of $`c\mathrm{\Delta }t_{\mathrm{obs}}`$, size, and $`D`$

The Monte Carlo method allows the accurate calculation of expected light curves for any emission region geometry and intrinsic source variability. The emission region is modeled in the jet frame, and is represented by $`N`$ “particles”, each of which emits precisely one “photon”. The number density of the particles models the geometry of the emission region, and each particle emits its photon at a time determined by the emission region geometry and variability model. The jet-frame 4-position of each photon emission event $`(t_i^{},x_i^{},y_i^{},z_i^{})`$, $`i=1,2,\mathrm{}N`$, is determined by model. The emission region moves in the $`x`$-direction with Lorentz factor $`\mathrm{\Gamma }`$, and the observer-frame 4-position of each photon emission event $`(t_i,x_i,y_i,z_i)`$, $`i=1,2,\mathrm{}N`$, is obtained by Lorentz transformation. The arrival time of each photon is calculated for a given viewing angle $`\theta `$ using Eq. 1, and is binned to give the lightcurve.

Two emission region geometries are considered: a solid cylinder, and a spheroidal 3D gaussian. We also consider four intrinsic jet-frame time distributions: a simultaneous burst (violates causality), a simultaneous gaussian pulse (violates causality) and simultaneous emission at plane shock. I shall investigate effects of varying the Doppler factor $`D`$, viewing angle $`\theta `$, jet Lorentz factor $`\mathrm{\Gamma }`$, and shock speed $`\beta _{\mathrm{shock}}`$. Finally, the effect of smoothing due to acceleration/radiation time delays is discussed.

### 3.1 Solid cylinder emission region geometry

If the emission region is “solid”, i.e. the emissivity is constant inside the emission region and zero outside, the lightcurve may be peaked and have finite duration to reflect the sharp edges of the emission region, and the shape of the lightcurve will change with with viewing angle. The Monte Carlo results for the cylinder, when plotted such that the time of observation is divided by the expected variability time $`\mathrm{\Delta }t_{\mathrm{obs}}`$ given by Eqs. 8 and 10, is shown in Fig. 2(a) for a simultaneous burst ($`\mathrm{\Delta }t^{}=0`$) , and in Fig. 2(b) for shock excitation. As can be seen, Eqs. 8 and 10 are verified by the Monte Carlo results. In both cases the shape of the lightcurve depends strongly on the viewing angle.

### 3.2 Gaussian emission region geometries

The density of “particles” representing the emission region is described by

$$\rho (x^{},y^{},z^{})\mathrm{exp}\left[\left(\frac{x_{}^{}{}_{}{}^{2}}{2\sigma _{x}^{}{}_{}{}^{2}}+\frac{y_{}^{}{}_{}{}^{2}+z_{}^{}{}_{}{}^{2}}{2\sigma _{r}^{}{}_{}{}^{2}}\right)\right]$$
(13)
with $`\sigma _x=\sigma _r`$ corresponding to a spherical gaussian distribution, $`\sigma _x<\sigma _r`$ to an oblate spheroidal gaussian distribution, and $`\sigma _x>\sigma _r`$ to an prolate spheroidal gaussian distribution. For the case of a simultaneous gaussian pulse, the probability of emission at jet-frame time $`t^{}`$ to $`(t^{}+dt^{})`$ is $`p(t^{})dt^{}`$, where $`p(t^{})=\mathrm{exp}(t_{}^{}{}_{}{}^{2}/2\sigma _{t}^{}{}_{}{}^{2})/\sqrt{2\pi }\sigma _t`$ is independent of position.

If the emission region has a gaussian shape the lightcurve will be smooth, and in many cases will also have a gaussian shape. This is true, for example, for the cases of a simultaneous gaussian pulse and excitation by a plane shock. The width of the lightcurve will depend on viewing angle in approximately the same way as for equivalent cylindrical emission volume. For example, $`r`$ and $`\mathrm{}^{}/2`$ of the solid cylinder should be related to $`\sigma _r`$ and $`\sigma _x`$, respectively, for the case of the spheroidal gaussian density. For the case of a simultaneous gaussian burst I find the standard deviation to be given by

$$\sigma =[(\sigma _r\mathrm{sin}\theta ^{}/c)^2+(\sigma _x\mathrm{cos}\theta ^{}/c)^2+\sigma _t^2]^{1/2}/D.$$
(14)
Note the similarity to Eq. 8 except that the terms are added in quadrature as the standard deviation is required instead of the maximum duration of the pulse which is, theoretically, infinite.

In the case of excitation by a plane shock wave with speed $`\beta _{\mathrm{shock}}c`$, the emission is from the plane where the shock cuts the spherical gaussian density. The emission region then has a surface density of emitting particles which is a two-dimensional gaussian, and the lightcurve reflects this distribution, and so is also gaussian. The standard deviation of the gaussian lightcurve depends on the viewing angle as a result of projection effects, and I find the standard deviation to be

$$\sigma =\frac{1}{cD}\left[(\sigma _r\mathrm{sin}\theta ^{})^2+(\sigma _x)^2\left(\frac{1}{\beta _{\mathrm{shock}}^{}}\mathrm{cos}\theta ^{}\right)^2\right]^{1/2}$$
(15)
Lightcurves for $`\mathrm{\Gamma }=10`$, $`\beta _{\mathrm{shock}}^{}=0.5`$ and various viewing angles are plotted in Fig. 3 and seen to lie on top of each other when plotted in units of $`\sigma `$, and to be a normal distribution. Equation 15 is valid for both forward and reverse shocks. For reverse shocks, $`\beta _{\mathrm{shock}}^{}`$ would be negative and give rise to a broader lightcurves (larger $`\sigma `$).

### 3.3 Acceleration/radiation delays

Instantaneous excitation by a shock wave is not a realistic approximation unless the time-scales for particle acceleration and radiation (cooling) are very short compared to the transit time of the shock through the emission region. If this is not the case, the lightcurve would be broadened and smoothed to reflect the time delays associated with particle acceleration and radiation. An equation of the form $`dN^{}/dt^{}=Q^{}(t^{})N^{}/t_{\mathrm{loss}}^{}`$ might describe the time-evolution of the number of particles, $`N^{}`$, radiating photons at the frequency corresponding to that observed, with $`Q^{}(t^{})`$ being the source term, and $`t_{\mathrm{loss}}^{}`$ representing the time-scale for particle losses (or time-scale over which the radiation is emitted). Then a simple smoothing function of the form

$$\varphi _{\mathrm{del}}(t_{\mathrm{del}}^{})=\{\begin{array}{cc}0\hfill & (t_{\mathrm{del}}^{}<0)\hfill \\ \frac{1}{\mathrm{\Delta }t_{\mathrm{gain}}^{}}\left[1\mathrm{exp}\left(\frac{t_{\mathrm{del}}^{}}{t_{\mathrm{loss}}^{}}\right)\right]\hfill & (0t_{\mathrm{del}}^{}<\mathrm{\Delta }t_{\mathrm{gain}}^{})\hfill \\ \frac{1}{\mathrm{\Delta }t_{\mathrm{gain}}^{}}\left[1\mathrm{exp}\left(\frac{\mathrm{\Delta }t_{\mathrm{gain}}^{}}{t_{\mathrm{loss}}^{}}\right)\right]\mathrm{exp}\left[\frac{(t_{\mathrm{del}}^{}\mathrm{\Delta }t_{\mathrm{gain}}^{})}{t_{\mathrm{loss}}^{}}\right]\hfill & (\mathrm{\Delta }t_{\mathrm{gain}}^{}t_{\mathrm{del}}^{})\hfill \end{array}$$
(16)
can be used for this purpose, where $`t_{\mathrm{del}}^{}`$ represents the delay between shock passage and emission by the radiating particles, $`\mathrm{\Delta }t_{\mathrm{gain}}^{}`$ represents the the duration of the acceleration following shock passage. In the case of a pre-existing population of thermal electrons shock heating is essentially instantaneous ($`\mathrm{\Delta }t_{\mathrm{gain}}^{}=0`$), and $`t_{\mathrm{loss}}^{}`$ represents the time-scale for cooling by thermal bremsstrahlung radiation. Similarly, for a pre-existing population of relativistic electrons, the passage of the shock will essentially instantly increase the magnetic field ($`\mathrm{\Delta }t_{\mathrm{gain}}^{}=0`$) and $`t_{\mathrm{loss}}^{}`$ would represent the time-scale for energy losses by synchrotron emission. If particle acceleration is required then $`\mathrm{\Delta }t_{\mathrm{gain}}^{}>0`$.

Fig. 4 shows the effect of acceleration/radiation delays on the lightcurve due to a plane shock ($`\beta _{\mathrm{shock}}^{}=0.3`$) exciting an oblate gaussian spheroidal emission region ($`\sigma _x/\sigma _r=0.2`$) with acceleration/radiation delays described by $`\mathrm{\Delta }t_{\mathrm{gain}}^{}=t_{\mathrm{loss}}^{}`$ for $`t_{\mathrm{loss}}^{}=0.1\sigma _r/c`$, $`\sigma _r/c`$ and $`10\sigma _r/c`$, and viewed at angle $`\theta ^{}=5^{}`$. In the case of no acceleration/radiation delays, the lightcurves would simply be normal distributions centred on $`t_{\mathrm{obs}}=0`$. Results are shown for the three $`t_{\mathrm{loss}}^{}`$-values, and $`t_{\mathrm{obs}}`$ is plotted in terms of $`\sigma _{\mathrm{tot}}`$ defined by

$$\sigma _{\mathrm{tot}}^2=\sigma ^2+\frac{t_{\mathrm{loss}}^{}{}_{}{}^{2}+\mathrm{\Delta }t_{\mathrm{gain}}^{}{}_{}{}^{2}}{4cD},$$
(17)
where $`\sigma `$ is given by Eqn. 15, such that $`\sigma _{\mathrm{tot}}`$ gives a crude measure of the expected duration of the lightcurve. These distributions can be obtained simply by convolution of a normal distribution with $`\varphi _{\mathrm{del}}`$ (Eq. 16), taking account of the fact that jet-frame times enter in Eq. 16, that for any point co-moving with the jet $`\mathrm{\Delta }t_{\mathrm{obs}}=\mathrm{\Delta }t^{}/D`$, and that $`t_{\mathrm{obs}}`$ is plotted in units of $`\sigma _{\mathrm{tot}}`$,

$$p\left(\frac{t_{\mathrm{obs}}}{\sigma _{\mathrm{tot}}}\right)=\frac{\sigma D}{\sqrt{2\pi }}_{\mathrm{}}^{\mathrm{}}\varphi _{\mathrm{del}}(t_{\mathrm{obs}}Dt)\mathrm{exp}\left(\frac{t^2\sigma _{\mathrm{tot}}^2}{2\sigma ^2}\right)𝑑t.$$
(18)
The solid curves, which agree with the Monte Carlo results (histograms), are obtained from Eqn. 18. As we see, if the time-scales for the acceleration/radiation process ($`\mathrm{\Delta }t_{\mathrm{gain}}^{}`$ and $`t_{\mathrm{loss}}^{}`$) are much less than the time-scale $`\sigma ^{}`$ associated with the shock passage and dimensions of the emission region then the light-curve will be symmetrical, and in the case of a gausian spheroidal emission region will be a gaussian distribution with standard deviation $`\sigma `$ (leftmost histogram). If the time-scales for the acceleration/radiation process are much larger than $`\sigma ^{}`$, then the light-curve will reflect that of the acceleration/radiation process, i.e. Eqn. 16 (rightmost histogram). The middle histogram shows an intermediate case.

Lightcurves of AGN show many different features and forms of variability. One example which appears to show variations reflecting the acceleration/radiation process is 3C 454.3. In Fig. 5, I show the lightcurve of 3C 454.3 at 37 GHz obtained with the Metsahovi and Crimea telescopes over ten years (Salonen et al. 1987, Terasranta et al. 1992). The flares appear non-symmetrical and have a shape similar to the rightmost histogram in Fig. 4 ($`\mathrm{\Delta }t_{\mathrm{gain}}^{}=t_{\mathrm{loss}}^{}\sigma `$). As an example, I have constructed a reasonably well-fitting lightcurve (solid curve) from a number of flares of the form given by Eqn. 16 with $`\mathrm{\Delta }t_{\mathrm{gain}}=t_{\mathrm{loss}}=0.75`$ yr (dotted curves) plus a background flux density of 3 Jy. One interpretation of these data would then be that the radio emission region in 3C 454.3, being modelled by the solid curve, has dimensions much less than $`c`$(0.75 yr)$`D=0.16D`$ pc, and that the energy-loss time-scale of the radiating electrons is $``$(0.75 yr)$`D`$.

Another example is the X-ray lightcurve of 1ES1959+65 from ARGOS/USA and RXTE/ASM (Giebels et al. 2002) shown in Fig. 6. In this case, the individual flares appear to be symmetrical in time, and have a roughly gaussian shape. I have constructed a reasonably well-fitting lightcurve (solid curve) from a number of flares of gaussian form with standard deviations $`\sigma =2.5`$ d (dotted curves) plus a background flux density of 2 mCrab. One could interpret this as indicating that the energy-loss time-scale of the radiating electrons is much less than $`2.5D`$ d, and that the emission region dimensions are $`c2.5D`$ d$`=2\times 10^3D`$ pc.

## 4 Dependence of energy density on dimensions of the blob

It is important to know the dimensions of the emission region for several reasons: (i) in some hadronic models the synchrotron photons are targets for photoproduction, (ii) in all models the synchrotron photons are targets for photon-photon pair production by $`\gamma `$-rays, (iii) in SSC models the synchrotron photon energy density determines Compton scattering, and (iv) knowing the systematic errors on photon energy density may help understand the so-called “Compton catastrophe” in IDV sources which have apparent brightness temperatures well in excess of the limit $`T_B<10^{12}`$ K imposed by Compton scattering (Kellermann and Pauliny-Toth 1969, Kardashev 2000) when the photon energy density in the emission region reaches the energy density in the magnetic field. One extreme example is PKS 0405-385 (Kedziora-Chudczer et al. 1997) which has $`T_B>5\times 10^{14}`$ K after correcting for interstellar scintillation (see also Walker 1998) requiring a Doppler factor of $`D=10^3`$ to satisfy the brightness temperature limit.

To illustrate how critically the energy density depends on the geometry, I shall consider the case of the jet-frame emissivity following a spheroidal gaussian density. Provided the emission is optically thin, as is almost certainly true for the optical–X-ray synchrotron hump in the SED of blazars, then it is straightforward to calculate the average energy density from the emission region geometry and the luminosity.

Assuming that the jet-frame luminosity, $`L^{}`$, is constant, we can estimate the average jet-frame photon energy density given the Doppler factor and emission region geometry. The simplest way of doing this, for any emission region geometry, is to use the Monte Carlo method to place $`N`$ points at positions $`\stackrel{}{r}_i,i=1,\mathrm{}N`$, distributed according to the emission region geometry, and to give each point a luminosity $`L^{}/N`$. Then at point $`\stackrel{}{r}_i`$ the energy density is

$`U_{\mathrm{phot}}^{}(\stackrel{}{r}_i){\displaystyle \frac{L^{}}{4\pi cN}}{\displaystyle \underset{ji}{}}r_{ij}^2`$ (19)
where $`r_{ij}=|\stackrel{}{r}_i\stackrel{}{r}_j|`$. Averaging over the emission region distribution we obtain

$`U_{\mathrm{phot}}^{}{\displaystyle \frac{L^{}}{4\pi cN^2}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{ji}{}}r_{ij}^2.`$ (20)
The result giving the jet-frame average photon energy density versus $`\sigma _x`$ is shown as the thick solid curve in Fig. 7 for the case of a gaussian spheroidal density.

In Fig. 7 we also plot the jet-frame average photon energy density that would be inferred if we assumed that the observed variability time scale $`\sigma `$ and an assumed or estimated Doppler factor $`D`$ gave the jet-frame radius of a spherical emission region, i.e. $`R^{}=\sigma Dc`$. This “observed” jet-frame average photon energy density is simply given by

$`U_{\mathrm{phot},\mathrm{obs}}^{}{\displaystyle \frac{L^{}}{4\pi (\sigma Dc)^2c}}.`$ (21)
Note that $`\sigma `$ (given by Eq. 14) depends on $`\theta ^{}`$, $`\sigma _t`$ and $`\sigma _x`$ so that the “observed” jet-frame average photon energy density depends also on $`\theta ^{}`$, $`\sigma _t`$ and $`\sigma _x`$ as well as $`D`$. The figure shows $`U_{\mathrm{phot},\mathrm{obs}}^{}`$ plotted against $`\sigma _x`$ and gives the range due to variation in $`\theta ^{}`$ (shaded) for various $`\sigma _t`$. The upper bound in each case gives the result for $`\theta ^{}=0`$ or $`\pi `$ and the lower bound is for $`\theta ^{}=\pi /2`$. We see that the “observed” value, $`U_{\mathrm{phot},\mathrm{obs}}^{}`$ can be several orders of magnitude higher or lower than $`U_{\mathrm{phot}}^{}`$ if the emission region is different from a sphere, or if the intrinsic variability time $`\sigma _t`$ is not small. For example, take the case of $`\sigma _t\sigma _r/c`$, if $`\sigma _x=10^2\sigma _r`$ and $`\theta ^{}=0`$ ($`\theta =0`$) then $`U_{\mathrm{phot},\mathrm{obs}}^{}10^3U_{\mathrm{phot}}^{}`$, whereas if $`\sigma _x=10^2\sigma _r`$ and $`\theta ^{}=90^{}`$ ($`\theta =5.74^{}`$ for $`\mathrm{\Gamma }=10`$) then $`U_{\mathrm{phot},\mathrm{obs}}^{}10^2U_{\mathrm{phot}}^{}`$. The above result has clear implications for both leptonic and hadronic models of AGN in which photons of the low-energy peak of the SED provide target photons for inverse-Compton scattering by electrons (leptonic models) or pion photoproduction by protons (hadronic models). Using the observed variability time together with assumed or estimated Doppler factor to estimate the emission region radius $`R^{}`$ can clearly lead to large errors in $`U_{\mathrm{phot}}^{}`$.

Although the radio emission in IDV sources is usually assumed to be optically thick, if this is not the case then the above result may also have implications for IDV sources as the photon energy density responsible for causing the brightness temperature limit may actually be a few orders of magnitude lower than estimated on the basis of time-variability, and in that case much lower Doppler factors would be required to avoid the Compton catastrophe. I explore this further in a separate paper (Protheroe 2002).

## 5 Variability due to photon pile-up in observation time

The simplest example of photon pile-up in observation time is a bent jet. Jets may bend if they pass through a stratified cold and high density region (Mendoza and Longair 2001). We approximate the trajectory of an emission region moving with speed $`\beta c=(11/\mathrm{\Gamma }^2)^{1/2}`$ along a bent section of jet by motion around a section of a circle of radius $`r`$ in the $`x`$$`y`$ plane: $`(t,x=r\mathrm{cos}(\omega t),y=r\mathrm{sin}(\omega t),z=0)`$, where $`\omega =\beta c/r`$. For an observer in the $`x`$$`z`$ plane at angle $`\eta `$ to the $`x`$-axis, the observation time is then given by

$`t_{\mathrm{obs}}=t(r/c)\mathrm{cos}(\omega t)\mathrm{cos}\eta `$ (22)
and the Doppler factor is

$`D=\{\mathrm{\Gamma }[1+\beta \mathrm{cos}\eta \mathrm{sin}(\omega t)]\}^1.`$ (23)
The Doppler factor raised to the 4th power is plotted against observation time for $`\mathrm{\Gamma }=10`$ and various observation angles in Fig. 8. I have plotted $`D^4`$ as it is appropriate for bolometric flux from a moving isotropic source; it also applies to the specific flux $`F_\nu `$ for $`F_\nu \nu ^1`$. We see that even for modest observation angles the lightcurve is strongly peaked, essentially a delta function when the emission region direction is closest to the line of sight. Of course the finite size of any emission region will broaden the distribution. For example, for a lab-frame emission region length along the jet $`\mathrm{}`$, the burst would have duration $`\mathrm{\Delta }t_{\mathrm{obs}}\mathrm{}/c`$.

### 5.1 Helical jet structures

VLBI observations show that helical jets or helical structures in jets may be fairly common in AGN (Rantakyro et al. 1998), and theoretical studies have shown that wave-like helical structures can occur as a result of jet precession (Hardee 2000). Several papers discuss helical jet models or the application of helical models to specific sources (Camenzind 1986, Rosen 1990, Tateyama et al. 1998, Qian et al. 1992, Schramm et al. 1993, Steffen et al. 1995, Villata and Raiteri 1999). Certainly helical jets or structures would be important in determining the lightcurve of $`\gamma `$-ray and neutrino emission from blazars, and various suggestions have been made about the mechanisms involved (Despringre and Fraix-Burnet 1997, Marcowith et al. 1995).

I consider a filamentary helical structure embedded in the jet with Lorentz factor $`\mathrm{\Gamma }`$ whose axis coincides with the jet axis and is excited by a plane shock travelling along the jet with jet-frame speed $`\beta _{\mathrm{shock}}^{}c`$. The helical structure could be, for example, a flux tube containing a relatively high magnetic field, or a tube of high plasma density arising from a density perturbation in the plasma entering the jet. Helical magnetic fields with an Archemedian spiral topology similar to the “Parker spiral” field of the heliosphere may well be expected in AGN jets.

The excitation of a filamentary helical structure is described by the jet-frame 4-vector $`(t^{},x^{},y^{},z^{})=(\lambda ^{}\varphi /\beta _{\mathrm{shock}}^{}c,\lambda ^{}\varphi ,r\mathrm{cos}\varphi ,r\mathrm{sin}\varphi )`$, where $`r`$ is the radius of the cylinder containing the helix, $`\lambda `$ is the helix wavelength, $`\lambda ^{}=\lambda /2\pi `$, $`\varphi =\beta _{\mathrm{shock}}^{}ct^{}/\lambda ^{}`$, and the jet is pointing in the $`x`$-direction. Lorentz transformation to the galaxy-frame gives $`(t,x,y,z)=[\mathrm{\Gamma }t^{}(1+\beta \beta _{\mathrm{shock}}^{}),\mathrm{\Gamma }ct^{}(\beta +\beta _{\mathrm{shock}}^{}),r\mathrm{cos}(\beta _{\mathrm{shock}}^{}ct^{}/\lambda ^{}),r\mathrm{sin}(\beta _{\mathrm{shock}}^{}ct^{}/\lambda ^{})]`$. The galaxy-frame speed of the location of the excited part of the helix can exceed $`c`$, but this does not violate causality as no particles or information propagates at this pattern speed which is

$`v_{\mathrm{pattern}}={\displaystyle \frac{\beta _{\mathrm{shock}}^{}+\beta }{1+\beta _{\mathrm{shock}}^{}\beta }}\left[1+\left({\displaystyle \frac{r\beta _{\mathrm{shock}}^{}}{\lambda ^{}\mathrm{\Gamma }(\beta _{\mathrm{shock}}^{}+\beta )}}\right)\right]c.`$ (24)
Observation in the $`x`$-$`z`$ plane at angle $`\theta `$ to the jet axis ($`x`$-axis) gives

$`t_{\mathrm{obs}}`$ $`=`$ $`\mathrm{\Gamma }t^{}[1+\beta \beta _{\mathrm{shock}}^{}(\beta _{\mathrm{shock}}^{}+\beta )\mathrm{cos}\theta ]+r\mathrm{sin}\theta \mathrm{sin}(\beta _{\mathrm{shock}}^{}ct^{}/\lambda ^{})/c`$ (25)
$`{\displaystyle \frac{dt_{\mathrm{obs}}}{dt^{}}}`$ $`=`$ $`\mathrm{\Gamma }[1+\beta \beta _{\mathrm{shock}}^{}\mathrm{cos}\theta (\beta _{\mathrm{shock}}^{}+\beta )]+r\mathrm{sin}\theta (\beta _{\mathrm{shock}}^{}/\lambda ^{})\mathrm{cos}(\beta _{\mathrm{shock}}^{}ct^{}/\lambda ^{}).`$ (26)
Whenever $`dt_{\mathrm{obs}}/dt^{}=0`$ the lightcurve will have a cusp. If cusps are possible, they will occur at times corresponding to

$`\mathrm{cos}\varphi `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }[1+\beta \beta _{\mathrm{shock}}^{}\mathrm{cos}\theta (\beta _{\mathrm{shock}}^{}+\beta )]}{r\mathrm{sin}\theta (\beta _{\mathrm{shock}}^{}/\lambda ^{})}}`$ (27)
provided the model parameters give $`\mathrm{cos}\varphi `$ in the range $`1\mathrm{cos}\varphi 1`$ ($`\mathrm{cos}\varphi `$ depends on the helix geometry, shock speed, jet Lorentz factor and viewing angle). If $`\mathrm{cos}\varphi =\pm 1`$ one cusp per helix wavelength will occur, and if $`1<\mathrm{cos}\varphi <1`$ multiple cusps occur, otherwise no cusps are present in the light curve. However, if $`\mathrm{cos}\varphi `$ is close to $`\pm 1`$ the lightcurve will be peaked. This is illustrated in Fig. 9 for $`\mathrm{\Gamma }=10`$, $`\beta _{\mathrm{shock}}^{}=0.5`$, and $`\lambda =r`$. Fig. 9(a) shows $`t^{}`$ plotted against $`t_{\mathrm{obs}}`$ for five viewing angles, and Fig. 9(b)–(f) shows the resulting lightcurve for each of the five viewing angles (the corresponding Doppler factors are also given).

Apart from the periodicity, the lightcurves shown in Fig. 9 are reminiscent of those of blazars. All of these lightcurves correspond to instantaneous emission from the point on the helical filament at the time of excitation. The cusps and peaks would in reality be smoothed to some extent by the finite width of any helical structure as well as by any delays associated with acceleration and radiation time scales. Furthermore, as the shock weakens while propagating down the jet the successive peaks/cusps in the lightcurve would decrease in height. One possible scenario could be that there is a succession of shocks at random intervals propagate down the jet, and because individual shocks weaken as they propagate, each shock would cause only one or two peaks (due to only one or two cycles of the helix). In this case, it may not be possible to distinguish between a simple bent structure in the jet and a helical structure – if appropriately aligned, a simple bent structure would cause a cusp in the lightcurve in exactly the same way as a helical structure. An alternative to a helical filamentary structure within the jet is a helical jet which would produce a qualitatively-similar lightcurve, but would have the Doppler factor varying with position along the jet as the viewing angle relative to the local jet direction changes.

### 5.2 Conical shocks

Lind and Bladford (1985) have considered the possibility that hotspots seen in VLBI images of radio jets may actually be relativistically moving conical shocks. They applied the relativistic shock jump conditions to a plane parallel flow entering a forward conical shock with cone angle $`\eta `$ to find the angle $`\zeta <\eta `$ with at which the flow initially diverges with respect to the jet axis. Defining the downstream region between the cones with angles $`\zeta `$ and $`\eta `$ to be the emission region, and taking account of Doppler boosting, they model the brightness distribution to simulate VLBI images.

The generation of conical shock structures is often seen to occur in simulations of relativistic jets after the introduction of a fast perturbation (Bowman 1994, Gomez et al. 1997). These structures, which are typically an alternating sequence of forward and reverse conical shocks, are stationary or slowly-moving in the galaxy frame, and have cone angles $`\eta 1/\mathrm{\Gamma }`$. Salvati et al. (1998) discuss emission by a conical shock. They consider the case of a density perturbation, confined to thin flat disk, travelling relativistically along the jet causing particle acceleration and emission where the disk cuts a stationary forward conical shock such as that illustrated in Fig. 10.

Salvati et al. (1998) assume that the time-scales for acceleration and emission are negligible. They show that their model can lead to highly peaked light curves whose shape depends on the viewing angle, and used these light curves to fit Markarian 421 TeV flare data. For the same input, and assuming there is no Doppler boosting of the emitted radiation, I am able to reproduce exactly their figure 2 which shows the observed flux for various viewing angles. However, bearing in mind that the shock is stationary and that the pile-up in observing time is already included, the bolometric flux emitted by part of the downstream flow will be Doppler boosted by $`D_{\mathrm{local}}^3`$ where $`D_{\mathrm{local}}=[\mathrm{\Gamma }_d(1\beta _d\mathrm{cos}\theta _{\mathrm{local}})]^1`$ where $`\mathrm{\Gamma }_d`$ is the Lorentz factor of the downstream flow, and $`\theta _{\mathrm{local}}`$ is the viewing angle with respect to the line of sight and the local downstream flow direction which varies around the shock as indicated in Fig. 10 ($`\theta _{\mathrm{local}}`$ ranges between $`\theta _a`$ and $`\theta _b`$). To obtain the Lorentz factor of the downstream flow, it is easiest to Lorentz transform in a direction parallel to the shock plane to a frame in which the flow is normal to the shock. For the case of cone angle $`\eta =\mathrm{sin}^1(1/\mathrm{\Gamma })`$ assumed by Salvati et al. (1998), i.e. an oblique shock at angle $`\eta `$ to the upstream flow, and using the relativistic equation of state, I find that for $`\mathrm{\Gamma }=10`$, $`\eta 7.25^{}`$ and $`\zeta 1.89^{}`$ and that the Lorentz factor of the downstream flow is related to that of the upstream flow by $`\mathrm{\Gamma }_d0.801784\mathrm{\Gamma }_u`$. My result for the observed flux for the same input as Salvati et al., but including the Doppler boosting taking into account the local downstream flow directions, is given in Fig. 11(a) and shows that the inclusion of Doppler boosting causes the peak at $`\theta =0`$ to be higher than that at $`\theta =0.9\eta `$, the opposite to that found by Salvati et al. However, since Salvati et al. (1998) used one viewing angle in their fits to Markarian 421 TeV flare data (their figure 3), and because the divergence of the downstream flow is rather small, their fits are still valid and their model remains an interesting mechanism for flare production. The same authors (Spada et al. 1999) have applied their model to IDV sources and are able to explain brightness temperatures up to $`3\times 10^{17}`$ K.

I wish to extend the work of Salvati et al. (1998) by including the reverse shocks, and ultimately a sequence of stationary reverse and forward shocks. In Fig. 11(b) I show the lightcurve for a reverse conical shock having identical parameters as the forward shock already discussed. Because the divergence of the downstream flow is rather small, the lightcurve of the reverse shock is approximately just the lightcurve of the forward shock reflected about $`t_{\mathrm{obs}}=0`$.

I shall consider next the case of a sequence of stationary reverse and forward conical shocks. Although in real AGN jets it may be possible that the downstream flow is re-accelerated to near the original value (depending on conditions external to the jet, and the jet production mechanism), I shall assume that each successive conical shock causes a reduction in the jet Lorentz factor by $`\mathrm{\Gamma }_d/\mathrm{\Gamma }_u=0.801784`$. I shall make an additional approximation that the divergence/convergence of the flow caused by the conical shocks can be neglected, and that the jet flow is always parallel to the jet axis. The lightcurve due to a thin, initially flat, disk travelling relativistically along the jet is then calculated by the Monte Carlo method, in which “particles” are placed over the surfaces of all of the cones. Each cone has the same number of particles, and they would appear uniformly distributed across the cross section of the jet when viewed along the jet axis. A side view showing the location of these particles is given in Fig. 12. An initially flat thin disk is launched along the jet with initial Lorentz factor $`\mathrm{\Gamma }=10`$, and each time a part of the disk crosses one of the shocks the Lorentz factor of that part of the disk drops appropriately, distorting the disk. The time at which each “particle” emits its photon is determined by the time at which the (distorted) disk reaches the “particle”. The resulting lightcurve for various viewing angles is shown in Fig. 13(a) where the emission is boosted using the Doppler factor corresponding to the Lorentz factor of the flow immediately downstream of each shock. If the Lorentz factor of the jet decreases at each shock crossing as assumed here then, because the Doppler factor is also reduced at each shock crossing, only the first few conical shocks would be prominent in the lightcurve resulting from a single jet perturbation. This is shown in Fig 13(b) which gives the lightcurve from Fig 13(a) for $`\theta =9^{}`$ on a linear scale.

As can be seen, the lightcurve is quasi-periodic. It would be strictly periodic if the downstream flow velocity were identical to the upstream flow velocity. This might be the case if the re-collimmation of the diverging flow from the conical shocks results in restoration of the flow Lorentz factor to near the upstream value, and in this case the flux from successive cycles would be at about the same level, as observed in Markarian 501 flares. The two flow velocities would also be roughly the same if the initial jet Lorentz factor were much higher than $`\mathrm{\Gamma }=10`$ used in Fig. 13(a) such that the intervals between flares due to a pair of conical shocks was roughly the same. However, in this case the flux of successive flares would diminish as the Doppler factor decreases. Turning this argument around, we may have a method of determining the minimum jet Lorentz factor from observational data, e.g. the $`23`$ day periodicity in the 1997 Markarian 501 data (Protheroe et al. 1998, Hayashida et al. 1998) may be used to put a lower limit on $`\mathrm{\Gamma }`$.

Fig. 14 shows an example lightcurve due to multiple jet perturbations. In this case, I have used an exponential distribution of times between the injection into the jet of a density perturbation with a mean interval of $`r/c`$. I have also used an exponential distribution for the strength of each perturbation. The resulting lightcurve looks, at least qualitatively, as good as any model for the variability of fluxes from blazars. Note that even though there is no strict periodicity resulting from the excitation of a sequence of shocks in this case, each flare episode has two strong peaks due to (for this viewing angle) the apex of the first reverse shock cone and the apex of the first forward shock cone, and in each flare the separation between these two peaks is identical. Fourier analysis would therefore show a strong peak at the frequency corresponding to the time interval between these two peaks in a single flare.

## 6 Conclusion

Many factors can influence observed variability time. The connection between $`\mathrm{\Delta }t_{\mathrm{obs}}`$, Doppler factor and emission region geometry is non-trivial, and so measuring $`\mathrm{\Delta }t_{\mathrm{obs}}`$ may give, at best, one of the dimensions of the emission region. Using $`R^{}=c\mathrm{\Delta }t_{\mathrm{obs}}D`$ and $`U_{\mathrm{phot}}^{}=L^{}/4\pi R_{}^{}{}_{}{}^{2}c`$ may then lead to over-estimation or under-estimtion of the jet-frame photon energy density by orders of magnitude. This is clearly of importance in any AGN model in which the low energy photons produced in the jet are targets for interaction of high energy particles or radiation, such as in SSC models and hadronic blazar models. Although not discussed in detail in the present paper, the escape of $`\gamma `$-rays from the emission region depends on the optical depth to photon-photon pair production interactions. This optical depth can be uncertain by orders of magnitude in the same way as the photon energy density, and will also depend on viewing angle. One must therefore be careful when using the observation of apparently unattenuated gamma-rays, and an observed variability time, to place limits on the Doppler factor.

The uncertainty in the jet-frame photon energy density discussed in this paper may also have implications for the high brightness temperature/Compton catastrophe problem of IDV sources. In this case, it is the energy density of target photons which limits the brightness temperature through the competition of inverse-Compton scattering with synchrotron radiation, and the target photon energy density may actually be lower than estimated if the emission region is non-spherical.

If the jet is bent or helical, or has some other favoured geometry (e.g. conical shocks) cusps in $`t^{}`$ vs. $`t_{\mathrm{obs}}`$, and/or a varying Doppler factor may cause narrow peaks in the observed lightcurve irrespective of other factors. Distinguishing between these cases from the observed lightcurve alone is likely to be difficult. One way of distinguishing whether a flare is due to (i) an emission region moving around a bent or helical path, or (ii) a shock is exciting a curved, conical or helical structure within a jet, is that the in the first case the flare is caused by a change in viewing angle with respect to the motion of the emission region leading to a change in Doppler factor, whereas in the second case the flare is due to a pile-up in observation times with no change in Doppler factor. Hence, in case (i) not only will the observed flux increase during a flare, but the photon energies also increase – increase in $`(\nu F_\nu )_{\mathrm{peak}}`$ by factor $`x^4`$ accompanied by shift in $`\nu _{\mathrm{peak}}`$ by factor $`x`$ ($`x`$ is ratio of final to initial Doppler factor). In case (ii), however, since there is no change in Doppler factor there should be no shift in $`\nu _{\mathrm{peak}}`$ accompanying an increase in $`(\nu F_\nu )_{\mathrm{peak}}`$.

Distinguishing between the excitation of a conical and a helical structure by a plane shock would be almost impossible – note the qualitative similarities between lightcurves depicted in Figs. 9(d) and 13(b). Relativistic jet simulations (Bowman 1994, Gomez et al. 1997) do show the presence of conical shocks, and these shocks appear after a large perturbation (e.g. from 4 to 10 in the simulations of Gomez et al. 1997) in Lorentz factor of the matter entering the jet through the nozzle. This can result in a sequence of quasi-stationary to superluminal reverse and forward conical shocks extending from the nozzle to the perturbation as it moves along the jet (Agudo et al. 2001). For the conical shock model (Salvati et al. 1998) discussed here to work, a subsequent perturbation would need to result in a plane shock, or plane thin perturbation of some kind, which could travel along the jet and excite the pre-existing conical shocks. As far as I am aware, whether or not this could occur has not been demonstrated. A similar uncertainty hangs over whether or not helical jet structures, which may themselves be shocks with a twisted ribbon topology, resulting from a perturbation entering the jet through the nozzle, could subsequently be excited by the passage of a plane shock or perturbation. Nevertheless, in both cases, if the viewing angle is favourable pile-ups in $`t_{\mathrm{obs}}`$, and hence flares, could occur simply as a result of the motion of the conical or helical patterns which may themselves be sites of enhanced emission. Note that in recent 3D relativistic jet simulations, the introduction of a 1 percent helical velocity perturbation at the nozzle results in a helical pattern propagating along the jet at nearly the beam speed (Aloy et al. 1999), and that the conical shocks resulting from a perturbation in jet Lorentz factor can range from being quasi-stationary to superluminal (Agudo et al. 2001).

In conclusion, in models for flaring in AGN in which the emission comes from a localized region (blob) co-moving with the jet, time variability is non-trivial to interpret in terms of emission region geometry and Doppler factor. A further complication is that flaring may arise instead due to curved or helical motion of a blob, even if the emission is constant in the istantaneous rest frame of the blob. In this case, apparent flaring is due to the change in viewing angle, and hence Doppler factor. Similarly, if the viewing angle is favourable, relativistic motion of curved or helical filaments or surfaces can lead to observation of flares. Excitation of curved or helical jet structures by shocks or perturbations can also lead to pile-ups in $`t_{\mathrm{obs}}`$, and hence large apparent increases in flux. Observations of time variability in AGN is therefore non-trivial to interpret and may lead to large systematic errors in estimated jet-frame photon energy density, Doppler factor and the physical parameters of the emission region.

## Acknowledgments

I thank Peter Biermann for helpful discussions and Alina Donea for a careful reading of the manuscript. I would like to thank the anonymous referees for their suggestions which have led to improvements to the paper. My research is supported by a grant from the Australian Research Council and a grant from the University of Adelaide.

## References

Agudo, I, Gomez, J.-L., Marti, J.-M., Ibanez, J.-M., Marscher, A.P., Alberdi, A., Aloy, M.-A., Hardee, P.E. 2001, ApJ, 549, L183

Aloy, M.-A., Ibanez, J.-M., Marti, J.M., Gomez, J.-L., Muller, E. 1999, ApJ, 523, L125

Bowman, M. 1994, MNRAS, 269, 137

Bednarek, W. & Protheroe, R.J. 1997, MNRAS, 292, 646

Bednarek, W. & Protheroe, R.J. 1999, MNRAS, 310, 577

Camenzind, M. 1986, A&A, 156, 136

Catanese, M. et al. 1998, ApJ, 501, 616

Chadwick, P.M. et al. 1999, ApJ, 513, 161

Despringre, V. & Fraix-Burnet, D. 1997, A&A, 320, 26

Gaidos, J.A. et al. 1996, Nature, 383, 319

Giebels, B. et al. 2002, ApJ, 571, 763

Gomez, L.L., Marti, J.M., Marscher, A.P., Ibanez, J.M, & Alberdi, A., 1997, ApJ, 482, L33

Hardee, P.E. 2000, ApJ, 533, 176

Hartman, R.C. et al. 1999, ApJSuppl, 123, 79

Hayashida, N. et al. 1998, ApJ, 504, L71

Kedziora-Chudczer, L., et al. 1997, ApJ, 490, L9

Kardashev, N.S., 2000, Astron. Reports, 44, 719

Kellerman, K.I., & Pauliny-Toth, I.I.K., 1969, ApJ, 155, L71

Kifune, T. 2002, PASA, 19, 1-4

Kniffen, D.A. et al. 1993, ApJ, 411, 133.

Konopelko, A. et al. 1999, Astropart. Phys., 11, 135

Lind, K.R & Blandford, R.D. 1985, ApJ, 295, 358

Mannheim, K. & Biermann, P.L., 1989, A&A, 221, 211

Mannheim, K., Protheroe, R.J. & Rachen, J., 2001, Phys. Rev. D 63, 023003.

Mücke, A., & Protheroe, R.J., 2001, Astropart. Phys., 15 121.

Mücke, A., Protheroe, R.J., Engel, R., Rachen, J. & Stanev, T., 2002, Astropart. Phys., Astropart. Phys., submitted.

Marcowith, A., Henri G. & Pelletier G., 1995, MNRAS, 277, 681.

Mendoza, S. & Longair, M.S. 2001, MNRAS, 324, 149

von Montigny C. et al., 1995, ApJ, 440, 525

Mukherjee R. et al., 1997, ApJ, 490, 116.

Protheroe, R.J. 1997, in Accretion Phenomena and Related Outflows, IAU Colloq. 163. ASP Conf. Ser. Vol. 121, ed. D.T. Wickramasinghe et al., p. 585

Protheroe, R.J. et al., 1998, Invited, Rapporteur, and Highlight Papers of the 25th Int. Cosmic Ray Conf. (Durban), eds. M.S. Potgieter et al., pub. World Scientific (Singapore), p. 317

Protheroe, R.J. 2002, MNRAS, submitted.

Punch M. et al., 1992, Nature, 358, 477

Qian S.-J. et al., 1992, Chinese A&A, 16, 137

Quinn J. et al. 1996, ApJ, 456, L83

Rantakyro F. T. et al, 1998, A&ASuppl, 131, 451

Rieger, E.M. & Mannheim, K. 2000, A&A, 359, 948

Rosen A., 1990, ApJ, 359, 296

Salonen, E., Terasranta, H., Urpo, S., Tiuri, M., Moiseev, I.G., Nesterov, N.S., Valtaoja, E., Haarala, S., Lehto, H., Valtaoja, L., Teerikorpi, P., Valotnen, M. 1987, A&ASuppl, 70, 409

Salvati, M., Spada, M. & Pacini, F 1998, ApJ, 495, L19

Schramm K.-J. et al, 1993, A&A, 278, 391

Spada, M., Salvati, M. & Pacini, F 1999, ApJ, 511, 136

Steffen W., Zensus, J.A., Krichbaum, T.P., Witzel, A., Qian, S.J., 1995, A&A, 302, 335.

Tateyama C. E. et al, 1998, ApJ, 500, 810

Terasranta, H., Tornikoski, M., Valtaoja, E., Urpo, S., Nesterov, N., Lainela, M., Kotilainen, J., Wiren, S., Laine, S., Nilsson, K., Valtonen, L. 1992, A&ASuppl, 94, 121

Urry, C.M. & Padovani, P. 1995, PASP, 107, 803

Villata, M. & Raiteri, C.M., 1999, A&A, 347, 30

Wagner, S.J. & Witzel, A., 1995, ARA&A, 33, 163

Walker, M.A., 1998, MNRAS, 294, 307


text:
# UWThPh-2002-29 Global analysis of neutrino oscillation data in four-neutrino schemes

## 1 INTRODUCTION

The neutrino oscillation interpretations of the solar  and atmospheric  neutrino data and the LSND experiment  require three neutrino mass-squared differences of different orders of magnitude. Since it is not possible to obtain this within the Standard Model framework of three active neutrinos it has been proposed to introduce a light sterile neutrino  to reconcile all the experimental hints for neutrino oscillations. Here we present an analysis of the global neutrino oscillation data in terms of four-neutrino mass schemes, including data from solar and atmospheric neutrino experiments, the LSND experiment, as well as data from short-baseline (SBL) experiments  and long-baseline reactor experiments  reporting no evidence for oscillations. This analysis updates the work presented at the meeting and is based on the data of summer 2002. We find that for all possible types of four-neutrino schemes different sub-sets of the data are in serious disagreement and hence, four-neutrino oscillations do not provide a satisfactory description of the global oscillation data including LSND. The details of our calculations can be found in Refs. .

Four-neutrino mass schemes are usually divided into the two classes (3+1) and (2+2), as illustrated in Fig. 1. We note that (3+1) mass spectra include the three-active neutrino scenario as limiting case. In this case solar and atmospheric neutrino oscillations are explained by active neutrino oscillations, with mass-squared differences $`\mathrm{\Delta }m^2sol`$ and $`\mathrm{\Delta }m^2atm`$, and the fourth neutrino state gets completely decoupled. We will refer to such limiting scenario as (3+0). In contrast, the (2+2) spectrum is intrinsically different, as there must be a significant contribution of the sterile neutrino either in solar or in atmospheric neutrino oscillations or in both.

## 2 NOTATIONS AND APPROXIMATIONS

Neglecting CP violation, in general neutrino oscillations in four-neutrino schemes are described by 9 parameters: 3 mass-squared differences and 6 mixing angles in the unitary lepton mixing matrix. Here we use a parameterization introduced in Ref. , which is based on physically relevant quantities: the 6 parameters $`\mathrm{\Delta }m^2sol`$, $`\theta sol`$, $`\mathrm{\Delta }m^2atm`$, $`\theta atm`$, $`\mathrm{\Delta }m^2lsnd`$, $`\theta lsnd`$ are similar to the two-neutrino mass-squared differences and mixing angles and are directly related to the oscillations in solar, atmospheric and the LSND experiments. For the remaining 3 parameters we use $`\eta _s,\eta _e`$ and $`d_\mu `$. Here, $`\eta _s(\eta _e)`$ is the fraction of $`\nu _s(\nu _e)`$ participating in solar oscillations, and ($`1d_\mu `$) is the fraction of $`\nu _\mu `$ participating in oscillations with $`\mathrm{\Delta }m^2atm`$ (for exact definitions see Ref. ). For the analysis we adopt the following approximations:

* We make use of the hierarchy

$$\mathrm{\Delta }m^2sol\mathrm{\Delta }m^2atm\mathrm{\Delta }m^2lsnd.$$
(1)
This means that for each data set we consider only one mass-squared difference, the other two are set either to zero or to infinity.
* In the analyses of solar and atmospheric data (but not for SBL data) we set $`\eta _e=1`$, which is justified because of strong constraints from reactor experiments .

Due to these approximations the parameter structure of the four-neutrino analysis gets rather simple. The parameter dependence of the different data sets solar, atmospheric and SBL is illustrated in Fig. 2. We see that only $`\eta _s`$ links solar and atmospheric data and $`d_\mu `$ links atmospheric and SBL data. All the other parameters are “private” to one data set.

## 3 (2+2): RULED OUT BY SOLAR AND ATMOSPHERIC DATA

The strong preference of oscillations into active neutrinos in solar and atmospheric oscillations leads to a direct conflict in (2+2) oscillation schemes. We will now show that thanks to the new SNO solar neutrino data  and the improved SK statistic on atmospheric neutrinos  the tension in the data has become so strong that (2+2) oscillation schemes are essentially ruled out.111Details of our analyses of the solar and atmospheric neutrino data can be found in Ref. . For an earlier four-neutrino analysis of solar and atmospheric data see Ref. . Indeed, latest solar neutrino data lead to the bound $`\eta _s0.45`$ at 99% CL, where $`\eta _s`$ is the parameter describing the fraction of the sterile neutrino participating in solar neutrino oscillations. In contrast, in (2+2) schemes atmospheric data imply $`\eta _s0.65`$ at 99% CL, in clear disagreement with the bound from solar data.

In Fig. 3 we show the $`\chi ^2`$ for solar data and for atmospheric combined with SBL data as a function of $`\eta _s`$. Furthermore, we show the $`\chi ^2`$ of the global data defined by

$$\overline{\chi }^2(\eta _s)\mathrm{\Delta }\chi ^2sol(\eta _s)+\mathrm{\Delta }\chi _{\mathrm{atm}+\mathrm{sbl}}^2(\eta _s).$$
(2)
From the figure we find that only if we take both data sets at the 99.95% CL a value of $`\eta _s`$ exists, which is contained in the allowed regions of both sets. This follows from the $`\chi ^2`$-value $`\chi _{\mathrm{PC}}^2=12.2`$ shown in the figure. In Ref.  we have proposed a statistical method to evaluate the disagreement of different data sets in global analyses. The parameter goodness of fit (PG) makes use of the $`\overline{\chi }^2`$ defined in Eq. (2). This criterion evaluates the GOF of the combination of data sets, without being diluted by a large number of data points, as it happens for the usual GOF criterion (for details see Ref. ). We find $`\chi _{\mathrm{PG}}^2\overline{\chi }_{\mathrm{min}}^2=23.5`$, leading to the marginal PG of $`1.3\times 10^6`$. We conclude that (2+2) oscillation schemes are highly disfavoured by the disagreement between the latest solar and atmospheric neutrino data. This is a very robust result, independent of whether LSND is confirmed or disproved.

## 4 (3+1): STRONGLY DISFAVOURED BY SBL DATA

It is known for a long time  that (3+1) mass schemes are disfavoured by the comparison of SBL disappearance data  with the LSND result. In Ref.  we have calculated an upper bound on the LSND oscillation amplitude $`\mathrm{sin}^22\theta lsnd`$ resulting from SBL and atmospheric neutrino data. From Fig. 4 we see that this bound is incompatible with the signal observed in LSND at the 95% CL. Only marginal overlap regions exist between the bound and LSND if both are taken at 99% CL. An analysis in terms of the parameter goodness of fit  shows that for most values of $`\mathrm{\Delta }m^2lsnd`$ SBL and atmospheric data are compatible with LSND only at more than $`3\sigma `$, with one exception around $`\mathrm{\Delta }m^2lsnd6`$ eV2, where the PG reaches 1%. These results show that (3+1) schemes are strongly disfavoured by SBL disappearance data.

## 5 COMPARING (3+1), (2+2) AND (3+0) HYPOTHESES

With the methods developed in Ref.  we are able to perform a global fit to the oscillation data in the four-neutrino framework. This approach allows to statistically compare the different hypotheses. Let us first evaluate the GOF of (3+1) and (2+2) spectra with the help of the PG method described in Ref. . We divide the global oscillation data in the 4 data sets SOL, ATM, LSND and NEV, where NEV contains the experiments KARMEN , CDHS , Bugey  and CHOOZ/Palo Verde , reporting no evidence for oscillations. Following Ref.  we consider

$$\begin{array}{ccc}\overline{\chi }^2& =& \mathrm{\Delta }\chi ^2sol(\theta sol,\mathrm{\Delta }m^2sol,\eta _s)\hfill \\ & +& \mathrm{\Delta }\chi ^2atm(\theta atm,\mathrm{\Delta }m^2atm,\eta _s,d_\mu )\hfill \\ & +& \mathrm{\Delta }\chi ^2nev(\theta lsnd,\mathrm{\Delta }m^2lsnd,d_\mu ,\eta _e)\hfill \\ & +& \mathrm{\Delta }\chi ^2lsnd(\theta lsnd,\mathrm{\Delta }m^2lsnd),\hfill \end{array}$$
(3)
where $`\mathrm{\Delta }\chi _X^2=\chi _X^2(\chi _X^2)_{\mathrm{min}}`$ ($`X`$ = SOL, ATM, NEV, LSND). In Tab. 1 we show the contributions of the 4 data sets to $`\chi _{\mathrm{PG}}^2\overline{\chi }_{\mathrm{min}}^2`$ for (3+1) and (2+2) oscillation schemes. As expected we observe that in (3+1) schemes the main contribution comes from SBL data due to the tension between LSND and NEV data in these schemes. For (2+2) oscillation schemes a large part of $`\chi _{\mathrm{PG}}^2`$ comes from solar and atmospheric data, however, also SBL data contributes significantly. This comes mainly from the tension between LSND and KARMEN , which does not depend on the mass scheme and, hence, also contributes in the case of (2+2). Therefore, the values of $`\chi _{\mathrm{PG}}^2`$ in Tab. 1 for (2+2) schemes are higher than the one given in Sec. 3, where the tension in SBL data is not included.

The parameter goodness of fit is now obtained by evaluating $`\chi _{\mathrm{PG}}^2`$ for 4 DOF . This number of degrees of freedom corresponds to the 4 parameters $`\eta _s,d_\mu ,\theta lsnd,\mathrm{\Delta }m^2lsnd`$ describing the coupling of the different data sets (see Eq. (3)). The best GOF is obtained in the (3+1) case. However, even in this best case the PG is only 0.56%. The PG of $`1.6\times 10^6`$ for (2+2) schemes shows that these mass schemes are essentially ruled out by the disagreement between the individual data sets.

Although we have seen that none of the four-neutrino mass schemes can provide a reasonable good fit to the global oscillation data including LSND, it might be interesting to consider the relative status of the three hypotheses (3+1), (2+2) and the three-active neutrino scenario (3+0). This can be done by comparing the $`\chi ^2`$ values of the best fit point (which is in the (3+1) scheme) with the one corresponding to (2+2) and (3+0). First we observe that (2+2) schemes are strongly disfavoured with respect to (3+1) with a $`\mathrm{\Delta }\chi ^2=17.8`$. The reason for the big change with respect to the value of $`\mathrm{\Delta }\chi ^2=3.7`$ found in Ref.  is the improved sensitivity of solar (SNO NC) and atmospheric (SK 1489-days) data to a sterile component. With this new data now (3+1) schemes are clearly preferred over (2+2): for 4 DOF a $`\mathrm{\Delta }\chi ^2=17.8`$ implies that (2+2) is ruled out at 99.87% CL with respect to (3+1). Further we find that (2+2) is only slightly better than (3+0), which is disfavoured with a $`\mathrm{\Delta }\chi ^2=20.0`$ with respect to (3+1).

## 6 CONCLUSIONS

Using latest solar and atmospheric data, we re-analyze the four-neutrino description of current global neutrino oscillation data, including the LSND evidence for oscillations. The higher degree of rejection for non-active solar and atmospheric oscillation solutions implied by the SNO neutral current result as well as by the latest 1489-day SK atmospheric neutrino data allows us to rule out (2+2) oscillation schemes. Using an improved goodness of fit method especially sensitive to the combination of data sets we obtain a GOF of only $`1.6\times 10^6`$ for (2+2) schemes. Further, we find that also (3+1) schemes are strongly disfavoured by the data. The disagreement between negative short-baseline experiments and LSND in (3+1) schemes imply a GOF of $`5.6\times 10^3`$. This leads to the conclusion that all four-neutrino descriptions of the LSND anomaly, both in (2+2) as well as (3+1) realizations, are highly disfavoured. Our analysis brings the LSND hint to a more puzzling status, and the situation will become even more puzzling if LSND should be confirmed by the up-coming MiniBooNE experiment .