text:
# Flare Heating in Stellar Coronae

## 1 Introduction

The source of heating of solar and stellar coronae still eludes understanding even after decades of study (e.g., Schrijver et al. 1999). Despite significant evidence that magnetic activity is the prime driver for transferring energy into the corona, the mechanism by which this transfer occurs is not established in either the case of the Sun or other stars (see e.g., Rosner, Golub, & Vaiana 1985, Narain & Ulmschneider 1996). Numerous heating mechanisms, such as acoustic wave dissipation (Stepien & Ulmschneider 1989), Alfvén wave dissipation (Cheng et al. 1979, Narain & Ulmschneider 1990), magnetic reconnection phenomena (Parker 1988, Lu & Hamilton 1991) have been proposed, all of which might play some role in the overall heating.

Recent work in the solar case has lent strong credence to the possibility of coronal heating being dominated by small-scale explosive events suggestive of Parker’s nanoflare model, which is based on magnetic reconnections releasing energies $`10^{24}`$ erg event-1. It has been well known that solar microflares and milliflares111 Because of the vast range of flare energies encountered, the energy ranges of the different flare types are not well defined. We adopt the convention (see e.g., Aschwanden et al. 2000) that milliflares cover the range $`E10^{2932}`$ ergs, microflares $`E10^{2629}`$ ergs, and nanoflares $`E10^{2326}`$ ergs. We consider all events down to the microflare regime to be ‘normal’ X-ray flare events, with similar origin, parameters, and effects, except for the differences in energy deposition.  are distributed in number as power laws of their energy output (Lin et al. 1984, Hudson 1991),

$$\frac{dN}{dE}=kE^\alpha $$
(1)
where $`E`$ is the energy of the flare and $`k`$ is a constant. This relation has been verified and extended to lower energies by various authors, but despite the near universal acceptance of the form of the function in Equation 1 (e.g., Crosby, Aschwanden, & Dennis 1993; see Kopp & Poletto 1993, Shimizu & Tsuneta 1997 for a different perspective), neither the index $`\alpha `$ nor the normalization $`k`$ are well determined. For instance, $`\alpha =1.61.8`$ in the HXR to microflare energy range, and is variously measured to lie in the range 1.8-2.9 at lower energies (Shimizu 1995 \[$`\alpha =1.51.6`$\], Porter, Fontenla, & Simnett 1995 \[$`\alpha =2.3`$\], Krucker & Benz 1998 \[$`\alpha =2.32.6`$\], Parnell & Jupp 2000 \[$`\alpha =2.02.6`$\], Aschwanden et al. 2000 \[$`\alpha =1.8`$\], Winebarger et al. 2001 \[$`\alpha =2.9\pm 0.1`$\], Veronig et al. 2002 \[$`\alpha =2.03\pm 0.09`$\]). Recently Aschwanden & Parnell (2002) have used scaling laws based on energy balance arguments to conclude that $`\alpha `$ must be $`1.6`$ on the Sun. The precise value of $`\alpha `$ is of considerable interest because if the power law is steep enough ($`\alpha >2`$), then in principle a multitude of small impulsive events would be sufficient to account for the energy output of the entire corona.

Here we reconsider in particular an outstanding question in stellar X-ray astronomy, which is the nature of the apparently quiescent emission from active stars: does this emission actually arise from a superposition of a multitude of impulsive events (such as milliflares and microflares), or from truly quiescent plasma? Previous work based on detecting flares in EUV data (see e.g., Audard et al. 2000) suggests that flare contribution is indeed an important factor. Further, correlations of quiescent X-ray flux with time-averaged U-band flare flux (Skumanich 1985, Doyle & Butler 1985) and the synchrotron radio luminosity (Güdel & Benz 1993), together with the similar correlations found in the solar case (Benz & Güdel 1994) strongly suggest a link between the apparently quiescent emission and flares. In addition, spectroscopic evidence for high temperature plasma ($`T10^7`$ K) during the quiescent phase (Butler et al. 1986, Kashyap et al. 1994, Drake 1996, Güdel et al. 1997, Güdel 1997) indicates that this quiescent emission could in fact be very similar to flare emission in origin. Thus, apparently quiescent coronae of active stars could be composed of a continuum of small unresolved flares, presumably distributed as power laws analogous to the Sun. This view is also supported by the double-peaked Differential Emission Measures (DEMs) that result when an ensemble of flaring, hydrodynamically evolving loops are modeled on active solar analogs (Güdel et al. 1997, Güdel 1997).

The possibility of stellar coronal heating due to small flares was considered by Ambruster, Sciortino, & Golub (1987) who searched for variability in Einstein data of active stars and discussed the contribution of low-level flaring to heating stellar coronae. They concluded that while flaring must contribute at some level, the evidence does not justify extending the solar power-law distributions to the stellar microflare case. Later studies of ensembles of strong stellar flares seen with EXOSAT and EUVE have shown these are distributed as power laws with index $`\alpha =1.61.8`$ (Collura et al. 1988, Pallavicini et al. 1990, Osten & Brown 1999), thus ruling out low-intensity flares as a significant contributor to the heating budget. In contrast, using a more sensitive method to detect fainter flares (see Crawford et al. 1970), Robinson et al. (1995, 1999, 2001) find that for stellar chromospheric and transition region events observed with the high-speed photometer and the imaging spectrograph on the HST, $`\alpha 1.762.17`$ in the chromosphere of the active dMe star CN Leo; $`\alpha =2.25\pm 0.1`$ in the chromosphere of the dMe star YZ CMi; and $`\alpha 2.22.8`$ in the transition region of the dM0e flare star AU Mic. (Note however that chromospheric and transition-region flare distributions have no known direct correspondence with the coronal case.) Applying a similar method to EUVE/DS data, and also correcting for overlaps in flares, Audard et al. (1999) find that for solar analogs EK Dra and 47 Cas, $`\alpha 2.2\pm 0.2`$. This analysis was further extended by Audard et al. (2000) to a larger sample of cool stars, and they find that $`\alpha `$ ranging from 1.5 to 2.6, with the majority of the measurements having $`\alpha >2`$. Similar results are obtained for AD Leo (Güdel et al. 2001,2002).

Note that the above studies are limited to relatively large flares ($`E10^{31}`$ ergs) because of instrument sensitivity, and also because the more numerous weaker flares are harder to detect in the presence of “contamination” by other weak flares. Thus, the low-energy end of the flare distribution is subject to large uncertainties.

We have developed a new method to model the undetectable stellar flares and thus derive estimates of flare indices covering the milliflare regime as well as directly estimating the flare contribution to the observed flux. We apply this method to active low-mass stars FK Aqr, V1054 Oph, and AD Leo. The datasets used are described in §2. The analysis method is detailed in §3 (a glossary of the terms used is given in Appendix A). The results are given in §4, and are summarized in §5.

## 2 Data

Here we analyze data from the Extreme Ultra-Violet Explorer satellite’s Deep Survey photometer (EUVE/DS) of 3 active low-mass stars.222 The intrinsic EUVE time resolution is 8 ms, and this is adequate to resolve the sources even at the maximum count rate seen in our observations (3 ct s-1). The EUVE/DS covers a useful spectral range of 52-246 Å, with a peak effective area of 28 cm-2 at 91 Å. This wavelength range includes many lines from highly ionized Fe XVIII to Fe XXIII normally found in the coronae of active stars, in addition to bound-free and free-free continua: plasma temperatures from $`1`$ to 30 MK are thus accessible for observation.  These stars are known to have significant flare activity, and do not undergo eclipses, and so are amenable to straightforward modeling:

##### FK Aqr

is a BY Dra type, spectroscopic, double-lined, non-eclipsing, low-mass, active binary (Table 1). Its flare energy output has remained steady over long intervals (8 years; Byrne et al. 1990), and optical modulation due to spot activity has been observed. Thus it is possible that flare heating could be a significant component of coronal emission for this star. Indeed, the EUVE/DS light curve shows evidence of a number of flare events (Figure 1) as well as large stochastic variability in the apparently quiescent emission.

##### V1054 Oph

is a low-mass astrometric binary (Table 2). Previous X-ray observations with ROSAT have confirmed the existence of almost continuous flaring (Kellet & Tsikoudi 1999) and the existence of high-temperature plasma that is responsible for most of the intensity variations (Giampapa et al. 1996). EUVE data show a large number of relatively weak flare events that blend into the variations in the quiescent emission (Figure 2). Such a dataset is very difficult to analyze by the traditional means of detecting and counting flares, but poses no difficulty to a direct modeling approach as is described below.

##### AD Leo

is a well studied low-mass single flare star (Table 3) with a high flare rate. A long duration exposure was obtained by Güdel et al. (2001,2002) that shows many large flares (Figure 3). The data were obtained in 6 segments, and since the first segment could not be optimally reduced, and the last segment suffered from a high background, we have ignored them in this analysis and have concentrated on the $`2^{nd}5^{th}`$ segments. In particular, we have carried out the analysis with the data grouped into two sets, segments 2 and 3 forming one set and segments 4 and 5 forming the other. The light curves show slightly different characters in the two parts, with the former part dominated by large flares while the latter part shows smaller identifiable flares (Figure 3).

In our analyses (see §3 below), we use the photon arrival times directly, and apply the dead-time and Primbsching corrections333 Primbsching refers to the photons lost due to telemetry bandwidth, and is measured by the ratio of the total counts incident in a quadrant of the detector (as determined on the spacecraft; the summed counts from all the quadrants is used to determine the instrument deadtime) to the number of events telemetered to the ground (see the EUVE Data Products Guide for a detailed description of its origin and correction). for the particular observation to the models (see §3.1) over good time intervals (GTIs) defined to exclude SAA passages and Earth blockages. That is, we process model light curves to produce a simulated set of photons and apply time windowing and statistical censoring (i.e., the instrument response in the time domain) to generate an event list that may be directly compared to the observed events. For the sake of brevity, these corrections will henceforth be referred to as “instrumental” corrections. The source photons are collected over a circle of radius $`4^{\prime \prime }`$ surrounding the source. In all cases, $`10\%`$ of the events are estimated to be due to background photons, except in the case of V1054 Oph where it is estimated to be 17%.

## 3 Analysis

Previous attempts to determine the values of the parameters in Equation 1 have concentrated on first detecting flares in the binned light-curve and then fitting power law expressions to the detected numbers. In contrast, we assume the reality of a power law distribution and set up a model to compare with the data. This model is described in §3.1, and the manner in which the model parameters are derived is described in §3.2. The applicability of the method, including verification, assumptions, advantages and disadvantages, are discussed in §3.3.

### 3.1 The Model

Because flares generally occur randomly (see §3.3.2), we cannot directly model the light curves, as it is not possible to “deconvolve” a complex light curve by specifying the location of each flare in a model (see, e.g., Figure 2). Instead, rather than match the flare locations in detail, we carry out a fitting process wherein only the number and intensity distributions of a set of model flares are matched with the data. This is accomplished by comparing the distributions of photon arrival-time differences. The assumptions made in defining the model described below are discussed in detail in §3.3.2.

We first generate a set of photon arrival times by simulation from a 3-parameter model

$$𝐌=\{\alpha ,r_F,r_C\},$$
(2a)
where $`\alpha `$ is the index of the power law as in Equation 1, $`r_F`$ is the average count rate due to flares, and $`r_C`$ is a constant component which is also expected to fully account for the background (see §3.3.2). It is also possible to use the average total flux $`r_T=r_F+r_C`$ ($`r_F<r_T`$) as the defining parameter instead of $`r_C`$,
$$𝐌^{}=\{\alpha ,r_F,r_T\},$$
(2b)
with equivalent results. The counts at time $`t`$ in an interval $`dt`$, $`C(t)Poisson[r(t)dt]`$ are Poisson-distributed according to the instantaneous rate $`r(t)`$. The rate $`r(t)`$ may be described as due to the sum of model counts due to a flare component $`f(t)`$ and a non-flare component $`r_C(t)`$, and a correction factor $`\varphi (t)`$ that takes into account Primbsch, dead-time, and GTIs,444 Note that $`0\varphi (t)1`$ determines the probability that any photons are collected at time $`t`$, and in particular is identically 0 outside the GTIs. , i.e.,

$$C(t)\varphi (t)Poisson[r_C(t)dt+f(t)dt],$$
(3a)
where $`r_C(t)`$ is taken to be constant unless stated otherwise. The flare component $`f(t)`$ may in turn be written as a superposition of numerous individual flares, i.e.,
$$f(t)=\underset{j=1}{\overset{N_f}{}}\mathrm{\Theta }(tt_j)F_je^{(tt_j)/\tau },$$
(3b)
where $`\tau `$ is the flare decay timescale,555 We assume $`\tau `$ to be fixed for purposes of simplicity. See §3.3.2 for a discussion of cases when $`\tau `$ may vary.  $`N_f`$ are the total number of flares, $`F_j`$ are the peak energy intensities of individual flares (the counts due to each flare, $`c=\tau F_j`$ are sampled from the distribution represented by Equation 1), and
$$\mathrm{\Theta }(x)=\{\begin{array}{cc}\hfill 0:& x<0\hfill \\ \hfill 1:& x0\hfill \end{array}$$

is a step function to represent flare onset.
Note that not only will the placement $`t_j`$ and peak intensity $`F_j`$ of the flares vary for each simulation, but so will the total number of flares $`N_f`$. Within the bounds of Poisson statistics, we expect that for any given simulation, the total energy due to all the simulated flares is determined by the average expected count rate due to flares and the total duration of the observation, $`\mathrm{\Delta }T`$, i.e.,

$$\underset{0}{\overset{\mathrm{}}{}}𝑑tf(t)=\underset{j=1}{\overset{N_f}{}}F_j\tau =r_F\mathrm{\Delta }T.$$
(3c)
Note that Equation 1 is written as a function of the energy deposited by a flare $`E`$, but assuming that the observed counts due to this flare $`c`$ is proportional to $`E`$ (see discussion in §3.3.2), i.e.,

$$\frac{dN}{dE}\frac{dN}{dc}=\kappa c^\alpha dc,$$
(4)
we can use the model parameter $`r_F`$ to fix the normalization $`\kappa `$ of the power law. By equating the total counts due to the flare component with the counts expected from the power law distribution, we get

$$r_F=\frac{\underset{c_{min}}{\overset{c_{max}}{}}𝑑cc\frac{dN}{dc}}{\mathrm{\Delta }T}.$$
(5)
The upper limit in the integration is defined by requiring that all the observed counts be due to a single model flare ($`c_{max}=\underset{j}{\mathrm{max}}\{F_j\tau \}=r_F\mathrm{\Delta }T`$), that is, no flare model may produce a flare with more counts than are observed. The lower limit in the integration is defined by requiring that each flare be assigned at least 2 counts (i.e., $`c_{min}=2`$; this is so that an arrival time difference may be determined even in the extreme case where the model may have just one weak flare – see §3.2 below). Thus, carrying out the integral in Equation 5 and rearranging the terms,

$`\kappa |_{(\alpha 2)}`$ $`=`$ $`{\displaystyle \frac{(2\alpha )r_F\mathrm{\Delta }T}{(r_F\mathrm{\Delta }T)^{2\alpha }2^{2\alpha }}}`$ (6a)
$`\kappa |_{(\alpha =2)}`$ $`=`$ $`{\displaystyle \frac{r_F\mathrm{\Delta }T}{\mathrm{ln}(r_F\mathrm{\Delta }T/2)}}.`$ (6b)
For $`\alpha >2`$, this implies that if $`r_Fr_T`$, then the adopted lower limit is very close to the theoretical lower limit to the extent of the power law distribution in order for it to account for all the observed counts (see Table 4).

In order to obtain the best values of the parameters that describe the data, and a confidence range on these parameters, we carry out a forward-fitting procedure based on a Bayesian formalism (see e.g., Loredo 1990 for a tutorial on the foundations of Bayesian probability theory): we compute the probability distribution of the model parameters given the data, which is a composite of whatever prior information we may have on the model parameter values and the likelihood of realizing the observed data for specified parameter values. That is, we derive the joint posterior probability $`p(𝐌|D,I)`$ of the model parameters conditional on the data,666 The expression $`p(x)`$ represents the probability that the logical statement “$`x`$” is true. In particular, conditional statements are written as “$`A|B`$”, i.e., $`p(A|B)`$ represents the probability that statement “$`A`$” is true given that statement “$`B`$” is true. These statements may be generalized to include models and parameter values; as an illustrative example, $`p(\{\alpha =2.1,r_F=0.1,r_T=0.3\}|D,I)=0.1`$ reads as “the probability is 0.1 that $`\alpha =2.1,r_F=0.1,\mathrm{and}r_T=0.3`$ given the data $`D`$ and supporting information $`I`$.”

$`p(𝐌|D,I)`$ $``$ $`p(\alpha |I)p(r_F|I)p(r_C|I)`$ (7a)
$`\times p(D|𝐌,I),`$
where $`D`$ represents the data, and $`I`$ represents assumptions necessary to solve the problem, including the effects of instrument characteristics. The first 3 factors on the right hand side of the equation are the a priori probability distribution functions of the model parameters, and the last factor, $`p(D|𝐌,I)`$ is the likelihood of obtaining the observed data given the model parameters. Computing the posterior probability distribution constitutes a complete solution to the inference problem for the specified model M. Note that the above is a simplified form of Bayes’ Theorem wherein the normalization factor $`p(D|I)`$ usually present on the right hand side of the equation is ignored (e.g., Kashyap & Drake 1998, van Dyk et al. 2001). If $`r_T`$ is used instead of $`r_C`$, the joint probability distribution takes the form
$`p(𝐌^{}|D,I)`$ $``$ $`p(\alpha |I)p(r_F|r_T,I)p(r_T|I)`$ (7b)
$`\times p(D|𝐌^{},I),`$
with $`p(r_F|r_T,I)=0`$ for $`r_F>r_T`$. In the following, we make no distinction between $`𝐌`$ and $`𝐌^{}`$. The Bayesian formalism allows us to determine the probability distributions of each of the parameters by marginalizing, i.e., integrating over the other parameters. Thus we can write for the probability distribution of $`\alpha `$ alone,

$$p(\alpha |D,I)=_{r_F}𝑑r_F_{r_C}𝑑r_Cp(𝐌|D,I)$$
(8)
and similarly for the other parameters (cf. Equation 11).

### 3.2 The Algorithm

The problem facing the modeling process is illustrated in Figure 4, where the light curve from the observation of FK Aqr is compared with selected simulated model light curves.777 We emphasize that these light curves are shown only for purposes of illustration, and that the analysis does not require binning the observed events.  It is easy to recognize that the model light-curve for $`\alpha =2.5`$ is the most similar in character (i.e., in the number and strength of discernible flares) to the data light-curve, but clearly the locations of the flares do not match. Normal fitting methods that rely on matching the expected model counts in a bin with the observed counts would fail on a problem such as this. We thus seek to employ a method which compares the interesting information between the data and the model without being misled by the obvious, though uninteresting, differences.

One method that would satisfy our requirements of simultaneously ignoring flare locations and yet be sensitive to flare numbers and intensities is to compare the distributions of photon arrival-time differences $`g(\delta t)`$ (where $`\delta t`$ is the interval between consecutive counts) between the data and the model.888 Other methods such as computing the fractal length of the events, multi-scale analyses of light curves, etc., may also be applied (see e.g., Vlahos et al. 1995). A detailed comparison of the benefits of one method versus another is beyond the scope of this article, but note that the results we derive here are robust within the regime of applicability of the adopted method (see §3.3).  In the absence of any variability, i.e., if the light curve is flat with expected rate $`r`$, the resulting set of $`\delta t`$ are distributed as an exponential,

$$g(\delta t)re^{r\delta t}.$$
(9a)
In the presence of variability in the expected rate, the observed distribution would be a superposition of many such distributions: if the fraction of time that a source spends at rate $`r_i`$ is given as $`\rho _i`$, then
$$g(\delta t)\underset{i}{}\rho _ir_ie^{r_i\delta t}.$$
(9b)
It is thus possible to distinguish between different magnitudes of variability. In particular, events generated from a model with low $`\alpha `$ (e.g., $`\alpha =1.8`$, where the light curve would be dominated by a few very large flares) would result in a distribution $`g(\delta t)`$ that is skewed to smaller values of $`\delta t`$, while those from a model with large $`\alpha `$ (e.g., $`\alpha =3.0`$, where the light curve would be composed of a large number of very small flares that overlap each other so finely that it would not be possible to distinguish it from a source with constant intensity) would approach the limiting case of Equation 9a above.

In Figure 5 we show the comparison between the distributions of arrival-time differences derived from the same datasets in Figure 4. The different curves may be compared using any of a number of statistical methods such as computing the $`\chi ^2`$, applying the Kolmogorov-Smirnoff test, etc. This approach succeeds in providing us with an objective measure of the “similarity” between the datasets; indeed of the 4 models considered in Figures 4 and 5, the one with $`\alpha =2.5`$ has the smallest $`\chi ^2`$ when compared with the data. Note that Figure 5 also illustrates a fundamental limitation of this method, viz., the method loses sensitivity for larger values of $`\alpha `$, which may be indistinguishable from sources with a constant intensity (see §3.3.1).

Because the model (Equation 3a) is stochastic, we use Monte-Carlo simulations to generate many realizations for each set of model parameters. The simulated distributions of arrival-time differences $`g_{MODEL}(\delta t)`$ are compared with the corresponding distribution derived from the data $`g_{DATA}(\delta t)`$ over a parameter grid. The likelihood is computed as the probability density of obtaining the observed $`\chi ^2`$ value for N degrees of freedom (see Eadie et al. 1971, their Equation 4.22):

$$p(D|𝐌)=\frac{\frac{1}{2}\left(\frac{\chi }{2}\right)^{\frac{N}{2}1}e^{\frac{\chi }{2}}}{\mathrm{\Gamma }(\frac{N}{2})}.$$
(10)
The a priori probability distributions for the parameters (Equation 7a) are taken to be non-informative, and thus flat, over the range of parameter values defining the grid. The basic steps in the algorithm we follow are outlined below:

1. From the data, derive the distribution of photon arrival-time differences, excluding the gaps in the data due to breaks in the GTIs.
2. For the specified values of model parameters $`𝐌`$ (Equation 2), obtain a realization of the photon event-list over the duration of the observation. This is done by first simulating a light curve incorporating all the flare events (Equation 3b), added to the base emission, and then deriving photon arrival times based on the instantaneous count rates.
3. Because we are to compare EUVE/DS event list data with simulated data, instrument effects must be taken into account. This is encoded in the factor $`\varphi (t)`$ (see Footnote 4). We apply Primbsching and dead-time corrections by discarding photons with probability $`1\varphi (t)`$ (i.e., sample a random number $`z`$ from a uniform distribution over the interval $`[0,1]`$, and discard the photon at time $`t`$ if $`z>\varphi (t)`$). The retained set of events is identical in its instrumental characteristics to the observed data.
4. From this set of simulated photon arrival times, derive the model distribution of arrival-time differences over a binning that maximizes the reduced $`\chi ^2`$.999 We do not know the optimal binning for $`g(\delta t)`$ a priori. The binning must be chosen such that differences between the two distributions being compared are highlighted to best advantage, at a scale that is determined by the datasets themselves. Because the total number of photons in the datasets being compared are approximately the same, very coarse binnings and very fine binnings both produce low values of the reduced $`\chi ^2`$, and we adopt as the optimal binning that which maximizes this value.
5. Compute the likelihood as in Equation 10 and the a posteriori probability at the specified grid point as in Equation 7a. Note that the number thus obtained is not normalized, and so must not be used to compare, for example, the relative probabilities of different types of models.

In practice, the above algorithm must be enhanced by some additional steps. For instance, the likelihood may be artificially reduced if a large model flare is fortuitously located coincident with large dead-time. We therefore shift the Primbsching and dead-time corrections by a random interval and recompute the arrival-times. This is mathematically equivalent to shifting the simulated events, but is done in this fashion because of the lower computational cost. Typically 3 such shifts are carried out for each simulated light curve, and the best comparison is chosen. In addition, in order to derive a robust estimator, we perform $`1520`$ simulations for each $`𝐌`$ (or $`𝐌^{}`$) and adopt the median value of the set as the final value of $`p(𝐌|D,I)`$. The grid of model parameters are chosen such that $`\alpha `$ is explored over the useful range of the algorithm ($`1.5\mathrm{}3.0`$; see §3.3.1); $`r_T`$ in a range within 3$`\sigma `$ of the average count rate; and $`r_C`$ and $`r_F`$ ranges from $`0`$ to the average count rate. There are typically $`15`$ grid points for each parameter.

The probability distribution along any of the axes is then obtained by summing over the other axes and normalizing, e.g.,

$$p(\alpha |D,I)=\frac{\underset{r_F}{}\underset{r_C}{}p(𝐌|D,I)}{\underset{\alpha }{}\underset{r_F}{}\underset{r_C}{}p(𝐌|D,I)}.$$
(11)
The derived probability distributions may be summarized by their means and variances,101010 Other well-known methods of summarizing probability distributions include noting the MAP (maximum a posteriori value; the mode of the distribution), as e.g., $`p(\alpha _{_{MAP}}|D,I)=\underset{\alpha }{\mathrm{max}}\{p(\alpha |D,I)\}`$, or a range corresponding to an integrated area under the curve, equivalent to some defined probability $`\pi =\underset{\alpha _{min}}{\overset{\alpha _{max}}{}}p(\alpha |D,I)`$, such that $`\alpha _{_{MAP}}[\alpha _{min},\alpha _{max}]`$, etc. Unless otherwise specified, we always report the mean values and $`1\sigma `$ errors.  e.g.,

$`\overline{\alpha }={\displaystyle \frac{\underset{\alpha }{}\alpha p(\alpha |D,I)}{\underset{\alpha }{}p(\alpha |D,I)}},`$ (12a)
$`var(\alpha )=\overline{\alpha ^2}(\overline{\alpha })^2.`$ (12b)
### 3.3 Applicability

As was demonstrated above (§3.2, Figures 4 & 5), the distributions of arrival-time differences $`g(\delta t)`$ provide an objective means to compare event lists that are dominated by stochastic events. Here, we verify that the algorithm gives reasonable results by creating simulated datasets and obtaining best-fit values for them (§3.3.1), then discuss the effects of some of the assumptions we have made in formulating the problem (§3.3.2), and detail the advantages and disadvantages of the adopted method (§3.3.3).

#### 3.3.1 Verification

We show here that the algorithm works as expected by simulating well defined datasets and then “fitting” to them. We generated flare models with specific values of $`\alpha `$ (labeled $`\alpha _{_{TRUE}}`$) ranging from 1.5 to 3, and for the sake of simplicity, held $`r_F=0.5`$ ct s-1 and $`r_C=0.1`$ ct s-1. Simulations were carried out over a time interval $`\mathrm{\Delta }T=100`$ ks and assuming a fixed flare decay timescale of $`\tau =3`$ ks. These values were chosen as being typical of EUVE/DS observations. The datasets to be fit to were chosen by first simulating 9 separate event lists and then choosing the one with the median number of counts, in order to avoid contaminating the verification process with extremal datasets. The fits were carried out as described in §3.2, and the resulting best-fit values $`\alpha _{_{FIT}}`$ are shown in Figure 6. The fitted values follow the true values closely, for all 3 parameters.

The method works well over the range of interesting $`\alpha `$, that is, spanning the value $`\alpha =2`$, and is capable of distinguishing between sources with $`\alpha `$ above or below this critical value (see §1). The sensitivity of the method naturally decreases as $`\alpha `$ increases to $`3`$ when the datasets become indistinguishable from that of a steady source. However, this is not a hard limit and can be extended for datasets with larger average flare component intensity $`r_F`$ and longer observations $`\mathrm{\Delta }T`$. The reliability also decreases as $`\alpha `$ decreases to $`1.5`$ since at these values of $`\alpha `$ the simulations are dominated by a very few but very large flares and are therefore subject to large fluctuations, and hence are not robust. A larger number of simulations at each grid point becomes necessary when a dataset exhibits smaller values of $`\alpha `$.

#### 3.3.2 Assumptions

We have made a number of simplifying assumptions in our analysis, and these are discussed below, with particular attention to the effect they have towards the robustness of the results. In general, all our assumptions are conservative, in the sense that they all act to generate a best-fit $`\alpha `$ and $`r_F`$ that is smaller than the true $`\alpha `$ and $`r_F`$. That is, the main results expounded here, that the flare distributions on active stars have $`\alpha >2`$, and that the apparently quiescent emission is dominated by flares, are not affected.

##### Power laws:

Energy release due to flare events in both solar and stellar environments has been well established to be highly intermittent and that the events are distributed as power laws spanning many decades in energy (e.g., Crosby et al. 1993, Güdel et al. 1997, Güdel 1997, Krucker & Benz 1998, Osten & Brown 1999, Audard et al. 2000, Aschwanden et al. 2000, Veronig et al. 2002). This suggests the absence of a characteristic scale for the intensity of a flare event, and is understood to arise from avalanche or SOC (self-organized critical) models (Lu & Hamilton 1991, Vlahos et al. 1995, Georgoulis & Vlahos 1998, Krasnoselskikh et al. 2001). Nevertheless, there is evidence from studies of solar flares that the index of the power law distribution does not remain the same in the microflare and nanoflare range (cf. Aschwanden et al. 2000 \[their Figure 10\], Winebarger et al. 2001). Studies of stellar flares (e.g., Ambruster et al. 1987) also suggest that more than one type of plasma instability may be present, and that the distribution may steepen or change at lower flare energies; for instance, compare $`\alpha 1.8`$ found by Osten & Brown (1999) with the generally larger values found by Audard et al. (2000) who include much weaker flares in their analysis (also see the results from AD Leo presented here in §4). Studies of flare distributions arising from the transition region also show a similar dichotomy (Robinson et al. 2001). Such changes in power law indices could arise from a variable driving mechanism (Georgoulis & Vlahos 1998). It is therefore possible that the true distribution departs from exact self-similarity in some complicated manner. However, present data, especially in the case of stellar flares, are insufficient to detect these departures (Audard et al. 1999,2000; also Güdel et al. 2001,2002). Here we assume that a single power law index is valid over at least 4 orders of magnitude (see Table 4) in flare energy. If the distribution does steepen for lower flare intensities, then note that first, the steepest parts, which approach the limiting case of constant emission (see Equation 9a), will contribute to enhancing the constant intensity $`r_C`$ (as described in §3.2). Hence the fitted $`r_F`$ will be a lower limit to the true value. Second, the fitted $`\alpha `$ will be a weighted average biased towards the high count rate, shallower distribution, and inasmuch as a “true” value of $`\alpha `$ may be said to exist, it would be greater than the fitted value.

##### Decay Timescales:

We analyze the data by fitting the 3 parameters $`\alpha ,r_F,`$ and $`r_C`$ of the model (Equation 2). We assume that the flare decay timescale $`\tau `$ (Equation 3b) is fixed (usually at 3 ks, as suggested by the detectable flares, and the radiative cooling timescales suggested by ROSAT observations; see e.g., Giampapa et al. 1996) and is the same for all flares. However, it is well known that $`\tau `$ varies for individual flares on active stars (e.g., Pallavicini et al. 1990) and especially so for RS CVn stars (Osten & Brown 1999). Note though that we confine ourselves to a specific class of active stars – low-mass main-sequence stars that flare frequently – and exclude RS CVn stars from our sample. There is evidence that flare duration scales with flare energy in various passbands (e.g., Crosby, Aschwanden, & Dennis 1993; also, Vlahos et al. 1995 for avalanche models, Temmer et al. 2001 for Hα flares, Georgoulis, Vilmer, & Crosby 2001 for deca-keV flares, Robinson et al. 2001 for chromospheric and transition region events, Veronig et al. for GOES SXR flares) such that more intense flares appear to last longer. Based on avalanche model simulations, Lu et al. (1993) find that on average, an event with decay timescale $`\tau `$ corresponds to an energy release $`E\tau ^{1.77}`$. (See also the discussion below about sympathetic flaring, which could affect measurements of flare durations.) But note that the decay timescales for soft X-ray flares (such as the ones we are concerned with) are primarily dependent on flare temperature and plasma density, and secondary effects such as changes in the heating rate come into play only for very large flares where the heating timescales approach or surpass the radiative cooling timescales. In a model that incorporates such variations of $`\tau `$, the fainter flares will have larger peak rates than in the regular model (in order to have the same total energy output), and the distribution of arrival-times $`g(\delta t)`$ will be skewed towards smaller $`\delta t`$. Thus, fitting to data that may be generated in this manner using a $`\tau `$ fixed by the higher intensity flares would result in $`\alpha _{_{FIT}}<\alpha _{_{TRUE}}`$ (see Figure 5) and therefore the fitted values would be lower limits to the true indices. Conversely, fitting a model where $`\tau `$ decreases with flare energy to a dataset that is not generated in this manner will result in $`\alpha _{_{FIT}}>\alpha _{_{TRUE}}`$. (See also discussion by Güdel et al. 2001,2002.) In order to test the sensitivity of our datasets to these variations in modeling, we have carried out fits to the data using models with an energy-dependent timescale, $`\tau E^\beta `$, in particular, $`\beta =\frac{1}{4}`$ (Güdel et al. 2002, based on fits to EXOSAT flare decay timescales of Pallavicini et al. 1990). We find that the best-fit value of $`\alpha `$ increases by $`0.30.4`$ for the complex model, e.g., for segment 3 of AD Leo $`\alpha `$ changes from 2.2 to 2.7 . We thus conclude that using the simpler model ($`\tau =constant`$) is preferable in that we do not overestimate $`\alpha `$, and further note that stronger dependences, e.g., that suggested by the theoretical study of Lu et al. (1993; $`\tau E^{0.56}`$) would result in even larger fitted values of $`\alpha `$. In addition, we have also explored the sensitivity of fits to the adopted value of $`\tau `$, by generating simulated data with small decay timescales $`\tau _{_{SIM}}`$ and then fitting to it a model with a larger decay timescale $`\tau _{_{FIT}}`$. We find that as expected the best-fit $`\alpha `$ is smaller than the true value, i.e., $`\alpha _{_{FIT}}|_{(\tau _{_{FIT}}>\tau _{_{SIM}})}<\alpha _{_{SIM}}|_{\tau _{_{SIM}}}`$.111111 This effect is large only for small values of $`\alpha `$. For example, $`\{\alpha _{_{SIM}}=1.7;\tau _{_{SIM}}=1000\}`$ is fit by $`\{\alpha _{_{FIT}}=1.5\pm 0.02;\tau _{_{SIM}}=3000\}`$, whereas $`\{2.1;500\}_{_{SIM}}`$ is fit by $`\{2.16\pm 0.17;3000\}_{_{FIT}}`$ (i.e., is indistinguishable from the true value), etc. Note that, as expected, this exercise also shows that when the data are dominated by multiple overlapping flare events, the results are insensitive to the precise value of the decay timescale, and effects such as stellar rotation may be ignored.  The reverse holds true for the opposite case, i.e., $`\alpha _{_{FIT}}|_{(\tau _{_{FIT}}<\tau _{_{SIM}})}>\alpha _{_{SIM}}|_{\tau _{_{SIM}}}`$, in which case the best-fit values imply flare distributions that are steeper than they should be. Thus, the relatively large value we adopt for $`\tau `$ (3 ks) provides a conservative estimate, and the fitted values may be considered lower limits to the true values.

##### Rise Times:

We model the flares as having an instantaneous rise and a slow exponential decay (Equation 3b). In general this is a good approximation, since the rise times are short compared to the flare duration (e.g., Reale, Peres, & Orlando 2001, Temmer et al. 2001). Further, our algorithm is insensitive to the exact sequence of the emission intensities (i.e., the process of forming $`g(\delta t)`$ destroys the sequential information in the event lists), and thus any discernible flare rise times will be indistinguishable in their effects from flares of small decay times. Thus, this approximation reduces to the problem discussed above, that of a sequence of flares that are “contaminated” by flares with smaller decay timescales, with similar effects on $`\alpha `$ and $`r_F`$.

##### Flare Waiting-times:

Much work has been done to characterize the time between flares within individual active regions. If we consider each flare to be an independent event then the waiting time between flares is a Poisson process (Rosner & Vaiana 1978). However, actual observations of solar Hard X-ray flares show that the waiting-time distribution (WTD) is a power law in intervals, with index ranging from $`2.16`$ to $`2.4`$ (Boffetta et al. 1999, Wheatland 2000, Lepreti, Carbone, & Veltri 2001), though Moon et al. (2001) find that for strong solar flares (strength greater than C1) the waiting times are indeed well characterized by a Poisson distribution. The power law distribution has been interpreted as due to sympathetic flaring (i.e., a cascade of small flares that follow a large flare, thus invalidating the assumption of event independence) by Wheatland, Sturrock, & McTiernan (1998), and adapted within SOC models as a non-stationary random process (Norman et al. 2001). However, since we can only observe disk-integrated flux in stellar data (it is not possible to monitor individual active regions), event independence is a better approximation. We thus assume in our modeling that the stellar flare WTD is Poisson. In such a case, the model undercounts the number of flares separated by short intervals (Wheatland et al. 1998). These flares would generally be associated with the stronger flares that set off a cascade, thus effectively increasing the decay timescale for large flares. As argued above, this causes $`\alpha _{_{FIT}}<\alpha _{_{TRUE}}`$.

##### Energy Deposition:

We have implicitly assumed that the observed counts track the energy deposited by the flares linearly (see Equations 1 and 4). That is, we assume that the energy deposition process that causes the flare event (generally considered to be magnetic reconnection – see Parker 1988) is, first, sufficiently energetic that a large fraction of the deposited energy goes towards thermal loading of the plasma and not into bulk motions (see Winebarger et al. 2001); second, that the resulting plasma temperatures for all flare energies lie near $`10^7`$ K; and third, that the temperature evolution of the plasma after the flare event is not drastic. These assumptions are supported by the emission measure analysis of several Yohkoh flares by Reale et al. (2001), who find that solar flares are dominated by emission at $`10`$ MK. Hydrodynamic modeling of an ensemble of flaring loops by Güdel et al. (1997) and Güdel (1997) also shows that the bulk of the flare DEM lies above 6 MK; indeed Güdel (1997) finds that the DEM is bimodal around 10 MK, attributable to the slightly smaller temperatures generated by the weaker flares (but which are still significantly hotter than the temperatures achieved by the quiet Sun). Thus, while it is reasonable to expect that flare emission would evolve from temperatures of $`20`$ MK to $`5`$ MK, and that flare events that are less energetic would heat the plasma to a lower temperature, in practice the effects of such variations are very little. Furthermore, because we model the distribution in counts space, i.e., $`\frac{dN}{dc}`$ rather than $`\frac{dN}{dE}`$ directly, we invariably obtain distributions that are shallower than the true distributions121212 Suppose we write the observed counts generated in a detector $`c_{obs}`$ due to emitted energy $`E_{true}`$ as a power law that deviates from linearity by a small amount, $`c_{obs}E_{true}^{1+\delta }`$, $`\delta >0`$. That is, as $`E_{true}`$ decreases, the counts produced in the detector decrease at a faster than linear rate, as would be expected to happen if the lower-energy flares result in lower plasma temperatures, and the detector has a smaller response to these temperatures. This is a reasonable approximation for broad-band instruments such as the EUVE/DS, though it may not be applicable for small-passband detectors such as TRACE (see Aschwanden & Charbonneau 2002). If $`\frac{dN}{dE_{true}}E_{true}^{\alpha _{true}}`$, then $`\frac{dN}{dc_{obs}}c_{obs}^{\alpha _{obs}}`$, where $`\alpha _{obs}=\frac{\alpha _{true}+\delta }{1+\delta }`$, and is always smaller than $`\alpha _{true}`$ for $`\delta >0`$ and $`\alpha _{true}>1`$. Conversely, if $`\delta <0`$, as may happen for extremely high flare energies and temperatures that lie above the best response of the detector, the observed distribution will be steeper than the true distribution, $`\alpha _{obs}>\alpha _{true}`$; however this case has little effect on our analysis results because of the wide temperature response of the EUVE/DS, and any flares that fall outside its range would be too few in number to affect the results.  (see also extensive discussion in Güdel et al. 2002). In the EUVE/DS, changes in temperature of this magnitude causes the observed counts to vary by a factor of $`2`$. Furthermore, the range of flare energies we can model (see Table 4) are much larger than the relatively low-energy “explosive events” that are characterized by large non-thermal velocities (Winebarger et al. 2001 and references therein), and which may not contribute significantly to the plasma heating. That is, we assume that lower energy depositions that may heat the plasma to lower temperatures would not be detectable by the EUVE/DS in any case. Note that the assumed abundances will also affect this to some extent, but its effect is minimal for a broad-band instrument such as is used here. We therefore assume that an observed EUVE/DS count corresponds to a photon of average energy $`1.7\times 10^{11}`$ ergs cm-2 ct-1 over the 59-250 Å range produced by a plasma at $`10^7`$ K (PIMMS v2.5).

##### Cut-offs:

We have adopted upper and lower limits to the extent of the power law distribution (see Equation 5) that are based strictly on numerical expediency: the upper limit is set by the requirement that the largest possible flare can produce no more than the observed number of counts, and the lower limit is defined by the minimum number of counts needed to define an interval. How do these limits compare with theoretical expectations? First note that the solar flare distribution has been explored to much lower energies, $`E10^{23}`$ ergs (Winebarger et al. 2001), than has been achieved for stellar observations. Recently Katsukawa & Tsuneta (2001) have estimated that the probable energy of Parker-type nanoflares is $`<10^{22}`$ ergs. In contrast, we find that if the observed EUVE emission is assumed to originate entirely in flares, then it is sufficient to extend the power law distributions to energies $`E10^{2829}`$ ergs (see Table 4). We suggest that the discrepancy is not due to a higher cut-off energy in the stellar case, but rather that our analysis is physically limited due to a combination of the lack of instrument sensitivity, possible power law changes, and systematic bias due to temperature effects in very small flares (see above). By modeling flares as driven dissipative avalanche systems, Lu et al. (1993) predict a high-energy rollover in the energy distribution whose magnitude depends on the size of the active region where the flares occur. Kucera et al. (1997) have applied this to SMM/HXRBS data and find empirically that the total energy of an event has a maximum, $`E_{cutoff}5\times 10^{28}A_{\mu \mathrm{hs}}^{5/4}`$ ergs, where $`A_{\mu \mathrm{hs}}`$ is the total area of active regions in units of solar micro-hemispheres, i.e., $`1\mu \mathrm{hs}=3.04\times 10^{18}`$ cm2. Thus, for active stars such as the ones we are considering, where active regions covering large fractions of the stellar surface are expected ($`A_{\mu \mathrm{hs}}>>10^4`$), total flare energies in excess of $`10^{33}`$ ergs are achievable, and thus our analyses are valid up to these energies.

##### Background Corrections:

All the stars in our sample are strong sources of EUV emission, and the background is generally small compared to the source strength. We do not subtract the background from the datasets, nor model it separately, but assume that the constant component in the flare model (Equation 3a) includes contamination due to the background. That is, we assume that the background does not vary on timescales smaller than the adopted decay timescale, and that any variations that exist in the instrumental and astrophysical background are small in magnitude and do not contribute to the flare component. Any departures from strict constancy in the background will result in a poorer determination of $`r_C`$ because of the larger spread in the distribution of arrival times, $`g(\delta t)`$ (see Figure 5). In cases where the above assumptions are invalid, the background data will contaminate the signal, but as estimated in Tables 1-3, this contamination will be $`10\%`$.

#### 3.3.3 Advantages and Disadvantages

Unlike previous methods that determine flare number distributions non-parametrically, that is, by directly counting the number of detected flares, we model the flare distribution as a power law, and generate photon arrival-times to compare with the observed event lists. Thus the validity of our results is directly dependent on the applicability of the model (see extensive discussion above in §3.3.2), i.e., the derived parameter values are only as good as the model. Thus, our results are correct subject to the caveat that we require that all flares are physically “similar” and that the flare distribution follows a power law. But the fact that we model the distribution means that we can usefully extend the analysis to the regime where the light curves are dominated by large numbers of small, undetectable flares, and hence this method is best suited to study continuous micro-flaring.

Because the model is stochastic in nature, any realization of the photon arrival times can be significantly different from the observed event list even if the parameter set matches exactly. Thus a large number of simulations must be carried out for each set of parameters, which is a time consuming process. On the other hand, the stochastic nature of the model renders unnecessary an exact match between the model realization and the observed light curve, allowing us to explore weak flare events.

As shown in §3.2 (Figure 5) and §3.3.1, the method loses sensitivity for $`\alpha 3`$ and loses stability for $`\alpha 1.5`$. However, over the range of $`\alpha `$ that is of scientific interest, i.e., spanning the critical value $`\alpha =2`$, the method is robust and can discriminate between the two cases where flare emission sl may be a significant contributor to the coronal emission budget ($`\alpha >2`$), or where the light curve may be dominated by large flares which do not contribute significantly to the energy budget ($`\alpha <2`$).131313Note that the mere fact that $`\alpha >2`$ does not guarantee that the coronal emission is dominated by flaring, but the normalization of the power law distribution must also be large. Our approach allows us to independently determine the flare contribution $`r_F`$ to the total count rate. The results (see §4) show that the normalization is such that the flare contribution is generally $`>50\%`$.

## 4 Results

We have applied the algorithm detailed above (§3) to EUVE/DS data on FK Aqr, V1054 Oph, and AD Leo (see §2). The results are summarized in Table 4. Below we comment on each dataset in detail.

### 4.1 FK Aqr

As anticipated above (see Figure 5), the maximum a posteriori value (MAP; see footnote 10) for FK Aqr are $`\alpha _{_{MAP}}2.5`$, $`r_{F_{_{MAP}}}0.22`$ and $`r_{T_{_{MAP}}}0.36`$, corresponding to a flare contribution of $`65\%`$ to the total emission. The joint probability distribution of $`\alpha `$ and $`r_F`$, marginalized over $`r_T`$, is shown in Figure 7, and the individual probability distributions of $`\alpha `$ and $`r_F`$, marginalized over the other parameters, are shown in Figure 8. The best-fit values are $`\alpha =2.60\pm 0.34`$, $`r_F=0.19\pm 0.12`$, and $`r_T=0.37\pm 0.01`$. The correlations between the parameters are evident in $`p(\alpha ,r_F)`$. Note that there is a small probability ($`0.1`$) that {$`r_F<<r_T`$ & $`\alpha <2`$}. It is worth pointing out the cause of the large values of the probability for small $`r_F`$. First note the presence of a secondary peak at ($`\alpha =2.2`$$`r_F=0.05`$) which contributes to a larger spread in the uncertainty in $`r_F`$, and incidentally also indicates the potential existence of multiple power-law components. Second, this effect is a measure of the stability of the model best-fit parameters to the data; if there are large parts of the parameter space which provide adequate (though not good) fits to the data, their contributions are enhanced due to the larger volume of the space they occupy. This is indeed the case here for combinations of small values of $`r_F`$ and $`\alpha `$, where the skewness induced in $`g_{MODEL}(\delta t)`$ by $`\alpha <2.5`$ is minimized due to the lower weight accorded to it because of the smaller values of $`\kappa `$. Data from Chandra or XMM-Newton, with their higher expected count rates, are necessary to explore the region of smaller values of $`\delta t`$ and thus better constrain the parameter ranges.

### 4.2 V1054 Oph

The light curve of V1054 Oph (Figure 2) shows considerable variability with some relatively weak flare-like events. This is a signature of a flare distribution with large values of $`\alpha `$, and indeed detailed analysis confirms this impression; we find $`\overline{\alpha }=2.74\pm 0.35`$, with an upper bound that is not well constrained (Figure 10). From the joint probability distribution $`p(\alpha ,r_F)`$, the most probable values are $`\alpha _{_{MAP}}2.5`$ and $`r_{F_{_{MAP}}}0.065`$, suggesting that it is likely that almost all the observed emission originates in flare-like events.

### 4.3 AD Leo

In general, the AD Leo data have been analyzed in two separate batches because of the large time intervals and large number of counts involved (which leads to very long computation times) and also because the character of the light curve (see Figure 3) appears to change from being dominated by an intense flare at the beginning (segments II+III) to being steadier, with weaker flare events (segments IV+V). The light-curves also suggest that the flare distribution is characterized by smaller values of $`\alpha `$ than are FK Aqr and V1054 Oph. Detailed analysis confirms the latter, with $`\overline{\alpha _{\mathrm{II}+\mathrm{III}}}=2.17\pm 0.03`$ and $`\overline{\alpha _{\mathrm{IV}+\mathrm{V}}}=2.31\pm 0.11`$ (see Figures 11,12) but the distributions $`p(\alpha )`$ overlap for the two segments, and we cannot rule out at the 10% confidence level that the $`\alpha `$’s are identical for the two datasets. However, the trend for the analyses of individual segments, ($`\overline{\alpha _{\mathrm{II}}}=2.03\pm 0.05`$, $`\overline{\alpha _{\mathrm{III}}}=2.22\pm 0.07`$, $`\overline{\alpha _{\mathrm{IV}}}=2.21\pm 0.06`$, and $`\overline{\alpha _\mathrm{V}}=2.31\pm 0.03`$) does suggest a gradual steepening of the flare distribution when large flares are absent from the dataset.

## 5 Summary

We have modeled the event arrival times from active stars FK Aqr, V1054 Oph, and AD Leo with particular attention to the component that arises from flare-like events. On the Sun, flares are known to be distributed as a power law in energy (Hudson 1991; also see Aschwanden et al. 2000 and references therein), and numerous studies have established that strong stellar flares also follow a power law distribution, with indices ranging from 1.5 - 2.5 (see Audard et al. 2000 and references therein). This is of considerable interest because if the power law index $`\alpha `$ is $`>2`$, then the coronal X-ray losses could in principle be ascribed to weak flare events that are nevertheless numerous enough to dominate the emission.

We consider active dMe stars with known variability in their light curves where numerous flares are seen. Audard et al. (2000) find that in general dF and dG stars tend to have $`\alpha >2`$ while dK and dM stars tend to have lower $`\alpha `$. In particular, they analyze an older dataset of AD Leo (from May 1996) and show that detectable flares have $`\alpha [1.18,2.35]`$ or $`\alpha =2.02\pm 0.28`$ using different methods. A more detailed analysis by Güdel et al. (2002) based on a larger sample of the dataset used here shows $`2.0<\alpha <2.5`$ .

We model the event arrival times using a simple two-component model comprising of a constant rate component and a statistical ensemble of flare components, with the flare energies distributed as a power law. In general, the simplifying assumptions we make (e.g., constancy of decay timescales, ignoring the rise times, assuming a constant counts-to-energy conversion factor, including the background directly in the model, etc.) are conservative, and tend to underestimate the value of $`\alpha `$. We find that all the stars in our sample clearly have $`\alpha >2`$: for AD Leo, $`\alpha `$ lies in the range 2.06 - 2.32; for FK Aqr, $`\alpha =2.60\pm 0.34`$; and for V1054 Oph, $`\alpha =2.74\pm 0.35`$. We thus conclude that coronal heating on these stars is dominated by impulsive energy release events whose energy output is $`23\times 10^{29}`$ ergs, reaching to the microflare range. These results are in contrast to the solar case, where over similar flare energy ranges the observed distribution of flares is shallower, with $`\alpha 1.8`$, i.e., below the critical value of 2.

Further, we directly estimate the contribution of the flare emission to the total observed count rate, as one of the parameters defining the model. Because the energy range over which the model is defined spans over 4 orders of magnitude, and the power law indices indicate steep distributions, we expect that the flare component should contribute significantly to the total emission. Indeed, we find this to be generally $`>50\%`$, and in some cases being $`>80\%`$. Thus, there appears to be no truly “quiescent” emission on some of these low-mass active stars, i.e., emission from apparently stable active region loops as on the Sun is not a dominant component of the observed emission.

We have also explored the possible dependence of the various model parameters on flare energies. The long observation of AD Leo suggests that $`\alpha `$ increases when strong flares are not evident in the data suggesting that the flare distribution steepens for flares of smaller energies, though we cannot rule out the possibility of a statistical fluctuation that mimics this trend. We have also searched for, but do not find, evidence of strong dependence of the decay timescale on flare energies.

We find that if the flare distributions extend to the microflare regime, the energy output due to these weak flares is quite sufficient to account for the entire coronal emission in the EUVE/DS passband. However, it must be noted that the error bars on the parameters derived for the fainter stars FK Aqr and V1054 Oph are quite large (for instance, values of $`\alpha <2`$ cannot be completely ruled out for FK Aqr, and a firm upper bound on $`\alpha `$ cannot be set for V1054 Oph) and it would be of considerable interest to verify and improve these results (and also to extend them to a larger sample of stellar types) using high-quality data such as those obtainable with Chandra and XMM-Newton. Data from these observatories are characterized by good time resolution and in general larger count rates, and will therefore allow us to explore the arrival time difference distribution functions $`g(\delta t)`$ at smaller values of $`\delta t`$, thereby extending the range of values of $`\alpha `$ that the method is sensitive to.

We would like to thank David van Dyk, Alanna Connors, Eric Kolaczyk, and Olivia Johnson for useful discussions. We also thank the referee, J. Linsky, for valuable comments that improved the paper. VLK was supported by NASA grants and the Chandra X-Ray Center during the course of this research. JJD was supported by the Chandra X-Ray Center NASA contract NAS8-39073. The PSI group acknowledges support from the Swiss National Science Foundation (grants 2100-049343 and 2000-058827). This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

## Appendix A Glossary of terms

Here we compile, for reference, all the symbols used in the text (Table 5).


format: text

text:
# Shot noise in generic quantum dots: crossover from ballistic to diffusive transport

## 1 Introduction

The discreteness of the electron charge $`e`$ causes time-dependent current fluctuations even at zero temperature. These fluctuations are known as shot noise and recently have received much attention in mesoscopic physics (for a review see ). One prominent class of mesoscopic systems are classically chaotic ballistic cavities (quantum dots). For incoherent transport through a chaotic quantum dot coupled to two leads, each having $`N`$ channels, the shot noise would assume the Poissonian value $`P_0=2eG_0V`$, where $`G_0=Ne^2/(2h)`$ and $`V`$ is the voltage difference between the leads. For low temperatures, correlations of electrons due to Fermi statistics suppress the noise $`P`$ by a factor of $`=P/P_0`$ relative to the value $`P_0`$ of uncorrelated electrons. In chaotic quantum dots, the suppression factor $`=_{\mathrm{ch}}=1/4`$ has been found to be universal, i.e., independent of the details of systems, by using random matrix theory  and semiclassical methods , which has been confirmed by an experiment . The origin of the low-temperature noise in chaotic quantum dots is the probabilistic nature of quantum transport: Each attempt to transmit a charge $`e`$ through the system succeeds with a finite probability.

Very recently, further reductions of shot noise below $`_{\mathrm{ch}}`$ due to residual signatures of classically deterministic scattering have been discussed by Agam, Aleiner, and Larkin . This non-universal suppression of the noise appears when the dwell time of electrons is smaller than the time scale of quantum diffraction (the Ehrenfest time in chaotic systems). This has also been verified in a recent experiment .

Generic ballistic quantum dots are not fully chaotic but possess a mixed classical phase space , where regular islands are separated from chaotic seas by impenetrable dynamical barriers. Signatures of the mixed phase space in quantum transport have been found in the conductance, which exhibits fractal fluctuations  and isolated resonances . In this proceeding paper, we address the effect of the mixed phase space and of disorder on the shot-noise suppression.

We study shot noise for generic quantum dots coupled to two leads, with varying the positions and the widths of the leads, and allow for an arbitrary strength of elastic diffractive impurity scattering inside the dots. In Ref.  this combination of effects was used to investigate the time scales of deterministic transport through the quantum dots. In the present paper we concentrate on the effect of impurity scattering itself. In absence of impurity scattering, the suppression factor $``$ is systematically reduced below the universal value $`_{\mathrm{ch}}`$ of fully chaotic systems, by an amount which varies with the parameters (positions and widths) of the attached leads. This reduction of $``$ is not only due to the deterministic nature of transport along short chaotic trajectories, but also because quantum diffraction is strongly reduced for transport through regular regions, as well, thus, the dependence of $``$ on the parameters of leads can be understood from the fact that the amount of deterministic processes coupled to the leads varies with the parameters. This feature becomes clarified after adding diffractive impurity scattering with the aid of the Poisson kernel  of random matrix theory. As the rate of diffractive scattering increases, the suppression factor $``$ becomes larger and eventually recovers the universal value $`_{\mathrm{ch}}`$ of fully chaotic systems.

In our numerical simulation, the annular billiard  is chosen as a model of a ballistic quantum cavity with mixed phase space. We will discuss the classical motion in the billiard in Sec. 2, and incorporate diffractive scattering in Sec. 3. In Sec. 4, the dependence of $``$ on the parameters of the leads will be discussed in absence and presence of diffractive scatterings, and the role of regular versus chaotic classical motion on shot noise will be clarified. In Sec. 5, finally, some concluding remarks will be given.

## 2 Annular billiards

An exemplary model for mixed regular and chaotic classical dynamics is the two-dimensional annular billiard , which consists of the region between two circles with radii $`R`$, $`r`$, and eccentricity $`\delta `$. Its schematic diagram is drawn in Fig. 1(a). To study shot noise, two leads (openings) of width $`W`$ are attached opposite to each other at an angle $`\theta [0,\pi /2]`$ with respect to the axis through the two circle centers. Two values of $`W`$ are selected as $`W=0.222R`$ and $`W=0.347R`$; the former corresponds to leads with six channels, while the latter to ten channels. The other parameters of $`r`$ and $`\delta `$ are fixed as $`r=0.6R`$ and $`\delta =0.22R`$, respectively, in terms of $`R`$ throughout this work.

### 2.1 Classical motion in annular billiards

The classical phase space of the annular billiard can be parameterized by the impact parameter $`s`$ and the transverse component of the momentum $`\mathrm{sin}\alpha `$ (with $`\alpha `$ the angle of incidence) of trajectories that are reflected at the exterior circle. The phase space of the closed annular billiard is shown in Fig. 1(b), which displays two regular whispering-gallery regions, a large regular island, neighboring satellite islands, and a chaotic sea. Figs. 1(c) and (d) show the phase space of the open annular billiard (with opening width $`W=0.222R`$), which only includes trajectories of particles that are injected into the billiard through leads attached at $`\theta =0`$ and $`\theta =\pi /2`$, respectively. The large island is well coupled to one of the openings for $`\theta =0`$, and is completely decoupled from both openings at $`\theta =0.5\pi `$ \[see the empty region at the center of Fig. 1(d)\]. The chaotic sea and whispering-gallery regions are well coupled for arbitrary position of the openings, while the coupling of the satellite islands depends on $`\theta `$ and $`W`$. In this way one can select regions in phase space by varying $`\theta `$ and $`W`$.

In the open billiard some classical trajectories entering into the chaotic sea escape without many billiard-boundary scatterings so that the chaotic sea appears less chaotic than in the closed billiard. One can see this modification of phase space, e.g., near $`s=\pi `$ in Fig. 1(c) and around the empty center in (d). The role of these short chaotic trajectories on shot noise will be discussed later.

### 2.2 Conductance and shot noise formula

We calculated the dimensionless conductance $`T`$ and the shot-noise suppression factor $``$ from the relations ,

$`T`$ $`=`$ $`\mathrm{tr}t^{}t,`$ (1)
$``$ $`=`$ $`(2/N)\mathrm{tr}t^{}t(1t^{}t).`$ (2)
Here, $`N`$ is the number of channels in each lead and the transmission matrix $`t`$ is a subblock of the scattering matrix

$$S=\left(\begin{array}{cc}r& t^{}\\ t& r^{}\end{array}\right),$$
(3)
where $`r`$ and $`r^{}`$ are reflection coefficients and the transmission matrix $`t^{}`$ contains the same information as $`t`$. The scattering matrix is obtained numerically by the method of recursive Green functions , for which space is discretized on a square lattice. In terms of the lattice constant $`a`$, we choose $`R=144a`$, $`r=86.4a`$, and $`\delta =31.7a`$. Energy $`E`$ will be measured in units of $`\mathrm{}^2/(2ma^2)`$ and time in units of $`2ma^2/\mathrm{}`$, with $`m`$ the mass of the charge carriers. In these units the mean level spacing $`\mathrm{\Delta }`$ is estimated as $`\mathrm{\Delta }=4\pi /A0.0003`$, where $`A`$ \[$`=\pi (R^2r^2)/a^2`$\] is the dimensionless area of billiard. We will work in the energy window $`E(0.408,0.433)`$. In this window, the Fermi wavelength $`\lambda _\mathrm{F}9.5a`$, resulting in $`N=\mathrm{Int}(2W/\lambda _\mathrm{F})=6`$ for $`W=0.222R=32a`$ and $`N=10`$ for $`W=0.347R=50a`$.

## 3 Elastic diffractive scatterings

We now introduce elastic diffractive scatterings with the aid of the Poisson kernel , a statistical ensemble of random matrix theory . In the Poisson kernel one averages the scattering matrix over an energy range $`E_{\mathrm{av}}`$, in this way eliminating the system-specific details of the dynamics with time scales longer than $`t_{\mathrm{av}}=\mathrm{}/E_{\mathrm{av}}`$, and replaces these by random dynamics of the same universality class as elastic diffractive impurity scattering. The effective mean free scattering time $`t_{\mathrm{av}}`$ can be tuned by changing the energy-averaging window $`E_{\mathrm{av}}`$.

In our numerical simulation, the elimination of the system-specific details on time scales larger than $`t_{\mathrm{av}}`$ is achieved by averaging the scattering matrix $`S`$ over an energy range $`[E_0E_{\mathrm{av}}/2,E_0+E_{\mathrm{av}}/2]`$ of width $`E_{\mathrm{av}}=\mathrm{}/t_{\mathrm{av}}`$ (which will be taken inside the total energy range \[0.408,0.433\] of our numerical simulation),

$$\overline{S}(E_{\mathrm{av}};E_0)=E_{\mathrm{av}}^1_{E_0E_{\mathrm{av}}/2}^{E_0+E_{\mathrm{av}}/2}𝑑ES(E).$$
(4)
Here the information on processes with longer time scales than $`t_{\mathrm{av}}`$ is lost, because it is encoded into the short-range energy correlations (fluctuations) of the scattering matrix, while the information on the dynamics on shorter time scales than $`t_{\mathrm{av}}`$ modulates the scattering matrix on larger energy scales and hence is retained. Diffusive motion then is introduced by coupling to an auxiliary system with scattering matrix $`S_0`$ taken from the appropriate circular ensemble of random matrix theory (observing the same symmetries as the original scattering matrix, as time-reversal or spatial parities ), resulting in

$`S^{}(E_{\mathrm{av}};E_0;S_0)=\overline{S}(E_{\mathrm{av}};E_0)+𝒯^{}(1S_0)^1S_0𝒯.`$ (5)
The ensemble of scattering matrices (5) is the so-called Poisson kernel , with $`\overline{S}`$ the so-called optical scattering matrix. The coupling matrices $`𝒯`$, $`𝒯^{}`$, and $``$ must be chosen such that $`S^{}`$ is a unitary matrix, but the invariance of the circular ensemble guarantees that results do not depend on their specific choice.

Then, the mean suppression factor $`(E_{\mathrm{av}})`$ (or the conductance $`T`$) can be obtained for fixed $`E_{\mathrm{av}}`$ (and hence fixed $`t_{\mathrm{av}}`$), first by averaging the noise (or the conductance) within each Poisson kernel (fixing also $`E_0`$), and then averaging these values over $`E_0(0.408,0.433)`$.

## 4 Results and discussions

### 4.1 Absence of diffractive scatterings: the dependence of shot noise on the positions and on the widths of the leads

The energy-averaged conductance $`T`$ and the energy-averaged suppression factor $``$ in the window of $`E(0.408,0.433)`$ are shown as the values at $`E_{\mathrm{av}}=0`$ of Figs. 2 and 3 for different $`\theta `$ and $`W`$. Our first observation for the shot noise is that for all values of $`\theta `$, the suppression factor is smaller than the universal value $`1/4`$ of fully chaotic motion. Moreover, the shot noise is more suppressed for the cases that one opening couples to the large regular island ($`0<\theta 0.3\pi `$) than for the cases in which the regular islands is decoupled from the openings ($`\theta 0.4\pi `$). This behavior indicates that electrons injected into the large regular island contribute less to the shot noise, which will be clarified in the analysis incorporating diffractive scatterings in Sec. 4.2.

The suppression factors are strongly reduced for $`W=0.347R`$ than for $`W=0.222R`$, as shown in Figs. 2 and 3. This feature will be also understood in Sec. 4.2 from the fact that the amount (per channel in the leads) of deterministic trajectories, which include those in the regular islands and those escaping the dots without many billiard-boundary scatterings in the chaotic sea, increases for the leads with broader widths.

We note that the dependence of conductance $`T`$ on $`\theta `$ shows a similar features to that of the suppression factor; conductance is smaller for the well-coupled cases to the regular islands ($`0<\theta 0.3\pi `$) than for the decoupled cases ($`\theta 0.4\pi `$). This behavior can be understood from the fact that for $`0<\theta 0.3\pi `$ the regular island is well coupled only to one opening, so that electrons entering the regular island do not contribute to transmission. This similarity, however, is accidental. As will be demonstrated in Sec. 4.2, deterministic processes will always suppress shot noise (when their area in phase space is larger than a Planck cell, so that they are well resolved), while they enhance or suppress conductance, depending on the shape of the quantum dot and the positions of the leads.

### 4.2 Role of regular versus chaotic classical motion in shot-noise suppression

The effect of the diffractive scattering on the conductance $`T`$ and the shot-noise suppression factor $``$ for different positions of the leads are shown in Figs. 2 and 3.

The values of $`T`$ and $``$ at $`E_{\mathrm{av}}=0`$ are those in absence of diffactive impurity scattering, because the Poisson kernel becomes a delta function at $`S(E_0)`$ as $`\overline{S}`$ approaches this unitary matrix for $`E_{\mathrm{av}}0`$. For increasing $`E_{\mathrm{av}}`$, $`T(E_{\mathrm{av}})`$ and $`(E_{\mathrm{av}})`$ approach their universal values of random-matrix theory, because $`\overline{S}=0`$ in this limit. Note that spatial symmetries such as up-down and left-right reflections affect the weak localization corrections to the universal values . The correction is enhanced (absent) for $`\theta =0`$ ($`0.5\pi `$), where there is up-down (left-right) reflection symmetry , compared with the other values of $`\theta `$ with no spatial symmetry.

The function $`(E_{\mathrm{av}})`$ is monotonically increasing almost everywhere (negative values of the slope are within the statistical uncertainty in the numerical simulation), indicating that the noise is enhanced when the system-specific deterministic details of the transport are replaced by the indeterministic transport of random matrix theory. Consequentially, one can define the probability distribution of deterministic processes with dwell time $`t_{\mathrm{dwell}}t_{\mathrm{av}}`$ by the rate of change of the shot-noise suppression $`(E_{\mathrm{av}})`$ which can be substantiated by a semiclassical analysis .

The crossover of shot noise from ballistic to diffusive transport occurs at different scattering rate $`\mathrm{}/E_{\mathrm{av}}`$, depending on $`\theta `$ and $`W`$. We will discuss first for the case of $`W=0.222R`$. In the three cases $`\theta =0`$, $`0.1\pi `$, $`0.2\pi `$, in which the leads couple to the large regular island in phase space, the crossover occurs already at very small $`E_{\mathrm{av}}`$, while it starts above an energy $`E_{\mathrm{indet}}0.002`$ for $`\theta =0.4\pi `$, $`0.5\pi `$ (where the regular region is decoupled). The former is due to the deterministic processes of very long dwell times in the large regular island, which can be affected by diffractive scatterings of very small scattering rate. On the other hand, the latter is due to the fact that the processes in the chaotic sea are already indeterministic for times larger than $`\mathrm{}/E_{\mathrm{indet}}`$, which can be interpreted as the Ehrenfest time.

For the leads with larger width $`W`$, more deterministic processes of regular regions can couple to the leads and more short deterministic trajectories in the chaotic sea appear. This feature entails that for the leads with $`W=0.347R`$ (Fig. 3), $`(E_{\mathrm{av}})`$ approaches to $`_{\mathrm{ch}}`$ more slowly than for the case of $`W=0.222R`$ (Fig. 2). Similarily, for the leads at $`\theta =0.5\pi `$ and with $`W=0.347R`$ the suppression factor $`(E_{\mathrm{av}})`$ starts to increases at earlier $`E_{\mathrm{av}}`$ than $`E_{\mathrm{indet}}`$ (which is obtained in the curves for $`W=0.222R`$), since small regular islands additionally couple to the leads with $`W=0.347R`$.

We note that in contrast to $`(E_{\mathrm{av}})`$, the sign of the slope of $`T(E_{\mathrm{av}})`$ depends on the system-specific geometry, e.g., on lead position $`\theta `$ \[see Figs. 2(a) and 3(a)\]. This shows that conductance can be either enhanced or suppressed by diffractive scattering (indeterministic processes).

## 5 Concluding remarks

The shot noise reduction due to deterministic transport suggests that shot noise can be used as a probe of classical phase space of generic quantum dots (conductance can not be used as such a probe). The crossover of shot noise from ballistic to diffusive transport, shown in Figs. 2(b) and 3(b), could be probed experimentally by measuring shot noise while tuning the disorder strength, which can be achieved by adjusting a gate voltage .

## References


text:
# Jet Quenching with Detailed Balance in a Thermal QCD Medium This work was supported by the NSFC under projects Nos. 19928511 and 10135030, and by the U.S. Department of Energy under Contract No. DE-AC03-76SF00098.

## 1 INTRODUCTION

Jet quenching in high-energy heavy-ion collisions suppresses large transverse momentum hadrons as compared to $`pp`$ collisions at the same energy . It provides a direct probe of the dense matter produced in high-energy heavy-ion collisions. Recent theoretical studies  of jet quenching have concentrated on gluon radiation induced by multiple scattering in medium. Because of the presence of thermal gluons in the hot medium the stimulated gluon emission and absorption by the propagating parton should also be taken into account consistently. Such a detailed balance is crucial for parton thermalization and should also be important for calculating the energy loss of an energetic parton in quark-gluon plasma (QGP). We will show that the partial cancellation between stimulated emission and thermal absorption results in a net reduction of parton energy loss induced by multiple rescattering. Though such a reduction is negligible for asymptotically large parton energy, it is still important for intermediate values of energy.

## 2 FINAL-STATE AND RESCATTERING-INDUCED ABSORPTION

In an axial gauge and in the leading-log approximation, the final-state radiation amplitude off a quark can be factorized from the hard scattering. Taking into account of both stimulated emission and thermal absorption in a thermal medium with finite temperature $`T`$, the probability of gluon radiation with energy $`\omega `$ can be written as

$$\frac{dP^{(0)}}{d\omega }=\frac{\alpha _sC_F}{2\pi }\frac{dz}{z}\frac{d𝐤_{}^2}{𝐤_{}^2}\left[N_g(zE)\delta (\omega +zE)+\left(1+N_g(zE)\right)\delta (\omega zE)\theta (1z)\right]P(\frac{\omega }{E}),$$
(1)
where $`N_g(|𝐤|)=1/[\mathrm{exp}(|𝐤|/T)1]`$ is the thermal gluon distribution, $`C_F`$ is the Casimir of the quark jet in the fundamental representation and the splitting function $`P_{gq}(z)P(z)/z=[1+(1z)^2]/z`$ for $`qgq`$. In Eq.(1) the first term is from thermal absorption and the second term from gluon emission with the Bose-Einstein enhancement factor.

To define the effective parton energy loss, we only consider gluon radiation outside a cone with $`|𝐤_{}|>\mu `$. Here $`\mu `$ should be a nonperturbative scale for parton fragmentation in the vacuum. In a hot QGP medium, it can be replaced by the gluon’s Debye screening mass. Assuming the scale of the hard scattering as $`Q^2=4E^2`$, one has the kinematic limits of the gluon’s transverse momentum, $`\mu ^2𝐤_{max}^24|\omega |(E\omega ).`$

Subtracting the gluon radiation spectrum in the vacuum, one then obtains the energy loss due to final-state thermal absorption and stimulated emission. For asymptotically large parton energy, $`ET`$, we have

$`{\displaystyle \frac{\mathrm{\Delta }E_{abs}^{(0)}}{E}}`$ $`=`$ $`{\displaystyle 𝑑\omega \frac{\omega }{E}\left(\frac{dP^{(0)}}{d\omega }\frac{dP^{(0)}}{d\omega }|_{T=0}\right)}`$ (2)
$``$ $`{\displaystyle \frac{\pi \alpha _sC_F}{3}}{\displaystyle \frac{T^2}{E^2}}\left[\mathrm{ln}{\displaystyle \frac{4ET}{\mu ^2}}+2\gamma _\mathrm{E}+{\displaystyle \frac{6\zeta ^{}(2)}{\pi ^2}}\right],`$
where $`\gamma _\mathrm{E}0.5772`$ and $`\zeta ^{}(2)0.9376`$.

During the propagation of the hard parton after its production, it will suffer multiple scattering with targets in the medium. Recently many theoretical studies  have investigated the rescattering-induced radiative energy loss. Here we will focus on the stimulated emission and thermal absorption associated with rescattering in a hot QGP medium.

Assume a hard parton is produced at $`𝐲_0=(y_0,𝐲_0)`$ inside the medium with $`y_0`$ being the longitudinal coordinate. The interaction between the jet and target partons is modeled by a static color-screened Yukawa potential which is initially proposed by Gyulassy-Wang (GW) ,

$$V_n=2\pi \delta (q^0)v(𝐪_𝐧)e^{i𝐪_n𝐲_n}T_{a_n}(j)T_{a_n}(n),v(𝐪_n)=\frac{4\pi \alpha _s}{𝐪_n^2+\mu ^2}.$$
(3)
Here $`𝐪_n`$ is the momentum transfer from a target parton $`n`$ at $`𝐲_n=(y_n,𝐲_n)`$, $`T_{a_n}(j)`$ and $`T_{a_n}(n)`$ are the color matrices for the jet and target parton.

Following the framework of opacity expansion developed by Gyulassy, Lévai and Vitev (GLV) and Wiedemann , we will only consider the contributions to the first order in the opacity expansion. It was shown by GLV that the higher order corrections contribute little to the radiative energy loss. Including stimulated emission and thermal absorption we can express the radiation probability at the first order in opacity as

$`{\displaystyle \frac{dP^{(1)}}{d\omega }}`$ $`=`$ $`{\displaystyle \frac{C_2}{8\pi d_Ad_R}}{\displaystyle \frac{dz}{z}\frac{d^2𝐤_{}}{(2\pi )^2}\frac{d^2𝐪_{}}{(2\pi )^2}v^2(𝐪_{})\frac{N}{A_{}}Tr\left[|R^{(S)}|^2+2Re\left(R^{(0)}R^{(D)}\right)\right]}`$ (4)
$`\times P({\displaystyle \frac{\omega }{E}})\left[\left(1+N_g(zE)\right)\delta (\omega zE)\theta (1z)+N_g(zE)\delta (\omega +zE)\right],`$
where $`C_2`$ and $`C_A`$ are Casimirs of the target parton in the fundamental representations in $`d_R`$ and $`d_A`$ dimension. Averaging over the longitudinal target profile is defined as $`\mathrm{}=𝑑y\rho (y)\mathrm{}`$ with $`\rho (y)=2\mathrm{exp}(2y/L)/L`$. We denote $`R^{(0)}`$, $`R^{(S)}`$ and $`R^{(D)}`$ as the radiation amplitude associated with zero, single and double Born rescattering with parton targets in the QGP, respectively. They can be obtained as

$`R^{(0)}`$ $`=`$ $`2igT_c{\displaystyle \frac{𝐤_{}ϵ_{}}{𝐤_{}^2}},`$ (5)
$`R^{(S)}`$ $`=`$ $`2ig\left(𝐇T_aT_c+𝐁_1e^{i\omega _0y_{10}}[T_c,T_a]+𝐂_1e^{i(\omega _0\omega _1)y_{10}}[T_c,T_a]\right)ϵ_{},`$ (6)
$`R^{(D)}`$ $`=`$ $`2igT_ce^{i\omega _0y_{10}}\left({\displaystyle \frac{C_F+C_A}{2}}𝐇e^{i\omega _0y_{10}}+{\displaystyle \frac{C_A}{2}}𝐁_1+{\displaystyle \frac{C_A}{2}}𝐂_1e^{i\omega _1y_{10}}\right)ϵ_{}.`$ (7)
Here $`ϵ_{}`$ is the transverse polarization of the radiative gluon, $`y_{10}=y_1y_0`$,

$`\omega _0`$ $`=`$ $`{\displaystyle \frac{𝐤_{}^2}{2\omega }},\omega _1={\displaystyle \frac{(𝐤_{}𝐪_{})^2}{2\omega }},`$ (8)
$`𝐇`$ $`=`$ $`{\displaystyle \frac{𝐤_{}}{𝐤_{}^2}},𝐂_1={\displaystyle \frac{𝐤_{}𝐪_{}}{(𝐤_{}𝐪_{})^2}},𝐁_1=𝐇𝐂_1.`$ (9)
The zero-temperature part of the parton energy loss corresponds to the radiation induced by rescattering which should be the same as obtained in the previous studies . We denote this part as $`\mathrm{\Delta }E_{rad}^{(1)}`$. The remainder or temperature-dependent part of energy loss induced by rescattering at the first order in opacity is

$`\mathrm{\Delta }E_{abs}^{(1)}`$ $`=`$ $`{\displaystyle 𝑑\omega \omega \left(\frac{dP^{(1)}}{d\omega }\frac{dP^{(1)}}{d\omega }|_{T=0}\right)}`$ (10)
$`=`$ $`{\displaystyle \frac{\alpha _sC_F}{\pi }}{\displaystyle \frac{L}{\lambda _g}}E{\displaystyle 𝑑z\frac{d𝐤_{}^2}{𝐤_{}^2}d^2𝐪_{}|\overline{v}(𝐪_{})|^2\frac{𝐤_{}𝐪_{}}{\left(𝐤_{}𝐪_{}\right)^2}N_g(zE)}`$
$`\times \left[P(z)Re(1e^{i\omega _1y_{10}})+P(z)Re(1e^{i\omega _1y_{10}})\theta (1z)\right],`$
where $`\lambda _g=C_F\lambda /C_A`$ is the mean-free-path of the gluon, the factor $`1\mathrm{exp}(i\omega _1y_{10})`$ reflects the destructive interference arising from the non-Abelian Landau-Pomeranchuck-Migdal (LPM)  effect, and $`|\overline{v}(𝐪_{})|^2`$ is the normalized distribution of momentum transfer from the scattering centers,

$$|\overline{v}(𝐪_{})|^2=\frac{1}{\pi }\frac{\mu _{eff}^2}{(𝐪_{}^2+\mu ^2)^2},\frac{1}{\mu _{eff}^2}=\frac{1}{\mu ^2}\frac{1}{q_{max}^2+\mu ^2},q_{max}^23E\mu .$$
(11)
In the limit of $`EL1`$ and $`E\mu `$, One can then get the approximate asymptotic behavior of the energy loss,

$`{\displaystyle \frac{\mathrm{\Delta }E_{rad}^{(1)}}{E}}`$ $``$ $`{\displaystyle \frac{\alpha _sC_F\mu ^2L^2}{4\lambda _gE}}\left[\mathrm{ln}{\displaystyle \frac{2E}{\mu ^2L}}0.048\right],`$ (12)
$`{\displaystyle \frac{\mathrm{\Delta }E_{abs}^{(1)}}{E}}`$ $``$ $`{\displaystyle \frac{\pi \alpha _sC_F}{3}}{\displaystyle \frac{LT^2}{\lambda _gE^2}}\left[\mathrm{ln}{\displaystyle \frac{\mu ^2L}{T}}1+\gamma _\mathrm{E}{\displaystyle \frac{6\zeta ^{}(2)}{\pi ^2}}\right].`$ (13)
Similarly to the final-state absorption, the minus sign in Eq. (13) indicates that the contribution from thermal absorption associated with rescattering is larger than that of stimulated emission, resulting in a net energy gain.

## 3 NUMERICAL RESULTS AND CONCLUSION

The significance of the thermal absorption relative to the induced radiation can be studied by evaluating Eqs. (2) and (10) numerically. From the perturbative QCD at finite temperature the Debye screening mass can be taken as $`\mu ^2=4\pi \alpha _sT^2`$ . In the GW model the mean-free-path for a gluon $`\lambda _g`$ is

$$\lambda _g^1=\sigma _{qg}\rho _q+\sigma _{gg}\rho _g\frac{2\pi \alpha _s^2}{\mu ^2}9\times 7\zeta (3)\frac{T^3}{\pi ^2},$$
(14)
where $`\zeta (3)1.202`$. With fixed values of $`L/\lambda _g`$ and $`\alpha _s`$, $`\mathrm{\Delta }E/\mu `$ is a function of $`E/\mu `$ only.

The ratios of the calculated radiative energy loss with and without stimulated emission and thermal absorption as functions of $`E/\mu `$ are illustrated in Fig. 1 for $`L/\lambda _g=3`$,5 and $`\alpha _s=0.3`$. The thermal absorption reduces the effective parton energy loss by about 30-10% for intermediate values of parton energy. The energy dependence of the effective parton energy loss becomes stronger. This will affect the shape of suppression of $`p_T`$ hadrons spectrum in the intermediate energy region. For partons with very high energy the effect of the gluon absorption is small and can be neglected.


text:
# 1 The zero temperature phase diagram of the dual [𝑈(1)]² Higgs model in the Bogomol’ny limit in the presence of the external electromagnetic field. The confinement and deconfinement phases are denoted as 𝐶 and 𝐷. The asymmetric confinement phases are denoted as 𝐴_𝑖 and 𝐴_{𝑗𝑘}, where subscripts indicate the components of the monopole field Φ_𝑖 which are condensed in these phases.

## Acknowledgments

The author is grateful to V.G. Bornyakov and Y. Koma for interesting discussions. The work is supported by the JSPS Fellowship P01023.