text:
# On the filling factor of emitting material in the upper atmosphere of ϵ Eri (K2 V)

## 1 Introduction

It is well known from direct observations that the solar chromosphere and transition region are far from homogeneous (see, for example, the uv images shown by Reeves 1976 and Gallagher et al. 1998). In spatially resolved observations of the Sun, transition region line emission tends to be dominated by contributions from the network boundaries. The solar-stellar analogy tells us that a similar effect must be expected on the surfaces of other main-sequence stars. Therefore, it is of interest to determine the filling factor of emitting material in solar-like stars and look for correlations with stellar parameters, such as has been done for surface (photospheric) magnetic filling factors by Montesinos & Jordan (1993). Such correlations may help to identify the physical processes which control coronal/chromospheric activity.

The current generation of high resolution spectrographic instruments are providing spectra of unprecedented quality in the ultraviolet (uv) and X-ray wavelength ranges. Unfortunately, such data cannot provide direct information about the spatial distribution of emitting material in main-sequence star atmospheres since it is beyond the scope of current instruments to resolve dwarf star surfaces, with the exception of the Sun. There is little prospect of overcoming this limitation in the near future.

In the absence of spatial information, semi-empirical modelling of main-sequence star spectra is effectively limited to the creation of “mean” 1-D atmospheric models whose parameters are those of a hypothetical homogeneous atmosphere which reproduces the observed spectrum (e.g. the $`ϵ`$ Eri models constructed by Kelch 1978 or Sim 2002). Although such models are useful, tighter constraints on the true atmospheric parameters, and therefore the nature of the physical processes which operate, could be achieved if the distribution of emitting material were known. (Information on the structure within the atmospheres of cool giant stars has been obtained through studies of molecular fluorescence, see McMurry & Jordan 2000.) The approach which could be taken is illustrated by the work of Cuntz et al. (1999). They have explored two-component chromospheres of K2 V stars, in which one component is heated by longitudinal tube waves and the other by acoustic waves. They simulate the supergranulation network using a flux tube model based on empirical correlations between observed surface magnetic fluxes and stellar rotation rates (Cuntz, Ulmschneider & Musielak, 1998) and calculate the flux in the Ca ii H and K lines expected from magnetoacoustic and acoustic waves. In the absence of observational constraints they use an average chromospheric network filling factor of 0.305. Thus a method which gives the filling factor in the transition region between the chromosphere and corona would provide valuable constraints on the above type of theoretical modelling.

Kopp (1972) discussed the energy balance in the solar transition region, in the light of early observations of the solar supergranulation structure. He introduced an area filling factor and by comparing emission measure distributions based on spatially integrated observations with the predictions of the energy balance equation he proposed that the transition region emission is restricted to an area of about 20 per cent of the solar surface area. In this case, the conductive flux in the mid-transition region is reduced and the lower transition region requires an additional source of heating. This type of structure was also explored by Gabriel (1976). The influence of the cross-section area and magnetic field expansion factor has also been recognized in the context of loop models of solar and stellar active regions (see Underwood, Antiochos & Vesecky 1981).

A method of deriving area factors has also been discussed by Jordan (2000), in which the total pressure is not kept constant but is allowed to vary according to hydrostatic equilibrium. Although Philippides (1996) used this approach to estimate the filling factor at $`T_e10^5`$ K in several main-sequence stars, the observational data available were not sufficient to find the variation of the filling factor with $`T_e`$. Here we apply the method to spatially unresolved stellar observations to derive, for the first time, the area filling factor as a function of $`T_e`$ in the upper transition region of $`ϵ`$ Eri.

The star selected for this analysis, $`ϵ`$ Eri (K2V), is a nearby dwarf that has moderately high levels of chromospheric/coronal activity and for which high quality spectroscopic data are available. This paper forms part of our study of the outer atmosphere of $`ϵ`$ Eri: see Jordan et al. (2001a), Jordan et al. (2001b), Sim & Jordan (2003) and Sim (2003, in preparation) for details.

In the next section the theory required for the analysis is discussed. Section 3 discusses the observations that are used to constrain the emission measure distribution (Section 4). Section 5 presents “first-cut” area filling factors deduced from the data and discusses possible sources of error. In Section 6, a self-consistent determination of the filling factors is presented. The results are discussed in Section 7 and conclusions are drawn in Section 8.

## 2 Theory

Spatially unresolved stellar observations can be used to constrain the apparent mean emission measure (see Pan & Jordan 1995) distribution defined by

$$Em^{0.3}(T_e)=_{T_e/\sqrt{2}}^{\sqrt{2}T_e}N_e(T)N_H(T)\frac{\text{d}r}{\text{d}T}G(r)f(r)\frac{A(r)}{A_{}(r)}\text{d}T$$
(1)
where $`T_e`$ is the electron temperature at radius $`r`$, $`N_e`$ and $`N_H`$ are the electron and hydrogen number densities respectively, $`f(r)=r^2/R_{}^2`$ ($`R_{}`$ is the photospheric stellar radius and $`f(r)=1`$ at the stellar surface). The dilution factor $`G(r)`$ specifies the fraction of photons (emitted in an optically thin line at radial coordinate $`r`$) that escape the star uninterrupted by the stellar surface. For spherical geometry,

$$G(r)=0.5\left(1\sqrt{1[f(r)]^1}\right).$$
(2)
The $`0.3`$ included in the notation for the emission measure in equation (1) is to remind the reader that the emission measure adopted is an integral over a logarithmic temperature range $`\mathrm{\Delta }\mathrm{log}T_e=0.3`$ dex. This temperature range is chosen since it is the typical temperature range across which individual spectral lines form in the transition regions of cool stars. The area filling factor $`A(r)/A_{}(r)`$ specifies the fraction of the total stellar surface area at radius $`r`$ ($`A_{}(r)=4\pi r^2`$) which is occupied by emitting material. The purpose of this paper is to show how values for $`A(r)/A_{}(r)`$ can be extracted from spatially unresolved data.

In the transition region where the pressure varies slowly as a function of the temperature, equation (1) can be written approximately as

$$Em^{0.3}(T_e)=G(r)f(r)\frac{A(r)}{A_{}(r)}\frac{\text{d}r}{\text{d}T_e}\frac{P_eP_H}{\sqrt{2}k_B^2T_e}$$
(3)
where $`P_e=N_ek_BT_e`$ and $`P_H=N_Hk_BT_e`$ are the electron and hydrogen pressures respectively, and are approximately constant over the region where an individual line is formed. The approximation that $`\text{d}r/\text{d}T_e`$ is constant in the region of line formation is less good, but is improved upon when finding final models by iteration.

As discussed by Jordan (2000), it is possible to explain the complete radiation losses from the upper transition region in terms of the divergence of the net conductive energy flux from the overlying corona. The justification for this is that the emission lines are formed over a pressure-squared scale height $`40,000`$ km, whereas the main heating of the quiet corona occurs at much greater heights. It cannot be ruled out that there is an additional source of heating in this part of the atmosphere, but it is certain that any additional heating must be small. The conductive flux can be expressed as

$$F_c=\kappa T_e^{5/2}\frac{\text{d}T_e}{\text{d}r}$$
(4)
where $`\kappa 1.1\times 10^6`$ ergs cm-2 s-1 K-7/2 at $`T_e10^6`$ K (Spitzer 1956). $`\kappa `$ increases at lower temperature (by about $`30`$ per cent by $`2\times 10^5`$ K), but it is taken as constant here for simplicity. Building on the simpler plane parallel approximation used by Jordan & Brown (1981), the logarithm of equation (3) can be differentiated and the temperature gradient replaced by the conductive flux using equation (4). By assuming an energy balance between conduction and radiation, the net conductive flux can be replaced by the radiation losses and the following equation obtained:

$$\begin{array}{c}\frac{\text{d}}{\text{d}\mathrm{log}T_e}\mathrm{log}\left[\frac{A_{}(r)Em^{0.3}}{A(r)}\right]=\frac{\text{d}\mathrm{log}G(r)}{\text{d}\mathrm{log}T_e}+2\frac{\text{d}\mathrm{log}f(r)}{\text{d}\mathrm{log}T_e}\hfill \\ \hfill +\frac{\text{d}\mathrm{log}A(r)/A_{}(r)}{\text{d}\mathrm{log}T_e}+2\frac{\text{d}\mathrm{log}P_g}{\text{d}\mathrm{log}T_e}+\frac{3}{2}\\ \hfill \frac{2P_{\text{rad}}(T_e)}{\kappa P_eP_HT_e^{3/2}G(r)^2f(r)^2}\left[\frac{A_{}(r)Em^{0.3}}{A(r)}\right]^2.\end{array}$$
(5)
$`P_{\text{rad}}(T_e)`$ is the radiative power loss function which describes the emissivity of a plasma of fixed composition at temperature $`T_e`$. $`P_g=N_gk_BT_e`$ is the total gas pressure.

This equation can be solved numerically to obtain the distribution of $`A_{}(r)Em^{0.3}/A(r)`$ as a function of $`T_e`$. By comparison with values of $`Em^{0.3}(T_e)`$ (deduced from observations) this allows the area filling factor $`A(r)/A_{}(r)`$ to be derived.

## 3 Observations

The mean emission measure distribution can be constrained using observed emission line fluxes. The details of this procedure are well known and have been discussed previously by various authors (e.g. Jordan & Brown 1981, Jordan et al. 1987; Griffiths & Jordan 1998).

In this paper, the mean emission measure distribution in the upper transition region of $`ϵ`$ Eri is constrained using the observed fluxes of emission lines of iron. For the most part, the observed fluxes used are from the analysis performed by Schmitt et al. (1996) using data recorded by the Extreme Ultraviolet Explorer (EUVE) satellite. The data and line flux measurements are discussed by Schmitt et al. (1996) and Laming, Drake & Widing (1996). Schmitt et al. (1996) give fluxes for Fe lines from ionisation stages ixxvi and xviiixxi. For the current work, only one line from each ionisation stage detected in the EUVE data is used. For Fe xii the three lines at around 193 Å have been added to give the total flux in the 4S –4P multiplet. In other cases where there are several lines from the same ion reported by Schmitt et al. (1996), the strongest line has generally been used. The lines used are listed in Table 1 together with the percentage errors in the measured fluxes given by Schmitt et al. (1996).

The observed EUVE fluxes have been converted to fluxes at the stellar surface using the stellar parameters for $`ϵ`$ Eri listed in table 1 of Jordan et al. (2001a). Subsequent to the work carried out by Schmitt et al. (1996), a more accurate determination of the hydrogen column density along the line of sight to $`ϵ`$ Eri has been made by analysis of interstellar absorption features in the Lyman $`\alpha `$ line profile (Dring et al. 1997). This gives a column density of $`\mathrm{log}N(\text{cm}^2)=17.875`$, which is smaller than that used by Schmitt et al. (1996) by close to a factor of 1.5. The wavelength dependent transmission of the interstellar medium with the new column density has been taken from calculations by Philippides (1996) and has been used to correct the observed fluxes for interstellar absorption.

The EUVE data are supplemented by our recent detection of the forbidden ultraviolet Fe xii lines in STIS spectra of $`ϵ`$ Eri. These data, and the measurement of the Fe xii line fluxes, have been presented elsewhere (Jordan et al. 2001a).

## 4 Emission measure loci

As discussed by Jordan & Brown (1981), each line flux can be used to construct a locus which places an approximate upper bound on the mean emission measure distribution:

$$Em^{0.3}(T_e)<F_{}/K(T_e)$$
(6)
where $`F_{}`$ is the stellar surface flux for the line under consideration and $`K(T_e)`$ is a contribution function defined by

$$K(T_e)=\frac{hc}{\lambda }\frac{N_u}{N_IN_e}\frac{N_I}{N_E}\frac{N_E}{N_H}A_{ul}$$
(7)
where $`\lambda `$ is the wavelength of the line, $`A_{ul}`$ is the Einstein coefficient for spontaneous emission in the line, and the $`N`$’s are number densities with the subscripts $`u`$, $`I`$ and $`E`$ referring to ions in the upper level in the transition, the ionisation stage in which the transition occurs and the element giving rise to the transition respectively. $`K(T_e)`$ has been calculated for each of the lines listed in Table 1, and for the Fe xii 1242-Å line. The CHIANTI atomic database (v4; Dere et al. 1997; Young et al. 2003) was used to compute $`N_u/N_I`$ as a function of $`T_e`$.The calculations were performed using the best value of the mean transition region electron pressure $`\mathrm{log}P_e(\text{cm}\text{-3}\text{ K})=15.68`$ obtained by Jordan et al. (2001b). Two sets of calculations for the ionisation fractions of iron have been used (Arnaud & Rothenflug 1985; Arnaud & Raymond 1992) and the stellar photospheric iron abundance \[Fe/H\]=-0.09 found by Drake & Smith (1993) has been adopted. (The solar iron abundance was taken from Grevesse & Sauval 1998). Fig. 1 shows the loci $`F_{}/K(T_e)`$ derived from the observed fluxes using the ionisation balance calculations of Arnaud & Rothenflug (1985) and Fig. 2 shows the loci obtained with the ionisation balance of Arnaud & Raymond (1992).

As expected, the emission measure loci suggest a smooth emission measure distribution which increases with temperature up to a peak at several million degrees. The temperature of the peak emission measure is difficult to determine since there are no Fe xvii lines in the data which would span the critical temperatures range. Comparing Figs. 1 and 2 shows that the choice of ionisation balance affects the form of the emission measure distribution. In particular there are significant differences in the ion populations of Fe ix, xvi and xviii which alter the apparent shape. The Arnaud & Raymond (1992) ionisation balance calculations suggest a higher peak temperature ($`\mathrm{log}T_e(\text{K})6.6`$) than those of Arnaud & Rothenflug (1985) ($`\mathrm{log}T_e(\text{K})6.4`$). In their analysis, Laming et al. (1996) adopted the Arnaud & Rothenflug (1985) ionisation balance (leading them to a coronal temperature $`\mathrm{log}T_e(\text{K})6.4`$) since solar data suggests that these calculations are more appropriate than those of Arnaud & Raymond (1992) for Fe ix (see the discussion by Laming, Drake & Widing 1995). In this paper, all calculations are carried out with both sets of ionisation fractions, and the differences between the results are assumed to be indicative of the errors in the ionisation balance.

Above the peak temperature the emission measure drops with increasing temperature. There appear to be some discrepancies between the loci for the high temperature ions. In particular the Fe xx locus appears to lie too high compared to those of Fe xix and Fe xxi. Material at temperatures which are significantly greater than that of the peak emission measure most likely occurs in stellar active regions, where it is magnetically confined in hot loops (such as discussed by Schmitt et al. 1996) and will not be discussed further here. Solar observations support this interpretation.

The locus from the 1242-Å STIS line lies above that from the EUVE Fe xii lines by about a factor of 1.5. This is significant but much smaller than the discrepancy reported by Jordan et al. (2001a) (they reported a factor of 3 difference between the STIS and EUVE results). Jordan et al. (2001a) computed the EUVE Fe xii locus using the flux reported by Laming et al. (1996) for the 195-Å line only. The EUVE Fe xii loci shown in Figs. 1 and 2 are computed using the fluxes for all three lines in the 4P – 4S multiplet from Schmitt et al. (1996). If the three fluxes reported by Schmitt et al. (1996) are used separately, the 195-Å line gives a locus that is significantly higher than that obtained by combining all three lines of the multiplet while the 192-Å line gives a locus that is significantly lower. The 193-Å line gives a locus that agrees well with that of the complete multiplet. It is plausible that the discrepancy (factor of 1.5) between the EUVE and STIS Fe xii loci is the result of systematic errors in the atomic data, but the measurement errors for the EUVE fluxes (see Table 1) are large enough to account for most of the difference. As suggested by Jordan et al. (2001a), there may be a real difference related to variations in coronal activity between the times at which the two sets of observations were made, but the uncertainty in the EUVE flux is too large to establish this with certainty.

## 5 Emission measure models

### 5.1 Models

As discussed in Section 2, equation (5) can be used to construct theoretical emission measure distributions under the assumption of an energy balance between thermal conduction and radiation in the upper transition region. To determine a unique emission measure distribution from equation (5) two boundary conditions are required: the pressure and emission measure must both be fixed at some temperature. The pressure at all other points is then determined by enforcing hydrostatic equilibrium. For this, a contribution to the pressure from turbulent motions is included, calculated using a most-probable turbulent velocity $`\xi 21`$ km s-1 which is appropriate for the temperature range under consideration (Sim & Jordan 2003). (Here, this turbulent pressure contribution is very small.) With the pressure so determined, equation (5) can be integrated numerically to determine the emission measure at all other points.

The first boundary condition adopted is that the electron pressure at $`\mathrm{log}T_e(\text{K})=5.3`$ should be consistent with the mean value determined by Jordan et al. (2001b). The effects of slightly varying the choice of fixed electron pressure at $`\mathrm{log}T_e(\text{K})=5.3`$ have been explored (see below).

It is a natural consequence of the signs of the terms in equation (5) that the computed emission measure distribution will pass through a minimum in the transition region when the magnitude of the last term in equation (5) is close to $`3/2`$. Therefore, the second boundary condition has been chosen to be that this minimum should occur at $`\mathrm{log}T_e(\text{K})=5.3`$, which is approximately the temperature at which the minimum in the emission measure distribution is observed to occur in the Sun (see e.g. Macpherson & Jordan 1999). The emission measure distribution for $`ϵ`$ Eri (Sim 2002) shows a decrease up to at least $`\mathrm{log}T_e(\text{K})=5.3`$.

Using these two boundary conditions, emission measure distributions ($`A_{}(r)Em^{0.3}/A(r)`$) have been computed using equation (5). At this stage, the third term on the right-hand side (RHS) of equation (5) which accounts for the variation in the area factor with $`T_e`$ has been ignored. The effect of this term is discussed in Section 6. The distributions have been calculated from an upper boundary temperature $`\mathrm{log}T_{\text{max}}=6.5`$ down to a lower boundary temperature $`\mathrm{log}T_{\text{min}}=5.3`$. Distributions of $`A_{}(r)Em^{0.3}/A(r)`$ were calculated using a range of starting values for $`A_{}(r)Em^{0.3}/A(r)`$ and the electron pressure at the top point in the model in order to find the distributions which satisfies the chosen boundary conditions (see above). Two of the computed distributions are shown in Fig. 3 (the solid and dotted lines). These distributions differ only in the choice of fixed pressure at the bottom end. Fig. 4 shows the electron pressure ($`P_e=N_eT_e`$) as a function of $`T_e`$ for the same models. The first model (solid line in the figures) has a pressure which is slightly lower than, but close to, the best pressure ($`\mathrm{log}P_e(\text{cm}\text{-3}\text{ K})15.68`$) while the other (dotted line) has a higher pressure, just within the error margins of the best pressure ($`\mathrm{log}P_e(\text{cm}\text{-3}\text{ K})15.78`$).

In both models, it can be seen that above the mid-transition region the theoretical emission measure increases roughly proportional to $`T^{3/2}`$, owing to the dominance of the penultimate term in equation (5). At higher temperatures, the emission measure begins to flatten off and eventually starts to decrease. This is the result of the negative contributions to the RHS of equation (5) from the 1st and 4th terms which become significant at high temperatures (the last term is small at high temperatures).

### 5.2 Filling factors

The theoretical distributions plotted in Fig. 3 can be used to estimate fractional filling factors for each of the observed iron lines. The apparent emission measure distribution can be used to compute approximate line fluxes via

$$F_{}K(T_e)Em^{0.3}(T_e)\text{d}\mathrm{log}T_e/0.3$$
(8)
where the integral runs over the whole temperature range of the model. If the apparent emission measure in equation (8) is replaced with the computed distribution of $`A_{}(r)Em^{0.3}/A(r)`$ then $`F_{}A_{}(r)/A(r)`$ is obtained instead of $`F_{}`$. This can be compared to the observed value of the line flux to deduce $`A(r)/A_{}(r)`$ in the region of line formation. The value so deduced is only an estimate because of the approximate nature of equation (8) and must be interpreted as a weighted average of the area filling factor over the region of line formation. It should be a good estimate provided that the physical extent of the region of line formation is reasonably small (as is the case for all the lines considered here).

The filling factors deduced from the distributions discussed above are given in Table 2. Filling factors are only given for individual lines that form in the upper transition region (Fe ix – Fe xvi) since the higher temperature lines probably form under different conditions in hot active region loops (see Section 4). For each line the approximate temperature of line formation (determined from the emission measure loci) and four filling factors are given. These correspond to the four possible combinations of the two models and the two sets of ionisation balance calculations (Arnaud & Rothenflug 1985; Arnaud & Raymond 1992). In addition to the filling factors deduced from the upper transition region lines, filling factors are given in Table 2 for the mid-transition region ($`\mathrm{log}T_e(\text{K})5.3`$) which were deduced by comparing the computed $`A_{}(r)Em^{0.3}/A(r)`$ with values of $`Em^{0.3}`$ from a mean transition region emission measure distribution constructed from STIS data (Sim 2002). This distribution will be discussed in a forthcoming paper (Sim 2003, in preparation).

The models suggest that the mid transition region filling factor is several tens per cent and that the filling factor starts increasing with temperature above $`\mathrm{log}T_e(\text{K})5.9`$, becoming close to 1.0 at coronal temperatures. To illustrate this point, Fig. 5 shows (open symbols) a plot of the area filling factors deduced from the model with higher pressure using the Arnaud & Raymond (1992) ionisation balance (the error bars indicate the random measurement errors but do not attempt to address the potential systematic errors which are discussed below). This trend with temperature is striking and consistent with what would be expected from the standard model of expanding magnetic funnels in the upper atmosphere (Gabriel 1976).

The filling factors are not strongly sensitive to the adopted boundary pressure (they are 20 – 40 per cent smaller at the larger pressure). As discussed above, however, several of the filling factors are quite sensitive to the adopted ionisation balance.

On physical grounds, the filling factor cannot be greater than 1.0. In all cases the filling factor calculated for Fe xvi is greater than 1.0 and, for all cases except the high pressure/Arnaud & Raymond (1992) combination, the filling factor is unphysical for some of the other lines. In most cases the random measurement errors are not sufficient to explain these unphysical results, and they point towards a systematic error in the calculations. Several possible sources of systematic errors are discussed below.

### 5.3 Systematic errors

#### 5.3.1 Ionisation balance

The quality of the atomic data is probably the most significant source of uncertainty, particularly those used in the ionisation balance. By comparing the filling factors deduced using Arnaud & Rothenflug (1985) and Arnaud & Raymond (1992) it is clear that for several of the ions (Fe ix, Fe xv and Fe xvi) the uncertainty in the ion balance may be comparable to, or even greater than, a factor of 2. For Fe xii and Fe xiii the uncertainty is smaller ($`40`$ per cent) and it is negligible for Fe x and Fe xi.

#### 5.3.2 Hydrogen column density

The results are sensitive to the adopted interstellar hydrogen column density since it determines the wavelength dependent transmission factor of the interstellar medium. This is most important for Fe xvi (the longest wavelength line). The error on the column density quoted by Dring et al. (1997) is only 20 per cent, however, and so this should not be a dominant contribution to the uncertainty.

#### 5.3.3 The abundance of iron

A potentially important source of systematic error is the adopted iron abundance. The quoted error in the photospheric abundance derived by Drake & Smith (1993) is only 0.05 dex, but it is possible that the iron abundance in the upper transition region is different from that in the photosphere. It has been suggested that the abundances of elements with low ($`<10`$ eV) first ionisation potentials (FIPs) are enhanced in the outer parts of cool star atmospheres relative to elements with high FIPs (Crawford, Price & Sullivan 1972 first identified this effect in the solar wind). This is known as the FIP-effect (see Jordan et al. 1998 and references therein for a full discussion). Laming, Drake & Widing (1996) were unable to conclusively detect an enhancement of the iron abundance in the corona of $`ϵ`$ Eri (as would be expected if the FIP-effect were acting), however they could not rule out an enhancement of similar magnitude to that in the Sun (increase by a factor of 2 – 3). If the upper transition region/coronal abundance is enhanced relative to the photospheric value adopted in the models, then the predicted area filling factors in Table 2 will all be systematically overestimated. (It can be shown from equation (5) that an increase in the adopted iron abundance by a factor $`x`$ would lead to a reduction in the derived area factors by approximately a factor of $`\sqrt{x}`$.)

#### 5.3.4 The stellar surface gravity

Errors in the assumed value for the stellar surface gravity $`g`$ will affect the calculated area filling factors in a systematic way. In particular, if the adopted gravity is too large, the area factors found for the high temperature lines will be too large. This is because the magnitude of the 4th term on the RHS of equation (5) is proportional to the adopted gravity. Since that term is negative, if the adopted gravity is too large it will mean that the peak in the calculated emission measure distribution will occur at too low a temperature (recall that the 4th term is important in determining where the peak occurs). To investigate this effect models have been computed using a lower surface gravity, $`\mathrm{log}g=4.65`$, which is consistent with the error bar on the value of the surface gravity found by Drake & Smith (1993) ($`\mathrm{log}g=4.75\pm 0.1`$). The emission measure distributions obtained with the lower gravity (using the same lower boundary pressures as before) are plotted in Fig. 6. These emission measure distributions have been used to compute area filling factors (using the Arnaud & Raymond 1992 ionisation balance) for comparison with those found in Table 2. The area factors deduced for the EUVE lines are given in Table 3.

Comparing the area factors in Table 3 with the appropriate entries in Table 2 shows that, as expected, the high temperature lines have the greatest sensitivity to the adopted gravity: for Fe ix – Fe xii the difference in the adopted gravity decreases the computed area factor by less than 10 per cent, but for Fe xvi the reduction is close to 50 per cent.

Therefore, it is possible that the unphysical area factors for the high temperature Fe xvi line could arise because the adopted surface gravity is too large. (Note that earlier work, e.g. by Kelch 1978 did adopt the significantly lower surface gravity of $`\mathrm{log}g=4.5`$.)

#### 5.3.5 The radiative power loss function, $`P_{\text{rad}}(T_e)`$

In all the calculations discussed here, the fit to the radiation loss function

$$P_{\text{rad}}(T_e)=1.25\times 10^{16}T_e^1\text{ergs cm}\text{3}\text{ s}\text{-1}$$
(9)
derived by Philippides (1996) has been adopted. Across the temperature range under consideration, this fit is accurate to within a factor of $`2`$. For future work, in which it is hoped that the other sources of systematic error may be eliminated, it would be valuable to obtain a more accurate power loss function using modern atomic data. For the present, however, this simple form is sufficiently accurate given the uncertainties in the ionisation fractions and iron abundance, and that the last term in equation (5) is only important at the low temperature end of the models.

#### 5.3.6 Discussion of systematic errors

To conclude this section, it seems probable that the uncertainties in the deduced filling factors are dominated by those in the adopted ionisation balance calculations. The uncertainty in the Fe xvi filling factor is likely to be the largest (a multiplicative factor of 2 – 3) since it is the most sensitive to errors in the adopted hydrogen column density, surface gravity and has amongst the most uncertain ionisation fraction. These sources of error are certainly large enough to explain the moderately unphysical filling factors deduced for this line.

The errors are likely to be smallest in the Fe x, Fe xi, Fe xii and Fe xiii filling factors since these are the least sensitive to the systematic errors discussed above.

The scale of the uncertainties in all cases is small enough that the filling factors presented in Table 2 can be regarded as reasonable estimates. That the higher pressure model with the Arnaud & Raymond (1992) ionisation balance gives the “most physical” set of filling factors could be interpreted as meaning that this is the best combination of pressure and ionisation balance, but the errors are large enough that this cannot be clearly established.

The detected trend with temperature is not expected to be the result of any of the errors discussed here and is almost certainly real: $`A(r)/A_{}(r)`$ varies systematically by more than a factor of 5 over the temperature range considered. This variation is larger than that expected for any of the uncertainties discussed above.

## 6 Self-consistent calculation

The estimates for $`A(r)/A_{}(r)`$ presented above are not self-consistent. This is because when equation (5) was used the 3rd term on the RHS was neglected, but the implied filling factors suggest a significant variation of $`A(r)/A_{}(r)`$ with $`\mathrm{log}T_e`$. To investigate this inconsistency, two more models have been constructed accounting for the previously neglected term.

To do this, the area filling factor must be specified as a function of temperature in the modelling. In principle self-consistency between the variation of the area factor adopted in the modelling and the derived area factors could be achieved iteratively, but since there are significant uncertainties (see above) a rigorous iterative solution is not warranted and only approximate self-consistency is pursued. To this end a simple form for $`A(r)/A_{}(r)`$ as a function of $`T_e`$ has been assumed in the calculation of the two new models, based on the approximate area factors presented in Table 2. Given the form of the relevant term in equation (5), it is computationally convenient to express the relationship as a power law. The relationship adopted in the calculations is

$$\mathrm{log}\frac{A(r)}{A_{}(r)}=\{\begin{array}{cc}0.7\hfill & \mathrm{log}T_e<5.8\hfill \\ 6.5+\mathrm{log}T_e\hfill & 5.8<\mathrm{log}T_e<6.5\hfill \end{array}.$$
(10)
This relationship crudely describes the form suggested by the results in the last column of Table 2: constant below $`\mathrm{log}T_e=5.8`$ and then increasing to approximately 1.0 in the corona. Distributions of $`A_{}(r)Em^{0.3}/A(r)`$ calculated from the two new models (which, like the first two models differ only in the adopted lower boundary pressure) are plotted in Fig. 3 (the dashed and dashed-dot lines). The electron pressures in these models are plotted in Fig. 4.

Comparing the distributions in Fig. 3 shows the effect of the area variation term: it steepens the predicted emission measure gradient before the peak (the kink in the dashed and dashed-dot lines is not physical but due to the discontinuity in equation 10). Also, above the peak the emission measure declines more rapidly. This is because the larger emission measures around the peak lead to a more rapid decline in the pressure with temperature, making the magnitude of the 4th term on the RHS of equation (5) greater.

Table 4 gives the area factors deduced from the new models in the same format as those given in Table 2. The area factors deduced from the new model with the higher pressure and the Arnaud & Raymond (1992) ionisation balance are shown in Fig. 5 (filled symbols). A curve showing the assumed $`A(r)/A_{}(r)`$ (equation 10) is over-plotted for comparison.

Given that there are several sources of systematic error, the agreement between the form of $`A(r)/A_{}(r)`$ used in the modelling and the values of $`A(r)/A_{}(r)`$ predicted by the model is sufficiently good that the higher pressure model can be regarded as self-consistent if the Arnaud & Raymond (1992) ionisation balance is adopted. It is noteworthy that the largest discrepancies are for Fe ix, Fe xv and, in particular, Fe xvi which are the ions that show the largest uncertainties in the ionisation balance. The area factors deduced from the model with lower pressure are slightly less self-consistent than those obtained from the higher pressure model, being systematically larger than those assumed in the modelling. Given the scale of the other sources of error, however, the agreement is adequate.

Comparing Tables 2 and 4, it is clear that the main effect of including the term describing the variation of the area factor with temperature is to reduce the implied area factors for most of the high temperature lines, leading to fewer unphysical ($`A(r)/A_{}(r)>1.0`$) values. This is the result of the steeper increase of emission measure with temperature. The only exception is Fe xvi for which the implied area factor is larger when the variation of the area factor is included in the model if the Arnaud & Raymond (1992) ionisation fractions are used. This is because Fe xvi forms at temperatures above the peak of the emission measure distribution, if the Arnaud & Raymond (1992) ionisation balance calculations are adopted.

In summary, the self-consistent treatment of the variation of the area factor presented in this section indicates that neglecting this variation in the modelling may cause the area factors in the upper transition region to be systematically overestimated to a small, but significant, extent. Accounting for the variation of the area factors does not, however, alter the apparent trend that the area factor increases across the upper transition region. For now it is concluded that the area factors presented in Table 4 (or equivalently described by equation 10) are the best current estimates of the area filling factors in the upper transition region of $`ϵ`$ Eri and that most of them are accurate to within a factor of two.

## 7 Discussion and prospects for future development

In Sections 5 and 6 it has been shown that the theory described in Section 2 can be used to provide estimates of the area filling factors of emitting material at various temperatures in stellar transition regions.

The results show (Fig. 5) that in $`ϵ`$ Eri the area filling factor increases with temperature markedly in the upper transition region, from $`0.2`$ at $`\mathrm{log}T_e<5.8`$ to $`1.0`$ in the corona. There is no evidence for significant variations across the mid-transition region $`5.3<\mathrm{log}T_e<5.8`$.

The variation of the area factor derived for $`ϵ`$ Eri is similar to that found from observations of the solar supergranulation network made with the Solar and Heliospheric Observatory (SOHO) (Gallagher et al. 1998). They find (see their figs. 5 and 8, in particular) that at the resolution of the Coronal Diagnostic Spectrometer (CDS) ($``$ 2 arcsec by 2 arcsec) the area occupied by the emitting material increases dramatically in the upper transition, but does not vary significantly across the mid-transition region. These results are in general agreement with earlier observations and the model by Gabriel (1976). When the method used here is applied to CDS observations of the average network, a filling factor of $``$ 0.2 at $`T_e2\times 10^5`$ K is deduced within the network (Smith & Jordan 2002), which typically occupies $``$ 50 per cent of the solar area at this temperature (Gallagher et al. 1998). The filling factors which we have derived are based on the assumption that the observed emission is dominated by the regions occupying the small areas. These could refer to structures within the supergranulation network of $`ϵ`$ Eri and may not be comparable to the total solar network area observed with the resolution of CDS. In this case, the area factors are larger than those found for the Sun.

Using the solar analogy, the area factor found at $`T_e2\times 10^5`$ K can be applied to the lower transition region. This increases the intrinsic radiation losses and since thermal conduction now plays no role below $`10^5`$ K, an additional source of heating is required. As found by Sim & Jordan (2003), if the observed non-thermal emission line widths are associated with an Alfvén wave (or fast-mode wave) flux, there is ample energy to heat both the corona and the lower transition region, although uncertainties concerning wave propagation and dissipation remain to be resolved.

The filling factor of $`0.2`$ found for the mid-transition region of epsilon Eri can be compared with the value of 0.3 assumed for the chromosphere by Cuntz et al. (1999). Since the filling factor of magnetic flux tubes should increase with height, our results suggest that a smaller chromospheric filling factor of $`0.2`$ would be more appropriate; this could be explored in future two-component models.

The accuracy of the results presented in this paper is limited by the various systematic sources of error which were discussed in Section 5.3. In particular the analysis of larger sets of lines, ideally covering a wide range of elements and ionisation species would help. This would reduce the systematic errors associated with uncertainties in the ionisation fractions for iron and possibly help to identify possible errors arising from the adopted elemental abundance. Suitable data for such analysis are gradually becoming available for a range of stars thanks to the current generation of X-ray satellite missions. X-ray data are also advantageous since their analysis is less sensitive to the assumed interstellar hydrogen column density. Improved calculations of atomic data for iron, and the radiation loss function $`P_{\text{rad}}`$ would also make the modelling more quantitatively reliable.

In future work this technique should be applied to a range of late-type stars with differing levels of activity. By doing so it should be possible to look for trends in the behaviour of the filling factors of emitting material with various stellar parameters. The ability to reproduce these trends will then be a new constraint on theoretical studies of the structure and heating of late-type stars, including the magnetic field geometry.

## 8 Conclusions

It has been shown that by assuming an energy balance between radiation and thermal conduction in the upper transition region, it is possible to use spatially unresolved stellar observations to estimate the fractional area filling factor of emitting material in stellar atmospheres.

The technique has been applied to $`ϵ`$ Eri and it has been found that the area filling factor is $`0.2`$ in the mid transition region and that it increases to $`1.0`$ in the corona. These filling factors are similar to those found from observations of the solar transition region, and from models in which the transition region consists of material in magnetic funnels whose radii are greatest is in the upper transition region/corona. Several potential sources of uncertainty in the calculations have been discussed and it is concluded that, although they are currently significant, they do not prevent a useful application of the methods adopted. Further observations and calculations of atomic data should be able to reduce the current uncertainties.

## Acknowledgments

This work was carried out while SAS was a research student in the Physics Department (Theoretical Physics) at the University of Oxford, supported by a PPARC studentship (PPA/S/S/1999/02862) and was completed while he was a PPARC funded PDRA in the Astrophysics Group at Imperial College London (PPA/G/S/2000/00032).

## Bibliography

Arnaud M., Raymond J., 1992, ApJ, 398, 394
Arnaud M., Rothenflug R., 1985, A&AS, 60, 425
Cuntz M., Ulmschneider P., Musielak Z. E., 1998, ApJ, 493,
L117
Cuntz M., Rammacher W., Ulmschneider P., Musielak Z.E.,
Saar S.H., 1999, ApJ, 522, 1053
Crawford H. J., Price P. B., Sullivan J. D., 1972, ApJL,
175, L149
Dere K. P., Landi E., Mason H. E., Monsignori Fossi B. C.,
Young P. R., 1997, A&AS, 125, 149
Drake J. J., Smith G., 1993, ApJ, 412, 797
Dring A. R., Linsky J., Murthy J., Henry R. C., Moos W.,
Vidal-Madjar A., Audouze J., Landsman W., 1997,
ApJ, 488, 760
Gabriel A. H., 1976, Phil. Trans. Roy. Soc. London, A291,
339
Gallagher P. T., Phillips K. J. H., Harra-Murnion L. K.,
Keenan F. P., 1998, A&A, 335, 733
Grevesse N., Sauval A. J., 1998, Space Science Reviews,
85, 161
Griffiths N. W., Jordan C., 1998, ApJ, 497, 883
Jordan C., 2000, Plasma Phys. Control. Fusion, 42, 415
Jordan C., Ayres T. R., Brown A., Linsky J. L., Simon T.,
1987, MNRAS, 225, 903
Jordan C., Brown A., 1981, in R. M. Bonnet & A. K.
Dupree eds., Solar Phenomena in Stars and Stellar
Systems, p. 199, Reidel. Dordrecht, Holland, NATO
ASIC, 68
Jordan C., Doschek G. A., Drake J. J., Galvin A. B.,
Raymond J. C., 1998, in ASP Conf. Ser. 154: Cool
Stars, Stellar Systems and the Sun, Vol. 10, p. 91
Jordan C., McMurry A. D., Sim S. A., Arulvel M., 2001a,
MNRAS, 322, L5
Jordan C., Sim S. A., McMurry A. D., Arulvel M., 2001b,
MNRAS, 326, 303
Kelch W. L., 1978, ApJ, 222, 931
Kopp R. A., 1972, Solar Physics, 27, 373
Laming J. M., Drake J. J., Widing K. G., 1995, ApJ, 443,
416
Laming J. M., Drake J. J., Widing K. G., 1996, ApJ, 462,
948
Macpherson K. P., Jordan C., 1999, MNRAS, 308, 510
McMurry A. D., Jordan C., 2000, MNRAS, 313, 423
Montesinos B., Jordan C., 1993, MNRAS, 264, 900
Pan H. C., Jordan C., 1995, MNRAS, 272, 11
Philippides D., 1996, D.Phil. Thesis, University of Oxford
Reeves E.M., 1976, Solar Physics., 46, 53
Schmitt J. H. M. M., Drake J. J., Stern R. A., Haisch B.
M., 1996, ApJ, 457, 882
Sim S. A., 2002, D.Phil. Thesis, University of Oxford
Sim S. A., Jordan C., 2003, MNRAS, 341, 517
Smith G. R., Jordan C., 2002, MNRAS, 337, 666
Spitzer L., 1956, Physics of ionized gases, Interscience, New
York
Underwood J. H., Antiochos S. K., Vesecky J. F., 1981,
in R. M. Bonnet & A. K. Dupree eds. Solar Phenomena
in Stars and Stellar Systems, p. 227, Reidel. Dordrecht,
Holland, NATO ASIC, 68
Young P. R., DelZanna G., Landi E., Dere K. P., Mason H.
E., Landini M., 2003, A&AS, 144, 135


format: text

text:
# Stabilizing the invisible axion in 3-3-1 models

## I Introduction

Recently, observations of the core mass distribution in the cluster of Galaxies Abell 2029 using the NASA’s Chandra X-ray Observatory suggest the existence of cold dark matter (CDM) . On the other hand, the Wilkinson Microwave Anisotropy Probe (WMAP) measurements of the cosmic microwave background temperature anisotropy and polarization are also consistent with CDM and a positive cosmological constant . Although the exact nature of the CDM is not known yet, candidates for this sort of matter are elementary particles such as neutralinos or the invisible axion . However, early invisible axion models  are unstable against quantum gravitational effects , which may generate a large axion mass and also spoil the value of the $`\overline{\theta }_{\mathrm{eff}}`$ parameter. One way to stabilize the axion is by considering large discrete gauge symmetries in the sense of Ref.  as was done in the multi-Higgs extension of the standard model , in the 3-3-1 model  or in the supersymmetric model . The search for dark matter is of course related to the search for new physics beyond the standard model which in turn is related to the existence of new fundamental energy scales. In the literature, the most easily recognized fundamental energy scales are those related to supersymmetry, the neutrino masses, grand unification, and superstring theory.

In this vein, it is worth recalling once more that it has been known for a long time that the measured value of the electroweak mixing angle $`\mathrm{sin}^2\theta _W(M_Z)=0.231131/4`$ appears to obey, at an energy scale $`\mu `$, an $`SU(3)`$ symmetry in such a way that $`\mathrm{sin}^2\theta _W(\mu )=1/4`$ . Hence, if the value of $`\mathrm{sin}^2\theta _W`$ is not an accident, it may be considered as an indication of a new fundamental energy scale of the order of a few TeVs. Notwithstanding, in models with $`SU(3)`$ electroweak symmetry there is trouble when we try to incorporate quarks. A solution to this issue is to introduce an extra $`U(1)`$ factor such as in 3-3-1 models , to embed the model in a Pati-Salam-like model , or even to embed it in theories of TeV gravity .

Independently of the axion or dark matter issues, 3-3-1 models are interesting possibilities, on their own, for physics at the TeV scale. At low energies they coincide with the standard model and some of them give at a least partial explanation of some fundamental questions that are accommodated but not explained by the standard model. For instance, i) in order to cancel the triangle anomalies the number of generations must be three or a multiple of three; ii) the model of Ref.  predicts that $`(g^{}/g)^2=\mathrm{sin}^2\theta _W/(14\mathrm{sin}^2\theta _W)`$; thus there is a Landau pole at the energy scale $`\mu `$ at which $`\mathrm{sin}^2\theta _W(\mu )=1/4`$, and according to recent calculations $`\mu 4`$ TeV ; iii) the quantization of the electric charge  and the vectorial character of the electromagnetic interactions  do not depend on the nature of the neutrinos, i.e., whether they are Dirac or Majorana particles; and iv) the model possesses $`𝒩=1`$ supersymmetry naturally at the $`\mu `$ scale . If right-handed neutrinos are considered to transform nontrivially, 3-3-1 models  can be embedded in a model with 3-4-1 gauge symmetry in which leptons transform as $`(\nu _l,l,\nu _l^c,l^c)_L(\mathrm{𝟏},\mathrm{𝟒},0)`$ under each gauge factor .

Models with $`SU(3)`$ (or $`SU(4)`$) symmetry may have doubly charged vector bosons. These types of bileptons may be detected by measuring the left-right asymmetries in Møller scattering , for instance, at the E158 SLAC experiments (which use 48 GeV polarized electrons scattering off unpolarized electrons in a liquid hydrogen target ); or in future lepton-lepton accelerators. It is interesting that the weak interaction’s parity nonconservation has never been observed in lepton-lepton scattering. Those asymmetries may also be used for seeking a heavy neutral $`Z^0`$ vector boson, which is also a prediction of these models, in $`e\mu `$ collisions . Singly and doubly charged vector bileptons may also be produced in $`e^{}\gamma `$  or $`\gamma \gamma `$  or hadron  colliders. New heavy quarks are also part of the electroweak quark multiplets in the minimal model representation. They are singlets under the standard model $`SU(2)_LU(1)_Y`$ group symmetry. In some versions their electric charge is different from the usual one, so that it can be used to distinguish such a model from their viable competitors. In fact, the $`p\overline{p}`$ production and decay of these exotic quarks at the energies of the Tevatron have been studied in Ref.  where a lower bound of 250 GeV on their masses was found. This sort of models is also predictive with respect to neutrino masses ; the models can implement the large mixing angle MSW solution to the solar neutrino issue , and also the almost bimaximal mixing matrix in the lepton sector .

Summarizing, from the present experimental data, say those from the CERN $`e^+e^{}`$ collider LEP, 3-3-1 models are safe if the symmetry breaking from 3-3-1 to 3-2-1 occurs at the level of TeVs; however, they have rich phenomenological consequences as we mentioned above. It will be interesting to search for some of the new particles that are present in these models, as extra Higgs scalars, exotic quarks and vector bileptons, at the energies of the upgrade DESY $`ep`$ collider HERA and Tevatron . The scalar sectors are equivalent to multi-Higgs-boson extensions of the standard model; for instance, under $`SU(2)_LU(1)_Y`$ the model with three triplets has two doublets and several non-Hermitian singlets, while the model with a sextet has three doublets, a complex triplet, and several complex singlets. In particular the neutral singlet ($`\chi ^0`$) is $`Z`$-phobic (its coupling with $`Z^0`$ vanishes when the scale of the $`SU(3)_L`$ symmetry goes to infinity) and for this reason it evades the LEP constraints. For a finite $`SU(3)_L`$ energy scale there are corrections that can be calculated by using the oblique S, T, and U radiative parameters which constrain the allowed masses for the leptoquarks and bileptons . These masses are of the same order of magnitude, a few TeV, as those allowed by the running of the coupling constants. Through the condition $`\mathrm{sin}^2\theta _W(\mu )=0.25`$, the running is sensitive to a new degree of freedom. Hence, the masses of exotic scalars and bileptons run from hundreds of GeV to a few TeV . We will return to this point later.

Turning back to the axion, the interesting point is that a Peccei-Quinn (PQ) symmetry  is almost automatic in the classical Lagrangian of 3-3-1 models. It is only necessary to avoid a trilinear term in the scalar potential by introducing a $`Z_2`$ symmetry . Unfortunately, even in this case the PQ symmetry is explicitly broken by gravity effects. In order to stabilize the axion, and at the same time automatically implement the PQ symmetry, we must introduce local discrete symmetries, $`Z_N`$. In fact, recently it was shown that in a version of the Tonasse and Pleitez 3-3-1 model  it is possible to implement both symmetries $`Z_{13}`$ and PQ automatically, thus the axion is naturally light and there is no domain wall problem .

We will consider in this work two 3-3-1 models in which only the known leptons transform nontrivially under the gauge symmetry, as in Refs. , but we add also right-handed neutrinos and exotic charged leptons transforming as singlets. In one model (model A) we consider a scalar sextet but it is possible to use only three scalar triplets (model B). Both models admit a large enough discrete $`Z_N`$ symmetry, implying a natural light invisible axion.

## II The axion in two 3-3-1 models

We will consider two versions of the 3-3-1 model of Ref. . In model A we use three scalar triplets and a sextet, while in model B we avoid the scalar sextet. In both models we introduce also a scalar singlet, $`\varphi (\mathrm{𝟏},\mathrm{𝟏},0)`$, and lepton singlets.

The quark and lepton sectors have the same representation content in both models. We have quarks transforming, under $`SU(3)_CSU(3)_LU(1)_N`$, as follows: $`Q_{mL}=(d_m,u_m,j_m)_L^T(\mathrm{𝟑},\mathrm{𝟑}^{},1/3),m=1,2`$; $`Q_{3L}=(u_3,d_3,J)_L^T(\mathrm{𝟑},\mathrm{𝟑},2/3)`$, and the corresponding right-handed components in singlets, $`u_{\alpha R}(\mathrm{𝟑},\mathrm{𝟏},2/3)`$, $`d_{\alpha R}(\mathrm{𝟑},\mathrm{𝟏},1/3)`$, $`\alpha =1,2,3`$; $`J_R(\mathrm{𝟑},\mathrm{𝟏},5/3)`$; $`j_{mR}(\mathrm{𝟑},\mathrm{𝟏},4/3)`$; the leptons are the known ones and transform as triplets $`(\mathrm{𝟑}_a,0)`$, $`\mathrm{\Psi }_{aL}=(\nu _a,l_a,l_a^c)_L^T`$; $`a=e,\mu ,\tau `$, and we also add right-handed neutrinos and a charged lepton in the singlets $`\nu _{aR}(\mathrm{𝟏},\mathrm{𝟏},0)`$, $`E_{L,R}(\mathrm{𝟏},\mathrm{𝟏},1)`$ . The scalar sector, in the minimal version, has only three triplets $`\chi =(\chi ^{},\chi ^{},\chi ^0)^T`$, $`\rho =(\rho ^+,\rho ^0,\rho ^{++})^T`$, $`\eta =(\eta ^0,\eta _1^{},\eta _2^+)^T`$, transforming as $`(\mathrm{𝟏},\mathrm{𝟑},1),(\mathrm{𝟏},\mathrm{𝟑},1)`$ and $`(\mathrm{𝟏},\mathrm{𝟑},0)`$, respectively, and a scalar singlet $`\varphi (\mathrm{𝟏},\mathrm{𝟏},0)`$.

With the quark and scalar multiplets above we have the Yukawa interactions

$`_Y^q`$ $`=`$ $`\overline{Q}_{iL}(F_{i\alpha }u_{\alpha R}\rho ^{}+\stackrel{~}{F}_{i\alpha }d_{\alpha R}\eta ^{})+\lambda _{im}\overline{Q}_{iL}j_{mR}\chi ^{}`$ (1)
$`+`$ $`\overline{Q}_{3L}(G_{1\alpha }u_{\alpha R}\eta +\stackrel{~}{G}_{1\alpha }d_{\alpha R}\rho )+\lambda _1\overline{Q}_{3L}J_{1R}^{}\chi `$ (2)
$`+`$ $`H.c.,`$ (3)
where repeated indices mean summation.

### A Model with a scalar sextet (Model A)

In this model we add a scalar sextet $`S(\mathrm{𝟏},\mathrm{𝟔},0)`$ with the following electric charge assignment:

$$S=\left(\begin{array}{ccc}\sigma _1^0\hfill & h_1^{}\hfill & h_2^+\hfill \\ h_1^{}\hfill & H_1^{}\hfill & \sigma _2^0\hfill \\ h_2^+\hfill & \sigma _2^0\hfill & H_2^{++}\hfill \end{array}\right),$$
(4)
and we will assume that only $`\sigma _2^0`$ gets a nonzero vacuum expectation value (VEV) in order to give the correct mass to the known charged leptons plus a mixing with the heavy leptons ($`K_a`$ and $`K_a^{}`$ terms below). The Yukawa interactions in the lepton sector are given by

$`_Y^l`$ $`=`$ $`H_{ab}^\nu \overline{\mathrm{\Psi }}_{aL}\nu _{bR}\eta +H_{ab}^l\overline{\mathrm{\Psi }}_{aL}S(\mathrm{\Psi }_{bL})^c+K_a\overline{\mathrm{\Psi }}_{aL}E_R\rho `$ (5)
$`+`$ $`K_a^{}\chi ^T\overline{E}_L(\mathrm{\Psi }_{aL})^c+G_E\overline{E}_LE_R\varphi +H.c.`$ (6)
where $`H_{ab}^l`$ is a symmetric matrix in the generation space; we have omitted $`SU(3)`$ indices. Neutrinos are strictly Dirac particles since the total lepton number will also be an automatically conserved.

Next we impose a $`Z_{13}`$ discrete symmetry under which the fields transform as $`Q_{iL}\omega _2^1Q_{iL}`$, $`Q_{3L}\omega _0Q_{3L}`$, $`u_{\alpha R}\omega _3u_{\alpha R}`$, $`d_{\alpha R}\omega _5^1d_{\alpha R}`$, $`J_R\omega _4J_R`$, $`j_{mR}\omega _6^1j_{mR}`$, $`\mathrm{\Psi }_L\omega _6\mathrm{\Psi }_L`$, $`E_L\omega _3E_L`$, $`\nu _R\omega _4^1\nu _R`$, $`E_R\omega _1E_R`$, $`\eta \omega _3^1\eta `$, $`\rho \omega _5\rho `$, $`\chi \omega _4^1\chi `$, $`S\omega _1^1S`$, $`\varphi \omega _2\varphi `$, where $`\omega _k=e^{2\pi ik/13},k=0,\mathrm{},6`$. Notice that if $`N`$ is a prime number the singlet $`\varphi `$ can transform under this symmetry with any assignment (but the trivial one), otherwise we have to be careful with the way we choose the singlet $`\varphi `$ to transform under the $`Z_N`$ symmetry. This symmetry implies that the lowest order effective operator that contributes to the axion mass is $`\varphi ^{13}/M_{\mathrm{Pl}}^9`$ which gives a mass of the order $`(v_\varphi )^{11}/M_{\mathrm{Pl}}^9`$ and also keeps the $`\overline{\theta }`$ parameter small as discussed in Ref. .

The most general scalar potential invariant under the gauge and $`Z_{13}`$ discrete symmetries is

$$V_{331}^{(A)}=V_\mathrm{H}+\left(\lambda _{\varphi 1}\varphi ϵ^{ijk}\eta _i\rho _j\chi _k+\lambda _{\varphi 2}\chi ^TS^{}\rho \varphi ^{}+\text{H. c.}\right),$$
(7)
where $`V_\mathrm{H}`$ denotes the Hermitian terms of the potential. This scalar potential has the correct number of Goldstone bosons and an axion field.

After imposing the $`Z_{13}`$ symmetry defined above we have that both the total lepton number $`L`$ and the PQ symmetry are automatic. The PQ charge assignment is as follows:

$`u_L^{}=e^{i\alpha X_u}u_L,d_L^{}=e^{i\alpha X_d}d_L,l_L^{}=e^{i\alpha X_l}l_L,`$ (8)
$`\nu _L^{}=e^{i\alpha X_\nu }\nu _L,j_L^{}=e^{i\alpha X_j}j_L,J_L^{}=e^{i\alpha X_J}J_L,`$ (9)
$`E_L^{}=e^{i\alpha X_{E_L}}E_L,E_R^{}=e^{i\alpha X_{E_R}}E_R,`$ (10)
and in the scalar sector we have the following PQ charges:

$`\eta ^0:`$ $`2X_u=2X_d=X_\nu X_{\nu _R},`$ (11)
$`\eta _1^{}:`$ $`(X_u+X_d)=X_u+X_d=X_lX_{\nu _R},`$ (12)
$`\eta _2^+:`$ $`(X_J+X_u)=X_j+X_d=(X_l+X_{\nu _R}),`$ (13)
$`\rho ^0:`$ $`2X_u=2X_d=X_lX_{E_R},`$ (14)
$`\rho ^+:`$ $`(X_u+X_d)=X_u+X_d=X_\nu X_{E_R},`$ (15)
$`\rho ^{++}:`$ $`(X_J+X_d)=X_j+X_u=(X_l+X_{E_R}),`$ (16)
$`\chi ^{}:`$ $`(X_u+X_J)=X_d+X_j=X_\nu +X_{E_L},`$ (17)
$`\chi ^{}:`$ $`(X_d+X_J)=X_u+X_j=X_{E_L}+X_l,`$ (18)
$`\chi ^0:`$ $`2X_J=2X_j=X_{E_L}X_l,`$ (19)
$`\varphi :`$ $`2X_j,`$ (20)
$`\sigma _1^0:`$ $`X_d+3X_j=2X_\nu ,`$ (21)
$`h_1^{}:`$ $`3X_jX_d=X_l+X_\nu ,`$ (22)
$`h_2^+:`$ $`4X_j=X_l+X_\nu ,`$ (23)
$`\sigma _2^0:`$ $`4X_j2X_d=0,`$ (24)
$`H_1^{}:`$ $`3(X_jX_d)=2X_l,`$ (25)
$`H_2^{++}:`$ $`5X_jX_d=2X_l.`$ (26)
From Eqs. (26) we obtain the relations $`X_j=X_J=\frac{1}{2}X_d=\frac{1}{2}X_u=\frac{2}{3}X_l=\frac{2}{3}X_{\nu _R}=\frac{2}{5}X_\nu =\frac{2}{5}X_{E_R}=2X_{E_L}=\frac{1}{2}X_\varphi `$.

Notice that the mass scale related to the exotic charged lepton $`E`$, up to an arbitrary dimensionless constant $`G_EO(1)`$, is related to the PQ energy scale since the requirement of the symmetries of the model imposes the Yukawa couplings in Eq. (6).

Notice also that, at the energy scale below the breakdown of the $`SU(3)_L`$ symmetry, this model has scalar multiplets transforming under $`SU(3)_CSU(2)_LU(1)_Y`$ as follows: two doublets $`(\rho ^+,\rho ^0)(\mathrm{𝟏},\mathrm{𝟐},+1)`$, $`(\eta ^0,\eta ^{})(\mathrm{𝟏},\mathrm{𝟐},1)`$ and a non-hermitian triplet $`(H^{},h_1^{},\sigma _1^0)(\mathrm{𝟏},\mathrm{𝟑},2)`$. With the lighter scalar multiplets, and the usual degrees of freedom, the energy scale at which $`\mathrm{sin}^2\theta _W=0.25`$ is 5.2 TeV. The doublets $`(\chi ^{},\chi ^{})(\mathrm{𝟏},\mathrm{𝟐},2)`$, $`(h_2^+,h^0)(\mathrm{𝟏},\mathrm{𝟐},+1)`$ and the extra vector bosons have masses proportional to $`v_\chi `$; the lepton singlet $`E`$ has a mass of the order of $`v_\varphi `$. More details will be given elsewhere.

### B Model with three scalar triplets (Model B)

In this model we do not introduce the scalar sextet and the Yukawa interactions are

$`_Y^l`$ $`=`$ $`H_{ab}^\nu \overline{(\mathrm{\Psi })}_{aL}\nu _{bR}\eta +H_{ab}^lϵ_{ijk}\overline{(\mathrm{\Psi })^c}_{iaL}\mathrm{\Psi }_{jbL}\eta _k`$ (27)
$`+`$ $`K_a\overline{\mathrm{\Psi }}_{aL}E_R\rho +K_a^{}\chi ^T\overline{E_L}(\mathrm{\Psi }_{aL})^c+G_E\overline{E}_LE_R\varphi ^{}`$ (28)
$`+`$ $`H.c.,`$ (29)
where $`H_{ab}^l`$ is now an antisymmetric matrix. In both Yukawa interactions above, a general mixing is allowed in each charge sector. As in the previous model, neutrinos are strictly Dirac particles. The charged leptons gain mass as in Ref. .

If we want to implement a given texture for the quark and lepton mass matrices we have to introduce more scalar triplets, and a larger $`Z_N`$ symmetry will be possible in the model.

Let us introduce a $`Z_4`$ symmetry with parameters denoted by $`\stackrel{~}{\omega }_0`$, $`\stackrel{~}{\omega }_1`$, $`\stackrel{~}{\omega }_1^1`$, and $`\stackrel{~}{\omega }_2\stackrel{~}{\omega }_2^1`$. $`u_{\alpha R}`$, $`Q_{iL}`$, and $`\nu _{aR}`$ transform with $`\stackrel{~}{\omega }_1`$; $`d_{\alpha R}`$, $`Q_{3L}`$, $`\mathrm{\Psi }_{aL}`$, $`E_R`$, $`\chi `$, and $`\varphi `$ transform with $`\stackrel{~}{\omega }_1^1`$, $`\eta `$ transform with $`\stackrel{~}{\omega }_2`$, and all the other fields remain invariant, i.e., transform with $`\stackrel{~}{\omega }_0`$. After $`Z_4`$ is imposed, the total lepton number $`L`$ and the PQ and $`Z_{11}`$ symmetries are all automatically implemented in the Yukawa sector and in the scalar potential. The most general scalar potential is then

$$V_{331}^{(B)}=V_\mathrm{H}+\left(\lambda \varphi ϵ^{ijk}\eta _i\rho _j\chi _k+\text{H. c.}\right).$$
(30)
The following $`Z_{11}`$ symmetry is automatically implemented in both the Yukawa interactions and in the scalar potential: $`Q_{iL}\omega _3Q_{iL}`$, $`Q_{3L}\omega _0Q_{3L}`$, $`u_{\alpha R}\omega _4u_{\alpha R}`$, $`d_{\alpha R}\omega _1^1d_{\alpha R}`$, $`J_R\omega _5^1J_R`$, $`j_{mR}\omega _3^1j_{mR}`$, $`\mathrm{\Psi }_L\omega _2\mathrm{\Psi }_L`$, $`E_L\omega _3E_L`$, $`\nu _R\omega _5^1\nu _R`$, $`E_R\omega _1E_R`$, $`\eta \omega _4^1\eta `$, $`\rho \omega _1\rho `$, $`\chi \omega _5\chi `$, $`\varphi \omega _2^1\varphi `$. It happens that, in addition to the $`Z_{11}`$ symmetry, the $`U(1)_{\mathrm{PQ}}`$ and the conservation of the total lepton number are also automatic i.e., a consequence of the gauge symmetry and renormalizability of the model, in the interactions in Eqs. (3), (29), and (30). The PQ charge assignments for the fermions in the model are as in Eq. (10); and in the scalar sector we have constraints equations as in the previous subsection. In this case proceeding as in the model A, we obtain the relations $`X_d=X_u=0`$ and $`X_l=X_{\nu _R}=X_\nu =X_{E_R}=\frac{1}{2}X_j=\frac{1}{2}X_J=\frac{1}{3}X_{E_L}=\frac{1}{4}X_\varphi `$. Notice that for leptons the PQ transformations are not a chiral symmetry. As in Model A, we see from Eq.(29) that the mass scale related to the singlet charged lepton $`E`$ is related to $`v_\varphi `$. Moreover, in this model we have that $`\mathrm{sin}^2\theta _W=0.25`$ at 4 TeV.

## III conclusions

We have built two invisible axion models in which the axion is naturally light and protected against quantum gravity effects. In model A, the $`Z_{13}`$ symmetry has to be imposed but in model B the $`Z_{11}`$ symmetry is automatically implemented in the classical Lagrangian after imposing a $`Z_4`$ symmetry. With a $`Z_{13}`$ symmetry the axion is protected from gravitational effects even if $`v_\varphi 10^{12}`$ GeV but, with a $`Z_{11}`$ symmetry, $`v_\varphi 10^{10}`$ GeV. In both models the PQ symmetry is automatically implemented in the classical Lagrangian in the sense that it is not imposed on the Lagrangian but is just a consequence of the particle content of the model, its gauge invariance, renormalizability, and Lorentz invariance.

We would like to stress the strong constraint put on model building by the approach proposed in Refs. . Once the symmetry $`Z_N`$ is used, automatic or imposed, there is no choice for new interactions. In this vein, in both models neutrinos are strictly Dirac particles, and for this reason both models will be ruled out if the neutrinos turn to be Majorana particles, say by observation of the neutrinoless double beta decay.

In the PQ solution to the strong CP problem the quark contributions to $`\overline{\theta }`$ are such that $`\overline{\theta }\overline{\theta }2\alpha _fX_f`$, where $`f`$ denotes any quark. In both models we have $`\overline{\theta }\overline{\theta }2\alpha X_j`$ and we have the domain wall problem related to $`Z_2U(1)_{\mathrm{PQ}}`$ for $`X_j=2`$.

Concerning the CDM, we would like to call attention to a new possible candidate: a light and stable scalar in nonsupersymmetric models . Although the mass of this scalar is in the range $`32\mathrm{GeV}m<M_Z`$, the model is still compatible with the LEP data since the lightest scalar is almost a singlet under $`SU(2)_LU(1)_Y`$, and it may be mistaken for a light, $`m_\chi 50`$ GeV , or for the usual, $`m_\chi \stackrel{>}{}50`$ GeV, neutralino. We recall that the latter bound comes from LEP2 searches for the corresponding chargino $`m_{\chi ^\pm }\stackrel{>}{}100`$ GeV.

This work was supported by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) and partially by Conselho Nacional de Ciência e Tecnologia (CNPq).