format: text

text:
# The Identification of Infrared Synchrotron Radiation from Cassiopeia A

## 1 Introduction

The energy distribution of $`0.110`$ GeV electrons in supernova remnants has been studied for decades through their radio synchrotron radiation. The distributions are typically power laws, with low energy spectral flattenings often seen due to absorption in the interstellar medium (Kassim, 1989), or perhaps even due to internal ionized gas in the case of Cas A (Kassim et al., 1995). The power law slopes (spectral indices) vary significantly from remnant to remnant (Green, 2001). Although it has long been accepted that diffusive shock acceleration is responsible for the GeV electrons (Bell, 1978), there are still no clear observational signatures for how this is regulated in different remnants. Spectral index variations within an individual remnant are often very small; this presents a dilemma because the asymmetries in remnant structures and dynamics should lead to large variations in spectral index if particle acceleration is due to first-order Fermi acceleration at shocks. In older remnants, such as the Cygnus Loop, spectral variations probably reflect local differences in absorption and confusion along the line of sight (Leahy & Roger, 1998). But in Cas A, the spectral variations across the remnant (Anderson & Rudnick, 1996) are most consistent with different power laws, indicating local variations in relativistic particle acceleration (Wright et al., 1999). Information about deviations from the power law shape could provide important information about the acceleration mechanisms and the possible modifications of the shocks by cosmic-ray backpressures (Berezhko & Ellison, 1999); only weak evidence for such deviations exists to date (Reynolds & Ellison, 1992).

The maximum energy to which cosmic rays can be accelerated in supernova remnants is another issue of current importance. Non-thermal X-ray synchrotron radiation from the supernova remnants SN1006 (Koyama et al., 1995; Reynolds, 1996), RX J1713.7-3946 (Koyama et al., 1997; Slane et al., 1999), Cas A (Allen et al., 1997), and RCW86 (Borkowski et al., 2001) suggest that relativistic electrons have energies up to 100 TeV. Gamma-ray detections of SN1006 (Tanimori et al., 1998) and RX J1713.7-3946 (Muraishi et al., 2000) have been interpreted as the Inverse Compton emission from these high energy electrons against the cosmic microwave background.

Although the energies of electrons that would produce synchrotron radiation in the infrared are considerably lower than those emitting X-rays and $`\gamma `$rays, they offer the possibility of a model-independent probe of electrons at much higher energies than seen in the radio. The possibility of near or mid-IR synchrotron emission was first suggested to us by R. Tuffs (private communication and Tuffs et al. (1997)), based on their impressions of the ISOCAM images (Lagage et al., 1996). The possibility was again raised by Gerardy & Fesen (2001), who noted the close similarity between their 2.2 $`\mu `$m (but not 1.2 $`\mu `$m) image and the radio images. This led us to seek verification of the infrared synchrotron origin through polarization measurements.

## 2 Observations

### 2.1 Choice of Field

The observations were made using NSFCAM on the IRTF in polarimetry mode on 1 August, 2001. Since NSFCAM has a field of view of about 75$`\mathrm{}`$ (0.3$`\mathrm{}`$ per pixel) and we would have time to observe only one field, we had to make a careful choice of location within the Cas A remnant. The 2.2 $`\mu `$m band image kindly provided by C. Gerardy and R. Fesen (Gerardy & Fesen, 2001) gave us information on the flux from Cas A at 2.2 $`\mu `$m . Our own radio images of Cas A (Anderson & Rudnick, 1996; Anderson, Keohane & Rudnick, 1995; Gotthelf et al., 2001) provided us with total intensity and polarized intensity maps that could be compared to the 2.2 $`\mu `$m band image.

Our choice of which region in Cas A to observe was influenced by a number of factors. First, there had to be a good correspondence between the radio synchrotron emission and the 2.2 $`\mu `$m emission. This was done to avoid small bright regions of emission at 2.2 $`\mu `$m that might be contaminated by emission lines unrelated to the radio synchrotron emission. Second, there had to be significant polarized intensity in the radio, which we hypothesized would correspond to similar polarization in the near infrared. Third, there had to be sufficient flux at 2.2 $`\mu `$m to make fractional polarization measurements of a few percent, a level similar to the observed radio polarization. Using these criteria, we chose a segment of Cas A’s bright ring in the northwest (azimuths of $`0to40`$ degrees) centered at RA 23h 23m 26s , Dec 58, 48$`\mathrm{}`$ 58.4$`\mathrm{}`$, (J2000), as shown in Figure 1.

### 2.2 Polarization Observations

The techniques for using NSFCAM as an imaging polarimeter are described in Jones (1997, 2000). Images are taken of the source and an off-source sky location at four positions of a half waveplate that can be rotated in the beam. This allows us to form sky subtracted intensity images at polarization position angles of 0, 45, 90 and 135 (waveplate angles of 0, 22.5, 45 and 67.5). We formed two independent estimates of the total intensity ($`I_0+I_{90}`$) and ($`I_{45}+I_{135}`$), and found that these were consistent with each other within the off-source noise errors. Photometric calibration was done using observations of HD 225023 from the Elias standards list (Elias et al., 1982). The 2.2 $`\mu `$m (K) intensity scale was set using the conversion $`0.0m_K=635`$ Jy. The final accuracy of our spectral indices is limited by the 2.2 $`\mu `$m extinction corrections, as discussed below.

From the images at four angles, we then computed Stokes images Q ($`I_0I_{90}`$) and U ($`I_{45}I_{135}`$) and the total intensity image $`\frac{1}{2}(I_0+I_{45}+I_{90}+I_{135})`$. Note that these Stokes parameter definitions differ by a factor of two from those commonly used in the radio. The details of our sequencing of source, sky and waveplate position differed somewhat from Jones (1997). For Cas A we rotated back and forth between two waveplate positions (either for Q or U) for several minutes, then repeated the same procedure on the sky position. Integrations at a single waveplate position were 10 seconds. It took about 2 seconds for the waveplate to move to a new position.

NSFCAM is not a true imaging polarimeter, but rather a camera retrofitted with some polarization optics, so we were limited by systematics in the precision of our polarimetry to about $`\pm `$0.3%. This limit is due primarily to fluctuations in atmospheric transmission between the time it takes to integrate at different positions of the waveplate. Since we are interested in polarizations of a few percent or more in Cas A, this instrumental limitation did not hinder our observations. In fact, the nebulosity in Cas A at 2.2 $`\mu `$m is sufficiently faint that our measurements are effectively photon noise limited by the background.

Calibration of the position angle was done using observations of AFGL 490 two nights earlier with the same observing procedure as for Cas A. AFGL 490 was assumed to have a polarization angle of $`\theta `$ = 115 at 2.2 $`\mu `$m (Kobayashi et al., 1978). The Q’ and U’ values quoted in Table 1 are given as observed; the PA contains a $`33`$ degree correction to place it in the proper sky frame. An extra check on the position angle calibration was done using observations on 1 August of the reflection nebulosity surrounding the young star RNO1 (Weintraub & Kastner, 1993). The centro-symmetric scattering pattern around this star produces all possible position angles and provides a very sensitive calibration. Calibration of the efficiency (which exceeds 98% at 2.2 $`\mu `$m ) was performed in previous observing runs using S1 in $`\rho `$ Oph which has a polarization of 1.9% at 2.2 $`\mu `$m . Instrumental polarization was checked during the night by observing an unpolarized standard star from the UKIRT faint standards list (FS 149) and was found to be too small to be measured.

The telescope was offset 100$`\mathrm{}`$ to the North for the sky frames. This was intended to move the telescope well off Cas A, but not so far that telescope movement and settling time seriously reduced our observing efficiency. In an effort to increase our time on Cas A itself, we reduced the number of background frames to about half the number of frames taken on the source. The total on source integration time was $`1`$ hr. Serious telescope drift during the long integrations on Cas A resulted in the movement of stars in the background frames relative to the source frame. This produced elongated stars in the background subtracted images (see Figure 1, 2.2 $`\mu `$m intensity). No correction has been made for these problems, and the effect on our spectral measurements is described below. Our polarization measurements are less affected by the star trails than the total intensity because the star polarizations are expected to be $`1\%`$ as discussed below, while Cas A’s polarization is $`510\%`$.

Since the time between source and background frames ranged over several minutes, fluctuations in sky background could cause a residual DC level to be present in the sky subtracted images. Although small compared to the net intensity of the brighter nebulosity in Cas A, this DC offset could be comparable to the polarized intensity, and produce a spurious net intensity in the Q and U images. These offsets were removed by first identifying several small regions of the sky-subtracted total intensity image with values fluctuating near zero for the intensity. We simply assume the true intensity is zero when averaged over these locations and subtract the respective means from the I, Q and U images to force them to be zero there. This procedure will remove the DC offset produced by sky background fluctuations in the Q and U images, but also would remove any real, but faint emission. Sensitive observations in the radio do show faint total intensity and polarized emission from these regions, at intensities of approximately 15% that of the bright ring; no correction has been made for this small bias.

### 2.3 Radio and Optical Images

The radio images were produced from observations at the Very Large Array 111The Very Large Array is an instrument of the National Radio Astronomy Observatory, a National Science Foundation facility operated under cooperative agreement by Associated Universities, Inc. using all four configurations over the time period from May, 1997 through March, 1998. Maps were made by using the data from 1.285 GHz for the 20 cm band, and from 4.64 GHz for the 6 cm band. The flux calibration scale was set using standard VLA procedures (VLA Calibration Manual, 2002) and the source 3C48. The polarization calibration also used standard VLA procedures with the angle fixed by assuming a position angle of $``$10 for 3C138. Other details of analysis and image creation are similar to those described by Anderson, Keohane & Rudnick (1995). The full resolution of these images was $`1.5\mathrm{}`$.

The optical emission in Cas A’s bright ring is dominated by clumpy emission lines on a wide range of scales, most recently studied in great detail by Fesen et al. (2001). The emission line material represents a distinct temperature component, separate from the synchrotron plasma, but which may dominate the dynamics driving the relativistic particles and magnetic fields. In particular, Fesen et al. (2001) noted from their high resolution WFPC2 observations, a very large number of filamentary structures that they suggested were Rayleigh-Taylor instabilities. Although there is significant small-scale scatter, we measured the approximate position angle of the filaments at the azimuths of the boxes used for the synchrotron measurements and list these in Table 1.

## 3 Data Analysis & Results

### 3.1 Image Registration

The registration of the 2.2 $`\mu `$m and radio images was done in two steps. The radio coordinates from the VLA are already in the FK4 reference frame to an accuracy of $`0.1\mathrm{}`$. To put the 2.2 $`\mu `$m image in the same frame, we first used the Sloan Digital Sky Survey to register the image of the entire Cas A field kindly provided by C. Gerardy & R. Fesen (Gerardy & Fesen, 2001). We then minimized the rms differences between our new small polarimetry field and the Gerardy/Fesen image, allowing for shifts in position, rotation, and scale. The final accuracy of this registration is $`<0.1\mathrm{}`$, and the uncertainties have no effect on the results discussed here.

### 3.2  Spatial averaging

The signal:noise on the $`0.3\mathrm{}`$ pixel level was not sufficient either for accurate spectral measurements or determination of the polarization. We therefore found the average brightness for I, Q, and U in each waveband in six boxes as shown in Figure 1, and listed in Table 1. These boxes were chosen to have sufficiently high polarized flux at 6 cm, to avoid regions with bright stars in the sky (“off”) frame, and to avoid the edges of the frame where the flat fielding is less accurate.

For display purposes, we also convolved the polarizations (Q and U) to a resolution of 10$`\mathrm{}`$. The resulting low resolution polarized intensities are shown on the right side of Figure 1. The middle image of the right side represents a prediction of the 2.2 $`\mu `$m polarizations using the radio data and 2.2 $`\mu `$m total intensities. It was constructed by taking the fractional polarization image at 6 cm and multiplying it by the total intensity image at 2.2 $`\mu `$m . The overall correspondence is good, but not perfect. The 2.2 $`\mu `$m polarized intensity map contains, e.g., a bright spot due to residuals from the bright star immediately to the SE of box 6. In addition, it has not been corrected for the noise bias, which is responsible for some of the structure in the image. The 6 cm polarization angles were used, after correction for Faraday rotation as described below, as a prediction for the 2.2 $`\mu `$m angles. The correspondence here is somewhat better than expected from the errors.

A number of corrections to the total and polarized intensities are needed before the 2.2 $`\mu `$m and radio intensities and polarizations can be quantitatively compared. First is the contamination of the 2.2 $`\mu `$m total intensity by foreground stars in the boxes. The fluxes of the brightest stars in each box were therefore subtracted yielding the “corrected” 2.2 $`\mu `$m intensities given in Table 1. We estimate that the remaining positive star contamination is less than a few percent.

More troublesome are the negative trails from stars in the sky frame, since the contribution from each individual star is spread across the box, and is therefore fainter and harder to recognize. Although the boxes avoid the brightest negative trails, there are certainly unrecognized trails that reduce the 2.2 $`\mu `$m intensities. In individual regions, such as the NW corner of the frame, we can look for negative trails in the absence of bright diffuse emission from Cas A. The distribution of pixel values in that region is consistent with a gaussian distribution around mean zero with a long negative tail. For this region, the negative tail can be recognized and ignored in determining the zero offset level. When negative star trails overlap bright 2.2 $`\mu `$m emission, however, they can reduce the intensity without being recognized. There is thus a negative bias in the 2.2 $`\mu `$m intensities, which we expect to be small on average ($`0.5\mu `$Jy/$`\mathrm{}\mathrm{}`$, or $`5\%`$ of the total intensity), but could possibly be as much as $`30\%`$ in an individual box, without being recognized. The 2.2 $`\mu `$m total intensity values quoted in Table 1 could therefore actually be somewhat larger, but because of the scientific significance of such larger values, we have taken a more conservative approach and not tried to make any corrections for this possible bias.

Another bias correction that we did not make to the data is due to the fact that the sensitive radio observations show that synchrotron emission is present everywhere in the 2.2 $`\mu `$m image, and therefore it is not really possible to define an appropriate zero, as discussed earlier. The characteristic intensity of the radio plateau is approximately 15% of the bright ring, so defining it to be zero at 2.2 $`\mu `$m results in spectral indices between 2.2 $`\mu `$m and 6 cm that are too steep by $`0.01`$.

Cas A’s high optical extinction ($`>4^m`$ Searle (1971) ) requires that corrections are needed even at 2.2 $`\mu `$m . Hurford & Fesen (1996) find reddenings of $`4.6<A_V<5.4`$ magnitudes for five FMKs and $`A_V`$ = 5.3 and 6.2 magnitudes for two QSFs, and provide an extensive discussion of the literature on Cas A’s extinction. In the west, $`A_V`$ may extend up to 8 magnitudes, according to the CO and HI analysis of Troland, Crutcher & Heiles (1985). In this paper, we use a value of $`A_V`$ = 5 magnitudes, which leads to $`A_{2.2\mu m}`$ = 0.55 magnitudes, adopting the conversion factors of Bessell, Castelli, & Plez (1998). Application of this extinction value results in the lines labeled “de-reddened” for I, Q, and U in Table 1, and the values plotted in Figure 2. We also show, in Figure 2, the range of corrections that would result for the range $`4.6<A_V<6.2`$ magnitudes.

For any reasonable value of the extinction, it thus appears that the 2.2 $`\mu `$m data consistently fall above the extrapolation from radio wavelengths, giving the spectra a concave shape, i.e., flatter at short wavelengths. There are several caveats to this conclusion. First, a mixture of different power-laws, as observed at cm wavelengths (Anderson & Rudnick, 1996; Wright et al., 1999), will result in a concave spectrum. The rms variations within our boxes are $`<=0.01`$, insufficient to cause a significant effect at 2.2 $`\mu `$m . However, it is likely that there is some contamination of the bright ring with emission from the steeper spectrum plateau, at least at radio wavelengths. The maximum contamination, for plateau emission with $`\alpha _{6cm}^{20cm}=0.9`$ would result in an increase of the 2.2 $`\mu `$m flux density of $`35\%`$. The actual 2.2 $`\mu `$m contamination is likely to be considerably lower, given our procedure for setting a zero background, as discussed earlier. However, in any case, this possible contamination is much less than the factor of two enhancement seen at 2.2 $`\mu `$m after correction for extinction.

A second concern is our calibration to the radio flux scale (VLA Calibration Manual, 2002), which itself is based on “absolute” calibrations of Cas A (Baars et al., 1977). These two steps are expected to be accurate to $`2\%`$, and so could each contribute another $`25\%`$ error to the 2.2 $`\mu `$m extrapolation. We therefore consider our detection of spectral concavity to be reasonable, but somewhat tentative until the extinction and flux scale issues can be fully resolved. At present, we are analyzing the 6 cm $``$ 2.2 $`\mu `$m spectra across the whole remnant, which should allow us to use differential measurements to separate real from calibration effects.

### 3.3 Polarization

The 2.2 $`\mu `$m Stokes parameter values in Table 1 are listed as Q, U; these represent the actual observations, which are only nominally Q and U prior to the calibration of the absolute polarization angle. The properly calibrated polarization position angles are listed in Table 1.

In either the (Q,U) or (Q,U) frame, the observed polarized intensity ($`\sqrt{(}Q^2+U^2)`$) represents a biased measure since it would give a positive result from noise, even in the absence of a true signal. As a simple correction for this bias, the fractional polarizations for 2.2 $`\mu `$m listed in Table 1 were calculated as $`\%P=[\sqrt{(}Q^2\sigma _Q^2+U^2\sigma _U^2)]/I`$ where $`\sigma _{Q,U}`$ were calculated as the rms scatter among boxes of approximately 35 pixels on a side, away from bright stars and away from the bright ring.

For the 2.2 $`\mu `$m polarization measurements, we also need to consider the possible contamination by the foreground interstellar polarization expected from the large extinction towards Cas A. Catalogs of optical polarimetry (e.g. Heiles (2000)) have no stars in the same field as Cas A, so a direct measurement of the interstellar polarization in the visual has not been made. We have measured the polarization of the bright star in the northwest of our image and find a 3$`\sigma `$ upper limit of 0.9% at 2.2 $`\mu `$m , although we do not know if it is as far away as Cas A. The expected 2.2 $`\mu `$m polarization for Cas A, given the visual extinction of 5-6 magnitudes, can be estimated using the mean trend in polarization with extinction analyzed in Jones (1989) and Jones, Klebe & Dickey (1991) of

$$<P_{2.2\mu m}(\%)>=2.2\tau _{2.2\mu m}=0.2A_V1\%.$$

Thus, we conclude that interstellar polarization is not a significant contaminant of our polarimetry of the nebulosity at 2.2 $`\mu `$m . After consideration of all these effects, and the faintness of the nebulosity being measured, the errors in the fractional polarizations and angles at 2.2 $`\mu `$m are dominated by the random contribution from sky background photons, as opposed to systematic effects. The measured percentage polarizations at 2.2 $`\mu `$m range from $`410`$ %, comparable to those at 6 cm.

To compare the 2.2 $`\mu `$m polarization angles to those at radio wavelengths, we must first correct for Faraday Rotation, which can arise both local to Cas A and along the line of sight. At 20 cm, the rotation is strong enough to actually depolarize the emission (Anderson, Keohane & Rudnick, 1995), and render rotation measure (RM) measurements at that wavelength unreliable for the higher frequencies. Kenney & Dent (1985) present RM measurements at low resolution but high frequencies, with a value of $`106\mathrm{rad}\mathrm{m}^2`$ for the northwestern part of the remnant. In order to look for variations in RM among the boxes used here, we created three different polarization maps at 4.41, 4.64 and 4.99 GHz, using the VLA D configuration data alone to produce a beam size of $`10\mathrm{}`$. The signal:noise was insufficient to detect significant variations among the boxes, but averaging over all boxes yields $`RM=105\pm 3\mathrm{rad}\mathrm{m}^2`$, consistent with the Kenney & Dent (1985) value. We therefore adopt the value of $`106\mathrm{rad}\mathrm{m}^2`$ for all boxes, and show the Faraday corrected position angles in Figure 1 and in the angle comparisons at the bottom of Table 1.

In Figure 3, we compare the position angles of the magnetic fields derived from the polarization vectors in each box with the approximate angles of the filamentary structures in the WFPC2 images of Fesen et al. (2001). All of the position angle measures change systematically with azimuth, with a small but significant clockwise offset from the local radial direction at least for the 2.2 $`\mu `$m and filament angles. It is not clear if the 6 cm angles show any offset, but better local rotation measures would be needed to make the comparison with 2.2 $`\mu `$m angles more accurately.

## 4 Scientific Implications

### 4.1 Maximum energy of accelerated electrons

Any mechanism to accelerate relativistic particles yields some maximum particle energy. The energy can be limited by a variety of factors including the time available for acceleration, the escape of the highest energy particles from the finite sized acceleration region, or the balance between energy gains and losses (Reynolds, 1998, 2001). In Cas A, a hard X-ray tail, which falls two orders of magnitude below the flux extrapolated from radio wavelengths, suggests that the synchrotron emission may extend to that regime (Allen et al., 1997). This conclusion has been called into question by Bleeker et al. (2001), who argue that the 10 $``$ 15 $`keV`$ emission is diffuse, and not concentrated at the outer shock as it should be if it were synchrotron. This argument is questionable, however, because the radio synchrotron radiation also does not show an enhancement at the location of the outer shock; Xray synchrotron emission is therefore still a possibility. Bremsstrahlung provides an alternative explanation for the hard X-ray emission (Laming, 2001), where lower hybrid waves in Cas A’s strong magnetic field scatter electrons up to energies of tens of keV. There is not yet any good test to distinguish between bremsstrahlung and synchrotron explanations for Cas A, so the nature of the 10 - 100 $`keV`$ emission is therefore still open to question. With the demonstration of a synchrotron origin for the 2.2 $`\mu `$m emission, it is therefore interesting to estimate the maximum energy of the accelerated particles.

We start by calculating the value of the field strength that minimizes the energy in both relativistic particles and magnetic field ($`B_{min}`$). We assume a pathlength through the emitting regions of 0.5 pc, with a filling factor of unity. We assume no energy in relativistic protons. These assumptions lead to a field strength of $`500\mu `$G for the regions in Table 1. This field strength can then be used to estimate the energy of the relativistic electrons that are radiating at 2.2 $`\mu `$m . In a $`500\mu `$G field, this leads to energies of $`150`$ GeV, or a relativistic gamma of $`3\times 10^5`$. As argued below, the synchrotron spectrum almost certainly extends a factor of 10 higher in frequency, which would result in particle energies(gammas) of $`450`$ GeV ($`9\times 10^5`$). This is still a factor of at least $`10^3`$ below energies to which SNRs are commonly expected to accelerate cosmic rays, but it is the first unambiguous demonstration of the extension of the same electron population to such high energies in a shell SNR.

### 4.2 Synchrotron loss limits

We use the same field strength calculation, and the lack of spectral steepening between 20 cm and 2.2 $`\mu `$m to place a limit on the time since the relativistic particles were last accelerated. The assumptions in the field strength calculation bias the lifetimes to longer values, so our lifetime estimates should be considered upper limits. At a factor of 10 in frequency below the “break” frequency, the spectrum should already have steepened by 0.1 (e.g., Leahy (1991), Figure 3.3a). Since, if anything, the spectra have actually flattened at high frequencies, we adopt 0.1 as a conservative upper limit to the steepening at 2.2 $`\mu `$m , yielding a “break” wavelength of $`<0.22\mu `$m. This leads to lifetimes for the electrons radiating in the infrared that are no longer than 80 years. It is thus likely that particle acceleration is an ongoing process, at the current stage of Cas A’s evolution. Based on our latest dynamical picture, this is the period when Cas A has probably swept up $`110`$ times its own mass (Gotthelf et al., 2001; Willingale et al., 2002; DeLaney & Rudnick, 2002) and during which the reverse shock is slowly coming to the end of its outward motion (in the sky frame). The reverse shock is still strong at this stage, and is encountering outward moving ejecta at speeds of $`50006000`$ km/s (Reed et al., 1995). The x-ray emitting filaments are rapidly decelerated to $`3500`$ km/s while those in the radio are slowed to $`1800`$ km/s (Koralesky et al., 1998; Vink et al., 1998; DeLaney & Rudnick, 2002); this provides at least one possible source of energy for the relativistic particle acceleration.

### 4.3 Non-linear shocks

If the acceleration of cosmic rays at SNR shocks is efficient, then there is a back reaction of the cosmic rays on the shock strength and structure (Drury & Völk, 1991). When the maximum cosmic ray momentum is limited by shock properties such as geometry, the compression ratio can become arbitrarily large (Eichler, 1984). The increasing gyro-radii at higher energies allows those particles to see a higher compression due to the sub-shock structure, so the resulting spectrum becomes concave (Bell, 1987; Berezhko & Ellison, 1999; Reynolds & Ellison, 1992), as observed here. The radio spectra can then be much steeper than the $``$0.5 expected in the strong shock test particle limit, and can flatten out considerably before they cut off, e.g., in the X-ray, due to radiative losses and escape. Reynolds & Ellison (1992) suggested that slight curvature might have been observed in the integrated radio spectra of Tycho and Kepler SNRs, although the signal was marginal.

In this work, we have shown that the 2.2 $`\mu `$m data, after correction for extinction, fall significantly above the power law extrapolation from the radio. The mean flattening between $`\alpha _{20cm}^{6cm}`$ and $`\alpha _{6cm}^{2.2\mu m}`$ is 0.06 (Table 1), equivalent to higher fluxes at 2.2 $`\mu `$m by a factor of $`2`$ . If the curvature were constant across this full range of wavelengths, then the local spectral index would change from $`0.75`$ at cm wavelengths to $`0.5`$ at 2.2 $`\mu `$m . This is essentially the same as the model shown by Berezhko & Ellison (1999). Using their Figure 5, we made a rough calculation of the flux expected at 2.2 $`\mu `$m for Box 1; this is shown as the nested square symbol in Figure 2. The agreement is quite good, suggesting that the signature of acceleration at modified shocks has been observed. In the context of the Berezhko & Ellison (1999) model, this implies a critical injection rate of relativistic particles of $`\stackrel{>}{}10^4`$ and, for most of the parameters they explored, efficiencies of 10 $``$ 100 %. Thus, with the caveats described earlier regarding the extrapolation of the radio spectrum, relativistic particles appear to be an important part of the dynamics of the shock and the overall evolution of the SNR.

### 4.4 Polarization

The fractional polarizations and the polarization angles are very similar at 6 cm and 2.2 $`\mu `$m . This provides confirming evidence not only for synchrotron radiation, but that we are looking at the same population of electrons at both wavelengths, validating the detailed spectral comparisons above. The 2.2 $`\mu `$m results confirm the expectation from the Kenney & Dent (1985) mm wave studies that fractional polarizations of $`510\%`$ represent real field disorder, as opposed to Faraday depolarization. Thus, although we speak about the “radial” magnetic fields in Cas A, this only refers to the small net ordered component field. Approximately half of the field energy must be in a component which is disordered on scales $`<1\mathrm{}`$ ($`<0.02`$ pc at a distance for Cas A of 3.4 kpc).

Gull (1975) first suggested that Rayleigh-Taylor instabilities in the decelerating bright ring of Cas A and other young SNRs would lead to radial magnetic fields. Detailed numerical simulations since that time confirm the prevalence of these filamentary structures and their role in the amplification and alignment of the magnetic fields (Jun & Norman, 1996; Jun & Jones, 1999). Approximately radial filamentary structures have been seen in the Crab SNR, and, more recently in the WFPC2 images of Cas A (Fesen et al., 2001). In the northwest region investigated here and shown in Fesen et al. (2001) Figure 8, the Cas A filaments emerge from the suggested location of the reverse shock in directions that are systematically non-radial. These non-radial motions could arise from large scale asymmetries in the ejecta, as is well established from the distribution of different elements seen in the X-ray (e.g. Hwang, Holt & Petre (2000)). Alternatively, gradients in the circumstellar medium could cause asymmetric deceleration (e.g., Koralesky et al. (1998)), again leading to large scale velocity shears. In any case, the close correspondence between the R-T finger orientations and the local magnetic field directions strongly supports a connection between the two. Blondin & Ellison (2001) show that the R-T instabilities generated at the reverse shock can extend all the way out to the outer shock front, providing radial magnetic fields instead of the tangential ones otherwise expected at the shock. The transition region between the radial and tangential regions may have been observed in Cas A (Gotthelf et al., 2001), although the signal-to-noise is quite low.

## 5 Conclusions

The extension of Cassiopeia A’s known synchrotron spectrum to 2.2 $`\mu `$m opens up a variety of important issues regarding the acceleration of relativistic particles. First, we have established that the acceleration mechanism, usually assumed to be diffusive acceleration at shocks, operates to energies of at least 150 GeV in Cas A. In addition, there is reasonable evidence that these shocks are efficient particle accelerators, arising from non-linear interactions between the particles and the shocks. This creates the opportunity to investigate these interactions in conjunction with the dynamical information about the several coupled plasmas (thermal and non-thermal) in the remnant. In turn, this will strongly constrain models for the high X-ray energy tail in Cas A’s spectrum and assist in the challenge of isolating the energy in relativistic protons. To make progress on these issues, a detailed investigation of the radio/2.2 $`\mu `$m spectral shape as a function of the radio spectral index variations across the remnant is underway. The association of magnetic fields with the local R-T instabilities opens the door for investigating the accompanying amplification of magnetic fields in the remnant, and the roles of both relativistic particles and fields in the further dynamical evolution of the remnant.

We appreciate the 2.2 $`\mu `$m images provided by C. Gerardy and R. Fesen, and a number of important discussions on the science with Tom Jones. Support for work on Cas A at the University of Minnesota was provided by National Science Foundation grant AST-96-19438 and for comparisons with HST data by NASA through grant number HST-AR-090537.01-A from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. The comments of the referee were very helpful in improving the paper.


text:
# 3D convection simulations of the outer layers of the Sun using realistic physics.

## 1 Introduction

It is now just possible to perform physically realistic 3D simulations of the surface layers of the Sun which take into account the complex interaction between radiative and convective energy transports (Stein & Nordlund 1998, Kim & Chan 1998 hereafter denoted SN and KC, respectively, Stein & Nordlund 2000). There have also been a number of 2D simulations of the surface layers (Steffen et al. 1990 and Gadun et al. 2000). To model solar convection realistically, requires a realistic equation of state, realistic opacities and proper treatment of radiative transfer in the shallow layers. SN and KC are the two most freqently cited 3D models of this type of stratified convection. As their approaches differ considerably, both in numerical methods and in input model physics, it is important to find out, what particular aspect of the respective simulations, caused their results to be different.

The aim of this paper is to describe solar surface convection which not only has the KC realistic physics, but also has a realistic (as is presently possible) geometry. To achieve this we had to increase both the depth and width of the original KC model, until the side walls and the horizontal boundaries had only a minimal effect on the flow. The simulations themselves model a region less than a few thousand kilometers in depth, as measured inwards from the visible solar surface. In the deep regions of the solar convection zone, the turbulent velocity is subsonic and the superadiabaticity is close to zero. However, within a few hundred kilometers of the solar surface, the convective flux starts to decrease. As the total flux is constant, the radiative flux must increase to offset the drop in the convective flux. This is achieved by a rise in the local temperature gradient $``$. This region of inefficient convection is called the superadiabatic layer (SAL). In the SAL, the superadiabaticity $`_{\mathrm{ad}}`$ is positive and of order unity (Demarque, Guenther & Kim 1997, 1999). As the buoyancy force is large in the SAL, the region is characterised by highly turbulent velocities and large relative thermodynamic fluctuations. The turbulent velocity also gives rise to a significant turbulent pressure (approximately 15 % of the gas pressure). This moves out the convection surface, modifying the SAL and the stratification. In one-dimensional (1D) models of the solar convection zone based on the mixing length theory (MLT)(Böhm-Vitense 1958), the velocity is set to zero above the convection boundary. However, three-dimensional (3D) numerical simulations described in Cattaneo et al. (1990), have shown that just above the top of the convection layer the turbulent velocities are still high.

There are several motivations for such simulations among stellar physicists. One is to understand the effects of turbulence on the structure of the outer solar layers as revealed by the observed frequencies of solar p-modes. Another is to explain the excitation mechanism of the p-modes. Still another is to derive more realistic surface boundaries for stars with convection zones from first physical principles, free of the arbitrary assumptions of the mixing length theory. And finally, such simulations may be of help in investigations of the solar dynamo.

Over the last few years, the science of helioseismology has provided some very precise measurements of the p-mode oscillation frequencies (to within one part in a thousand) in the surface layers of the Sun \[Harvey et al. 1996\]. The discrepancy between the observed p-mode frequencies and those calculated from solar models, is known to be primarily due to the approximations made in modelling the surface layers, where turbulent and radiative losses are significant (Balmforth 1992; Guenther1994). In the 3D simulation described in KC, the turbulent pressure pushed the convective boundary radially outwards from its original position (that was computed using the MLT). This situation was mimicked in the 1D solar model by tweaking the opacity in the outer layers \[Demarque, Guenther & Kim 1999\]. This resulted in improved p-mode frequencies for low and intermediate degrees. However, Demarque et al. also showed that the mixing length prescription of Canuto & Mazzitelli (1991), which had a completely different SAL structure than KC, could produce a similar improvement of the p-mode frequencies for the same $`l`$-values. Full details of the different approaches, the contrasting SAL structures and the resulting p-modes, are described in Demarque et al. (1997, 1999). Later, Rosenthal et al (1999) used another approach to compute the p-mode frequencies. These authors patched the mean stratification (horizontal average) of a 3D hydrodynamical simulation, onto a 1D MLT envelope model. In order to get a smooth 1D model, they adjusted the mixing length and the amplitude of the turbulent pressure to match the 3D simulation. The computed frequencies were closer to the observed frequencies than if a standard solar model is used, thus showing the importance of including turbulence in modeling the outer layers of the Sun. More recently, Li et al. (2002) found similar results to Rosenthal et al. by inserting the averaged turbulent pressure and turbulent kinetic energy directly into the 1D stellar model. Their method is applied to two of the simulations described in this paper (see section 4.3). As the turbulent kinetic energy and turbulent pressure were very small at the base of the 3D model, they were set to zero in regions of the 1D stellar model that lay below the 3D model domain. This required the usual adjustment of the mixing length parameter and the helium abundance to calibrate the perturbed stellar model. No other adjustable parameters were employed. The improvement in the eigenfrequencies was found to be primarily due to the inclusion of turbulent kinetic energy flux in the 1D stellar model.

After extensive testing, we found that our simulations are in good agreement with other numerical studies of the surface layers (e.g. Rosenthal et al. 1999, Asplund et al. 2000). In addition, by incorporating the computed 3D turbulence into a 1D stellar model, we were able to produce solar surface eigenfrequencies (p-modes) that were very close to the observed frequencies. As our eventual goal is to simulate the SAL in stars other than the Sun, it is vital to be sure we are modelling the Sun as correctly as possible.

## 2 Modelling realistic solar surface convection

To model surface layer convection in the Sun as realistically as possible, we take the following approach:

1. Using a stellar evolution code (YREC e.g. see Guenther et al. (1992)), we compute a standard stellar model from which the initial density $`\rho `$ and internal energy $`e`$, required by the 3D simulations, are derived. From an arbitrary initial velocity field $`𝐯`$ , we then compute $`\rho `$, $`E(=1/2\rho v^2+e)`$, $`\rho v_x`$, $`\rho v_y`$ and $`\rho v_z`$. These are the dependent variables of the governing equations. The horizontal directions are $`x`$ and $`y`$, and $`z`$ is radially outwards.
2. Using the same tables for the equation of state and the opacities as in the stellar model, we then compute the pressure $`P(\rho ,e)`$, temperature $`T(\rho ,e)`$, Rosseland mean opacity $`\kappa (\rho ,e)`$, specific heat capacity at constant pressure $`c_\mathrm{p}(\rho ,e)`$, adiabatic gradient $`_{\mathrm{ad}}(\rho ,e)`$ and some thermodynamic derivatives.
3. The radiation flux is then computed using the diffusion approximation in the optically thick regions and the 3D Eddington approximation in the optically thin layers.
4. We then integrate the Navier-Stokes equations over one time step to compute a new set of dependent variables and return to (ii).

### 2.1 Realistic initial conditions: steps (i) and (ii)

The solar model uses the same realistic physics as described in Guenther & Demarque (1997). In particular, the low temperature opacities of Alexander & Ferguson (1994) and the OPAL opacities and equation of state were used \[Iglesias & Rogers 1996\]. Hydrogen and helium ionisation, and the diffusion of both helium and heavy elements are included.

In the initial model, the atmospheric layers, which are in radiative equilibrium, are assumed to be gray. Deeper in, in the convectively unstable region, the thermal structure is described by the MLT, which prescribes the temperature gradient $``$. In the original KC simulation, the Eddington approximation T($`\tau `$) relation was used in the atmosphere. In this case, the values of the parameters, X, Z and $`\alpha `$, in the calibrated Standard Solar Model (SSM) are $`(X,Z,\alpha )=(0.7385,0.0181,2.02)`$, where X and Y are the hydrogen and helium abundances by mass, and $`\alpha `$ is the ratio of mixing length to pressure scale height in the convection zone, required to match precisely the solar radius. It should be mentioned that although the values of X and Y in the calibrated SSM depend little on the treatment of the surface layers, the value of $`\alpha `$ is sensitive to the choice of T($`\tau `$) relation in the atmosphere (Guenther et al. 1992).

We note that for simulations discussed in this paper (those listed in the appendix), the empirical Krishna-Swamy $`T(\tau )`$ relation for the Sun \[Krishna Swamy (1966)\] was used instead of the Eddington approximation. The same model as KC, but constructed with the Krishna-Swamy T($`\tau `$) relation, yields $`(X,Z,\alpha )=(0.7424,0.01706,2.1319)`$ (this is model KC2 in the appendix). Because of the inclusion of helium and heavy element diffusion during the evolution, the initial X and Z (denoted by X0,Z0) were slightly different. For KC they were (X0,Z0)=(0.7066, 0.0201), and for KC2 they were (X0,Z0)=(0.710,0.019). The difference in initial atmospheric structures between simulations KC and KC2 provides us with the opportunity to verify that that the final state in the simulations is not strictly dependent on the choice of T($`\tau `$) relation in the initial model. This important test is discussed in the appendix.

For each timestep of the numerical integration we need to solve the complex equation of state (step (ii)). This adds considerable computational time compared to an ideal simulation, such as with a perfect gas. The 3D hydrodynamical simulations use identical opacities and equation of state as used in the 1D reference SSM, which served as the initial model. The side boundaries are periodic, while the top and bottom boundaries are stress free. A constant heat flux flows through the base, and the top is a perfect conductor. To ensure that mass, momentum and energy are fully conserved, we use impenetrable (closed) top and bottom boundaries.

A particular model is specified by $`g`$ (the surface gravity) and $`T_{\mathrm{eff}}`$ (the effective temperature). Aside from the viscosity coefficients there are no other free parameters.

### 2.2 Radiative transfer: step (iii)

In the SAL, the photon mean free path may not be small enough to use the diffusion approximation. Consequently one is forced to either solve the full radiative transport equation or consider the three-dimensional Eddington approximation (Unno & Spiegel 1966), which is a higher order approximation than the diffusion approximation, and is valid in the optically thin regions. Computationally, to solve the full radiative transport equations is formidable, and only a few ray directions are currently used in this approach (SN). In our simulations, we have chosen to use the three-dimensional Eddington approximation.

In the deeper part of the domain ($`\tau >10^4`$), we use the diffusion approximation,

$$Q_{\mathrm{rad}}=\left[\frac{4acT^3}{3\kappa \rho }T\right],$$
(1)
where $`\kappa `$ is the Rosseland mean opacity, $`a`$ is the Boltzmann constant and $`c`$ is the speed of light. In the shallow region $`Q_{\mathrm{rad}}`$ is computed as

$$Q_{\mathrm{rad}}=4\kappa \rho (JB)$$
(2)
where the mean intensity $`J`$ was computed by using the generalized three-dimensional Eddington approximation (Unno & Spiegel),

$$\left(\frac{1}{3\kappa \rho }J\right)\kappa \rho J+\kappa \rho B=0,$$
(3)
where $`B`$ is the Planck function. This formulation is exact for isotropic radiation in a gray atmosphere, and without requiring local thermodynamic equilibrium, the Eddington approximation describes the optically thick and thin regions exactly (Rutten 1995). To study spectral line profiles and spectral energy distribution, requires frequency dependent radiative transfer. However, in a study of the SAL, to serve as a surface boundary condition for stellar models or for comparison with results of helioseismology, we found a gray atmosphere to be adequate.

### 2.3 Hydrodynamics: step (iv)

For deep ($`v<<c_\mathrm{s}`$, where $`v`$ is the flow velocity and $`c_\mathrm{s}`$ is the isothermal sound speed) and efficient ($`_{\mathrm{ad}}`$ just above zero) convection, Chan & Sofia (1989) showed that the MLT is a very good approximation to the real situation. However, in the SAL, both $`v/c_\mathrm{s}`$ and $`_{\mathrm{ad}}`$, can be of order unity. In this case, the MLT is unlikely to apply. The validity of this argument is confirmed by helioseismology. The run of the sound speed in the SAL derived by inversion of the helioseismic data does in fact disagree with the MLT model (Basu & Antia 1997). In such an environment, the governing hydrodynamic equations are the fully compressible Navier-Stokes equations (see for example Kim et al. 1995).

$`\rho /t`$ $`=`$ $`\rho 𝐯`$ (4)
$`\rho 𝐯/t`$ $`=`$ $`\rho \mathrm{𝐯𝐯}P+𝚺\text{ }+\rho 𝐠`$ (5)
$`E/t`$ $`=`$ $`[(E+P)𝐯𝐯𝚺\text{ }+𝐟]`$
$`+\rho 𝐯𝐠+Q_{\mathrm{rad}}`$
where $`E=e+\rho v^2/2`$ is the total energy density and $`\rho ,𝐯,P,e`$ and $`𝐠`$, are the density, velocity, pressure, specific internal energy and acceleration due to gravity, respectively. $`Q_{\mathrm{rad}}`$ is the energy transferred by radiation (see previous section) and $`𝐟`$ is the diffusive flux. Ignoring the coefficient of bulk viscosity, the viscous stress tensor for a Newtonian fluid is $`\mathrm{\Sigma }_{ij}=\mu (v_i/x_j+v_j/x_i)2\mu /3(𝐯)\delta _{ij}`$, where $`\mu `$ is the dynamic viscosity and $`\delta _{ij}`$ is the delta function.

In fully developed turbulence the ratio of the length of the largest eddy to the dissipation length is $`\mathrm{Re}^{3/4}`$, where Re is the Reynolds number, (Landau & Lifshitz 1987). In the Sun, Re is the order of $`10^{12}`$ which means about $`10^9`$ scales per dimension. A 3D direct numerical simulation of the Sun would thus require about $`10^{27}`$ grid points !

The Large Eddy Simulation (LES) approach assumes that the small scales are independent of the resolved scales (large eddies) and can be parameterised as a diffusion process. In this case $`\mu `$ is an eddy viscosity defined in terms of the resolved velocity (Smagorinsky 1963),

$$\mu =\rho (c_\mu \mathrm{\Delta })^2(2𝝈\mathbf{:}𝝈)^{1/2}.$$
(7)
The colon inside the brackets denotes tensor contraction of the rate of strain tensor $`\sigma _{ij}=(_iv_j+_jv_i)/2`$. The sub-grid scale (SGS) eddy coefficient $`c_\mu `$, is set to 0.2, the value for incompressible turbulence, and $`\mathrm{\Delta }=(\mathrm{\Delta }_x\mathrm{\Delta }_y)^{1/2}\mathrm{\Delta }_z`$ is an estimate of the local mesh size. To handle shocks, $`\mu `$ is multiplied by $`1+C(𝐯)^2`$, where the constant $`C`$ is made as small as possible, while still maintaining numerical stability. As $`\mu `$ is dependent on the velocity divergence, any large velocity gradients are smoothed out by the increased viscosity. If a shock occurs it is not resolved, but smeared out by a local increase in viscosity.

The diffusive flux $`𝐟=(\mu /\mathrm{Pr})TS`$, where the horizontal mean of the entropy gradient $`S0`$ i.e. the convection zone, and $`𝐟=(\mu c_\mathrm{p}/\mathrm{Pr})T`$ where the horizontal mean of $`S0`$ i.e. the radiation zone (the Prandtl number Pr is defined below). In the convection layer the SGS diffusive flux tends to smooth out entropy fluctuations and make the layer close to adiabatic (in analogy with turbulent mixing). The change in the form of the diffusive flux above the convection boundary is necessary because the SGS’s should continue to transport heat radially outwards though the top of the box. Away from the horizontal boundaries, $`𝐟`$ is close to zero. At the base, $`𝐟`$ is equal to $`\sigma T_{\mathrm{eff}}^4`$, where $`\sigma `$ is the Stefan-Boltzmann constant. The Prandtl number Pr $`=\nu /\kappa `$, where $`\nu `$ is the kinematic viscosity and $`\kappa `$ is the thermal diffusivity. In the simulations Pr = 1/3. Due the inclusion of radiative energy transport the effective Pr is actually much smaller and not constant.

## 3 Numerical integration: obtaining accurate statistics

To simulate the highly stratified SAL of the Sun we need to relax the initial layer and then compute accurate statistics. The former requires a long computation, while the latter a small timestep. In compressible hydrodynamics however, with an explicit numerical method, the timestep must be less than the time for a sound wave to traverse two adjacent grid points. This is known as the Courant-Friedrichs-Levy (CFL) stability criterion. Because of these considerations the simulations were done in two stages.

Firstly, using an implicit code in which the time step is restricted by the flow speed rather than the sound speed, the initial hydrostatic layer was allowed to adjust its thermodynamical structure until it was close to hydrodynamic equilibrium. This thermal adjustment phase took at least 5 hours of solar surface convection time.

Secondly, using a second order accurate explicit code, quantities were averaged over a time that was long enough for the averages to be independent of the integration time. As the explicit timestep is about five times smaller than the implicit timestep, prior to statistical averaging, the code was run for a few thousand timesteps. This allowed the simulation to adjust to the new timestep. The statistical convergence took at least an hour of solar surface convection time.

### 3.1 Thermal adjustment

The implicit code was the Alternating Direction Implicit Method on a Staggered grid (ADISM) developed by Chan & Wolff (1982). This code was used to relax the fluid to a self consistent thermal equilibrium.

The entire layer was assumed to be relaxed when

* the energy flux leaving the top of the box was within 5 $`\%`$ of the input flux at the base;
* the horizontally averaged vertical mass flux was less than $`10^4`$ g/$`\mathrm{cm}^2/\mathrm{s}`$ at every vertical level;
* the overall thermal structure did not change much over time;
* and the maximum velocity in the box was roughly constant.

These criteria must be satisfied, before any useful statistical data can be gathered.

### 3.2 Statistical convergence

A second order explicit method (Adams-Bashforth time integration) gathered the statistics of the time averaged state. This code is much more accurate than ADISM, for example the mean energy flux leaving the top of the box was within 1 % of the input flux at the base. On a 667 MHz Alpha processor, each integration step on a $`80\times 80\times 80`$ grid, required about 5 seconds of CPU time.

The time required for statistical convergence depends on the particular quantity being averaged. Conserved quantities converge very fast. For example, the horizontally averaged vertical mass flux was less than $`10^5`$ g/$`\mathrm{cm}^2/\mathrm{s}`$ after 5 minutes of solar time integration. On the other hand, rms velocities converged in a few eddy turnover times (at least 30 minutes of solar time). While second order turbulent quantities, such as the horizontal Reynolds stress, took even longer to converge (at least 80 minutes of solar time).

#### 3.2.1 Some statistical definitions

In a turbulent fluid a quantity $`q`$ can be split into a mean and a fluctuating part,

$$q=\overline{q}(z)+q^{}(x,y,z,t).$$
(8)
The overbar represents a combined horizontal and temporal average, i.e.

$$\overline{q}(z)=\frac{1}{t_2t_1}\underset{t_1}{\overset{t_2}{}}\left(\frac{1}{(L_xL_y)}q𝑑x𝑑y\right)𝑑t.$$
(9)
$`t_1`$, is a time after the system has reached a self-consistent thermal equilibrium (the thermal adjustment time). $`L_x`$ and $`L_y`$ are the horizontal widths of the box in the x and y direction respectively. The time required for statistical convergence is $`t_2t_1`$.

The rms value of a quantity $`q`$ is defined as

$$q^{\prime \prime }=\overline{q^2}\overline{q}^2,$$
(10)
while the correlation coefficient of two quantities $`q_1`$ and $`q_2`$, is defined as

$$C[q_1^{}q_2^{}]=\frac{\overline{q_1q_2}\overline{q_1}\overline{q_2}}{q_1^{\prime \prime }q_2^{\prime \prime }}x.$$
(11)
As the simulations have periodic side boundaries, symmetry requires that

* $`C[v_x^{}v_y^{}]=0`$
* $`v_x^{\prime \prime }=v_y^{\prime \prime }`$

The run of $`v_x^{\prime \prime }`$ and $`v_y^{\prime \prime }`$ after 80 minutes of time integration is shown in Fig. 1. The closeness of the horizontal velocities confirms that the simulation is close to statistical convergence.

By examining many different simulations we found that the run of $`C[v_x^{}v_y^{}]`$ is generally a much stricter test of convergence. Fig. 2 shows $`C[v_x^{}v_y^{}]`$ measured for 4 different integration times. Even after 80 minutes, $`C[v_x^{}v_y^{}]`$ is still about 0.1 near the bottom. Convergence is faster in the upper layers because when the convecting fluid elements move up through the stratification, their rapid expansion smoothes out small scale fluctuations.

## 4 Main results

### 4.1 Effect of domain size

After substantial numerical testing which is described in detail in the appendix, we found that :

1. An undesirable effect of the impenetrable horizontal boundary at the bottom was too speed up the overall flow. This produced an artificially high convective flux. Just below the surface the radiative flux is a significant fraction of the total flux. As the total energy flux is fixed, to accommodate the increased convective flux, the radiative flux had to reduce. This was achieved by an (unphysical) drop in the temperature gradient $``$ in the SAL region.

To avoid this, the lower boundary has to be positioned far enough away from the surface, so that the velocity at the base is small and uncorrelated from the motions near the surface.
2. If the width of the box was too small then the turbulent kinetic energy of the granules was artificially small. This is because the movement of the larger granules was restricted by the walls of the box.

To avoid this, the aspect ratio was doubled until the turbulent kinetic energy was unaffected by a further increase in aspect ratio. For the Sun this required a width of 2.7 Mm.

### 4.2 Comparison with previous 3D numerical simulations

#### 4.2.1 Vertical velocity

Using the SN code, Asplund et al. (2000) computed a series of simulations of solar surface convection in a domain with a depth of about 4 Mm and width of 6 Mm. The best resolution in Asplund et al. was $`200^2\times 82`$. The rms vertical velocity and the mean velocity are shown in Fig.s 3 and 4 in their paper. Corresponding velocity plots from our best model (model C in the appendix) are shown in Fig. 3. Note that we have radially upwards as the positive direction, the opposite to the Asplund paper. Despite all the differences between the SN and KC approaches, away from the boundaries, the rms and mean vertical velocities in our best model are very similar to those in the best Asplund et al. simulation. The most noticeable differences occur at the top and bottom because in our simulations the vertical velocity is forced dropped to zero. In Asplund et al. the transmitting boundaries allowed the velocity to decrease more gently.

#### 4.2.2 Superadiabaticity

The superadiabaticty $`_{\mathrm{ad}}`$ for the MLT, and models KC2 and C, are plotted in Fig. 4. The abscissa is the radius of the simulation divided by the radius of the Sun. For model C, the superadiabaticity has a maximum of about 0.6, which is close to the value given by the SN code (see Fig. 3 in Rosenthal et al.). The position of the top of the convection layer (as determined by the Schwartzchild criterion) was pushed out further in KC2 than in C. This is because of the higher turbulent pressure.

When the SAL was moved outwards the convective efficiency was reduced and radiation was forced to carry more of the total flux. This resulted in an increase in the height of the SAL in C (triple dot dashed line), compared to the MLT (crosses). However, in KC2 (dashed line) the convection was speeded up by the lower boundary and thus is (incorrectly) more efficient than the MLT. This resulted in a drop in the height of the peak of the SAL compared to the MLT.

### 4.3 Comparison with observational results: p-mode oscillation frequencies

#### 4.3.1 Implementation of 3D turbulence into 1D stellar models

By analogy with the work of Lydon & Sofia (1995) on magnetic effects, the 3D turbulence is parameterised in terms of two quantities (Li et al. 2002), the turbulent kinetic energy per unit mass,

$$\chi =\frac{1}{2}v_{}^{\prime \prime }{}_{}{}^{2}$$
(12)
where $`v_{}^{\prime \prime }{}_{}{}^{2}=v_{x}^{\prime \prime }{}_{}{}^{2}+v_{y}^{\prime \prime }{}_{}{}^{2}+v_{z}^{\prime \prime }{}_{}{}^{2}`$ is a dimensional quantity, and an anisotropy parameter,

$$\gamma =1+2(v_z^{\prime \prime }/v^{\prime \prime })^2.$$
(13)
The $`z`$ direction is parallel to the radial direction.

The turbulence is included by incorporating 2 new variables $`\gamma `$ and $`\chi `$ into the 1D stellar model. To understand how this is works, consider a perfect gas in which the ratio of the specific heats $`\gamma =c_\mathrm{p}/c_\mathrm{v}`$, the internal energy $`e=c_\mathrm{v}T`$ and the gas pressure $`P=\rho RT`$. The quantities $`c_\mathrm{v}`$ and $`R`$ are for unit mass of a gas and $`T`$ is the temperature. The previous three equations can be expressed as,

$$\gamma =1+P/(\rho e)$$
(14)
If we replaced gas quantities by turbulent quantities, i.e. $`P`$ by $`P_{\mathrm{turb}}(=\rho v_{z}^{}{}_{}{}^{\prime \prime }{}_{}{}^{2})`$ and $`e`$ by $`\chi `$ and rearrange, then we would get:

$$P_{\mathrm{turb}}=(\gamma 1)\rho \chi .$$
(15)
where $`\gamma `$ is defined in terms of turbulent quantities (i.e. rms velocities). Including $`\chi `$ and $`\gamma `$ is equivalent to including $`\chi `$ and $`P_{\mathrm{turb}}`$ in the 1D stellar model. The two variables $`\gamma `$ and $`\chi `$ are included in the mathematical system in a self consistent manner (see Li et al. for the full mathematical treatment)

For example, the equation of state for the 1D model becomes,

$`\rho =\rho (P_T,T,\chi ,\gamma )`$
where $`P_T=P_{\mathrm{gas}}+P_{\mathrm{rad}}+P_{\mathrm{turb}}`$.
The continuity equation and the equation of transport of energy by radiation remain the same regardless of turbulence. In terms of $`\gamma `$ and $`\chi `$, the hydrostatic equilibrium is,

$$\frac{P_T}{M_r}=\frac{GM_r}{4\pi r^4}\frac{2(\gamma 1)\chi }{4\pi r^3}.$$
(16)
where $`M_r`$, $`G`$ and $`r`$ have their usual standard stellar model (SSM) definitions. The energy conservation equation,

$$\frac{L_r}{M_r}=ϵT\frac{dS}{dt},$$
(17)
is also affected by the inclusion of $`\chi `$, as

$$TdS=dU+(P_TP_{\mathrm{turb}})d(1/\rho )+d\chi ,$$
(18)
where $`d\chi `$ represents the work done by turbulence. The quantities $`L_r`$ and $`ϵ`$ have their usual SSM definitions.

The convective flux in the 1D model $`F_{\mathrm{conv}}=\rho TDSV_{\mathrm{conv}}`$, now includes the turbulent kinetic energy flux $`\chi V_{\mathrm{conv}}`$ (without turbulence $`\chi =0`$, so $`d\chi =\chi `$). $`V_{\mathrm{conv}}`$ is the MLT convective velocity and DS is the entropy excess. Equation 18 is the ‘non-standard form’ of the 2nd law as described in Lydon & Sofia. As turbulence is in general anisotropic it is wrong to treat the work done by turbulence as $`P_{\mathrm{turb}}dV`$. This is essentially because $`P_{\mathrm{turb}}`$ is a tensor and $`P_{\mathrm{gas}}`$ is a scalar.

The equation of energy transport by convection, does not change in form, but $``$ is different from that without turbulence. The equations that govern envelope integrations are changed accordingly (see Li et al.).

#### 4.3.2 Solar p-mode oscillation frequencies

To investigate the effect of turbulence on the solar p-mode oscillation frequencies, the runs of $`\chi `$ and $`\gamma `$ for models KC2 or C were incorporated into the 1D stellar model using the Li et al method. For each run, the p-mode frequencies for $`l=0,1,2,3,4,10\mathrm{}100`$, were computed using Guenther’s pulsation code, (1994), under the adiabatic approximation. When comparing the computed and observed eigenfrequencies, one should really include the proper modeling of radiative gains and losses by computing the non-adiabatic frequencies (Guenther 1994). However, to help isolate the effects of turbulence on the p-mode frequencies, we computed the adiabatic frequencies. The difference between observed and computed adiabatic p-mode frequencies for the standard solar model is shown in Fig. 7. The frequency difference is scaled by the mode mass $`Q_{nl}`$ (e.g., Christensen-Dalsgaard & Berthomieu 1991).

The $`Q_{nl}`$ weighting attempts to correct the p-mode frequency differences, by removing the dependency on mode inertia. The mode differences with high mode inertia get more weight than those with lower inertia. This tends to give greater weight to the low $`l`$, deep penetrating modes and less significance to the very high $`l`$ shallow modes. The result is the removal of the $`l`$ dependence in the frequency differences due to perturbations in the near surface regions. We are left with only $`n`$ dependence, which is equivalent to frequency dependence. The $`Q_{nl}`$ weighting enables us to see that the discrepancy between the observed and the computed adiabatic frequencies is worse at the high frequencies, thus pointing to a problem in the surface layers. Without the $`Q_{nl}`$ weighting, the mismatch would appear to varying degree, at all frequencies and therefore it would not be as easy to claim that the outer layers are responsible for the mismatch.

In the statistically steady state (after the flow is thermally relaxed), the turbulent kinetic energy flux should be proportional to the energy input rate at the base. The turbulent kinetic energy $`\chi `$ should be independent of the geometry of the computational domain. If the box is too shallow or the width too narrow, then $`\chi `$ can be wrong. To illustrate the effect of the geometry on the frequencies, $`\gamma `$ and $`\chi `$ from the KC2 model was used to compute a stellar model. The computed adiabatic frequencies with turbulence derived from KC2 are shown in Fig. 7. The resulting frequencies are much worse than the standard solar model. However, if we use model C instead of KC2, then the derived frequencies were improved considerably (Fig. 7). The difference between the computed and the observed frequencies is less than 5 $`\mu `$ Hz. The inclusion of turbulent kinetic energy in 1D models in addition to turbulent pressure, has a similar effect on the adiabatic frequencies to the inclusion of radiative losses/gains (computation of non-adiabatic frequencies). Both reduce the difference between the adiabatic frequencies and the observed frequencies by an order of magnitude. This mean that the turbulent correction to the p-modes is as important as radiative losses/gains.

### 4.4 Characteristics of the solar granules

#### 4.4.1 Size

Fig. 8 shows a set of contours of the instantaneous vertical velocity. From top to bottom the frames depict depths of -0.14 Mm, 0.03 Mm, 0.2 Mm and 1.0 Mm. The depth was measured positively inward from the visible solar surface. The contours themselves have been derived from a simulation of $`80^3`$ grid points, in a domain with a horizontal area of 3.75 Mm $`\times `$ 3.75 Mm, and a depth of 2.5 Mm. At this instant in time there appear to be between 3 and 4 granules near the surface. The thick black lines represent the strong downflows which occur at the sides of the granules. The lighter regions denote upflowing fluid or weak downflows. The granular pattern seems to persist over the first three frames, suggesting that the granule structure is correlated over most of the SAL (between 0 Mm and 0.25 Mm). By the fourth frame there is not much sign of granulation. As the last contour is still more than 2 pressure scale heights above the base, the break up of the granules is probably not due to the impenetrable bottom boundary. In other words the bottom boundary is far enough away from the surface. Furthermore, as the shape of the granules does not appear to be influenced by the side walls, the width of the box is also big enough.

The thin closely spaced dark vertical/slanted parallel lines that appear in the first contour, are signs of two grid waves. This is an indication of insufficient damping and is only seen at the very top of the box. It is due to the low SGS viscosity, which is proportional to density ((equation 7). This is is an undesirable numerical effect. But as the oscillations occur only at the very top of the box, where the density is very small, their effect on the convection below is minor. Before starting the time integration, we were able to damp out such oscillations, by slightly increasing the viscosity near the top. As these contours only represent one instant in time, such a picture can only provide a limited idea of the nature of the solar granulation. If the same contours are computed at a later time, a completely different picture would be seen. For example the third frame of figure 22 is the same simulation, but the contour was computed 2 solar minutes later.

Clearly any useful data can only be derived from long time averaged statistics. One characteristic vertical length scale associated with convective turbulence is the half-width of the 2-point vertical velocity correlation C\[$`v_z^{}v_z^{}]`$. If we assume convective eddies to have an aspect ratio of unity, then this length scale gives us an idea of the size of the eddies in a turbulent fluid. Using simulation D (described in the appendix), we computed C\[$`v_z^{}v_z^{}]`$ at a series of depths from the top of the box inwards. The results are shown in Fig. 9. From about 0 Mm to 1 Mm the eddy size remains close to 600 km. From the begining to the end of the plateau, both density and pressure have increased by an order of magnitude. Over this range the eddy does not seem to be affected by stratification. However, over the next Mm the eddy size increases to about 1100 km. After that the lower boundary starts to influence the eddies. The bends at the far left and far right of the plot signal the approach of the upper and lower boundaries.

The velocity contours and the velocity correlation length both suggest that the granules have some coherent vertical structure. The contours show intergranule lanes that do not shift much between 0 and 0.2 Mm, i.e. the granules are like cylinders. The half-width suggests a vertical size of about 600 km, between 1 Mm and the surface. This length scale does not seem to be affected by the stratification. The strong vertical coherence implies that the convection is dominated by strong downflows that originate close to the surface. SN showed that solar granulation is primarily driven by radiative cooling at the surface. When ascending fluid approaches the surface, it looses heat so rapidly that very strong downflows are created. The downflows are not deflected by the surrounding fluid until they are about 1 Mm below the surface, by which point they have weakened enough to be deflected/broken up by the surrounding fluid motion.

#### 4.4.2 Heat transport

In Chan & Sofia (1987,1989) the 2-point vertical velocity and temperature correlations were found to scale with pressure scale height rather than density scale height. This study was for deep efficient convection. In shallower layers Kim et al. (1995) found that while the vertical velocity correlation scaled with both density and pressure scale height, the temperature could not be scaled with either. The difference from Chan & Sofia’s results, was caused by the inclusion of the coupling of the partial ionisation with convection (they were treated separately by Chan & Sofia). As a fluid parcel moves upwards through regions of decreasing ionisation, it liberates ionisation energy which increases the buoyancy of the fluid parcel. However, because radiation was modelled by the diffusion approximation (which is known to break down at some point in the SAL), the Kim et al. simulations could only include the lower half of the SAL.

The present models include all of the SAL. Figs. 12 and 12 show $`C[v_z^{}v_z^{}]`$ and $`C[S^{}S^{}]`$ at depths of 0.06 Mm, 0.2 Mm, 0.3 Mm, 0.47 Mm and 1 Mm, denoted by solid, small-dashed, dot-dashed, triple-dot dashed and long-dashed lines, respectively. The curves have been centered about their respective maxima. As the plots do not coincide, neither quantity scales with the pressure scale height. Furthermore, if density is used rather than pressure, there was no noticeable improvement (not shown).

The width of the $`C[v_z^{}v_z^{}]`$ distribution decreases between 0.06 Mm and 1 Mm. This is because the width of the geometric distribution is roughly constant (i.e. the half-width is always about 600 km, as shown in the previous figure), but the scale height increases inwards. The entropy correlation looks very different. For $`\mathrm{\Delta }\mathrm{ln}P<0`$ the entropy correlation for the three most shallow depths drops rapidly, while for $`\mathrm{\Delta }\mathrm{ln}P>0`$ it shows a smoother decay. This implies that ascending parcels near the surface, lose heat very quickly, while descending parcels maintain their heat content appreciably longer. Over the last Mm before the surface, the eddies may have about the same diameter, but they transport less and less heat as the surface is approached.

This behaviour can be understood by examining the correlation $`C[v_z^{}S^{}]`$ for upflows and downflows. Fig. 12 shows that the granules are more efficient at transporting entropy downwards than upwards. At the peak of the SAL, $`C[v_z^{}S^{}]`$ is 0.9 for downflows and 0.6 for upflows. This explains the one-sidedness of $`C[S^{}S^{}]`$. Consider a fluid parcel that is initially just below the surface (say at a depth of 0.1 Mm). If that parcel moves downwards it will maintain its original entropy for about 1 PSH, whereas if the parcel moves upwards it it will loose its original entropy in 0.5 PSH. At greater depths $`C[v_z^{}S^{}]`$ is similar in the upflows and downflows, hence $`C[S^{}S^{}]`$ is more symmetric deeper down.

#### 4.4.3 A variable mixing length ?

Over the last 1000 km before the surface, the eddies seem to have a nearly constant diameter of about 600 km. This would seem to imply a mixing length ratio (the correlation length divided by the pressure scale height) which increases towards the surface. However, as the surface is approached, radiation plays a greater and greater role in heat transport, so the mixing length ratio should approach zero at the convection surface. This paradox can be partly resolved by noticing that as the surface is approached, the eddies become less and less efficient at transporting heat upwards. This is shown in the plot of $`C[v_z^{}S^{}]`$ for the upflows. A more reasonable candidate for the mixing length might be the product of $`C[v_z^{}S^{}]`$ and the velocity correlation half width.

## 5 Summary and implications for solar models

As we intend to use our code to model the SAL in other stars, this paper is an important benchmark for future studies. From these simulations of the Sun, we found that:

1. The impenetrable lower boundary needs to be far enough away from the SAL so that by the time the fluid reaches the boundary, the velocities have become both weak and uncorrelated from velocities in the SAL. If the lower boundary is too close to the SAL, the kinetic energy will be overestimated.
2. The horizontal cross section needs to be large enough so that the side-walls do not restrict the movement of the granules. If the box width is too small then the kinetic energy will be underestimated.
3. There is a region close to the surface in which the vertical correlation length remains constant, even though the density and pressure vary by an order of magnitude.
4. In that region the mixing length theory, which assumes the correlation length to be a constant multiple of the pressure scale height, will not work.
5. The final equilibrium state is not strictly dependent on the initial model atmosphere (see appendix, section 6.1).

While these results have been found in a simulation of the Sun, it is reasonable to assume that similar criteria would apply to the convection-radiation transition layers in other stars. The effect of the boundaries should certainly be considered when simulating other convection-radiation layers. The correlation length (half-width) seems to be a very robust feature of the solar granules and could easily be measured in other computations. For example, in a recent proceedings (Robinson et al 2002) we describe the application of this model to the Sun at the sub-giant and the start of the red giant branch. Preliminary results indicate that, near the surface of each convection zone, the ratio of the half-width to the stellar radius depends directly on the surface gravity.

This research is supported in part by NASA grant NAG5-8406 (FJR, PD, SS). YCK is supported by a Department of Astronomy, Yonsei University grant (2001-1-0134). DBG’s research is supported in part by a grant from NSERC of Canada. We also would like to acknowledge computer support from the Centre for High Performance Computing at Saint Mary’s University in Canada. Finally, we thank H.-G. Ludwig for helpful comments.

## 6 Appendix : Numerical tests

Numerical simulations of solar granulation need to be robust. Part of this means that the walls of the computational box should not control the radiative hydrodynamics of the interior (e.g. the size of the granules, the structure of the SAL, the turbulent pressure, etc.). Nor should any small changes in the initial conditions affect the final equilibrium. The following simulations were designed to minimise these uncertainties.

The features of the individual simulations are summarised in Table 1. In this table the geometric width in both the x and y directions is given in Mm, while the depth is given in Mm and as the number of pressure scale heights (PSH). Models A and B were created by vertically extending model KC2, while C, D and E were made by periodically extending model B. After periodic extension, the horizontal resolution in models D and E was halved. Each time a new model was created, it was allowed to return to thermal equilibrium before the statistical integration was begun.

For each model we computed three non-dimensional statistical quantities:

1. The turbulent pressure divided by the mean gas pressure,

$$P_{turb}^{}{}_{}{}^{}=\overline{\rho }v_{z}^{\prime \prime }{}_{}{}^{2}/\overline{P},$$
(19)
where the ‘\*’ denotes a non-dimensional quantity.
2. The turbulent kinetic energy per unit mass divided by the isothermal sound speed squared,

$$\chi ^{}=\frac{1}{2}\frac{v_{}^{\prime \prime }{}_{}{}^{2}}{c_{\mathrm{s}}^{}{}_{}{}^{2}}$$
(20)
where $`v_{}^{\prime \prime }{}_{}{}^{2}=v_{x}^{\prime \prime }{}_{}{}^{2}+v_{y}^{\prime \prime }{}_{}{}^{2}+v_{z}^{\prime \prime }{}_{}{}^{2}`$ and $`c_\mathrm{s}=\sqrt{\overline{P}/\overline{\rho }}`$ is the isothermal sound speed.
3. The superadiabaticity

$$_{\mathrm{ad}}=\overline{\frac{\mathrm{ln}T}{\mathrm{ln}P}}\overline{_{\mathrm{ad}}}$$
(21)
where $`_{\mathrm{ad}}`$ is computed using the OPAL equation of state

From now on we will drop the ’\*’ on $`P_{\mathrm{turb}}`$ and $`\chi `$, for convenience. All quantities in this section are non-dimensional.

### 6.1 The influence of varying the initial conditions

Here we compare the original KC simulation with simulation KC2 listed in Table 1. Both simulations have the same input physics (opacity tables, equation of state). As noted in section 2.1, the only difference between KC2 and KC, is that the 1D initial models were based on slightly different model stellar atmospheres \[i.e. the Eddington T($`\tau `$) relation in KC, vs. the empirical Krishna-Swamy $`T(\tau )`$ relation in KC2\]. As the Krishna-Swamy relation is from observations of the Sun, it was closer to the final state of the solar simulation.

Fig. 13 shows the initial hydrostatic (based on the MLT) and relaxed hydrodynamic runs of Log P vs Log T for the two simulations. The dotted and solid lines are the initial plots for KC and KC2, respectively. The horizontal and temporal averages of the relaxed states, are denoted by diamonds and crosses for KC and KC2, respectively. To demonstrate statistical convergence, the data of model KC2 was initially averaged over about 600 seconds of solar surface convection, whereas in KC, the averaging time was 2240 seconds. After a sufficient amount of time, KC2 (crosses) converged to KC (diamonds). The MLT provides the initial thermal structure, but the hydrodynamical turbulence shapes the final structure. Provided the initial conditions are not too different, the hydrodynamics converges to the same equilibrium state.

This is an important preliminary result. Even if the initial atmospheres are slightly different, the hydrodynamics dictates the final equilibrium. The thermodynamics of the resultant 3D simulation does not depend on the precise details of the MLT. The final equilibrium is not determined by the exact value of the mixing length ratio in the stellar model, on which the initial state of the simulation is based. The choice between the two initial conditions in the atmosphere, does not affect the turbulent pressure, turbulent kinetic energy, or the superadiabatic temperature gradient. These quantities will depend on the geometry and grid resolution in the box.

### 6.2 The influence of the lower impenetrable boundary

The structure of the SAL in KC2 (or KC which is identical) is at odds with the SAL constructed with the Canuto-Mazzitelli (1991) approach (see Fig 1. in Demarque, Guenther & Kim (1999)) or presented in the numerical simulation by Rosenthal et al. The most obvious difference is that KC2 predicts a broader SAL with a maximum of about 0.4, while Canuto-Mazzitelli predicts a peak of about unity and Rosenthal et al a peak of 0.6. The height of the peak is related to the convective efficiency with respect to radiative heat transport.

Models KC2, A and B differ only in their domain depths. Simulation A was constructed by adding 3 extra pressure scale heights to the base of model KC2 (see Table 1). We first thermally relaxed the lower layer, because with only hydrostatic support (no turbulent pressure), the overlying layer collapsed, and the simulation crashed. Also, when joining the two layers, the turbulent viscosity was temporarily increased. This smoothed out the fluctuations that were produced by suddenly removing the lower boundary. Model B was made by adding a hydrostatic layer computed using the MLT to model A. This increased the depth from 2.1 Mm to 2.8 Mm. As the turbulent pressure near the bottom of model A was small (about 1 % of the gas pressure), in this case we did not need to relax the hydrostatic lower layer separately.

#### 6.2.1 Effect on kinetic energy and turbulent pressure

To demonstrate the effect of the lower boundary on the turbulent flow, we computed the ratio of the horizontal kinetic energy to the gas pressure. This is equivalent to the horizontal turbulent Mach number squared, i.e.

$$\chi _{\mathrm{horiz}}=\frac{1}{2}\overline{\rho }(v_{x}^{\prime \prime }{}_{}{}^{2}+v_{y}^{\prime \prime }{}_{}{}^{2})/\overline{P}=\frac{1}{2}(v_{x}^{\prime \prime }{}_{}{}^{2}+v_{y}^{\prime \prime }{}_{}{}^{2})/c_{s}^{}{}_{}{}^{2}$$
(22)
Fig. 14 shows $`\chi _{\mathrm{horiz}}`$ for models KC2, A and B. In KC2 the magnitude of $`\chi _{\mathrm{horiz}}`$ shoots up at the base. The net effect of the fast downflows striking the lower boundary was to speed up the overall flow. This increased the total turbulent kinetic energy $`\chi `$ throughout the box (Fig. 15).

The upturn near the bottom of each layer, is clearly reduced as the boundary is moved deeper, and is almost eliminated when the depth is 2.8 Mm. This problem is much less severe at the top boundary because the underlying region is sub-adiabatic (the radiative zone). The vertical velocity is already small when the flow hits the top boundary. The effect on $`P_{\mathrm{turb}}`$ of moving down the lower boundary is shown in Fig. 16. In model KC2, $`P_{\mathrm{turb}}`$ drops sharply as the impenetrable bottom is approached. While for A and B it has a much smoother decay.

#### 6.2.2 Effect on superadiabaticity

The run of $`_{\mathrm{ad}}`$ for the MLT, and models KC2, and C was shown previously in Fig. 4. When the SAL is moved outwards by turbulent pressure, the convective efficiency is reduced and radiation is forced to carry more of the total flux. This should result in an increased height of the SAL compared to the MLT. However, in KC2 the convection is overly turbulent and is thus (incorrectly) more efficient than convection computed via the MLT. This results in a drop in the height of the SAL compared to the MLT.

### 6.3 The influence of the upper impenetrable boundary

As the top boundary is also impenetrable we need to ensure that its position does not effect the convection either. We need to prove that the decrease of vertical velocity between the convection surface (where the depth = 0 Mm) and the top of the box, is not due to the top boundary. Rather it should be due to the stable layer at the top, i.e. convection-radiation losses.

To address this issue we damped the horizontal velocity at the top by replacing the stress free boundary with a no-slip top. The flow is then relaxed and statistics are gathered as usual. Fig. 17 shows $`_{ad},v_x^{\prime \prime }`$ and $`v_z^{\prime \prime }`$ for simulations with slip and no-slip top boundaries. The velocities are non-dimensionalised by the local isothermal sound speed. This conveniently enables them to be on the same axis as $`_{ad}`$. Apart from the top 100 km, where the two $`v_x^{\prime \prime }`$’s diverge because of the different top boundary conditions, the rms vertical and horizontal velocities are nearly the same for the no-slip and stress free top boundaries. This implies the drop in vertical velocity near the top is primarily due to to convection to radiation losses, rather than the top boundary. Furthermore, the horizontal velocity in the upper atmosphere does not affect the SAL structure (i.e. $`_{ad}`$) much.

### 6.4 The influence of the domain width

Models B, C, D and E have widths of 1.35 Mm, 2.7 Mm, 2.7 Mm and 5.4 Mm. To judge the effect of width on simulations without changing the grid spacing, we should compare B to C and D to E. All four models have similar $`P_{\mathrm{turb}}`$ and $`_{\mathrm{ad}}`$ (Figs. 19 and 19). Increasing the width seems to have a minor effect on either $`P_{\mathrm{turb}}`$ or $`_{\mathrm{ad}}`$. In general it is essential to resolve turbulent motions inside the SAL region. As this is only about 250 km thick, the grid spacing in the SAL need to be very small. To show that all our simulations have resolved the SAL, we have plotted the individual vertical grid points on the figure.

The variation of $`\chi `$ with domain width is more interesting. This is shown in Fig. 20. When the width is increased from 1.35 Mm to 2.7 Mm, $`\chi `$ increases (especially near the top). However, when the width is increased from 2.7 Mm to 5.4 Mm there is only a very small change in $`\chi `$. A box width of 2.7 Mm seems to be a sufficient to resolve $`\chi `$ for solar granules.

To provide further evidence that 2.7 Mm is a large enough box width, we computed $`v_x^{\prime \prime },v_y^{\prime \prime }`$ and $`v_z^{\prime \prime }`$ for models D and E. Fig. 21 shows all three velocity components for both D and E. Doubling the width has only a minor effect on any particular velocity component. The small differences in the deeper part, are due to insufficient convergence. As the domain has periodic lateral boundaries, eventually the horizontal velocities should all be the same.

Fig. 22 shows contour plots of the instantaneous vertical velocity, for a horizontal cross-section, for models of different widths. The uppermost frame shows that not even a single granule can fit in the box when the width is 1.35 Mm, while about 2 fit easily in to the box in the second frame. In the final frame about 9 or 10 larger cells fit in to the box. Due to computational restrictions the last figure is computed on a coarser mesh than the previous three frames so the granules are not clearly depicted.

### 6.5 The influence of numerical resolution

While the main point of these tests is to show that the walls of the computational box can have a significant effect on granular convection, we can also partially address the effect of changing the horizontal grid spacing. Models C and D produced very similar $`P_{\mathrm{turb}}`$ and $`_{\mathrm{ad}}`$, while $`\chi `$ differed only slightly. As $`\chi `$ depends on the horizontal as well as the vertical velocity, this probably reflects the sensitivity of the horizontal component of kinetic energy to horizontal grid resolution. This is particularly noticeable near the top where the flow is mostly horizontal.


text:
# Neutrino mean free path in neutron stars

## 1 INTRODUCTION

The interaction of neutrinos with baryons has been mostly studied in connection with the stability of nuclei. But it also plays a crucial role in the thermal evolution of supernovae and neutron stars, where the nuclear medium is far from stability. A large effort has been devoted to study the production of neutrinos via direct or modified URCA processes  and, more recently, via the bremsstrahlung from nucleons in the strong magnetic field of neutron stars. The propagation of neutrinos in neutron matter has received less attention. Being the experimental information on this topic still insufficient, the theoretical predictions have to be developed further, extending the microscopic description of many-body systems to a wide range of density, temperature and isospin asymmetry. This means extending the study of the nuclear equation of state far away from the saturation point.

In this note we report on the calculation of the response function to weak interaction coupling in neutron matter. It is based on the equation of state (EOS) predicted by the Brueckner theory, including relativistic and finite-size effects (nucleon resonances). Recently it has made a step forward as to reproducing the empirical saturation properties of nuclear matter. The effective NN interaction is extracted from the Brueckner theory and cast in terms of the Landau parameters. Then, the response function is calculated in RPA along with the neutrino mean free path in neutron matter for various conditions of neutron density and temperature.

## 2 NUCLEAR MATTER EOS FROM THE BRUECKNER THEORY

The EOS of nuclear matter has been investigated from several microscopic approaches. Their common feature is that of adopting as $`V_{NN}`$ interaction a realistic potential, i.e. fitted on the experimental phase shifts of the nucleon-nucleon scattering in the vacuum. Among them the Brueckner theory has proved to be very powerful for its capability of incorporating both relativistic and nucleon finite-size effects . This feature is illustrated in Fig. 1, where a calculation with the interaction Argonne $`V_{18}`$ is plotted. In the first order of the hole-line expansion (BHF approximation), the saturation density turns out to be too large, when the NN interaction is just a two-body force (2BF). Extending the BHF calculation to include a three-body force (3BF) effects one improves the agreement with the empirical saturation properties of nuclear matter. The origin of 3BF is twofold: the $`\overline{N}N`$ excitations missing from the non relativistic Brueckner theory, amount to a special 3BF as shown by G.E. Brown et al.; low-lying nucleonic excitations as intermediate states in the interaction between two nucleons are pure 3BF. As shown in Fig.1 both effects play a important role for the reproduction of the empirical saturation properties of nuclear matter.

The predictions of the EOS can be extended far beyond the saturation density of cold and symmetric nuclear matter. A variety of physical situations can be studied, including neutron-rich and spin-polarized matter at finite temperature (here we limit ourselves to zero temperature). For such cases there are less physical constraints: the isospin-symmetry energy at the saturation point is constrained by the empirical value of $`30MeV`$ extracted from the nuclear mass table: finally, the spin-symmetry energy can be probed by spin-flip excitations of nuclei.

A set of Brueckner calculations have been performed in different states of spin and isospin polarization of nuclear matter, defined by

$$\delta =\frac{A_{}A_{}}{A},\beta =\frac{NZ}{A},$$
(1)
respectively. Here $`A_{}`$($`A_{}`$) are the number of nucleons with spin polarized upward (downward). In both cases the Brueckner calculations predict a quadratic law so that the isospin and spin symmetry energies describe completely the properties of isospin and spin polarized nuclear matter, respectively. In Fig. 2 the isospin symmetry energy and in Fig. 3 the spin symmetry energy are plotted as a function of density.

## 3 LANDAU PARAMETERS FROM THE BRUECKNER THEORY

The Brueckner theory basically provides us with the $`G`$-matrix, which is nothing else than the in-medium NN interaction. The latter is not suitable for calculations, but it can be cast in the form of effective interaction within the quasiparticle description of Fermi systems . Alternatively, for the present purpose one first determines from the EOS of nuclear matter fundamental equilibrium properties such as the compression modulus, symmetry energy and spin susceptibility and then calculates the corresponding Landau parameters.

In the Landau theory of Fermi liquids one writes the particle-hole residual interaction as

$$N(0)𝒱_{NN}=F+F^{}\stackrel{}{\tau }\stackrel{}{\tau }^{}+G\stackrel{}{\sigma }\stackrel{}{\sigma }^{}+G^{}\stackrel{}{\tau }\stackrel{}{\tau }^{}\stackrel{}{\sigma }\stackrel{}{\sigma }^{},$$
(2)
where $`F,F^{},G,G^{}`$ are the Landau parameters. Each of them is expanded in a Legendre series, for instance $`F=F_lP_l(cos\theta )`$ where $`\theta `$ is the Landau angle. The physical quantities related to the static response functions can be expressed in a simple way in terms of the Landau parameters. The response to density fluctuation, i.e. the compression modulus is directly related to $`F_0`$; the response to charge variation, i.e. the symmetry energy, is related to $`F_0^{}`$; the response to spin variation, i.e. the magnetic susceptibility, is related to $`G_0`$; finally, the response to spin-isospin variation, which is the isospin susceptibility, is related to $`G_0^{}`$.

From the Brueckner calculations we found that the energy per nucleon of nuclear matter as a function of $`\delta _n`$ and $`\delta _p`$ (Eq. (1) for neutron and proton, separately) can be fairly approximated by the quadratic form

$$E_A(\delta _n,\delta _p)E_A(0,0)=\mathrm{\Delta }_{nn}\delta _n^2+\mathrm{\Delta }_{pp}\delta _p^2+2\mathrm{\Delta }_{np}\delta _n\delta _p.$$
(3)
This equation is the basis of our calculation of the zero-order Landau parameters, which are plotted in Fig. 4 and Fig. 5. Notice that $`F_0`$ becomes less than -1 (negative compression modulus) at low density which corresponds to the unphysical spinodal region of the EOS, signaled by the multifragmentation events of heavy ion collisions.

## 4 NEUTRINO MEAN FREE PATH

Let us consider the propagation of neutrinos in pure neutron matter. The neutrino mean free path $`\lambda `$ is derived from the transport equation. One obtains

$$\frac{c}{\lambda _\stackrel{}{k}}=W_{fi}(kk^{})(1n_\stackrel{}{k}^{})+W_{fi}(k^{}k)n_\stackrel{}{k}^{},$$
(4)
where $`n_\stackrel{}{k}`$ is the occupation number of neutrinos and $`W_{fi}`$ is the neutrino-neutron scattering probability via the weak neutral current. The medium effects are incorporated into $`W_{fi}`$ via the dynamical form factors,

which, in turn, are related to the imaginary part of the dynamical response function. In the low-frequency limit $`\omega /(kv_F)1`$ the latter is expressed in a simple way in terms of the Landau parameters, which is the approximation adopted in the present calculation. The Landau parameters were extracted from the Brueckner prediction of the EOS of neutron matter. They are listed in Table 1 in comparison with other predictions. In Fig. 6 (left panel) our prediction for $`\lambda _\nu `$ in interacting neutron matter is reported as a function of the density in comparison with that in a free neutron gas. One can see that the large medium effects are due to the large value taken by both Landau parameters $`F_0`$ and $`G_0`$ in the interacting system. For densities in the range of neutron-star density $`\lambda _\nu `$ is by one order of magnitude larger than in a free Fermi gas. In the right panel of Fig. 6 the Hartree-Fock predictions from two Skyrme forces and Gogny force taken from Ref. are reported for comparison. In our case $`\lambda _\nu `$ is increasing much faster with density due to the more repulsive character of our 3BF. Fig. 7 shows the dependence of $`\lambda _\nu `$ on the temperature T. The expected T dependence of the Landau parameters is neglected. The relative weight of each elementary process giving rise to the neutrino propagation is plotted in Fig. 8. The role played by the Landau zero mode appearing as a singularity of the response function, is negligible.

In summary, we have presented a prediction of mean free path of neutrinos in neutron matter. Since the strong nuclear correlations enhance the propagation of neutrino one may predict a faster cooling of a neutron stars via neutrino escape than in a degenerate neutron gas. The neutron star temperature plays also important role in the neutrino propagation. Due to the large value of $`T`$ involved in proton-neutron stars, the present predictions should be improved by introducing $`T`$ dependent Landau parameters.

This work has been supported in part by the Chinese Academy of Science, within the one Hundred Person Project, the Knowledge Innovation Project of CAS under No. KJCX2-SW-N02, and the Important Pre-research Project of the Ministry of Science and technology under No.2002CCB00200, China.