format: text

text:
# The Optical Spectrum of the SN 1006 Supernova Remnant Revisited

## 1. INTRODUCTION

Among the known galactic supernova remnants (SNRs), SN 1006 stands out as one of the few with low extinction ($`E(BV)=`$0.11, Schweizer & Middleditch 1980), a well determined age and a reasonably well known distance ($`D`$2 kpc, Kirshner, Winkler & Chevalier 1987 (hereafter KWC87), Long et al. 1988, Laming et al. 1996, Winkler & Long 1997). This SNR was probably produced by a Type Ia explosion (Schaefer 1996) and appears as a limb brightened shell 30 in diameter at X-ray (Willingale et al. 1996, Winkler & Long 1997) and radio (Reynolds & Gilmore 1986) wavelengths. Using the Hopkins Ultraviolet Telescope (HUT), Raymond et al. (1995) first detected ultraviolet emission lines of He II $`\lambda `$1640, C IV $`\lambda \lambda `$1548, 1551, N V $`\lambda \lambda `$1239, 1243 and O VI $`\lambda \lambda `$1032, 1038 in SN 1006. Raymond et al. (1995) observed equal velocity widths for the emission lines detected by HUT. Combining this information with the relative strengths of the He II, C IV, N V and O VI lines, Laming et al. (1996) concluded that the lines observed by HUT were due to collisional excitation in 2600 km$`\mathrm{s}^1`$ nonradiative shocks with little or no electron-ion temperature equilibration at the shock front. In this case the high shock speed and low density ($`n`$ 1 cm-3) result in negligible radiative cooling downstream (hence the term ‘nonradiative shock’). All of the line emission is produced in a narrow ionization zone close to the shock front.

In the optical, SN 1006 appears as a faint network of Balmer-dominated filaments (van den Bergh 1976, Long, Blair & van den Bergh 1988, Winkler & Long 1997). These filaments emit an optical spectrum dominated by the Balmer lines of H, produced in the same ionization zone as the ultraviolet lines. Balmer line profiles consist of two components: a narrow component produced when cold ambient neutrals are overrun by the shock and a broad component produced when charge exchange with protons generates a hot neutral population behind the shock (Chevalier & Raymond 1978, Chevalier, Kirshner & Raymond 1980). The broad component width is directly proportional to the postshock proton temperature, while the ratio of broad to narrow flux is sensitive to the electron-proton temperature equilibration and the preshock neutral fraction. These diagnostics are invaluable tools for measuring shock speeds in supernova remnants and studying the efficiency of collisionless heating processes in high Mach number, nonradiative shocks.

## 2. SPECTROSCOPIC OBSERVATIONS

We obtained a deep long-slit spectrum of the NW filament in SN 1006 from the 4-m Blanco telescope at CTIO on 1998 June 24 (UT). The RC Spectrograph was used with the 527 line mm-1 grating KPGL3, the Blue Air Schmidt camera, and the Loral 3K CCD. This combination covered the spectral range 3760$``$7480 Å at a dispersion of 1.22 Å pixel-1. For our slit width of 2$`\stackrel{}{\mathrm{.}}`$0, the instrumental resolution was approximately 4.5 Å (205 km$`\mathrm{s}^1`$ at H$`\alpha `$) along the entire spectrum. The plate scale for this configuration was 0$`\stackrel{}{\mathrm{.}}`$51 pixel-1 and the slit length was 5. We rotated the slit to PA = 63.9 and placed its center at coordinates 15h02m13$`\stackrel{\mathrm{s}}{\mathrm{.}}`$5, $``$414522 (2000), as shown in Fig. 1. We chose this position to cut obliquely across the NW filament in order to capture as much of the bright region as possible while giving enough background sky, away from the filament, for a good sky subtraction. It is virtually identical to the position used by KWC87. We obtained 4 exposures at this location, with a total exposure time of 8400 s, which we combined into a single two-dimensional spectrum following initial reduction.

The data were reduced using standard IRAF111IRAF is distributed by the National Optical Astronomy Observatories, which is operated by the AURA, Inc. under cooperative agreement with the National Science Foundation procedures of bias subtraction, flat fielding and illumination correction. Spectra of a HeNeAr lamp taken at the beginning and end of the series of object frames were used for wavelength calibration. We subtracted the night sky contribution from the two-dimensional spectrum using emission from the two ends of the slit. The H$`\alpha `$ line profiles (both broad and narrow components) were nearly constant along the observed filament; therefore, we integrated the emission along the slit to obtain a single high signal-to-noise spectrum from a 51 section of the filament (Figs. 2 and 3). Exposures of several spectrophotometric standard stars from Hamuy et al. (1994) were used for flux calibration. We estimate an absolute photometric accuracy of 20% for our quoted emission line fluxes.

## 3. DETECTED EMISSION LINES

Among the most interesting new features in our optical data are the detection of broad H$`\beta `$ and He I $`\lambda `$ 6678. To our knowledge, this is the first detection of the $`\lambda `$6678 line in a pure nonradiative shock. The neutral He atoms passing downstream are unaffected by the electromagnetic turbulence at the shock front. Therefore, neutral He atoms remain cool throughout the shock and produce narrow emission lines unresolved in our data. On the other hand since He+ is an ion, it is heated at the shock front by the same collisionless processes which heat the electrons and ions. Therefore, we expect the He II $`\lambda `$4686 line to be broad and more difficult to detect than the He I $`\lambda `$6678 line. There is a hint of the He II $`\lambda `$4686 line at a low (1.5$`\sigma `$) statistical significance. The surface brightnesses of the detected emission lines from NW SN 1006 appear in Table 1. We have used $`E(BV)=0.11`$ (Schweizer & Middleditch 1980) to obtained dereddened line intensities from our measured spectra.

## 4. SPECTRAL LINE FITS

We fit the line profiles using a standard nonlinear $`\chi ^2`$ reduction routine. We assumed Gaussian profile shapes for all lines and fitted the H$`\alpha `$, He I, H$`\beta `$, H$`\gamma `$ and He II lines separately. Each Balmer line profile was fit with two Gaussians, one each for the narrow and broad components. The He lines were fit with single Gaussians. To reduce the required degrees of freedom, we set the widths of the narrow component Balmer lines and He I line equal to the instrumental resolution during the fitting process.

Due to the faintness of the $`\lambda `$4686 emission, we cannot estimate the width of this line from the fits themselves. Therefore, we set the velocity width of the He II emission line width equal to that of the H$`\alpha `$ line. This assumption is supported by the HUT analysis of Raymond et al. (1995), who found that to within the errors, the velocity widths of the He II $`\lambda `$1640 and broad Ly $`\beta `$ lines are equal (i.e., electron-ion and ion-ion temperature equilibration is inefficient in SN 1006). The best fit for the He II $`\lambda `$4686 line profile appears in Fig. 4. We note that due to the large width and faintness of the He II line, the flux predicted by the line fit is quite sensitive to the baseline level: the lower the baseline, the higher the predicted He II flux. Both the statistical uncertainty and the systematic baseline uncertainty contribute to the large error bars on the flux. For the purpose of our subsequent analysis, the most important quantity here is the upper limit on the He II line flux.

Our best fit for the H$`\alpha `$ broad component width is 2290$`\pm `$80 km$`\mathrm{s}^1`$, agreeing very well with the value of 2310$`\pm `$210 km$`\mathrm{s}^1`$ measured by Smith et al. (1991). The line center shift between the broad and narrow component lines is 29$`\pm `$18 km$`\mathrm{s}^1`$, less than the uncertainty in broad component width, indicating a nearly edge-on viewing angle ($`90^{}\mathrm{\Theta }`$2). The H$`\alpha `$ broad-to-narrow ratio is 0.84$`{}_{0.01}{}^{}{}_{}{}^{+0.03}`$, lying slightly above the extreme upper limit quoted by Smith et al. (1991), 0.73$`\pm `$0.06. The H$`\beta `$ broad-to-narrow ratio is somewhat larger, 0.93$`{}_{0.16}{}^{}{}_{}{}^{+0.18}`$. The H$`\beta `$ broad-to-narrow ratio is larger than that of H$`\alpha `$ due to the lower conversion efficiency of Ly $`\gamma `$ photons into H$`\beta `$ in the narrow component (c.f. Ghavamian et al. 2001). The broad-to-narrow ratio of the H$`\gamma `$ line (0.63$`{}_{0.25}{}^{}{}_{}{}^{+0.53}`$) is poorly constrained due to the faintness of the broad component. The observed H$`\alpha `$ surface brightness of the Balmer filament in SN 1006 is (3.2$`\pm `$0.7)$`\times `$10-16 ergs cm-2 s-1 arcsec-2 (including the 20% photometric uncertainty), in reasonable agreement with the value quoted by KWC87.

## 5. ANALYSIS AND MODELS

To model the Balmer line profiles, we used the one-dimensional, plane-parallel shock code described in Ghavamian et al. (2001) to compute a grid of numerical models. Briefly, the code predicts the density and temperature of electrons, protons and hot neutrals behind a nonradiative shock. Broad to narrow flux ratios are computed by a Monte Carlo simulation which follows collisionally excited Ly $`\beta `$ and Ly $`\gamma `$ photons until they are absorbed or escape from the grid. Collisional excitation and ionization of H by both electrons and protons is included in the models. Free parameters of the model are the preshock neutral fraction $`f_{H^0}`$, shock speed, $`v_S`$, and fractional equilibration immediately behind the shock. The latter quantity denotes the degree of electron-ion equilibration by plasma turbulence at the shock front, where $`f_{eq}`$ = 0 represents no (pure Coulomb) equilibration, while $`f_{eq}`$ = 1 represents full equilibration. The electron and proton temperatures are then defined by the relation

$$T_e=\frac{3}{16}\frac{m_pv_S^2}{k}(\mu f_{eq}+\frac{m_e}{m_p}(1f_{eq}))$$
(1)
$$T_p=\frac{3}{16}\frac{m_pv_S^2}{k}(\mu f_{eq}+(1f_{eq}))$$
(2)
where $`\mu `$ is the mean molecular weight ( = 0.6 for Galactic abundances, Allen 1973) and $`\frac{3}{16}\mu m_pv_S^2`$ is the mean postshock temperature. For an input $`f_{eq}`$, the subsequent evolution of $`T_e`$ and $`T_p`$ by Coulomb collisions is followed as a function of position behind the shock. We have updated the code to include the calculation of both the He ionization balance and the contribution of He ionization to the postshock electron density. To increase the signal-to-noise of our simulations, we have also increased the number of excitations followed by the Monte Carlo simulation from 10,000 used in the earlier paper (Ghavamian et al. 2001) to 50,000.

Using the same method outlined in Ghavamian et al. (2001), we computed a grid of models in equilibration and shock velocity, combining each $`f_{eq}`$ with the $`v_S`$ required to match the observed FWHM of the broad component. We sampled the subset $`f_{eq}`$ = (0, 0.03, 0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 1), corresponding to $`v_S`$ = (2865, 2880, 2895, 2920, 2965, 3100, 3235, 3395, 3580) km$`\mathrm{s}^1`$, respectively. The H$`\alpha `$ and H$`\beta `$ broad-to-narrow ratios predicted by our models are shown in Fig. 5. It is immediately obvious that multiple combinations of $`f_{eq}`$ and $`I_b/I_n`$ are consistent with the observed broad-to-narrow ratios. This makes the unambiguous, simultaneous determination of $`f_{H^0}`$ and $`f_{eq}`$ from Fig. 5 impossible. In an earlier calculation, Laming (2000) obtained similar shapes for the broad-to-narrow ratio curves, the main difference being that our curves flatten out and then decline slightly at high equilibrations. The reason for the disagreement is that unlike our models, the Laming (2000) calculation maintains a constant shock speed while varying $`f_{eq}`$. The slight downturn in our computed narrow ratio curves is caused by the decline of the H+$``$H0 charge exchange rate at high equilibrations (i.e., high shock speeds). As a result, the hot neutral density peaks at progressively smaller values behind the shock as $`v_S`$ increases beyond 3000 km$`\mathrm{s}^1`$.

Fortunately, we can use the He I/He II and He I/H$`\alpha `$ flux ratios to break the degeneracy in Fig. 5. Since every He atom in the postshock ionization zone will spend some time as He+ (the preshock He++ fraction typically amounts to only a few percent), the predicted He II $`\lambda `$4686 flux is insensitive to the preshock He ionization fraction. Furthermore, as demonstrated by Hartigan (1999) and Laming et al. (1996), the He II $`\lambda `$4686 flux is also insensitive to the degree of electron-ion temperature equilibration behind the shock (a result that is confirmed by our models). Therefore, we can use the He I/He II ratio to constrain the preshock He0 fraction. Once this quantity is known, we can then use the He I/H$`\alpha `$ ratio to constrain the preshock H0 fraction. This will enable us to determine the appropriate choice of broad-to-narrow ratio from Fig. 5. From our shock models we estimate that regardless of the assumed equilibration, 90% of the He II emission is produced within a zone $`<`$1 thick (assuming a total preshock density of 1 cm-3 (Winkler & Long 1997), shock speed 3000 km$`\mathrm{s}^1`$ and a distance of 2 kpc). This is smaller than both the pixel scale (0$`\stackrel{}{\mathrm{.}}`$51 pixel-1) and the slit width (2$`\stackrel{}{\mathrm{.}}`$0). Therefore, we believe that all of the He II emission layer was contained within the slit.

In most interstellar shocks where He I emission is observed (mainly radiative shocks), the $`\lambda `$5876 line is brighter than the $`\lambda `$6678 line. The $`\lambda `$6678 emission line arises from the singlet transition $`{}_{}{}^{1}P`$(1s2p) $``$$`{}_{}{}^{1}D`$(1s3d) and is most efficiently excited from the ground state (i.e., by collisional excitation). The $`\lambda `$5876 line, however, arises from a triplet transition and is most efficiently excited by cascades from high lying energy levels (i.e., from recombination). Therefore, since recombination is negligible in SN 1006 shocks, we expect the former line to be significantly brighter than the latter (indeed, $`\lambda `$5876 emission is not detected in our optical spectrum). In the following analysis we limit our calculations to the $`\lambda `$6678 line intensities.

To estimate the He I/H$`\alpha `$ ratio for SN 1006, we computed the total H$`\alpha `$ and He I $`\lambda `$6678 fluxes through the shock front using the combinations of $`f_{eq}`$ and $`v_S`$ listed above. Calculation of the He I $`\lambda `$6678 flux in SN 1006 is relatively straightforward because proton excitation and ionization of He is negligible (Laming et al. 1996). In addition, the only significant contribution to the He I flux other than $`{}_{}{}^{1}S`$(1s2) $``$$`{}_{}{}^{1}D`$(1s3d) excitation are cascades from the 4f level via $`{}_{}{}^{1}D`$(1s3d) $``$$`{}_{}{}^{1}F`$(1s4f). We have included both direct excitation and cascade contributions to 3d in our He I $`\lambda `$6678 flux calculation. We obtained the electron excitation cross sections from the convergent close coupling (CCC) results of Fursa & Bray (1995). At high energies, we approximated the cross sections with a Born-like dependence, scaled to match the CCC values at 0.9 keV. Our He ionization balance calculation uses the electron ionization rates of Younger (1981) for He0 and He+. We also utilized the CCC cross sections of Fursa and Bray (1995) to compute the collisional excitation rate of He II $`\lambda `$4686 electrons (proton excitation is negligible, see Laming et al. 1996). Given the low density, the order of magnitude lower abundance of He relative to H, and the large line width of He II, we assumed Case A conditions in the He II line, i.e., the radiative excitation of He II $`\lambda `$4686 by He II $`\lambda `$243 (Ly $`\gamma `$) absorption is negligible behind the shock.

The model F(6678)/F(4686) ratios are shown in Fig. 6 for a range of preshock He0 fractions and temperature equilibrations. The dashed lines indicate error bars on the observed line ratio after dereddening. Since the ratio of the He I excitation rate to the He I ionization rate depends steeply on electron temperature for $`T_e`$105K, the He I/He II ratio declines sharply at low values of $`f_{eq}`$ ($``$0.01). It is clear from this figure alone that matching the predicted He I/He II ratio with the observed value requires a preshock He neutral fraction $``$0.5 and low electron-ion equilibration.

The predicted He I/H$`\alpha `$ ratios are displayed in Fig. 7 for a range of preshock H neutral fractions. For simplicity, we assumed $`f_{He^0}=`$1 in these models (varying this quantity merely moves the curves up or down). To facilitate comparison with the broad-to-narrow ratios, we used the same set of neutral fractions and combinations ($`f_{eq}`$, $`v_S`$) employed by the models shown in Fig. 5.

(1) Case A in both the broad and narrow components, and (2) Case A in the broad component, Case B in the narrow component. Since the inferred Ly $`\beta `$ and Ly $`\gamma `$ optical depths behind Balmer-dominated shocks are often $``$1 in the narrow component line (Chevalier, Kirshner & Raymond 1980, Ghavamian 1999, Ghavamian et al. 2001), the Case A and Case B line ratio curves mark the largest and smallest values of He I/H$`\alpha `$ attainable for each combination of preshock H and He ionization fraction.

Comparing Figs. 6 and 7, we find that (1) the observed He I/He II and He I/H$`\alpha `$ ratios are simultaneously matched for $`f_{H^0}`$0.1 and $`f_{He^0}`$0.7 in the preshock gas, and (2) only the lowest equilibration models ($`f_{eq}`$0.05) can match the observations. Returning to Fig. 5 and choosing the $`I_b/I_n`$ curves for $`f_{H^0}=`$0.1, we find a match between the observed and modeled H$`\alpha `$ and H$`\beta `$ broad-to-narrow ratios for $`f_{eq}`$0.1, or $`T_e/T_p`$0.07. This corresponds to a shock speed of 2890$`\pm `$100 km$`\mathrm{s}^1`$ (including the measurement uncertainty of the broad component). This result is in good agreement with the shock speed estimated by Laming et al. (1996) from ultraviolet line analysis of HUT data: 2600$`\pm `$300 km$`\mathrm{s}^1`$. It also agrees with the low equilibration inferred by Laming et al. (1996): $`T_e/T_i`$0.05.

Given the intrinsic faintness of the He II emission, it is useful to compare the measured He II $`\lambda `$4686 surface brightness with the value inferred from observation of the He II $`\lambda `$1640 line by HUT (Raymond et al. 1995). Assuming Case A conditions in the He II lines (reasonable due to the low density of He), the predicted He II $`\lambda `$4686 flux in the RC Spectrograph slit is

$$F(4686)=\frac{\nu _{4686}}{\nu _{1640}}\frac{\gamma _{4686}}{\gamma _{1640}}\frac{F_{RCSpec}(H\alpha )}{F_{HUT}(H\alpha )}F_{HUT}(1640)$$
(3)
where $`\gamma _{4686}`$ is the number of $`\lambda `$4686 photons per incident He+ ion and $`\gamma _{1640}`$ is the corresponding quantity for the $`\lambda `$1640 line. We have used our H$`\alpha `$ image (the star-subtracted version of Fig. 1) to scale from the He II $`\lambda `$4686 flux expected in the HUT aperture to the value expected in the RC Spectograph slit. This scaling is valid if the H neutral fractions are roughly equal in the two filaments (reasonable given their close proximity and comparable surface brightnesses). From the H$`\alpha `$ image we estimate $`F_{RCSpec}(H\alpha )/F_{HUT}(H\alpha )`$0.13. Next, from the calculations of Hartigan (2000) we estimate $`\gamma _{4686}`$/$`\gamma _{1640}`$0.1 for 10$`{}_{}{}^{5}KT_e`$107 K. Inserting these values into Eqn. 3 and using $`F_{HUT}(1640)`$ from Raymond et al. (1995, who also assumed $`E(BV)=`$0.11), we find $`F(4686)`$ = (1.1$`\pm `$0.2)$`\times `$10-15 ergs cm-2 s-1. This agrees well with the $`\lambda `$4686 flux we derive from our data, (9.4$`\pm `$6.4)$`\times `$10-16 ergs cm-2 s-1.

## 6. DISCUSSION

In their detailed imaging study of SN 1006, Winkler & Long (1997) found that X-ray emission in the NW peaks approximately 12 behind the optical filaments. Using nonequilibrium ionization models convolved with the ROSAT HRI response function, they found that a total preshock density $``$1 cm-3 was required to match the spatial distribution of the X-ray emission. Combined with our estimated H preshock neutral fraction, this implies a preshock neutral density $`n_{H^0}`$0.1 cm-3 with $`\tau _{Ly\beta }`$0.5 in the narrow component. The inferred value of $`\tau _{Ly\beta }`$ will be even lower if Winkler & Long (1997) overestimated the total preshock density. The dereddened narrow Balmer decrement is 3.35$`\pm `$0.14, significantly lower than the value of 4.3$``$6 observed in Tycho’s SNR (Ghavamian 2001). The absorption of narrow component Lyman photons within the shock structure scales linearly with $`f_{H^0}`$, as long as $`f_{H^0}`$ is small ($``$0.5). This enhances the H$`\alpha `$ emission more than H$`\beta `$, steepening the narrow component Balmer decrement. Therefore, the modest decrement we observe in SN 1006 is consistent with our low derived Lyman optical depth.

### 6.1. The Preshock Ionization Fractions

The detection of the He I $`\lambda `$6678 line indicates that a significant fraction of the preshock He is neutral. The collisional excitation of He II $`\lambda `$304 behind the shock certainly produces some photoionization and heating of the ambient gas (see, for example, Tycho’s SNR, Ghavamian et al. 2000). Since the photoionization rate of He0 by $`\lambda `$304 photons is nearly an order of magnitude higher than the photoionization of H0 (Ghavamian et al. 2000), some additional mechanism is required which can ionize ambient H0 more effectively than He0. There are several possibilities which we discuss in turn below.

#### 6.1.1 Photoionization From the Reverse Shock?

One process which may generate high H ionization fractions in the ambient medium while leaving most He neutral is photoionization by O II$``$O VIII emission from the outermost layer of shocked ejecta (e.g., Hamilton & Fesen 1988, hereafter HF88). The strongest of these lines have energies lying between the H and He0 ionization thresholds. Using the calculations of HF88 and the He II $`\lambda `$304 photoionization models of Ghavamian et al. (2000), we can estimate the ratio of the H photoionization rate to the He photoionization rate in the interstellar medium due to radiation from the blast wave and shocked ejecta. Since each O atom entering the reverse shock will pass through all ionization stages between O++ and O+7, the ratio of photoionizations $`\alpha _{H/He}(R)`$ at a distance $`R`$ from the remnant center is

$$\alpha _{H/He}(R)\frac{_i\gamma (\lambda _i)\sigma ^H(\lambda _i)+\eta \gamma (304)\sigma ^H(304)}{_i\gamma (\lambda _i)\sigma ^{He}(\lambda _i)+\eta \gamma (304)\sigma ^{He}(304)}$$
(4)
where

$`\eta ={\displaystyle \frac{\frac{N_{bw}(He)}{(RR_{bw})^2}}{\frac{N_{ej}(O)}{(RR_{ej}(O))^2}}}`$
Here, $`\gamma (\lambda _i)`$ is the number of ionizing photons per ion of wavelength $`\lambda _i`$ (the Bethe parameter as tabulated by HF88 for the reverse shock), $`\gamma (304)`$ is the Bethe parameter for He II $`\lambda `$304 photons from the blast wave ($``$1.1 for $`v_S>`$2000 km$`\mathrm{s}^1`$, Ghavamian et al. 2000), $`\sigma ^H`$ is the H0 photoionization cross section from Osterbrock (1989) and $`\sigma ^{He}`$ is the corresponding cross section for He0 (Reilman & Manson 1979). The current O-layer and blast wave radii are $`R_{ej}(O)`$ and $`R_{bw}`$, respectively, while $`N_{ej}(O)`$ is the number of shocked O atoms in the ejecta and $`N_{bw}(He)`$ is the number of He atoms swept up by the blast wave. In our calculation we have assumed that the ionization state of the ISM is negligibly affected by recombination. HF88 found that the relative masses of Fe, Si and O predicted by carbon deflagration model CDTG7 (Woosley & Weaver 1987) matched the observed Fe II $`\lambda \lambda `$2586, 2599 absorption profile in SN 1006 best. Therefore, using the shocked O ejecta mass of 0.41 $`M_{\mathrm{}}`$ and blast wave mass of 4.3$`M_{\mathrm{}}`$ from model CDTG7, $`R_{ej}(O)=`$6 pc and $`R_{bw}=`$8.3 pc from HF88 and neglecting radiative transfer effects on the photoionizing radiation, we find $`\alpha _{H/He}`$0.1 for a parcel of interstellar gas just outside the blast wave radius. This result suggests that reverse shock photoionization cannot produce our derived low H neutral fraction. The smallness of $`\alpha _{H/He}`$ is mainly due to that fact that (1) there are nearly an order of magnitude more He atoms in the blast wave than there are O atoms in the reverse shock, and (2) since $`R_{ej}(O)<R_{bw}`$, a given parcel of interstellar gas sees a larger emitting area from the blast wave than from the reverse shock.

A more detailed, time dependent calculation would be desirable. However, we note that the reverse shock has already propagated through the entire O-layer, while the blast wave continues to accumulate interstellar material (and hence He atoms). Therefore, the ratio $`M_{ej}(O)/M_{bw}`$ is decreasing as the supernova remnant evolves and the blast wave is becoming the increasingly dominant photoionization source for the ISM around SN 1006. Finally, we note that even if we had included photoionization from the Si-rich and Fe-rich layers in Eqn. 4, our conclusions would remain unchanged since most of the EUV/X-ray lines from Si and Fe will ionize He much more effectively than H.

#### 6.1.2 Relic Ionization from the Supernova Explosion?

An interesting source of photoionization is the EUV flash produced when photons escape from the detonating white dwarf (Type Ia SN). Using the EUV burst models presented by Blinnikov & Sorokina (2000) and Sorokina & Blinnikov (2001, in preparation), we estimate $``$10$`{}_{}{}^{55}`$1056 ionizing photons (E$`>`$13.6 eV) are released over the integrated Type Ia SN light curve. The precise number can vary by a factor $``$5, depending on which detonation model is used in the light curve calculation (see Blinnikov & Sorokina (2000) for model descriptions). Due to the steepness of the EUV spectrum, the H0 ionization rate is much larger than the He0 ionization rate, a desirable feature. However, assuming an average ISM density of 0.06$``$1.0 cm-3 around SN 1006 and distance of 2 kpc, there were $``$5$`\times `$1057 H atoms contained within the current blast wave radius at the time of the SN explosion. Even if these H atoms were already partially ionized beforehand, the total ionizing flux from a Type Ia explosion still falls an order of magnitude short of the value needed to effectively ionize the gas at the current blast wave radius. The disparity is worsened if we include radiative transfer loss of the released photons as they (1) propagate through the ISM, and (2) possibly encounter a circumstellar envelope left by a nondegenerate donor star prior to the explosion (Cumming et al. 1996). Therefore, regardless of whether the EUV emission ionizes H more effectively than He, there are not enough photons produced in the explosion to explain our inferred neutral fractions. This result also agrees with the SN flash ionization calculations of Cheng & Bruhweiler (1996).

#### 6.1.3 Photoionization from Galactic Background Radiation?

If neither the supernova event nor the postshock emission from the supernova remnant can account for our deduced H0 and He0 fractions, we must consider nonlocal sources of ionizing radiation. Given the lack of late O and early B stars in the vicinity of SN 1006, the most promising source is the Galactic ionizing background. Reynolds (1984) observed the H$`\alpha `$ recombination line toward a sample of high Galactic latitude pulsars ($`b>`$5) at distances of 2$``$3 kpc. He found that Lyman continuum radiation from O stars in the Galactic plane was the most likely source of ionization in the diffuse ISM. Since the Lyman continuum is produced by recombination, it declines nearly exponentially beyond each ionization edge (Brussard & Van De Hulst 1962). Therefore, incident Lyman continuum photons will ionize H far more effectively than He, a desirable property.

The remaining question is whether the Galactic Lyman continuum flux is high enough at the position of SN 1006 ($`b=`$+14.6) to produce the necessary amount of H ionization. Assuming a density of 1 cm-3 around SN 1006, the recombination rate for 5000 K$``$T$``$10,000 K is $``$5$`\times `$10-13 s-1 (Osterbrock 1989). Therefore, obtaining an H neutral fraction $``$0.1 requires a photoionization rate $``$5$`\times `$10-12 s-1. Dove, Shull & Ferrara (2000) estimated a Lyman continuum flux $``$5$`\times `$107 photons cm-2 s-1 from O stars in the Galactic disk. In their time dependent radiative transfer calculation for the ionizing continuum, they estimated that roughly 3%$``$6% of the Lyman continuum photons escape through superbubbles and Galactic chimneys into the halo of the Milky Way, producing the “Reynolds Layer” of diffuse H$`\alpha `$ emission (Reynolds 1984). If the Lyman photons propagating into the halo suffer no further attenuation, the corresponding H ionization rate in the vicinity of SN 1006 is approximately $`F_{Ly}\sigma ^H`$(912) $``$1$``$2$`\times `$10-11 photons s-1, more than enough to produce the necessary H ionization around SN 1006.

If the Galactic Lyman continuum background pre-ionizes most of the ambient H on the NW side of SN 1006, the question arises of why the ambient H on the eastern half of Tycho’s SNR (also the remnant of a Type Ia explosion) was found to be mostly neutral (Ghavamian et al. 2000). First, the eastern half of Tycho’s SNR is likely interacting with the outer edge of a warm H I cloud (Reynoso et al. 1999, Ghavamian et al. 2000). Second, the presence of large structures in the Galactic disk (H I clouds, superbubbles, molecular cloud complexes, etc.) produces large variations in the Lyman continuum flux from position to position and time to time (Dove, Shull & Ferrara 2000). Since Tycho’s SNR is close to the Galactic plane ($`b=`$+1.4), it is possible that this remnant is shielded from the Lyman continuum emission by the many dense clouds in the Galactic disk. SN 1006, on the other hand, lies nearly 500 pc above the Galactic plane. Therefore, neutral H near SN 1006 may “see” a larger time averaged Lyman continuum than the neutral H around Tycho’s SNR.

### 6.2. Comparison With Ionization Estimates of the Local Bubble

An interesting comparison can be drawn between our derived H and He neutral fractions and values estimated for the local ($``$200 pc) ISM. For example, using EUVE spectra of 29 DA white dwarfs, Wolff et al. (1999) estimated $`f_{He^+}`$ 0.4 in the Local Bubble, with no evidence of a gradient or other general pattern. They interpreted this result as evidence for a lack of He-ionizing sources in any particular location. On the other hand, Wolff et al. (1999) found a positive gradient in the H ionization fraction toward the Canis Major cavity, implicating the B-stars Adara and Mirzam as the photoionizing sources. In some respects the ionization properties of the Local Bubble and the ISM near SN 1006 are similar: the H ionization is strongly influenced by EUV radiation, while He remains largely unaffected. The main difference between the two cases, of course, is the lack of a local H-ionizing source near SN 1006. Our derived He ionization fraction ($``$0.3) for SN 1006 is lower than the average value of 0.4 measured by Wolff et al. (1999), but there are individual cases from their WD sample with He ionization fractions similar to ours.

### 6.3. Comparison With Earlier Results

Our optical analysis lends strong support to the picture developed by earlier X-ray, UV and optical studies of SN 1006. Multi-epoch optical imagery of SN 1006 has yielded a proper motion of 0.3$`\pm `$0.04 yr-1 (Long, Blair & van den Bergh 1988) for the Balmer filaments in the NW. Combined with our estimated shock speed this gives a distance of 2.0$`{}_{0.3}{}^{}{}_{}{}^{+0.4}`$ kpc, agreeing well with the value 1.8$`\pm `$0.3 obtained by Laming et al. (1996). The SN 1006 blast wave is the fastest shock we have thus far modeled with our Balmer-dominated shock code. In modeling Balmer spectra from the Cygnus Loop ($`v_S`$300 km$`\mathrm{s}^1`$), RCW 86 ($`v_S`$600 km$`\mathrm{s}^1`$) and Tycho’s SNR ($`v_S`$ 2000 km$`\mathrm{s}^1`$), we noted that the equilibration $`f_{eq}`$ declined from near 1 for the slowest shocks to $``$0.2 for the fastest shocks (Ghavamian 1999, Ghavamian et al. 2001). The shocks in SN 1006 are closer to 3000 km$`\mathrm{s}^1`$ and are matched by an even lower $`f_{eq}`$, further evidence of an inverse correlation between shock speed and equilibration. The case for low equilibration at high Mach numbers has recently received further support from X-ray observations of SN 1987A. In their analysis of high resolution HETG data from Chandra, Michael et al. (2001) found that a high blast wave speed (3500 km$`\mathrm{s}^1`$) along with a low electron temperature ($`f_{eq}<`$0.1) were required to match the X-ray profiles of the SN 1987A spectrum.

### 6.4. Implications for Collisionless Shock Physics

There are few numerical simulations available of high Mach number, collisionless shocks. Cargill & Papadapoulos (1988) modeled the electron-ion heating in a shock moving perpendicular to the ambient magnetic field, with $`M_A=`$50 ($`M_Av_S/(c_s^2+v_A^2)^{1/2}`$ is the Alfvénic Mach number, where $`c_s`$ is the sound speed and $`v_A`$ is the Alfvén speed of the preshock gas). In their hybrid simulation the electrons were treated as a resistive fluid, while the ions were treated as individual particles. Cargill & Papadapoulos (1988) found that two-stream plasma instabilities produced by ions reflected upstream resulted in $`T_e/T_i`$0.25 at the shock front. Shimada & Hoshino (2000) tried an approach where both electrons and ions were treated as particles. At the highest Mach numbers covered by their simulations, they predicted that the equilibration would increase with Mach number. In their fastest simulated shock ($`M_A`$ = 20), they found $`T_e/T_i`$ 0.2. Their curve of $`T_e/T_i`$ vs. $`M_A`$ flattens considerably as $`M_A`$20; extrapolating this curve to $`M_A`$200 leads to significantly larger equilibration than implied by the ultraviolet (Laming et al. 1996) and optical (this paper) observations of SN 1006.

During the last few years, a considerable effort has been made in modeling collisionless heating processes in solar wind shocks ($`M_A`$20). Hull & Scudder (2000) combined an electron fluid treatment with constraints from ISEE 1 data of Earth’s bow shock to model the electron-ion equilibration process. For the first time, they were able to match the observed electron temperature jumps in fast shocks propagating both parallel and perpendicular to the magnetic field. These models emphasize the role of coherent (DC) electric forces in the collisionless heating, rather than wave-particle interactions. The advantage of the Hull & Scudder (2000) treatment is that it draws upon in situ observations of real collisionless shocks rather than computer simulations. Both these models and the observations (Schwartz et al. 1988) predict $`T_e/T_i`$0.25 for shocks in the solar wind. Although we are unable to perform similar detailed observations of supernova shocks, extending the models of Hull & Scudder (2000) to much higher Mach numbers may be an interesting direction for future work on collisionless shocks.

## 7. SUMMARY

The optical spectrum presented here is the deepest yet obtained of the nonradiative filament along the northwest rim of SN 1006. Lines of H$`\alpha `$, H$`\beta `$, H$`\gamma `$, H$`\delta `$, He I $`\lambda `$6678 and possibly He II $`\lambda `$4686 have been detected$``$more lines than in the spectrum of any other purely nonradiative SNR filament. We have detected both broad and narrow component emission lines of H$`\alpha `$, H$`\beta `$ and H$`\gamma `$. The broad and narrow components result, respectively, from fast neutrals produced by charge exchange and slow neutrals that are collisionally excited.

The broad component widths, the ratios of broad to narrow flux in H$`\alpha `$ and H$`\beta `$, and the Balmer line strengths relative to He I and He II enable us to determine the important parameters for the shock. We find a shock speed of 2890$`\pm `$100 km$`\mathrm{s}^1`$ and a very low degree of electron-ion equilibration, $`T_e/T_p`$0.07, both consistent with earlier measurements. In addition, our models can match the observations only if the preshock H is mostly ionized ($``$90%) and preshock He is mostly neutral ($``$70%). In seeking an explanation for this ionization structure, we have concluded that neither radiation from the blast wave nor radiation from the reverse shock can effectively ionize the ambient H surrounding SN 1006. We also find that the EUV flash from the supernova explosion is incapable of producing the necessary H ionization. The most plausible explanation for the high H ionization fraction seems to be photoionization from the Galactic Lyman continuum. Finally, comparing our derived electron-ion equilibration with values derived for other nonradiative SNRs reinforces the earlier conclusion that the fastest shocks have the least equilibration. A thorough understanding of this trend remains as a challenge for the study of collisionless shock physics.

The authors gratefully acknowledge the support staff and general hospitality of CTIO, where various parts of this work were completed. We would like to thank P. Lundqvist for helpful correspondence on the EUV bursts from SNe and E. I. Sorokina for providing us with EUV fluxes from her emission calculations of Type Ia SNe in advance of publication. P. G. would like to thank J. P. Hughes for valuable discussions, and Marc Hemsendorf for valuable advice on optimizing the codes used in this paper. The authors wish to acknowledge their sources of financial support: P. G., through Chandra grants GO0-1035X and GO1-2052X to Rutgers University; P. F. W., through NSF grant AST-9618465 and NASA grant NAG 5-8020; J. C. R., through FUSE grant NAG5-10352; and K. S. L., through Chandra grants GO0-1120X and GO1-2058A.


text:
# 1 Introduction

## 1 Introduction

In a unified treatment of the semi-leptonic electroweak interactions in nuclei, like lepton-nucleus scattering , $`\beta `$-decay modes , charged-lepton nuclear capture ($`\mu ^{}`$ capture, etc.) , neutrino-nucleus induced reactions , exotic semi-leptonic processes in nuclei (double beta-decay , $`\mu ^{}e^{}`$ conversion, etc. ), at first, one takes advantage of the well known electromagnetic interactions to probe the nucleus and determine accurately the nuclear charge, convection current and spin magnetization distributions . Afterwards, having the nuclear structure uncertainties reduced to a minimum, one can predict the cross sections for a variety of semi-leptonic weak interaction processes.

As is well known, the reaction rates $`\mathrm{\Gamma }_{if}`$ between the initial $`|i`$ and a final $`|f`$ discrete nuclear states of any of the aforementioned processes is written in terms of the matrix elements of a specific effective Hamiltonians $`H_{eff}`$ as

$$\mathrm{\Gamma }_{if}\left|f|H_{eff}|i\right|^2$$

where the $`H_{eff}`$ can be assumed to be constructed within the current-current interaction hypothesis by first writing down the leptonic ($`j_\lambda `$) and hadronic ($`J_\lambda `$) currents.

The effective hadronic current-density operator $`\widehat{𝒥}_\lambda `$ which enters the $`H_{eff}`$, contains only strong isoscalar ($`T=0`$) and isovector ($`T=1`$) pieces and, in general, it has both vector (V) and axial vector (A) pieces as

$$[\widehat{𝒥}_\lambda ]^{TM_T}=\beta _V^{(T)}[\widehat{J}_\lambda ]^{TM_T}+\beta _A^{(T)}[\widehat{J}_\lambda ^5]^{TM_T}.$$
(1)
where $`\beta _V^{(T)}`$ and $`\beta _A^{(T)}`$ are polar-vector and axial-vector coefficients involving the relevant couplings (the purely electromagnetic interactions involve only vector current density, but the weak interactions involve in addition axial vector current). In neutral-current processes $`M_T=0`$ for both isoscalar and isovector parts, while in charge changing interactions $`T=1`$ with $`M_T=\pm 1`$. One, usually, assumes that for a nuclear target (denoted as $`(A,Z)`$ with Z, A the proton- and mass-number respectively) the states $`|i`$ and $`|f`$ are characterized by definite angular momentum and parity ($`|\alpha |J^\pi `$). In neutral-current reactions (lepton scattering, $`\mu e`$ conversion, etc.) the final state $`|f`$ is, in general, an excited state of the nucleus $`(A,Z)`$ while in charged-current processes $`|f`$ denotes an excited state of the $`(A,Z1)`$ nucleus.

The nuclear calculations of $`f|H_{eff}|i`$ proceed usually by decomposing the Fourier transform of the hadronic current-density operator into irreducible tensors in strong isospin space . Then, spherical Bessel functions $`j_l(x)`$ in central or non-central interaction components are obtained. The nuclear transition rates are subsequently determined by the matrix elements of these operators. The possibility to obtain closed expressions for the matrix elements of the basic one-body operators describing the above mentioned reactions has been examined in previous theoretical studies with special emphasis on the use of harmonic oscillator wave functions .

In Refs. - a method was developed which provides explicit compact expressions for some basic nuclear form factors (those needed for the $`\mu ^{}e^{}`$ conversion in nuclei, scattering of cold dark matter candidates off nuclei, etc.). The advantage of this method is that it allows the easy calculation of the required reduced matrix elements by separating the geometrical coefficients from the kinematical parameters of the reaction. In this way, the transition matrix elements for every value of the momentum transfer q can be evaluated immediately. It is worth remarking that such a compact formalism provides a very useful insight to those authors who wish to adopt a phenomenological approach and fit the nuclear transition strengths.

In the present work, we extend the method of Ref. - so as to provide expressions for the reduced matrix elements of any of the seven basic single-particle tensor operators describing the above-mentioned semi-leptonic processes . We especially deal with the operators involving the differential operator $``$ (nabla), which have not been studied previously . To this end, we first treat the corresponding radial matrix elements which contain derivatives with respect to the coordinate r. Finally, we discuss some representative examples as applications of our formalism focusing on the improvement of the previous expressions giving the radial part of the nuclear moments in the neutrinoless double-beta decay -.

## 2  Multipole decomposition of the hadronic currents

From a nuclear physics point of view, the hadronic current ($`J_\mu `$) which enters the weak and electromagnetic interaction Hamiltonian $`H_{eff}`$ is of primary interest. In our convention the four-current operator $`\widehat{J}_\mu `$ is written as

$$\widehat{J}_\mu (𝐫)=(\widehat{\rho }(𝐫),\widehat{𝐉}(𝐫))$$
(2)
where $`\widehat{\rho }(𝐫)`$ the density and $`\widehat{𝐉}(𝐫)`$ the three-current operators. The standard multipole expansion procedure (see Appendix A) , applied on the matrix elements of $`\widehat{J}_\mu (𝐫)`$, leads to spherical tensor operators which are given in terms of projection functions involving spherical Bessel functions $`j_L(r)`$ and spherical Harmonics $`Y_M^L(\widehat{r})`$ or vector spherical Harmonics $`𝐘_M^{(L,1)J}(\widehat{r})`$ as:

$$M_M^J(q𝐫)=\delta _{LJ}j_L(qr)Y_M^L(\widehat{r}),$$
(3)
$$𝐌_M^{(L1)J}(q𝐫)=j_L(qr)𝐘_M^{(L1)J}(\widehat{r}),$$
(4)
where

$$𝐘_M^{(L1)J}(\widehat{r})=\underset{m,q}{}Lm1q|JMY_m^L(\widehat{r})\widehat{e}_q,$$
(5)
$`\widehat{e}_q`$ is a unit vector in the direction of the three-momentum transfer q. The magnitude of $`𝐪`$ ($`q=|𝐪|`$) is given from the kinematics of the studied process (see e.g. Ref. ).

Using the projections functions Eqs. (3) and (4), in the case of the polar-vector current ($`\widehat{J}_\lambda `$) we obtain the multipole operators

$$\widehat{M}_{JM;TM_T}^{coul}=𝑑𝐫M_M^J(q𝐫)\widehat{\rho }(𝐫)_{TM_T},$$
(6)
$$\widehat{L}_{JM;TM_T}=i𝑑𝐫\left(\frac{1}{q}M_M^J(q𝐫)\right)\widehat{𝐉}(𝐫)_{TM_T},$$
(7)
$$\widehat{T}_{JM;TM_T}^{el}=𝑑𝐫\left(\frac{1}{q}\times 𝐌_M^{JJ}(q𝐫)\right)\widehat{𝐉}(𝐫)_{TM_T},$$
(8)
$$\widehat{T}_{JM;TM_T}^{mag}=𝑑\mathrm{𝐫𝐌}_M^{JJ}(q𝐫)\widehat{𝐉}(𝐫)_{TM_T}.$$
(9)
which have well-defined parity and are called as: Coulomb, longitudinal, transverse-electric, and transverse-magnetic multipoles, respectively. The first three operators have parity $`()^J`$ (normal parity operators), while the parity of $`\widehat{T}_J^{mag}`$ is $`()^{J+1}`$ (abnormal parity operator). The analogous axial-vector multipoles are

$$\widehat{M}_{JM;TM_T}^5=𝑑𝐫M_M^J(q𝐫)\widehat{\rho }(𝐫)_{TM_T}^5,$$
(10)
$$\widehat{L}_{JM;TM_T}^5=i𝑑𝐫\left(\frac{1}{q}M_M^J(q𝐫)\right)\widehat{𝐉}(𝐫)_{TM_T}^5,$$
(11)
$$\widehat{T}_{JM;TM_T}^{el5}=𝑑𝐫\left(\frac{1}{q}\times 𝐌_M^{JJ}(q𝐫)\right)\widehat{𝐉}(𝐫)_{TM_T}^5,$$
(12)
$$\widehat{T}_{JM;TM_T}^{mag5}=𝑑\mathrm{𝐫𝐌}_M^{JJ}(q𝐫)\widehat{𝐉}(𝐫)_{TM_T}^5.$$
(13)
The first three axial-vector multipoles have parity $`()^{J+1}`$ while $`\widehat{T}_J^{mag5}`$ is a normal parity operator. For a conserved vector current (CVC) like the electromagnetic, the longitudinal multipoles can be written in terms of the Coulomb ones as

$$\widehat{L}_{JM_J;Tm_T}(q)=\frac{q_0}{q}\widehat{M}_{JM_J;TM_T}(q)$$
(14)
($`q_0`$ is the time component of the four-momentum transfer, $`q_\mu =(q_0,𝐪)`$). In this case, the number of independent operators resulting from the projection procedure is reduced to seven.

The isospin dependence of the operators (6)-(13) includes the operator $`I_T^{M_T}`$ given by

$$I_T^{M_T}=\{\begin{array}{cc}1,T=0,M_T=0\hfill & \\ \tau _0=\tau _3,T=1,M_T=0\hfill & \\ \tau _\pm =\frac{1}{\sqrt{2}}(\tau _1\pm \tau _2),T=1,M_T=\pm 1\hfill & \end{array}$$
(15)
The exact form of $`I_T^{M_T}`$ is determined by the specific semi-leptonic reaction in question (see the Introduction).

The matrix elements of the seven basic operators of Eqs. (6)-(13) involve isospin dependent form factors $`F_X^{(T)}`$ (see Appendix A). For this reason, we define seven new operators of which the spin-space parts are written (in the first quantization) as :

$`T_1^{JM}M_M^J(q𝐫)=\delta _{LJ}j_L(\rho )Y_M^L(\widehat{r}),`$ (16)
$`T_2^{JM}\mathrm{\Sigma }_M^J(q𝐫)=𝐌_M^{JJ}𝝈,`$ (17)
$`T_3^{JM}\mathrm{\Sigma }_{}^{}{}_{M}{}^{J}(q𝐫)=i\left\{{\displaystyle \frac{1}{q}}\times 𝐌_M^{JJ}(q𝐫)\right\}𝝈,`$ (18)
$`T_4^{JM}\mathrm{\Sigma }_{}^{\prime \prime }{}_{M}{}^{J}(q𝐫)=\left\{{\displaystyle \frac{1}{q}}M_M^J(q𝐫)\right\}𝝈,`$ (19)
$`T_5^{JM}\mathrm{\Delta }_M^J(q𝐫)=𝐌_M^{JJ}(q𝐫){\displaystyle \frac{1}{q}},`$ (20)
$`T_6^{JM}\mathrm{\Delta }_{}^{}{}_{M}{}^{J}(q𝐫)=i\left\{{\displaystyle \frac{1}{q}}\times 𝐌_M^{JJ}(q𝐫)\right\},`$ (21)
$`T_7^{JM}\mathrm{\Omega }_M^J(q𝐫)=M_M^J(q𝐫)𝝈{\displaystyle \frac{1}{q}}.`$ (22)
(hereafter we use the shorter notation for the projection functions of Eqs. (4) and (5) as $`𝐌_M^{LJ}`$ and $`𝐘_M^{LJ}`$, respectively). For the reader’s convenience, in Eqs. (16)-(22) we adopt the usual notation . Using properties of the nabla operator ($``$), Eqs. (18), (19) and (21) can be rewritten as:

$$T_3^{JM}\mathrm{\Sigma }_{}^{}{}_{M}{}^{J}=[J]^1\left\{J^{1/2}𝐌_M^{J+1J}+(J+1)^{1/2}𝐌_M^{J1J}\right\}𝝈,$$
(23)
$$T_4^{JM}\mathrm{\Sigma }_{}^{\prime \prime }{}_{M}{}^{J}=[J]^1\left\{(J+1)^{1/2}𝐌_M^{J+1J}+J^{1/2}𝐌_M^{J1J}\right\}𝝈,$$
(24)
$$T_6^{JM}\mathrm{\Delta }_{}^{}{}_{M}{}^{J}=[J]^1\left\{J^{1/2}𝐌_M^{J+1J}+(J+1)^{1/2}𝐌_M^{J1J}\right\}\frac{1}{q},$$
(25)
(throughout this work we use the common symbol $`[J]=(2J+1)^{1/2}`$). By glancing at Eqs. (16)-(25), we conclude that, the seven basic single-particle operators of definite parity are built with the aid of four operators

$$𝒪_1^{JM}=M_M^J,𝒪_2^{JM}=𝐌_M^{LJ}𝝈,$$
(26)
$$𝒪_3^{JM}=𝐌_M^{LJ}\frac{1}{q},𝒪_4^{JM}=M_M^J𝝈\frac{1}{q}.$$
(27)
which contain the projection functions $`M_M^J`$ and $`𝐌_M^{LJ}`$ or their products with the nucleon-spin $`𝝈`$ and/or the $``$ operator.

## 3 The Single-particle matrix elements

Many quantities of the semi-leptonic electroweak processes are expressed (to a good approximation) in terms of single-particle nuclear matrix elements of the one-body operators $`T_i^{JM},i=1,2,\mathrm{}7`$ . These matrix elements, by applying the Wigner-Eckart theorem, are written as:

$$j_1m_1|T_i^{JM}|j_2m_2=()^{j_1m_1}\left(\begin{array}{ccc}j_1& J& j_2\\ m_1& M& m_2\end{array}\right)j_1T_i^Jj_2,$$
(28)
In our notation, the single-particle wave functions are labelled (see Appendix B) with the quantum numbers ($`nlsjm_j`$). The reduced matrix elements of the operators $`𝒪_i^J`$, $`i=1,2,3,4`$, after applying the re-coupling relations, can be written in the closed forms shown below.

1. For the operators $`𝒪_i^J`$, with $`i=1,2`$, the reduced matrix elements $`j_1𝒪_i^{(L,S_i)J}j_2`$ have been compactly written as

$$j_1𝒪_i^{(L,S_i)J}j_2=(l_1Ll_2)𝒰_{LS}^Jn_1l_1|j_L(\rho )|n_2l_2,i=1,2$$
(29)
where the symbols $`(l_1Ll_2)`$ and $`𝒰_{LS}^J`$ are

$$(l_1Ll_2)()^{l_1}\frac{1}{\sqrt{4\pi }}[l_1][L][l_2]\left(\begin{array}{ccc}l_1& L& l_2\\ 0& 0& 0\end{array}\right),$$
(30)
$$𝒰_{LS}^J[j_1][j_2][J](S+1)^{1/2}(S+2)^{1/2}\left\{\begin{array}{ccc}l_1& l_2& L\\ 1/2& 1/2& S\\ j_1& j_2& J\end{array}\right\}.$$
(31)
In the case of the Fermi-type operator $`𝒪_1^J`$, we must put $`S=0`$ in the 9-j symbol of Eq. (31) while in the case of the Gamow-Teller operator $`𝒪_2^J`$, $`S=1`$.

2. The reduced matrix elements of $`𝒪_3^J`$ after some manipulation can be cast in the form

$$j_1𝐌^{LJ}(q𝐫)\frac{1}{q}j_2=\underset{\alpha }{}𝒜_L^\alpha (j_1j_2;J)n_1l_1|\theta _L^\alpha (\rho )|n_2l_2,\alpha =\pm $$
(32)
where the coefficients $`𝒜_L^\pm `$ are

$`𝒜_L^\pm (j_1j_2;J)`$ $`=`$ $`\pm ()^{l_1+L+j_2+1/2}[j_1][j_2][J]\left({\displaystyle \frac{2l_2+11}{2}}\right)^{1/2}(l_1Ll_21)`$ (33)
$`\times 𝒲_6(l_1,j_1,1/2,j_2,l_2,J)𝒲_6(L,1,J,l_2,l_1,l_21),`$
with $`𝒲_6`$ representing the common 6-j symbol

$$𝒲_6(l_1,j_1,1/2,j_2,l_2,J)\left\{\begin{array}{ccc}l_1& j_1& 1/2\\ j_2& l_2& J\end{array}\right\}.$$
(34)
3. Similarly, for the reduced matrix element of $`𝒪_4^J`$ we can write

$$j_1M^J(q𝐫)𝝈\frac{1}{q}j_2=\underset{\alpha }{}_L^\alpha (j_1j_2;J)n_1l_1|\theta _J^\alpha (\rho )|n_2l_2,\alpha =\pm $$
(35)
where

$`_L^\pm (j_1j_2;J)`$ $`=`$ $`\pm \delta _{j_2,l_21/2}[j_1][j_2](l_1J2j_2l_2)𝒲_6(l_1,j_1,1/2,j_2,2j_2l_2,J).`$ (36)
From Eqs. (29), (32) and (35) we notice that all basic single-particle reduced matrix elements required for our purposes rely on the following three types of radial integrals:

$$n_1l_1|\theta _l^\alpha (\rho )|n_2l_2𝑑rr^2R_{n_1l_1}^{}(r)\theta _l^\alpha (\rho )R_{n_2l_2}(r),\alpha =\mathrm{\hspace{0.17em}0},\pm $$
(37)
with

$$\theta _l^0(\rho )=j_l(\rho ),\theta _l^\pm (\rho )=j_l(\rho )\left(\frac{d}{d\rho }\pm \frac{2l_2+1\pm 1}{2\rho }\right).$$
(38)
The argument $`\rho `$ in Eqs. (38) is equal to $`\rho =qr`$. We note that the derivatives with respect to $`\rho `$ appeared in the radial matrix elements come from the application of the gradient formula in the matrix elements of Eqs. (32) and (35).

### 3.1 Expressions for the Radial Integrals

It is well known that, for single-particle wave functions with arbitrary radial dependence it is not easy to perform analytically the integrations over r in Eq. (37). These integrals, however, can be simplified in the case when harmonic oscillator basis (see Appendix B) is used . Then they take the elegant expressions shown below.

(i) For the operator $`\theta _l^0(\rho )`$ in Ref.  it was proved that

$$n_1l_1|j_L(\rho )|n_2l_2=e^yy^{L/2}\underset{\mu =0}{\overset{n_{max}}{}}\epsilon _\mu ^Ly^\mu ,y=(qb/2)^2$$
(39)
$$n_{max}=(N_1+N_2L)/2,$$

where $`N_i=2n_i+l_i`$ represent the harmonic oscillator quanta of the $`i_{th}`$ level. The coefficients $`\epsilon _\mu ^L(n_1l_1n_2l_2)`$ are given by

$$\epsilon _\mu ^L(n_1l_1n_2l_2)=G\frac{\pi ^{1/2}}{2}\underset{m_1=\varphi }{\overset{n_1}{}}\underset{m_2=\sigma }{\overset{n_2}{}}n!\mathrm{\Lambda }_{m_1}(n_1l_1)\mathrm{\Lambda }_{m_2}(n_2l_2)\mathrm{\Lambda }_\mu (nL),$$
(40)
with

$$n=m_1+m_2+(l_1+l_2L)/2.$$

Also, $`G=b^3N_{n_1l_1}N_{n_2l_2}/2`$, where $`N_{n_1l_1}`$ is defined in Appendix B. The other symbols of Eq. (40) are explained in Ref. .

(ii) The formulation of the radial matrix elements which include differential operators $`\theta _l^\pm `$, which constitute the main task of our present paper, proceeds in a similar manner to that of Eq. (39). This leads to the expressions

$$n_1l_1|j_L(\rho )\left(\frac{d}{d\rho }\pm \frac{2l_2+1\pm 1}{2\rho }\right)|n_2l_2=e^yy^{(L1)/2}\underset{\mu =0}{\overset{n_{max}}{}}\zeta _\mu ^\pm (L)y^\mu ,$$
(41)
where the geometrical coefficients $`\zeta _\mu ^\pm (n_1l_1n_2l_2;L)`$ are given in terms of those of Eq. (40), i.e.

$$\zeta _\mu ^{}(L)=\frac{1}{2}\{\begin{array}{cc}(n_2+l_2+3/2)^{1/2}\epsilon _\mu ^L(n_1l_1n_2l_2+1)+n_2^{1/2}\epsilon _\mu ^L(n_1l_1n_21l_2+1),\hfill & 0\mu <n_{max}\hfill \\ (n_2+l_2+3/2)^{1/2}\epsilon _{n_{max}}^L(n_1l_1n_2l_2+1),\mu =n_{max}\hfill & \end{array}$$
(42)
$$\zeta _\mu ^+(L)=\frac{1}{2}\{\begin{array}{cc}(n_2+l_2+1/2)^{1/2}\epsilon _\mu ^L(n_1l_1n_2l_21)+(n_2+1)^{1/2}\epsilon _\mu ^L(n_1l_1n_2+1l_21),\hfill & 0\mu <n_{max}\hfill \\ (n_2+1)^{1/2}\epsilon _{n_{max}}^L(n_1l_1n_2+1l_21),\mu =n_{max}\hfill & \end{array}$$
(43)
The value of the index $`n_{max}`$ in the latter two cases of radial integrals is determined by

$$n_{max}=(N_1+N_2L+1)/2.$$

At this point, it should be noted that, the explicit and general formulas of Eqs. (40), (42) and (43) hold for every combination of the levels $`(n_1l_1)j_1`$, $`(n_2l_2)j_2`$, and give the geometrical (momentum independent) coefficients of the polynomials defined in Eqs. (39) and (41). These closed expressions have been constructed by inverting properly the multiple summations involved in the corresponding matrix elements (see Ref. ), so as the final summation is performed over the harmonic oscillator quantum number $`N=2n+l`$.

In order to show the advantages of the above formalism we discuss below some applications.

## 4 Applications

As mentioned before, our formalism is applicable for an infinite single-particle basis and can be applied to any semi-leptonic electroweak process which take place in the field of nuclei . In a realistic case one is forced to truncate the basis set. Thus, e.g. the description of the nuclear moments for the neutrinoless double beta decay $`{}_{}{}^{48}Ca^{48}`$Ti can be done using a model space including the seven orbitals of the major shells 2$`\mathrm{}\omega `$ and 3$`\mathrm{}\omega `$ .

In Tables 1-3 we list the coefficients $`\epsilon _\mu ^L`$ and $`\zeta _\mu ^\pm `$ needed to evaluate the radial matrix elements $`n_1l_1|\theta _L^\alpha (\rho )|n_2l_2`$ in the above (1s-0d-1p-0f) model space. For some additional geometrical coefficients $`\epsilon _\mu ^L`$, which determine the transition matrix elements for scattering of dark matter candidates off nuclei in the model space including the orbitals 2p-0h and 0i-1g, the reader is referred to Ref. .

By exploiting the analytic expressions of the Secs. 2 and 3, in the next subsection we simplify the formalism of Ref.  constructed for the calculation of the nuclear moments in neutrinoless double-beta decay.

### 4.1 Nuclear moments for the neutrinoless double-beta decay

In the description of the nuclear moments for the neutrinoless double-beta decay - the following radial integrals are needed:

$$_L^\kappa (n_1l_1n_2l_2;q)=n_1l_1|j_L(qr)r^\kappa |n_2l_2$$
(44)
$$_{L_1L_2}^\kappa (n_1l_1n_2l_2,n_1^{}l_1^{}n_2^{}l_2^{};\omega )=q^{2+\kappa }\upsilon (q;\omega )_{L_1}^0(n_1l_1n_2l_2;q)_{L_2}^0(n_1^{}l_1^{}n_2^{}l_2^{};q)𝑑q$$
(45)
where $`\upsilon (q;\omega )`$ are functions describing the energy of the intermediate neutrinos in various gauge models (light-neutrino mixing, etc.) . The parameter $`\omega `$ is related to the excitation energy of the intermediate nucleus. Using the compact formalism presented in the previous section the above integrals are simplified and written as follows:

a. For the radial moments $`_L^\kappa (n_1l_1n_2l_2;q)`$ working as in Eq. (39) we find

$`_L^\kappa (n_1l_1n_2l_2;q)=b^\kappa y^{L/2}e^y{\displaystyle \underset{\mu =0}{\overset{n_{max}}{}}}\epsilon _\mu ^Ly^\mu `$ (46)
which implies that the nuclear moments $`_L^\kappa (n_1l_1n_2l_2;q)`$ are simply obtained with the aid of the coefficients $`\epsilon _\mu ^L(n_1l_1n_2l_2)`$ of Eq. (40). Note that for $`\kappa =0`$ Eq. (46) reduces to Eq. (39).

b. By replacing the radial integrals in the right-hand side of Eq. (45) with Eq. (39) and manipulating their product we take

$`_{L_1}^0(n_1l_1n_2l_2;q)_{L_2}^0(n_1^{}l_1^{}n_2^{}l_2^{};q)=e^{2y}y^{(L_1+L_2)/2}{\displaystyle \underset{\mu =0}{\overset{n_{max}}{}}}c_\mu y^\mu `$ (47)
where

$$n_{max}=n_{1max}+n_{2max},$$

with

$$n_{1max}=n_1+n_2+(l_1+l_2L_1)/2,n_{2max}=n_1^{}+n_2^{}+(l_1^{}+l_2^{}L_2)/2$$

The coefficients $`c_\mu `$ in Eq. (47) are obtained from those of $`\epsilon _\nu ^L`$ as

$$c_\mu =\underset{\lambda =0}{\overset{\mu }{}}\epsilon _\nu ^{L_1}\epsilon _{\mu \nu }^{L_2}$$
(48)
In the last summation $`\epsilon _\nu ^{L_1}=\mathrm{\hspace{0.17em}0}`$, for $`n_{1max}<\nu n_{max}`$, and $`\epsilon _{\mu \nu }^{L_2}=\mathrm{\hspace{0.17em}0}`$, for $`n_{2max}<\mu \nu n_{max}`$. By inserting Eq. (47) into Eq. (45) and putting $`q=\mathrm{\hspace{0.17em}2}y^{1/2}/b`$ \[see Eq. (39)\], we have

$$_{L_1L_2}^\kappa (n_1l_1n_2l_2,n_1^{}l_1^{}n_2^{}l_2^{};\omega )=\frac{2^{\kappa +2}}{b^{\kappa +3}}\underset{\mu =0}{\overset{n_{max}}{}}c_\mu ^\lambda (\omega )$$
(49)
where

$$\lambda =\mu +(L_1+L_2+\kappa +1)/2$$

The quantities $`^\lambda (\omega )`$ represent the integrals

$$^\lambda (\omega )=_0^{\mathrm{}}\upsilon (y;\omega )e^{2y}y^\lambda 𝑑y.$$
(50)
which can be computed after providing the exact form of $`\upsilon (y;\omega )`$ . Equation (49) is much simpler than Eq. (D.1) of Ref. , since the latter does not include a simple summation as Eq. (49), but a double one over similar integrals to those of Eq. (50).

## 5 Summary and Conclusions

Our main concern in the present paper was the general formalism needed for calculating the weak and electromagnetic processes in nuclei. Starting from the decomposition of the hadronic current into tensor multipole operators, we investigated the possibility of constructing analytic formulas for the matrix elements of the principal semi-leptonic operators. We payed special attention on the multipole matrix elements taken between single-particle states and restricted ourselves on the formulation of their radial part.

By utilizing harmonic oscillator basis, we have achieved all types of radial moments related to the seven basic one-body operators (involved in semi-leptonic reactions), to be written in the form of closed analytic expressions, i.e. as products of an exponential times a polynomial of even powers in the momentum transfer q. The coefficients of the polynomials found are in general simple numbers (for diagonal matrix elements they are always rational numbers) and they can be calculated with simple codes. We applied the above compact expressions in order to simplify the previous formalism  giving the nuclear moments in the interesting neutrinoless double-beta decay process. This helped us to see the advantages of the formalism described in Secs. 3 and 4.

## 6 Acknowledgements

This work has been supported in part by the Research Committee of the University of Ioannina.

Appendix A

1. The identities which are needed in the multipole decomposition procedure are

$$exp(i𝐪.𝐫)=\mathrm{\hspace{0.17em}4}\pi \underset{JM_J}{}(i)^JM_M^J(𝐫)Y_{JM}^{}(\widehat{e}_{q_0})$$
(51)
$$\widehat{e}_{q_\lambda }exp(i𝐪.𝐫)(2\pi )^{1/2}\underset{J1}{}(i)^J(2J+1)^{1/2}[\lambda 𝐌_M^{JJ}(𝐫)+\frac{1}{q}\times 𝐌_M^{JJ}(𝐫)]$$
(52)
where $`\lambda =\pm `$ and the unit vector $`\widehat{e}_{q_\lambda }`$ has the spherical components

$$\widehat{e}_{q_0}\widehat{e}_z=𝐪/q$$
(53)
$$\widehat{e}_{q_\pm }\frac{1}{\sqrt{2}}(\widehat{e}_x\pm i\widehat{e}_y)$$
(54)
The components of the rank-one spherical tensor $`_\mu `$ are

$$_0=\frac{}{z}$$
(55)
$$_\pm =\left(\frac{}{x}\pm i\frac{}{y}\right)$$
(56)
2. For a single free nucleon (using Lorentz covariance, conservation of parity, time-reversal invariance, and isospin invariance) the matrix elements of the vector and axial-vector currents, in spin-isospin space, are written as

$`𝐤^{}\lambda ^{};1/2m_t^{}|J_\mu (0)_{TM_T}|𝐤\lambda ;1/2m_t`$ $`=`$ $`i\overline{u}(𝐤^{}\lambda ^{})\left[F_1^{(T)}\gamma _\mu +F_2^{(T)}\sigma _{\mu \nu }q_\nu +iF_S^{(T)}q_\mu \right]u(𝐤\lambda )`$ (57)
$`\times `$ $`1/2m_t^{}|I_T^{M_T}|1/2m_t`$
$`𝐤^{}\lambda ^{};1/2m_t^{}|J_\mu ^5(0)_{TM_T}|𝐤\lambda ;1/2m_t`$ $`=`$ $`i\overline{u}(𝐤^{}\lambda ^{})\left[F_A^{(T)}\gamma _5\gamma _\mu iF_P^{(T)}\gamma _5q_\mu F_T^{(T)}\gamma _5\sigma _{\mu \nu }q_\nu \right]u(𝐤\lambda )`$ (58)
$`\times `$ $`1/2m_t^{}|I_T^{M_T}|1/2m_t`$
The plane-wave single-nucleon states are labelled with 3-momenta $`𝐤(𝐤^{})`$, helicities $`\lambda (\lambda ^{})`$ and isospins $`1/2m_t(1/2m_t^{})`$. The single-nucleon form factors $`F_X^{(T)}=F_X^{(T)}(q_\mu ^2)`$, with $`T=0,1`$, and $`X=1,2,S,A,P,T`$ (vector (Dirac), vector (Pauli), scalar, axial, pseudoscalar,and tensor) are all functions of the momentum transfer $`q_\mu ^2`$.

Appendix B

The radial part of the wave functions in a (three-dimensional isotropic) harmonic oscillator potential is written as

$`R_{nl}=N_{nl}x^le^{x^2/2}_n^{l+1/2}(x^2),`$ (59)
where $`x=r/b`$, with $`b`$ is the harmonic oscillator size parameter and $`N_{nl}`$ the normalization factor

$`N_{nl}^2={\displaystyle \frac{2n!}{b^3\mathrm{\Gamma }(n+l+3/2)}}.`$ (60)
$`\mathrm{\Gamma }(x)`$ denotes the known gamma function and $`_n^{l+1/2}`$ represent the Laguerre polynomials defined by:

$`_n^{l+1/2}(x)={\displaystyle \underset{m=0}{\overset{n}{}}}\mathrm{\Lambda }_m(nl)x^m={\displaystyle \underset{m=0}{\overset{n}{}}}{\displaystyle \frac{()^m}{m!}}\left(\begin{array}{c}n+l+1/2\\ nm\end{array}\right)x^m.`$ (61)
In Sect. 3, by writing $`|n(l1/2)jm_j`$, we mean $`R_{nlj}(r)[Y_l(\mathrm{\Omega }_r)\chi _{1/2}]_{m_j}^j`$ with $`\chi _{1/2}`$ the Pauli spinor. We note that for the harmonic oscillator it holds $`R_{nlj}(r)=R_{nl}(r)`$. The adopted, in the present work, sequence of single-particle states is: 0s1/2, 0p3/2, 0p1/2, …