format: text

text:
# Search for Relativistic Curvature Effects in Gamma-Ray Burst Pulses

## Introduction

The temporal structure of gamma-ray bursts (GRBs) varies drastically, with no apparent pattern among bursts. Of the $`>`$2704 GRBs detected by the BATSE instrument onboard the Compton Gamma Ray Observatory (CGRO) (Fishman et al., 1994), less than 10 $`\%`$ of the detected light-curve profiles can be categorized as being similar in overall shape or morphology, in that they comprise coherent structures, or pulses of radiation. It is generally believed that these pulses represent the fundamental constituent of GRB time profiles (light curves), and appear as asymmetric pulses with a fast rise and a slower decay, often denoted FRED ”Fast Rise and Exponential Decay” (see Figure 1). Each pulse is assumed to be associated with a separate emission episode, with complex GRBs being superpositions of several such episodes, of varying amplitude and intensity. Although the duration and amplitude of the FRED pulses vary considerably, the shape is the only recurring pattern that can be distinguished among the vast range of complex GRB light curves that have been observed.

The observed $`\gamma `$-ray pulses are believed to be produced in a highly relativistic outflow, an expanding and collimated fireball, based on the large energies and the short time scales involved. There are several possible mechanisms underlying the emission, although the commonly assumed scenario is that individual pulses are created, when shocks internal to the relativistic outflow, drain the kinetic energy and accelerate leptons which radiate. In this paper we base our discussion on this standard, fireball model, in which the $`\gamma `$-rays arise from the internal shocks at a distance of $`R10^{13}10^{17}`$ cm from the initial source (Piran, 1999). In the context of the fireball model, the episodic nature of the outflow causes inhomogeneities in the wind to collide and thus create the shocks. The dominant emission mechanisms are most probably non-thermal synchrotron (Tavani, 1996; Lloyd & Petrosian, 2000)) and/or inverse Compton emission (Liang, Kocevski, Boettcher, 2003), but there have been other suggestions, for instance, thermal, saturated Comptonization (Liang, 1997). The simplest scenario here is to assume an impulsive heating of the leptons and a subsequent cooling and emission. Therefore, the rise phase of the pulse is attributed to the energization of the shell and the decay phase reflects the cooling of the energized particles.

As many authors have pointed out, this cooling interpretation does not work well when applied to non-thermal synchrotron and/or inverse Compton emission because the radiative cooling timescales alone are typically much too short to explain the pulse durations. For example, if we interpret the average break energy of a GRB spectrum ($`E_{pk}`$) as the characteristic energy of synchrotron self-absorption, then the resulting magnetic field must be extremely high, about $`1\times 10^7`$ to $`1\times 10^8`$ Gauss. Similar estimates are obtained if we assume equipartition conditions between the lepton energy and the magnetic field. In either case, such a high $`B`$ field would create a synchrotron cooling timescale of the order of $`10^5`$ seconds in the comoving frame (Wu & Fenimore, 1996). One resolution to this problem is to introduce relativistic effects. If the relativistic fireball expands with a Lorentz factor of $`\mathrm{\Gamma }`$ $`>`$ 100, then the geometry of the shell would make radiation emitted off the line of sight delayed and affected by a varying Doppler boost (see Figure 2). This would cause the GRB spectra to evolve to lower energies and produce decay profiles much longer than the microscopic cooling timescale. These relativistic curvature effects would produce a signature decay profile that can be obtained analytically and searched for in a sample of GRB pulses. Thus it is of great importance to characterize the individual pulse profiles within GRBs and study their parameter distributions.

Several investigations along the lines of pulse modelling have perviously been made. Initially, the pulse profiles were modelled by ”stretched exponential” functions, both for the rise phase and for the decay phase (Norris et al., 1996; Lee et al., 2000a, b):

$$F(t)=F_0e^{(|tt_{\mathrm{max}}|/\tau _{\mathrm{r},\mathrm{d}})^\nu },$$
(1)
where $`t_{\mathrm{max}}`$ is the time of the maximum flux, $`F_0`$, of the pulse, $`\tau _{\mathrm{r},\mathrm{d}}`$ are the time constants for the rise and the decay phases, respectively, and $`\nu `$ is the peakedness parameter111For $`\nu >1`$, equation (1) is, strictly speaking, a compressed exponential.. Such a function is very flexible which makes it possible to describe the whole shape of most pulses, and to quantify the characteristics of the pulses for a statistical analysis, most notably their location, amplitude, width, rise phase, transition phase, and decay phase. Norris et al. (1996) studied a sample of bursts observed by BATSE and found that the decay generally lies between a pure exponential ($`\nu =1`$) and a Gaussian ($`\nu =2`$). Lee et al. (2000a, b) studied approximately 2500 pulse structures, in individual energy channels, using the high time resolution BATSE TTS data type. They confirmed the general behavior, namely, that pulses tend to have shorter rise times than decay times.

The stretched exponential was introduced because of its very flexible nature, although recently Ryde & Svensson (2000, 2002) have proposed that the decay phase would be better described by a power law, which in terms of energy flux is :

$$F(t)=F_0\left(1+\frac{t}{T(d1)}\right)^d,$$
(2)
See §2 for the derivation of equation (2) and definition of $`d`$. The motivation for this type of shape is entirely based on empirical relations describing the evolution of the GRB spectra during the decay phase of individual pulses.

Ryde & Petrosian (2002) have recently shown that the form expressed in equation (2) can be produced through simple relativistic kinematics when applied to a spherical shell expanding at extreme relativistic velocity. The curvature of a relativistic fireball would make the photons emitted off the line of sight (LOS) delayed and affected by a varying Doppler boost due to the increasing angle at which the photons were emitted with respect to the observer (See Figure 2 for an illustration). They show that this Doppler boosting of off-axis photons can reproduce the two well known empirical correlations observed in the GRB spectra, namely the hardness-intensity and hardness-fluence correlations (HIC and HFC, respectively), in a manner that is highly predictive of the resulting pulse profile.222 See Kobayashi, Ryde, & MacFadyen (2001) for the effects on the morphology of a complex GRB light curve. Following Ryde & Petrosian (2002), we can derive the expected emission profile from a spherical shell that radiates for an infinitesimal period of time at a peak energy $`E_{pk}^{}`$ in the comoving frame. The Lorentz boosting factor for transformations from the comoving frame to the observers frame of photons emitted from different locations on the surface of a spherical shell is given by

$$𝒟(\mathrm{\Gamma },\mu )=\frac{1}{\mathrm{\Gamma }(1\beta \mu )}=\mathrm{\Gamma }(1+\beta \mu ^{}),$$
(3)
where the angle $`\theta \mathrm{arccos}\mu `$ shown in figure 2. If we define the point where the flow velocity is parallel to the line of sight (LOS) to be $`\theta =0`$, then the difference in light travel time between photons emitted along the LOS and photons emitted at an angle $`\theta `$ is given by $`\mathrm{\Delta }t=R_0(1\mu )/c`$, which gives $`\mu =1c\mathrm{\Delta }t/R_0`$. Substituting this into the Lorentz boost factor, we find

$$𝒟(\mathrm{\Gamma },\mathrm{\Delta }t)=\frac{1}{\mathrm{\Gamma }(1\beta +\beta c\mathrm{\Delta }t/R_0)}.$$
(4)
For extremely relativistic outflows $`\beta =1(2\mathrm{\Gamma }^2)^1`$, and the boost factor becomes

$$𝒟(\mathrm{\Gamma },\mathrm{\Delta }t)=2\mathrm{\Gamma }(1+\mathrm{\Delta }t/\tau _{\mathrm{ang}})^1,\tau _{\mathrm{ang}}R_0/(2\mathrm{\Gamma }^2c).$$
(5)
The outcome of such a Doppler profile is that if the emitted spectra from different parts of the shell are identical, then the observed spectra will be gradually redshifted in time by a factor of $`𝒟`$ as photons from different parts of the shell are received by the observer. The peak energy of the spectra will be observed to evolve as

$$E_{\mathrm{pk}}(t)=E_{\mathrm{pk}}^{}𝒟=\frac{E_{\mathrm{pk},0}}{(1+t/\tau _{\mathrm{ang}})}$$
(6)
Where $`E_{\mathrm{pk}}^{}`$ is the peak energy in the comoving frame. By similar arguments, Ryde & Petrosian (2002) show that the resulting bolometric energy flux should evolve as

$$F_{\mathrm{bol}}(t)=F_0𝒟^2=\frac{F_0}{(1+t/\tau _{\mathrm{ang}})^2}$$
(7)
Therefore, the decay of the resulting light curve, produced purely by relativistic effects, should exhibit a distinct profile, given by equation (7). This result becomes more complicated when one considers the three primary timescales relevant for the shape of the pulse: the shell crossing time, the angular spreading time, and the cooling time. The first two timescales are comparable and dominate over the third and compete with each other to produce the observed results. The authors find that for the case of ”pure” spherical curvature, where the curvature timescale dominates over the intrinsic light-curve profile (the shell crossing time), the convolved light curves indeed reach the behavior of equation (7) fairly quickly. However, if the crossing timescale and the curvature timescale are comparable, then the resulting light curve will have some convolved shape and the index value of 2 will only be manifested at late times in the pulse decay.

The Ryde & Petrosian (2002) treatment assumes that the variations in the low- and high-energy power law indices ($`\alpha `$ and $`\beta `$) that are commonly used to describe the GRB spectra tend to be small over the length of the pulse, which may not necessarily be the case. In order to account for this, Fenimore et al. (1996) and Fenimore & Sumner (1997) used a photon number spectrum which they described by a single power law index equal to the average value of $`(\alpha +\beta )/2`$. They examined the resulting light curve and spectral evolution that would be expected if the FRED profile were due primarily to relativistic effects. They conclude that the decay phase of FRED profiles should be independent of intrinsic spectral variations and should scale roughly as $`T^{\alpha 1}`$.

Therefore, one of the primary goals of this paper is to model the FRED light curves using an analytic function derived from physical first principles and relations describing the spectral evolution of GRBs to obtain parameters that uniquely quantify the shape of smooth pulses and compare those values to the predicted signatures of relativistic kinematics. Preliminary results of this analysis have been presented in Kocevski & Liang (2002b) and Ryde, Kocevski, & Liang (2003), where we derived analytic models and tested for a small number of pulses. This paper expands on these previous results to include a much larger sample and a more detailed analysis of the distribution of pulse attributes. In $`\mathrm{\S }1`$, we review the empirical hardness-intensity and hardness-fluence correlations which are crucial to the derivation of an analytic pulse model. In $`\mathrm{\S }2`$, we derive an analytic function based on the above-mentioned spectral correlations. In $`\mathrm{\S }3`$, we define a sample of FRED pulses observed by CGRO to which we fit our function to obtain a distribution of model parameters, most notably the decay power law index and pulse asymmetries. In $`\mathrm{\S }4`$, we discuss our measured parameter distributions and compare them to the values predicted by relativistic geometry. We discuss several mechanisms that can skew the distributions away from the expected values.

## 1 Spectral Correlations

Gamma-ray burst continuum spectra have a well-known evolution as the burst proceeds, both over the entire, often complex light curves, and over individual pulses. The latter evolution is often more spectacular. This evolution is generally characterized by an overall softening of the spectra to lower energies with time. Two specific correlations between observable quantities have been found that describe this softening in a qualitative manner. The first is the hardness-intensity correlation or HIC, which relates the instantaneous hardness of the spectra and the instantaneous energy flux $`F_E`$, within individual pulses. For the decay phase of a pulse the most common behavior of the HIC is

$$F_E=F_0(E_{\mathrm{pk}}/E_{\mathrm{pk},0})^\eta ,$$
(8)
where $`E_{\mathrm{pk},0}`$ and $`F_0`$ are the initial values of the peak of the $`EF_E`$ spectrum and the energy flux at the beginning of the decay phase of each pulse, respectively. $`\eta `$ is the power law index333We follow the notation introduced in Ryde & Petrosian (2002). Note that Golenetskii et al. (1983) used $`\gamma `$ for the power law index, which unfortunately is also a notation for the Lorentz factor. This is a general form for non-thermal energy flux spectra that is largely independent of the emission mechanism which gives $`F`$ an implicit time dependence because the peak energy evolves with time. The original study by Golenetskii et al. (1983) found the power law index to vary between $`1.51.7`$ over the whole GRB. Moreover, Borgonovo & Ryde (2001) studied a sample of 82 GRB pulse decays and found them to be consistent with a power law HIC in, at least, 57% of the cases and for these found $`\eta =2.0\pm 0.7`$. Several physical emission processes are described by such a relation, with the most prominent example being single particle synchrotron emission with a constant magnetic field and a variable electron Lorentz factor $`\mathrm{\Gamma }`$. In this particular case, $`E_{pk}\mathrm{\Gamma }^2`$ while $`F\mathrm{\Gamma }^2`$, so the HIC index is $`\eta =1`$.

The second correlation is the hardness-fluence correlation or HFC (Liang & Kargatis, 1996) which describes the observation that the instantaneous hardness of the spectra decays exponentially as a function of the time-integrated flux, or fluence, of the burst:

$$E_{pk}=E_0e^{\mathrm{\Phi }/\mathrm{\Phi }_0},$$
(9)
where $`\mathrm{\Phi }(t)`$ is the photon fluence integrated from the start of the burst and $`\mathrm{\Phi }_0`$ is the exponential decay constant. Liang & Kargatis (1996) found the HFC to be valid for 35 of 37 smooth GRB pulses. These findings were later confirmed by Crider et al. (1999) who studied a sample of 41 pulses in 26 GRBs and found that the spectral evolution in the majority of the pulses could be described by the HFC. The physical interpretation of the HFC is not immediately clear, but upon differentiation of equation (9), it can be seen that it simply states that the rate of change in the hardness is proportional to the luminosity of the radiating medium (or, equivalently, to the energy density);

$$dE_{pk}/dt=\frac{F_NE_{pk}}{\mathrm{\Phi }_0}\frac{F_E}{\mathrm{\Phi }_0}.$$
(10)
Here the last equivalency is only approximately true. The energy flux is actually defined as the integral of the photon energy times the photon spectrum, $`F_E=EN(E)𝑑E`$, where in this case the integral would be carried over the BATSE energy band. We make the assumption that if $`E_{pk}`$ falls well within that band, then the energy flux is approximately $`E_{pk}`$ times the instantaneous photon flux, $`F_NE_{pk}F_E`$. Since $`E_{pk}`$ is defined to be the peak of the $`EF_E`$ spectrum, it represents the energy of most of the arriving photons and therefore this equivalency is robust as long as $`E_{pk}`$ does not evolve outside of the BATSE window. We also make the assumption that the evolution of the low-energy power law index $`\alpha `$ should not effect our results except for pulses with an extreme amount of evolution. This assumption is supported by the findings of Crider et al. (1999) and Borgonovo & Ryde (2001) who avoided making the assumption that $`F_NE_{pk}F_E`$ by directly using photon flux to test the HFC and HIC relations. Using independent methods, they show the approximation to be valid for a large sample of GRBs.

Liang and Kargatis originally proposed that the HFC could be the result of a confined plasma with a fixed number of particles cooling via $`\gamma `$-radiation. This type of exponential decay of the break energy with photon fluence that is seen in the HFC is expected if the average energy of the emitted photons is directly proportional to the average emitting particle energy such as in thermal bremsstrahlung or multiple Compton scattering, although this interpretation is not unique.

## 2 Analytic Light Curve Profile

We now follow Ryde & Svensson (2000) to find an analytical description of the energy-flux decay-profile. Note that in their original description they used the photon flux, and we will here reformulate their results in terms of the energy flux instead, which makes a physical interpretation easier. Combining the empirical relations, the HIC and HFC, gives the following differential equation governing the spectral evolution

$$\dot{E}=\frac{F_0}{\varphi _0E_0^\eta }E^\eta $$
(11)
The solution gives an expression for the energy flux during the decay phase of a pulse

$$F(t)=\{\begin{array}{cc}F_0\left(1+\frac{(\eta 1)t}{T}\right)^{\eta /(\eta 1)}\hfill & \text{if }\eta 1\text{ }\hfill \\ F_0e^{t/T}\hfill & \text{if }\eta =1\text{ ,}\hfill \end{array}$$
(12)
and correspondingly the peak of the $`EF_E`$ spectra follow

$$E_{\mathrm{pk}}(t)=\{\begin{array}{cc}E_{\mathrm{pk},0}\left(1+\frac{(\eta 1)t}{T}\right)^{1/(\eta 1)}\hfill & \text{if }\eta 1\text{ }\hfill \\ E_{\mathrm{pk},0}e^{t/T}\hfill & \text{if }\eta =1\text{ ,}\hfill \end{array}$$
(13)
where $`T\mathrm{\Phi }_0E_{\mathrm{pk},0}/F_0`$. Note that this gives the possibility to measure $`\mathrm{\Phi }_0`$ directly from the light curve with the knowledge of $`E_{\mathrm{pk},0}`$. Introducing $`d\eta /(\eta 1)`$ ($`d`$ as in the asymptotic decay of the energy flux) we can describe the $`E_{pk}(t)`$ and the $`F(t)`$ decays as

$$F(t)=F_0\left(1+\frac{t}{T(d1)}\right)^d$$
(14)
$$E_{\mathrm{pk}}(t)=E_{\mathrm{pk},0}\left(1+\frac{t}{T(d1)}\right)^{1d}$$
(15)
for the $`\eta 1`$ case. Note that this representation is completely model-independent and is based solely on the empirical HIC and HFC relations that describe GRB’s spectral evolution. As will be shown in the following sections, specific models, such as relativistic kinematics, make specific predictions regarding the value of the HIC power-law index and hence the value of decay constant $`d`$.

### 2.1 Inclusion of the Rise Phase

The analytic profile derived above describes only the decay phase of a pulse, since that is the portion of the GRB time history for which the HIC is often a power law (eq. ) and the HFC follows equation (4). This introduces an ambiguity when fitting equations (14) and (15) to the data, since the power law index $`d`$ and the time constant $`T`$ are coupled and are therefore not well constrained by fitting. Furthermore, choosing the moment when the pulse changes from the rise to the decay phase is somewhat subjective. In order to eliminate these ambiguities, a rise phase must be supplemented to equations (9) and (10) in order to produce a complete description of the FRED profile. This also allows for additional parameters such as the pulse rise time and asymmetry values to be measured.

We first note the fact that in most physical models both the peak of the energy spectrum and the luminosity are proportional to the random Lorentz factor of the shocked electrons to some power (Rybicki & Lightman, 1979):

$$E_{\mathrm{pk}}(t)\mathrm{\Gamma }_\mathrm{r}^a(t)g(t),$$
(16)
and

$$F(t)\mathrm{\Gamma }_\mathrm{r}^b(t)h(t).$$
(17)
The functions $`g(t)`$ and $`h(t)`$ parameterize the unknown time dependencies on particle densities, optical depth, magnetic field and kinematics, etc. The correlation between hardness and the energy flux can thus be described as

$$F(t)E_{\mathrm{pk}}^\eta ^{}f(t)$$
(18)
with $`\eta ^{}=b/a`$ and $`f(t)=h(t)g(t)^\eta ^{}`$. During the decay phase of a pulse the power law relation dominates the HIC, which has been elucidated in previous studies. However, the HIC during the rise phase will be dominated by the unknown function $`f(t)`$.

We will now assume two simple representations of $`f(t)`$ and recapitulate the analytical discussion in the previous section to find a description for the shape of the entire pulse.

#### 2.1.1 Power-Law Rise

The simplest assumption one can make is that the rise phase is dominated by a power-law increase of flux. This leads to a simple prescription of the function

$$f(t)t^r$$
(19)
where the power-law index, $`r`$, stands for the rise. In the case of optically-thin synchrotron emission the energy flux is given by $`F_EE_{\mathrm{pk}}\tau _TB`$, so the HIC exponent in equation (8) is unity, $`\tau _T`$ is the Thompson depth, and $`B`$ is the magnetic field. In this case the explicit time dependence arises from the evolution of the Thompson depth and the magnetic field strength. Our initial assumption gives rise to a general form that should also emerge when considering other non-thermal emission mechanisms.

To fully model the FRED light curve we have to define the time at which the light curve peaks as $`t=t_{\mathrm{max}}`$ or $`t_\mathrm{m}`$ for brevity. This allows us to define several new parameters, $`E_{\mathrm{pk}}(t=0)=E_0`$ and $`E_{\mathrm{pk}}(t=t_\mathrm{m})=E_\mathrm{m}`$ and $`F(t=t_\mathrm{m})=F_\mathrm{m}`$ and that $`\tau _\mathrm{m}F_\mathrm{m}/E_\mathrm{m}\mathrm{\Phi }_0`$. The HIC can then be written as:

$$F=F_\mathrm{m}\left(\frac{E_{\mathrm{pk}}}{E_\mathrm{m}}\right)^\eta \left(\frac{t}{t_\mathrm{m}}\right)^r$$
(20)
which combined with the HFC now gives the differential equation

$$\dot{E_{\mathrm{pk}}}=\frac{F_\mathrm{m}}{\mathrm{\Phi }_0E_\mathrm{m}^\eta t_\mathrm{m}^r}E_{\mathrm{pk}}^\eta t^r$$
(21)
Now, we define the exponent $`d`$ (for finite $`r`$) to describe the asymptotic power law behavior of the light curve, $`F(t)t^d`$ as $`t\mathrm{}`$ which gives $`\eta =(d+r)/(d1)`$. We can now better define the peak of the pulse to be at $`t=t_\mathrm{m}`$ by letting $`\mathrm{d}F/\mathrm{d}t(t=t_\mathrm{m})=0`$ and get the condition that for $`t_\mathrm{m}0`$,$`\tau _\mathrm{m}=t_\mathrm{m}(d+r)/[(d1)r]`$. Therefore, for a simple power-law rise, the final pulse shape and $`E_{pk}`$ evolution can be described by:

$$F(t)=F_\mathrm{m}\left(\frac{t}{t_\mathrm{m}}\right)^r\left(\frac{d}{d+r}+\frac{r}{d+r}\left(\frac{t}{t_\mathrm{m}}\right)^{r+1}\right)^{\frac{r+d}{r+1}}$$
(22)
$$E_{\mathrm{pk}}(t)=E_\mathrm{m}\left(\frac{d}{d+r}+\frac{r}{d+r}\left(\frac{t}{t_\mathrm{m}}\right)^{r+1}\right)^{\frac{d1}{r+1}}$$
(23)
This four parameter model fully describes a single pulsed GRB event with direct measurements of an observed light curve yielding the values for $`F_\mathrm{m},t_\mathrm{m},r,d`$. This then gives $`\tau _\mathrm{m}`$, which in its turn, gives the product $`E_\mathrm{m}\mathrm{\Phi }_0`$. As $`E_\mathrm{m}`$ can be measured from the $`E_{\mathrm{pk}}`$ decay, $`\mathrm{\Phi }_0`$ can be deduced. Example profiles with varying decay parameters for the power-law rise-model are shown in Figure 3.

It must be noted that in this scenario (rise+decay), the HIC power law index, $`\eta ^{}`$, will correspond to both the implicit power law behavior, $`E^\eta `$ as well as the explicit behavior, $`t^r`$. Here $`\eta ^{}`$ denotes the power law in the second description. The two descriptions coincide as $`r`$ tends to $`0`$ and as $`E_{\mathrm{pk},0}E_\mathrm{m}`$ which gives

$$F=\frac{F_\mathrm{m}}{\left(1+\frac{t}{(d1)\tau _\mathrm{m}}\right)^d}$$
(24)
and thus $`T=\tau _\mathrm{m}`$ by identification. Note that $`t_\mathrm{m}`$ tends to $`0`$ as $`r`$ approaches $`0`$. It is also readily seen that the asymptotic behaviors as $`t`$ approaches $`\mathrm{}`$ are the same.

#### 2.1.2 Exponential Rise

It is reasonable that the rise phase is connected to some transient process, e.g., an initial decrease in optical depth, increase in the number of energized particles, or the merging (or crossing) of two shells. After this initial phase, the original decay behavior as described by equations (14) and (15) should emerge. We therefore try the following prescription

$$f(t)1e^{t/\tau _\mathrm{r}}$$
(25)
where the time constant $`\tau _\mathrm{r}`$ now represents the rise phase. This corresponds to $`h(t)g(t)^\eta `$ as $`t>>\tau _\mathrm{r}`$. We therefore have

$$F=F_\mathrm{m}\left(\frac{E_{\mathrm{pk}}}{E_\mathrm{m}}\right)^\eta \frac{1e^{t/\tau _\mathrm{r}}}{1e^{t_\mathrm{m}/\tau _\mathrm{r}}}$$
(26)
which gives the differential equation

$$\dot{E_{\mathrm{pk}}}=\frac{F_\mathrm{m}E_{\mathrm{pk}}^\eta (1e^{t/\tau _\mathrm{r}})}{\mathrm{\Phi }_0E_\mathrm{m}^\eta (1exp(t_\mathrm{m}/\tau _\mathrm{r}))}$$
(27)
In the same manner as before we define the decay index $`d`$ by requiring $`F(t)t^d`$ as $`t\mathrm{}`$. This gives $`\eta =d/(d1)`$. Solving the differential equation (27) and combining with equation (26) we find

$$F(t)=\frac{A_0(1e^{t/\tau _\mathrm{r}})}{(1+(t+\tau _\mathrm{r}e^{t/\tau _\mathrm{r}})/\tau _\mathrm{d})^d}$$
(28)
$$E_{\mathrm{pk}}(t)=\frac{A_1}{(1+(t+\tau _\mathrm{r}e^{t/\tau _\mathrm{r}})/\tau _\mathrm{d})^{(d1)}}$$
(29)
where $`A_0`$ and $`A_1`$ are analytical functions of $`d,\tau _\mathrm{d},\tau _\mathrm{r},\mathrm{\Phi }_0,F_\mathrm{m},E_\mathrm{m}`$ and $`t_\mathrm{m}`$. There is no analytical expression for the peak time $`t_\mathrm{m}`$, but it can be solved for numerically if the other parameters are known. Equation (28) describes the pulse by four parameters: $`A_0,\tau _\mathrm{r},\tau _\mathrm{d},`$ and $`d`$. This description coincides with the decay function (eqs.  and ) when $`\tau _\mathrm{r}0`$ or $`t\mathrm{}`$ with $`\tau _\mathrm{d}(d1)T`$ (and $`E_\mathrm{m},F_\mathrm{m}E_{\mathrm{pk},0},F_0`$). Note that in this description there are two explicit time constants $`\tau _\mathrm{r}`$ and $`\tau _\mathrm{d}`$, while in equation (22) the decay constant is an expression of $`r`$ and $`t_\mathrm{m}`$. For a sharp rise phase the explicit time dependence will disappear. This is, however, not the case for equation (22).

## 3 Data Analysis

For the purposes of this study, we utilize 64ms count data provided by the BATSE instrument onboard the CGRO spacecraft. These data were gathered by BATSE’s Large Area Detectors (LADs) which provide discriminator rates with 64 ms resolution from 2.048 s before the burst to several minutes after the trigger (Fishman et al., 1994). The discriminator rates are gathered in four broad energy channels covering approximately 25-50, 50-100, 100-300, and 300 to about 1800 keV allowing for excellent count statistics since the photons are collected over a wide energy band. For our analysis, we combine the data from all four channels in order to study the ”bolometric” light curve profile. We chose to study the count flux light curves rather then directly using the energy flux because of the spectral fitting that is trespassory to obtain energy flux measurements. Although BATSE’s spectroscopic detectors (SDs) can gather enough data to produce spectra on a 0.128 s timescale, the data usually has to be integrated to get an acceptable fit to a given spectra model. This reduces the resolution of the energy flux time history to a level that is unacceptable for a four or five parameter fit. Furthermore, Borgonovo & Ryde (2003) have shown that, empirically, there is a power-law relation between the count flux and the energy flux, with an index clustering around $`1`$, which supports the use of the count light-curves for our study of the pulse shapes. The background is, in general, fitted to a quadratic polynomial using 1.024 s resolution background data that are available from 10 minutes before the trigger to several minutes after the burst. The data along with the background fit coefficients were obtained from the CGRO Science Support Center (GROSSC) at Goddard Space Flight Center through its public archives.

We selected bursts from the entire BATSE catalog with the criteria that the peak flux be greater than 1.0 photons cm-2 s-1 on a 256 ms timescale. Of these, we selected bursts that exhibited clean, single-peaked events, or in the case of multi-peaked bursts, pulses that were well distinguished and separable from each other. We limited the bursts to events with durations longer than 2 s (full width half max - FWHM), primarily because this study focuses on the properties of long GRB events. There are strong indications that short bursts with $`T_{90}<2.6`$ s have different temporal behaviors compared to long bursts and that they may actually constitute a different class of GRBs (Norris, Scargle, & Bonnell, 2001). Therefore, in this study we are primarily interested in the shape and temporal properties of long GRB events. In the search we did not assume any functional description of the pulse profile, i.e. pulses that did not look like FREDs were also selected, in order to eliminate any preconceived idea of what the fundamental pulse shape should be. We note that we do not know what the clean, ”generic” pulse actually looks like, and some noise, or substructure, in the comoving light curve can be expected, for instance, if the surface brightness of the shell is non-homogeneous and/or if the density profile of the shell, through which the shock is traversing, has a complicated radial profile. Therefore, we adopt the following criteria. If a significant substructure on top of a smoothly fitted pulse is within a time interval of 1/4 the FWHM of the main pulse peak and it contains more than 15 $`\%`$ of the total flux in the main pulse, the pulse structure is rejected. If the substructure occurs in a time interval of more than twice the FWHM of the primary pulse, then the pulse is treated as an independent event and is fitted separately. We believe that this limits the possible contamination of the sample with pulses that are actually a composite of overlapping emission episodes. The final sample consisted of 76 pulses within 68 bursts. This sample differs from those examined in previous studies such as Norris et al. (1996) because they primarily represent separable pulses that exhibit long decay tails that can be followed to background. Therefore, no attempt was made to perform deconvolution fits of overlapping pulses, for the purpose of obtaining robust measurements of the asymptotic decay of each pulse, which is not possible for overlapping regions of activity. For more details on the sample selection, see F. Ryde $`\&`$ D. Kocevski (in prep.). These bursts are presented in Table 1, and are denoted by both their GRB names and BATSE trigger numbers. The pulses within multi-peaked bursts are distinguished by their $`t_{max}`$ value.

All of the background-subtracted light-curves were fitted with the power-law rise- (eq. 22) and the exponential rise- (eq. 28) functions via $`\chi ^2`$ minimization. To aid in this, we developed a graphical IDL fitting routine which allowed the user to set the initial pulse attributes manually before allowing the fitting routine to converge on the best fit model. The convergence routine utilized the Marquardt technique to minimize the difference between the model and the time profile by adjusting the model parameters and tracking the resulting effects on the $`\chi ^2`$. In the case of the power-law rise-function, the free parameters were: peak flux $`f_{max}`$, position of peak flux $`t_{max}`$, the rise index $`r`$, and the decay index $`d`$. Most of the GRBs had significant emission prior to the BATSE trigger time, so we introduced a fifth parameter $`t_0`$ which measures the offset between the start of the pulse and the trigger time. The initial value of $`t_0`$ was typically set to where the light curve would rise above 10$`\%`$ of the background, helping to eliminate the ambiguity associated with the start of the pulse. In the case of the exponential rise-function, the free parameters were: pulse amplitude $`A_0`$, rise time constant $`\tau _r`$, decay time constant $`\tau _d`$, the asymptotic decay power law index d, and the trigger offset $`t_0`$.

## 4 Results

Presented in Figure 4 are four type-fits to the power-law rise-function (eq. 22). These pulses represent a good example of the clean separable FRED pulses chosen for this analysis and the degree of acceptability of $`\chi ^2`$ for our fits. Overall the power-law rise-function appears to give a better fit to most of the pulses compared to the exponential rise-function (eq. 28). This can be seen in the spread of the respective $`\chi ^2`$ distribution for our entire sample shown in Figure 5, with the power-law rise-model giving a mean $`\chi ^2`$ of 1.16 with a standard deviation of 0.19. This is most likely due to the fact that the power-law function allowed for the modelling of both concave and convex rise profiles, whereas the nature of the exponential function is limited to modelling concave shapes. This is a major constraint for the latter model since many GRB pulses that are arguably single emission episodes are seen to have concave rise profiles, at least in the initial stage.

The narrow distribution of $`\chi ^2`$ values also indicates that a power law is sufficient to model the decay portion of the pulse light curve. According to Fenimore et al. (1996), this would not be the case if evolution of the low energy spectral index $`\alpha `$ greatly affected the resulting light curve. As discussed in $`\mathrm{\S }1`$, they show that the decay phase should scale roughly as the power of the low energy spectral index $`T^{\alpha 1}`$. If $`\alpha `$ exhibited a large range of evolution during the decay phase, then a single power law would no longer be adequate to model the temporal behavior of the light curve. The $`\chi ^2`$ distributions support the assumption that the evolution of the low energy spectral index does not significantly alter our results.

### 4.1 Rise and Decay Timescales

The primary pulse attributes that can be derived from fitting the time profiles are related to the pulse rise and decay rates. In the case of the power-law rise-function, these consist of the asymptotic rise and decay power-law indices $`r`$ and $`d`$ described in $`\mathrm{\S }`$2.3. The resulting $`r`$ and $`d`$ parameter distributions for our entire sample are shown in Figure 6. The rise index appears to be the tighter of the two, with a mean value of $`r=1.49\pm 0.25`$ and a standard deviation of 0.63. The distribution of the decay index is somewhat broader with a mean value of $`d=2.39\pm 0.12`$ with a standard deviation of 0.76. The exponential-rise function has no corresponding rise index, but rather the rise timescale is associated with the time constant $`\tau _r`$ in equation (28). The fit results show that the rise timescale covers a very broad range of values with a median value of $`\tau _\mathrm{r}=2.40\pm 0.10`$ and a standard deviation of 2.50. The asymptotic decay index is similar to that found with the power law rise model except that it has a much broader distribution with a median value of $`d=2.75\pm 0.12`$ with a standard deviation of 1.50.

### 4.2 Pulse Asymmetry

A measure of the average pulse time-asymmetry can be easily constructed using the fit parameters obtained from both models. Several authors have studied burst asymmetry properties in the past (Norris et al., 1996; Fenimore et al., 1996) \- and found interesting correlations between the pulse asymmetries and other temporal properties. Norris et al. (1996) measured asymmetry values for some 400 fitted pulses and reported a general trend that narrower pulses tend to be more symmetric and have harder spectra. This is often referred to as the GRB pulse paradigm. Here we define the asymmetry as the fraction of the FWHM belonging to the rise phase of the pulse to the fraction belonging to the decay phase (see Figure 7 for an illustration). This description is independent of the model parameters used to fit the pulses and hence allows the asymmetry to be measured for both of the functions utilized in this analysis in a comparative manner. The resulting asymmetry distribution for the power-law rise-model is remarkably tight with a median value of 0.47 $`\pm `$ 0.08 and a standard deviation of 0.09. The exponential rise model provides similar results yielding a median value of 0.45 $`\pm `$ 0.07 and a standard deviation of 0.07.

Alternatively, the asymmetry can be defined using the $`T_{90}`$ of the pulse instead of the FWHM. The rise phase would then be defined as the interval from when where 5$`\%`$ of the pulse counts are accumulated till the time of peak flux as measured by the model parameter $`t_{max}`$. Similarly, the decay phase is simply the time between the peak flux till when 95$`\%`$ of the pulse counts are collected. Note that here we are using $`T_{90}`$ for the individual $`pulses`$, not the entire burst which is the common practice. Although this description is straightforward, it is somewhat more problematic when dealing with multi-pulsed bursts because subsequent pulses will routinely begin before the signal from the first pulse falls to within 5$`\%`$ of the background. It is presumed that the decay phase of the pulse is hidden under the rise phase of the following pulse, therefore in these cases the decay time is only a lower limit to the true length of the decay phase of the pulse. This occurred in 10$`\%`$ of the pulses observed, so the effect on the overall distribution should be minimal. The resulting median asymmetry values using the $`T_{90}`$ definition for both the power law and exponential models are 0.30 $`\pm `$ 0.04 and 0.21 $`\pm `$ 0.03 respectively. The resulting asymmetry distributions for the power-law rise-function for the FWHM definition is plotted in Figure 8.

Surprisingly, we find no evidence of a connection between the pulse asymmetry and the pulse duration ($`T_{90}`$) for the individual FRED pulses. Likewise, we find no significant correlations between the pulse asymmetry and the pulse width (FWHM). The lack of any significant correlation can be seen in Figure 9 and measured quantitatively by the use of the Spearman rank-order statistic which provides a robust and convenient means of evaluating the statistical significance of a given correlation. A Spearman rank order analysis on the data shown in Figure 9 yields a correlation coefficient of $`R=0.057`$, where $`R=1`$ represents a 1:1 correspondence. Both of these conclusions are contrary to the pulse paradigm referred to earlier and the disagreement is most likely due to the difference in samples used in the two studies. The Norris et al. (1996) study examined complex overlapping bursts for a majority of their data sample, while the majority of our sample consist of single-peaked FRED bursts. This leads to the conclusion that most of the simple FRED pulses that were examined in our sample are extremely self-similar, independent of their duration and amplitude.

### 4.3 Correlations

We now examine several correlations between temporal properties measured from the pulse profiles such as the FWHM, rise and decay times and fluence values. First, Figure 10 shows the pulse rise time versus the pulse width, as measured by the FWHM, for our entire sample. The resulting linear correlation has a slope of 0.323 $`\pm `$ 0.01 with an associated Spearman rank order correlation coefficient of $`R=0.981`$ and is a direct consequence of the tight asymmetry distribution discussed above. This same effect is also observed between the time constant $`\tau _r`$ derived from the exponential rise model and the pulse width. Here again, the rise time constant, which is an alternative way of measuring the characteristic rise timescale, is linearly correlated to the FWHM. Note that the FWHM values plotted in Figure 10 are not normalized to any timescale such as the pulse duration, reasserting that the rise time is a constant fraction of the overall pulse width independent of the duration of the pulse.

Another correlation born from the data is the relationship between the normalized rise time and the decay power-law index, which is shown in Figure 11 and has a correlation coefficient of $`R=0.799`$. Here the rise time is normalized to the FWHM of the pulse, so the longer the rise time is, relative to the FWHM, the steeper the resulting power-law decay of the pulse envelope becomes. This is very interesting since the rise time is thought to be proportional to the dynamic or shell-crossing time of the pulse and is generally believed to be independent of the subsequent fall of the FRED profile. A similar effect can be seen in Figure 12 where the pulse asymmetry data are plotted versus the asymptotic decay power-law index. The degree of asymmetry is seen to be proportional ($`R=0.831`$) to the decay index in such a way that pulses with steep decay profiles (large $`d`$) have a larger rise-to-decay ratio. The latter correlation is to be expected since, as $`d`$ tends to 0, the decay time approaches $`\mathrm{}`$ and hence the asymmetry tends to 0. A similar correlation is not seen with the rise power-law index.

## 5 Discussion

The fit results presented above indicate that the analytic functions derived in §3 are very successful in describing the overall profile of individual GRB pulses. These descriptions differ from those used in previous modelling surveys as they are motivated by physical first principles within the general fireball model with the requirement that the previously observed empirical relations for the spectral evolution must hold. These descriptions also build on the original Ryde & Svensson (2002) analysis by including a description of the rise phase of the pulse. Of the two assumed rise profiles presented in §2, a simple, power-law rise seems to better describe the majority of the FRED pulses examined. This is one of the simplest assumptions that can be made regarding the dynamic (or shell crossing) phase of a GRB pulse and may point to the explicit time-dependance of the parameters internal to the shock such as Thompson depth or the magnetic field. Hence, the measured pulse parameters have the potential to be used to diagnose pulse characteristics, such as bulk Lorentz factors, $`\mathrm{\Gamma }`$, shell radii, and thickness.

The decay power-law index appears to cover a wide range of values with a median centered around 2.39 $`\pm `$ 0.12, which is 3.25 $`\sigma `$ steeper than the analytically derived value of $`d=2`$ for pulse profiles created principally by relativistic effects due to spherical symmetric shells. Overall, 30 of the 77 pulses that we analyzed had decay indices that were consistent (within 1$`\sigma `$) with the predicted $`d=2`$ value, constituting 39$`\%`$ of our sample. We note that the distribution of HIC indices, which were shown in $`\mathrm{\S }2`$ to be directly related to the decay index, found by Borgonovo & Ryde (2001) also has a broad distribution, with several cases differing substantially from $`\eta =2`$. There are several mechanisms that can cause the decay envelope to deviate from the predicted spherical scenario, but the most obvious would be the breaking of local spherical symmetry. Both prolate and oblate shell geometries as seen by the observer would create light curves with resulting power-law indices that differ from the spherical case. In $`\mathrm{\S }1`$ we derived the shape that the light curve is expected to take if the pulse profile were created purely by a spherical shell. In the case of an elliptical geometry, the time delay between the arrival of off axis photons will be given by $`\mathrm{\Delta }t=(R_0R\mu )/c`$, where $`R_0`$ is the radius of the shell at $`\theta =0`$. The radius $`R`$ of the shell is no longer independent of the angle $`\theta `$, and is given by

$$R=a\sqrt{\frac{1e^2}{1e^2\mu ^2}}$$
(30)
where e is the eccentricity of the shell, a is the semi-major axis, and $`\mu =cos(\theta )`$. Substituting this back into $`\mathrm{\Delta }t`$, we find that

$$\mathrm{\Delta }t=aa\sqrt{\frac{1e^2}{1e^2\mu ^2}\mu }$$
(31)
Solving for $`\mu `$ and substituting into the Doppler boost expression (equation 3) gives the Doppler profile (Doppler boost as a function of time) for elliptical shells

$$𝒟(\mathrm{\Gamma },t)=\frac{1}{\mathrm{\Gamma }(1\frac{\beta (t+a)}{\sqrt{e^2t^22e^2taa^2}})}$$
(32)
Following from the discussion in $`\mathrm{\S }1`$, the resulting evolution of the spectral break-energy and light curve of the energy-flux should follow $`E_{\mathrm{pk}}(t)=E_{\mathrm{pk}}^{}𝒟`$ and $`F_{bol}(t)=F_0𝒟^2`$, so the profile of the varying Doppler boost should directly give rise to the time dependence of the $`E_{pk}`$ evolution and the pulse light curve. The Doppler profile, expressed in equation 32, is shown for several different prolate and oblate shock fronts in Figure 13, with the corresponding shell geometry plotted as an inset. It shows that the prolate shells (solid lines) have steeper Doppler profiles resulting in asymptotic power law slopes which are greater than the spherical case (thick solid line) whereas the opposite is true for the oblate geometries (dashed lines). The extreme limit of the oblate scenario is such that the shell becomes a parallel slab, in which case for an unresolved source the Doppler profile approaches zero and any pulse evolution is directly due to the shell’s intrinsic emission profile. Because of this effect, the observed range of power-law decay indices may be explained by a simple distribution of shell geometries, where the median decay index of 2.39 $`\pm `$ 0.12 would indicate that, in the context of relativistic curvature, the majority of the analyzed pulses were produced by shells with a degree of curvature greater than that exhibited by a spherical shell. Such conditions are not unreasonable in the context of the fireball model and a simple angular dependance of the bulk Lorentz factor can easily produce an elliptical shock front that evolves as the shock propagates. The evolution of the shock geometry, from oblate at early times to prolate at later times, and a variation in the time at which the shell ”ignites” during this evolution, could account for the distribution of decay indices.

Fenimore et al. (1996) and Fenimore & Sumner (1997) provide several other arguments in support of the breaking of local spherical symmetry, the most prominent being that the pulse photon light curves appear to evolve uncoupled to the GRB spectra. They find that the GRB spectra also evolve faster than the predicted rate of evolution due to the simple boost factor that is expected from a spherically symmetric shell. According to the investigation performed by Ryde & Petrosian (2002), this result is consistent with our findings that the average light curve profile also decays faster than predicted. In the context of their model, a range of relative sizes between the comoving and the curvature timescales defined in $`\mathrm{\S }1`$ can result in a variety of pulse shapes. Light curves that exhibit a substantially different shape than that predicted by spherical curvature are attributed to a scenario where the comoving timescale dominates within a burst that already has an intrinsically fast decay rate ($`d>2`$) of the comoving light curve. This, according to their model, must also be manifested in the spectral evolution in a corresponding manner. Therefore, an analysis of the spectral evolution rates and $`\alpha `$ values of the pulses that vary the furthest from the predicted scenario of $`d=2`$ may test how much these values affect the resulting light curve profiles. This research is currently being pursued (F. Ryde,& D. Kocevski in prep).

The asymmetry results are rather surprising and the reason for the tight asymmetry distribution among such a variety of GRB pulses is not immediately clear. It would indicate that the dynamical and angular timescales are not fully independent and that the observed rate at which the shell becomes active is dependant on the decay timescale and hence the curvature of the shell. A situation can be envisioned where a short comoving rise profile is delayed and boosted to the same degree as the shell’s intrinsic cooling profile. If both of these timescales are extremely short as compared to the angular spreading timescale, then a relationship can be produced between the observed rise and decay profile. This effect could produce a narrow asymmetry distribution and would also predict a correlation between the rise timescale and decay power-law indices, which is indeed observed in our sample. Ryde & Petrosian (2002) have discussed extensively the effect of varying intrinsic emission profiles and angular spreading timescales on the resulting light curve that is observed. They find that the resulting pulse asymmetry relies heavily on the convolution of the intrinsic emission profile and the angular spreading timescale. Therefore, future modelling of pulse timescales and shell geometries, with the constraint that they produce the observed asymmetry distribution, may yield clues to the intrinsic profile of the emission episode.

We are grateful to the GROSSC at NASA/GSFC for providing the HEASARC Online Service. We should also like to thank Dr. Vahé Petrosian for his valuable discussions and insights. D.K. wishes to express his gratitude to the Department of Physics at Stanford University for hospitality and also acknowledges the NASA GSRP fellowship program for their support. F.R. acknowledges financial support from the Swedish Foundation for International Cooperation in Research and Higher Education (STINT) and the Ludovisi Boncomagni, née Bildt, foundation.


text:
# DFT calculation of the intermolecular exchange interaction in the magnetic Mn4 dimer

## Abstract

The dimeric form of the single-molecule magnet \[Mn4O3Cl4(O2CEt)3(py)3\]2 recently revealed interesting phenomena: no quantum tunneling at zero field and tunneling before magnetic field reversal. This is attributed to substantial antiferromagnetic exchange interaction between different monomers. The intermolecular exchange interaction, electronic structure and magnetic properties of this molecular magnet are calculated using density-functional theory within generalized-gradient approximation. Calculations are in good agreement with experiment.

Single-molecule magnets (SMMs), such as \[Mn12O12(CH3COO)16(H2O)4\]$``$2(CH3COOH)$``$4(H2O) (hereafter Mn12)LIS80  and \[Fe8O2(OH)12(tacn)6\]Br8$``$9(H2O) (hereafter Fe8)WIEG84  have received tremendous attention due to macroscopic quantum tunnelingVILL94  and possible use as nanomagnetic storage devices. Hysteresis loop measurements on the SMMs Mn12 and Fe8 showed magnetization steps at low temperatures upon magnetic field reversal.SESS93  This is due to quantum tunneling between spin-up states and spin-down states despite a large effective spin $`S`$$`=`$10 for each molecule. The resonant tunneling fields in these systems are primarily determined by the magnetomolecular anisotropy. Recently a dimerized single-molecule magnet \[Mn4O3Cl4(O2CEt)3(py)3\]2 (hereafter Mn4 dimer) where Et=CH2CH3 and py=NC5H5, has been formedHEND92 ; WERN02-NAT  which exhibited qualitatively different tunneling behavior: quantum tunneling prior to magnetic field reversal and an absence of quantum tunneling at zero field in contrast to other SMMs such as Mn12 and Fe8.WERN02-NAT  To understand the basis for the qualitative deviation we have calculated both the magnetomolecular anisotropy and the intermolecular exchange interaction in the Mn4 dimer using density-functional theory. Our results confirm that there exists an appreciable antiferromagnetic exchange interaction between monomers and that tunneling fields in this dimer are strongly influenced by the presence of the monomer-monomer exchange interaction. This interaction produces a bias field that encourages monomeric magnetic-moment reversal below zero field and prevents two monomers from simultaneously flipping their magnetic moments at zero field. We determine that the origin of the exchange interaction is not dominated by either kinetic or exchange-correlation terms and that the total ”exchange” interaction is in fact an order of magnitude smaller than the kinetic contribution. For Mn12 and Fe8, the intermolecular exchange interaction has not been observed experimentally and it is generally accepted that the overlap between neighboring molecules is negligible.

In this work, we discuss calculations on the Mn4 dimer which is formed by inversion of the three-fold symmetric monomer shown in Fig.1. The magnetic core of the Mn4 monomer consists of three ferromagnetically coupled Mn3+ ($`S`$$`=`$2) ions coupled antiferromagnetically to the remaining Mn4+ ($`S`$$`=`$3/2) ion leading to a total ground-state spin of $`S`$$`=`$2$`\times `$3$``$3/2$`\times `$1$`=`$9/2 (refer to Fig.1). The core has a similar cubane structure as the inner core of the SMM Mn12, although there are four Mn4+ for Mn12. We investigate the electronic structure and magnetic properties of this SMM Mn4 using density-functional theory (DFT). We calculate optimized geometries for the Mn4 monomer and dimer, their binding energy, the monomeric magnetic anisotropy barrier (MAE), and the exchange coupling constant between monomers. Results are compared with experiment.

Our DFT calculationsKOHN65  are performed with the all-electron Gaussian-orbital-based Naval Research Laboratory Molecular Orbital Library (NRLMOL).PEDE90  Here we use the Perdew-Burke-Ernzerhof (PBE) generalized-gradient approximation (GGA).PERD96  Before discussing energetics and magnetic phenomena we discuss two structural issues. First, we have considered monomers and dimers that are terminated by both H and by the CH2CH3 radicals found in the experimental structure. Second we have considered structures based on two conformers of the monomeric unit. While a complete vibrational analysis will be discussed in a later publication, all indications are that both conformers are stable. The conformers have slightly different arrangements of the pyridine ligands. The first conformer was identified by a density-functional-based geometry optimization of the hydrogenated monomer. The second conformer was identified by improvements on the monomer deduced from the experimental x-ray data. In the remainder of the paper we refer to these monomers as the computationally determined conformer (CDC) and the experimentally determined conformer (EDC).

Each Mn4 monomer has threefold symmetry so there are 26 inequivalent atoms to consider. The number of inequivalent atoms is reduced to 20 when the CH2CH3 radical is replaced by H. A pyridine ring is initially constructed to lie in the plane defined by the vector connecting Mn3+ and neighboring N and the sum of the two vectors connecting Mn3+ with two closest Cl’s (refer to Fig. 1). The geometries for the pyridine ring and the cubane were first optimized separately to generate an initial geometry for the DFT calculations on the full monomer. The initial geometry for the monomer was relaxed using NRLMOL with the Cl atom fixed to reproduce the experimental Cl-Cl distance (3.86 Å) upon dimerization (i.e. adding inversion symmetry). Relaxation continues until forces exerted on all atoms become $``$0.001 hartree/bohr. The CDC dimer is then obtained by inversion of the CDC monomer with the fixed value of $`d=3.86`$ Å (marked as dotted in Fig. 1). For the case of the x-ray deduced experimental geometry, the C-H bond lengths are underestimated (0.71 Å to 0.96 Å) in comparison to standard hydrogen bond lengths, which yields self-consistent forces on hydrogen atoms as large as 0.8 hartree/bohr. To improve the experimental geometry, all hydrogen positions were first moved to create C-H bond lengths as 1.1 Å, and then additional optimization of the experimental geometry was performed with the fixed Cl-Cl distance. The experimental geometry without corrected hydrogen positions was 53 eV higher in energy than that of the structure with corrected hydrogen positions. Hereafter unless we specify, the EDC monomer refers to the optimized experimental geometry with corrected hydrogen positions.

We have used full basis sets for all six different atoms and fine mesh.PORE99  Charges and magnetic moments for Mn’s from the CDC monomer agree well with those from the EDC monomer. For example, a sphere with a radius of 2.23 Bohr captures charges of 23.4 and 23.7, and magnetic moments of $`3.6\mu _B`$ and $`2.5\mu _B`$ for Mn3+ and Mn4+ respectively. The total magnetic moment for the monomer is $`9\mu _B`$ in good agreement with experiment. The HOMO-LUMO gap for majority (minority) spin is 1.02 eV (2.42 eV). The energy difference between the minority (majority) LUMO and the majority (minority) HOMO is 1.17 eV (2.28 eV), which ensures that the system is stable with respect to the total magnetic moment. As clearly seen in Fig. 2, right below the Fermi level for majority spins the projected Mn3+(3$`d`$) DOS is dominant over the projected Mn4+(3$`d`$) DOS, while for minority spin the opposite trend is observed. This confirms the experimental picture of three Mn3+ spins antiferromagnetically coupled to a Mn4+ spin.

We calculate the binding energy by subtracting the dimer energy from twice the monomer energy. We find that the dimer is stable for both the CDC and EDC. For the CDC (EDC), the binding energy is about 0.16 eV (0.78 eV). The magnitude of the binding energy suggests attractive electrostatic interactions between different monomers. The discrepancy between the binding energy for the CDC and that for the EDC may be attributed to our substitution of ethyl for hydrogen in the CDC and/or the fact that the plane where a pyridine ring sits is different for both geometries. To check the former possibility, we calculate the binding energy of the EDC terminated by hydrogen, and obtain 0.45 eV. We have also verified that the conformation of a pyridine ring for the EDC is slightly different from that for the CDC. Thus, the discrepancy arises from both reasons.

We have calculated the monomeric MAE in zero magnetic field for both the CDC and EDC with the assumption that spin-orbit coupling is a major contribution to the MAE. For this calculation, we follow the procedure developed in Ref.PEDE99, . Our calculations show that the Mn4 monomer has uniaxial anisotropy along the threefold axis (the bond between Mn4+ and Cl in the cubane), in agreement with experiment.HEND92 ; WERN02-NAT  For uniaxial systems, the energy shift $`\mathrm{\Delta }`$ due to the spin-orbit interaction can be simplified to $`\gamma _{zz}S_z^2`$ up to constant terms independent of $`S_z`$ if the $`z`$ axis is assigned as the easy axis.PEDE99  Then the classical barrier (MAE) to be overcome to monomer magnetization reversal $`M_z`$$`=`$+9/2 to $`M_z`$$`=`$$``$9/2 is $`\gamma _{zz}((9/2)^2(1/2)^2)`$. For the CDC (EDC) monomer, the MAE is 11.3 K (11.6 K), which is close to that for the hydrogenated EDC monomer. As shown in Table 1, all these numbers are close to the experimental value of 14.4 K. The difference between our estimated MAE and the experimental value might be ascribed to other effects on the barrier such as spin-vibron coupling.PEDE02

To calculate the exchange coupling constant $`J`$ between monomers, we assume that a monomer is an ideal $`S=9/2`$ object and that its effective spin is aligned along the easy axis and of Ising type (either $`M_z`$$`=`$+9/2 or $`9/2`$). Then we calculate self-consistently energies of ferromagnetic (parallel monomeric spins) and antiferromagnetic configuration (antiparallel monomeric spins) of the dimer, and take a difference $`\delta `$ between the two energies. We find that the antiferromagnetic configuration is favored. The antiferromagnetic exchange constant $`J`$ is determined from $`\delta =2J(9/2)^2`$. For the CDC, the energy difference is 31 microhartree so that $`J`$$`=`$0.24 K, while for the EDC, $`J`$$`=`$0.27 K. These can be compared to the experimentally measured value of $`J`$$`=`$0.1 K.WERN02-NAT  The numerical uncertainty in the total-energy difference for our DFT calculations is at most 5 microhartree, which can be translated to the uncertainty in the exchange $`J`$ as 0.04 K. We achieve high-accuracy in the total-energy difference, because we use exactly the same optimized dimer geometry with the same parameter values for a self-consistent approximation except for the effective spin configurations of monomers. Although our DFT estimated value of $`J`$ is somewhat higher than the experimental value, this may be acceptable considering the assumptions we made and the fact that DFT calculations often overestimate exchange interactions. In some cases, the PBE generalized-gradient approximation may not fully cancel the self-interaction in the Coulomb potential. Therefore, the electrons in our calculations are slightly more diffuse, which should lead to overestimated exchange interaction.

It is interesting to examine whether the exchange interaction varies significantly with the monomer-monomer separation. We consider the case that each monomer is displaced toward or away from the center of mass of the dimer along the easy axis. Then we calculate the exchange constant $`J`$ for the EDC dimer with three different monomer-monomer distances from the experimentally measured value. The monomer-monomer distance is varied by changing the two central Cl-Cl distance with a monomer geometry fixed. If the central Cl-Cl bond length increases by 1 Bohr, then $`J`$ decreases down to 0.10 K. If the bond length decreases by 0.5 Bohr (1 Bohr), $`J`$ increases to 0.47 K (0.81 K). Table 1 summarizes the separation dependence of $`J`$ and of the monomeric MAE. As shown in Table 1, the monomeric MAE does not depend on the exchange interaction between monomers, because the monomer geometry has not changed during this process. Figure 3 shows that $`J`$ increases exponentially with decreasing the separation distance. This tells us how quickly the overlaps of neighboring wavefunctions decrease with increasing the distance. We have decomposed the J values into kinetic, coulombic and exchange-correlation contributions. The kinetic contribution is an order of magnitude larger than the total value of J and it is significantly cancelled by the exchange-correlation contributions to the J value.

Since we estimated the anisotropy barrier and exchange constant, we can construct a model Hamiltonian for the dimer according to

$``$ $`=`$ $`\gamma _{zz}(S_{1z}^2+S_{2z}^2)+J\stackrel{}{S}_1\stackrel{}{S}_2`$ (1)
where the uniaxial anisotropy parameter $`\gamma _{zz}=0.58`$ K and $`J=0.27`$ K. To determine whether our values of $`\gamma _{zz}`$ and $`J`$ can reproduce the experimental values of the resonant tunneling fields (Fig. 4 in Ref.WERN02-NAT, ), we calculate these fields using exact diagonalization of Hamiltonian (1). Although it is crucial to include some small transverse terms in the Hamiltonian (such as transverse anisotropy and transverse fields) for calculations of tunnel splittings, the transverse terms do not affect the resonant fields much. Table 2 summarizes the resonant fields for some low-energy states for three different values of $`\gamma _{zz}`$ and $`J`$: (1) the experimental values $`\gamma _{zz}=0.72`$ K, $`J=0.1`$ K; (2) $`\gamma _{zz}=0.58`$ K, $`J=0.27`$ K; (3) slightly modified version of our DFT results $`\gamma _{zz}=0.58`$ K, $`J=0.1`$ K. Let us focus on two tunnelings which were prominent in the experimental measurements: $`|M_1=9/2,M_2=9/2|M_1=9/2,M_2=9/2`$, and $`|9/2,9/2|9/2,9/2`$, where $`M_1`$ and $`M_2`$ are the eigenvalues of the spin operator projected along the easy axis for each monomer. For these two tunnelings, the resonant fields are solely determined by $`J`$ and are independent of $`\gamma _{zz}`$: $`B_{\mathrm{res}}9J/(2g\mu _B)`$. Therefore, model Hamiltonian (1) with our estimated values will not quantitatively reproduce the experimental resonant fields. However, in this case (when $`\gamma _{zz}`$ becomes comparable to $`J`$), we notice that the hysteresis loop exhibits richer features such as more magnetization steps before magnetic field reversal. Since DFT often overestimates exchange interactions, we also calculate the resonant fields with $`J`$ decreased to 0.1 K and $`\gamma _{zz}`$ fixed to examine if agreement with experiment improves. We find that some resonances agree with experiment and some do not agree.

In summary, we have calculated optimized geometries for a monomer and dimer of the SMM Mn4 using DFT. For both the CDC and EDC, we calculated binding energy, monomeric MAE, and the exchange interaction between monomers. The binding interaction between monomers is electrostatic. Our calculated anisotropy barrier is close to the experimental value. The exchange interaction between monomers is twice or three times larger than the experimental value. Overall, our DFT calculations are in qualitative accord with experiment.

KP and SLR were funded by W. M. Keck Foundation, MRP was supported in part by ONR and the DoD HPC CHSSI program, and NA and GC were supported by NSF Grants Nos. CHE-0123603 and DMR-0103290.


text:
# References

Quantum Interference Effects in the Detection Probability of Charged Leptons Produced in Charged Current Weak Interactions

J.H.Field

Département de Physique Nucléaire et Corpusculaire, Université de Genève

24, quai Ernest-Ansermet CH-1211Genève 4.

E-mail: john.field@cern.ch

## Abstract

The phases of quantum interference effects in charged lepton production by neutrinos (neutrino oscillations) following pion decay at rest and in flight, muon decay and nuclear $`\beta `$decay at rest as well as for ‘muon oscillations’ following pion decay at rest and in flight are calculated. The same phase is found for neutrino and muon oscillations following pion decay at rest. The results found disagree with the conventionally used value: $`\varphi _{12}=\mathrm{\Delta }m^2L/(2P)`$. Differences to previous treatments are briefly commented on.

PACS 03.65.Bz, 14.60.Pq, 14.60.Lm, 13.20.Cz
Keywords ; Quantum Mechanics, Neutrino Oscillations.

This letter sketches briefly a derivation of spatially dependent interference effects (conventionally called ‘neutrino oscillations’) in the detection probability of charged leptons produced in neutrino interactions following the weak decays of unstable ‘source’ particles: pions, muons and $`\beta `$radioactive nuclei. Similar effects for the detection probability of the decay products of muons produced in pion decay (‘muon oscillations’) are also considered. The calculations are based on Feynman’s path amplitude formulation of quantum mechanics  in which the probability of transition from a set of initial states $`I=_li_l`$ to a set of final states $`F=_mf_m`$ is given by the relation:

$$P_{FI}=\underset{m}{}\underset{l}{}\left|\underset{k_1}{}\underset{k_2}{}\mathrm{}\underset{k_n}{}f_m|k_1k_1|k_2\mathrm{}k_n|i_l\right|^2$$
(1)
where $`k_j,j=1,n`$ are (unobserved) intermediate quantum states. Only the essential elements of the calculation and the results are presented here. Full details are given elsewhere .

The application of Eqn(1) to pion decay at rest is illustrated in the ideal experiment shown in Fig.1. A $`\pi ^+`$ comes to rest in a stopping target T at time $`t_0`$ as recorded by the counter CA (Fig.1a)). The pion at rest constitutes the initial state of the path amplitudes. In Fig.1b) and Fig.1c) are shown two alternative histories of the stopped pion. In Fig.1b) the pion decays at time $`t_1`$ into the neutrino mass eigenstate $`|\nu _1>`$, of mass $`m_1`$, and in Fig.1c) into the neutrino mass eigenstate $`|\nu _2>`$, of mass $`m_2`$, at the later time $`t_2`$. If $`m_1>m_2`$, then, for a suitable choice of the times $`t_1`$ and $`t_2`$, interference between the path amplitudes corresponding to the different physical processes shown in Fig.1b) and Fig.1c) will occur when a neutrino interaction $`(\nu _1,\nu _2)ne^{}p`$ takes place at time $`t_D`$ and distance $`L`$ from the pion decay point(Fig.1d)). The final state of the neutrino interaction event is also that of the path amplitudes. For this experiment, the path amplitudes corresponding to the two alternative histories of the decaying pion are:

$`A_i`$ $`=`$ $`<e^{}p|T_R|n\nu _i>U_{ei}D(x_fx_i,t_Dt_i,m_i)U_{i\mu }`$ (2)
$`<\nu _i\mu ^+|T_R|\pi ^+>e^{\frac{\mathrm{\Gamma }_\pi }{2}(t_it_0)}D(0,t_it_0,m_\pi )i=1,2`$
Here $`U_{\alpha i}`$ is the Maki-Nakagawa-Sakata (MNS)  matrix element giving the charged current coupling strength of the charged lepton of flavour $`\alpha `$ ($`\alpha =e,\mu ,\tau `$) to the neutrino mass eigenstate $`\nu _i`$. The ‘reduced’ matrix elements $`<e^{}p|T_R|n\nu _i>`$ and $`<\nu _i\mu ^+|T_R|\pi ^+>`$ are given by setting $`U_{\alpha i}=1`$ in the charged current, so that, for example, $`<e^{}p|T|n\nu _i>=U_{ei}<e^{}p|T_R|n\nu _i>`$. Since the purely kinematical effects of non-vanishing neutrino masses are expected to be negligible, then, to a very good approximation, $`|\nu _1>`$ and $`|\nu _2>`$ may be replaced, in the reduced matrix elements, by the massless neutrino wavefunction $`|\nu _0>`$. The former are then independent of neutrino flavour. In the case of a $`2\times 2`$ MNS matrix for the first two generations of leptons, the unitarity of the matrix implies that it is completely defined by a single real parameter, conventionally chosen to be an angle, $`\theta `$, such that:

$`U_{e1}`$ $`=`$ $`U_{1e}=U_{\mu 2}=U_{2\mu }=\mathrm{cos}\theta `$ (3)
$`U_{e2}`$ $`=`$ $`U_{2e}=U_{\mu 1}=U_{1\mu }=\mathrm{sin}\theta `$ (4)
In Eqn(2), $`m_\pi `$ and $`\mathrm{\Gamma }_\pi `$ are the pion pole mass and decay width and $`D`$ is the Lorentz invariant configuration space propagator  of a neutrino or the pion. In the limit of large time-like separations, or of on-shell particles, appropriate to the experiment shown in Fig.1, $`D\mathrm{exp}[im\mathrm{\Delta }\tau ]`$ where $`m`$ is the pole mass of the particle and $`\mathrm{\Delta }\tau `$ the increment of proper time corresponding to the path. In the following, the additional functional dependence $`(m/\mathrm{\Delta }\tau )^{\frac{3}{2}}`$ of $`D`$ in the asymptotic region (leading to solid angle correction factors) is neglected. With these approximations:

$`D(\mathrm{\Delta }x,\mathrm{\Delta }t,m)`$ $``$ $`\mathrm{exp}[im\sqrt{(\mathrm{\Delta }t)^2(\mathrm{\Delta }x)^2}]`$ (5)
$`=`$ $`\mathrm{exp}[im\mathrm{\Delta }\tau ]`$
$``$ $`\mathrm{exp}[i\mathrm{\Delta }\varphi ]`$
The phase increments, $`\mathrm{\Delta }\varphi `$, corresponding to the paths of the neutrinos and the pion in the amplitudes $`A_i`$ are:

$`\mathrm{\Delta }\varphi _i^\nu `$ $`=`$ $`m_i\mathrm{\Delta }\tau _i={\displaystyle \frac{m_i^2}{E_i}}\mathrm{\Delta }t_i={\displaystyle \frac{m_i^2}{P_i}}L`$ (6)
$`\mathrm{\Delta }\varphi _i^\pi `$ $`=`$ $`m_\pi (t_it_0)=m_\pi (t_Dt_0){\displaystyle \frac{m_\pi L}{v_i}}`$ (7)
$``$ $`m_\pi (t_Dt_0)m_\pi L\left\{1+{\displaystyle \frac{m_i^2}{2P_0^2}}\right\}`$
where $`m_i`$, $`E_i`$, $`P_i`$ and $`v_i`$ are the mass, energy, momentum and velocity of the mass eigenstate $`\nu _i`$ and

$$P_0=\frac{m_\pi ^2m_\mu ^2}{2m_\pi }=29.8\mathrm{MeV}$$
(8)
The neutrinos are assumed to follow classical rectilinear trajectories such that $`\mathrm{\Delta }t_i=L/v_i=E_iL/P_i`$ and the time dilatation formula $`\mathrm{\Delta }t=\gamma \mathrm{\Delta }\tau =E\mathrm{\Delta }\tau /m`$ has been used in Eqn(6). In Eqn(7) the neutrino velocityaaaUnits with $`\mathrm{}=c=1`$ are used. $`v_i`$ is expressed in terms of the neutrino mass to order $`m_i^2`$. Using Eqns(5-7), neglecting terms of O($`m^4`$) and higher, and replacing the massive neutrino wavefunctions in the reduced matrix elements by those of massless neutrinos, the path amplitudes of Eqn(2) may be written as:

$`A_i`$ $`=`$ $`<e^{}p|T_R|n\nu _0>U_{ei}U_{i\mu }<\nu _0\mu ^+|T_R|\pi ^+>`$ (9)
$`\mathrm{exp}[i\varphi _0{\displaystyle \frac{\mathrm{\Gamma }_\pi }{2}}(t_Dt_0t_i^{fl})]\mathrm{exp}i\left[{\displaystyle \frac{m_i^2}{P_0}}\left({\displaystyle \frac{m_\pi }{2P_0}}1\right)L\right]i=1,2`$
where the neutrino times-of-flight $`t_i^{fl}=t_Dt_i`$ have been introduced and

$$\varphi _0m_\pi (Lt_D+t_0)$$
(10)
Using now Eqn(1) to calculate the transition probability, and integrating over the detection time $`t_D`$ , gives, for the probability to observe the reaction $`(\nu _1,\nu _2)ne^{}p`$ at distance $`L`$ from the decay point:

$$P(e^{}p|L)=C_N(\nu ;\pi )\mathrm{sin}^2\theta \mathrm{cos}^2\theta (1F^\nu (\mathrm{\Gamma }_\pi )\mathrm{cos}\varphi _{12}^{\nu ,\pi })$$
(11)
where

$$\varphi _{12}^{\nu ,\pi }=\frac{\mathrm{\Delta }m^2}{P_0}\left(\frac{m_\pi }{2P_0}1\right)L=\frac{2m_\pi m_\mu ^2\mathrm{\Delta }m^2L}{(m_\pi ^2m_\mu ^2)^2}$$
(12)
$$F^\nu (\mathrm{\Gamma }_\pi )=\mathrm{exp}\left(\frac{\mathrm{\Gamma }_\pi m_\pi }{2m_\mu ^2}\varphi _{12}^{\nu ,\pi }\right)$$
(13)
and Eqns(3) and (4) have been used. Here $`C_N(\nu ;\pi )`$ is an $`L`$ independent normalisation factor and $`\mathrm{\Delta }m^2m_1^2m_2^2`$. The first and second terms in the second member of Eqn(12) are the contributions of $`\mathrm{\Delta }\varphi _i^\pi `$ and $`\mathrm{\Delta }\varphi _i^\nu `$ to the interference phase. The latter is a factor of two larger than in the conventional value  $`\varphi _{12}=\mathrm{\Delta }m^2L/(2P)`$. $`\mathrm{\Delta }\varphi _i^\pi `$ gives a numerically large ($`m_\pi /2P_0=2.34`$) contribution to the oscillation phase.

The above calculation is readily repeated for the case of detection of muon decay $`\mu ^+e^+(\nu _1,\nu _2)(\overline{\nu }_1,\overline{\nu }_2)`$ at distance $`L`$ from the $`\pi ^+`$ decay point. The phase increments analagous to Eqns(6) and (7) are:

$`\mathrm{\Delta }\varphi _i^\mu `$ $`=`$ $`{\displaystyle \frac{m_\mu ^2L}{P_0}}\left[1+{\displaystyle \frac{m_i^2E_0^\mu }{2m_\pi P_0^2}}\right]`$ (14)
$`\mathrm{\Delta }\varphi _i^{\pi (\mu )}`$ $`=`$ $`m_\pi (t_Dt_0){\displaystyle \frac{m_\pi L}{v_0^\mu }}\left[1+{\displaystyle \frac{4m_i^2m_\pi ^2m_\mu ^2}{(m_\pi ^2m_\mu ^2)^2(m_\pi ^2+m_\mu ^2)}}\right]`$ (15)
where

$$E_0^\mu =\frac{m_\pi ^2+m_\mu ^2}{2m_\pi }\mathrm{and}v_0^\mu =\frac{m_\pi ^2m_\mu ^2}{m_\pi ^2+m_\mu ^2}$$

The result found for the time-integrated decay probability is:

$$P(e^+\nu \overline{\nu }|L)=C_N(\mu ;\pi )(1F^\mu (\mathrm{\Gamma }_\pi )\mathrm{sin}2\theta \mathrm{cos}\varphi _{12}^{\mu ,\pi })$$
(16)
where

$$\varphi _{12}^{\mu ,\pi }=\frac{m_\mu ^2\mathrm{\Delta }m^2}{2P_0^3}\left(1\frac{E_0^\mu }{m_\pi }\right)L=\frac{2m_\pi m_\mu ^2\mathrm{\Delta }m^2L}{(m_\pi ^2m_\mu ^2)^2}=\varphi _{12}^{\nu ,\pi }$$
(17)
$$F^\mu (\mathrm{\Gamma }_\pi )=\mathrm{exp}\left(\frac{\mathrm{\Gamma }_\pi m_\pi }{(m_\pi ^2m_\pi ^2)}\varphi _{12}^{\mu ,\pi }\right)$$
(18)
The first and second terms in the second member of Eqn(17) are the contributions of $`\mathrm{\Delta }\varphi _i^{\pi (\mu )}`$ and $`\mathrm{\Delta }\varphi _i^\mu `$ to the interference phase. It is interesting to note that neutrino and muon oscillation phases are the same for given values of $`\mathrm{\Delta }m^2`$ and $`L`$. For oscillation phases $`\varphi _{12}^{\nu ,\pi }=\varphi _{12}^{\mu ,\pi }=1`$, the damping factors of the oscillation term, due to the non-vanishing pion lifetime take the values $`F^\nu (\mathrm{\Gamma }_\pi )=11.58\times 10^{16}`$ and $`F^\mu (\mathrm{\Gamma }_\pi )=14.4\times 10^{16}`$. This damping effect is thus completely negligible in typical neutrino oscillation experiments with oscillation phases of order unity.

Formulae like Eqn(11) have been derived in a similar manner  for ‘$`\overline{\nu }_\mu \overline{\nu }_e`$’ oscillations bbb If the neutrinos are massive, the ‘lepton flavour eigenstates’ $`\overline{\nu }_e`$ and $`\overline{\nu }_\mu `$ are unphysical. However the phrase ‘$`\overline{\nu }_\mu \overline{\nu }_e`$ oscillations’ is still a useful and compact way of describing experiments where anti-neutrinos, produced in association with a muon, give rise to a detection event containing an electron. The amplitudes for all physical processes contain, however, only the physical states $`\nu _i`$. following muon decay at rest, detected via the process $`(\overline{\nu }_1,\overline{\nu }_2)pe^+n`$ and ‘$`\overline{\nu }_e\overline{\nu }_e`$’ oscillations following nuclear $`\beta `$decay, detected via the same process. The results found, for the time-integrated decay probabilities, are:

$`P(e^+n,\mu |L)`$ $`=`$ $`C_N(\overline{\nu };\mu )\mathrm{sin}^2\theta \mathrm{cos}^2\theta \left[1\mathrm{cos}{\displaystyle \frac{\mathrm{\Delta }m^2}{P_{\overline{\nu }}}}\left({\displaystyle \frac{m_\mu }{2P_{\overline{\nu }}}}1\right)L\right]`$ (19)
$`P(e^+n,\beta |L)`$ $`=`$ $`C_N(\overline{\nu };\beta )\left[\mathrm{sin}^4\theta +\mathrm{cos}^4\theta +2\mathrm{sin}^2\theta \mathrm{cos}^2\theta \mathrm{cos}{\displaystyle \frac{\mathrm{\Delta }m^2}{P_{\overline{\nu }}}}\left({\displaystyle \frac{E_\beta }{2P_{\overline{\nu }}}}1\right)L\right]`$ (20)
In these formulae, $`P_{\overline{\nu }}`$ is the momentum of the detected $`\overline{\nu }_1`$, $`\overline{\nu }_2`$ , $`E_\beta `$ is the total energy release in the $`\beta `$decay process, and damping corrections due to the finite lifetimes of the decaying particles are neglected.

Finally, in Ref. , the cases of neutrino and muon oscillations following the decay in flight of ultrarelativistic pions were considered. The oscillation probability formulae are the same as Eqns(11) and (16) respectively, with the phases:

$`\varphi _{12}^{\nu ,\pi }(\mathrm{in}\mathrm{flight})`$ $`=`$ $`{\displaystyle \frac{m_\mu ^2\mathrm{\Delta }m^2L}{(m_\pi ^2m_\mu ^2)E_\nu \mathrm{cos}\theta _\nu }}`$ (21)
$`\varphi _{12}^{\mu ,\pi }(\mathrm{in}\mathrm{flight})`$ $`=`$ $`{\displaystyle \frac{2m_\mu ^2\mathrm{\Delta }m^2(m_\mu ^2E_\pi m_\pi ^2E_\mu )L}{(m_\pi ^2m_\mu ^2)^2E_\mu ^2\mathrm{cos}\theta _\mu }}`$ (22)
The neutrinos or muons are detected at distance $`L`$ from the decay point, along the direction of flight of the parent pion, with decay angles $`\theta _\nu `$ and $`\theta _\mu `$ respectively.

Because the parent pion and the daughter muon in the decay process $`\pi \mu \nu `$ are unstable particles, their physical masses $`W_\pi `$ and $`W_\mu `$ have distributions depending, through Breit-Wigner amplitudes, on their pole masses $`m_\pi `$ and $`m_\mu `$ and decay widths $`\mathrm{\Gamma }_\pi `$ and $`\mathrm{\Gamma }_\mu `$. Energy-momentum conservation in the pion decay process then leads to a distribution of path amplitudes with different neutrino momenta, an effect neglected in the above discussion. In Ref. corrections resulting from the (coherent) momentum smearing due to $`W_\mu `$ and (incoherent) smearing due to $`W_\pi `$ are calculated using a Gaussian approximation for the Breit-Wigner amplitudes. The resulting damping corrections due to $`W_\pi `$ to the interference terms in Eqns(11) and (16) are found to be vanishingly small. For $`\mathrm{\Delta }m^2=`$ (1eV)2 and $`L=30`$m (typical of the LNSD  or KARMEN  experiments) the corresponding damping factors are found to be $`11.3\times 10^{29}`$ (for neutrino oscillations) and $`15.2\times 10^{30}`$ (for muon oscillations) . The damping corrections from the variation of $`W_\mu `$ are even smaller.

Corrections to Eqns(11) and (16) due to thermal motion of the decaying pion and finite target and detector sizes have also been evaluated in Ref.. For neutrino oscillations with $`\mathrm{\Delta }m^2=`$ (1eV)2 and $`L=30`$m and a room-temperature target, the damping factor of the interference term is found to be $`16.7\times 10^{10}`$ and the shift in the oscillation phase to be $`1.2\times 10^9`$rad. In summary, all known sources of damping of the interference terms in Eqn(10) and (15) are expected to be completely negligible in any forseeable neutrino or muon oscillation experiment.

The first published calculation of the neutrino oscillation phase  gave a result (only the contribution of the neutrino propagators was taken into account) in agreement with Eqn(6) above. A later calculation  assumed instead that the phase of the neutrino propagator evolves as $`Et`$, i.e. according to the non-relativistic Schrödinger equation. The assumptions were also made that the different neutrino mass states are produced at the same time and have equal momenta. The first assumption leads to the ‘standard’ formula for the oscillation phase: $`\varphi _{12}=\mathrm{\Delta }m^2L/(2P)`$  which has subsequently been used for the analysis of all neutrino oscillation experiments. The Lorentz-invariant phase of Eqn(6) $`m^2t/E`$ evidently agrees with the result of Ref. in the non-relativistic limit where $`Em`$, but such a limit is clearly inappropriate to describe experiments with ultra-relativistic neutrinos.

The most important difference in the treatment given in the present paper to previous ones that have appeared in the literature is allowing the possibility for the different neutrino mass eigenstates to be produced at different times. Only in this way can the constraints, of both space-time geometry (the detection event is at a unique space-time point), and exact energy-momentum conservation in the decay process , be satisfied.

The other new feature is the inclusion of the important contribution to the oscillation phase from the propagator of the decaying particle, a necessary consequence of the different production times of the different mass eigenstates.

In the standard treatment, the common production time of the mass eigenstates follows from the assumption that a ‘neutrino flavour eigenstate’, that is a coherent superposition of mass eigenstates, is produced in the decay of the source particle. It has been shown that the production of such states in pion decay is incompatible with the measured branching ratio $`\mathrm{\Gamma }(\pi e\nu )/\mathrm{\Gamma }(\pi \mu \nu )`$ and the values of the MNS matrix elements deduced from atmospheric and solar neutrino oscillations . Actually, as first pointed out by Shrock , as a consequence of the Standard Model structure of the charged lepton current:

$$J_\mu (CC)^{lept}=\underset{\alpha ,i}{}\overline{\psi }_\alpha \gamma _\mu (1\gamma _5)U_{\alpha i}\psi _{\nu _i}$$
(23)
the different neutrino mass eigenstates are produced incoherently in different physical processes. This is the basic assumption of the calculations presented above. Clearly, in this case, the different mass eigenstates may be produced at different times in the (classically) alternative histories corresponding to the different path amplitudes. It has also been recently demonstrated  that the factor of two difference in the contribution of neutrino propagation to the oscillation phase between the calculation presented above and the standard formula is a necessary consequence of the assumption of equal production times for all mass eigenstates and hence equal space-time (as contrasted to kinematical cccThe ‘kinematical’ velocity, $`v_{kin}`$ is defined as $`v_{kin}=p/E`$.) velocities. In the same paper  it is also shown that different kinematical hypotheses that lead to different, but unequal, kinematical velocities (equal momenta, equal energies, or exact energy-momentum conservation in the production process) differ in their predictions for the oscillation phase only by terms of O($`m_i^4`$), and so are all equivalent at O($`m_i^2`$).

It may be remarked that the physical interpretation of ‘neutrino oscillations’ provided by the path amplitude description is different from the conventional one in terms of ‘lepton flavour eigenstates’. In the latter, the amplitudes of different lepton flavours in the neutrino are supposed to vary harmonically as a function of time. In the amplitudes for the different physical processes in the path amplitude treatment there is, instead, no variation of lepton flavour in the propagating neutrinos. Only in the detection process itself, the associated MNS matrix elements give different coupling strengths between the different neutrino mass eigenstates and the final state charged lepton and the interference effect occurs between the different path amplitudes that is described as ‘neutrino oscillations’. In the case of the observation of the recoil muons no such ‘lepton flavour projection’ occurs in the detection event, but exactly similar interference effects are predicted to occur. As previously emphasised , the ‘flavour oscillations’ of neutrinos, neutral kaons and b-mesons are just special examples of the universal phenomenon of quantum mechanical superposition that is the physical basis of Eqn(1).

The most important practical conclusion of the work presented here is that identical information on neutrino masses is given by the observation of either neutrinos or muons from pion decay at rest. In view of the possibility of detecting muons simply and with essentially 100 $`\%`$ efficiency, in contrast to the tiny observable event rates of neutrino interactions, the recent indications for ‘$`\overline{\nu }_\mu \overline{\nu }_e`$’ oscillations with $`\mathrm{\Delta }m^2`$(1eV)2 following $`\mu ^+`$ decay at rest  could be easily checked by a search for muon oscillations following $`\pi ^+`$ decay at rest. For $`\mathrm{\Delta }m^2`$(1eV)2, the first absolute maximum of the interference term in Eqn(16) occurs at $`L8`$m. Note that Eqn(16) is valid for any muon detection process that does not distinguish between events where $`\nu _1`$ or $`\nu _2`$ are produced.


text:
# A bird’s eye view of Matrix Distributed Processing

## 1 Introduction

Matrix Distributed Processing (MDP)  is a collection of classes and functions written in C++ to be used as components for fast development of efficient parallel algorithms. Typical algorithms include solvers for partial differential equations, mesh-like algorithms and various types of graph-based problems. These algorithms find frequent application in many sectors of physics, engineering, electronics and computational finance.

MDP components can be divided into two main categories:

* Non parallel components: Linear Algebra components (class mdp\_complex, class mdp\_array, class mdp\_matrix) and Statistical Analysis components (class Measure, class Jackboot)
* Parallel components: (class mdp\_lattice, class mdp\_site, class mdp \_field, etc.)

In this paper we will focus exclusively on the Linear Algebra and the Parallel components111The parallel components can interoperate with other third party C/C++ linear algebra packages and can be used to parallelize existing applications with minimal effort..

MDP is based on MPI and can be used on any machine with an ANSI C++ and support for the MPI communication protocol. No specific communication hardware is required but a fast network switch is suggested. MDP has been tested on Linux PC clusters, SUN workstations and a Cray T3E.

The best way to introduce MDP is to write a program that solves a typical problem:

Problem: Let’s consider the following differential equation:

$$^2\phi (x)=f(x)$$
(1)
where $`\phi (x)`$ is a field of $`2\times 2`$ Complex matrices defined on a 3D space (space), $`x=(x_0,x_1,x_2)`$ limited by $`0x_i<L_i`$, and

$`L`$ $`=`$ $`\{10,10,10\},`$ (2)
$`f(x)`$ $`=`$ $`A\mathrm{sin}(2\pi x_1/L_1),`$
$`A`$ $`=`$ $`\left(\begin{array}{cc}1\hfill & i\hfill \\ 3\hfill & 1\hfill \end{array}\right)`$ (5)
The initial conditions are $`\phi _{initial}(x)=0.`$ We will also assume that $`x_i+L_i=x_i`$ (torus topology).

Solution: In order to solve eq. (1) we first discretize the Laplacian ($`^2=_0^2+_1^2+_2^2`$) and rewrite it as

$$\underset{\mu =0,1,2}{}\left[\phi (x+\widehat{\mu })2\phi (x)+\phi (x\widehat{\mu })\right]=f(x)$$
(6)
where $`\widehat{\mu }`$ is a unit vector in the discretized space in direction $`\mu `$. Hence we solve it in $`\phi (x)`$ and obtain the following a recurrence relation

$$\phi (x)=\frac{_{\mu =0,1,2}\left[\phi (x+\widehat{\mu })+\phi (x\widehat{\mu })\right]f(x)}{6}$$
(7)
The following is a typical MDP program that solves eq. (1) by recursively iterating eq. (7). The program is parallel but there are no explicit call to communication functions:

```
00    #include "mdp.h"
01
02    void main(int argc, char** argv) {
03       mdp.open_wormholes(argc,argv);    // open communications
04       int L[]={10,10,10};               // declare volume
05       mdp_lattice      space(3,L);      // declare lattice
06       mdp_site         x(space);        // declare site variable
07       mdp_matrix_field phi(space,2,2);  // declare field of 2x2
08       mdp_matrix       A(2,2);          // declare matrix A
09       A(0,0)=1;  A(0,1)=I;
10       A(1,0)=3;  A(1,1)=1;
11
12       forallsites(x)                    // loop (in parallel)
13          phi(x)=0;                      // initialize the field
14       phi.update();                     // communicate!
15
16       for(int i=0; i<1000; i++) {       // iterate 1000 times
17          forallsites(x)                 // loop (in parallel)
18             phi(x)=(phi(x+0)+phi(x-0)+
19                     phi(x+1)+phi(x-1)+
20                     phi(x+2)+phi(x-2)-
21                     A*sin(2.0*Pi*x(1)/L[1]))/6;  // equation
22          phi.update();                  // communicate!
23       }
24       phi.save("field_phi.mdp");        // save field
25       mdp.close_wormholes();            // close communications
26    }

```

Notes:

* Line 00 includes the MDP library.
* Lines 03 and 25 respectively open and close the communication channels over the parallel processes.
* Line 04 declares the size of the box $`L=\{L_0,L_1,L_2\}`$
* Line 05 declares a lattice, called space, 3-dimensional, on the box $`L`$. MDP supports up to 10-dimensional lattices. By default a lattice object is a mesh with torus topology. It is possible to specify an alternative topology, boundary conditions and any parallel partitioning for the lattice. Notice that each lattice object contains a parallel random generator.
* Line 06 declares a variable site, called x, that will be used to loop over lattice sites (in parallel).
* Line 07 declares a field of $`2\times 2`$ matrices, called phi, over the lattice space. MDP is not limited to fields of matrices. It is easy to declare fields of any user-defined structure or class.
* Lines 08 through 10 define the matrix A.
* Lines 12 and 13 initialize the field phi. Notice that phi is distributed over the parallel processes and forallsites is a parallel loop.
* Line 14 performs communications so that each process becomes aware of changes in the field performed by other processes (*synchronization*).
* Lines 16 through 24 perform 1000 iterations to guarantee convergence. In real life applications one may want to implement some convergence criteria as stopping condition.
* Line 17 loops over all sites in parallel.
* Lines 18 through 21 implement eq. (7). Notice the similarity in notation. Here phi(x) is a $`2\times 2`$ complex matrix
* Line 22 performs *synchronization*.
* Line 24 saves the field. Notice than any field, including the user defined ones, inherit methods save and load from a basic class mdp\_field.
* It should also be noted that all MDP classes and functions are both type and exception safe. Moreover MDP components can be used without knowledge of C pointers and pointer arithmetics.

## 2 Linear Algebra

MDP includes a Linear Algebra package. The basic classes are:

* class mdp\_real, that should be use in place of float or double.
* class mdp\_complex, (just another implementation of complex numbers).
* class mdp\_array, for vectors and/or multidimensional tensors.
* class mdp\_matrix, for any kind of complex rectangular matrix.

The most notably difference between our linear algebra package and other existing packages is its natural syntax.

For example:

```
mdp_matrix A,B;
A=Random.SU(3);
B=exp(inv(A))*hermitian(A+5);

```

reads like

| $`A`$ and $`B`$ are matrices |
| --- |
| $`A`$ is a random $`SU(3)`$ matrix |
| $`B=e^{(A^1)}(A+5\mathrm{𝟏})^H`$ |

(8)
Notice that each matrix can be resized at will and is resized automatically when a value is assigned.

## 3 Lattice, Site and Field

An mdp\_lattice is a container for *topology* and *partitioning* information about the sites. In more abstract terms a lattice is any collection points (vertices) embedded in a multi-dimensional space and connected with directional links. The set of links determines the lattice topology and the boundary conditions. The term partitioning refers to the function that assigns each site (vertex) to one of the parallel process. A lattice, by default, is a mesh.

Each lattice is partitioned over the parallel processes at runtime. There is a default topology and default partitioning but it is possible pass any topology and partitioning functions to the mdp\_lattice constructor.

A lattice also contains a parallel random number generator: each site of each lattice has its own independent random number generator.

On each lattice it is possible to allocate one or more fields. Some fields are built-in, for example: mdp\_complex\_field, mdp\_vector\_field, mdp\_matrix\_field, etc. All of them extend (inherit from) mdp\_field$`<`$mdp\_complex$`>`$.

Class mdp\_complex can be used to declare any type of field. For example:

```
class W {
public: int w[10];
};
int L[]={30,30};
mdp_lattice plane(2,L);
mdp_field<W> psi(plane);

```

declares a $`30\times 30`$ lattice (plane) and a field (psi), that lives on the plane. The field variables of psi, psi(x) assuming x is an mdp\_site of plane, belong to class W.

Each user-defined field can be saved:

```
psi.save("filename");

```

loaded

```
psi.load("filename");

```

and synchronized

```
psi.update();

```

as any of the built-in fields.

## 4 Optimization Issues

Once an mdp\_lattice object is declared the constructor of class mdp\_lattice performs the following operations:

* Declares a parallel random number generator associated to each site (it uses the Marsaglia random number generator).
* Builds tables containing topology and partitioning information that will be used by the fields to optimize (minimize) communication. Basically each site determines which other sites are its neighbors and where they are located (on which parallel process). When a new field is created on the lattice the field will use these tables to create buffers for the communications.

Once a field object is declared the constructor of class mdp\_field (or derived field) performs the following operations:

* Each process allocates memory to store the local sites (i.e. sites that will be managed by the process itself).
* Each process loops over every other process and determines if the other process allocated sites that are neighbors of the local sites (this information is already stored in the tables maintained by the lattice object). If this occurs the two processes are said to *overlap*: they have sites in common that need to be synchronized.
* Each process allocates buffers to store copies of the sites that are not local but are neighbors of the local ones and need to be synchronized with the overlapping processes (in this paper we are assuming only next-neighbor synchronization but actually MDP supports also extended synchronization such as next-to-next-neighbor and more). Buffers are created according with some conditions: sites synchronized with the same overlapping process are stored contiguously in memory so that communication can be performed in a single *send/receive*. These buffers are created independently by each field.

Notice that two processes may be overlapping in respect to a given lattice and not overlapping in respect to a different lattice in the same program.

Every time a field changes, for example in a parallel loop such as

```
forallsites(x) phi(x)=0;

```

the program notifies the field that its values have been changed by calling

```
phi.update();

```

The method update performs all required communication to copy site variables that need to be synchronized between each couple of overlapping processes. These communications are optimal in the sense that:

* Each process, at each one time, is involved only in one send and one receive.
* Two different processes communicate only if they are overlapping in respect of the lattice associated to the field.
* If two different processes are overlapping, they perform a single send/receive of all sites variables that are synchronized between the two.
* Only the sending process needs to create a temporary buffer. The receiving process receives the site variables in the same buffer where they are normally stored without reordering (and without need for a temporary buffer).

We will refer to our set of communication rules as a “communication policy” (it is possible, in principle, to change this policy to deal with non-standard network solutions). Although our communication policy does not overlap communication with computation it has the advantage of minimizing network jam and calls to send/receive. Hence this communication policy is almost insensitive to network latency and is dominated by network bandwidth. Benchmarks are application dependent since parallel efficiency is greatly affected by the lattice size, by the amount of computation performed per site, processor speed and type of interconnection. In many typical applications, like the one described in the preceding example, the drop in efficiency is less than 10% up to 8 nodes (processes) and less than 20% up to 32 (our tests are usually performed on a cluster of Pentium 4 PCs (2.2GHz) running Linux and connected by Myrinet).

## 5 Conclusions

MDP is a powerful and reliable tool for developing efficient parallel numerical applications. Even if, on the one side, MDP is still undergoing development, on the other side, all of the features here described are fully functional and have been tested in real-life applications. For example MDP constitutes the core of the FermiQCD project  developed by the University of Southampton (UK) and Fermilab (Department of Energy). FermiQCD is collection of parallel algorithms for Quantum Chromo Dynamics computations. The typical FermiQCD problem is equivalent to solving iteratively a system of stochastic differential equations in a 4-dimensional space. Typical field variables are vectors of complex matrices. FermiQCD programs are used in production runs in parallel on 8 or more nodes.

We believe MDP could be a useful tool for scientists developing parallel numerical applications. MDP version 2.0 (current) is open source and is free for research and educational purposes.

Project web pages:

* http://www.pheonixcollective.org/mdp/mdp.html (license and source code)
* http://www.fermiqcd.net (Lattice QCD applications)