text:
# Time dependent diffusion in pulsating white dwarf stars: Asteroseismology of G117-B15A

## 1 Introduction

Asteroseismology has become a powerful method to disentangle the internal structure and evolution of stars by means of the study of their oscillatory pattern. This technique, very sophisticated in the case of the Sun, has also undergone a strong development in other stars, in particular in variable white dwarf (WD) stars (see e.g. Brown & Gilliland 1994; Gautschy & Saio 1995; 1996).

Pulsating WDs show multiperiodic luminosity variations in three ranges of effective temperatures ($`T_{\mathrm{eff}}`$ ) corresponding to the currently called DOV, DBV and DAV (see, e.g., the review by Winget 1988). Of interest in this work are the DAVs (hydrogen - dominated atmospheres), or ZZ Ceti, that pulsate in the instability strip corresponding to 11000 K $``$ $`T_{\mathrm{eff}}`$ $``$ 13000 K. The periodicities in the light curves of pulsating WDs are naturally explained in terms of nonradial g-modes of low harmonic degree ($`\mathrm{}2`$), driven by the “$`\kappa `$ mechanism” working in a partial ionization region near to the stellar surface (Dolez & Vauclair 1981, Winget et al. 1982). The periods ($`P`$) are found within a range of $`100sP1200s`$ and photometric amplitudes reach up to 0.30 magnitudes.

Most asteroseismological studies performed on WDs to date rely on stellar models constructed under some simplifying hypotheses. One of the most relevant is related to the profile of the internal chemical composition in the interface zones. In such regions, the equilibrium diffusion in the trace element approximation has been widely employed to infer the profile of the chemical distribution (see Tassoul, Fontaine & Winget 1990). The main motivation for considering this approximation is to avoid the solution of time dependent element diffusion as the WD evolves. In the frame of such an approximation, the profile of the interface region is very simple: its functional form is a power law. The transition zone is separated into two parts: an upper one in which one element is dominant and the other is considered as a trace and a lower region in which the role of the respective elements is reversed. Because these two power law solutions are matched for fulfilling the condition of conservation of mass of each element, a discontinuity in the derivative occurs just at the matching point. The exponent of the power law solution is directly related to the state of ionization of the stellar plasma. Thus, the structure of the whole interface zone, in the frame of this standard treatment, can be modified only if the plasma suffers from a modification in the state of ionization as a result of stellar evolution. Calculations of asteroseismology of WDs in the frame of such a standard treatment for the chemical interfaces are those of, e.g., Bradley (1996, 1998, 2001), Bradley & Winget (1994), Brassard et al. (1991, 1992ab), Fontaine et al. (1992), Montgomery & Winget (1999), Metcalfe, Nather & Winget (2000), Metcalfe, Winget & Charbonneau (2001), Montgomery, Metcalfe & Winget (2001).

From an evolutionary point of view, the shape of the chemical interfaces may not be critical111An obvious and important exception is when the tail of chemical profiles are subject to nuclear burning., but in pulsational studies they provide a non - negligible contribution to the shape of the Ledoux term of the Brunt-Väisälä frequency (Brassard et al. 1991). Thus, we can, in principle, expect differences between studies with equilibrium diffusion and those in which other more physically sound treatment is performed. This is a very important point to be made in connection with the aim of the present work.

Since the pioneering work of Iben & MacDonald (1985), we know that element diffusion modifies the chemical abundance distribution within a WD star even during evolutionary stages corresponding to the ZZ Ceti domain (see Iben & MacDonald 1985, particularly their Fig. 4). Few calculations exist in the literature in which the evolution of WDs is addressed in a self consistent way with time dependent element diffusion. Amongst them, we mention the study of MacDonald, Hernanz & José (1998) aimed at studying the carbon pollution in cool WDs. Also, Dehner & Kawaler (1995) used non - equilibrium diffusion together with evolutionary calculations to study WDs with helium envelopes. Finally, Althaus, Serenelli & Benvenuto (2001) have recently shown that diffusion induces the occurrence of thermonuclear flashes in helium - core WDs, causing the evolution of such kind of WDs to occur on timescales significatively shorter than predicted by models without diffusion. This has been particularly important in solving the discrepancy (Van Kerkwijk et al. 2000) about the age of binary systems containing a millisecond pulsar and a helium WD.

As far as we are aware, the only work aimed at exploring the role played by diffusion in the period and its rate of change for $`g`$-modes in pulsating DA WDs is that of Córsico et al. (2001b). In particular, the authors found that the differences in the shape of the chemical profiles at the interface zones induce appreciable changes in the periods, as compared to the case of equilibrium diffusion in the trace element approximation. Also, there are noticeable changes in the period derivative, which are due in part to the evolution of the chemical profile during the cooling of the WD across the ZZ Ceti instability strip.

Since sometime ago there have been several works available in the literature in which the observed period structure of a particular object is fitted to theoretical predictions. Such works show that, in principle, information about the stellar mass and the stratified outer layer structure can be inferred (see, e.g., Winget et al. 1991; Fontaine et al. 1992; Pfeiffer et al. 1996; Bradley & Kleinman 1997; Bradley 1998, 2001; Bradley & Winget 1994). As these studies are based on a simplified treatment of diffusion, we believe that it is worth revisiting this problem on the basis of more detailed models of WD structure. In particular, in view of the results found in Córsico et al. (2001b), we expect to find differences in the period fitting to a particular object, as compared with the situation in which the standard treatment is used. The only way to find how important can be such differences is by performing a detailed asteroseismological study choosing a well studied object.

We consider G117-B15A as an optimal target for our study. G117-B15A is an otherwise typical DA WD, the variability of which was discovered by Mc Graw & Robinson (1976) and, since then, it has been monitored continuously. The mass and in particular the effective temperature of this star have been the subject of numerous spectroscopic redeterminations (see Gautschy, Ludwig & Freytag 1996 for a summary), in particular values of 0.59 $`M_{}`$and 11620 K, respectively, have been derived by Bergeron et al. (1995). More recently, Koester & Allard (2000) (hereafter KA) have suggested a somewhat lower value for the mass of 0.53 $`M_{}`$and a higher effective temperature of $`T_{\mathrm{eff}}`$ =$`11900\pm 140`$ K.

G117-B15A has periods of oscillation of 215.2, 271 and 304.4 s (Kepler et al. 1982). Notably, for the 215.2 s mode it has been possible to find a value of its temporal derivative (Kepler et al. 2000) to be $`\dot{P}=(2.3\pm 1.4)\times 10^{15}`$ s s-1. Interestingly, the 215.2 s mode present in G117-B15A is the most stable optically registered oscillation with a stability comparable to that of the most stable millisecond pulsars. As G117-B15A is a well know oscillator, it has motivated the interest of several researchers in its internal structure. For example, Bradley (1998) has found the best fitting to the period pattern with a model of $`0.6`$ $`M_{}`$and depending on the identification of the modes (see below for an explanation of its meaning), the favoured value for the hydrogen envelope mass fraction $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ is $`7<\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})<6`$ or $`5<\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})<4`$, and $`\mathrm{log}(\mathrm{M}_{\mathrm{He}}/\mathrm{M}_{})`$ =-2. More recently, Córsico et al. (2001a) have found, on the basis of WD models in which diffusion is neglected, the best fitting to the period pattern with a 0.55 $`M_{}`$WD model with a carbon - oxygen interior, $`\mathrm{log}(\mathrm{M}_{\mathrm{He}}/\mathrm{M}_{})`$ =-2 and $`\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})`$ =-4 as predicted by stellar evolution. Notably, such a fitting is in nice agreement with one of those proposed by Bradley (1998).

It is the aim of this work to perform a fitting to the period structure present in G117-B15A by computing the non-radial eigenmodes in the frame of linear and adiabatic approximation, and evolutionary WD models in which time dependent element diffusion is properly accounted for. About the way we shall perform such a fitting, some words are in order. In handling models like those we shall employ it is quite obvious that we are not in a position to employ techniques like those used in the so - called genetic algorithm (see its application for the case of the DBV GD358 in Metcalfe et al. 2000). In the context of ZZ Ceti stars a similar approach though to a less extent has been repeatedly employed by numerous investigations (e.g the complete and detailed study performed by Bradley 1998, 2001). Up to now, the shape and thickness of the relevant chemical interfaces are usually treated as a free parameter. However, a more physically sound treatment can be performed when account is made of time dependent element diffusion in evolutionary models of WDs. This aspect of WD evolution is one of the most important when an attempt is made to compare observations with theoretical expectations from pulsating WDs. Unfortunately, there are currently some uncertainties (e.g., the treatment of convection, the rate of the critical nuclear reaction $`{}_{}{}^{12}\mathrm{C}(\alpha ,\gamma ){}_{}{}^{16}\mathrm{O}`$, wind mass loss episodes, etc.) that prevent us, provided a mass value and pre-WD composition, from predicting a whole definite internal structure. This is especially true for $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ . Thus, in view of these facts, we shall consider the value of $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ as a free parameter. We assumed $`\mathrm{log}(\mathrm{M}_{\mathrm{He}}/\mathrm{M}_{})`$ =-2 throughout this paper.

The remainder of the paper is organized as follows: In Section 2 we discuss the physical ingredients we have employed together with the computational strategy we have employed. Section 3 is devoted to presenting the details of the evolutionary models we have constructed down to the conditions relevant for the WD we are investigating. In section 4 we present the asteroseismological analysis we performed for the particular case of G117B-15A. Section 5 is devoted to the discussion of the results we have found and finally in Section 6 we give some concluding remarks and insights for future investigations.

## 2 Details of the Computations

For this work we have employed the same code and physical ingredients as in Córsico & Benvenuto (2002) and Córsico et al. (2001ab). In particular for the diffusion processes we considered gravitational settling, chemical and thermal diffusion as in Althaus & Benvenuto (2000).

In our code, evolutionary and pulsational calculations are performed in a fairly automatic way. After selecting a starting stellar model we choose an interval in $`P`$ and $`T_{\mathrm{eff}}`$ . The evolutionary code computes the model cooling until the hot edge of the $`T_{\mathrm{eff}}`$ -interval is reached. Then, the program calls the set of pulsation routines beginning the scan for modes. When a mode is found, the code generates an approximate solution which is iteratively improved to convergence (of the eigenvalue and the eigenmode simultaneously) and stored. This procedure is repeated until the period interval is covered. Then, the evolutionary code generates the next stellar model and calls pulsation routines again. Now, the previously stored modes are taken as initial approximation to the modes of the subsequent model and iterated to convergence. Such a procedure is automatically repeated for all evolutionary models inside the chosen $`T_{\mathrm{eff}}`$ \- interval. The computational strategy described above has been successfully applied in fitting the observed periods of G117-B15A to impose constraints on the mass of axions (Córsico et al. 2001a) in the frame of WD models less detailed than those we shall employ in the present work. We refer the reader to Córsico & Benvenuto (2002) and Córsico et al. (2001ab) and references therein for further details. Let us quote that most of our evolutionary models have been divided in more than 2000 meshpoints, whereas for mode calculations we employed 5000 meshpoints.

In order to start our evolutionary calculation we have employed an artificial heating technique detailed in Althaus & Benvenuto (2000). Let us quote that in this technique, the stellar model is forced to undergo some unphysical evolution by introducing an artificial energy release that, when the star gets bright enough, is switched off smoothly. Then, the star relaxes to the physical cooling branch after few tens of models.

In employing this artificial technique, we must be very careful in avoiding generating an artificial chemical profile. Thus, in order to verify that the star not only relaxes to the physical cooling branch but simultaneously relaxes to a correct internal chemical profile, we have employed, for a fixed $`\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})`$ value, different shapes of the initial hydrogen and helium profiles. We found that, if the initial model is far away enough from the evolutionary stage we are interested in, the fine details of the initial profile are irrelevant because diffusion largely evolves it to a well defined one. As a matter of fact, the structure and the chemical profiles are physically plausible far before the star attains the effective temperature range relevant for G117-B15A.

A better approach for constructing pulsating WD models would obviously be to start computing evolutionary models from much earlier in the star history. This possibility, better from a physical point of view, again would require an enormous amount of computing time for a definite sequence. Thus trying to fit a period pattern in such a way is, in our opinion, far beyond what seems reasonable at the present state of the art. In any case, it would be interesting to explore the pulsational properties of models resulting from a full evolutionary calculation. We plan to present such a study in a future paper.

To be specific, we have started the calculations with models of $`\mathrm{log}(\mathrm{L}/\mathrm{L}_{})`$ $``$ 2.05, $`\mathrm{log}T_{\mathrm{eff}}`$ $``$4.46 and masses of 0.50, 0.55, and 0.60 $`M_{}`$(in performing the fitting of the models to observations we have also computed supplementary models of 0.525 and 0.5375 $`M_{}`$). The mass fraction embraced by the outermost hydrogen layer has been assumed in a way that, after evolving along the cooling branch down to $`T_{\mathrm{eff}}`$ =12500 K the values of $`\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})`$ present in the star are those given in Table 1. At the starting conditions we have assumed the Salaris et al. (1997) core composition, whilst for the outer layers the CNO abundances were those corresponding to the pre-WD evolution. In order to get different values of $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ , we have started from the model with the largest $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ and simply replaced 1H by 4He at the base of the hydrogen envelope. Let us quote that the initial model has internal temperatures low enough that prevent the occurrence of nuclear reactions up to the moment at which the WD model has relaxed to the physical branch. We should mention that models become physically adequate (after the unphysical transitory) for describing a WD at $`\mathrm{log}(\mathrm{L}/\mathrm{L}_{})`$ $``$ 1.4, $`\mathrm{log}T_{\mathrm{eff}}`$ $``$4.83, far from the conditions relevant for DAV WDs. Also, the values of $`\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})`$ at the instability strip are somewhat lower than those included in Table 1 simply because nuclear reactions are still operating at such evolutionary stages. However this effect is noticeable only for the case of the thickest hydrogen envelopes on the most massive objects considered here.

Table 1. Values of $`\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})`$ for our computed models

at $`T_{\mathrm{eff}}`$ =12,500 K    $`0.50M_{}`$ $`0.55M_{}`$ $`0.60M_{}`$ $`0.525M_{}`$ $`0.5375M_{}`$    -3.815 -3.862 -3.941 -3.831 -3.843  -4.193 -4.224 -4.249 -4.198   -4.679 -4.684 -4.692    -5.175 -5.177 -5.180    -5.671 -5.671 -5.672    -6.160 -6.158 -6.156 -6.159   -6.640 -6.633 -6.700 -6.634   -7.071 -7.047 -7.028    -7.423 -7.385 -7.349

To our notice, this is the first opportunity in which carbon - oxygen, intermediate mass WD models with a set of values in $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ are evolved in the frame of time dependent element diffusion. Prior to this work, the kind of calculations carried out here were performed only for the values of $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ predicted by stellar evolution theory (Iben & MacDonald 1985; 1986).

## 3 Evolutionary Results

For each initial model, we have computed its evolution down to the conditions corresponding to the DAV instability strip. In doing so, we computed about ten thousand models for each sequence. Perhaps the most relevant characteristic of the models for the kind of study we are performing here is the internal chemical profile. Results representative for the case of the models corresponding to the DAV instability strip are presented in Fig. 1. In this figure we depict the internal chemical profiles for a 0.55 $`M_{}`$WD model at $`\mathrm{log}T_{\mathrm{eff}}`$ $``$ 4.096 for $`\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})`$ =-3.862. The adopted effective temperature is slightly higher than the last estimation for G117B-15A ($`\mathrm{log}T_{\mathrm{eff}}`$ $``$ 4.076). Thus, these profiles should be representative of the conditions expected for the interior of G117B-15A in the frame of our set of models.

Let us discuss the main characteristics of these profiles. Before discussing them, we have to warn the reader that the way we have modified the value of $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ is ad-hoc. In fact, it is not clear at present if there exists, for a fixed stellar mass value, an evolutionary mechanism capable to produce a distribution of $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ values differing each other by orders of magnitude. In addition, the outward tail of the distribution of initial heavy element abundances may be different from the results we shall present below. However, we should remark that the profile of the light elements present in the outer parts of the models should be essentially correct.

With regard to the chemical profiles of our models we want to mention that for $`\mathrm{log}(1\mathrm{M}_\mathrm{r}/\mathrm{M}_{})>1.6`$ the carbon - oxygen distribution is that predicted by pre-WD evolutionary calculations and it is the same for all of our evolutionary sequences. However this is not the case for the chemical structure of the outer layers, which are, in turn, crucial for the pulsational properties of the models. The main variations are related to the whole distribution of 1H, 3He and the outwards tail of 4He. While a variation in the assumed value of $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ obviously implies a modification in the profiles, we find a very noticeable variation in the slope of the profiles for the isotopes present in the outer layers of our WD models. In the model with $`\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})`$ =-3.815 it is remarkable the large slope in the tail of the hydrogen distribution. In spite of this steep profile, degeneracy prevents hydrogen from diffusing further inwards. In models with a smaller hydrogen content, the tail of the hydrogen distribution is smoother. Note also that the shape of the tail of the hydrogen distribution at the evolutionary conditions relevant for DAV WDs is strongly dependent on the value of the assumed $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ .

Regarding the tail of the 4He distribution, this is largely due to gravitational settling that forces 4He to sink down. In particular our models also present a zone in which the chemical composition is largely dominated by helium. Also, it is of interest the distribution of 3He. It is clearly noticeable that most of the 3He is located in stellar layers at which 1H and 4He are present. This is a consequence of the interplay between nuclear reactions that form 3He by means of the first two reactions of the PPI chain (Clayton 1968)

$`{}_{}{}^{1}\mathrm{H}+{}_{}{}^{1}\mathrm{H}{}_{}{}^{2}\mathrm{H}+e^++\nu `$

$`{}_{}{}^{2}\mathrm{H}+{}_{}{}^{1}\mathrm{H}{}_{}{}^{3}\mathrm{He}`$

and gravitational settling. Notice that the outer tail of the 3He distribution is similar to the corresponding to 4He which is, in part, due to the operation of the reaction

$`{}_{}{}^{3}\mathrm{He}+{}_{}{}^{3}\mathrm{He}{}_{}{}^{4}\mathrm{He}+2{}_{}{}^{1}\mathrm{H}`$.

As the chemical profile of the hydrogen in surface layers is one of the key ingredients in determining the g-modes of DAV WDs, we show in Fig. 2 the base of the hydrogen envelopes for our suite of 0.55 $`M_{}`$models at the ZZ Ceti domain. In addition, we show the profiles predicted by the standard treatment of equilibrium diffusion in the trace element approximation for the same stellar mass and $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ values. It is clear that time dependent element diffusion predicts hydrogen profiles that are somewhat different from those given by the standard treatment. In the standard treatment, the slope of the profile is determined solely by the state of ionization of the plasma (see Eqs. 26 - 29 of Tassoul et al. 1990). Thus, if the state of ionization is the same, the slope should remain essentially unaltered with further evolution. We have verified that, even in the case of the thinnest hydrogen envelopes we have considered, the base of the hydrogen envelope is located at layers with temperatures above $`10^6`$K, which ensures complete ionization for hydrogen and helium and this is so even at the effective temperatures corresponding to the DAV instability strip. Also, notice that the sudden change of the slope at the matching point of the two solutions of the standard treatment (denoted by circles along the curves) has a direct consequence on the enhancement of the trapping of some modes in the hydrogen envelope (see, e.g., Brassard et al. 1992a, Bradley 1996). Thus, our results show the necessity of allowing for a more detailed treatment of the chemical evolution at performing asteroseismology of WDs.

Now, let us discuss the impact of our chemical profiles on the Brunt-Väisälä frequency, a key quantity entering as a coefficient in the equations of non-radial pulsations. This is defined as (Unno et al. 1989)

$$N^2=g\left(\frac{1}{\mathrm{\Gamma }_1}\frac{d\mathrm{ln}P}{dr}\frac{d\mathrm{ln}\rho }{dr}\right)$$
(1)
where all the symbols have their usual meaning222In this Section $`P`$ stands for the pressure but in the rest of the paper it stands for the period of oscillation. and $`\mathrm{\Gamma }_1`$ is the first adiabatic exponent. However, in order to avoid numerical noise, we use the expression given by Brassard et al. (1991):

$$N^2=\frac{g^2\rho }{P}\frac{\chi _T}{\chi _\rho }\left[_{ad}+B\right].$$
(2)
The Ledoux term $`B`$, for the case of a M-component plasma, is given by

$$B=\frac{1}{\chi _T}\underset{i=1}{\overset{M1}{}}\chi _{X_i}\frac{d\mathrm{ln}X_i}{d\mathrm{ln}P}.$$
(3)
where

$$\chi _\rho =\left(\frac{\mathrm{ln}P}{\mathrm{ln}\rho }\right)_{T,\{X_i\}}$$

$$\chi _T=\left(\frac{\mathrm{ln}P}{\mathrm{ln}T}\right)_{\rho ,\{X_i\}}$$

$$\chi _{X_i}=\left(\frac{\mathrm{ln}P}{\mathrm{ln}X_i}\right)_{\rho ,T,\{X_{ji}\}}.$$

In Fig. 3 we show the squared Brunt-Väisälä frequency and the Ledoux term for models of 0.55 $`M_{}`$with $`\mathrm{log}T_{\mathrm{eff}}`$ $``$4.096 and all the considered values of $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ . We also show in an inset the Brunt-Väisälä frequency at the hydrogen - helium interface for the case of the thickest hydrogen envelope according to our computations and to the prediction of equilibrium diffusion (thick and thin lines, respectively). In order to show the correspondence of the structure of these functions with the internal chemical profiles of the star we have also included the curves corresponding to the abundances of 4He and 12C.

The Ledoux term $`B`$ has a structure that remains largely invariant for the stellar core under changes in $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ . From the centre outwards, the three peaks are due to steep slopes in the carbon profile. The first two are due to the structure of the carbon - oxygen interface whereas the other is due to the helium - carbon interface. In contrast to the internal behaviour of the Ledoux term, the other peak, due to the hydrogen - helium transition, is largely modified, not only in position but also in its height. Notice that the lower $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ the thinner the peak (see Brassard et al. 1991 and Tassoul et al. 1990 for comparison of the shape of the $`B`$ term and Bradley 1996 for the profile of the Brunt-Väisälä frequency).

It is worthwhile to comment on the fact that because of our realistic chemical profiles, we would expect to find differences in the mode trapping properties as compared with previous studies. This important issue surely deserves a careful exploration, which will be reported in a subsequent communication.

## 4 Asteroseismological Results

An important aspect for the present work is the so-called mode identification. That is, the identification of the $`\mathrm{}`$ and $`k`$ values corresponding to each observed period. In the case of G117-B15A, the results of Robinson et al. (1995) indicate that the 215.2 s period corresponds to a dipolar mode ($`\mathrm{}=1`$). Following the work by Bradley (1998) we shall assume that the other two modes cited above are also dipolar (Brassard et al. 1993; Fontaine & Brassard 1994), and also that the other periodicities present in the light curve of the star are not associated with actual eigenmodes but are actually due to non-linear effects in the envelope (Brassard et al. 1993). With regard to the radial order of the modes, there exist two possible identifications: following Clemens (1994) the observed periods are dipolar modes with $`k=2`$ (215.2 s), $`k=3`$ (271 s) and $`k=4`$ (304.4 s). On the other hand, Fontaine, Brassard & Wesemael (1994) identified these modes as dipolar with $`k=1`$ (215.2 s), $`k=2`$ (271 s) and $`k=3`$ (304.4 s). According to these facts, in computing the g-modes of our WD models, we have considered, as previously mentioned, only $`\mathrm{}=1`$ and radial order $`k=1,\mathrm{},4`$. Higher radial orders correspond, for our stellar mass values, to periods long in excess compared to those observed and thus are not considered here.

In order to look for a fit to observations we shall consider a function $`\psi `$ defined as

$$\psi =\mathrm{log}\underset{i=1}{\overset{3}{}}\left[1\frac{P_i^C}{P_i^O}\right]^2,$$
(4)
where $`P_i^O`$ and $`P_i^C`$ are the i-th observed and calculated periods respectively. For the case of G117B-15A, $`P_i^O`$= 215.2, 271 and 304.4 s. As a result of the pulsational calculations we have the values of the periods for $`k=1,\mathrm{},4`$. Then, assuming the two reasonable possible identifications of the observed modes (i.e. $`P_i^O`$ corresponding to $`k=1,2,3`$, hereafter case A, and $`P_i^O`$ corresponding to $`k=2,3,4`$, hereafter case B) we can compute the function $`\psi `$ in terms of the effective temperature of the star for each of the considered sequences. Obviously in the case that the fitting were exact, $`\psi \mathrm{}`$. Thus the criterium we shall assume in selecting the best model for representing G117B-15A is simply to look for the minima of the $`\psi `$ function considering the whole set of models. In such a way we shall be able to determine not only the internal structure and mass of G117B-15A, but also its asteroseismological effective temperature.

In Fig. 4 we show $`\psi `$ as a function of the effective temperature for the case of the set of models with 0.50 $`M_{}`$. In this figure the left panel labeled as A corresponds to the identification of the observed modes with $`k=1,2,3`$ whereas the right panel B corresponds to the other identification of modes we shall consider: $`k=2,3,4`$ modes. As it has been found by Bradley (1998), for cases A and B the best fitting to observations are found for models with values of $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ very different. While assuming case A we find the best fit with rather thin hydrogen envelopes ( $`\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})`$ = -6.16) and in the case B the best fit is found for the case of the thickest considered hydrogen envelope: $`\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})`$ = -3.81. It is remarkable that the minimum found for case A is a very deep one, corresponding to a mass and effective temperature ($`T_{\mathrm{eff}}`$ $``$11400K) rather lower than those corresponding to the last determination of these parameters for G117B-15A.

If we increase the stellar mass to 0.55 $`M_{}`$, the behaviour of the function $`\psi `$ is largely changed (see Fig. 5). Despite that in case A the best fitting is for $`\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})`$ = -6.63, it corresponds however to an effective temperatures far lower than that suggested by KA (and outside the DAV instability strip). By contrast, in case B the best fit is again provided by the model with the thickest hydrogen envelope, but again it occurs at effective temperatures somewhat too low.

Finally, in Fig. 6 we show the behaviour of function $`\psi `$ for the case of the models constructed assuming a stellar mass value of 0.60 $`M_{}`$. It is clear that in both cases of mode identifications the fits are far poorer than those for the other, lower mass values. Thus, the results presented up to this point strongly suggest that the asteroseismological mass of G117B-15A should be in between 0.50 $`M_{}`$and 0.55 $`M_{}`$. Moreover, the only viable mass for the case A seems to be 0.50 $`M_{}`$, because larger mass values give rise to poorer fits.

Let us analyse the results corresponding to case B more carefully, centering our attention in models with the thickest hydrogen envelopes. From Fig. 4b, it can be noticed that in the case of 0.50 $`M_{}`$model, the fit is better for increasing effective temperatures. However, notice that (Fig. 5b) for the 0.55 $`M_{}`$model the situation is just the opposite: the fit improves for lower effective temperatures. Guided by this fact, it can be expected a minimum in $`\psi `$ for some intermediate value of the stellar mass. In order to explore this possibility we computed models with 0.525 $`M_{}`$restricting ourselves to values of $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ similar to those that provided good fits (see Table 1). For this mass value, the identification of modes in case A shows a good fit for $`\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})`$ = -6.63 at $`T_{\mathrm{eff}}`$ $``$11400 while case B presents an interesting minimum at stellar mass and effective temperature values very similar to those determined by KA for this object by spectroscopic analysis. Finally, in order to perform a finer analysis, we compute a sequence of models with 0.5375 $`M_{}`$and a thick hydrogen envelope (see Table 1). The results for the whole set of masses and thickest hydrogen envelopes for case B is shown in Fig. 7. The labels of each curve correspond to the values of $`\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})`$ presented in Table 1.

The results presented in this Section indicate that, if we expect the asteroseismological masses and effective temperatures to be similar to those determined spectroscopically, the best fit is found in the frame of case B for a model with a stellar mass of 0.525 $`M_{}`$. For the sake of completeness, we also show in Fig. 8 the periods of oscillation for radial order $`k=2,3,4`$ for the case of models with the thickest hydrogen envelope for each considered value of the stellar mass.

## 5 Discussion

From the results presented in the preceding Section, we conclude that the best fit to observations is found with stellar mass and $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ values that are very sensitive to the adopted mode identification. This has been previously found by Bradley (1998). If we restrict ourselves to analyse the case of G117B-15A only in the frame of asteroseismological models we find no way to be conclusive about the actual values for its stellar mass, $`\mathrm{M}_\mathrm{H}/\mathrm{M}_{}`$ , and effective temperature. However, if we incorporate the best available spectroscopic determinations of stellar mass and effective temperatures to this analysis and expect them to be compatible with the asteroseismological results, we can discriminate between the identification of modes and also on the other values for the characteristics of the models.

In our opinion, the best fit found in the frame of case A mode identification, though being very good, ($`\psi `$ shows a very deep minimum) occurs however at a temperature and mass values too low compared to those determined by KA. Nevertheless, for case B we find also a very deep minimum at an effective temperature compatible to that of KA and a mass value almost identical to that given in KA. If we accept that this coincidence is not by chance, we have to conclude that this is the best fit to the structure and evolutionary conditions for G117B-15A. Thus, our results strongly suggest that the hydrogen envelope present in G117B-15A is a thick one with $`\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})`$ $``$ -3.8.

The result regarding the actual amount of hydrogen present in G117B-15A should be considered as a very interesting one. In fact this is a very thick hydrogen envelope. Notice that this is the upper limit we considered for this quantity. Because of this fact, in improving our fit we may have considered even thicker envelopes (see Figs. 4-6). However, notice that performing asteroseismology in the way we have done, it is not clear if such a procedure would produce stellar models plausible from a physical point of view. Notice that even in the hypothetical case that our artificial procedure for generating the stellar models were able to handle such models, there is no way to be confident that previous nuclear burning were effective enough to prevent the existence of such structures in Nature. Because of this reason we think it would be a more safe procedure to compute the structure of objects of the mass of G117B-15A starting from the hydrogen burning main sequence throughout all the evolutionary stages previous to the WD stage. This is the only way we envisage in order to assess the possible existence of hydrogen envelopes thicker than those we have considered here. In any case, we should point out that previous detailed calculations of the evolution of objects that finally form a WD predict values near $`\mathrm{M}_\mathrm{H}=10^4M_{}`$ for 0.60 $`M_{}`$objects at the start of the cooling branch. However a serious study of this problem would carry us too far afield and is beyond the scope of the present paper.

Finally, we can revisit the estimation of the parallax of G117B-15A as done in Bradley (1998). Because the effective temperature is higher (11790K vs. 11600K) and the stellar mass lower (0.52 $`M_{}`$vs. 0.60 $`M_{}`$, and thus a larger radius), the luminosity ($`\mathrm{log}(\mathrm{L}/\mathrm{L}_{})`$ =-2.44 vs. $`\mathrm{log}(\mathrm{L}/\mathrm{L}_{})`$ =-2.52) is higher than predicted by Bradley (1998) and for the same visual apparent magnitude $`m_V=11.50`$ and bolometric correction of -0.611 mag (Bergeron, Wesemael & Beauchamp 1995) we find that G117B-15A should have a parallax of 15.89 mas whereas Bradley determined it to be 16.5, which is still too high to be compatible (to a $`1\sigma `$ level) with the determination of $`10.5\pm 4.2`$ mas by van Altena, Lee & Hoffleit (1994). In view of the amount of compatible evidence on the value of the stellar mass and effective temperature we expect the observed parallax of G117B-15A to be underestimated. It seems difficult to adjudicate the discrepancy in parallax as due to an overestimation of the value of the stellar mass. Were this the case, it should be significantly lower than found in this work.

Table 2. Summary of results for the best fitting

| Quantity | Fitted model | G117-B15A |
| --- | --- | --- |
| $`M_{}/M_{}`$ | 0.525 | 0.53 (1) |
| $`\mathrm{log}g`$ | 7.85 | 7.86 $`\pm `$ 0.14 (1) |
| $`T_{\mathrm{eff}}`$ \[K\] | 11,790 | 11,900 $`\pm `$ 140 (1) |
| $`\mathrm{log}(M_\mathrm{H}/M_{})`$ | -3.83 | $`\mathrm{}`$ |
| $`\mathrm{log}(M_{\mathrm{He}}/M_{})`$ | -2.00 | $`\mathrm{}`$ |
| Parallax \[mas\] | 15.89 | 10.5 $`\pm `$ 4.2 (2) |
| $`P(k=2)`$ \[s\] | 208.80 | 215.2 (3) |
| $`P(k=3)`$ | 278.85 | 271.0 (3) |
| $`P(k=4)`$ | 308.57 | 304.4 (3) |
| $`\dot{P}(k=2)`$ \[$`10^{15}\mathrm{s}\mathrm{s}^1`$\] | 4.43 | $`2.3\pm 1.4`$ (4) |
| $`\dot{P}(k=3)`$ | 3.22 | $`\mathrm{}`$ |
| $`\dot{P}(k=4)`$ | 5.76 | $`\mathrm{}`$ |

References: (1) Koester & Allard (2000), (2) van Altena et al. (1994), (3) Kepler et al. (1982), (4) Kepler et al. (2000)

## 6 Conclusions

In this paper we have studied the structural characteristic of the variable DA white dwarf (WD) G117B-15A by applying the methods of asteroseismology. In doing so we have constructed models of WD evolution considering updated and detailed physical ingredients. It should be remarked than we have included several processes responsible for the diffusion of elements in the WD interior. In particular, we considered gravitational settling, chemical and thermal diffusion. Starting from an artificial model, we have considered several thickness for the outermost hydrogen layer, whereas for the inner helium-, carbon- and oxygen-rich layer we considered the profiles predicted by Salaris et al. (1997). The range of stellar masses we have considered is $`0.50M_{}/M_{}0.60`$. As far as we are aware, this is the first study in which evolutionary models of DA WDs considering element diffusion have been constructed for a set of values of hydrogen mass fraction.

The evolution of each model sequences were followed down to effective temperature of 12500K from where on, we considered the evolution coupled to the oscillations. We considered dipolar, adiabatic g-modes with radial order $`k=1,\mathrm{},4`$. After constructing the full set of models, we considered the location of the minima of a function conveniently defined for the purpose of fitting. We considered two possible mode identifications according to Bradley (1998): that observed modes 215.2, 271 and 304.4 s correspond to $`k=1,2,3`$ and $`k=2,3,4`$. We find that the employment of asteroseismology does not provide an univocal answer about the correct identification of modes present in G117B-15A. However, if we use them together with data deduced from spectroscopy it is found that the identification of the observed periods of 215.2, 271 and 304.4 s with dipolar modes with $`k=2,3,4`$ is clearly the best.

The results presented in this work (summarised in table 2) strongly suggest that G117B-15A is a DA WD with a stellar mass of 0.525 $`M_{}`$, a hydrogen mass fraction $`\mathrm{log}(\mathrm{M}_\mathrm{H}/\mathrm{M}_{})`$ =-3.83 and effective temperature $`T_{\mathrm{eff}}`$ $``$11800K. Notably, the values of the effective temperature and stellar mass are in very nice agreement with those predicted by spectroscopic analysis by Koester & Allard (2000). This represents the main result of the present paper.

While the favoured value for the mass fraction of hydrogen is the maximum we have considered, we have not tried to study the case of slightly (25% say) thicker hydrogen layers because our artificial starting evolutionary technique prevents us from being confident with the possibility that such structures can be actually be formed in Nature. The only way to answer this question is to perform full evolutionary calculations from the initial stages of evolution (hydrogen main sequence) to the WD stage.

Thus, performing full evolutionary studies coupled to oscillations is, in principle, interesting in two senses. First, it can provide significative improvements to the fitting of the oscillatory modes we have found for G117B-15A in this paper. Second, and more importantly, it can provide us with a stringent test of the quality of detailed theoretical evolutionary models of WD stars by comparing the theoretical pattern of non-radial oscillations against accurate observational determinations of the oscillatory pattern of WD stars. Work in this sense is in progress and will be published elsewhere.

## 7 acknowledgments

O.G.B. warmly thanks R. E. Martinez and E. Suarez for their help in getting some key information for the present paper.


format: text

text:
# 1 The parameter 𝐶 appearing in the exponential suppression factor of Voloshin as a function of the number of extra dimensions, 𝛿. The upper(lower) curve on the right-hand side of the plot is for the GT(DL) scheme.

Theories with extra dimensions and a low effective Planck scale($`M_{}`$) offer the exciting possibility that the production rate of black holes(BH) somewhat more massive than $`M_{}`$ can be quite large at future colliders. For example, cross sections of order 100 pb at the LHC, and even larger ones at the VLHC, have been advertised in the analyses presented by Giddings and Thomas(GT) and by Dimopoulos and Landsberg(DL). Although in practice the actual production cross section critically depends on the BH mass, the exact value of $`M_{}`$ and the number of extra dimensions, following the analysis of the authors in Ref., one finds very large rates over almost all of the interesting parameter space. These earlier analyses and discussions of the production of BH at colliders have been elaborated upon by several groups of authors and the production of BH by cosmic rays has also been considered. The most important question to address is whether or not the BH cross sections are actually this large or, at the very least, large enough to lead to visible rates at future colliders.

The basic idea behind the original collider BH papers is as follows: we consider the collision of two high energy Standard Model(SM) partons which are confined to a 3-brane, as they are in both the models of Arkani-Hamed, Dimopoulos and Dvali(ADD) and Randall and Sundrum(RS). In addition, we imagine that gravity is free to propagate in $`\delta `$ extra dimensions with the $`4+\delta `$ dimensional Planck scale assumed to be $`M_{}1`$ TeV. The curvature of the space is assumed to be small compared to the energy scales involved in the collision process so that quantum gravity effects can be neglected. When these partons have a center of mass energy in excess of $`M_{}`$ and the impact parameter for the collision is less than the Schwarzschild radius, $`R_S`$, associated with this center of mass energy, a $`4+\delta `$-dimensional BH is formed with reasonably high efficiency. It is expected that a very large fraction of the collision energy, i.e., $`\sqrt{s}`$, actually goes into the BH formation process so that $`M_{BH}\sqrt{s}`$. The subprocess cross section for the production of a non-spinning BH is thus essentially geometric for each pair of initial partons:

$$\widehat{\sigma }ϵ\pi R_S^2;$$
(1)
where $`ϵ`$ is an factor that accounts for finite impact parameter and angular momentum corrections and which is expected to be $`\begin{array}{c}<\hfill \\ \hfill \end{array}1`$. (We will assume that $`ϵ=1`$ in our calculations below.) We note that the $`4+\delta `$-dimensional Schwarzschild radius, $`R_S`$, scales as

$$R_S\left[\frac{M_{BH}}{M_{}^{2+\delta }}\right]^{\frac{1}{1+\delta }},$$
(2)
apart from an overall $`\delta `$\- and convention-dependent numerical prefactor. This convention dependence has been found to be somewhat confusing for unprepared readers of the literature and particularly to experimenters investigating the BH production phenomena. The confusion can be traced back to two different expressions for the $`4+\delta `$ dimensional Schwarzschild radius used by the original two sets of authors GT and DL. Explicitly, the two sets of authors employed distinct relationships between their $`4+\delta `$-dimensional Planck masses, $`M_{GT,DL}`$, and the associated generalized Newton’s constant, $`G_{4+\delta }`$: one finds that $`M_{DL}^{2+\delta }=G_{4+\delta }^1`$, while $`M_{GT}^{2+\delta }=(2\pi )^\delta /4\pi G_{4+\delta }`$. Depending on how the input parameters are chosen, i.e., which of the above masses is identified with $`M_{}`$, this scheme dependent prefactor can turn out to be relatively important. In particular, it is found to lead to an apparently very different $`\delta `$-dependence for the BH production cross section in the two cases. In the DL case the $`\delta `$-dependence of the numerical prefactor is rather weak whereas is it somewhat stronger in the GT analysis. Naively, for the same input value of $`M_{BH}`$ one finds the ratio of the cross sections obtained by the two sets of authors to be

$$\frac{\widehat{\sigma }_{GT}}{\widehat{\sigma }_{DL}}=\left[\frac{(2\pi )^\delta }{(4\pi )}\right]^{\frac{2}{1+\delta }}\left[\frac{M_{DL}^2}{M_{GT}^2}\right]^{\frac{2+\delta }{1+\delta }},$$
(3)
which is always greater than unity for $`\delta 2`$ and grows as $`\delta `$ increases if one assumes $`M_{}=M_{GT}=M_{DL}`$ as an input. (The later equality is, of course, invalid.) It is important to note, however, that when the differences in the definitions of the Planck scale are accounted for, both cross sections lead to the same numerical results. This can be accomplished by using a common value of $`G_{4+\delta }`$ as input in both cases. Written in terms of $`G_{4+\delta }`$ the above ratio of cross sections is unity. Unfortunately, since $`G_{4+\delta }`$ is a rather cumbersome object to use as an input parameter, in our discussion below we will display results for both sets of notation, i.e., we will identify $`M_{}`$ as either $`M_{DL}`$ or $`M_{GT}`$.

The approximate geometric subprocess cross section expression is claimed to hold by both GT and DL when the ratio $`M_{BH}/M_{}`$ is “large”, i.e., when the system can be treated semi-classically and quantum gravitational effects are small; one may debate just what “large” really means, but it most likely means “at least a few”. Certainly when $`M_{BH}/M_{}`$ is near unity one might expect curvature and stringy effects to become important and even the finite extent of the incoming partons associated with this stringy-ness would need to be considered. Clearly caution must be applied when $`M_{BH}M_{}`$ in interpreting cross sections evaluated in this parameter space region.

In order to obtain the actual cross section at a collider, one takes the geometric parton-level result, folds in the appropriate parton densities, sums over all pairs of possible partons and then integrates over the relevant kinematic variables. The resulting total cross section for BH with masses $`M_{BH}^{min}`$ is then given by the expression

$$\sigma =_{M_{BH}^{min2}/s}^1𝑑\tau _\tau ^1\frac{dx}{x}\underset{ab}{}f_a(x)f_b(\tau /x)\widehat{\sigma }(M_{BH}),$$
(4)
where we have summed over all possible pairs of initial state partons with their associated densities $`f_i(x)`$. These parton densities are evaluated at a scale $`Q^2=M_{BH}^2`$ in the numerical results presented below.

Voloshin has recently provided several arguments which suggest that an additional exponential suppression factor, $`S=e^{I_E}`$, must be included which seriously damps the pure geometric cross section for this process even in the semi-classical case. This implies we should rescale the purely geometric result as $`\widehat{\sigma }S\widehat{\sigma }`$ in the equations above. Here, $`I_E`$ is the Euclidean Gibbons-Hawking action for the BH generalized to the case of $`4+\delta `$ dimensions. Then $`S`$ is given explicitly by

$$S=exp[4\pi R_SM_{BH}/(1+\delta )(2+\delta )]=exp\left[C\left(\frac{M_{BH}}{M_{}}\right)^{\frac{2+\delta }{1+\delta }}\right],$$
(5)
where $`C`$ is a relatively small, though $`\delta `$-dependent, constant whose detailed form is also dependent upon whether the GT or DL schemes are followed. Fig. 1 shows the value of $`C`$ as a function of $`\delta `$ for both schemes. In the DL(GT) scheme, for example, when $`\delta =3`$ one obtains $`C=0.474(0.999)`$; similarly for $`\delta =5`$ one obtains $`C=0.233(0.706)`$. Note that in either scheme we see that $`C`$ is rather small and decreases substantially as $`\delta `$ is increased. Again, when the differences in notation in the two schemes are accounted for , i.e., using the same value of $`G_{4+\delta }`$, the value of the suppression factor is identical. The smallness of $`C`$ thus seriously reduces the numerical impact of the exponential suppression. From the $`\delta `$-dependence of $`C`$ and $`\delta `$-dependence of the power of the ratio $`M_{BH}/M_{}`$ in the exponential, it is clear that the Voloshin suppression will become less effective as the number of extra dimensions increases.

While the possibility of the existence of exponential suppression remains controversial, and strong arguments have been made on either side of the argument, for purposes of this discussion we will assume that it is indeed present. (However, we warn the reader that the jury is still out on this issue. In either case we anxiously await the resolution of this important argument.) If Voloshin’s criticisms of the geometrical cross section are valid one may certainly worry that the resulting exponentially suppressed rates for heavy BH production will possibly be too small to be observable at the LHC; as we will see below this need not be so. This is partly due to the smallness of $`C`$ and partly due to the expected huge size of the unsuppressed production rates for BH at the LHC.

Just how do the suppressed and unsuppressed cross sections at the LHC compare? As can be seen in Fig. 2 for the case $`M_{}=1`$ TeV, the unsuppressed rates for BH production at the LHC are quite large over a wide range of masses and numbers of extra dimensions using either set of authors’ cross section expressions. (It is important to remember that in this figure and the others below we appropriately identify $`M_{}`$ as either $`M_{GT}`$ or $`M_{DL}`$ depending on which set of predictions are under discussion.) Note that the results of Giddings and Thomas always appear larger than those of Dimopoulos and Landsberg due to the different definitions used for the Planck scale and that these apparent differences between the two sets of predictions increases as $`\delta `$ increases as was discussed above. This result is of course just a simple artifact of the known notational difference. We also see that Fig. 2 shows the effects of the suppression predicted by Voloshin in the two cases. From these results we make the important observation that for at least for some ranges of parameters BH will still be produced at rates that are large enough to be observable at the LHC even when the Voloshin suppression is active. For example, assuming that $`M_{BH}^{min}=5`$ TeV, we see that cross sections can easily be in the 100-1000 fb range. Although this is not a huge cross section the associated rates at the LHC will be quite large given an integrated luminosity of order 100 $`fb^1`$/yr. Note that the Voloshin suppression factor apparently modifies the two sets of predictions in quite different manners due to the two different expressions used for $`R_S`$. Since $`(R_S)_{GT}>(R_S)_{DL}`$ for all $`\delta 2`$, assuming the same input values for $`M_{GT}`$ and $`M_{DL}`$, the GT results are found to appear more suppressed than are those of DL. Again, if we take the notational differences in the two schemes into account, identical cross sections are obtained in both cases.

What happens as we vary $`M_{}`$? Figs. 3 and 4 show the effects of increasing $`M_{}`$ from 1 TeV to 2 and 3 TeV, respectively. As expected the unsuppressed rates for any fixed value of $`M_{BH}`$ decreases but we also see that the Voloshin suppression becomes less effective. This is also to be expected since the ratio $`M_{BH}/M_{}`$ in the exponent of the factor $`S`$ has been decreased for fixed $`M_{BH}`$. Again we see that for BH in the 5-6 TeV range it is relatively likely that the production cross section can quite easily be in excess of 100 fb.

We remind the reader that once produced these BH essentially decay semi-classically, mostly on the brane, via Hawking radiation into a reasonably large number $`25`$ or more final state partons with energies of order 50-100 GeV in a highly spherical pattern. Hadrons will be seen to dominate over leptons by a factor of order 5-10 for such final states. These unusual signatures would not be missed at either hadron or lepton colliders. (We note that an alternative decay scenario has been advocated by Casadio and Harms.) These features are sufficiently unique that BH production above conventional backgrounds should be observable at the LHC even if the cross sections are substantially smaller than the original estimates.

In summary, we have examined the production of BH at the LHC assuming that the exponential suppression of the geometric cross section predicted by Voloshin is realized. We have found that even when this suppression is significant the resulting BH production rates are still quite large for a wide range of model parameters given an integrated luminosity of order 100 $`fb^1`$. If the scale $`M_{}`$ is $``$ TeV then BH production should provide an exciting signature at the LHC and open a door to even more exciting physics.

Acknowledgements

The author would like to thank S. Thomas, G. Landsberg and J. Hewett for discussions related to this work.


text:
# Incommensurate Magnetism around Vortices and Impurities in High-𝑇_𝑐 Superconductors

## Abstract

By solving self-consistently an effective Hamiltonian including interactions for both antiferromagnetic spin-density wave (SDW) and $`d`$-wave superconducting (DSC) orderings, a comparison study is made for the local magnetic structure around superconducting vortices and unitary impurities. To represent the optimally doped regime of cuprates, the parameter values are chosen such that the DSC is dominant while the SDW is vanishingly small. We show that when vortices are introduced into the superconductor, an oscillating SDW is induced around them. The oscillation period of the SDW is microscopically found, consistent with experiments, to be eight lattice constants ($`8a_0`$). The associated charge-density wave (CDW) oscillates with a period of one half ($`4a_0`$) of the SDW. In the case of unitary impurities, we find a SDW modulation with identical periodicity, however without an associated CDW. We propose neutron scattering experiments to test this prediction.

The common feature of all high-transition-temperature (high-$`T_c`$) cuprates is the proximity between antiferromagnetic (AF) and d-wave superconducting (DSC) phases controlled by the doping—At low temperatures, the AF order of the parent compounds is replaced by a DSC order upon doping. The main theme of high-$`T_c`$ superconducivity research centers on how to establish the connection between these two kinds of orderings, i.e., whether they exclude each other, or coexist microscopically. In one of the strategies, the interplay between the AF and DSC ordering is explored by weakening the superconductivity in the optimally or slightly overdoped regime. The inelastic neutron scattering (INS) measurements on both optimally doped and underdoped La2-xSrxCuO4 (LSCO) samples by Lake and co-workers provide the first evidence of the field-induced magnetic order Lake01 ; Lake02 . The field-induced spin fluctuations have a spatial periodicity of $`8a_0`$ with the wave vector pointing along the Cu-O bond directions. The principle magnetization oscillations in the vortex state are coherent over a distance $`L_M>20a_0`$, substantially longer than the superconducting coherence length $`\xi _0`$ ($`5a_0`$). The field-induced enhancement of the Bragg peak intensity was also observed by Khaykovich et al. Khaykovich01  in the elastic neutron scattering (ENS) measurement on a related material, La2CuO4+y (LCO), but with $`L_M>100a_0`$, indicating the existence of the field-induced static AF order of $`8a_0`$ periodicity. Strong AF fluctuations have also been observed by Mitrovic et al. Mitrovic01  in a high-field nuclear magnetic resonance (NMR) imaging experiment on near-optimally doped YBa2Cu3O7-x (YBCO). More recently, the scanning tunneling microscopy imaging by Hoffman et al. Hoffman02  has revealed the quasiparticle states around the vortex cores in slightly overdoped Bi2Sr2CaCu2O8+x (BSCCO) a Cu-O bond-oriented “checkboard” pattern with $`4a_0`$ periodicity. The periodicity ($`4a_0`$) of charge modulation is one half of that ($`8a_0`$) of the field-induced SDW modulation. Theoretically, Demler and co-workers, by focusing their attention on the far-field regions outside the vortex core, proposed a phenomenological model Demler01  that when the superconductivity is weakened by the circulating currents caused by the vortices, a coexisting SDW plus DSC phase appears surrounding the core; more recently they also argued Zhang01  that the halved periodicity of the static CDW modulation is associated with the “Friedel oscillation of the spin gap”. As an extension of earlier work within the SO(5) theory Zhang97 ; Arovas97 , originally suggesting an AF insulating region inside the core, Hu and Zhang Hu01  showed that the vortex-induced AF region can be greater than $`\xi _0`$, due to the light effective mass of the dynamic AF fluctuations at optimal doping. In these two cases, the $`8a_0`$ periodicity of the magnetic modulation is not fully understood due to the phenomenological nature of the models. Distinctly, stripe models Poilblanc89 ; Zaanen89 ; Emery90 ; Zachar98  predict that the spin modulation of wavelength $`\lambda `$ in cuprates should be associated with the charge modulation of wavelength $`\lambda /2`$. However, two dimensional modulation is apparent in the STM image Hoffman02 . One of us and C.S. Ting Zhu01  have applied an effective microscopic mean-field model accounting for the competition between the AF and DSC orderings Martin00  to successfully explain the peak-split structure around the Fermi surface in the local density of states (LDOS) at the vortex core center Maggio95 ; Pan00 ; Hoogenboom01 . Nevertheless, due to the use of the simplest band structure parameter values, the spin and charge modulation is too weak for the determination of the periodicity. In this work, using realistic model parameter values, we present a comprehensive study of the spin/charge structure around the vortices and nonmagnetic unitary impurities in optimally doped high-$`T_c`$ superconductors. We find, for the first time, microscopically: (i) For the vortex case, the oscillation periods of the SDW and the associated CDW are indeed $`8a_0`$ and $`4a_0`$, respectively, in good agreement with the INS Lake01  and STM Hoffman02  experiments. (ii) Around the nonmagnetic unitary impurity, the modulation of the SDW order still exhibits $`8a_0`$ periodicity. However, the charge density shows only the Friedel oscillation with a period of the Fermi wavelength due to the strong potential scattering from the impurity.

Consider a minimal model defined on a two dimensional (2D) lattice square , in which the on-site repulsion is solely responsible for the antiferromagnetism while the nearest neighbor attraction causes the $`d`$-wave superconductivity Martin00 . With the application of an external magnetic field and/or in the presence of nonmagnetic impurities, the effective mean-field Hamiltonian can be diagonalized by solving self-consistently the Bogoliubov-de Gennes equation Zhu01 :

$$\underset{j}{}\left(\begin{array}{cc}_{ij,\sigma }& \mathrm{\Delta }_{ij}\\ \mathrm{\Delta }_{ij}^{}& _{ij,\overline{\sigma }}^{}\end{array}\right)\left(\begin{array}{c}u_{j\sigma }^n\\ v_{j\overline{\sigma }}^n\end{array}\right)=E_n\left(\begin{array}{c}u_{i\sigma }^n\\ v_{i\overline{\sigma }}^n\end{array}\right).$$
(1)
Here $`(u_{i\sigma }^n,v_{i\overline{\sigma }}^n)`$ is the quasiparticle wavefunction corresponding to the eigenvalue $`E_n`$, the single particle Hamiltonian $`_{ij,\sigma }=t_{ij}e^{i\phi _{ij}}+(m_{i\overline{\sigma }}+ϵ_i\mu )\delta _{ij}`$, where $`t_{ij}=t`$ for the nearest neighbor hopping while $`t_{ij}=t^{}`$ for the next-nearest neighbor hopping, $`ϵ_i`$ is the single site potential describing the scattering from impurities, and the Peierls phase factor $`\phi _{ij}=\frac{\pi }{\mathrm{\Phi }_0}_{𝐫_j}^{𝐫_i}𝐀(𝐫)𝑑𝐫`$ with $`\mathrm{\Phi }_0=hc/2e`$ the superconducting flux quantum in the presence of an externally applied magnetic field. Notice that the quasiparticle energy is measured with respect to the Fermi energy. The self-consistency conditions read:

$`m_i`$ $`=`$ $`Un_i=U{\displaystyle \underset{n}{}}|u_i^n|^2f(E_n),`$ (2a)
$`m_i`$ $`=`$ $`Un_i=U{\displaystyle \underset{n}{}}|v_i^n|^2[1f(E_n)],`$ (2b)
and

$$\mathrm{\Delta }_{ij}=\frac{V}{4}\underset{n}{}(u_i^nv_j^n+v_i^nu_j^n)\mathrm{tanh}\left(\frac{E_n}{2k_BT}\right),$$
(3)
where $`U,V`$ are the strength of the on-site repulsion and the nearest neighbor attraction, respectively, and the Fermi distribution function $`f(E)=1/[e^{E/k_BT}+1]`$. Here the summation is over the eigenstates with both positive and negative eigenvalues Zhu01 . We report results below for two cases at zero temperature. For the Abrikosov vortex state, the magnetic field effects enter through the Peierls phase factor $`\phi _{ij}`$ and no impurities are introduced ($`ϵ_i=0`$). For the effects of a single impurity, we set $`𝐀=0`$ so that $`\phi _{ij}=0`$. Hereafter we measure the length in units of the lattice constant $`a_0`$ and the energy in units of the hopping integral $`t`$. To mimic a hole-like Fermi surface, as relevant to the hole-doped cuprates, we take $`t^{}=0.2`$. As a model calculation, the on-site repulsion and pairing interactions are taken to be, $`U=2.5`$ and $`V=1.0`$. In addition, the filling factor $`n_f=_{i,\sigma }n_{i\sigma }/N_xN_y`$ is fixed to be 0.84, which corresponds to an optimal hole doping $`n_h=1n_f=0.16`$. Here $`N_x,N_y`$ are the linear dimensions of the unit cell under consideration. We use an exact diagonalization method to solve the BdG equation (1) self-consistently. In the absence of magnetic field and impurities, we recover all results reported in Martin00 . In the optimal doping regime ($`n_h=0.16`$), the AF SDW order is absent while the DSC order is homogeneous. When the DSC is weakened by the application of an external magnetic field or by introducing impurities, the AF SDW will be nucleated in the region where the DSC order is depressed.

Field-induced SDW and CDW. When an external magnetic field is applied perpendicular to the 2D Cu-O plane, $`𝐇=H\widehat{𝐳}`$ ($`H_{c1}HH_{c2}`$), an Abrikosov vortex state is formed. As the vortex core is approached, the DSC order parameter vanishes topologically due to the circulating supercurrent surrounding the core. It is assumed that the superconductor is in the extreme type-II limit where the Ginzburg-Landau parameter $`\kappa =\lambda _0/\xi _0`$ goes to infinity so that the screening effect from the supercurrent is negligible. We choose a Landau gauge to write the vector potential as $`𝐀=(Hy,0,0)`$, where $`y`$ is the $`y`$-component of the position vector r, so that the Peierls phase factor $`\phi _{ij}`$ is uniquely determined. By taking the strength of magnetic field, $`H=2\mathrm{\Phi }_0/N_xN_y`$, such that the flux enclosed by each unit cell is twice $`\mathrm{\Phi }_0`$, we solve self-consistently the BdG equation (1) with the aid of the magnetic Bloch theorem as given by Eq. (3) in Ref. Zhu01 . For the calculation, we consider the magnetic unit cell of size $`N_x\times N_y=48\times 24`$, which gives a square vortex lattice. The numerics shows that each unit cell accommodates two superconducting vortices each carrying a flux quantum $`\mathrm{\Phi }_0`$, which conforms to the above prescription for the magnetic field strength. The unit cell is equally partitioned between these two vortices, each located at the center of area $`\frac{N_x}{2}\times N_y`$ sites. Typical results on the structure around one vortex core are displayed in Fig. 1, where the left column is three-dimension plots and the right column is contour plots. As shown in Fig. 1(a), the DSC order parameter vanishes at the core center (dark-blue site) and approaches its zero-field value, which is about $`\mathrm{\Delta }_0=0.08`$ for the chosen parameter values, away from the core center. Notice that the DSC order parameter is not uniform beyond the distance $`\xi _0`$ away from the core center. Instead, it is weakly modulated. The maximum modulation amplitude is less than $`0.05\mathrm{\Delta }_0`$. This modulation is closely related to the appearance of the field-induced SDW, as will be discussed immediately. Fig. 1(b) displays the spatial distribution of the staggered magnetization of the local SDW order defined as $`M_s=(1)^i(n_in_i)`$. Clearly, the maximum strength of $`M_s`$ is pinned at the vortex core center. This AF SDW order exhibits a modulation pattern with the satellite peaks (yellow spots) and valleys (dark-blue spots) regularly spaced throughout the unit cell along the Cu-O bond directions. The modulation pattern of the SDW implies a much longer magnetic correlation length as compared to $`\xi _0`$, which is consistent with the INS measurements Lake01 . The appearance of the SDW order around the vortex core also strongly affects the electron density $`n_i=_\sigma n_{i\sigma }`$. As shown in Fig. 1(c), at the vortex core center, where the SDW amplitude reaches the global maximum, the electron density is strongly enhanced. In addition, the charge density also exhibits regular modulation. By comparing the spatial distribution of the DSC order parameter, the field-induced SDW as well as the associated CDW, one finds that as long as the absolute amplitude ($`|M_s|`$) of SDW reaches a local maximum (both yellow and dark-blue spots in Fig. 1(b)), the associated CDW also reaches a local maximum (green spots in Fig. 1(c)), that is, a local minimum in the hole density, while the DSC order parameter has a local minimum (shallow red spots in Fig. 1(a)). The one to one correspondence between the SDW and associated CDW orderings at zero field, was discussed in the context of stripes Poilblanc89 ; Zaanen89 ; Emery90 ; Zachar98 ; Martin00 . To make a closer inspection of the periodicity of the SDW and associated CDW modulation, we perform a Fourier transform of these two quantities. As shown in Fig. 2(a), the strongest spectral intensity of the SDW modulation occurs at the wave vectors $`𝐤=2\pi (\frac{3}{24},0)`$ and $`𝐤=2\pi (0,\frac{3}{24})`$, which gives unambiguously the period of $`8a_0`$ along the Cu-O bond directions. To verify this point, we have performed a calculation on the unit cell of different size $`N_x\times N_y=52\times 26`$, and found that the period of the SDW modulation remains $`8a_0`$, indicating convincingly that the periodicity of the field-induced SDW modulation is an intrinsic property. Therefore, our result explains very well the INS measured value of the SDW modulation Lake01 . We notice that, in the optimal doping regime, the circulating supercurrent around the vortex, which favors the isotropy between $`x`$ and $`y`$ directions, overcomes the one-dimensionality of stripes. In the underdoping regime, the stripe behavior becomes dominant Ichioka01 .

Impurity-induced SDW and CDW. The $`d`$-wave superconductivity can also be suppressed locally by a strong nonmagnetic impurity. The effects of impurities on superconductors have been of theoretical and experimental interest even in its own right for a long time. In view of the recent observations of the AF magnetism around the vortex core, it is important to see whether a similar magnetic structure can also be induced around an impurity. We perform the numerical calculation on a unit cell of $`N_x\times N_y=36\times 36`$ sites. The single-site potential is taken to be $`ϵ_0=100`$ at the unitary impurity and zero on other sites. Since $`\phi _{ij}=0`$ in this case, the conventional Bloch theorem is used. The results on the spatial distribution of the three orderings is displayed in Fig. 3. The DSC order parameter is depressed dramatically at the impurity site and approaches the bulk value at the scale $`\xi _0`$. Out of this range, almost no modulation of the DSC order parameter is seen. The induced staggered moment of the SDW is zero at the impurity site and has maxima on the four nearest neighbor Cu sites of the impurity. Away from the impurity, the induced SDW shows a modulation very similar to the vortex case, although with a much weaker amplitude. However, the electron density, which is zero at the impurity, only exhibits a Friedel-like oscillation within a limited range around the impurity. The Fourier transform, as shown in Fig. 4, gives a clear picture: The main modulation period of the SDW is $`8a_0`$; the CDW is modulated with a period of $`2a_0`$, which is roughly the Fermi wave length. Therefore, the correspondence between the SDW and CDW is absent in this impurity case. The reason lies in the fact that the CDW is very sensitive to the strong impurity scattering. This result is consistent with the STM measurements in BSCCO Pan00b  that show no evidence for CDW modulation around the impurities. However, SDW appears to be a robust feature induced by the impurity, which should be observable in neutron scattering experiments.

Nuclear magnetic resonance (NMR) measurements have shown that when a Cu2+ in the Cu-O plane is substituted by a strong nonmagnetic impurity, such as Zn2+, an effective magnetic moment can be induced on the Cu sites around the impurity site Alloul91 ; Ishida96 ; Julien00 ; Bobroff01 . More recent NMR measurements Bobroff01  show for the first time that near optimal doping, the Kondo screening effect observed above the superconducting transition temperature, is strongly reduced in the superconducting state. This indicates the stabilization of the magnetic moments. Our result of the strongest staggered magnetic moments on the four nearest-neighbor Cu sites of the impurity is consistent with the above experiment Bobroff01 . This leads us to speculate that the magnetism around vortices and around impurities may share a common origin.

In conclusion, we have studied the magnetism around vortices and nonmagnetic unitary impurities. In the case vortices, the experimentally observed $`8a_0`$ period of the SDW modulation is explained for the first time based on a microscopic model. The correspondence between the SDW and CDW modulations has also been established. Around the unitary impurity, we have also shown the existence of the SDW modulation with a period of $`8a_0`$. The existence of such modulation can be tested by neutron scattering experiments.

Acknowledgments: We wish to thank A.V. Balatsky, J.C. Davis, B. Khaykovich, G. Ortiz, S. Sachdev, C.S. Ting for useful discussions. This work was supported by the US Department of Energy through the Los Alamos National Laboratory.