text:
# Close frequency pairs in Delta Scuti stars

## 1 Introduction

The $`\delta `$ Scuti stars are variable stars situated in the Classical Instability Strip, which pulsate with a large number of simultaneously excited radial and nonradial modes, which makes them well-suited for asteroseismological studies. The amplitudes of the more dominant modes in the typical $`\delta `$ Scuti star are a few millimag, which is much higher than found in the Sun. It is now possible for ground-based telescopes to detect a large number of simultaneously excited modes with millimag amplitudes in stars other than the Sun. These studies require hundreds of telescope hours at observatories spread around the world.

If it can be established that $`\delta `$ Scuti stars contain close frequencies of nonradial pulsation as close as 0.05 cd-1, this has important asteroseismological implications. Such close frequencies cannot be explained by first-order rotational splitting since almost all $`\delta `$ Scuti stars rotate faster than 10 kms-1. The existence of extreme splitting asymmetries, which may be produced in rapidly rotating stars (e. g., Goupil et al. 2000), may produce a close doublet in the power spectrum: the observational evidence for such theoretically predicted behavior is still missing. Similar frequencies are also predicted in the asymptotic case if the two modes belong to different nonradial degrees with $`\mathrm{}`$ values separated by 2. In fact, the small frequency separations are used as an important asteroseismological tool. Furthermore, it is presently debated whether the large number of closely spaced modes of mixed p and g character predicted by theory are indeed found in evolved stars (e.g., see Breger & Pamyatnykh 2002). Close modes might provide some evidence for such modes.

It therefore appears prudent to examine whether close frequencies of pulsation exist in $`\delta `$ Scuti stars and if they do, to look at their properties in order to determine their origin.

It is important to emphasize that the presence of double peaks in the power spectrum of a pulsator, or the emergence of a second peak after prewhitening the first peak at a similar frequency, do not necessarily mean that there exist two pulsation modes with close frequencies. A number of different explanations for such a behavior need to be tested. These range from instrumental effects, methods to find the frequencies of pulsation, to other effects intrinsic to the star. Before one can conclude the presence of closely spaced double modes, the possible problems associated with the available data sets need to be investigated in detail, in particular:

(i) The problem of frequency resolution: A crucial requirement for the discovery of close frequencies is sufficient frequency resolution. Loumos & Deeming (1978) derive a frequency resolution of 1.5/$`\mathrm{\Delta }`$T, where $`\mathrm{\Delta }`$T is the total span of observations. They make the point that the generally used figure of 1/$`\mathrm{\Delta }`$T is not really correct from a theoretical point of view. Furthermore, data are not equidistantly spaced and longer time gaps can occur. This can decrease the resolution. In an extreme case, the resolution of two months of observing spaced one year apart corresponds to that of a single month. Poor frequency resolution, alone, does not lead to the detection of a double mode. Nevertheless, the combination of poor frequency resolution with one of the problems listed below, can cause severe problems in correctly interpreting the power spectrum.

(ii) The problem of incorrect prewhitening: In the analysis to deduce the multiple excited modes present in photometric data, prewhitening detected frequencies is a common tool to detect further frequencies. Prewhitening with an incorrect frequency value leaves a spurious peak in the power spectrum with a frequency close to the prewhitened frequency, i. e., a doublet is found which does not correspond to two separate modes in the star. Consequently, it is very important to examine every detected frequency doublet whether a single frequency with a better value can explain the doublet.

(iii) Artifacts of amplitude variability of a single mode: Amplitude variability on the time scale of years (or even months) is a common feature in $`\delta `$ Scuti stars (e.g., see Breger 2000a, Arentoft et al., 2001). Some $`\delta `$ Scuti stars are known to exhibit very strong amplitude variability, e. g., the star 4 CVn with an average annual amplitude variability of 12$`\%`$ for the different pulsation modes and 40$`\%`$ over a decade (Breger 2000a). The techniques used to detect the frequencies of pulsation usually rely on constant amplitudes. The amplitude variability shows up in the power spectrum as multiple peaks. Strictly speaking, amplitude variability associated with a single frequency does not lead to a double peak, but to a more complex structure. However, with amplitude variability the exact value of the pulsation frequency can sometimes not be determined, so that the problem of incorrect prewhitening also occurs. In some cases, with amplitude variability only a second peak with sufficient power is noticed, while the other features blend into the noise present in the power spectrum.

(iv) Observational problems: These range from small systematic time errors in the data from one of the multiple observing sites to different effective wavelengths for nominally similar filters at different sites. The latter is caused by the fact that the amplitudes of $`\delta `$ Scuti stars are extremely wavelength dependent. These problems should be carefully checked for by the observers during data reduction, since at a later stage it is difficult to discover such effects.

(v) Effects of aliasing: In a multiple-mode star, aliasing may result in combinations of peaks with incorrect frequencies. Prewhitening such peaks, again, leaves artifacts in the power spectrum, which can be interpreted as double modes. For small-amplitude modes, we regard this effect as one of the most dangerous possibilities. This needs to be examined carefully with statistical programs specializing in the discrimination between various multifrequency solutions, such as PERIOD (Breger 1990) and PERIOD98 (Sperl 1998).

## 2 $`\delta `$ Scuti stars with reported close frequencies

A number of $`\delta `$ Scuti stars have been reported to have close frequency pairs with separations less than 0.1 c/d. In order to estimate whether this is a common occurrence, we have examined the literature for all stars listed in the $`\delta `$ Scuti star catalogue (Rodríguez et al. 2000). For most of these stars, very little information is available, so that the question of close frequencies cannot be addressed. Consequently, we selected those stars which were well-studied so that close frequencies (if they exist) could be found. Due to the inhomogeneity of the data, no rigorous, consistent definition is possible. Therefore we selected stars with more than 200 hours of photometric observations, annual observing runs of 20 days or longer, as well as no obvious problems of analysis, measurements, and the distribution of the observations in time. 17 stars could be selected. Of these, seven stars had reported close frequencies. This underestimates the true occurrence because many $`\delta `$ Scuti stars have very small amplitudes so that only a few modes are observed. If we further restrict the sample to those stars for which more than ten frequencies of pulsation had been found, every one of the seven stars had at least one reported close frequency pair with separations of less than 0.06 c/d. This is surprising, since the frequencies typically cover a wide range between 5 and 15 c/d. Due to inhomogeneity of the data and the small number of well-studied stars, one should be careful not to overinterpret these numbers. A safe conclusion would be that the majority of well-studied $`\delta `$ Scuti stars show close frequency pairs.

Table 1 lists the seven stars with reported close frequency pairs with separations of 0.06 c/d or less.

Let us mention two of these stars: Alvarez et al. (1998) reported two close frequency pairs for the star BW Cnc in the Praesepe clusters. The pairs at 138.7/139.1 and 260.6/261.5 $`\mu `$Hz have separations of 0.4 and 0.9 $`\mu `$Hz, respectively. This translates to 0.035 and 0.078 cd-1, respectively. The length of the campaign of 26d suggests a frequency resolution of 0.058 cd-1, according to the Deeming & Loumos criterion so that at least one of the reported doublets has sufficient frequency resolution. Another star is BI CMi, which was studied extensively by the Delta Scuti Network for a number of years so that a frequency separation as small as 0.005 cd-1 could still be resolved.

Some of these studies are quite extensive. Due to the large amounts of data it appears inappropriate to dismiss all reports of close frequencies in $`\delta `$ Scuti stars as due to observational difficulties. It is important to examine a star in close detail to test the different interpretations.

We have chosen the star BI CMi for these detailed examinations. The choice was motivated by the fact that the star has extensive observations available and that we are intimately familiar with the data, since they were obtained by the Delta Scuti Network located in Vienna.

After the pioneering study of BI CMi by Mantegazza & Poretti (1994), from 1997 to 2000, the Delta Scuti Network obtained more than 1000 hours of photometric data and deduced more than 25 frequencies of pulsation (Breger et al. 2002). The 1998/1999 and 1999/2000 observing seasons covered 125d and 109d, respectively. This provides a frequency resolution of 0.003 cd-1 for the two years. One of the important and puzzling results of the study concerns the presence of five features in the power spectrum which can not be explained as single pulsation modes with constant amplitudes. Their unusual nature modes was discovered when prewhitening the frequency of each mode revealed an additional peak in the power spectrum with a frequency very close to that of the prewhitened mode and/or when different amplitudes were found for different years. We note that the detection of these features is statistically very significant (amplitude signal/noise ratio as high as 6), so that their reality will not be questioned here. The frequency spectrum of BI CMi is shown in Fig. 1, where the location of these features is also shown.

Before the nature of these features is examined in details, one should eliminate the possibility of their origin due to observational uncertainties, complex aliasing effects from other modes or incorrect prewhitening. While such problems can never be eliminated completely, except for feature 5 (11.26 cd-1) we consider them extremely unlikely as an explanation for the following reason: The data were subdivided into various subgroups based on year and observatory and a large number of new multifrequency solutions with a variety of hypothetical scenarios computed. The result was always the same: two modes with close frequencies or amplitude variability. In the case of feature 5, we note the existence of a pulsation mode at 10.24 cd-1 with a similar amplitude, so that interaction through 1 cd-1 aliasing presents a potential influence. Again, no obvious alternate explanation could be derived, but for this feature we would like to advise some caution.

Two explanations remain: pairs of close frequencies and amplitude variability. The next section will concentrate on distinguishing between these possibilities.

## 3 Close frequency pairs or just amplitude variability?

In principle, the choice between the hypotheses looks simple: one needs to consider both models and calculate optimum fits of the relevant parameters to the data. The models with the lowest residuals between fit and observations would then be chosen. This valid method has one serious drawback: in a multiperiodic star, the two models may lead to only slightly different residuals. Furthermore, the solution with the lowest residuals may not be the correct one due to complex aliasing effects (see above) and the different number of free parameters fitted.

It appears prudent to adopt a more classical, analytical approach. The two hypotheses both lead to a variable amplitude, but to very different phasing behavior: two beating frequencies appear as a single frequency with variable phase. The most extreme case occurs for two frequencies with the same amplitudes: there is a half cycle phase jump at every new cycle of the beat period. Even if the amplitudes are not identical, the phase jump still exists, but is progressively smeared out in time as the amplitudes become more and more unequal. A numerical modeling of the amplitude ratios and phase shifts can be seen in Breger (1981), where the two-frequency model could be ruled out for the cepheid HR 7308.

To distinguish between the two hypotheses, we will examine the amplitude and phase variability for each mode of BI CMi and compare the one- and two-frequency models in the time, amplitude and the time, phase domains. All computations of amplitudes and phases of the mode to be tested for duplicity were carried out simultaneously with all the other frequencies or frequency pairs known for BI CMi as listed in Breger et al. (2002). The simultaneous calculation eliminates any bias caused by spectral leakage from other modes which were assumed to have constant annual amplitudes. Numerical simulations showed, however, that such bias is very minor to start with. The available data have been obtained through the Stromgren $`y`$ and $`v`$ filters. We have used both filters to reduce the noise. For the analysis, the $`v`$ amplitudes were multiplied by an experimentally determined factor of 0.632 to agree with the $`y`$ amplitudes and any small shifts between the times of maximum for the two filters ignored. This is a safe procedure because in the test we are looking for phase jumps of up to half a cycle. Furthermore, the results are insensitive to the chosen value for the amplitude ratio because of the similar coverage within each for year for the two colors. The data, which were obtained at different observatories with different detectors and filters, were examined for systematic differences (which might lead to errors in the power spectra). No systematic effects or problems were found.

The test relies mainly on the phase shifts of the assumed single frequency. The results can be represented in different graphical forms based on the uncertainties of the calculated points and the number of measured data points required for each phase range considered. Consequently, we have chosen slightly different graphical forms for demonstrating the results for the different frequencies and to strengthen the conclusions.

Table 2 lists the four best candidates for double modes in BI CMi and shows the relevant parameters for both the two-frequency and single-frequency hypotheses.

### 3.1 The mode(s) near 10.43 cd-1

The pulsation mode near 10.43 cd-1 shows strong amplitude variability on the timescale of about 132d. To apply the phase test we have assumed the single (best) frequency of 10.43142 cd-1 and computed separate solutions of the amplitudes and phases in 10 to 15 d intervals. This could be done for all three observing seasons. The separate results were then plotted together utilizing a beat frequency of 132d. The results are shown in Fig. 2.

We consider the results as surprisingly unambiguous and find:

* The rapid phase change is found. The fact that the phase is almost constant outside the times of the rapid switch indicates that we have chosen a correct trial period. The rapid phase change is incompatible with the single frequency/variable amplitude hypothesis and a beating phenomenon of two close periods with similar amplitudes may be involved.
* The fast phase changes occur near the times of minimum amplitude. This is another requirement of the two-frequency hypothesis.
* Both the 1998/1999 and 1999/2000 observing seasons show exactly the same behavior. Apart from confirming our conclusions in favor of the two-frequency hypothesis this result indicates that the amplitudes of the two pulsation modes are not strongly variable on an annual basis.
* The 7 nights of data from 1997 spanning 15 days can provide only one point in Fig. 2. This point is in agreement with the other data and supports the two-frequency model.
* We have computed the best values of the two frequencies, amplitudes and phases by a least-squares fit. In fact, the two-frequency model successfully predicts both the observed amplitudes and phases for all three observing seasons.

We conclude that the model with a single frequency and variable amplitudes can be ruled out from the observations and that the model with two close frequencies provides excellent fits to the data.

### 3.2 The mode(s) near 8.65 cd-1

The successful analysis for the previous (double) mode was repeated for the 8.65 cd-1 mode. Here, the shorter beat cycle of $``$ 56d makes the test more difficult than in the previous example, because the time bins for the examination have to be shorter in length and thereby contain fewer points. Fig. 3 lists the important phase information for the two extensive 1998/1999 and 1999/2000 observing seasons. The data were plotted on the same diagram with different symbols by adopting an experimentally determined best shift of 335d (6 beat cycles).

As for the previously examined mode, there exist regular phase jumps of $``$ half a cycle. The beat cycle is short enough for each observing season to contain two of these phase discontinuities. These jumps occur at minimum amplitude.

We conclude that the mode near 8.65 cd-1 is actually double with two close frequencies and that the phase jumps are incompatible with the single mode/variable amplitude hypothesis. As was previously also found for the 10.43 cd-1 pair, the amplitudes associated with the two modes are constant or nearly constant on an annual basis.

### 3.3 The mode(s) near 4.8 cd-1

The one or two modes near 4.8 cd-1 are even more difficult to analyze because of the short beat cycle of 28.3 d and the small amplitudes near 1 mmag. Consequently, it was no longer feasible to examine different beat cycles separately. The 1998/1999 and 1999/2000 data were combined. The short data set from 1997 was omitted in order to avoid potential errors accumulated from the long time base.

To apply the phase test we have assumed the single (best) frequency of 4.7826 cd-1 and computed separate solutions of the amplitudes and phases for each for the 20 phase bins.

The results are shown in Fig. 4. We can see systematic shifts in the phasing of the assumed single frequency. These shifts correspond to the shifts calculated from the two-frequency hypothesis. Because of the different amplitudes associated with the two modes, a sudden phase jump is not expected in this case.

We conclude that in BI CMi there exist two close frequencies near 4.8 cd-1.

### 3.4 The mode(s) near 12.35 cd-1

The Fourier analyses and multiple least squares solutions show two frequencies at 12.3499 and 12.3268 cd-1. This leads to a beat period of 43d. As before, we have divided the data into 17 bins (each with $``$ 1000 data points) grouped according to the phasing of the beat period. However, because of the short beat period and the low amplitudes, we were not able find an optimum single frequency to apply the test to. Statistically, the best result always turned out to be one of the frequencies of the frequency pair with a variable amplitude, as would occur by the effect of a neighboring mode or true amplitude variability. We have chosen the frequency with the best residuals, viz., 13.50 cd-1. Fig. 5 shows that the phasing test is not conclusive for this case. Both the size of the residuals of the multiperiodic solution and the power spectra favor the two-frequency hypothesis over the single-frequency hypothesis. Nevertheless, a convincing proof is lacking for this frequency pair.

### 3.5 The mode(s) near 11.26 cd-1

The mode at 11.2608 cd-1 shows strong amplitude variability from year to year. During the 1998/99 observing season, covering 120d, amplitude variability including a minimum was found. On the other hand, no significant variability was found during the 110d of the 1999/2000 observing season. These amplitudes, assuming a single frequency, are shown in Table 3, where the formal uncertainties in the amplitudes and phasing are also shown. As before, the formulae for the computation of these uncertainties assume random errors. The true uncertainties may therefore be somewhat larger.

The amplitude variability indicates a long beat period. Extensive two-frequency modelling shows a beat period between 500 and 600d, or $`\mathrm{\Delta }`$f $``$ 0.002 cd-1 between the two (hypothetical) modes. The modelling showed that the two-frequency hypothesis is indeed a strong possibility, but that a unique, reliable solution with two close frequencies could not be found. Including the 1991 data by Mantegazza & Poretti ($``$ 2 mmag amplitude) did not improve the situation because of the long observing gap between 1991 and 1997.

The data in Table 3 shows why the phasing test of the two-frequency hypothesis has to be unreliable: the test relies on the phase change near minimum amplitude. This phase change is detected, viz. 0.12 $`\pm `$ 0.07 cycles, but is not statistically significant in the light of the above discussion on the uncertainties.

We conclude that due to the long beat cycle of the 11.26 cd-1 mode, we cannot distinguish between the variable-amplitude and two-frequency hypotheses.

## 4 The strange case of 4 CVn

We have seen that one of the signatures of two close modes is the phase shift of up to half a cycle near minimum amplitude. Such a phase shift has been seen before in another $`\delta `$ Scuti star: the 7.37 cd-1 mode of 4 CVn. This mode exhibits the strongest long-term amplitude variability seen so far. Details can be found in Breger (2000a, b: Figs. 11 and 12). The amplitude dropped from 15 mmag in 1974 to 4 mmag in 1976 and 1 mmag in 1977. In subsequent years, the amplitude has been increasing again. The amplitude variability, by itself, is not unusual: the fascinating change is the phase shift of 0.48 $`\pm `$ 0.02 cycles between 1976 and 1977, i.e., at minimum amplitude. The determination of accurate phases became possible when the extensive 1996/7 data of 4 CVn led to accurate frequency values, which also fit the earlier data well.

Although the half cycle shift was suggestive of a double mode, in our previous discussions of 4 CVn we preferred the hypothesis of a single mode which was re-excited at a new epoch. The reason was that numerical simulations of a double mode did not lead to an improvement of the fit to the photometric data relative to that of the single-frequency variable-amplitude hypothesis. At this stage, we cannot present evidence proving the correctness of either the single-mode or double-mode hypothesis, but note that in the light of the new BI CMi results, the double-mode hypothesis should also be kept in mind for 4 CVn. The reason for the less definite conclusions concerning the 7.37 cd-1 feature in 4 CVn is the relative lack of data available for the crucial 1974–1978 time period and the long time gap up to the 1983/4 data. We also note that 4 CVn contains other close frequency pairs such as the 6.404 and 6.440 cd-1 pair. However, if the double-mode hypothesis can be shown to apply to the 7.37 cd-1 feature, then the best two-frequency solution gives a separation of $`\mathrm{\Delta }`$f = 0.0015 cd-1. This would then be the smallest separation observed in $`\delta `$ Scuti stars so far.

## 5 Possible reasons for close frequency pairs

There are a number of possible explanations for the close frequency pairs:

(i) Mixed modes: Theoretical models predict a very large number of excited nonradial modes in evolved $`\delta `$ Scuti stars due to the dual nature of these modes (p modes in the envelope and g modes in the core). This dense spectrum of modes leads to small frequency spacing. The pulsation models of 4 CVn (see Breger & Pamyatnykh 2002) predict 554 unstable modes with radial orders $`\mathrm{}`$ = 0 to 2 over a 7 cd-1 range. Since BI CMi is in a similar stage of evolution, the number of modes can be regarded as an estimate for this star as well. The average spacing of $``$ 0.01 cd-1 is of the same order as the separation of the close frequencies found in this paper. If the theoretically predicted modes are indeed unstable, the close frequencies found may find a natural explanation.

We note here that the large number of closely spaced modes is predicted by more than one group using independent codes, e. g., the models for evolved $`\delta `$ Scuti stars by Guzik, Bradley & Templeton (2000) show a similar large number of predicted mixed modes.

The problem with this explanation lies in the fact that the large number of theoretically predicted modes of mixed character have not (yet) been discovered. The $`\delta `$ Scuti stars select only a small number of the possible modes to pulsate with detectable amplitudes and these sometimes show relatively large amplitudes. It is not obvious why the (unknown) mode selection mechanism favors so many close modes (but see below).

(ii) Trapped modes: Breger & Pamyatnykh (2002) have explored the possibility that the modes with low kinetic energy are preferentially selected since they are excited more easily. These are the modes trapped in the outer envelope and the analogs of pure acoustic modes. In this hypothesis, only the modes with the lowest kinetic energy would be selected out the many possible modes. The effect is very strong for $`\mathrm{}`$ = 1 and weaker for $`\mathrm{}`$ = 2 modes. This scenario has been successful in explaining the observed $`\mathrm{}`$ = 1 frequency spacing in the star 4 CVn and is promising for the $`\mathrm{}`$ = 2 modes. Again, the arguments could also apply to BI CMi because of its similar evolutionary status.

In spite of the severely reduced number of excited low-degree modes in the trapped mode hypothesis, two modes with close frequencies are still possible. Calculations of the kinetic energy associated with all the possible modes show that two almost identical modes with neighboring frequencies can have the same low kinetic energy. Both modes might, therefore, be excited.

The average spacing of the mixed-mode hypothesis applies here as well, so that the observed frequency separations can be matched. However, in contrast to the mixed-mode hypothesis, close frequencies would be predicted to occur only for trapped modes, i.e., at a relatively small number of predictable frequencies. We are, at present, engaged in testing this hypothesis with pulsation calculations.

(iii) Rotational splitting: Close frequency pairs can, in principle, be explained by rotational splitting with only two of the split modes visible. A simple calculation independent of the details of stellar models appears to rule out this hypothesis. The transformation from a coordinate system rotating with the star to the observer’s coordinate system leads to a frequency separation near $`m\mathrm{\Omega }`$, where $`\mathrm{\Omega }`$ is the frequency of rotation. Since almost all $`\delta `$ Scuti stars rotate faster than 10 km s-1, the predicted frequency splitting is too large. However, this calculation may be too simple because of the extreme splitting asymmetries which can be produced in rapidly rotating stars (e. g., Goupil et al. 2000). Such asymmetries produce rotational splitting which can be very much larger or smaller than that given by the $`m\mathrm{\Omega }`$ term. An observational test of this explanation and subsequent theoretical modelling requires reliable observational identification of the the $`m`$ values associated with each mode. These are not available.

(iv) The Small Spacing: Two modes with only a small frequency difference are expected if their $`\mathrm{}`$ values differ by 2, at least in the asymptotic case (e. g., Tassoul 1980). This is shown by

$$\nu _n\mathrm{}(n+\mathrm{}/2+1/4+\alpha )\mathrm{\Delta }\nu ,$$
(1)
where $`n`$ is the radial order and $`\mathrm{}`$ the degree of the pulsation mode, $`\mathrm{\Delta }\nu `$ the (large) spacing between successive orders, and $`\alpha `$ depends on the phase change at the stellar surface and is, in general, a slowly varying function of frequency. This is related to the so-called Small Spacing, which is expected to be an important asteroseismological tool for $`\delta `$ Scuti stars (see Christensen-Dalsgaard 2000). The fact that many $`\mathrm{}`$ = 0 and 2 as well as $`\mathrm{}`$ = 1 and 3 pairs are separated by only very small frequency differences could provide an explanation for the observed close frequencies. In $`\delta `$ Scuti stars the pulsations are of lower radial order, so that the asymptotic case does not really apply. Nevertheless, models show that if this explanation for the close modes applies, the observed frequencies would provide an important input into the pulsation models.

An important prediction of this explanation would be that the observed close frequency pairs consist of component modes with $`\mathrm{}`$ values differing by 2. Because of photometric cancellation effects across the surface, the only pairs that we could expect to observe photometrically would be (0, 2) as well as (1, 3). This can be checked if mode identifications are available. For the close frequency pair in BI CMi with the highest amplitudes, the mode identifications from phase differences (see Breger et al 2002) suggest $`\mathrm{}`$ = 2 for 10.429 cd-1 and $`\mathrm{}`$ = 0 (or maybe 1) for 10.437 cd-1. The available mode identifications for the other modes and the observed frequency spacings between the modes (i.e., the $`\mathrm{\Delta }\nu `$ value from equation 1) are at this stage not uniquely determined. Until more data become available, the $`\mathrm{\Delta }\mathrm{}`$ = 2 explanation cannot be excluded.

(v) Mode coupling: Frequency shifts can also be produced by nonlinear mode coupling by modes with the same or even different $`\mathrm{}`$ values (see Buchler, Goupil & Serre 1995, and Goupil 2000). These shifts could lead to close frequency pairs. This effect has had very little observational study and additional discussion is beyond the scope of this paper.

## 6 Conclusion

A literature search has shown that close frequency pairs with separations less than 0.06 cd-1 are common among well-studied $`\delta `$ Scuti stars. One of the best candidates for a more detailed analysis of this phenomenon is the star BI CMi. The variable-phase technique, which examines the phase jumps near the times of minimum amplitude of an assumed single frequency, was applied to test whether theses pairs are caused by separate excited stellar pulsation modes, single modes with variable amplitudes, or observational problems.

It was shown that at least three features are indeed pairs of separate pulsation modes beating with each other: at times of minimum amplitude the predicted phase jumps were observed. Both the observed amplitude and phase variations were predicted correctly by assuming two separate modes of nearly equal frequencies. Artifacts caused by observational error, insufficient frequency resolution or variable amplitudes could be ruled out. A fourth pair has a probable origin in two excited modes, while a 5th case is inconclusive due to long time scales of variability and small amplitudes.

The existence of close frequencies needs to be taken into account in planning the lengths of earth-based as well as space campaigns so that sufficient frequency resolution is obtained.

###### Acknowledgements.

It is a pleasure to thank Alosza A. Pamyatnykh for many interesting discussions. This investigation has been supported by the Austrian Fonds zur Förderung der wissenschaftlichen Forschung, project number P14546-PHY.


format: text

text:
# Complex shock structure in the western hot-spot of Pictor A

## 1 Introduction

The FR II radio galaxy Pictor A ($`z=0.0342`$) consists of two lobes dominated by hot-spots separated by $`4^{}`$ ($`210h^1\mathrm{kpc}`$) from the flat spectrum nucleus. Near the western hot-spot, there is a faint, $`24^{\prime \prime }`$ ($`10.8h^1\mathrm{kpc}`$) long radio and optical filament, “upstream of the hot-spot” first detected by Röser & Meisenheimer (1987). Wagner et al. (2000) found the filament to be highly polarised, confirming that it radiates via optical synchrotron emission.

The filament has been studied in detail at several radio frequencies by Perley et al. (1997). The radio emission at 8415MHz, figure 1, which bears a striking correspondence to the optical emission, also reveals a knot further upstream from the filament. This is of particular interest since at various timesteps throughout our simulation such a morphology is captured.

The discovery of optical synchrotron radiation, emitted by electrons with radiative life-times of about 100 yr (for magnetic fields of the order of equipartition), in a filament which extends several thousand light years perpendicular to the jet strongly indicates that the electrons are accelerated in situ. The crossing time from the hot-spot to the filament is much longer than the cooling time-scale, and efficient acceleration to optical emitting energies is required.

Our objective here is to ascertain whether purely hydrodynamical simulations can qualitatively reproduce the radio and optical morphology. Indeed on the basis of our simulations we suggest below that the filament is the result of a quasi-periodic pulsing behaviour of the terminal jet shock.

## 2 Summary of the Observations

Before discussing the simulations, we summarise relevant aspects of the observations and in the following sections we point out the relevance of the features in the simulations to the observational data.

A total intensity radio map (Perley et al. 1997), displayed in figure 1 shows the high-surface brightness western hot-spot; a bar-shaped region (the filament) east of the hot-spot, extends perpendicular to the inferred direction of the jet. The region between the hot-spot and the ridge line of the bar shows a “pedestal” of lower surface brightness.

High S/N observations were obtained with the VLT (UT1) during commissioning of FORS1. Figure 2 shows total intensity emission at 450 nm. Wagner et al. (2000) presented polarization maps derived from these data, showing up to $`60\%`$ polarization in the hot-spot, and up to $`40\%`$ polarization in the filament. Such a high degree of polarization cannot be the result of scattering and strongly suggests that the filament radiates optical synchrotron emission. Moreover, in the optical image, the bar-shaped filament is quite bright, and spatially coincides with a similar feature in the radio images. The eastern edge of the bar is very pronounced, but the surface brightness gradient in the western direction is more gradual. The filament is brightest at its ends and near the centre, where it overlaps the inferred jet direction. In the optical image there are faint extensions from the hot-spot extending in the direction of the ends of the filament.

## 3 Simulation Methodology

### 3.1 Hydrodynamics

The earliest simulations of astrophysical jets were published by Norman et al. (1982) and Smith et al. (1985). Their simulations revealed hot-spot morphologies that were different from the simple one-shock structure envisaged by Blandford & Rees (1974). These papers were seminal in leading to a good physical understanding of the properties of supersonic extragalactic jets.

In this work we have carried out similar simulations with a view to addressing in detail the time variability of hot-spot morphology. In particular we address the question of whether the Pictor A morphology could arise naturally in avlight supersonic jet. Our simulations were conducted using the VH-1 code (Blondin & Lukfin 1993) which is an implementation of the Piecewise Parabolic Method (PPM) (Colella & Woodward 1984). An advantage of PPM for this type of simulation is its excellent resolution of shocks.

The simulations were computed on an axisymmetric grid, with $`600\times 300`$ cells in the axial and radial directions. We performed calculations with nine different cases of the jet parameters, i.e. jet Mach number $`M=5,10,50`$ and the ratio of jet to background densities, $`\eta =10^2,10^3,10^4`$. In each case the jet density and pressure were normalized by the respective values, $`(\rho _0,p_0)`$, in the background medium. If we take $`T=10^7T_7\mathrm{K}`$ to be the temperature of the external medium, then the normalising value of the velocity is $`v_0=\sqrt{p_0/\rho _0}=\sqrt{kT/\mu m_p}3.65\times 10^7T_7^{1/2}\mathrm{cm}\mathrm{s}^1`$. The unit of length, $`x_0`$, is defined by the value of the $`\text{jet radius}=0.25`$, in numerical units. If we take the observed western hot-spot $`\text{FWHM}1.2^{\prime \prime }2.7\times 10^{21}\mathrm{cm}`$, for $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, as indicative of the jet diameter, then $`x_09.0\times 10^{21}\mathrm{cm}`$ and the unit of time, $`t_0=x_0v_0^19.4\times 10^6T_7^{1/2}\mathrm{yr}`$. As with simulations of this type, the exact physical parameters of velocity and density are not recovered. However, a high Mach number and low density ratio serve to reproduce, at least qualitatively, many of the features found in real jets.

Table 1 lists the jet parameters, time resolution and left boundary condition of our simulations. We define the grid to be $`10`$ units in the $`z`$ direction and $`5`$ units in the radial direction. Table 2 indicates the jet’s total mass flux, force and power for each choice of the basic parameters, in terms of fiducial values of the temperature, $`T_7`$, and the number density of the background medium, $`n_3=n/10^3\mathrm{cm}^3`$. For a jet with Mach number $`M`$, density contrast $`\eta `$, adiabatic index $`\gamma `$ (assumed equal to the ambient medium), and cross-sectional area $`A=\pi r_j^2`$, the initial flow velocity, mass flux, force and kinetic power are

$$v_\mathrm{j}=v_0\sqrt{\frac{\gamma }{\eta }}M,$$
(1)
$$\dot{M}_\mathrm{j}=A\rho _0v_0\sqrt{\gamma \eta }M$$
(2)
$$F_\mathrm{j}=A(\rho v^2+p)=A\rho _0v_0^2(\gamma M^2+1),\text{ and}$$
(3)
$$P_\mathrm{j}=A\left(\frac{1}{2}\rho v^3+\frac{\gamma }{\gamma 1}pv\right)=A\rho _0v_0^3\left(\frac{1}{2}M^2+\frac{1}{\gamma 1}\right)M\sqrt{\frac{\gamma ^3}{\eta }}.$$
(4)
These are subrelativistic expressions, although the jets of Pictor A may involve relativistic flows of plasma. One may relate the $`(\eta ,M)`$ parameters of a classical jet to an equivalent pair of parameters for a relativistic jet with equivalent thrust or power, e.g. Rosen et al. (1999), Carvalho & O’Dea (2001). However we are not modelling the system to determine the Lorentz factor of the Pictor A jet to great precision.

Each simulation is initialised with the high-velocity, low-density jet material extruding a small distance into the ambient medium. At every subsequent time step the conditions of the initial jet are recreated in the cells at the base of the jet on the left (inner $`z`$) boundary of the grid, and thus the jet’s injected fluxes of mass, momentum and energy remain constant in time. We performed simulations with two different choices for the condition applied at the left boundary outside the radius of the jet. (1) A reflecting (“closed”) boundary condition represents the effects of the symmetry plane through the nucleus, felt by a jet that is close to its origin and counterjet. (2) An “open” boundary allows free outflow of mass, momentum and energy, thereby representing systems where the hot-spot and bow shock are distant from the nucleus and not directly influenced by the opposite lobe. The outer $`z`$ and $`r`$ boundaries are open to outflow.

### 3.2 Ray-tracing

In our simulations, a scalar variable $`\phi `$, as introduced in Saxton et al. (2001), is passively advected with the gas. Cells occupied purely by jet material are initially assigned the value $`\phi =1`$; and the ambient medium initially has $`\phi =0`$. Thus $`\phi `$ traces the local fraction of matter originating in the jet.

We adopt the approximation of using the pressure field to trace the emissivity of the relativistic particles within jet plasma (e.g. Saxton et al. 2001). This approximation is based on the fact that the synchrotron emissivity, $`j_\nu `$ is proportional to $`pB^{1+\alpha }`$ where $`p`$ is the pressure, $`B`$ is the magnetic field and $`\alpha `$ is the spectral index. If the magnetic pressure, $`B^2/8\pi `$ tracks the particle pressure, then $`j_\nu p^{(3+\alpha )/2}`$, and $`\alpha 0.6`$ typically. Even if the magnetic pressure does not track the pressure faithfully, it is likely that regions of high pressure still correspond to regions of high emissivity.

We revolve the hydrodynamic data frames of $`\phi `$ and $`p`$ to form 3D cylindrical structures. Using $`\phi `$ to distinguish synchrotron-emitting jet material from gas originating in the background medium, the weighted emissivity is then $`j_\nu \phi p^{1.8}`$. By integrating $`j_\nu `$ along rays through the 3D volume we generate surface brightness images of the jet and cocoon at a selected orientation. The resulting images are compared with observations in §5.

## 4 Simulation Results

### 4.1 General features

The density distribution and velocity fields of representative frames from six of the simulations are shown in Figure 3. At the beginning of each simulation, a bow shock is quickly established around the jet. The regions just downstream and inside the bow shock are dominated by the shocked dense interstellar medium (ISM). The innermost regions immediately surrounding the jet are occupied by a cocoon of low-density gas moving in a turbulent backflow opposite the jet direction. Flow velocities in some regions are comparable to the initial jet velocity $`v_\mathrm{j}`$, and are much larger than the bow shock velocity (see Figure 4). Turbulent velocities are typically lower in the outer parts of the cocoon, further from the jet. The scalar tracer $`\phi `$ (see Figure 5) reveals that jet-derived matter dominates the inner regions of the backflow, whereas the outer regions just downstream of the bow shock are entirely composed of background gas.

The evolution of the backflow depends upon the left boundary condition. When the left boundary of the grid is open to outflow, the backflow occupies a cylindrical “sheath” surrounding and flowing opposite the jet, ultimately flowing off the grid. This “sheath” has a radius a few times that of the jet. It increases gradually as the jet progresses across the grid, but is insensitive to $`\eta `$. Very little of the surrounding gas mixes with backflowing jet plasma: most of the fast-moving backflow still exhibits high concentrations of jet plasma ($`\phi 0.5`$) when it exits the left boundary (see right column of Figure 5).

In simulations with a closed, reflecting boundary, cocoon eddies accumulate inside a roughly conical volume, which is widest near the left boundary, and tapers towards the head of the jet. This “wake” cocoon expands laterally as the jet progresses forwards, and jet material accumulates and mixes with the shocked ISM. The eddies are typically smaller and have lesser concentrations of jet matter (lower $`\phi `$) with increasing distance from the head of the jet (see upper panels of Figure 5).

The density contrast of the jet, $`\eta `$, affects the gross structure of the surrounding radio emitting regions. If the left boundary is closed, jets with lower $`\eta `$ produce wider bow shocks and wider “wake” cocoons (see left panels of Figure 3). However in cases with an open boundary the high-$`\phi `$, fast backflowing “sheath” cocoon has approximately the same width, $`4r_\mathrm{j}5r_\mathrm{j}`$, for all the values of $`\eta `$ we have studied: $`\eta =10^2,10^3`$ and $`10^4`$ (see right panels of Figure 5). In such systems the bow shock widens marginally with decreasing $`\eta `$.

The amount of mixing between jet and ISM also shows a dependence upon $`\eta `$. For the open boundary mixing per unit length is more pronounced in the backflow for the larger values of $`\eta `$ (see the right panels of Figure 5).

There is a similar dependence for the case of the closed boundary with some additional features. For $`\eta =10^2`$, the cocoon widens almost linearly from the head of the jet (radius $`3r_\mathrm{j}`$) and the degree of mixing in the backflow increases more or less gradually with distance from this point. However in cases with $`\eta =10^4`$, the head of the jet is persistently surrounded by a broad (radius $`6r_\mathrm{j}`$), low-density, high-$`\phi `$ ellipsoidal region that is essentially the first turbulent cell in the backflow. (see upper-left panels of Figure 35). There is a lot of mixing at the trailing edge of this cavity and average $`\phi `$ values decline rapidly to the left of that region. Cases with $`\eta =10^3`$ show an intermediate behaviour: the cavity is more extended in the $`z`$ direction and includes more fingers of entrained dense gas.

In the ultra-light ($`\eta =10^4`$) jets with an enduring cavity, the head of the jet surges back and forth along the axis between the front and back ends of the cavity. Similar surging behaviour occurs in most of the simulations, to different extents and at different frequencies. This cyclic process begins when turbulence in the backflow briefly compresses and constricts the jet at one point (see Figure 7, panels 2 & 3). A portion of the front of jet is cut off, and its remnant mixes into the front of the cocoon. Jet plasma accumulates upstream of the obstruction (Figure 7, panel 4), until it is able to surge through to a point much closer to the bow shock (Figure 7, panels 5 & 6). Thus the terminal shock of the jet is not a constant feature: it may occur close to the bow shock (when the jet has just undergone a surge), or it may occur many jet diameters further upstream (when the jet is being disturbed by turbulence). The extent and frequency of surging depends on the jet parameters and the condition of the left boundary. In the closed-boundary, $`\eta =10^4`$ cases, the strong, persistent, frontal eddy frequently pinches off the jet at the back end of the cavity, where the backflow converges back towards the surface of the jet. Cases with $`(\eta ,M)=(10^4,50)`$ exhibit unusually extreme surging behaviour, because the relatively violent turbulence throughout the entire cocoon can pinch this ultra-light jet at almost any point along its length.

### 4.2 Shocks & emission regions

As discussed above (§3.2), the regions of strongest emission are high pressure shock features with a significant concentration of jet plasma (i.e. large $`\phi `$). In agreement with past studies (Norman et al. 1982), we find several types of shocks characterising particular regions of the jet and its surroundings. The plots of $`|p|`$ in Figure 16 show the shocks appearing as intense lines, in several regions of the flow.

We classify the different shocked regions as follows:

1. Firstly there is the bow shock driven by and encompassing the jet and its cocoon. The bow shock is persistent but is clearly not involved in the radio emission.
2. The strongest but most ephemeral shocks occur near the head of the jet where the jet deposits a large part of its energy and momentum. This feature pulsates back and forth along the axis over distances several times the jet diameter, sometimes disappearing and reforming.
3. Shocks also propagate in the cocoon and either originate from near the head of the jet or are the result of local interactions between the jet, regions of the turbulent backflow and the surrounding layer of dense gas. Annular shocks are seen propagating both forwards and backwards through the backflow.
4. The jet itself is criss-crossed by diamond shocks caused by momentary and local compression by the turbulent backflow. As implied by the classical perturbation analysis of Birkhoff & Zarantonello (1957), the diamond shocks are more widely separated in jets with greater $`M`$.

The distribution of $`\phi `$ informs us which of the shocks produce significant synchrotron emission. Diamond shocks are always bright because they occur within pure jet plasma. For the same reason, the jet’s terminal shock yields a prominent hot-spot in most frames of the simulations. However wider luminous structures are restricted to the extent of the high-$`\phi `$ “sheath” cocoon in the systems with an open boundary, or the relatively unmixed frontal regions of the “wake” cocoons arising in the systems with a closed boundary. Thus, depending on the jet parameters and the boundary condition, there are no bright features exceeding $`2r_\mathrm{j}7r_\mathrm{j}`$ in radius (see Figure 5).

In summary, allowing for the axial symmetry of our simulations, the kinds of bright three-dimensional structures that are possible include: bright point-like, disk-like and cone-like shocks within the jet; in the backflow there may occur a diversity of rings and circular ribbons, with potentially flat, curved or solid cross-sections. In the following section, we consider how these structures may appear in two-dimensional projections onto the sky, and how these images may compare with observations.

## 5 Relationship to Pictor A

In this section we compare rendered images with the morphology of bright structures evident in the observations of Pictor A. The optical observations (Figure 2) are a guide to the current positions of the shocks. Because of the longer lifetime of radio emitting electrons the radio image (Figure 1) is indicative of the overall distribution of jet plasma within the cocoon.

In the rendering of our simulations we seek physical configurations that show a bright hot-spot at the front of the jet and a transverse bar upstream of the hot-spot. As refinements, we prefer images in which the bar appears brightened at its edges. The optical image (Figure 2) also suggests faint arcs connecting the ends of the bar to the hot-spot. We also have in mind (in the optical data) the bright knot of emission at the intersection of the jet axis and the bar, and the faint tail of emission leading from that knot to the east.

We generated ray-traced images of each frame of each simulation, with the jet oriented at a variety of inclinations, $`\theta =45^{},60^{},70^{},80^{},90^{}`$ ($`\theta `$ is the angle between the jet and the line of sight). We inspected and compared all the resulting images and sorted them somewhat subjectively according to their consistency with the topology, proportions and brightness ratios seen in the Pictor A observations. In the following we consider the effects of orientation upon the quality of the matching rendered images, and the effects of the jet parameters. We also inspect the underlying physical nature of the bright features and their temporal evolution.

### 5.1 Orientation

The orientation of the jet has a substantial effect on the projected appearance of the axially symmetric luminous features surrounding the jet.

The hot-spot and knot are not decisive characteristics, since diamond and terminal shocks appear at numerous instances and their appearance is virtually independent of orientation. However the transverse bar in the Pictor A observations provides a tight constraint on the selection of appropriate rendered images. This elongated feature must be the projection of at least one luminous structure that is three-dimensionally extended. Therefore it is necessary but not sufficient to select snapshots that contain extended high-pressure features. When projected onto the sky at some orientation, only a subset of these configurations produce images similar to the observations.

When the jet is at a large angle to the plane of the sky, e.g. $`\theta =45^{}`$, very few frames show the desired topology and brightness: typically a few times $`0.1\%`$ of frames have a subjective resemblance to the images. The best images of jets at $`45^{}`$ inclinations are shown in Figure 8. In these cases, the appearance of a straight “bar,” if any, is a coincidence due to the visual superposition of the near and far edges of unconnected rings in the backflow; it is not a single physical structure. One of the positive features of the better renderings at $`\theta =45^{}`$ is the abundance of snapshots in which the hot-spot is projected in front of an edge of a bright ring in the backflow, thereby resembling the faint arcs connecting the hot-spot to the bar in the optical observations. However the occurence of conspicuously elliptical features throughout the image tends to ruin the overall resemblence, and the accidental bars appearing by superposition are never much wider than the hot-spot. Indeed the least-bad matches have small features compared to the initial jet, $`r2r_\mathrm{j}`$, and we speculate that the resemblences might be even less convincing if these structures could have been resolved more finely by the simulations.

When the system is viewed at inclinations $`70^{}\theta 90^{}`$, the matches are more abundant and better in quality. Approximately $`4\%`$ of all frames rendered at $`\theta =90^{}`$ show features qualitatively similar to the Pictor A observations. Approximately $`3\%`$ of the frames at $`\theta =70^{}`$ show the same quality of resemblance. Inclinations of $`\theta =90^{}`$ produce images with a straight, bright transverse bar behind the hot-spot. Lower inclinations, progressively down to $`\theta 70^{}`$, yield more images in which arcs appear to connect the ends of the bar with the hot-spot. However the trade-off is that the apparent bar may be less straight and distinct.

At inclinations as low as $`\theta =60^{}`$, frames that discernibly resemble observation are as rare as $`0.6\%`$. Some of the best examples are presented in Figure 9, but they are not compelling matches. Like the best morphological matches at $`\theta =45^{}`$, their bar radii are $`2r_\mathrm{j}`$, and they tend to involve coincidental superposition of the far side of an annulus behind the hot-spot. There are instances showing a favourable “pedestal” or “arcs” near the hot-spot, but the best simulated “bar” features are plainly elliptical rather than straight.

At a given inclination, good morphologies appear with roughly equal abundance for all choices of the jet parameters, with the exception of the extreme case $`(\eta ,M)=(10^4,50)`$ where the turbulent backflow and the jet behave abnormally (see §5.2). For example renderings, see Figures 10-15, and a selection of corresponding plots of $`|p|`$ in Figure 16, which indentify the principal intense linear features with shocks. We display more examples for cases where $`\eta =10^4`$ but this is only because these simulations produced more frames of data, and the emission features evolve more rapidly compared to the jet’s advance across the grid (see §5.3 below).

As implied by the distribution of $`\phi `$ and the commonplace occurrence of annular shocks throughout the backflow, the rendered snapshots include bright bars of diverse widths, (up to approximately $`6r_\mathrm{j}`$ in cases with $`\eta =10^4`$). The observed filament (Figures 12) has a radius four times that of the hot-spot. If the observed hot-spot corresponds to the width of the jet, then all of our simulations are in principle able to produce bright shocks in the backflow that are consistent in size with the observed filament.

### 5.2 Brightness distribution

Let us now consider the overall brightness distribution of various features in the image. This includes the brightness distribution of the filament, its brightness relative to the hot-spot and the overall brightness of the hot-spot and bar relative to the lobe. Note that in the observations, the brightest part of the filament is approximately $`7\%`$ of the peak brightness in the hot-spot, in the radio and about $`2\%`$ in the optical.

Some frames rendered from the simulations have shapes that are consistent with the observations of Pictor A, in terms of qualitative topology and approximate proportions, but may differ from the observations in terms of brightness. For instance, in some imperfect examples, the bar is dim on the edges rather than being edge-bright, e.g. $`(\eta ,M,t,\theta )=(10^4,10,0.6216t_0,90^{})`$ (in centre panel, upper row of Figure 11, and right column, Figure 23) and $`(\eta ,M,t,\theta )=(10^4,5,0.6221t_0,90^{})`$ (centre panel, lower row of Figure 11). In these cases, the bar is the image of a disk-like (rather than annular) shock that cuts the jet and a high-$`\phi `$ transverse stream diverging from the jet.

In some cases in which the morphology is approximately correct, the ratio of filament to hot-spot brightness is larger than observed: e.g. $`(\eta ,M,t,\theta )=(10^2,10,0.296t_0,90^{})`$ shown in the central panel of Figure 14; and $`(\eta ,M,t,\theta )=(10^2,10,0.8212t_0,90^{})`$ (right column, Figure 22), in which the wide bar is a shock cutting through both the jet and cocoon. However, we cannot rule out cases with overly bright bars since relativistic beaming may be important in the hot-spot. There are also examples in which the ratio of filament to hot-spot brightness is approximately correct, e.g. the three examples in the uper row of Figure 15. Here the brightness ratios are approximately $`12\%`$, $`8\%`$ and $`12\%`$ (ordered from left to right). The appropriate ratios are preferentially found in the highest Mach number cases. A filament at the $`7\%`$ level is generally more difficult to find in a low Mach number lobe because of the lack of contrast with the background cocoon. The physics of the contrast between hot-spot and cocoon is discussed further below.

The observed filament is extended more to the south of the inferred jet than the north. Because of the axisymmetry of the simulations, this feature is not replicated. However, an axisymmetric approximation is justified in view of the very straight jet evident in the X-ray image (Wilson et al. 2001). Nevertheless, in most of the matches that we have presented in Figures 10-15, the simulated filament is seen to be brighter towards the ends. This is the result of the annular shocks arising in the region of the hot-spot.

A feature of the observations that is replicated in approximately $`40\%`$ of the matches in Figures 10-15 is a small “knot” of emission at the centre of the bar. Such a feature is often the result of the superposition the bright annulus with a section of the jet that is also bright, such as a shock upstream of and related to the terminal shock of the hot-spot.

Yet another observational feature that appears in almost every panel in Figures 10-15 is the radio knot south-east of the bar along the inferred direction of the jet. These knots correspond to diamond shocks in the jet; their exact location depends stochastically upon the details of the turbulence in the cocoon.

In the radio observations (Figure 1), diffuse emission from the radio lobe in regions near the bar and hot-spot is $`0.5\%`$ of the maximum surface brightness of the hot-spot. As Perley et al. noted, such a large contrast is highly unusual for a radio galaxy. Ideally, our ray-traced images should produce a similar contrast between the brightness of bar and hot-spot and the wispy rings and dim emission pervading the rest of the cocoon.

Ideally, higher Mach number jets produce hot-spots with a higher contrast to the lobe. In simple terms, if the jet terminates in a single steady shock, the pressure contrast between the hot-spot and the cocoon is given by $`p_{\mathrm{hs}}/p_{\mathrm{cocoon}}M^2`$. This aspect of jet physics lead Perley et al. to suggest that the Mach number may be as high as $`40`$. In terms of a realtivistic jet, this is equivalent to a Lorentz factor of 24. However, the inferred parameters of the hot-spot (especially the pressure) could be significantly affected by beaming. The beaming pattern may be quite complex in view of the velocity field produced by the terminal oblique shock (see §5.4).

Ray-traced images from the simulations with parameters $`(\eta ,M)=(10^4,5)`$ (see Figure 10) give bar and hot-spot features that are typically only $`10`$ times brighter than the background cocoon. The brightness contrast is greater in simulations with $`(\eta ,M)=(10^4,10)`$ (see Figure 11). In this case the contrast varies bewteen $`10^{1.3}`$ and $`10^{1.7}`$. Simulations of very fast jets, (e.g. case C2L, with parameters $`(\eta ,M)=(10^2,50)`$) can yield brightness contrasts in the range, i.e $`10^{22.5}`$ – consistent with the radio images of Perley et al. The tendency towards increasing contrast in the hot-spot and bar can be seen in the example contour maps of surface brightness for $`\eta =10^2`$ and for the lower-$`M`$ cases of $`\eta =10^4`$ in Figure 17.

However this trend breaks in the case C4L where $`(\eta ,M)=(10^4,50)`$; the brightness contrast is exceptionally low – typically $`30`$ – compared to the more sharply defined features in the $`(\eta ,M)=(10^2,50)`$ and $`(\eta ,M)=(10^3,50)`$ simulations. This is an extreme manifestation of the surging behaviour that we have already noted. In this case the backflow turbulence is particularly vigorous and remains so at locations well behind the front of the bow shock. Strong eddies pinch off the extremely light jet at numerous locations along its length, and the jet’s penetration into the cocoon varies dramatically within a dynamical timescale. Consequently, jet energy is deposited widely throughout the cocoon, rather than being concentrated at a quasi-stable head; moderately bright but transient rings and hot-spots occur throughout a large volume. This aspect may well be a feature of the restriction to axisymmetry causing the two-dimensional turbulence within the cocoon to be focussed on the jet.

### 5.3 Temporal variation

Physical configurations producing the desired morphology of surface brightness are ephemeral products of the jet’s pulsating head and turbulent motions within the cocoon. It is therefore interesting to consider how these transient features vary in time. These considerations are also motivated by the potential of using hot-spot brightness ratios as a measure of relativistic beaming. Although these calculations are not relativistic we can examine the intrinsic distribution of hot-spot intensities. If it were to be the case, for example, that order of magnitude fluctuations in the ratio of hot-spot luminosities were to be likely, then the use of luminosity ratios as a measure of relativistic speeds within hot-spots would be limited. We address this point below following consideration of the structure function of individual hot-spots.

The dynamical timescale of structures in the jet material, such as the hot-spot and shocks in the jet and backflow, characteristically depends upon the jet’s diameter and internal sound speed,

$$t_{\mathrm{dyn}}\frac{2r_\mathrm{j}}{c_{\mathrm{s},\mathrm{j}}}0.39\sqrt{\eta }t_0.$$
(5)
Thus $`t_{\mathrm{dyn}}0.04t_0`$ for jets with $`\eta =10^2`$, and $`t_{\mathrm{dyn}}0.004t_0`$ for jets with $`\eta =10^4`$. Since the brightest feature is usually a terminal hot-spot or some other similarly compact structure, we expect the instantaneous peak surface brightness to vary in time on a time scale on the order of magnitude of $`t_{\mathrm{dyn}}`$. Indeed this is the case, as shown by the structure functions (Simonetti et al. 1985) of the time-variation of the peak brightness of images rendered at $`\theta =90^{}`$ (see Figure 18. The structure functions were calculated from a time series based upon the peak intensity in each frame divided by the long term mean. (Frames containing initial transients were omitted.) Hence the amplitude of the structure function gives one an estimate of the square of the relative variabiliity on the relevant timescale.

On time scales longer than $`t_{\mathrm{dyn}}`$ the behaviour of most structure functions is effectively white noise (power index $`0`$); whereas the behaviour is like “red noise” on shorter time-scales. The transition timescale between these behaviours is of order $`0.10.5t_{\mathrm{dyn}}`$ and depends little on $`M`$. However jets with greater $`M`$ generally show greater relative variability.

The one exceptional case is the $`(\eta ,M)=(10^4,50)`$ simulation with the closed boundary. In this simulation there is long term evolution (evident in the upturned curve at long time-scales in the upper left panel of Figure 18) related to the highly turbulent cocoon upon which we have already commented in §5.2.

In all cases bar the exceptional case, the left boundary condition has negligible effect on the variation of peak surface brightness (compare left and right columns of Figure 18). This is the result of the proximity of the brightest structures to the terminal jet shock, rather than parts of the flow directly affected by the left boundary.

In physical units, with $`t_010^7T_7^{1/2}\mathrm{yr}`$, the dynamical timescale is $`t_{\mathrm{dyn}}4\times 10^{5(4)}\mathrm{yr}`$ for a jet with $`\eta =10^{2(4)}`$. These time-scales are much shorter than the time taken for light or jet plasma to traverse the $`2.1\times 10^2/h\mathrm{kpc}=6.4\times 10^5/h\mathrm{ly}`$ projected distance between the nucleus and either hot-spot. Therefore the intrinsic flickering of two opposite hot-spots are causally independent. Figure 19 shows the cumulative distribution function of the ratio of the peak brightness of one spot with respect to the other, $`R=I_1/I_2`$, calculated from simulation output according to the derivation in Appendix A. The ratio has a more extended distribution in cases with larger Mach number, consistent with the larger amplitude of the structure function. For a given choice of jet parameters $`(\eta ,M)`$, an open boundary condition gives rise to substantially greater variability. For instance, the probability that $`R>4`$ or $`R<\frac{1}{4}`$ (see Table 3) is approximately $`20\%`$ for cases with $`M=50`$ and an open boundary but only ranges from $`413\%`$ in the cases we studied with a closed boundary. Likewise for small Mach number cases, for example with $`(\eta ,M)=(10^3,5)`$ the probability that $`R>4`$ or $`R<\frac{1}{4}`$ ranges from $`37\%`$ with an open boundary but only $`0.41.5\%`$ with a closed boundary. Thus in cases where the jet has a low Mach number, the relative brightnesses of hot-spots on different sides of a source is a reasonable measure of non-intrinsic effects, including relativistic beaming, however intrinsic variability may be significant if the Mach number is large ($`M10`$). In gemeral terms these distributions provide a means of quantifying the uncertainty in using hot-spot ratios in various applications.

### 5.4 Underlying structure and activity

The reproduction of images resembling observations of the western hot-spot of Pictor A requires a strong shock on the jet axis, and some other, wider shock structure in high-$`\phi `$ regions of the backflow. However more than one kind of physical situation can produce an appropriate image. Here we consider the particular instantaneous arrangements of shocks and jet plasma that yield some of the example images selected and displayed in Figures 10-15. Examples of the related dynamics are given in Figures 20-24, in which the panels show simulated intensity maps along with corresponding snapshots of pressure, $`\phi `$ and density with overlaid velocity vectors. (In referring to these figures we use L for left and R for right panels respectively.)

In practice, a bar and hot-spot may arise from different underlying structures at any stage of the jet’s surging behaviour. A hot-spot is commonly the termination shock at the head of the jet. However amongst our $`\eta =10^4`$ selections, there are some cases such as Figures 20L, 21L, 21R, and the detailed case Figures 25-27 where the head of the jet is at the back of the leading cavity, and not directly connected to the hot-spot. In such cases the forward hot-spot consists of a high pressure region of moderately mixed jet and thermal gas, or is an isolated blob of high-$`\phi `$ matter previously separated from the jet.

A bar may result from a variety of different kinds of annular or disk-like structures. In some cases (Figures 20L, 20R, 24L) the bar is a broad shock at the head of the jet, or a body of jet plasma near the head. In some cases, it is the shocked rim of a funnel-shaped or disk-shaped volume of plasma that has split off the head of the jet but yet remains largely unmixed with the thermal gas in the cocoon (see Figures 21R, 22R, 23R, 25-27). In other instances, the bar is a wide annular region in a part of the backflow well separated from the head of the jet, when such a region either preserves a moderate or high concentration of jet plasma, or experiences a strong shock locally (see Figures 22R, 21L ).

There are cases in which an ideal bar and hot-spot are formed, but additional features spoil the resemblance to Pictor A. For example, in one selection at $`t=0.4085t_0`$ (Figure 23L; or upper right panel of Figure 11) the bar and hot-spot are high-$`\phi `$ shocks in the vicinity of the jet’s head, but an additional strong reverse shock in the jet appears as another bright feature to the left of the bar. Such reverse shocks are an occasional by-product of the activity that produces a hot-spot as a forward shock and annular shocks propagating into the backflow. Milder versions of the reverse shock may give the appearance of the “knot” or “pedestal” observed on the jet axis.

It is not possible within the space of one paper to demonstrate the physical details underlying the evolution of all the physical configurations that give the appearance of a bar and hot-spot. In order to describe one particular process for the formation of a bright hot-spot and bar, we show in Figure 25, a time sequence of linearly scaled pressure images from the closed-boundary, $`(\eta ,M)=(10^4,5)`$ simulation with fine time resolution (simulation c4V). The first frame is at $`t=0.9952t_0`$ and subsequent frames are at intervals of $`0.0006t_0`$ numerical time units. Figure 26 shows corresponding frames of the distribution of $`\phi `$. Figure 27 shows rendered images of the surface brightness at $`\theta =90^{}`$. Figure 28 shows corresponding flow velocity fields superimposed upon $`\mathrm{log}\rho `$ greyscale images.

At the beginning of the sequence (frame 1), we see the usual haze of emission from weak shocks in the $`\phi 0.5`$ cavity at the head of the cocoon. A bright hot-spot appears at the front of the cavity and this is a high-pressure, shocked region with $`\phi 1`$, left behind by the head of the jet when it was previously at the foremost extremity of its surging motion. In this frame, the head of the jet is towards the rear of the cavity. The jet splits at a terminal shock, dividing into a funnel-shaped structure consisting of unmixed jet plasma moving in a radially outwards, positive-$`z`$ direction, at speeds comparable to $`v_\mathrm{j}`$. The outermost edge of this structure is curled back upon itself, and trails off into the backflow closest to the surface of the jet. The region between this backflow and the inner surface of the funnel is underpressured compared to average conditions in the cocoon, however the outer lip of the funnel is a high pressure surface.

In the second and third frames, the edge of the funnel has advanced outwards. In the third frame it collides with a flow of $`\phi 0.5`$ matter moving in the negative-$`z`$ direction near the outer surface of the cocoon. The flow of gas around the funnel causes a strong annular shock which also results in low presures at smaller $`r`$ and $`z`$. The shock produces significant emission, and the result is a bright bar visible in images rendered at $`\theta 90^{}`$.

During the activity that creates the bright ring, the diamond and transverse shocks moving in both directions within the jet produce several brighter knots of emission. The knots move and fluctuate in brightness, but they maintain an approximately constant separation. Frames 3 to 5 also show a hazy skirt of emission surrounding the jet to the left of the brightest knot; it is a backflowing, high-$`\phi `$ body of fluid that previously came off the edge of the funnel.

In frames 5 and 6 the front hot-spot fades, as it spreads laterally, loses pressure and undergoes a reduction in $`\phi `$ as a result of mixing. Meanwhile, the bright bar that persisted from frames 3 to 4 has also faded. The interaction between the funnel’s edge and the well-mixed, outer-front backflow has progressed. The shocks of the initial interaction have separated and spread into a wider and complicated shock structure. The outer edge of the high-$`\phi `$ funnel is swept in the negative-$`z`$, radially outwards direction by the well-mixed backflow from the head of the cocoon. In later frames (not shown), the funnel structure is greatly disrupted.

## 6 Discussion

### 6.1 Orientation

In this paper, we have presented a detailed analysis of simulations of light jets with a wide choice of parameters with the main aim of determining the cause of the filament upstream of the western hot-spot in Pictor A. The dominant physical principle underlying this analysis is that complex terminal shock structure is manifest in simulations of light supersonic jets and that, in projection, such a shock complex may resemble the Pictor A radio and optical observations of the region near the western hot-spot. Similar complexity has been realised in three dimensional simulations of light supersonic jets by Tregillis et al. (2001). One of the advantges of restrictio to a two dimensional simulation, at this stage, is the ability to map out a larger region of parameter space and to utilise higher spatial resolution. In particular, we have carried out simulations with a density ratio of $`10^4`$; these would be quite expensive in three dimensions.

One of the most important constraints on the simulations is that the bar is a thin structure and if it is three dimensional, it would appear to be almost edge-on. Our rendering of synthetic radio images at different orientations bears this out. Adequate morphologies are not reproduced at low inclinations. Another decisive requirement that implies an almost edge-on structure, is reproduction of a filament with low curvature.

At inclinations $`\theta <70^{}`$, geometric effects cannot convincingly produce a long, straight bar frequently enough . The projection of an annular structure with a ribbon-like cross-section makes one side brighter than the other, but the image is still significantly curved. The accidental coincidence of the near and far sides of separate rings in projection may give the appearance of a straight feature, but these tend to be narrower than the observed filament. A structure with a very large radius of curvature might appear straight, but luminous structures of this type occur in our simulations in regions where the coocon is no more than $`2r_\mathrm{j}`$ in radius. Three dimensional hydrodynamical effects may admit more complicated explanations for bar-like morphology. However extreme departures from axial symmetry in the backflow seem unlikely, given the strightness of the jet as indicated by the X-ray observations. On the basis of non-relativistic hydrodynamics, we therefore think it likely that the inclination angle, $`\theta >60^{}`$.

If there is any easy way to produce the appearance of a bar in a basically axisymmetric backflow viewed at $`\theta <70^{}`$, it probably requires relativistic effects. A ring or disk moving relativistically appears more edge-on than if it were at rest (see Bicknell & Begelman 1996). This would require the filament to be transient structure moving with a pattern speed larger than the expected hot-spot speed $`0.10.2c`$. Relativistic beaming could also affect the brightness distribution if the matter passing through the filament and hot-spot shocks moves relativistically, brightening an annular structure on whichever side has a post-shock flow directed towards our line of sight. However the significance of these relativistic effects, for hot-spots in general and for Pictor A in particular, remains to be demonstrated, and is beyond the scope of our present study.

### 6.2 Implications for physical parameters

Although all choices of the jet parameters, jet density ratio, $`\eta `$, and Mach number, $`M`$, yield instants with appropriate morphology, the brightness characteristics of the resulting images may distinguish the most plausible values of the Mach number. Jets with larger Mach number typically produce terminal shock features with higher peak intensity. Perley et al. (1997) had deduced a high Mach number ($`M40`$) for the jet in Pictor A based on the large ratio of the minimum pressures in the hot-spot and lobes. We have confirmed this in our simulations that show, for example, a brightness contrast $`100`$ for $`(\eta ,M)=(10^2,50)`$ (see §5.2). Nevertheless, the unusual $`(\eta ,M)=(10^4,50)`$ simulation produces a lower contrast and widely distributed rings in the backflow, and these features are clearly ruled out by observation. This may mean that the real jet parameters are less extreme, or else the powerfully turbulent backflow that occurs in this case may be less disruptive in three dimensions.

We do not find such a clear diagnostic of the density contrast of the jet. Cases with lower $`\eta `$ tend to produce broader cocoons, especially the rounded cavity in $`\eta =10^4`$ simulations. For jets with larger $`\eta `$, the jet plasma and shocked thermal gas mix more rapidly in the cocoon. However all of the parameter sets we investigated permit the appearance of ephemeral bar-like sructures with a range of widths, including some with $`4r_\mathrm{j}`$ radius as required by the observations.

Many of our simulations imply relativistic bulk flow, $`\beta 1`$ (see Table 2). One aspect of the observations where relativistic effects may be relevant concerns the large contrast between the hot-spot and radio lobe empshasised by Perley et al. (1997). On the basis of hydrodyamics alone, this would indicate a Mach number of approximately 40-50, or equivalently a Lorentz factor of approximately 25. This has motivated the $`M=50`$ simulations described above. Nevertheless, relativistic beaming may substantially contribute to the observed hot-spot brightness. This aspect may need to be taken into account with future simulations with a fully relativistic code.

### 6.3 The orientation of the jet and filament in relation to radio and X-ray observations

In their analysis of the Chandra X-ray Observatory data on Pictor A, Wilson et al. (2001) have suggested that an inclination $`\theta 23^{}`$ is viable. Larger angles are possible if jet magentic fields considerably less than equipartition are allowed. An inclination angle as small as $`23^{}`$ would be in direct contradiction to our conclusion that the Pictor A radio structure is viewed at an angle of $`\theta 70^{}90^{}`$. It is therefore worthwhile to examine the question of orientation in some detail.

The first part of the Wilson et al. (2001) argument relates to the arm-length ratio of the radio source. Assuming a velocity of advance for the hot-spot, $`\beta _{\mathrm{HS}}0.11`$ (Arshakian & Longair 2000) and the measured arm-length ratio of the source, Wilson et al. (2001) derive $`\theta 23^{}`$ but caution that this estimate by itself should not be taken too seriously for an individual source. The main limitation is probably the following: The statistical result for hot-spot advance includes sources in which the momentum of the jet is spread over a wide area as a result of the “Dentist Drill” effect (Scheuer 1982), thereby slowing its average advance. However, as we have remarked in justifying an axisymmetric simulation, the western jet in Pictor A appears to be very straight. That is, it does not appear to be substantially deflected and therefore the advance of the hot-spot may be faster than the average Arshakian & Longair (2000) value. Moreover, in contrast to the western lobe, the eastern lobe contains two hot-spots indicating some spreading of the jet force over a larger area than each hot-spot and a consequent slower expansion rate and larger arm-length ratio.

Now consider the effect of orientation on the X-ray emission and jet brightness ratio. As Wilson et al. (2001) remark, the X-ray emission resulting from scattering of the cosmic microwave background radiation by the electrons in the jet is proportional to $`\delta ^{4+2\alpha }`$ where $`\delta `$ is the Doppler factor and $`\alpha 0.9\pm 0.5`$ is the spectral index (Begelman & Sikora 1987; Dermer 1995). The ratio of jet to counterjet X-ray brightnesses, $`R10`$ implies that $`\beta \mathrm{cos}\theta 0.2`$ for $`\alpha =0.9`$. The radio estimate of $`R3`$ does not provide as stringent a limit – $`\beta \mathrm{cos}\theta 0.14`$.

The main constraint on the angle $`\theta `$ is that it should not be so close to $`90^{}`$ that the western jet is Doppler dimmed, in absolute terms, irrespective of the ratio of jet to counterjet fluxes. That is, $`\delta `$ should not be less than approximately $`0.7`$ for an order of magnitude decrease in the X-ray flux.

For a given ratio, $`R`$, of jet to counterjet X-ray fluxes, we have

$$\beta \mathrm{cos}\theta \frac{R^{1/p}1}{R^{1/p}+1}$$
(6)
where $`p=4+2\alpha `$ for inverse Compton emission off the microwave background. In Figure 29, we show plots of the Doppler factor against inclination for values of $`R=10,20`$ and $`50`$. An inclination of $`70^{}`$ is consistent with $`\delta 0.7`$ for $`R=10,20`$ but not for $`R=50`$. If $`\theta =70^{}`$ then Wilson et al. (2001) have shown that the magnetic field is about 0.03 times the minimum energy value. This constitutes no difficulty in principle since one does not necessarily expect minimum energy conditions to apply in the jet even though turbulent processes may be sufficient to produce minimum energy conditions in the lobes. Moreover, one possible estimate of the magnetic field in the hot-spot is an order of magnitude below equipartition. For reference, a jet inclination of $`\theta =70^{}`$ and $`R=10(20)`$ implies a jet velocity of $`\beta 0.57(0.74)`$. These estimates of velocity are consistent with other estimates in powerful sources (e.g. Wardle & Aaron 1997). We conclude therefore that an orientation of $`\theta 70^{}`$ is consistent with the X-ray and radio constraints. We note, however, that the steep decline in the Doppler factor at around $`\theta =70^{}80^{}`$ indicates that the constraints implied by the X-ray data are only just satisfied if $`\theta =70^{}`$.

A jet velocity of $`\beta 0.60.7`$ has other interesting ramifications since gas would then emerge from oblique shocks in the hot-spots with this velocity. Some of this gas would be moving directly towards us with $`\delta 3.1`$ (for $`\beta 0.6`$). The flux from that region of the hot-spot would be enhanced by a factor of 60 and this may be another explanation for the large ratio of hot-spot surface brightness to lobe surface brightness (Perley et al. 1997).

If an almost side-on orientation proves to be untenable, then an alternative is that the limitations of an axisymmetric solution do not allow for a partial ring. That is, the ring could be brighter in one part that another as a result of three dimensional hydrodynamic effects. This may be consistent with the faint structures north of the filament seen in the greyscale image of the Perley et al. (1997) radio data, the Roeser (1989) optical image (see Figure 3 of Wilson et al. 2001) and the VLT optical image (Figure 2). These features may indicate a non uniformly bright ring viewed at a more acute angle. We do not favour this explanation because of the straightness of the observed filament.

Another point in favour of $`\theta 90^{}`$ is that the western hot-spot appears, in projection, at the edge of a wide lobe. If the source were to viewed at, say $`\theta 20^{}`$, then the hot-spot would most likely appear closer to the middle of the lobe. In addition such a small viewing angle makes the total extent of the source quite large ($`1.3\mathrm{Mpc}`$). Radio galaxies this large are not unknown, but are fairly exceptional.

### 6.4 Structure functions and hot-spot ratios

The construction of cumulative distribution functions for hot-spot ratios has been incidental to the main purpose of this paper. However, it is of general interest and quantifies the role of intrinsic variability in using such ratios for the estimation of say, relativistic effects.

The calculation of structure functions in jet simulations could have wide applicability. For example, the structure function of internal shocks could be applied to blazars.

## 7 Acknowledgements

This work was supported by an Australian Research Council Large Grant, A69905341 and grants of computing time from the ANU Supercomputer Facility.

## Appendix A Probability distribution of hot-spot ratio

In each simulation the bright features in the vicinity of the hot-spot are unsteady, dynamic structures. The intrinsic variability of the peak intensity, $`I`$, of the hot-spot is easily quantified as a cumulative probability distribution function, $`G(I)`$. This function is numerically approximated by sorting the instantaneous $`I`$ values from a representative period $`t_0tt_0+\mathrm{\Delta }t`$ of the simulation, in which $`\mathrm{\Delta }tt_{\mathrm{dyn}}`$. The sample must be representative in the sense that it begins and ends during the time in which the cocoon is well established and both the head of the jet and the bow shock are not close to any of the boundaries of the computational grid. If the intensity samples $`\{I_n\}`$ are ranked in ascending order, $`_n`$, for $`n=1,\mathrm{},n_{\mathrm{max}}`$, then the cumulative distribution function of the peak intensity is

$$G(I)\mathrm{Pr}(I^{}I)=\frac{(I)}{n_{\mathrm{max}}}$$
(7)
where $`I`$ is one of the sample values. If $`I`$ is not one of the sample values but has neighbouring values are, $`I_{}<I<I_+`$, then $`G(I)`$ can be interpolated from $`G(I_{})`$ and $`G(I_+)`$. Then we may define a probability density function,

$$g(I)=\frac{dG}{dI},$$
(8)
which may be evaluated for any $`I`$ by numerical differentiation and/or interpolation as appropriate.

If the two jet and counter-jet have the same physical parameters and the respective hot-spots are separated by a distance exceeding the $`t_{\mathrm{dyn}}c`$ then the hot-spots are causally independent. The probability distribution, $`f(R)`$, of the ratio, $`R`$, of hot-spot peak intensities can be derived from the intrinsic distributions of one hot-spot.

$$f(R)\mathrm{Pr}\left(RI_2<I_1<(R+dR)I_2\right)\frac{dF}{dR}=_0^{\mathrm{}}g(RI_2)g(I_2)I_2𝑑I_2.$$
(9)
Alternatively, the cumulative distribution function $`F(R)`$ may be calculated by taking a double sum over the sample of intensity values,

$$F(R)\mathrm{Pr}\left(R^{}R\right)=\frac{1}{n_{\mathrm{max}}^{}{}_{}{}^{2}}\underset{i,j}{}H(RI_jI_i),$$
(10)
where $`H`$ is the Heavyside step function ($`H(x)=1`$ for $`x1`$ and $`H(x)=0`$ otherwise).