text:
# Tracing the power-law component in the energy spectrum of black hole candidates as a function of the QPO frequency

## 1 Introduction

In 1992, the WATCH all-sky monitor onboard GRANAT discovered the X-ray source GRS 1915+105 as a bright transient (Castro-Tirado et al. 1992). Mirabel & Rodríguez (1994) reported the observation of radio relativistic expansions interpreted in terms of collimated emissions of matter, which made GRS 1915+105 the first galactic object with superluminal jets, for which the term “microquasar” was coined. This source was originally supposed to host a black hole because of its high X-ray luminosity, and for its similarity with another superluminal source, GRO J1655-40, for which it was possible to obtain a dynamical estimate of the mass (M=6.3$`\pm `$0.5M) of the compact object (Bailyn et al. 1995; Greene et al. 2001). No optical counterpart of GRS 1915+105 was found because of the high Galactic extinction and only an IR counterpart was detected (Mirabel et al. 1994). Recently, an orbital period of 33.5 days was discovered with IR spectroscopic observations, leading to a dynamical estimate of 14$`\pm `$4 M for the compact object (Greiner et al. 2001, see also Borozdin et al. 1999). The numerous observations of GRS 1915+105 made with the instruments onboard the Rossi X-ray Timing Explorer (RXTE) satellite show an extremely complex variability in the X-ray band, which Belloni at al. (1997a, 1997b) interpreted as an oscillation of the inner region of the accretion disk caused by a thermal-viscous instability. Belloni et al. (2000) introduced a model-independent classification of the different types of X-ray variability in three basic spectral states. Similarly to more conventional Black Hole Candidates (BHC), the energy spectrum of GRS 1915+105 can be fitted with the superposition of a soft thermal component, modeled as a disk-blackbody (Mitsuda et al. 1984), plus a harder component, modeled as a power law with a high energy cutoff (Tanaka & Lewin 1995). The state transition can be attributed to changes in the relative contributions of these variable components, again similarly to what observed in other BHCs (Belloni 1998).

Detailed timing analysis of the RXTE data (Morgan et al. 1997) revealed three types of Quasi Periodic Oscillations (QPO): a 1-10 Hz QPO with variable centroid frequency similar to that observed in other BHCs (see van der Klis 1995), two quasi-stable QPOs with centroid frequencies $``$40 and $``$67 Hz (Morgan et al. 1997; Belloni et al. 2001; Strohmayer et al. 2001a), and variable low frequency QPOs ($`10^31`$ Hz), which can be identified with transitions between the states of Belloni et al. (2000). The presence of the 1-10 Hz is always associated with a spectral state C (Belloni et al. 2000). The 1-10 Hz QPO centroid frequency is well correlated with both disk and power-law spectral components (Muno et al. 1999, Reig et al. 2000, Rodriguez et al. 2002) and with the X-ray flux (Chen et al. 1997; Markwardt et al. 1999).

Several theoretical models associate the origin of this QPO with the innermost region of the disk (Miller et al. 1998; Nowak & Wagoner 1991; Titarchuk et al. 1998; Stella & Vietri 1998; Stella et al. 1999; Varnière et al. 2002, Psaltis & Norman 2002). Di Matteo & Psaltis (1999) made use of the fact that the fastest variability timescale at any radius around a compact object is the Keplerian orbital frequency to obtain upper limits on the radius at which the 1-10 Hz QPO is produced. As it is apparent in their work, there is a general positive correlation, considering different systems, between the photon index of the power-law component in the energy spectrum and the centroid frequency of the QPO. Sobczak et al. (2000a) in their analysis of RXTE data of two transient BHCs, XTE J1550-564 and GRO J1655-40, examined the correlation between the spectral parameters, inner disk radius and power-law photon index, with the QPO centroid frequency. Surprisingly, they find correlations of opposite sign for the two sources. Rodriguez et al. (2002) analyzed the correlation between QPO centroid frequency and inner disk radius (determined from X-ray spectral fits).

In this paper, we concentrate on the correlation between power-law photon index and QPO centroid frequency. We present the results of timing and spectral analysis of a number of RXTE observations of GRS 1915+105, focusing on this correlation. We also consider published data relative to four other systems: GRO J1655-40 (Remillard et al. 1999a; Sobczak et al. 2000a), XTE J1550-564 (Sobczak et al. 2000a,b; Homan et al. 2001), XTE J1748-288 (Revnivtsev et al. 2000), 4U 1630-47 (Trudolyubov at al. 2001; Tomsick & Kaaret 2000). We then compare the results in term of possible theoretical models.

## 2 Data analysis

Our study of GRS 1915+105 is divided in two parts, according to the different classes of observations that we considered. In the first part, we analyze the so-called plateau observations (Fender et al. 1999; Reig at al. 2000), which consist in long intervals of persistent state C (class $`\chi `$ in Belloni et al. 2000). The second part is dedicated to observations that show large variability in the X-ray flux and at least one state-C interval. Again following the classification of Belloni et al. (2000), we selected observations corresponding to classes $`\alpha `$, $`\beta `$, and $`\nu `$.

### 2.1 Plateau observations

Plateau intervals show persistent QPOs that are very easily detectable in the Power Density Spectrum (PDS). We based our choice of appropriate observations upon the extensive analysis of Reig et al. (2000), which allowed us to identify the larger possible range of centroid frequencies of the QPO. We selected 7 RXTE/PCA observations of GRS 1915+105 (for a total of 15 satellite orbits), corresponding to class $`\chi `$ in the classification of Belloni et al. (2000). The observation log is reported in Table 1.

The plateau observations are relatively straightforward to analyze from the point of view of timing analysis, because of the stable positive correlation between the QPO frequency and count rate (see Reig et al. 2000). Therefore, no special timing selection is necessary and it is possible to accumulate all the data to produce one PDS and one energy spectrum. To produce the PDS, for each observation we used high-time resolution data in the 2-13 keV range. Some of the observations have a time resolution of 1/512 s, limited to 1/128 s for the others (see Table 1). We divided each observation interval (orbit of the satellite) in intervals 16 s long and produced an average PDS per observation summing the PDS of the individual segments. We subtracted the contribution due to Poissonian statistics and the Very Large Event window (Zhang et al. 1995, Zhang 1995) from each PDS and renormalized them to squared fractional rms per Hertz (see Belloni & Hasinger 1990). Every PDS was then logarithmically rebinned and fitted to find the centroid frequency of the 1-10 Hz QPO, which is very sharp and often shows higher harmonics and even a sub-harmonic component (as shown in Fig. 1 for three different observations). A model consisting of the superposition of a number of Lorentzians (from 5 to 9 depending on the observation) was used to fit every feature present in the power spectrum (see Belloni et al. 2002). Because of the relative strength of the QPO peak, the centroid frequency did not change appreciably when using different models for the continuum noise: the QPO frequencies we found with this analysis are in close agreement with those from Reig at al. (2000). The resulting QPO centroid frequencies are listed in Table 1.

In Class $`\chi `$ observations, the source is in a long steady state C (Belloni et al. 2000). As shown by Belloni et al. (2000), this corresponds to an undetectable thermal component in the PCA energy spectrum. For each observation interval we extracted PCA energy spectra in the range 3–25 keV, averaged over the whole interval, from spectral (Standard2) data. For each spectrum, we subtracted the background estimated with pcabackest v. 2.1e, and created an appropriate detector response matrix with pcarsp v. 2.43. A systematic error of 1% was added to the spectra to account for remaining uncertainties in the detector calibration. Our starting model for the spectral fits consisted of a power law, corrected for interstellar absorption, plus a Gaussian emission line to take into account an excess at 6.4 keV. With this model, it was not possible to lower the reduced $`\chi ^2`$ below 2. We modified the model, adding a high-energy cutoff to the power law. The resulting spectral parameters are shown in Table 1. This model is consistent with that used by Trudolyubov et al. (1999) and Muno et al. (1999). Notice that Trudolyubov et al. (1999) found rather high cutoff energies ($``$100 keV), while Muno et al. (1999) found also much lower values, even below 10 keV. Our values are in the range 20–30 keV.

Only 2 among the selected plateau observations (10408-01-27-00 and 20402-01-52-00) were also analyzed by Muno et al. (1999). They used an often-used model for Black Hole Candidates, which consists of a superposition of a multicolor disk-blackbody (thermal component) and a power law (without high energy cutoff). Their spectral analysis results in a different power law photon index, $`0.7`$ grater than the values we found. However, they derive an extremely high inner disk temperature ($`4`$ keV) and an inner disk radius consistent with zero ($`2\pm 2`$ km). These parameters do not represent a realistic physical picture of the source and therefore suggest that the inclusion of a disk-blackbody component is not physically necessary. Our choice of a simple power law with a high-energy cutoff is motivated by two reasons: the need of a simple model for efficiently performing a large number of fits, and the need to use a common model that allows comparison with values in the literature. In order to check that for high values of the spectral index $`\mathrm{\Gamma }`$ this model yields meaningful values, we fitted our spectra also with the bmc model in XSPEC, which has been used by other authors to describe the energy spectra of Black Hole Candidates in the framework of Comptonization models. These fits give values of $`\mathrm{\Gamma }`$ consistent with those from a simple power-law model, with a slight systematic excess of 3%.

From Table 1, one can see that there are variations in the best fit values of the column density ($`\mathrm{N}_\mathrm{H}`$) (see also Belloni et al. 2000; Klein-Wolt et al. 2002).

In Fig. 2 we plot the power-law photon index versus the QPO centroid frequency. A positive correlation is evident.

### 2.2 Class $`\alpha `$, $`\beta `$ and $`\nu `$ observations

In order to examine the correlation between QPO centroid frequency and power-law photon index, and in particular to extend it to higher values of the QPO frequency, we also analyzed observations in which GRS 1915+105 shows large variability on time scales longer than 1 second. While belonging to different classes, they all have in common the presence of at least one state-C interval of more than 700 seconds (see Belloni et al. 2000). The log of selected observations is reported in Table 2. After a close inspection of the light curves produced from Standard1 data (corresponding to the full PCA energy band 2–60 keV), we identified the start and end times of the state-C intervals to analyze. Flux values during these periods typically range from 3000 to 40000 cts/s. For each observation, we checked the choice of intervals by producing a dynamical PDS (a series of PDS as a function of time), with which we verified that a 1–10 Hz QPO was present throughout the whole interval. The dynamical PDS also confirmed that the X-ray intensity of the source is strongly correlated with the QPO frequency (Chen et al. 1997, Markwardt et al. 1999, Reig et al. 2000). Therefore, following the PCA count rate, the QPO frequency first decreases rather fast, then starts increasing at a lower rate (see Fig. 4).

To perform timing analysis, we inspected the same type of data (high-time resolution and low energy range) as for plateau observations. In this case however, we divided every state-C interval in segments with a length of 4 s rather than 16, in order to better follow the rapid variations of the QPO centroid frequency, and produced a PDS for every segment. After subtracting the contribution from the Poissonian statistics and the Very Large Event window, we summed these PDS in sets of four, rebinned them logarithmically and fitted the result to find the QPO centroid frequency. As in the previous case, a model consisting of the superposition of a number of Lorentzians was used (see Belloni et al. 2002). We performed the summing process because we wanted to correlate timing and spectral results (which have a minimum time bin of 16 s), and this is possible only if we analyze intervals with temporal coincidence.

Spectral parameters were obtained extracting energy spectra from Standard2 data, with an integration time of 16 seconds, subtracting the background and creating the detector response matrix as for the plateau observations. We fitted every spectrum using the “standard” model for BHCs: the superposition of a multicolor disk-blackbody and a power law for the high energy tail. We corrected for the interstellar absorption, fixing $`\mathrm{N}_\mathrm{H}=7\times 10^{22}`$ cm<sup>-2</sup> (Klein-Wolt et al. 2002) and added a Gaussian emission line with central energy fixed to 6.4 keV. With respect to the model adopted for the plateau observation, we needed to include a disk-blackbody (see Belloni et al. 1997b, 2000) and we had to fix $`\mathrm{N}_\mathrm{H}`$ due to the limited statistics of the 16s spectra. Again, we added in quadrature a systematic error of 1%. The fits give acceptable results, with a reduced $`\chi ^2`$ usually around 1. Therefore, we can trace the evolution (with a time step of 16 s) of the spectral parameters during the state-C intervals. We can then compare the spectral and timing evolution by correlating the resulting parameters. As in the case of plateau observations we fitted selected spectra with the bmc model (see Sect. 2.2), finding consistent results for $`\mathrm{\Gamma }`$.

Examples of both PDSs and energy spectra (with their best fit models) can be seen in Fig. 3. The correlation between QPO centroid frequency and power-law index for classes $`\beta `$ and $`\nu `$ is plotted in Fig. 5 (upper panel). The behaviour is similar to $`\chi `$ observations, although there is a spread due to the lower statistics available from spectra averaged over only 16 seconds. At high QPO frequencies it is possible to notice a “turnoff” in the correlation for observations 15 and 16, while observation 19 (diamonds) does not show this feature. Figure 5 (lower panel) shows the same correlation from classes $`\alpha `$ and $`\nu `$. While $`\nu `$ data are compatible with those in the upper panel, $`\alpha `$ data show a correlation similar, but offset to lower photon-index values.

## 3 Other sources

In order to compare our results for GRS 1915+105 with other systems, we searched the literature for similar data from other black hole candidates. We identified four systems for which values for the QPO centroid frequency and the power-law index are available, and for which there are at least a few different values. Notice that, as mentioned above, we do not consider parameters like inner disk radius, which are difficult to compare between different sources. Providing the adopted spectral models are compatible, the values of $`\mathrm{\Gamma }`$ are directly comparable, and so of course are the QPO centroid frequencies.

### 3.1 GRO J1655-40

This source is a BHC with a dynamical estimate for the mass of the black hole of $`6.3\mathrm{M}_{}`$ (Bailyn et al. 1995; Greene et al. 2001). The system, a bright X-ray transient which appeared for the first time in 1994 and showed a subsequent outburst in 1996/1997, showed superluminal radio jets (Hjellming & Rupen 1995). GRO J1655-40 is therefore another of the few objects belonging to the microquasar class. This object shows quite a few different QPO peaks in a range between 0.1 and 450 Hz, (see Remillard et al. 1999a; Strohmayer 2001, Remillard et al. 2002b). Among them, we are interested in the variable 14–22 Hz QPO (see Sobczak et al. 2000a). We consider here PCA/RXTE data from an intensive observational campaign of the source during its second outburst in 1996–1997. The timing parameters (i.e., the QPO centroid frequencies) are taken from Table 3 of Remillard et al. (1999a), and the spectral parameters (i.e., the power law photon indexes) from Table 1 of Sobczak et al. (1999a). The energy spectra in the 2.5–20 keV energy range were fitted using a model consisting of a multicolor disk-blackbody plus a power law; a detailed description of the spectral analysis can be found in Sobczak et al. (1999a). The power density spectrum for each observation was analyzed in the 2–30 keV energy range, obtaining the QPO central frequency fitting with a Lorentzian function plus a power law for the local continuum (see Remillard et al. 1999a). The plot of power-law index versus QPO centroid frequency for GRO J1655-40 is shown in Fig. 6 (squares). There is a clear anti-correlation between these two parameters, as already shown by Sobczak et al. (1999a).

### 3.2 XTE J1550-564

This BHC is the brightest X-ray transient observed with RXTE (the 1998–1999 flare reached 6.8 Crab in the 2–10 keV band). In the PDS, it shows both low-frequency (0.8–18 Hz) and high-frequency (180–280 Hz) QPOs (Cui et al. 1999; Remillard et al. 1999b; Homan et al. 2001; Remillard et al. 2002a,b). We consider here the low-frequency one (see Remillard et al. 2002a). The observations presented here are the whole set made with PCA/RXTE before 1999 March 22 (Gain Epoch 3). The power-law photon indexes obtained from the spectral analysis of these observations are taken from Table 1 of Sobczak et al. (1999b), while the QPO centroid frequencies from the PDS analysis are listed in Table 1 of Sobczak et al. (2000a). The spectral model used by Sobczak et al. (1999b) to fit the 2.5–20 keV energy spectra consists again of the superposition of a multicolor disk-blackbody and a power law, with corrections for interstellar absorption and a Gaussian emission line at $`67`$ keV. The centroid frequency of the QPO was found fitting the peak with a Lorentzian function and the continuum with a power law.

Data obtained from the literature for XTE J1550-564 are plotted in Fig. 6 (circles). They show a positive correlation between power law photon index and QPO centroid frequency very similar to the behaviour of GRS 1915+105 (see Fig. 2 and 5). However, while the correlation is very good at low frequency, at higher QPO values there is a larger spread. Surprisingly, the points from GRO J1655-40, although showing a negative correlation, are placed on the natural continuation of the upper points for XTE J1550-564. As one can see from Cui et al. (1999), Remillard et al. (2002b) and Homan et al. (2001), the low-frequency QPO in XTE J1550-564 does not maintain the same properties throughout the whole outburst. It is therefore natural to examine its behaviour in separate intervals of the light curve.

We subdivide the light curve in four separate sections (see Fig. 7): a ‘rise’ interval from the first observation up to the peak flux, a ‘down’ interval from the peak following the smooth decrease in flux until the first small flare, a ‘flat’ interval where the flux remains rather stable, with small flares, and a ‘decay’ interval where the flux decays at the end of the first part of the outburst. After the last observation examined here, no QPO was detected for a long time, until the flux rose again to the second part of the outburst (see Homan et al. 2001).

In Fig. 8 the behaviour of XTE J1550-564 in the QPO-$`\mathrm{\Gamma }`$ plane is shown with different symbols for the four different intervals. One can see that the first two intervals trace a very smooth curve in Fig. 8, the second interval tracing back the path of the first. The scattered points correspond to the remaining two intervals, and it is interesting to note that all these points are located below the others. What is different between the first two and the second two classes of observations that have a different behaviour in this plot? From Sobczak et al. (2000b) one can see that a clear dividing line between the two groups is the ratio of the power-law flux over the total flux, estimated from the observed PCA spectra. The first observations have a lower value of the ratio. Also notice that the points that contribute to the scatter in Fig. 8 are the ones from the ‘decay’ interval, which contains QPOs classified as A/B in Homan et al. 2001, while the QPO observed in the previous intervals is classified as C (notice that these QPO classes have nothing to do with the A/B/C classes defined in Belloni et al. 2000).

### 3.3 XTE J1748-288

This BHC shows superluminal radio jets and belongs to the class of galactic microquasars, for the same reasons as GRS 1915+105 and GRO J1655-40 (Smith et al. 1998; Hjellming et al. 1998; Revnivtsev et al. 2000; Kotani et al. 2000; Miller et al. 2001). During the 1998 outburst the source was observed to display spectral and timing properties typical for BHCs systems (Revnivtsev, et al. 2000). A QPO feature was found in the PDS at 0.5 Hz and 20–30 Hz. Here we consider spectral and timing parameters given in Table 2 and 3 of Revnivtsev et al. (2000) who reported an analysis of the PCA/RXTE observations during the 1998 outburst. They fitted the energy spectra between 3 and 25 keV using the standard model for BHCs consisting of a multicolor disk-blackbody plus a power law, with low energy interstellar absorption and a Gaussian emission line at $`6.5`$ keV. For the seven observations considered here, Revnivtsev et al. (2000) fitted the PDS using the sum of a flat-topped band limited noise component, a power law component, and a Lorentzian function to describe the QPOs. The resulting points in the QPO-$`\mathrm{\Gamma }`$ plot are shown as empty diamonds in Fig. 9. There is a clear positive correlation between $`\mathrm{\Gamma }`$ and QPO, but it is not identical to the same correlations of the three sources previously discussed. Notice that the QPO frequencies are the highest considered here.

### 3.4 4U 1630-47

This X-ray transient is considered a BHC for which the optical thin radio emission, observed during the 1998 outburst, suggests the presence of radio lobes (Hjellming et al. 1999). We obtained timing (QPO centroid frequency) and spectral (power-law index) parameters of the observations made during the decay of the 1998 outburst from Table 4 and 5 of Tomsick & Kaaret (2000). The timing analysis was made in the 2–21 keV energy band fitting the PDS with a flat-top component for the continuum and a Lorentzian function for the QPO feature. The centroid frequency was between 0.2 and 3.4 Hz. The energy spectra were fitted in the 2.5–20 keV range using a power law, a multicolor disk-blackbody soft component, a Gaussian emission line at $`7`$ keV, the interstellar absorption and a broad iron absorption edge. For the peak of the outburst, we used the results of Trudolyubov et al. (2001). Here timing analysis was performed from data in the 2–13 keV energy range and the PDS were fitted with Lorentzian models plus a low frequency power-law. The PDS were rather complex: as centroid frequency we used the peak at the lowest frequency in the 1–20 Hz range. The 3-20 keV PCA energy spectra were fitted with a disk-blackbody plus a power law, although the presence of an excess around 6–8 keV was recognized.

The points corresponding to this source are shown in Fig. 9. It is possible to recognize a behaviour similar to those shown by GRS 1915+105 and XTE J1550-564, although in this case the correlation is not as smooth. The points corresponding to the end phase of the outburst have a rather constant $`\mathrm{\Gamma }`$ around 1.7–1.8, while the QPO frequency decreases from 3 Hz down to 0.2 Hz (data from Tomsick & Kaaret 2000). Notice that the points corresponding to the lowest frequencies (0.2 to 0.5 Hz) do match rather well the correlation followed by XTE J1550-564. At higher QPO frequencies there is a correlation between the two quantities, with a rather large scatter and a possible indication of a turn-off.

## 4 Summary and discussion

Our results can be summarized as follows. For the $`\chi `$ class of GRS 1915+105, when the energy spectrum is rather hard and consists exclusively of the power-law component and the power density spectrum shows a clear 1–10 Hz QPO, there is a strong positive correlation between QPO centroid frequency and photon index of the power law. At the highest values of the QPO centroid frequencies, this correlation flattens and appears to start reversing. Published data from other sources indicate that this general correlation exists for a number of sources, and indeed in some cases there is a turnoff in the relation for high QPO frequencies. The correlation is not exactly the same for all sources and in the case of XTE J1550-564 there is evidence for different turnoffs during different intervals of the source outburst.

The turnoff in the correlation is difficult to understand. If the QPO frequency is associated to a specific radius in the accretion disk, only the rising part of the correlation can be explained with a simple Comptonization model. If the QPO-radius relation also shows a reversal (see for instance Rodriguez et al. 2002 and Varniére et al. 2002), then a turnoff like that presented here can be produced.

The source for which we have more data, and the largest range in QPO frequencies is XTE J1550-564. From Fig. 8 one can see that not only the points from the first two sections of the outburst (see Fig. 7) follow a very narrow correlation, but also that all subsequent points lie below this correlation. This suggests that the correlation followed in the first part, when the overall PDS was in the form of a flat-top noise component plus a QPO (Cui et al. 1999; Remillard et al. 2002b), is perturbed in a very specific way when the source changes state and a powerful disk component appears in the energy spectrum (Sobczak et al. 1999b,2000a,b). There is increasing evidence that the power-law component observed in the energy spectra of BHCs in the Low/Hard state is physically different from that seen in the Very High State (see e.g. Zdziarski et al. 2001). It is therefore possible that the turnoff is due to a gradual transition between these two regimes.

It is interesting to discuss the position of the turnoff in the QPO-$`\mathrm{\Gamma }`$ relation. Figure 10, where all the sources are shown together, suggests that: (a) XTE J1550-564 and GRO J1655-40 have a similar turnoff frequency around 10 Hz; (b) XTE J1748-288 has a higher turnoff frequency, if any, larger than 30 Hz; (c) GRS 1915+105 (see Fig. 5) has a lower turnoff frequency at about 5 Hz. It is tempting to compare these values with the relative values of the estimated dynamical masses for the compact objects in the systems. Excluding XTE J1748-288, for which no dynamical mass estimate is available, we see that GRS 1915+105 has a mass roughly twice and a turnoff frequency half that of GRO J1655-40, as one could expect with a simple mass scaling. However, the recent mass estimate for XTE J1550-564 is around 10 (Orosz et al. 2002), from which we would expect a turnoff frequency lower by a factor of $``$2 with respect to that of GRO J1655-40, which is evidently not observed.

In conclusion, we have found that in GRS 1915+105 and a number of other BHCs, the centroid frequency of the low-frequency QPO and the photon index of the power-law component in the energy spectra follow a very similar kind of correlation. In particular, for some sources there is evidence for the presence of a turnoff in this correlation. If the QPO frequency is associated to a specific radius of the accretion disk, this provides an important clue to the possible origin of the power-law component. More data are needed in order to understand the physical nature of the turnoff, but these results stress the importance of such simple spectral/timing correlations for black hole candidates.

###### Acknowledgements.

TB acknowledges the hospitality of MIT and thanks the Cariplo Foundation for financial support.


format: text

text:
# Calculations of the Far-Wing Line Profiles of Sodium and Potassium in the Atmospheres of Substellar-Mass Objects

## 1 Introduction

Absorption lines of the neutral alkali metals sodium (Na), potassium (K), cesium (Cs), rubidium (Rb), and lithium (Li) are prevalent in the spectra of L dwarfs, T dwarfs, and irradiated giant planets (Kirkpatrick et al. 1999,2000; Burgasser et al. 1999,2000,2002; Burrows, Marley, and Sharp 2000 (BMS); Burrows et al. 2001,2002; Liebert et al. 2000; Charbonneau et al. 2001; Martín et al. 1999). The atmospheres of such substellar mass objects (SMOs) are cool enough that the neutral alkali atoms predominate. Furthermore, since aluminum, magnesium, iron, silicon, and calcium are sequestered in grains that have settled out, a cool atmosphere in the gravitational field of such a planet or star is depleted of these more refractory metals (Burrows and Sharp 1999; Lodders 1999; BMS). As a consequence, the less-refractory neutral alkali metals assume an importance in the optical and near-infrared spectra of SMOs and in their atmospheres from $``$800 K to $``$2000 K that is unique among astronomical objects.

In standard stellar atmospheres, atomic lines are superposed on a background continuum and the concepts of equivalent width and curve of growth make conceptual and practical sense. An individual line is but a perturbation on the local spectrum. However, due to the rainout of metals in cool atmospheres and the consequent paucity of continuum and alternate opacity sources between 0.4 $`\mu `$m and 1.0 $`\mu `$m , the wings of the strong resonance doublets centered at $``$7700 Å (K I) and 5890 Å(Na-D) assume the role of the continuum throughout most of this spectral range (BMS; Tsuji, Ohnaka, and Aoki 1999). In particular, the red wing of the $`4s^2S_{1/2}4p^2P_{3/2}`$ transitions of K I provides the pseudo-continuum in cool molecular atmospheres all the way from 0.77 $`\mu `$m to $``$1.0 $`\mu `$m and the Na-D doublet, centered as it is in the middle of the visible, determines the true color of brown dwarfs (magenta/purple; Burrows et al. 2001).

Hence, whereas in traditional stellar atmospheres the Lorentzian core and Gaussian wings of a line are not generally of relevance beyond $``$20 Å detunings ($`\mathrm{\Delta }\lambda `$ from the line core), in cool substellar atmospheres the relevant reach of the Na I and K I resonance lines can be thousands of Å. Given this, to achieve accurate spectral fits for brown dwarf, L dwarf, T dwarf, and hot giant planet atmospheres, the shapes of the far wings of these alkali lines as a function of pressure and temperature must be ascertained.

In this paper, we perform ab initio calculations of the energy shifts of the ground and excited states of sodium and potassium immersed in H2\- and helium-rich atmospheres and obtain the resulting opacity profiles of the red and blue wings of the Na-D and K I resonance lines. In §2, we review the general theory of line profiles in the quasi-static limit using the Unified Franck-Condon (UFC) formalism (Szudy and Baylis 1975,1996). In §3, we describe our use of the quantum chemical code GAMESS (Schmidt et al. 1993) to obtain the interaction potentials as a function of H2 and helium perturber distance and H2 orientation. The derived potential curves are presented and described in §4, where we also compare to related work in the literature. The resultant cross sections as a function of wavelength are given in §5. The latter section contains our major results. We compare our new theory with the approximate theories found in Burrows, Marley, and Sharp (2000) and Burrows et al. (2002) and find both similarities and important differences. In §6, we summarize our major conclusions.

## 2 Quasi-Static Theory of Absorption Spectra in the Far Wings of Atomic Spectral Lines

The energy levels of an atom of sodium or potassium immersed in a sea of molecular hydrogen will be perturbed by the potential field of the diatomic hydrogen. At a distance of $`R_i`$, the molecular hydrogen will shift both the ground and the excited states of the alkali metal atom. The difference of these level shifts is the amount by which the corresponding absorption line will be shifted (Griem 1964; Breene 1957,1981):

$$h\nu =h\nu _0+V_m(R_i)V_n(R_i)=h\nu _0+\mathrm{\Delta }V_{mn}(=\mathrm{\Delta }V),$$
(1)
where $`\nu `$ is the photon frequency, $`\nu _0`$ is the unperturbed photon frequency, $`m`$ and $`n`$ indicate the excited and ground states, and for an H2 perturber the potential shifts $`V`$ are a function of the orientation (angle) of the molecule and the distance. Note that $`V_m(R_i)`$ and $`V_n(R_i)`$ are the potential shifts with respect to the asymptotic shifts at infinity. $`\mathrm{\Delta }V`$ is a central quantity in the quasi-static theory of line profiles.

The quasi-static theory of absorption strength assumes that the perturbers are arrayed about the transitioning atom in a static, thermal distribution with a density $`N_p`$ given by its partial pressure and temperature ($`T`$) (Holtzmark 1925; Holstein 1950). A full quantum mechanical treatment of line shapes has been developed, but it has been applied to only a limited number of systems and is significantly more involved (Szudy and Baylis 1996). Fortunately, at least for the Na-He system, Pontius and Sando (1983) have shown that the semi-classical/quasi-static and the quantal methods yield very similar results.

In the quasi-static theory one imagines a shell of radius $`R_i`$ and thickness $`dR_i`$ centered at the atom. The strength of the absorption line in the wings at the frequency given by eq. (1) is then proportional to the number of perturbers in that shell:

$$4\pi R_i^2N_pdR_i=\frac{4\pi R_i^2}{d\mathrm{\Delta }V/dR_i}N_pd\mathrm{\Delta }V.$$
(2)
The quantity $`d\mathrm{\Delta }V/dR_i`$ has been extracted to show $`d\mathrm{\Delta }V`$, which is proportional to the differential frequency ($`d\nu `$; eq. (1)). Eq. (2) indicates that the absorption cross section at frequency $`\nu `$ is proportional to:

$$\frac{4\pi R_i^2}{d\mathrm{\Delta }V/dR_i}N_pe^{\frac{V_0(R_i)}{k_BT}},$$
(3)
where $`k_B`$ is Boltzmann’s constant, $`V_0(R_i)`$ is the ground state energy shift by a perturber at distance $`R_i`$, and the density of perturbers has been corrected by the corresponding thermal Boltzmann factor. Inverting eq. (1) to obtain $`R_i(\nu )`$ at a given absorption frequency yields an expression for the absorption cross section in the line wings at that frequency. Each solution for $`R_i`$ at a given $`\nu `$ is called a Condon point, where the subscript $`i`$ is now revealed to be the Condon point index. In principle, since $`\mathrm{\Delta }V`$ can be multiple-valued, at a given frequency there can be multiple Condon points. Expressions like eq. (3) must be summed for all such points to obtain the total cross section.

Expression (3) is not only missing the dipole transition moment (oscillator strength $`f_{mn}`$), but is not normalized by the dipole (Thomas-Reiche-Kuhn) sum rule to yield the correct cross section $`\sigma (\nu )`$. Furthermore, it must be corrected at small detunings (large $`R_i`$). The correct expression for the absorption cross section in the wings of a spectral line in the quasi-static approximation is given by the Unified Franck-Condon formula (Szudy and Baylis 1975,1996):

$$\sigma (\nu )=\underset{i}{}\sigma _i(\nu )+\underset{i,j=i+1}{}\sigma _{ij}(\nu ),$$
(4)
where

$$\sigma _i(\nu )=A\frac{R_i^2}{|\mathrm{\Delta }V^{}(R_i)|}\left(\pi z_i\right)^{1/2}\left(L(z_i)+\frac{M(z_i)}{z_i}\right)e^{\frac{V_0(R_i)}{k_BT}}$$
(5)
and

$$\sigma _{ij}(\nu )=2A\frac{R_i^2}{|\mathrm{\Delta }V^{}(R_i)\mathrm{\Delta }V^{}(R_j)|^{1/2}}\left(\pi Z_{ij}\right)^{1/2}\left(L(Z_{ij})\frac{M(Z_{ij})}{Z_{ij}}\right)e^{\frac{V_0(R_i)}{k_BT}}.$$
(6)
$`A`$ contains the transition oscillator strength and perturber number density and is given by:

$$A=12\pi N_p\frac{\pi e^2}{m_ec^2}f_{mn}.$$
(7)
The arguments $`z_i`$ and $`Z_{ij}`$ are given by the expressions:

$$z_i=\frac{1}{2}\left(\frac{\mu }{k_BT}\right)^{1/3}\left(\frac{\mathrm{\Delta }V^{}(R_i)}{\mathrm{}}\right)^2\left|\frac{\mathrm{\Delta }V^{\prime \prime }(R_i)}{\mathrm{}}\right|^{4/3}$$
(8)
and

$$Z_{ij}=\frac{1}{2}\left(\frac{\mu }{k_BT}\right)^{1/3}\left(\frac{\mathrm{\Delta }V^{}(R_i)\mathrm{\Delta }V^{}(R_j)}{\mathrm{}^2}\right)\left|\frac{\mathrm{\Delta }V^{\prime \prime }(R_i)\mathrm{\Delta }V^{\prime \prime }(R_j)}{\mathrm{}^2}\right|^{2/3},$$
(9)
where

$$\mathrm{\Delta }V^{}(R_i)=\frac{\mathrm{\Delta }V}{R_i}$$
(10)
and $`\mu `$ is the reduced mass of the interacting pair. $`L(z_i)`$ and $`M(z_i)`$ are functions defined through the integrals:

$$L(z_i)=_0^{\mathrm{}}𝑑\zeta \zeta ^2Ai^2(z_i\zeta )e^{\zeta ^3}$$
(11)
and

$$M(z_i)=_0^{\mathrm{}}𝑑\zeta \zeta ^3Ai^2(z_i\zeta )e^{\zeta ^3},$$
(12)
where $`Ai(x)`$ is the Airy function and $`Ai^{}(x)`$ is its derivative. The asymptotic behavior of $`L(z_i)`$ and $`M(z_i)`$ for positive $`z_i`$ is:

$$L(z_i)=\left(36\pi z_i\right)^{1/2}$$
(13)
and

$$M(z_i)=\left(\frac{z_i}{36\pi }\right)^{1/2}.$$
(14)
Eq. (5) is the expression that corresponds to eq. (3), the sum is over the Condon points ($`i`$ and $`j`$), and $`\sigma _{ij}(\nu )`$ is the interference term between the Condon points. The interference term is important only when $`R_iR_j`$, at which point $`\mathrm{\Delta }V^{}(R_i)`$ is formally zero. This singularity (seen in eqs. (3) and (5)) yields what is referred to as a “rainbow” or “satellite” feature (Beuc and Horvatic 1992) in the wings of an absorption spectrum. We expect such features in the Na I (5890 Å) and K I (7700 Å) absorption cross sections (§4 and §5). Note that at large $`z_i`$, the expression $`\left(\pi z_i\right)^{1/2}\left(L(z_i)+\frac{M(z_i)}{z_i}\right)`$ in eq. (5) becomes a constant and eq. (5) reduces to eq. (3), modulo constants and the oscillator strength. Moreover, at large $`Z_{ij}`$, the expression $`\left(\pi Z_{ij}\right)^{1/2}\left(L(Z_{ij})\frac{M(Z_{ij})}{Z_{ij}}\right)`$ in eq. (6) goes to zero.

The exponential terms in eqs. (5) and (6) provide natural cutoffs at positive (red) and negative (blue) detunings and depend only upon the ground-state potential shifts and temperature. Such cutoffs have been shown to be necessary for reasonable spectral fits to T dwarf spectra longward of $``$1.0 $`\mu `$m (Burrows et al. 2002) and are naturally provided by the quasi-static UFC theory.

The line cores (not the subject of this paper) are determined by distant encounters and are handled by assuming a van der Waals interaction potential with an adiabatic impact theory (Weisskopf 1933; Ch’en and Takeo 1957; Dimitrijević and Peach 1990), as in BMS and Nefedov, Sinel’shchikov, and Usachev (1999). The result is a pressure-broadened core of Lorentzian shape that applies out to only tens of wavenumber (cm-1) detunings. The major contribution to the alkali line profile shape is determined by the statistical UFC theory, which applies beyond the core to thousands of wavenumber (and thousands of Angstrom) detunings. The procedure for obtaining the wing absorption profiles is first to calculated the potential curves for the relevant levels of Na and K (§3 and §4), then to calculate the Condon points, and finally to employ eq. (4). We do this for Na+H2, Na+He, K+H2, and K+He interaction pairs, and for the states that contribute to both the red and the blue wings of the 5890 Å and 7700 Å resonance features. For perturbations by H2, the potential interactions are functions of angle (orientation), introducing yet another degree of complexity.

## 3 Technique for Calculating the Interaction Potentials and Line Shifts

The Na-D doublet of neutral sodium involves the $`3s^2S_{1/2}3p^2P_{3/2,1/2}`$ transitions, where the D2 transition (5890 Å) is from the $`1/2`$ state to the $`3/2`$ state with a multiplicity of 4 ($`=2j+1`$) and the D1 transition (5896 Å) is from $`1/2`$ to $`1/2`$ with a multiplicity of 2. Hence, the ratio of the D2 and D1 line strengths is 2. The analogous transitions for potassium are $`4s^2S_{1/2}4p^2P_{3/2,1/2}`$ at 7665 Å and 7699 Å, respectively, with a similar strength ratio. Since the upper states of the D2 transitions are $`j=3/2`$ states, they include $`p_x`$ and $`p_y`$ states that in molecular orbital theory can form bonding states with H2 (oriented perpendicularly). The $`p_z`$ state associated with the $`{}_{}{}^{2}P_{1/2}^{}`$ state is anti-bonding to H2. Figure 6 depicts the orientations of the interacting H2 and Na pair for C2v symmetry (H$``$H line perpendicular to the Na-H2 line). The potassium $`4p`$ states can be depicted in a perfectly analogous way. The $`{}_{}{}^{2}A_{1}^{}`$ states depicted in Fig. 6 are anti-bonding and the interaction potential is purely repulsive at all distances $`R_i`$. The $`{}_{}{}^{2}B_{1,2}^{}`$ states are attractive at large distances and repulsive at small distances 111Note that in order to connect to both communities, we use both the spectroscopist’s ($`A`$ and $`B`$) and the physicist’s ($`P_{1/2,3/2}`$) nomenclature for the states.. The ground states of both the Na+H2 and the K+H2 systems are repulsive for all distances $`R_i`$. The upshot is that the energy difference between the 2B state of the D2 transition and the ground state decreases with decreasing $`R_i`$ until $`R_i`$ is quite small. Hence, the D2 line, though it is blueward of the D1 line for atoms in isolation, provides the red wing of the alkali line. Correspondingly, the D1 line provides the blue wing.

A schematic of the excited-state minus ground-state energy shifts of the Na-D and K I resonance lines is given in Fig. 6. Figure 6 depicts the mapping between energy shift, line core, red wing, and blue wing. Note that $`\mathrm{\Delta }V`$ on the blue wing has a maximum (extremum), which results in a rainbow satellite feature. The absence of a corresponding extremum on the red wing means that there is no corresponding satellite feature redward of the line core.

To calculate the interaction potentials, we employ the multi-configurational self-consistent-field (MCSCF) variant of the quantum chemical code GAMESS (Schmidt et al. 1993). The Hartree-Fock (HF) method is a self-consistent field approach that employs iteration to find a consistent solution to the multi-electron, multi-nucleus wavefunction (itself a Slater determinant), but assumes that the potential is a smoothed average and does not explicitly include electron-electron interactions. MCSCF is a post-HF method that includes electron-electron interactions and exchange forces, the so-called configuration interaction (CI), but assumes that the CI term is small. It allows variation of not only the mixing coefficients of the various HF configurations, but also of the coefficients of the basis functions in the constituent molecular orbitals. In this way, the HF orbitals from which the total wavefunction is constructed are optimized simultaneously. The MCSCF method requires care in the selection of the basis set, but generally achieves a good solution for problems having low-lying excited states, as in the problem at hand. We use the Complete Active Space Self-Consistent Field (CASSCF) version of MCSCF and the 6-31G∗∗ split-valence basis set that incorporates Gaussian Type Orbitals (GTO) (Frisch et al. 1998). In practice, we first obtain a solution with the HF method, and then use this solution as a first guess in the full MCSCF/CASSCF calculation. This two-step approach has often proven necessary to achieve converged solutions.

We use the above methods and the GAMESS code to solve the many-electron Schrödinger equation for the wavefunctions and the energy levels of the various electronic states as a function of distance $`R_i`$ for the composite system of alkali metal atom and either H2 or helium. For H2, this is done for a variety of orientations (from C2v \[$`\theta =90^{}`$\] to C \[$`\theta =0^{}`$\] (parallel) symmetries) and we keep the H–H bond length constant (Botschwina et al. 1981). We have included electron correlation effects successfully for He perturbations and for H2 perturbations germane to the red wings of both the Na and K lines. However, when calculating the blue wings for orientation angles other than 0 or 90 we had difficulty obtaining converged solutions with the electron correlation effects (CI) included. Hence, for the blue wings of the Na/K+H2 interaction pairs we have settled for the pure Hartree-Fock solutions. This approximation is good for large $`R_i`$, but begins to break down at small $`R_i`$ and large detunings.

For all the relevant $`{}_{}{}^{2}A`$ and $`{}_{}{}^{2}B`$ states we calculate the interaction potentials as a function of both distance and (for H2) angle of orientation. The calculations are done at hundreds of distances and (for H2) for five angles ($`\theta =0^{},20^{},45^{},70^{},`$ and $`90^{}`$). Then, we calculate the energy shifts ($`\mathrm{\Delta }V`$) of the resonance transitions of Na and K. Though we reproduce the wavelengths of the isolated D2 and D1 lines of Na to within $``$14% (0.3 eV) and for K to within $``$21% (0.34 eV), we shift our calculated energy shifts by a constant amount that ensures that the line shift at very large $`R_i`$ is zero. This puts the line core exactly at the measured wavelength by construction. Furthermore, we do not calculate the dipole transition matrix elements as a function of $`R_i`$, but use the published oscillator strengths of the unperturbed Na and K resonance doublets in the calculation of the line profiles (eq. 7). These incorporate the 2:1 strength ratio due to the spin degeneracy of the final states in the unperturbed D2 and D1 transitions.

## 4 Potential Curves for the Na+H2, K+H2, and Na/K+He Systems

Figure 6 depicts the energy levels as a function of separation of the various states (relative to the ground state at infinite separation), for both Na and K perturbed by both helium and H2. For H2, the run of energy shift with $`R_i`$ for the four states is given for various orientation angles. As Fig. 6 shows, the ground state interaction is always repulsive and is most repulsive for H2 for colinear orientations (C symmetry). This is the $`V_0(R_i)`$ term to be used in eqs. (5) and (6) in calculating the cutoff exponential. The $`{}_{}{}^{2}B`$ states associated with the $`p_x`$ and $`p_y`$ atomic orbitals of sodium or potassium show a slight attraction for large angles (close to perpendicular) and large distances, but then become universally repulsive as the hard core is approached (Botschwina et al. 1981). Figure 6 portrays the resulting transition energies for the D lines of both Na and K as a function of separation and angle. As is clear from the figure, the energy differences for the transitions for which the excited states are $`{}_{}{}^{2}B`$ states decrease with decreasing separation and, hence, are associated with the red wings of the spectral lines. Similarly, since the ground-state to $`{}_{}{}^{2}A_{1}^{}`$-state transitions have the opposite behavior, they are associated with the blue wings. Furthermore, the peaks seen at the bottom of Fig. 6 at separations of $`R_i2.42.5`$Å and $`R_i3.13.2`$Å, respectively, for Na and K perturbed by H2 imply that their blue wings have satellites at wavelengths of $``$0.49-0.5 $`\mu `$m and $``$0.66–0.68 $`\mu `$m . The K satellite might be near the Li feature at 6708 Å. Fig. 6 schematically shows how the panels of Fig. 6 merge to make the corresponding absorption profiles.

It is useful to compare these potential calculations with previous related work using different methods. Rossi and Pascale (1985, RP) have done $`l`$-dependent pseudopotential calculations for the Na+H2 and K+H2 systems and Botschwina et al. (1981) have done RHF-SCF and PNO-CEPA calculations for the Na+H2 system. The PNO-CEPA potential curves of Botschwina et al. are similar to those of RP. Neither did the calculations for the arbitrary H2 orientations and small separations necessary to perform the line profile calculations that motivated this paper. Rather, they concentrated solely on separations greater than 2.0 Å and on the C and C2v symmetries. Nevertheless, their potential curves for these symmetries are similar to those found in Fig. 6. For C2v symmetry and the Na+H2 system, RP derive that the potential minimum for the $`{}_{}{}^{2}B_{1}^{}`$ state is at 2.4 Å and has a depth of 0.12 eV; we obtain 2.45 Å and 0.1 eV. For the $`{}_{}{}^{2}B_{2}^{}`$ state, RP derive corresponding numbers of 2.12 Å and 0.35 eV; we obtain 2.27 Å and 0.2 eV. For the K+H2 system, RP derive for the $`{}_{}{}^{2}B_{2}^{}`$ state a minimum at 3.0 Å with a depth of 0.49 eV; we obtain 3.0 Å and 0.5 eV. RP estimate that the change in the potential of the $`{}_{}{}^{2}A_{1}^{}`$ excited state of the Na+H2 system with C2v symmetry in going from 4.0 Bohr radii (2.12 Å) to 8.0 Bohr radii (4.23 Å) is 0.55 eV; we obtain 0.7 eV. Furthermore, Kleinekathöfer et al. (1996), who use a surface integral method to derive the ground-state potential curve of the Na+He system, derive V($`R_i`$ = 3.2 Å) = 0.10 eV; we obtain 0.08 eV at the same distance. All in all, our results compare favorably with those in previous studies, with the Na+H2/$`{}_{}{}^{2}B_{1}^{}`$ and K+H2/$`{}_{}{}^{2}B_{2}^{}`$ potential curves corresponding most closely. However, the correspondences are not perfect and we sometimes see discrepancies of as much as $``$0.1 eV at a given separation. We suspect that these discrepancies are traceable either to our choice of basis sets or to the MCSCF assumption that the CI term is small.

Since the second derivative of the potential ($`V^{\prime \prime }`$) is needed when using eqs. (5) and (6) to obtain the satellite profile shapes, small errors in the potential calculations when $`V^{}=0`$ can be amplified at the rainbow positions. Ch’en and Wilson (1961) measured the “violet band separation” for the satellite features on the blue wings of the Na+He and the K+He systems to be 1161$`\pm `$15 cm-1 and 838$`\pm `$5 cm-1, respectively, which translate into satellite positions at 0.55 $`\mu `$m and 0.72 $`\mu `$m . These contrast with the satellite positions we derive of 0.515 $`\mu `$m and 0.69 $`\mu `$m , respectively, for perturbation by helium. Not unexpectedly, the centers of the satellite features are the most poorly determined quantities in our theory. We now turn to the calculation of the opacity profiles in the wings of the D2 and D1 lines of Na and K due to perturbations by H2 and helium using the potential curves represented in Figs. 6 and 6.

## 5 The Absorption Cross Sections on the Red and Blue Wings of the Na (5890Å) and K (7700Å) Resonance Lines

Using the UFC/quasi-static eqs. (5) and (6), we obtain the absorption spectra as a function of photon wavelength for the D2 and D1 lines of Na and K centered near 5890 Å (0.589 $`\mu `$m ) and 7700 Å(0.77 $`\mu `$m ), respectively. These are depicted in Fig. 6 at a temperature of 1000 K and a total pressure of one atmosphere. The contributions at various orientation angles of H2, the angle-integrated spectra (solid), and the results for the Na+He and K+He systems are shown. The Lorentzian cores are not included on this plot. Since the ground-state potential shift increases quickly at small separations (Fig. 6), and the associated exponential terms in eqs. (5) and (6) introduce a large cutoff, Condon points at such small separations make only a marginal contribution to the overall opacities. As a consequence, though there can be more than two Condon points (Fig. 6), we include at most the outer two in our theory. For the same reason, we ignore possible satellites on the very far red wings.

Figure 6 encapsulates the central results of our study. The cutoff on the red wing of the potassium feature due to the exponential term (eq. 5) containing the ground-state interaction potential of the K+H2 system is seen to be situated between 0.95 $`\mu `$m and 1.0 $`\mu `$m , close to the 0.98 $`\mu `$m used in Burrows et al. (2002). The cutoff of the red wing of the Na feature is near 0.8 $`\mu `$m , but is a bit more gradual. The corresponding cutoffs for the alkali-He systems are more abrupt, but of less importance due to the lower abundance of helium. It should be noted that for the temperatures in cool atmospheres the precise form and position of this cutoff is a sensitive function of the ground-state potential. Errors in the potential of $`\stackrel{<}{}`$0.1 eV that may arise at small $`R_i`$ for the alkali-H2 systems translate directly into shifts in the position of the turnover and in the magnitude of the cross section at the turnover. Such an effect is demonstrated by the difference between the cross sections with and without electron-correlation/configuration-interaction terms (depicted below in Fig. 6). Figure 6 also makes clear that the red and blue wings are asymmetrical, as one would expect from the discussion in §4 and Fig. 6.

The satellite features identified on the blue wings of both the Na and K features by the asterisks on Fig. 6 have been alluded to before, but we note once again that our calculations are too imprecise to derive the detailed shape of the actual rainbow structures (which for the alkali-H2 systems is a composite of the corresponding curves for the various angles). Nevertheless, the existence and approximate wavelengths of the satellites remain clear.

Figure 6 portrays the temperature and pressure dependences of the line profiles. Three curves (at 800, 1000, and 1500 K) show the temperature dependence at one atmosphere pressure and one curve at 1000 K and 10 atmospheres shows, by comparison with the other 1000 K curve, the pressure dependence. Due to the $`N_p`$ dependence in eq. (7), higher pressures lead to higher wing strengths. The Thomas-Reiche-Kuhn sum rule is still satisfied by a corresponding diminution in the strength in the Lorentzian core. Due to the exponential term, increasing the temperature raises the cross section in the far wings. However, due to the constraint of constant pressure, increasing the temperature slightly lowers the cross sections in the near wings.

It should be emphasized that, unlike the theory of BMS, and the line profile formalism used in Burrows et al. (2002), there are no free parameters in our theory. The cutoffs and shapes arise naturally because the physics of the line-broadening process has been treated. Limitations on the accuracy of the current work are in the accuracy of the quantum chemistry calculations of the potential curves and in the use of the semiclassical formalism to treat the broadening. The earlier theory of BMS was constructed to fill the theoretical and experimental void concerning the far wing profiles of the alkali metal lines perturbed by molecular hydrogen. Burrows et al. (2002) followed with a simpler theory that used a Gaussian-truncated Lorentzian profile. Previously in astrophysics the character of the far wings of the K I and Na D resonance features perturbed by H2 had not been needed.

Given our new calculations, it is useful to compare them with the stopgap theories of BMS and Burrows et al. (2002). Figure 6 is a comparison of the BMS (red) and Burrows et al. (2002) (green) potassium cross sections with the corresponding cross sections using the theory of this paper (blue). For the BMS prescription, we arbitrarily set the associated $`q`$ parameter equal to 0.5. A temperature of 1000 K and a pressure of one atmosphere were assumed. The BMS and the Burrows et al. (2002) curves include the Lorentzian core, while the newly-calculated quasi-static curve does not. Also shown are the new cross sections calculated without the configuration interaction (CI) (blue, dashed). At the inner red wing, the theory of Burrows et al. (2002), in which the pressure-broadened Lorentzian was truncated with a Gaussian near 0.98 $`\mu `$m , is within a factor of three of the more detailed calculation, but clearly deviates from it. However, on the red wing the simple theory of BMS is actually surprisingly close to the new result. On the far red wing, the cross sections of the new theory plummet more precipitously than either BMS or Burrows et al. (2002). This emphasizes even more strongly the existence of the steep cutoff inferred by Burrows et al. (2002) to be necessary in order to reproduce the $`Z`$ ($`1.0`$ $`\mu `$m ) and $`J`$ ($`1.2`$ $`\mu `$m ) band fluxes observed for the T dwarfs (Burgasser et al. 2002). Note, however, that the new cross sections on the far blue wing are as much as a factor of four higher than the corresponding cross sections of the old theories. The same result obtains for the blue wing of the Na-D line. In sum, though the BMS theory was good in the near red wings, in places the contrast between both the BMS and the Burrows et al. (2002) formulations and our parameterless UFC calculations is large. As a result, we recommend that in the future the new theory be used for detailed T and L dwarf spectral calculations.

## 6 Conclusions

Using the quantum chemical code GAMESS, we have calculated as a function of distance the interaction potentials of the excited and ground states of sodium and potassium perturbed by H2 and helium. For H2, we also calculated the dependence of the interaction potentials on orientation angle. Using these potentials and the Unified Franck-Condon formalism, we derived the absorption line profiles in the wings of the Na-D doublet and the K I resonance doublet at detunings of thousands of Angstroms. These lines and their wings are of central importance in the atmospheres and spectra of brown dwarfs, T dwarfs, L dwarfs, and hot Jupiters. Our theory has no free parameters and naturally accounts for the cutoffs in the far red wings inferred to be present to explain the high $`J`$ and $`Z`$ band fluxes observed in T dwarfs (Burrows et al. 2002). However, the new cross sections differ from those used previously by having sharper cutoffs, stronger blue wings, and satellite features. Furthermore, the red/blue asymmetry is pronounced.

From our calculations using the full configuration interaction at the subset of useful H2 orientation angles for which we were able to obtain convergence for the $`{}_{}{}^{2}A_{1}^{}`$ states (of relevance for the calculation of the blue wings) and from a comparison of the Hartree-Fock and CI results for perturbations by helium, we have determined that the HF and CI numbers on the blue wings differ little, though the associated potential differences on the red wings can differ by $``$0.1 eV. This gives us some confidence that, despite the fact that we use the HF approach on the blue wings for perturbations by H2, our blue wing results are satisfactory. Fortunately, for all the red wing results, we were able to obtain full convergence using the MCSCF/CI formalism at all H2 angles.

The need for line wing strengths at such large detunings has few precedents in astrophysics. The next step is to incorporate the new alkali line opacities into a spectral synthesis code and to identify the differences between the fits using the old and the new algorithms. The line cores will be unchanged, but the blue wings and the spectra from 0.8 $`\mu `$m to 0.9 $`\mu `$m on the red wings of the K I line will be slightly altered, as else being equal. The major virtue of our new algorithm is that it is parameterless. As a result, we reduce by one (or even a few) the number of free parameters needed in fits to the increasing number of optical and near-IR spectra being measured for T dwarfs (in particular), thereby making the extraction of physical quantities such as effective temperature and gravity more robust and reliable.

We are happy to thank Ludwik Adamowicz, Jonathan Lunine, Christopher Sharp, David Sudarsky, and Bill Hubbard for help, advice, and useful discussions. This work was supported in part by NASA grants NAG5-10760, NAG5-10629, NAG5-7499, and NAG5-7073 and with the help of JPL contract # 1241086 under the NASA Navigator program. Tables of our potential curves are available from the first author upon request.


text:
# Axions: Past, Present, and Future

## 1 The Old and New U(1) Problems

Consider quantum chromodynamics (QCD) with color group SU(3) and three flavors of massless quarks. The lagrangian is invariant under an additional global flavor group, U$`(3)_\mathrm{L}\times `$U$`(3)_\mathrm{R}`$. Although we will use it anyway, this “left–right” terminology is actually somewhat misleading; in four dimensions, a right-handed fermion field is the hermitian conjugate of a left-handed fermion field, and so we could (and will) adopt the convention that all fundamental fermion fields are left-handed (with, of course, right-handed hermitian conjugates). In this case, there would be three left-handed fields transforming as $`\mathrm{𝟑}`$’s of SU$`(3)_{\mathrm{color}}`$, and three left-handed fields transforming as $`\overline{\mathrm{𝟑}}`$’s of SU$`(3)_{\mathrm{color}}`$; thus the flavor group would be better named U$`(3)_\mathrm{𝟑}\times `$U$`(3)_{\overline{\mathrm{𝟑}}}`$. \[We note in passing that if quarks were in a real or pseudo-real representation of the color group, such as $`\mathrm{𝟑}`$ of O$`(3)_{\mathrm{color}}`$ or $`\mathrm{𝟐}`$ of SU$`(2)_{\mathrm{color}}`$, the lagrangian flavor group would be U(6), since in this case there is no distinction between left-handed quarks and left-handed antiquarks.\] For now we ignore the effects of the quantum anomaly, which will of course play a key role later.

QCD exhibits dynamical breaking of the flavor symmetry; this can be understood as the formation of a quark condensate,

$$0|\chi _{\alpha ia}\stackrel{~}{\chi }^\beta {}_{\overline{ȷ}b}{}^{}|0=\frac{1}{6}\mathrm{\Lambda }^3\delta _\alpha {}_{}{}^{\beta }\delta _{i\overline{ȷ}}^{}\epsilon _{ab}.$$
(1)
Here $`\chi _{\alpha ia}`$ is a left-handed fermion field transforming as a $`\mathrm{𝟑}`$ of SU$`(3)_{\mathrm{color}}`$ with color index $`\alpha =1,2,3`$, “left” flavor index $`i=1,2,3`$, and left-handed spinor index $`a=1,2`$; $`\stackrel{~}{\chi }^\beta _{\overline{ȷ}b}`$ is a left-handed fermion field transforming as a $`\overline{\mathrm{𝟑}}`$ of SU$`(3)_{\mathrm{color}}`$ with color index $`\beta =1,2,3`$, “right” flavor index $`\overline{ȷ}=1,2,3`$, and left-handed spinor index $`b=1,2`$; $`\epsilon _{ab}`$ is the antisymmetric invariant symbol of SU(2); and $`\mathrm{\Lambda }`$ is a parameter with dimensions of mass. Under a U$`(3)_\mathrm{L}`$ transformation, $`\stackrel{~}{\chi }`$ is unchanged, and $`\chi _{\alpha ia}L_i{}_{}{}^{j}\chi _{\alpha ja}^{}`$, where $`L`$ is a unitary matrix; under a U$`(3)_\mathrm{R}`$ transformation, $`\chi `$ is unchanged, and $`\stackrel{~}{\chi }^\beta {}_{\overline{ı}b}{}^{}(R^{})_{\overline{ı}}{}_{}{}^{\overline{ȷ}}\stackrel{~}{\chi }_{}^{\beta }_{\overline{ȷ}b}`$, where $`R^{}`$ is an independent unitary matrix (the complex conjugation is a notational convention). The condensate is unchanged only by transformations in the “vector” subgroup U(3$`)_\mathrm{V}`$ specified by $`R=L`$. Thus, U$`(3)_\mathrm{L}\times `$U$`(3)_\mathrm{R}`$ is spontaneously broken down to U(3$`)_\mathrm{V}`$. \[I leave it as an exercise for the reader to show that the unbroken subgroup would be O$`(6)`$ if quarks were $`\mathrm{𝟑}`$’s of O$`(3)_{\mathrm{color}}`$, and Sp(6) if quarks were $`\mathrm{𝟐}`$’s of SU$`(2)_{\mathrm{color}}`$.\] The nine broken generators lead us to expect nine Goldstone bosons; these can be thought of as long wavelength excitations of the condensate,

$$0|\chi _{\alpha ia}(x)\stackrel{~}{\chi }^\beta {}_{b\overline{ȷ}}{}^{}(x)|0=\frac{1}{6}\mathrm{\Lambda }^3\delta _\alpha {}_{}{}^{\beta }\epsilon _{ab}^{}U(x)_{i\overline{ȷ}},$$
(2)
where $`U`$ is a unitary matrix field. Under a U$`(3)_\mathrm{L}\times `$U$`(3)_\mathrm{R}`$ transformation,

$$U(x)_{i\overline{ȷ}}L_i{}_{}{}^{k}(R^{})_{\overline{ȷ}}^{}{}_{}{}^{\overline{\mathrm{}}}U(x)_{k\overline{\mathrm{}}},$$
(3)
or equivalently $`ULUR^{}`$. We can write an effective lagrangian for $`U`$; it must of course be U$`(3)_\mathrm{L}\times `$U$`(3)_\mathrm{R}`$ invariant. There are no allowed terms with no derivatives, and two with two derivatives:

$$L=\frac{1}{4}f_\pi ^2\mathrm{Tr}^\mu U^{}_\mu U\frac{1}{18}(f_1^2\frac{3}{2}f_\pi ^2)^\mu (detU^{})_\mu (detU)+\mathrm{},$$
(4)
where $`f_\pi `$ and $`f_1`$ are parameters with dimensions of mass. If we write

$$U(x)=\mathrm{exp}\left[i_{a=1}^8\lambda ^a\pi ^a(x)/f_\pi i\pi ^9(x)/f_1\right],$$
(5)
where the Gell-Mann $`\lambda `$ matrices are hermitian and normalized via $`\mathrm{Tr}\lambda ^a\lambda ^b=2\delta ^{ab}`$, then the hermitian Goldstone fields $`\pi ^a(x)`$ ($`a=1,\mathrm{},9`$) have canonical kinetic terms, $`L=\frac{1}{2}^\mu \pi ^a_\mu \pi ^a`$. The parameter $`f_\pi `$ is called the pion decay constant, and it has the measured value of $`92.4`$MeV. (It is measured via the $`\pi \mu \nu `$ decay rate; see Ref. \[\] for details.)

In the real world, the three light quarks have small masses:

$$L_{\mathrm{mass}}=M^{i\overline{ȷ}}\epsilon ^{ab}\chi _{\alpha ia}\stackrel{~}{\chi }^\alpha {}_{b\overline{ȷ}}{}^{}+\mathrm{h}.\mathrm{c}.,$$
(6)
where $`M`$ is in general a complex matrix with no particular symmetries; $`M^{}M`$ has eigenvalues $`m_u^2`$, $`m_d^2`$, and $`m_s^2`$, where $`m_{u,d,s}`$ are the three light-quark masses. $`M`$ can be made diagonal with positive real entries $`m_u`$, $`m_d`$, and $`m_s`$ via a U$`(3)_\mathrm{L}\times `$U$`(3)_\mathrm{R}`$ transformation that leaves the rest of the lagrangian unchanged. In terms of the effective lagrangian, it is easy to see from eqs. (1) and (6) that

$$L_{\mathrm{mass}}=\mathrm{\Lambda }^3\mathrm{Tr}(MU+\mathrm{h}.\mathrm{c}.).$$
(7)
Expanding in the Goldstone fields, we get

$$L_{\mathrm{mass}}=\frac{1}{4}\mathrm{\Lambda }^3\mathrm{Tr}[M(\lambda ^a\lambda ^b+\lambda ^b\lambda ^a)]\pi ^a\pi ^b/f_\pi ^2,$$
(8)
where $`\lambda ^9(f_\pi /f_1)I`$. The key conclusion we want to obtain is reached most easily by taking the exact isospin limit $`m_u=m_dm`$; then the eigenvalues of the mass-squared matrix in eq. (8) are $`m_{1,2,3}^2=4\mathrm{\Lambda }^3m/f_\pi ^2`$ (these are the three pions) and $`m_{4,5,6,7}^2=2\mathrm{\Lambda }^3(m+m_s)/f_\pi ^2`$ (these are the four kaons). The $`\pi ^8`$ and $`\pi ^9`$ fields have a mass-squared matrix that must be diagonalized; taking $`mm_s`$, the eigenvalues are $`2\mathrm{\Lambda }^3m_s(\frac{4}{3}f_\pi ^2+f_1^2)`$ and $`12\mathrm{\Lambda }^3m/(f_\pi ^2+\frac{4}{3}f_1^2)`$. This implies that there is a fourth pion-like particle with mass less than $`\sqrt{3}m_\pi `$, where $`m_\pi =m_{1,2,3}=135`$MeV is the neutral pion mass. There is, in nature, no such particle. This discrepancy between theory and experiment is the “U(1) problem” or (nowadays) the “old U(1) problem”.

The old U(1) problem is solved by the axial anomaly: while the lagrangian is invariant under U$`(3)_\mathrm{L}\times `$U$`(3)_\mathrm{R}`$, the measure over the fermion fields in the functional integral is not; the measure transforms nontrivially under the “axial” U$`(1)_\mathrm{A}`$ subgroup $`L=R^{}=e^{i\alpha }I`$. Furthermore, the presence of instanton solutions of the euclidean field equations of QCD allow us to see explicitly the physical effects of the anomaly in semiclassical calculations. Thus, when we formulate an effective lagrangian for the low-energy fields, we have no reason to expect it to be invariant under this subgroup. This allows us to add terms involving $`detU`$ (and its derivatives) to the lagrangian:

$$L_{\mathrm{anom}}=f_\pi ^4_{n=1}^{\mathrm{}}c_n(detU)^n+\mathrm{h}.\mathrm{c}.+\mathrm{}.$$
(9)
In the absence of quark masses, the functional integral of QCD is CP invariant. In the effective theory, CP exchanges $`U`$ and $`U^{}`$; CP invariance then implies that the $`c_n`$ coefficients in eq. (9) must all be real. The main effect of $`L_{\mathrm{anom}}`$ is to give the $`\pi ^9`$ field a large mass, $`m_9^2=18Cf_\pi ^4/f_1^2`$, where $`C=_{n=1}^{\mathrm{}}c_nn^2`$; the $`\pi ^9`$ field is then essentially removed from the effective theory, and the mass of the $`\pi ^8`$ field (physically, the $`\eta `$ meson) is given by $`\frac{8}{3}\mathrm{\Lambda }^3m_s/f_\pi ^2`$.

A shadow of the $`\pi ^9`$ field remains, however; in the presence of $`L_{\mathrm{anom}}`$, we cannot remove the overall phase of the quark mass matrix $`M`$ via a U$`(3)_\mathrm{L}\times `$U$`(3)_\mathrm{R}`$ transformation. If $`detM=m_um_dm_se^{i\theta }`$, and if we then make a U$`(3)_\mathrm{L}\times `$U$`(3)_\mathrm{R}`$ transformation on $`U`$ that brings $`M`$ to the form $`\mathrm{diag}(m_u,m_d,m_s)`$, the $`c_n`$’s change: $`c_nc_ne^{in\theta }`$. This would result in CP violating effects in hadron interactions; the most visible of these would be an electric dipole moment for the neutron. Experimental limits on the electric dipole moment of the neutron imply $`\theta <10^9`$. This, then, is the “strong CP problem”, or the “new U(1) problem”: why is $`\theta `$ so small?

One possibility is that $`detM=0`$ because $`m_u=0`$; then there is no phase to remove from the quark mass matrix. It turns out that this possibility cannot be ruled out solely by experimental evidence (since higher-order corrections in the $`d`$ and $`s`$ quark masses can mimic a nonzero $`u`$ quark mass), but it does require some severe theoretical contortions (such as the invalidity of the large-$`N_{\mathrm{color}}`$ expansion to all orders); see Ref. \[\] for a discussion.

## 2 The Peccei–Quinn Solution

Let us choose to work in a field basis in which the $`c_n`$’s are real and the quark-mass matrix $`M`$ is diagonal, but contains the unremovable phase $`\theta `$: $`M=\mathrm{diag}(m_ue^{i\theta },m_d,m_s)`$. Now the potential energy is minimized at nonzero values of $`\pi ^3`$, $`\pi ^8`$, and $`\pi ^9`$. For small $`\theta `$, this minimum energy works out to be

$$V_{\mathrm{min}}(\theta )=\frac{m_um_dm_s}{m_um_d+m_dm_s+m_sm_u}\mathrm{\Lambda }^3\theta ^2.$$
(10)
This suggests a possible resolution of the strong-CP problem. If $`\theta `$ was some sort of dynamical variable, then it would want to relax to zero, in order to minimize the energy. This is the essential physics of the solution proposed by Peccei and Quinn. In their model, $`\theta `$ is effectively replaced by $`a(x)/f_a`$, where $`a(x)`$ is a new scalar field, the axion, and $`f_a`$ is a new parameter, the axion decay constant. The observed value of $`\theta `$ is now zero, since that is what minimizes the energy.

Of course, this new field provides us with a new particle as well. Its mass is given by

$`m_a^2`$ $`=`$ $`{\displaystyle \frac{1}{f_a^2}}{\displaystyle \frac{d^2V_{\mathrm{min}}(\theta )}{d\theta ^2}}|_{\theta =0}`$ (11)
$`=`$ $`{\displaystyle \frac{2m_um_dm_s}{m_um_d+m_dm_s+m_sm_u}}{\displaystyle \frac{\mathrm{\Lambda }^3}{f_a^2}}`$
$`=`$ $`\left[{\displaystyle \frac{m_um_dm_s/(m_u+m_d)}{m_um_d+m_dm_s+m_sm_u}}\right]{\displaystyle \frac{f_\pi ^2}{f_a^2}}m_\pi ^2,`$
where we have used $`m_\pi ^2=2(m_u+m_d)\mathrm{\Lambda }^3/f_\pi ^2`$ to get the last line. The factor in brackets has the numerical value $`0.225\pm 0.005`$. We see that if $`f_af_\pi `$, then the axion is very light.

Peccei and Quinn did not present their solution in this language. They were concerned with constructing a renormalizable extension of the Standard Model of quarks and leptons that did not suffer from the strong CP problem. In their model, they used two Higgs fields, instead of one: one to give mass to the up quarks, and the other to give mass to the down quarks (and, somewhat incidentally, the charged leptons). The axion then appears as the relative phase of these two fields, and it is massless (classically) if the lagrangian is invariant under an extra U(1) symmetry, now known as Peccei-Quinn symmetry, that is spontaneously broken by the Higgs vacuum expectation values. In working out the physics, Peccei and Quinn assumed that the anomaly would give a large mass to the axion; the fact that it has a small mass was noticed by Weinberg and Wilczek. In the Peccei-Quinn model, the axion decay constant is given by $`f_a=\frac{1}{6}\mathrm{sin}(2\beta )v_{\mathrm{EW}}`$, where $`v_{\mathrm{EW}}=(\sqrt{2}G_\mathrm{F})^{1/2}=246`$GeV is the electroweak scale, and $`\mathrm{tan}\beta `$ is the ratio of the two Higgs-field VEVs. The couplings of the axion to quarks, leptons, and photons are all proportional to $`1/f_a`$, and hence of typical weak-interaction strength. Thus the axion should be roughly as visible as a neutrino. It soon became apparent that there is no such particle in nature.

## 3 Invisible Axions

From our presentation in the last section, it should be clear that there is no fundamental principle that relates $`f_a`$ to $`v_{\mathrm{EW}}`$; this is simply a consequence of the specific model explored by Peccei and Quinn. Making $`f_a`$ larger would make the axion lighter, but also more weakly coupled to ordinary matter of all kinds, and hence harder to produce or detect.

The first successful theoretical proposal to detach $`f_a`$ from $`v_{\mathrm{EW}}`$ came from Kim. His model is just the Standard Model (with the usual one Higgs doublet), plus an additional sector consisting of a Higgs field that is an SU(3)$`\times `$SU(2)$`\times `$U(1) singlet, with a Yukawa coupling to a new quark field that is an SU(2) singlet and carries an essentially arbitrary hypercharge. (Thus, the electric charge of this quark is equal to its hypercharge.) The Peccei-Quinn symmetry of this model involves only these new fields. The singlet Higgs field gets a vacuum expectation value $`v`$, breaking the PQ symmetry. Since this Higgs VEV does not break any gauge symmetries, it can be arbitrarily large. If $`vv_{\mathrm{EW}}`$, then the phase of this field becomes the axion, and $`f_av`$. The new quark gets a mass given by $`v`$ times the value of the Yukawa coupling. Kim’s paper was submitted to Physical Review Letters in February 1979 and published in July 1979.

Next came the proposal of Zhitnitskiĭ; his model consisted of the Peccei-Quinn model (with its two Higgs doublets), plus an extra singlet Higgs field with (as in Kim’s model) a large VEV $`v`$. This Higgs singlet is coupled to the two doublets in such a way as to preserve the Peccei-Quinn symmetry of the original lagrangian. The large singlet VEV spontaneously breaks the PQ symmetry, and (as in Kim’s model) this results in the phase of this field becoming the axion, with $`f_av`$. Zhitnitskiĭ’s paper was submitted to Yadernaya Fizika in May 1979 and published in February 1980. The English translation (by Jonathan Rosner) appeared in the Soviet Journal of Nuclear Physics about a year later.

Neither of these papers attracted much attention at the time. I will now play amateur sociologist and attempt to guess why.

The title of Kim’s paper, “Weak-Interaction Singlet and Strong CP Invariance,” gives little clue as to its content. The abstract, however, is clear and succinct: “Strong CP invariance is automatically preserved by a spontaneously broken chiral U(1$`)_\mathrm{A}`$ symmetry. A weak-interaction singlet heavy quark $`Q`$, a new scalar meson $`\sigma ^0`$, and a very light axion are predicted. Phenomenological implications are included.” In the introduction, Kim takes a different tack, beginning with a classification scheme for theories of CP violation (his categories are hard, dynamical, and spontaneous), and remarking that “a simple theory should belong to one of these classifications for completeness of the gauge theory of weak and electromagnetic interactions.” He then presents his model. However, Kim does not emphasize the arbitrariness of $`f_a`$; when he needs a numerical value, he uses $`f_a=100`$TeV, and he takes 100 GeV for the mass of the new quark. He considers the new quark (and the heavy radial excitation of the new Higgs, $`\sigma ^0`$) to be important aspects of the model’s phenomenology. He concludes: “What are the possible experiments to prove the present scheme? Probably high-precision experiments of the axion search will do. But the easier verification of the weak-interaction singlets $`Q`$ and $`\sigma ^0`$ in $`pp`$ or $`\overline{p}p`$ machines ($`e^+e^{}`$ annihilation machine also if $`Q`$ is charged) will shed light on the whole idea of the spontaneously broken chiral U(1$`)_\mathrm{A}`$ invariance and the multiple vacuum structure of QCD.” A new quark at the electroweak scale with nonstandard electroweak interactions that could not be fit into grand unification schemes must have seemed to many at the time to be an unattractive proposition.

Zhitnitskiĭ’s paper, titled “On the possible suppression of axion-hadron interactions,” is clear, succinct, and correct throughout. His abstract reads, “A possible mechanism for strong suppression of the axion-fermion interaction is considered. Two models in which this mechanism is realized are described in detail.” The emphasis is on the axion coupling to quarks; the couplings to leptons and photons (which are also suppressed) are not mentioned. Also, the axion mass is mentioned only in the paper’s last sentence: “Here the mass of the axion is $`m_a100\delta `$keV,” where $`\delta v_{\mathrm{EW}}/v`$. In addition to the model discussed above, Zhitnitskiĭ also presents a model based on the gauge group SU(3)$`\times `$SU(2$`)_\mathrm{L}\times `$SU(2$`)_\mathrm{R}\times `$U(1); in this model, no new Higgs fields are needed other than the ones already needed to break SU(2$`)_\mathrm{R}`$. The emphasis on a left–right symmetric model probably contributed to this paper’s initially limited impact; these models were by then waning in popularity.

Next came the paper entitled “Can confinement ensure natural CP invariance of the strong interactions?” by Shifman, Vainshtein, and Zakharov (SVZ), which was submitted to Nuclear Physics B in July 1979 and published in April 1980. SVZ point out that the usual argument for the physical effect of the $`\theta `$ parameter rests on the properties of long-range instanton field configurations, and it is not immediately obvious how to treat these in the presence of color confinement. Could it be that color confinement renders $`\theta `$ unobservable after all? Their answer is no; my presentation in Sect. 1 of low-energy physics with $`\theta 0`$ essentially follows their argument. They show, in particular, that nonzero $`\theta `$ results in a nonzero amplitude for $`\eta \pi \pi `$ decay. So the effect of $`\theta `$ is really there, and the strong CP problem is really a problem. This analysis makes up the bulk of the paper.

So, what to do? SVZ begin the last section of their paper by stating that, “In fact, the problem of the $`\theta `$-term cannot be solved by QCD alone since it is intimately related to the origin of the quark masses and the mass generation is considered to be a prerogative of the weak interactions.” They then note than even the phase of an arbitrarily heavy quark would contribute to $`\theta `$ (which is just the phase of the determinant of the quark mass matrix), and that this would even apply to the mass of a Pauli-Villars regulator field! “Since $`\theta `$ is directly observable (see above),” they write, “it seems to be a unique example of how a laboratory experiment can shed light on regulator properties. Alternatively, it is very suspicious.”

Their suspicions aroused, they then proceed to construct a model identical to Kim’s. They conclude their paper with this paragraph: “If $`v`$ becomes arbitrarily large, the axion interaction with normal hadrons vanishes. (The same is true for its mass.) The new quark also becomes unobservable. This theoretical phantom still restores the natural P and T invariance of QCD. Although the model discussed is evidently a toy one, one might hope that something of the kind happens in unified theories.”

Here, then, is a clear presentation of the key point, by well known and respected physicists, published in the most widely read particle-physics journal of the time. Still, the initial impact was minimal. Why?

One clue can be found in this paragraph, near the top of the paper’s second page: “The overwhelming reaction to this observation is that the example of the one-instanton solution does prove CP violation in the presence of the $`\theta `$-term. Our final result complies with the conclusion on CP violation so that most readers can, justifiably, lose interest in the paper at this point.” Presumably, most readers did just that!

The fourth discovery of the arbitrariness of $`f_a`$ came in the paper “A simple solution to the strong CP problem with a harmless axion,” by Dine, Fischler, and myself, submitted to Physics Letters B in May 1981 and published in August 1981. We had been working on (very messy) models of supersymmetric technicolor, and had found an axion with a large decay constant in one of them. We began wondering whether this phenomenon could be duplicated in a simple Higgs model, and quickly constructed one (which turned out to be identical to Zhitnitskiĭ’s). We were excited at first, then wary. This was too easy! If it really worked, surely others would have thought of it before.

So we began trying to search the literature (without, of course, a computer database). We found Kim’s paper pretty quickly. We were disappointed, but luckily our specific model was different: we didn’t have an extra quark, only extra Higgs fields. We felt this would be considered more palatable. Also, we could fit our model into the GUT paradigm of the day. Our abstract reflects these points: “We describe a simple generalization of the Peccei-Quinn mechanism which eliminates the strong CP problem at the cost of a very light, very weakly coupled axion. The mechanism requires no new fermions and is easily implemented in grand unified theories.” As the paper was being typed (by a secretary, on a typewriter), we continued to poke around in the library, because we were still worried. Eventually we came across the SVZ paper. We were relieved to find that their model was the same as Kim’s. In our text, we talked about “Kim’s model,” and obviously we needed to change this to “the model of Kim and Shifman, Vainshtein, and Zakharov”. Too late! Our paper was typed, and re-typing was out of the question. We were at least able to add “see also SVZ” at the end of the reference to Kim. It came out looking as though we thought the SVZ paper was somehow inferior to the Kim paper. I don’t know if anyone ever noticed this, but I would like to take this occasion of Arkady’s sixtieth birthday to apologize to him and his collaborators!

Much to our surprise, our paper received a great deal of attention. Wise, Georgi, and Glashow constructed an explicit SU(5) model, and made the important point that solving the strong CP problem in this way did not require any extra fine-tuning in the Higgs sector. They also coined the phrase “invisible axion,” which has stuck. Nilles and Raby constructed a similar supersymmetric model. Many more papers followed. The two original models became known as the KSVZ and DFS models. Zhitnitskiĭ’s paper was still generally unknown and unreferenced.

In 1986, I had a (wrong, as it turned out) idea about the electric dipole moment of the neutron, and consulted Shalabin’s 1983 review. In it was a reference to Zhitnitskiĭ, right where I expected a reference to DFS. Wondering what this paper could be, I looked it up, began reading, and felt a flush of embarrassment wash over me as I realized that (a) it had a model identical to ours and (b) it was written and published long before ours. I began to reference it (and to ask others to reference it in axion papers that I refereed). I don’t know if this was the catalyst or not, but eventually the DFS model became the DFSZ model.

## 4 Axion Astrophysics and Cosmology

Detailed discussions of the axion couplings to hadrons, electrons, and photons can be found in Ref. \[\]. All are proportional to $`1/f_a`$, and the strongest limits on them come from astrophysical processes: axion emission from various objects (red giants, stars in globular clusters, supernova 1987A, etc.) must be slow enough to avoid changing the physics. The results are model dependent, but roughly require $`f_a>6\times 10^9`$GeV, or $`m_a<0.01`$eV. For a review, see Ref. \[\].

The axion can play a cosmological role for larger values of $`f_a`$. At high temperatures, the instanton effects which give rise to the axion mass go away. In the very early universe, then, the value of $`\theta `$ is undetermined, and presumably varies slowly from place to place. Inflation would pick out a particular value in our horizon volume. If subsequent reheating is to a temperature below $`f_a`$, this value is frozen in place until much later, when the universe cools and the axion mass appears. The axion field then begins coherent oscillations around its minimum. If the initial value of $`\theta `$ is $`O(1)`$, then the energy stored in these oscillations would overclose the universe if $`f_a`$ is larger than about $`10^{12}`$GeV. If $`f_a`$ is near this value, this energy (which is in the form of axions with near zero momentum) could be the cold dark matter.

Things are more complicated if there is no inflation, or reheating after inflation to a temperature above $`f_a`$. In this case a network of cosmic-string defects forms in the axion field, and the evolution of this network is a complicated numerical problem. The axions from string decay contribute at least as much dark matter as the axions from the initial misalignment, and possibly much more; see Ref. \[\] for more details.

The possibility that axions form the cold dark matter leads to the exciting possibility that we might be able to detect them. The axion-photon interaction lagrangian is

$$L_{a\gamma \gamma }=\frac{\alpha }{2\pi f_a}(CC^{})a𝐄𝐁,$$
(12)
where $`C`$ and $`C^{}`$ are numerical constants. This interaction allows an axion to convert to a photon in a magnetic field, and Sikivie proposed using this effect to search for dark-matter axions. The value of $`C`$ is computed in terms of the Peccei-Quinn charges of the fermions, and $`C^{}`$ arises from axion-pion and axion-eta mixing; it is model independent, and given by

$`C^{}`$ $`=`$ $`{\displaystyle \frac{2}{3}}{\displaystyle \frac{m_um_d+4m_dm_s+m_sm_u}{m_um_d+m_dm_s+m_sm_u}}`$ (13)
$`=`$ $`1.93\pm 0.04.`$
Unfortunately, $`C`$ is typically positive, so there is a cancelation. In grand-unifiable models (that is, models where the fermions come in SU(5) multiplets with the same PQ charge; whether the extra heavy particles implied by unification exist or not is irrelevant) such as the DFSZ model,

$$C_{\mathrm{DFSZ}}=\frac{8}{3}.$$
(14)
In KSVZ models with a single heavy quark of charge $`Q`$,

$$C_{\mathrm{KSVZ}}=6Q^2.$$
(15)
If axions are the dark matter, $`m_a1`$ to $`100\mu `$eV is the most interesting mass range ($`1\mu \mathrm{eV}=10^6`$eV). Experiments to search for them are currently underway, and one has ruled out KSVZ axions (with $`Q=0`$) in the mass range $`2.8\pm 0.5\mu `$eV as a significant component of the local dark matter. The experiments are continuing; for a review, see Ref. \[\].

## 5 The Future of Axion Physics

Over twenty years have passed since the invention of the invisible axion, and we still do not know whether or not this is the correct solution to the strong CP problem. It is vitally important that the current searches continue until they have covered at least the most plausible mass and coupling ranges; it would be a great shame if such important physics surrounded us, and we left it undiscovered. Of course, even if these searches do not find dark matter axions, this means only that the dark matter is not axions, and not that axions do not exist. There are many cosmological scenarios for this possibility.

Also, it should be noted that axions arise rather naturally in superstring models, and it may be that any ground state of string/M theory that resembles the Standard Model always includes an axion with, say, $`f_a(v_{\mathrm{EW}}M_{\mathrm{Planck}})^{1/2}10^{11}`$GeV. However at present we seem a long way from being to reach this sort of conclusion.

The road ahead for axion physics is thus likely to be a hard one, both theoretically and experimentally. But the reward for a successful traversal will make the journey worthwhile.

I would like to thank the organizing committee for Arkadyfest—Keith Olive, Mikhail Shifman, and Mikhail Voloshin—for the opportunity to present this talk, and Arkady Vainshtein for providing the occasion and the inspirational physics. I also acknowledge the NSF for financial support.