text:
# Irradiated Planets

## 1 Introduction

Many indirect detections of substellar objects in close orbit around stars of spectral type later than F have been made since 1996. However, currently there are few observations that can be useful for constraining the large parameter space (e.g. chemical composition, albedo, age, and inclination) for any of these so called extra-solar giant planets (EGPs). The transits observed for HD209458a by Henry et al. (2000) and Charbonneau et al. (2000) have helped constrain the gross physical parameters of HD209458a (radius, mass, mean density, etc.) but not any of the atmospheric properties. There is great hope that the handful of ambitious space and ground based projects planned for the next decade will be capable of making direct photometric observations of EGPs and possibly measurements of EGP spectral features. Until then, we must rely on atmospheric modeling to provide insight into the basic properties of these objects to guide the observers while we wait patiently for their observations.

An early work by Saumon et al. (1996) investigated the properties of EGPs for various masses and ages near primaries of different spectral type but approximated the reflected and thermal flux as gray bodies. More recently there have been several radiative equilibrium models produced for the purpose of predicting certain observables for 51 Pegasi B Seager & Sasselov (1998) and $`\tau `$ Boo Goukenleuque et al. (2000). Also, a broad range (100 – 1700$`\mathrm{K}`$) of EGP models were studied and loosely classified by Sudarsky et al. (2000) using ad hoc temperature – pressure profiles. In this paper, we present equilibrium models for several scenarios well within the known parameter space for EGPs. We have investigated the variations of the thermal structures and emergent flux as functions of both the spectral type of the primary and the orbital separation. We also address the importance of spherical versus plane parallel radiative transfer. Our emphasis in this paper is not on any one particular EGP, but instead on the basic understanding of these objects and the effects of their close proximity to a stellar companion.

## 2 Model Construction

We have used our multi-purpose atmosphere code PHOENIX (version 10.9) to generate the models discussed below. Most details of the radiative transfer method may be found in Hauschildt & Baron (1999), but, for clarity, we repeat some of the basic features and discuss a few changes needed for irradiated models. PHOENIX solves either the full spherically symmetric or plane parallel radiative transfer equation (PPRTE) using an operator splitting (OS) technique. For the majority of the calculations presented here, we have chosen plane parallel geometry and assumed hydrostatic equilibrium. However, we will discuss the importance of spherical geometry in section 4. All models are subject to an energy conservation constraint such that the total Flux (convective and radiative) is constant at each layer. Each model atmosphere spans a range of optical depth ($`\tau _{std}`$ defined at $`1.2\mathrm{\mu m}`$) from $`0`$ at the top of the atmosphere down to $`100`$ at the deepest layer. Convection is treated according to the Mixing Length Theory from the onset of the Schwarzschild criterion with mixing length parameter, $`\alpha =1`$.

The usual boundary conditions for an isolated star are that the inward directed flux at the surface should be zero ($`\mathrm{I}_\nu ^{}(\tau _{std}=0,\mu )=0`$, where $`1\mu =\mathrm{cos}(\theta )0`$) and that the diffusion approximation holds at the bottom of the atmosphere. For a close binary, the situation is clearly different. At the surface of the secondary, the boundary condition on $`\mathrm{I}_\nu ^{}`$ is determined by the incident flux ($`\mathrm{F}_\nu ^{\mathrm{inc}}`$) given by:

$$2\pi \underset{1}{\overset{0}{}}\mathrm{I}_\nu ^{}(\mu )\mu 𝑑\mu =\mathrm{F}_\nu ^{\mathrm{inc}}(\tau _{std}=0)$$
(1)
where

$$\mathrm{F}_\nu ^{\mathrm{inc}}(\tau _{std}=0)=\left(\frac{R^{}}{a}\right)^2\mathrm{F}_\nu ^{}$$
(2)
In the equations above, $`\mathrm{I}_\nu ^{}(\mu )`$ refers to the inward directed intensities along direction $`\mu `$, $`R^{}`$ is the radius of the primary, a is the surface to surface primary-secondary separation, and $`\mathrm{F}_\nu ^{}`$ is the monochromatic flux from the primary. For $`\mathrm{F}_\nu ^{}`$, we use a synthetic spectrum taken from a previous PHOENIX calculation Hauschildt et al. (1999); Allard et al. (2000). For the models presented below, we have made the simplifying assumption that the impinging radiation field is isotropic, meaning that $`\mathrm{I}_\nu ^{}(\mu )`$ at the surface is the same for all $`\mu `$ (i.e., $`\mathrm{I}_\nu ^{}(\mu )=\mathrm{I}_\nu ^{}`$). We have also assumed that the flux is not globally redistributed over the planet’s surface. A more detailed discussion of these last two issues may be found in section 4.

As was done in the work on irradiated M Dwarfs by Brett & Smith (1993), all of the incident radiation from the primary is re-radiated outward by the secondary in the form of reflected flux ($`\mathrm{F}^{\mathrm{ref}}`$) and as a contribution to the thermal flux ($`\mathrm{F}^{\mathrm{therm}}`$). This constraint is required by energy conservation and implies that the integrated flux at the surface is equal to $`\sigma \mathrm{T}_{\mathrm{int}}^4+\mathrm{F}^{\mathrm{inc}}`$. Throughout this paper, $`\mathrm{T}_{\mathrm{int}}`$ refers to the effective temperature of the planet in the absence of irradiation and $`4\pi \mathrm{R}_{\mathrm{p}}^{}{}_{}{}^{2}\sigma \mathrm{T}_{\mathrm{int}}^4`$ equals the planet’s intrinsic luminosity where $`\mathrm{R}_\mathrm{p}`$ is the planet’s radius. The intrinsic luminosity is an age dependent quantity which represents the energy released by the planet as it cools and contracts. $`\mathrm{T}_{\mathrm{int}}`$ also relates irradiated planets to isolated planets in which case $`\mathrm{T}_{\mathrm{int}}`$ is identical to the more commonly used $`\mathrm{T}_{\mathrm{eff}}`$.

The “effective temperature” ($`\mathrm{T}_{\mathrm{eff}}`$), which is customarily defined as the temperature a black body would have to have in order to radiate the total flux, has important physical and observational significance for isolated stellar and substellar objects. However, for irradiated planets (and stars) $`\mathrm{T}_{\mathrm{eff}}`$ loses some of its connection to the fundamental properties of the planet because it is difficult to separate, by observation, those photons which are thermally radiated by the planet from those which originated from the primary and are merely reflected by the planet. We shall only use $`\mathrm{T}_{\mathrm{eff}}`$, with its customary definition, to describe non-irradiated objects. When describing irradiated objects, we will define another quantity which describes the equilibrium temperature of the planet’s day side;

$$\sigma \mathrm{T}_{\mathrm{eq}}^4=\sigma \mathrm{T}_{\mathrm{int}}^4+(1A_B)\mathrm{F}^{\mathrm{inc}},$$
(3)
where $`A_B`$ is the Bond albedo. Since the reflected flux is not directly related to the equilibrium thermal structure of the atmosphere, it has been omitted from equation 3 and, therefore, distinguishes $`\mathrm{T}_{\mathrm{eq}}`$ from $`\mathrm{T}_{\mathrm{eff}}`$. It is also important to realize that $`\mathrm{T}_{\mathrm{eq}}`$, as defined above, represents the equilibrium state at a given age and allows for the possibility that the intrinsic luminosity has not reached zero. $`\mathrm{T}_{\mathrm{int}}`$ will only be important for young (or more massive) planets when the primary is a solar type star. However, for planets orbiting M dwarfs, $`\mathrm{T}_{\mathrm{int}}`$ can be a significant contribution to $`\mathrm{T}_{\mathrm{eq}}`$. While $`\mathrm{T}_{\mathrm{eq}}`$, which likely varies across the planet’s day side, tells us little about the interior or the non-irradiated face of the planet it does provide a useful measure of the effects on the thermal structure of the planet’s day side. It is also important to realize that $`\mathrm{T}_{\mathrm{int}}`$ is not necessarily the equilibrium temperature of the planet’s night side since energy may be carried over from the day side.

To produce an irradiated atmospheric model, we chose the structure of a converged non-irradiated model taken from Allard et al. (2001) as the initial structure of an irradiated atmosphere located $`5\mathrm{A}\mathrm{U}`$ from the primary. At this distance, the impinging radiation produces only small changes in the structure of the upper layers and convergence is achieved after only a few iterations. This new model is then used as the initial structure for an irradiated atmosphere located at $`2.5\mathrm{AU}`$ from the primary. This process of moving the planet closer to the primary is repeated until the desired orbital separation is reached. Each intermediate calculation is iterated until the changes to the temperature structure are less than 1 $`\mathrm{K}`$ at every depth point and energy conservation is satisfied to within a few percent. Models obtained in this way will have structures, at the deepest layers, similar to the initial non-irradiated model chosen from Allard et al. (2001), and thus depend on $`\mathrm{log}(g)`$, and the assumed value of $`\mathrm{T}_{\mathrm{int}}`$.

The opacity setups used for the EGPs are identical to the “AMES-Cond” and “AMES-Dusty” models of Allard et al. (2001) which refer to two limiting cases. AMES-Cond refers to the situation when dust forms in the atmosphere at locations determined by the chemical equilibrium equations, but has been entirely removed from the atmosphere by efficient gravitational settling. Dust formation, therefore, acts only to remove refractory elements and reduce the number of certain molecules but does not contribute to the overall opacity. The opposite case, AMES-Dusty, ignores settling altogether. Dust forms based on the same criteria as in the AMES-Cond models yet remains present to contribute to the opacity. The Cond model may be thought of as “clear skies”, and the Dusty model as “cloudy skies”.

The opacities include H2O and TiO lines by Partridge & Schwenke (1997), CH4 lines from the HITRAN and GEISA databases Rothman et al. (1992); Husson et al. (1992), and H2, N2, Ar, CH4, and CO2 collision induced absorption (CIA) opacities according to (Borysow et al., 1997a, b; Borysow & Frommhold, 1986a, 1987a, b; Borysow & Tang, 1993; Samuelson et al., 1997; Borysow & Frommhold, 1986c, 1987b; Gruszka & Borysow, 1997, and references therein). To include dust grains, we have assumed an interstellar size distribution with diameters ranging from 0.00625 to 0.24$`\mathrm{\mu m}`$ and the chemical equilibrium equations incorporate over 1000 liquids and crystals Allard et al. (2001). We follow the prescriptions of Grossman (1972) and use Gibbs free energies of formation from the JANAF database Chase et al. (1985) to compute the number densities for each grain species. The condensation equilibrium is directly incorporated into the chemical equilibrium equations to account, self-consistently, for the depletion of refractory elements as a function of gas temperature and pressure. With this approach, any changes to the structure brought on by irradiation will be automatically accounted for in the chemical equilibrium equations.

## 3 Results

We have divided our study into two separate cases: one for which the primary star is a dM5 ($`\mathrm{T}_{\mathrm{eff}}=3000\mathrm{K}`$, see Leggett et al. (2000); Allard et al. (2000)) and another where the primary is a G2 solar type star ($`\mathrm{T}_{\mathrm{eff}}=5600\mathrm{K}`$, see Hauschildt et al. (1999)). For the former case, one finds few observed objects in the literature (eg Gl876B). However, surveys for additional objects have begun Delfosse et al. (1999) . In the last decade, most of the observations have been spent searching for objects that fall into the second group, and, to date, there are roughly 50 planets known to be in orbit around solar type stars. For both of these cases, we consider different effective temperatures and several orbital separations. Below, we present the structures and spectra for each case and discuss the effects caused by the impinging radiation.

### 3.1 EGPs Around M Dwarfs

Gl876B was the first planetary companion found orbiting an M dwarf. This was an important discovery given that M dwarfs constitute nearly 70% of all stars in the galaxy Henry (1998). Gl876B is also one of the nearest EGPs (only 5 pc away) making it a promising candidate for direct imaging by future adaptive optics and interferometry missions Marcy et al. (1998). Our first group of irradiated models is intended to represent objects, like Gl876B, which have cool stellar primaries. We begin with a planet having $`\mathrm{T}_{\mathrm{int}}=500\mathrm{K}`$, $`\mathrm{R}_\mathrm{p}=1\mathrm{R}_{\mathrm{jup}}`$ and $`\mathrm{log}(g)=3.5`$ located very close to a dM5. Both planet and primary have solar compositions.

In Figure Irradiated Planets, we show the temperature structure of the non-irradiated planet compared to irradiated thermal structures at 0.1 and 0.05AU for both the AMES-Cond and AMES-Dusty cases. The outermost regions are significantly altered by the radiation from the dM5 resulting in a generally flatter temperature profile. Without dust opacity (the AMES-Cond model), the impinging radiation is capable of heating the lower layers of the planet thereby reducing the temperature gradient. Above the photosphere ($`\tau _{std}<10^4`$), the largest temperature increase occurs and a very slight temperature inversion forms at 0.05AU. Dust opacity generally produces a hotter atmosphere at all depths with a smoother spectral energy distribution. For more details on the effects of dust in non-irradiated atmospheres, see Allard et al. (2001). It is important to stress that the solution of the chemical equilibrium equations is based on the final temperature structure and thus incorporates the irradiation effects.

In the AMES-Dusty models, Fe, Mg2SiO4, MgSiO3 and CaMgSi2O6 and MgAl2O4 are the dominant dust species which form a cloudy region extending from roughly $`\tau _{std}=1.0`$ to the top of the atmosphere. Below this region ($`\tau _{std}>1.0`$) is a complicated mixture of various other condensates. In Figure Irradiated Planets, the concentrations of the dominant dust species are shown for the irradiated and non-irradiated planets. There is little change in the abundance of the condensates brought on by the impinging radiation from the dM5.

The spectra for the irradiated and non-irradiated planets are shown in Figure Irradiated Planets. The isolated AMES-Dusty planet produces a very smooth spectrum in the optical as is expected for an atmosphere dominated by grain opacity. In the IR, several distinct molecular bands are present: CH4 at 1.6$`\mathrm{\mu m}`$, 2.1 – 2.5$`\mathrm{\mu m}`$, and 3.0 – 4.0$`\mathrm{\mu m}`$, and H2O between 2.5 and 3.0$`\mathrm{\mu m}`$. The reflective properties of the dust between the near-UV and the near-IR is clearly shown in the spectrum of the irradiated AMES-Dusty planet. Even at 0.1AU, the planet reflects a considerable amount of radiation in the optical bands compared to the non-irradiated planet. However, the ratio ($`ϵ`$) of the planet flux to that of the incident flux is quite small. Averaged over 4500 to 5200Å, $`ϵ`$ is about $`10^7`$ for the AMES-Dusty at 0.05AU. In the IR, the flux also increases and the molecular bands become progressively shallower as the planet approaches the primary. The increase in flux at infrared wavelengths is not due to reflection, but is entirely a thermal effect.

The effects of irradiation on the spectrum of the 500$`\mathrm{K}`$ AMES-Cond planet are also shown in Figure Irradiated Planets. Distinct molecular bands, primarily due to water and methane, are clearly visible for $`\lambda >1\mathrm{\mu m}`$. These bands are sensitive to the temperature at various depths. The methane absorption (at $`1.8\mathrm{\mu m}`$, 2.2 – 3.0$`\mathrm{\mu m}`$ and 3.1 – 4.0 $`\mathrm{\mu m}`$) probes the deeper layers and indicates that, up to 0.1AU, little change in the temperature has taken place in the photosphere. Once 0.05AU separation has been reached, the irradiation has produced a much larger temperature rise resulting in nearly an order of magnitude increase in the emerging flux at 3.5$`\mathrm{\mu m}`$. Similarly, the water bands (at $`0.9\mathrm{\mu m},1.1\mathrm{\mu m}`$ and $`1.4\mathrm{\mu m}`$) probe the upper layers and are sensitive to the irradiation even at 0.1AU. The AMES-Cond model also produces a reflected component in the UV and optical regions, though less dramatic than in the Dusty model when compared to the non-irradiated planet’s intrinsic flux. In this case, $`ϵ=3.6\times 10^6`$ at 0.05AU, which is roughly 10 times that seen in the Dusty model.

It is quite possible that many of the young EGPs have $`\mathrm{T}_{\mathrm{int}}`$ closer to that of an L or T dwarf. To investigate this possibility, we have modeled an EGP with $`\mathrm{T}_{\mathrm{int}}=1000\mathrm{K}`$ (all other parameters are the same as for the above $`500\mathrm{K}`$ planet) irradiated by the same dM5 used above. The results are qualitatively similar but weaker, for both the Cond and Dusty situations, than those of the previous cooler planets. So, in order to explore the regime where irradiation effects become important, we must study these planets at even smaller orbital separations. Figures Irradiated Planets and Irradiated Planets demonstrate the effects of irradiation on the structure and spectrum, respectively, for the 1000$`\mathrm{K}`$ AMES-Dusty and AMES-Cond models. As with the cooler dusty planet, the inner regions are not affected by the irradiation even at 0.005 AU ($`1\mathrm{R}_{}`$), a separation more common in cataclysmic variables. The AMES-Dusty spectrum (Fig. Irradiated Planets) is featureless in the IR for separations $`<0.01\mathrm{AU}`$, but many of the absorption lines present in the dM5 are now reflected in the optical spectrum of the planet. The AMES-Cond planet experiences a significant temperature increase at all photospheric layers (Fig. Irradiated Planets) and a suppression of the convective zone also seen in the cooler AMES-Cond planet. The spectra (Fig. Irradiated Planets) show several orders of magnitude increase in the flux for the methane and water bands at 0.01AU with many reflected atomic features in the optical when the Cond planet is at 0.005AU from the primary.

### 3.2 EGPs Around Solar Type Stars

Many EGPs have been discovered orbiting solar type stars (e.g., 51 Peg B), and have already received much attention in the literature in both observational and theoretical studies. Recent works include models of 51 Peg B Goukenleuque et al. (2000), HD209458a Seager & Sasselov (2000) and a study of EGP photometric and polarization properties Seager et al. (2000). We have taken a slightly more general approach; we do not focus on any one particular known object. We have chosen, instead, a planet with $`\mathrm{T}_{\mathrm{int}}=500\mathrm{K}`$ and a young planet with $`\mathrm{T}_{\mathrm{int}}=1000\mathrm{K}`$ orbiting an object similar to the Sun as average representatives of most currently known EGPs orbiting F, G, and K stars.

Figure Irradiated Planets demonstrates the effects of impinging radiation from a G2 primary on the $`500\mathrm{K}`$ AMES-Cond and AMES-Dusty atmospheric structure for several separations. As in the dusty cases presented above, the inner part of the AMES-Dusty atmosphere is unaffected by the irradiation. However, regions around $`\tau _{std}10^4`$ show a dramatic rise in temperature as the planet is moved from 1.0 to 0.3AU. At 0.3AU, the Dusty model reaches $`\mathrm{T}_{\mathrm{eq}}=767\mathrm{K}`$ and the temperature in the upper atmosphere ($`\tau _{std}<10^4`$) has more than doubled. Also at 0.3AU, the AMES-Cond model reaches $`\mathrm{T}_{\mathrm{eq}}=735\mathrm{K}`$ and, for regions near $`\tau _{std}=1`$, the temperature has risen by over 150$`\mathrm{K}`$. Even the deepest layers feel the presence of the primary. The temperature at $`\tau _{std}=10`$ increases by $`100\mathrm{K}`$, and the boundary between the radiative and convective zone, which was well above $`\tau _{std}=1`$ at 1.0AU, has retreated to nearly $`\tau _{std}=50`$. Large concentrations of dust species are still capable of forming in the upper atmosphere. In the AMES-Dusty model, at 0.15AU, there still remains a thick cloudy region with concentrations similar to those seen in the $`500\mathrm{K}`$ planet near a dM5 (see Fig. Irradiated Planets). The only major difference being Fe which extends to the upper regions where $`\mathrm{log}(\mathrm{P}_{\mathrm{gas}})1.0`$ dynes cm-2.

The optical and IR regions of the spectrum for each orbital separation are shown in Figure Irradiated Planets where the differences between the AMES-Dusty and AMES-Cond models can clearly be seen. The dusty models produce very smooth spectra except for the reflected features in the optical bands. The dust is entirely responsible for the reflection effects in this case. As the dusty planet is brought closer to the G2, the few absorption features (primarily CH4) seen in the non-irradiated planet completely disappear. At 0.5AU, very little change has taken place for $`\lambda >7000\text{Å}`$ in the Cond model spectrum and the water and methane bands remain strong. However, for $`\lambda 7000\text{Å}`$, a large amount of reflected light is already present even at 0.5AU. As the planet is brought closer to the primary, the reflected light around $`5000\text{Å}`$ steadily increases. Nearly all ($`95\%`$) of the light reflected by the Cond model is due to Rayleigh scattering by the two most abundant species, H2 and He. At 0.3AU nearly $`10^7`$ times more light emerges from the planet around $`5000\text{Å}`$, however it remains very faint compared to the incident radiation with $`ϵ=2\times 10^7`$ between 4500 and 5200Å.

The younger $`\mathrm{T}_{\mathrm{int}}=1000\mathrm{K}`$ planet behaves in a similar manner as in the cases near a dM5 presented in section 3.1. The structure of the AMES-Cond planet, displayed in Figure Irradiated Planets, shows a significant suppression of the convective zone (below $`\tau _{std}=10`$ at 0.05AU) and the temperature has increased by a factor of 3 at $`\tau _{std}=10^4`$. The Dusty model also shows a large temperature inversion above $`\mathrm{log}(\mathrm{P}_{\mathrm{gas}})=3.0`$ dynes cm-2. As can be seen in Figure Irradiated Planets, the Dusty model at 0.15AU exhibits a complex mixture of cloud species throughout the atmosphere. Fe, Mg2SiO4, and MgSiO3 are the most prominent species except for the deeper layers where CaMgSi2O7 begins to dominate. The spectra, shown in Figure Irradiated Planets, display similar results as in the previous cases. However, at 0.05AU, the molecular bands have nearly disappeared in the Cond model and the once broad Na I and K I doublets (5890 and 7680Å) are extremely weak. Figure Irradiated Planets shows, in more detail, the steady reduction in equivalent width of both lines for the AMES-Cond models. The decrease in equivalent width is almost entirely due to the changes in thermal structure shown above. At 0.05AU, the Cond planet has reached $`\mathrm{T}_{\mathrm{eq}}=1752\mathrm{K}`$.

## 4 Discussion

### 4.1 Importance of Self-consistent Models

Many previous studies have used the structures of a non-irradiated planets with $`\mathrm{T}_{\mathrm{eff}}=\mathrm{T}_{\mathrm{eq}}`$ and simply computed a spectrum (or albedo) which included the incident flux and neglected the effects on the structure and chemical composition Marley et al. (1999). This procedure will, however, result in gross errors in the emergent flux from the optical to infrared. In Figure Irradiated Planets, we compare the thermal spectrum of the self-consistent irradiated (AMES-Cond, $`\mathrm{T}_{\mathrm{int}}=1000\mathrm{K}`$) planet at 0.065AU from a G2 star to that of a non-irradiated AMES-Cond model with an equal amount of thermal flux (i.e., $`\mathrm{T}_{\mathrm{eff}}=\mathrm{T}_{\mathrm{eq}}`$). Based on the irradiated spectrum without the reflected component, we estimate that $`\mathrm{T}_{\mathrm{eq}}=1560\mathrm{K}`$. The non-irradiated model significantly underestimates the flux in the water and methane bands and overestimates the flux in the regions outside the molecular bands. Also, the non-irradiated model overestimates the amount of reflected light blueward of 6000Å by $`35\%`$. Atomic features are also affected as can be seen for the Na I D line (Fig. Irradiated Planets) where the equivalent width is overestimated in the non-irradiated model even when the reflected light is included in the spectrum calculation. The differences between the structures (see Fig. Irradiated Planets) are just as striking. As might be expected, the non-irradiated model is still cooler (by 350$`\mathrm{K}`$) than the irradiated model in the upper atmosphere ($`\tau _{std}10^4`$). At deeper layers the two structures actually intersect and, at $`\tau _{std}10`$, the non-irradiated model is hotter than the irradiated model by roughly 350$`\mathrm{K}`$. An even more apparent difference is the location of the boundary between the radiative and convective regions. In the non-irradiated model, the convective zone reaches layers above $`\tau _{std}=1.0`$ while in the irradiated case, the convective zone has retreated to layers below $`\tau _{std}=10.0`$. These differences, between a non-irradiated structure with $`\mathrm{T}_{\mathrm{eff}}=\mathrm{T}_{\mathrm{eq}}`$ and a structure based on a self-consistent inclusion of the impinging radiation, will have significant consequences for any interior and evolution calculations of irradiated planets.

### 4.2 Limiting Effects of Dust

The question of how much energy is redistributed by ”weather patterns” is extremely important, and affects the upper boundary condition for irradiated models. To answer this question, 3D dynamical models of grain growth and diffusion would be needed. In this work, we explore the limiting effects of these patterns with Cond models, corresponding to clear skies, and Dusty models corresponding to cloudy skies. We find that irradiated Cond and Dusty models yield systematically different thermal spectra. The infrared spectra of Dusty and Cond planets around M dwarfs are affected by the impinging radiation only for extremely small orbits. However at separations around 0.05AU from a G2 star, the sudstellar point (the point on the planet closest to the star) will have a nearly featureless spectrum even in the Cond limit. The lack of spectral features is a consequence of the structure becoming nearly isothermal and is not a result of any significant decrease in abundances.

The emergent flux from both our irradiated and non-irradiated models indicates that Cond atmospheres are brighter at optical wavelengths than Dusty atmospheres. In general, Cond models are much brighter than the Dusty cases simply because the dust opacity blocks most of the thermal radiation (see Allard et al. (2001) for more details). When incident radiation is present, the dust grains reflect large amounts of light while continuing to block most of the intrinsic optical flux. Irradiated Cond atmospheres remain very transparent and allow large amounts of flux to emerge from the deep hotter layers. When this intrinsic thermal radiation is combined with the reflected flux, the irradiated Cond model appears brighter than the Dusty model.

The fraction of the stellar light reflected by our Cond model (at 0.05AU from a G2) around 4900Å is less than $`5\times 10^6`$. This is well below the results published by Cameron et al. (1999) and the upper limit published by Charbonneau et al. (1999). Also, the strong color dependence of Rayleigh scattering causes the reflected optical light (scattered by H2 and He) in the Cond models to be considerably different from the optical spectrum of the primary. However, the reflected light in the Dusty atmospheres is due almost entirely to Mie scattering which is a fairly grey process. The result is a near reflected copy of the stellar optical spectrum. A comparison between the stellar light and the reflected optical spectrum for a $`\mathrm{T}_{\mathrm{int}}=1000`$K planet located at 0.15AU from a G2 primary can be seen in Figure Irradiated Planets. For wavelengths less the $`4500\text{Å}`$, the reflected light in the Dusty model matches closely the stellar light. At redder wavelengths, the two spectra differ not only in the slope of the continuum, but also in the depth of H$`\alpha `$ and the Na I doublet at $`5890\text{Å}`$. The Cond model has a very different optical spectrum with the redder wavelengths being dominated by an extremely broad Na I doublet ($`5890\text{Å}`$). In general, the two spectra have little in common. Recent attempts have been made Cameron et al. (1999); Charbonneau et al. (1999) to observe a Doppler shifted copy of the stellar light reflected by the planet orbiting $`\tau _{std}`$ Boo. Figure Irradiated Planets suggests that planets in the Cond limit would be poor candidates for such techniques. Dusty planets, however, are clearly better choices and observations at shorter wavelengths may be more fruitful. In either case, observations similar to those of Cameron et al. (1999) and Charbonneau et al. (1999) could help determine the Cond or Dusty nature of EGPs.

### 4.3 Energy Redistribution

Clearly the fraction of the irradiated face seen by an observer depends on orbital phase and inclination. In addition, the planet’s atmospheric structure will likely vary across the planet’s surface. Guillot et al. (1996) claimed that one could assume that the radiation received by the planet’s day side is quickly redistributed over the entire planet surface. However, a recent calculation by Guillot (2000) indicates that this redistribution may take place over longer time-scales ($`10^5`$ seconds) than previously thought. This would imply that there exists, in those planets which are tidally locked, a large temperature difference between the day and night sides. Furthermore, if such a difference existed, the structure of the irradiated face would vary as a function of latitude and longitude and should approach the structure of the non-irradiated face for regions near the terminator111The terminator is the line dividing night and day on a planet.. Therefore, it is unlikely that an irradiated planet can be characterized by a single 1-D plane parallel model atmosphere. Though it is reassuring that our results agree qualitatively with those of previous studies Seager & Sasselov (1998); Seager et al. (2000); Goukenleuque et al. (2000), to accurately predict the spectrum and reflected light as a function of inclination and phase will require, at the very least, many 1-D models each accounting for the different amounts of energy deposited on the planet’s day side at various longitudes and latitudes. We are currently working on a sequence of non-isotropic irradiated models which include only the incident flux that a certain latitude222For this purpose, we assume that the symmetry axis is a line connecting the planet and primary. In this context, lines of constant latitude would refer to concentric rings about the substellar point. would receive. With such a sequence, we would be capable of modeling the reflected light for any phase and inclination.

Many of the recent studies Marley et al. (1999); Seager & Sasselov (2000); Seager et al. (2000); Goukenleuque et al. (2000); Sudarsky et al. (2000) have made the assumption that the planet receives an amount of energy equal to $`\pi \mathrm{R}_\mathrm{p}^2\mathrm{F}^{\mathrm{inc}}`$ which is quickly redistributed Guillot et al. (1996) over the entire planet surface. As a result, these previous works have essentially included the minimum amount of external flux ($`\frac{1}{4}\mathrm{F}^{\mathrm{inc}}`$) that could be received by a planet at a given separation. If the more recent calculation by Guillot (2000) is correct and global redistribution is inefficient, then the models in these previous studies would be valid only for separations larger than what the authors had originally intended and only for a specific ring around the substellar point. In our calculations, we have assumed no redistribution and therefore have included a much larger amount of external flux ($`\mathrm{F}^{\mathrm{inc}}`$). Despite the assumption of an isotropic external radiation field, we feel that our models adequately represent the substellar point where this assumption is more reasonable.

### 4.4 Transmission Spectra

In the recent study by Seager & Sasselov (2000) an attempt was made to model the transmitted flux of HD209458a (recently shown to be an eclipsing system) through the upper regions of the planet’s atmosphere. Seager & Sasselov (2000) assumed that the structure at the poles (and along the terminator) was the same as that given by their fully irradiated model. However, the planet receives only a small amount of incident flux at the terminator and only in the upper atmosphere. In the absence of strong energy redistribution, the terminator would have a structure more closely resembling that of the non-irradiated face.

For a transmission study, one must calculate the stellar flux passing through a thin region at the top of the planet’s atmosphere ($`0.01\mathrm{R}_p`$ thick) which includes latitudes above $`82^{}`$. While, in general, EGPs are adequately described by plane parallel geometry, as one approaches the limb, the plane parallel assumption becomes increasingly less accurate. Using plane parallel geometry to calculate the transmitted spectrum assumes that the stellar flux passes through a region of constant thickness and height (a slab). In reality, the flux is passing through a section of a sphere encompassing the limb, in which case, rays entering at different latitudes will pass through different amounts of the planet’s atmosphere. A simpler and more accurate way of calculating the transmitted flux is to solve the spherically symmetric radiative transfer equation (SSRTE) and take, as the transmitted flux, the average intensities for direction cosines which pass through this band. The main advantage over the plane parallel solution is that the correct geometry is automatically built into the SSRTE. When solving the SSRTE, the atmosphere is modeled as a discrete number of concentric shells surrounding the interior (or core). The transfer equation is then solved along characteristics which are divided into two categories: those which reach the core (core intersecting) and those which do not (tangential). Along a given characteristic, the direction cosine is now a function of depth whereas in the plane parallel case it is constant. Also, the diffusion approximation must hold at the inner boundary for the core intersecting characteristics but not for the tangential characteristics. The solution along the tangential characteristics which pass through the outer most shells is all that is needed to calculate the transmitted flux and will already account for the curvature of the limb.

Consider again the 1000$`\mathrm{K}`$ AMES-Cond planet located 0.05AU from a G2 primary and imagine that we observe the planet at (or near) inferior conjunction. Under such conditions one would observe a brighter planetary limb due to the transmitted stellar flux. In a spherically symmetric geometry, the effect emerges from the model in a very natural and physical manner. In Figure Irradiated Planets, we compare the transmitted flux for two different structures: our fully irradiated AMES-Cond structure ($`\mathrm{T}_{\mathrm{eq}}=1633\mathrm{K}`$) and the non-irradiated model ($`\mathrm{T}_{\mathrm{eff}}=1000\mathrm{K}`$). In both cases, the planet leaves its mark on the stellar flux as it passes through the planetary limb in the form of much stronger absorption features than seen in the unpolluted stellar spectrum. The transmitted spectrum based on the irradiated structure peaks at $`3\mathrm{\mu m}`$, but the non-irradiated planet (see Fig. Irradiated Planets) is still 100 times brighter in this region and over $`10^5`$ times fainter than the primary. In both cases, the planet’s limb is $`10^3`$ times brighter between 4500 and $`5000\text{Å}`$ than the planet’s night side. However, the disk of the G2 dominates the optical spectrum and is $`10^6`$ times brighter than the planet’s limb.

## 5 Conclusions

In this paper we have presented atmospheric models of planets in the presence of strong impinging radiation from primaries of different spectral type and at various orbital separations. We have also studied the effects of irradiation in two limiting cases: efficient settling (accounting for the depletion of elements by condensation) and complete cloud coverage.

Irradiation has only small effects on the atmospheres of planets with dM primaries except for extremely close orbits. For an object like Gl876B, which orbits a dM4 at 0.2AU, it is unlikely that any irradiation effects will be observed. However, for planets in the Cond limit orbiting G2 primaries, the effects are non-negligible. The upper layers of the atmospheres are significantly heated by the impinging radiation even at 0.5AU in both the Cond and Dusty limits. The inner layers of the Cond models also experience considerable heating and a suppression of the radiative-convective boundary. The innermost layers of the Dusty models are essentially unaffected by the irradiation even for close orbits. These results will have strong implications for interior and evolution calculations. The heating of the inner layers in the Cond limit can bring the temperatures at the bottom of the atmosphere close to those in the Dusty limit. This would suggest that, in certain situations, the interior models would no longer depend on the Dust-Cond uncertainty published by Chabrier et al. (2000). However, it is still important that detailed interior and evolution models be calculated which use irradiated atmospheres to set upper boundary condition.

The existence of “weather patterns” on EGPs in the form of strong zonal winds and clear, cloudy or partly cloudy skies, will greatly influence the effects caused by irradiation. Our models indicate that an EGP with significant cloud coverage will reflect a copy of the stellar optical light while EGPs with clear skies will have very different reflected spectra. We have also found that young EGPs with clear skies are actually brighter in the optical than cloudy EGPs despite the more efficient reflective properties of dust. In general, our models indicate that observations would need to be sensitive to variations in the planet+star spectrum on the order of $`10^5`$ to $`10^6`$ in order to disentangle the reflected or transmitted flux.

It is apparent from observations of Jupiter and the other Jovian planets that clouds are not homogeneously distributed over the surface and that complex weather patterns persists. There is also evidence (from Galileo) that, unlike our simplified dusty models, clouds form in thin decks and any condensation occurring at the uppermost layers would soon “rain” out onto the lower atmosphere possibly instigating a cascade of diffusive settling Encrenaz (1999). A self-consistent treatment of this behavior is currently being added to PHOENIX and results for both irradiated and non-irradiated atmospheres will be presented in future papers. Until these results are available, the models presented above are representatives of clear (AMES-Cond) and completely cloudy (AMES-Dusty) skies. It is possible to combine our Dusty and Cond models to investigate the intermediate “partly cloudy” cases.

Though it is unlikely that any of the models produced thus far are exact representations of a currently known EGP, one may hope that our predictions and those of others will aid observers in making the much needed detailed observations of EGPs.

We would like to thank Tristan Guillot, Adam Showman and Jean-Pierre Caillault for providing valuable comments and suggestions which greatly improved this paper. Travis Barman would especially like to thank everyone at CRAL, where a large part of this work was carried out. Their hospitality was greatly appreciated. This research was supported by the CNRS and NASA ATP grant NAG 5-3018, NAG 5-8425, LTSA grant NAG 5-3619 to the University of Georgia and LTSA grant NAG5-3435 and NASA EPSCor grant to Witchita State. Some of the calculations presented in this paper were performed on the IBM SP2 at the UGA UCNS and on the IBM SP2 of the San Diego Supercomputer Center (SDSC), with support from the National Science Foundation, and on the IBM SP2 of the NERSC with support from the DoE. We thank all these institutions for a generous allocation of computer time.


format: text

text:
# I Introduction

## I Introduction

Instantons play an important role in the QCD vacuum from the breaking of the axial $`U_A(1)`$ symmetry to chiral symmetry breaking and the low energy hadron spectrumSchafer and Shuryak (1998). Lattice studies support many of the theoretical predictionsPerez (2001). For example, the Witten-Veneziano formula relates the pure gauge topological susceptibility to the masses of the $`\eta ,\eta ^{}`$ and K mesons predicting $`\chi _{\mathrm{}}^{1/4}=180`$MeV. The topological susceptibility on pure gauge lattices has been measured by several groups using different methods obtaining consistent results $`\chi _{\mathrm{}}^{1/4}=205(5)`$MeVHasenfratz and Nieter (1998); de Forcrand et al. (1998); Alles et al. (1997).

Light fermions suppress the topological susceptibility. Chiral perturbation theory predictsLeutwyler and Smilga (1992)

$$\chi =\frac{m_q}{n_f}\mathrm{\Sigma }+O(m_q^2)=\frac{f_\pi ^2m_\pi ^2}{2n_f}+O(m_\pi ^4),$$
(1)
where $`m_q`$ is the quark mass, $`n_f`$ is the number of fermion flavors and $`\mathrm{\Sigma }`$ is the chiral condensate per fermionic flavor. In this normalization $`f_\pi =92`$MeV. Several recent lattice studies measured $`\chi `$ on dynamical configurations. Calculations with Wilson and clover fermions cover a fair range of lattice spacing $`a0.080.20`$fm and pion masses $`m_\pi r_01.32.5`$ ($`r_0`$ is the Sommer parameter Sommer (1994)) Hart and Teper (2000a, b, c); Bali et al. (2001); Khan et al. (2000). The results appear controversial. UKQCD uses clover fermionsHart and Teper (2000a, b, c). Their data on lattices with $`a0.1`$fm is basically consistent with eqn.(1). SESAM/T$`\chi L`$ uses unimproved Wilson fermions. Their topological susceptibility at similar lattice spacing does not decrease with the pion mass though the statistical errors are too large to claim inconsistency with theoretical expectationsBali et al. (2001). CP-PACS published data using clover fermions and improved gauge action at lattice spacing $`a0.17`$fmKhan et al. (2000). Their conclusion is the same as SESAM/T$`\chi L`$. The situation with staggered fermions is not much better. Only the Pisa group measured the topological susceptibility with two and four flavors of staggered fermions at lattice spacing $`a0.090.17`$fmAlles et al. (2000a, b). They do not see the reduction of the topological susceptibility with decreasing quark mass either.

Can we understand what is going on with the dynamical simulations? Eqn.(1) is valid only in case of $`n_f`$ light *chiral* fermions creating $`n_f^21`$ light pions. Both Wilson and staggered fermions violate this assumption. Wilson fermions break chiral symmetry explicitly. The addition of the clover term reduces the symmetry breaking and improves chiral behavior. Staggered fermions have only a residual $`U(1)`$ chiral symmetry and only one true Goldstone boson, the other pion-like states can be heavy. Fat link fermions considerably improve flavor symmetry and consequently chiral behavior. Both fermionic formulations should show the expected chiral behavior in the continuum limit, clover and fat link fermions sooner than the unimproved ones, but it is not clear when scaling in the topological susceptibility actually sets in. In a recent paper Durr (2001) it was suggested that the reason the topological susceptibility from available lattice data does not follow the expected theoretical behavior is a combination of three effects: too large lattice spacing, too small volumes and too small Leutwyler-Smilga Leutwyler and Smilga (1992) parameter. While all three conditions could indeed be important, we feel that the non-chiral behavior of the fermionic actions is the main cause of the problem. Chiral symmetry violation of the fermionic action is a scaling violation effect and is covered by the first condition of Durr (2001). However, it is not a universal quantity, it can strongly depend on the fermionic action.

Why are the topological properties of the vacuum important? Phenomenological instanton models predict that the low energy hadron spectrum is governed by the near-zero eigenvalue modes of the Dirac operator, which, in turn, are related to instantons. If the fermion-instanton interaction is different from the continuum one either because of chiral symmetry breaking or other lattice artifacts, the Dirac spectrum and consequently the low energy hadron spectrum can also be different. Recently, using chiral symmetric overlap fermions on the lattice, we showed that with light quarks the first few modes of the Dirac operator saturate the low lying hadron propagators on quenched $`a0.12`$ lattices, just as the phenomenological models predictDeGrand and Hasenfratz (2000, 2001). In a subsequent publication contradictory conclusion was reached using Wilson fermionsHorvath et al. (2001). For us that implies that Wilson fermions, at least at large to moderate lattice spacings, interact differently with instantons than chirally symmetric fermions.

In this paper we investigate the topological susceptibility on both $`n_f=2`$ and $`n_f=4`$ staggered dynamical fermion configurations. The $`n_f=2`$ configurations are $`16^332`$ thin link staggered fermion lattices created by the Columbia and MILC collaborations and downloaded from the NERSC archiveNER (http://qcd.nersc.gov/). The $`n_f=4`$ lattices are smaller, $`8^324`$ configurations used in the study of fat link fermions in Ref. Knechtli and Hasenfratz (2000). Two of the $`n_f=4`$ sets were generated with thin link staggered fermions and one with N=3 level APE blocked fat link fermions. The latter action has about an order of magnitude better flavor symmetry than thin link actions at similar parameter values. The thin link staggered fermion results with $`a0.1`$fm are more or less consistent with eqn.(1). The thin link $`a0.17`$fm data shows clear deviation, the topological susceptibility is consistent with the quenched value, independent of the pion mass. On the other hand, the fat link data at the same coarse lattice spacing is in perfect agreement with eqn.(1), suggesting that improved flavor/chiral symmetry indeed has a strong effect on the topology.

To determine the topological charge of the configurations we used a topological charge density operator constructed from hypercubic blocked (HYP) fat linksHasenfratz and Knechtli (2001a). Hypercubic blocking was introduced in a recent paper as an alternative to repeated APE smearingAlbanese et al. (1987). Hypercubic blocking mixes links only from the hypercubes that attach to a given link, thus the fat link is very compact yet the configuration is as smooth as after three levels of APE blocking. To avoid the distortion effect of extended operators, we consider only one to three levels of hypercubic blocked operators (HYP1, HYP2 and HYP3). We have calculated the additive and multiplicative renormalization factors of the topological susceptibility for these operators. After two to three levels of HYP blocking the renormalization factors turn out to be consistent with their tree level values and, after correction, all three HYP topological susceptibility measurements are consistent.

In Sect. 2. we describe the hypercubic topological operator and illustrate the measurement of the renormalization factors on pure gauge Wilson plaquette configurations. Sect. 3 contains our results for two- and four-flavor staggered fermions. Sect. 4 is a short discussion on flavor symmetry and the summary of our results.

## II Hypercubic Topological Charge Operator

Most large scale simulations use pure gauge observables to measure the topological susceptibility with some discretized version of the continuum $`F\stackrel{~}{F}`$ as the lattice charge density operator $`q_L(x)`$. The relation between the continuum and lattice topological susceptibilities contains both a multiplicative and an additive renormalization factor. The lattice charge density operator is related to the continuum one through a multiplicative renormalization factor

$$q_L(x)=Za^2q(x).$$
(2)
In addition the correlator of two topological density operators have an additive correction term as well due to the mixing of $`q_L(x)`$ with other lattice operators

$$q_L(x)q_L(0)=Z^2a^4q(x)q(0)+m(x),$$
(3)
thus the topological susceptibility on the lattice is

$$\chi _L=q(x)q(0)d^4x=Z^2a^4\chi +M,$$
(4)
$$M=m(x)d^4x.$$
(5)
The renormalization factors Z and M depend on the lattice parameters $`\beta `$ and $`m_q`$. The renormalization constants in principle can be determined non-perturbatively using the heating method proposed by the Pisa groupChristou et al. (1996).

Due to lattice artifacts, mainly dislocations, on typical lattices $`Z0.25`$ and $`M/\chi _L1`$ for thin link topological density operators which makes the direct determination of $`\chi `$ almost impossible. Local smearing or cooling removes most of these lattice artifacts moving Z towards 1 and M to 0. Repeated smoothing gives $`Z1`$ and $`M/\chi 1`$ leading to $`\chi _La^4\chi `$. Since the non-perturbative calculation of Z and M introduces statistical and systematical errors, a topological density operator where the renormalization constants are small or can be neglected is desirable. However, repeated smearing and cooling methods have their drawbacks as well. While removing vacuum fluctuations and lattice artifacts, both methods remove topological objects, mainly small instantons and close-by pairs. In addition the size of the remaining objects change as well, as can be demonstrated by monitoring individual instantons during the smearing processHasenfratz and Nieter (1998).

In this paper we construct the topological charge density operator $`q_L(x)`$ using hypercubic blocked links and the improved thin link charge operator of Refs. DeGrand et al. (1998) and Hasenfratz and Nieter (1998). Hypercubic blocking was introduced in Ref. Hasenfratz and Knechtli (2001a) as a localized alternative to repeated APE blocking. HYP fat links mix original links from the hypercubes that are attached to the fat link only, yet they create configurations that are as smooth as the ones obtained after three levels of APE blocking.

Topological charge measurement are difficult because of the presence of dislocations: short distance vacuum fluctuations that can be mistaken for small instantons. Smearing the links in the topological density operator reduces this problem in part by sharpening the transition between the topologically non-trivial $`Q=d^4xq(x)=1`$ sector and the trivial $`Q=0`$ sector. Figure 1 shows the topological charge of a smooth instanton measured with the thin link, HYP1, HYP2 and HYP3 fat link operators as the function of the instanton radius .

The charge measured with the HYP operators rises sharply at instanton radius $`\rho /a1`$.

Even though the HYP operators remove most dislocations, they still can have non-trivial renormalization factors. We have measured the renormalization factors Z and M following the heating method of Ref. Christou et al. (1996). To measure M one has to heat a trivial $`Q=0`$ configuration. The configuration thermalizes fast at short distances, and since the origin of M is local contact terms, short distance thermalization is sufficient to measure M. Before heating creates non-trivial topological objects the thermalization must be terminated and restarted with a different random seed. It is important to make sure the measurement of M is done on $`Q=0`$ configurations only. Fortunately with the HYP3 and frequently with the HYP2 operators it is easy to separate the trivial configurations from the occasional $`Q0`$ ones. Figure 2/a shows the topological charge on 500 heated $`Q=0`$ configurations measured with HYP2 operator. Ordered $`8^4`$ configurations were heated with 10-50 heat bath steps using $`\beta =6.0`$ pure gauge Wilson action. Most configurations have very small topological charge, a few has $`|Q|1`$ and only 2-3 configurations have topological charge whose interpretation is not clear. I chose, based on Figure 1 but somewhat arbitrarily, a cut $`Q_{max}=0.3`$ to separate the $`Q=0`$ and $`Q0`$ configurations. The final results are fairly insensitive to the precise choice of this cut. The topological charge measured with the HYP1 operator fluctuates more as is illustrated in Figure 2/b. Again, a configuration is accepted as $`Q=0`$ if its topological charge is $`|Q_L|<0.3`$ as measured with the HYP2 operator. The multiplicative renormalization constant Z is calculated similarly on $`Q=1`$ configurations that contain a single smooth instanton of size $`\rho /a=3.0,2.5`$ or 2.0. Z can be different on different instanton size backgrounds, though in none of the measurements we performed was the difference more than a few per cent. We chose Z on the background where the smooth instanton’s size was closest to the expected average instanton size of the Monte Carlo configurations, $`\rho /a0.3`$fm$`/a`$.

Results for $`\beta =6.0`$ pure gauge Wilson action are summarized in Table 1 where Z, M, $`\chi _L`$, $`\chi a^4=(\chi _LM)/Z^2`$ and $`\chi r_0^4`$ with Sommer parameter $`r_0=5.37(1)`$ are given both for the HYP1 and HYP2 operator, and, for comparison, one level APE smeared operator with parameter $`\alpha =0.75`$ as well. The topological susceptibility was measured on 220 $`16^332`$ $`\beta =6.0`$ configurations from the NERSC archiveKilcup et al. (1997). The results for $`a^4\chi `$ are consistent for all three operators, in physical units $`\chi ^{1/4}=196(5)`$MeV. The renormalization factors, on the other hand, are quite different. The background term M is about 17% of $`\chi _L`$ for the APE1 operator, 2% for the HYP1 operator and less than 0.1% for the HYP2 operator. The renormalization factors for the HYP2 operator are negligible, for all practical purposes we can use $`Z=1`$ and $`M=0`$.

Results are similar at $`\beta =5.7`$ though there one needs three levels of HYP blocking to reduce M to 0 and Z to 1 as is illustrated in Table 2. Here the topological susceptibility was measured on 350 $`8^324`$ lattices and the renormalization constants are obtained from 500 heated $`8^4`$ lattices. The predictions for the continuum topological susceptibility are consistent from all three operators, $`\chi ^{1/4}=188(9)`$MeV. Within statistical errors it is also consistent, though about 15% lower, then the value obtained on the $`\beta =6.0`$ data set. Whether the difference is only statistical fluctuation or the $`a0.17`$fm lattice spacing at $`\beta =5.7`$ is too coarse to support all the topological objects on the lattice cannot be decided with the present statistics.

## III The Topological Susceptibility on dynamical configurations

### III.1 Two flavors of staggered fermions

We have analyzed several two-flavor configuration sets using HYP1-HYP3 operators. The renormalization factors have to be calculated independently at every parameter value. The results for both Z and M were very similar to the quenched values at similar lattice spacing. We found that on configurations with lattice spacing $`a0.1`$fm the renormalization constants of the HYP2 operator could be neglected, while configurations with $`a0.17`$fm the HYP2 operator had a few percent correction from the renormalization constants, the HYP3 operator had none.

Table 3 collects our results. All four data sets are from the NERSC archive, the first three were generated by the Columbia group, the last one by MILCBrown et al. (1991); Bernard et al. (1999). All lattices are $`16^332`$ and use standard thin link staggered fermions with plaquette Wilson gauge action. $`am_\pi `$ was measured in the original studies, $`r_0/a`$ for the MILC lattice is from Ref. Tamhankar and Gottlieb (2000). For the Columbia lattices we have measured $`r_0/a`$ using HYP blocked Wilson loops. The HYP potential has greatly reduced statistical errors making it possible to obtain reliable values even from 33-83 configurationsHasenfratz and Knechtli (2001a). The lattice spacing in the fourth column was obtained using $`r_0=0.5`$fm and is listed for future reference. Since the different HYP topological charge operators give consistent results, only the continuum value $`a^4\chi `$ is listed. The topological susceptibility for the first data set has been measured in Ref. Hasenfratz (2000) using 10-40 APE smeared operators. The values for $`a^4\chi `$ with APE smeared operators are consistent with the present result.

In Figure 3/a $`\chi r_0^4`$ is plotted as the function of $`(m_\pi r_0)^2`$. The filled octagons correspond to the Columbia data sets with $`a0.1`$fm. The filled diamond is the MILC lattice data with $`a0.17`$fm. The filled square at the lowest $`m_\pi r_0`$ value is also from MILC. It is on a $`24^364`$, $`\beta =5.6`$, $`am_q=0.01`$ configuration set and was measured by C. DeTar using 10-20 level APE blocked topological operatorsDeTar (1998). The lattice spacing on these configurations is $`a0.1`$fm, $`r_0/a=4.99`$. The dashed line in Figure 3 is the leading order theoretical prediction from eqn.(1) using $`f_\pi =92`$MeV. The three lowest $`m_\pi r_0`$ data points on the finer lattices with $`a0.1`$fm are consistent, though a bit higher then the leading order theoretical curve. However, the data point with $`a0.17`$fm is very different. The topological susceptibility on those lattices is consistent with the quenched result even though the pion mass is fairly small. Since on quenched configurations the topological susceptibility can be measured successfully even on coarse lattices, this discrepancy is not likely to be the consequence of the gauge system. Rather, it appears that the instantons do not feel the effect of light staggered fermions on coarse lattices. Flavor symmetry violation of staggered fermions can explain these findings. In Figure 3/b published Wilson/clover data is added to the $`a0.1`$fm staggered data. Octagons correspond to the clover fermion simulations of UKQCD, diamonds to the Wilson fermion simulations of SESAM/T$`\chi `$L. All Wilson/clover data has lattice spacing $`a0.1`$fm. The data from UKQCD is consistent with the staggered results and with the leading order theoretical curve as well. However the topological susceptibility with unimproved Wilson fermions does not follow the expected behavior at small pion mass. It remains large, consistent with the quenched result, the same behavior we saw with staggered fermions at lattice spacing $`a0.17`$fm. It appears that Wilson fermions have a very different effect on instantons even on finer, $`a0.1`$fm lattices than clover or even unimproved staggered fermions. This is likely the consequence of chiral symmetry breaking of the Wilson action.

### III.2 Four flavors of staggered fermions

Unfortunately we could not find any large $`n_f=4`$ staggered data sets to use. In Ref. Knechtli and Hasenfratz (2000), where we proposed a dynamical fat link fermion update, we generated three $`8^324`$ $`n_f=4`$ data sets to study flavor symmetry breaking of thin and fat link actions. The first two sets were generated with standard thin link staggered fermions, the third one with N=3 APE blocked fat link fermions. All three sets have lattice spacing $`a0.17`$fm as shown in Table 4. The pion mass values in the table are also from Ref.Knechtli and Hasenfratz (2000) but we have re-analyzed the potential data using HYP Wilson loops and the $`r_0/a`$ values listed in Table 4 are more reliable and slightly different from the ones used in Ref. Knechtli and Hasenfratz (2000)

The difference in the topological susceptibilities between the thin and fat link actions is striking. There is no Wilson/clover fermion data available at $`n_f=4`$ to compare the staggered result. In Figure 4 we plot $`\chi r_0^4`$ as the function of $`(m_\pi r_0)^2`$ for the three data sets from Table 4. The dashed line is the leading order theoretical prediction from eqn.(1) using the experimental value $`f_\pi =92`$MeV. The fat link action data shows the expected behavior while the thin link action topological susceptibility values are consistent with the quenched result, independent of the dynamical quark mass. For the thin link action this is the same behavior we saw with the $`n_f=2`$ data at similar lattice spacing.

## IV Discussion and Summary

The interaction between light quarks and instantons changes the QCD vacuum substantially. One of the easiest way to get information about this effect is through the topological susceptibility that is expected to scale with the square of the pion mass in the small quark mass limit.

Our results indicate that thin link staggered fermions at lattice spacing $`a0.17`$fm do not have the expected continuum effect on instantons, a lattice spacing $`a0.1`$fm is needed to recover the proper chiral behavior. In contrast to that, fat link staggered fermions show the expected behavior even on coarser $`a0.17`$fm lattices. Can this be understood as the effect of flavor symmetry violation? In Refs. Knechtli and Hasenfratz (2000); Hasenfratz and Knechtli (2001a) we studied flavor symmetry violation on quenched lattices both with the thin link, N=3 APE blocked and hypercubic blocked fat link fermions. The two fat link formulations have about the same level of flavor symmetry violations as measured by the relative mass splitting $`\mathrm{\Delta }_\pi =(m_\pi m_G)/m_G`$ between the Goldstone pion $`m_G`$ and the other pion like objects, $`m_\pi `$. For thin link fermions at $`(m_\pi r_0)^22.0`$ we found $`\mathrm{\Delta }_\pi 0.7`$ at $`a0.17`$fm, $`\mathrm{\Delta }_\pi 0.2`$ at $`a0.1`$fm for the lightest non-Goldstone pions. The corresponding values for fat link fermions are $`\mathrm{\Delta }_\pi 0.09`$ and $`\mathrm{\Delta }_\pi 0.01`$, respectively. QCD with two(four) light flavors should have 3(15) light pions. Apparently when there is only one light pion and the other pion-like objects are 70% or more heavier than the Goldstone pion, the vacuum does not look like a two(four)-flavor QCD vacuum. The flavor symmetry breaking has to be reduced below 20% to get acceptable results. Fat link fermions can do that easily even on coarse lattices.

Available lattice data for Wilson/cover fermions can be understood similarly. The topological susceptibility with clover fermions are consistent with the theoretical predictions at lattice spacing $`a0.1`$fm but the Wilson fermion data at the same lattice spacing indicates that unimproved Wilson fermions, that have much larger chiral symmetry violations, do not have the expected continuum like interaction even at $`a0.1`$fm. Results from CP-PACS at lattice spacing $`a0.17`$fm indicate that at that lattice spacing not even clover fermions can reproduce the continuum topological behavior. Fat link clover fermions improve chiral symmetry and could provide a better alternative to thin link clover fermions. It would be interesting to find a parameter similar to $`\mathrm{\Delta }_\pi `$ that correlates with the level of chiral symmetry breaking for the Wilson like fermions and compare it with the topological susceptibility.

One should be concerned about dynamical simulations where the topological susceptibility is not reproduced correctly, since that indicates that the vacuum at those simulations is more like the quenched vacuum rather than the expected dynamical one. Improving chiral symmetry can have a profound effect. This point is further underscored in a forthcoming publication about the finite temperature phase diagram obtained with fat link fermionsHasenfratz and Knechtli (2001b).

###### Acknowledgements.

This work was largely inspired by the collaboration with F. Knechtli on fat link fermions and hypercubic blocking. G. Bali convinced me that in order to compare the topological susceptibilities from different simulations one needs to use a reliable scale, like $`r_0`$. That prompted the potential measurements with hypercubic blocked links. S. Gottlieb helped me out with some of the MILC staggered fermion data. The numerical calculations of this work were performed on the Colorado-HEP $`\alpha `$ cluster.


text:
# PREDICTION OF SUPERSOFT X-RAY PHASE, HELIUM ENRICHMENT, AND TURN-OFF TIME IN RECURRENT NOVA CI AQUILAE 2000 OUTBURST

## 1 INTRODUCTION

The second recorded outburst of CI Aquilae was discovered in April 2000 by Takamizawa (2000) since the first recorded outburst in 1917 (Reinmuth, 1925; Duerbeck, 1987). Thus, CI Aql has been recognized as a member of the recurrent novae (Yamaoka, Ayani, & Shirakami, 2000). After optical brightness reached its maximum $`V9`$ in early May 2000, it rapidly decreased to $`V13`$ in about 50 days. A plateau phase follows; the brightness levelled off at $`V13.5`$. CI Aql is still as bright as $`V14`$ until now about 300 days after maximum, being the longest record of on-time among the recurrent novae (e.g., VSNET archives, http://www.kusastro.kyoto-u.ac.jp/vsnet/).

Recurrent novae constitute a small class of cataclysmic variables, the outbursts of which repeat in a decade to a century. Ejecta of the recurrent novae are not metal-enriched but similar to the solar value ($`Z0.02`$), indicating that the WDs in the recurrent novae are not eroded but have grown in mass (e.g., Starrfield, Sparks, & Truran, 1985). Masses of white dwarfs (WDs) in the recurrent novae are estimated to be so massive as the Chandrasekhar mass limit from the light curve fitting (e.g., Hachisu & Kato, 1999; Hachisu et al., 2000a; Hachisu & Kato, 2000a, b, 2001). These features strongly suggest that the recurrent novae are an immediate progenitor of Type Ia supernovae (SNe Ia) if the WD consists of carbon and oxygen (Nomoto, Thielemann, & Yokoi, 1984). Evolutionary paths to SNe Ia via recurrent novae have been pointed out by Hachisu, Kato, & Nomoto (1999a).

Recurrent novae are divided into three subclasses depending on their orbital periods: T Pyx subclass with a dwarf companion, U Sco subclass with a slightly evolved main-sequence (MS) or subgiant companion, and RS Oph subclass with a red giant companion (e.g., Warner, 1995). CI Aql has an orbital period of 0.618 days, belonging to the U Sco subclass. It is therefore likely that the light curve of CI Aql develops like U Sco.

The 1999 outburst of U Sco was densely observed from the rising phase to the final decay phase through the mid-plateau phase. The detection of supersoft X-rays during the plateau phase (Kahabka et al., 1999) strongly supports the theoretical description of nova outbursts — massive wind phase, steady hydrogen burning phase (no winds), and final cooling phase, — which corresponds to the rapid decline, plateau, and final decline phases in the light curve (e.g., Hachisu et al., 2000a), respectively. Because the development of CI Aql outbursts is much slower than that of U Sco, there is still a chance to observe supersoft X-rays. Supersoft X-rays are direct evidence of hydrogen shell burning on a WD so that they provide us valuable information on recurrent novae and SN Ia progenitors.

In this Letter, we determine the mass and the growth rate of the WD as accurately as possible by theoretically analyzing light curves. Then, (1) we predict the turn-off time of nuclear burning and the duration of supersoft X-ray source (SSS) phase of the present outburst and (2) examine the possibility whether or not CI Aql will explode as an SN Ia. In §2, we describe our model of CI Aql in quiescence and estimate the inclination angle, the size and thickness of the accretion disk, and so on. In §3, the 1917 outburst is numerically analyzed and we determine the mass of the WD. In §4, the 2000 outburst is also numerically analyzed and the duration of wind is estimated. In §5, we predict the duration of the SSS phase. Discussion follows in §6.

## 2 BASIC MODEL IN QUIESCENCE

CI Aql is an eclipsing binary. We have considered a binary model as illustrated in Figure 1 which consists of (1) the WD photosphere, (2) the MS photosphere which fills its Roche lobe, and (3) the accretion disk (ACDK). We assume that the size and thickness of the ACDK are given by two parameters, $`\alpha `$ and $`\beta `$, as

$$R_{\mathrm{disk}}=\alpha R_1^{},$$
(1)
and

$$h=\beta R_{\mathrm{disk}}\left(\frac{\varpi }{R_{\mathrm{disk}}}\right)^\nu ,$$
(2)
where $`R_{\mathrm{disk}}`$ is the outer edge of the ACDK, $`R_1^{}`$ the effective radius of the inner critical Roche lobe for the WD component, and $`h`$ the height of the surface from the equatorial plane, $`\varpi `$ the distance on the equatorial plane from the center of the WD. Here, we adopt $`\varpi `$-square law ($`\nu =2`$) to mimic the effect of flaring-up at the rim of the ACDK (e.g., Schandl, Meyer-Hofmeister, & Meyer, 1997).

The accretion luminosity of the WD (e.g., Starrfield, Sparks, & Shaviv, 1988), the viscous luminosity of the ACDK, and the irradiations by the WD (e.g., Schandl et al., 1997) are also included. Here, we assume that the surfaces of the MS and of the ACDK emit photons as a black-body at a local temperature of each surface patch area (see Fig. 1) which is heated by radiation from the WD. The efficiencies of irradiation of the MS and of the ACDK are assumed to be $`\eta _{\mathrm{MS}}=0.5`$ and $`\eta _{\mathrm{DK}}=0.5`$, respectively. The numerical method is essentially the same as that in Hachisu et al. (2000b) and fully described in Hachisu & Kato (2001). A circular orbit is assumed. Light curves of CI Aql in quiescence has been obtained by Mennickent & Honeycutt (1995). They determined the orbital period and ephemeris to be JD 2,448,412.167(25) + 0.618355$`\times E`$ at eclipse minima.

Figure 2 shows the averaged (smoothed) observational points (open circles) together with our calculated $`V`$ light curves for the suggested WD mass of $`M_{\mathrm{WD}}=1.2M_{}`$ ($`R_{\mathrm{WD}}=0.0072R_{}`$) and mass accretion rate of $`\dot{M}_{\mathrm{acc}}1\times 10^7M_{}`$ yr-1 (see the following sections). To fit our theoretical light curves with the observational points, we have calculated $`V`$ light curves by changing the parameters of $`\alpha =0.5`$—1.0 by 0.1 step, $`\beta =0.05`$—0.50 by 0.05 step, $`T_{\mathrm{ph},\mathrm{MS}}=5000`$—8000 K by 100 K step, $`T_{\mathrm{ph},\mathrm{disk}}=3000`$—6000 K by 500 K step, and $`i=70`$$`90\mathrm{°}`$ by $`1\mathrm{°}`$ step and seek for the best-fit model by minimizing RMS of the residuals. The light curves are also calculated for four companion masses of $`M_{\mathrm{MS}}=0.8`$, 1.1, 1.5 and $`2.0M_{}`$. Since we obtain similar light curves for all of these four cases, we adopt $`M_{\mathrm{MS}}=1.5M_{}`$ because the mass transfer is as high as $`1\times 10^7M_{\mathrm{}}`$ yr-1 and this indicates a thermally unstable mass transfer (van den Heuvel et al., 1992), which requires $`M_{\mathrm{MS}}/M_{\mathrm{WD}}1.1`$ (Webbink, 1985). A 1.5 $`M_{\mathrm{}}`$ MS expands to fill the Roche lobe of $`1.7R_{\mathrm{}}`$ after a large part of hydrogen at the center is consumed (e.g., Hachisu et al., 1999b).

In the best-fit model, a rather blue color index of $`(BV)_c0.14`$ is obtained as shown in Figure 3. This suggests a large color excess of $`E(BV)=(BV)_o(BV)_c0.86`$ with the observed color of $`(BV)_o1.0`$ (Mennickent & Honeycutt, 1995). This color excess is very consistent with $`E(BV)=0.85`$ estimated by Kiss et al. (2001).

## 3 THE 1917 OUTBURST

First we reproduce the light curve of the 1917 outburst. Its Tycho $`B`$ magnitude light curve has recently been reported by Williams (2000), which is shown in Figure 4 together with our calculated $`B`$ light curves. Our numerical method for calculating light curves is fully described in Hachisu & Kato (2001).

It has been established that the WD photosphere expands to a giant size at the optical maximum in nova outbursts. Since the separation of CI Aql is $`a4R_{\mathrm{}}`$ (Greiner & Alcala, 1996), the WD envelope engulfs the MS companion so that the early phase light curve is calculated only from the WD photosphere (e.g., Kato & Hachisu, 1994; Kato, 1999). The accretion disk and the companion do not contribute to the $`B`$ light until it decreases to $`B13`$ mag. We determine the WD mass to be $`M_{\mathrm{WD}}=1.2\pm 0.05M_{\mathrm{}}`$ from the early phase light curve as seen in Figure 4.

Assuming $`M_{\mathrm{WD}}=1.2M_{\mathrm{}}`$, we estimate the envelope mass at the optical maximum to be $`\mathrm{\Delta }M=5.8\times 10^6M_{}`$ for the hydrogen content $`X=0.7`$. If the duration of quiescence before the 1917 outburst is almost the same as the quiescent phase between the 1917 and 2000 outburst, this envelope mass indicates a mass accretion rate of $`\dot{M}_{\mathrm{acc}}0.7\times 10^7M_{}`$ yr-1 if no WD matter has been dredged up. In our wind model, about 81% ($`4.7\times 10^6M_{}`$) of the envelope mass has been blown in the wind while the residual 19% ($`1.1\times 10^6M_{}`$) has been left and added to the helium layer of the WD. The net mass-increasing rate of the WD is $`\dot{M}_{\mathrm{He}}=0.13\times 10^7M_{}`$ yr-1. Thus, CI Aql may become an SN Ia if the donor is massive enough.

## 4 THE 2000 OUTBURST

CI Aql bursted in April 2000. Its optical maximum ($`V9`$ mag) was reached on May 5, 2000 (HJD 2451669.5, Kiss et al., 2001). As shown in Figure 5, the visual brightness quickly decreased to 13.5 mag in about 60 days. Then, it stays at $`V14`$ mag, i.e., a plateau phase. This mid-plateau phase is very similar to that of U Sco and can also be explained by the irradiation of the accretion disk (Hachisu et al., 2000a). Therefore, we have reproduced the light curve in the plateau phase by assuming similar model parameters to those of U Sco. Theoretical light curves in Figure 5 are calculated for a pair of $`1.2M_{}`$ WD $`+`$ $`1.5M_{}`$ MS. In the plateau phase of the light curve, i.e., when the WD photosphere becomes much smaller than the binary size, the light curve is determined mainly by the irradiations of the ACDK and of the MS. The unheated surface temperatures are fixed to be $`T_{\mathrm{ph},\mathrm{MS}}=6300`$ K for the MS and to be $`T_{\mathrm{ph},\mathrm{disk}}=5300`$ K at the disk rim, which are the same as the basic model in quiescence.

The luminosity of the accretion disk depends on both the thickness $`\beta `$ and the size $`\alpha `$. We have examined a total of 160 cases for the set of $`(\alpha ,\beta )`$, which is the product of 16 cases of $`\alpha =`$ 0.5—2.0 by 0.1 step and 10 cases of $`\beta =`$ 0.05—0.50 by 0.05 step. Here, we have adopted a set of $`\alpha =1.4`$ and $`\beta =0.2`$ during the wind phase and a set of $`\alpha =0.7`$ and $`\beta =0.2`$ after the wind stopped.

The decline of the early phase ($`t0`$—70 days) hardly depends on the hydrogen content $`X`$ of the WD envelope but the ends of the wind phase and of the hydrogen shell-burning phase depend sensitively on the hydrogen content $`X`$ as shown in Figure 5. Matsumoto & Kato (2000) reported a sharp $`1`$ mag drop of $`R`$ magnitude in late November 2000 and Kiyota (VSNET archives) also observed a $`1.5`$ mag drop of $`I`$ magnitude as shown in Figure 5. If we attribute this drop to the end of wind phase, the hydrogen content is somewhere between $`X=0.35`$ and $`X=0.5`$, more closer to $`X=0.5`$. We here adopt $`X=0.5`$. In the case of $`X=0.5`$, the continuous wind stopped in early December 2000 about 220 days after maximum and the steady hydrogen shell-burning will end in early August 2001 about 450 days after maximum as shown in Figure 5.

## 5 SUPERSOFT X-RAY SOURCE PHASE

We here strongly encourage X-ray observations until the mid August 2001, that is, the end of steady hydrogen shell-burning. No supersoft X-rays are observed in the early decline phase because the photospheric radius of the WD is rather large and the photospheric temperature is relatively low. After the recurrent nova enters the plateau phase, the photospheric temperature of the WD increases to high enough to emit supersoft X-rays. During the massive wind phase, however, we do not expect supersoft X-rays because they are self-absorbed by the wind itself. U Sco was observed as a luminous supersoft X-ray source in the mid-plateau phase of the 1999 outburst just after the massive wind stopped (Kahabka et al., 1999; Hachisu et al., 2000a). Our analysis suggests that the massive wind stopped in early December 2000 and that supersoft X-rays are now very likely detectable until August 2001.

## 6 DISCUSSION

The calculated distance modulus of the 1917 outburst is about $`(mM)_B=14.5`$. Therefore, the distance to CI Aql is estimated to be 1.6 kpc for the reddening of $`E(BV)=0.85`$ and the absorption of $`A_B=3.5`$ and $`A_V=2.7`$. The apparent distance modulus of $`(mM)_V=13.7`$ from the 2000 outburst as shown in Figure 5 is also consistent with the distance of 1.6 kpc and $`A_V=2.7`$.

On the other hand, we obtain an apparent distance modulus of $`(mM)_V=12.86`$ in the quiescent phase as shown in Figure 2, which corresponds to the distance of 1.1 kpc. This is not consistent with our estimation from the outburst light curves. To examine this difference, we check other observations and found that three others indicate $`V16.2`$ in quiescence, which is 0.7 mag darker than $`V15.5`$ in Figure 1 of Mennickent & Honeycutt (see Table 2 of Mennickent & Honeycutt, 1995; Szkody & Howell, 1992). The $`V`$ magnitudes in their Figure 1 were obtained with a fully automated telescope and the errors in the zero point (formally 0.2 mag as in their text) could be larger than expected (private communication with Mennickent). If we adopt $`V=16.2`$ in quiescence, then the distance modulus is $`(mM)_V=13.56`$ and the distance to CI Aql becomes $`d=1.5`$ kpc, being consistent with the distance estimate from the outbursts. Thus, we determine the distance to CI Aql to be $`d1.6`$ kpc with $`A_V=2.7`$.

The lacks of flickering and of Balmer emission lines in quiescence lead Mennickent & Honeycutt (1995) to conclude that CI Aql is not interacting at least at the present time. However, the inner part of the accretion disk cannot be seen from the Earth as easily understood from Figure 1. This is the reason for the lacks of emission lines and of flickering even when the system includes a mass-accreting WD. This is essentially the same situation as in U Sco (Schaefer, 1990; Johnston & Kulkarni, 1992; Schaefer & Ringwald, 1995; Hachisu et al., 2000b).

We are very grateful to the VSNET members who observed CI Aql and sent their valuable data to the VSNET. We also thank R. Mennickent for his comment on their photometry and the anonymous referee for many kind comments that help us to improve the manuscript. This research has been supported in part by the Grant-in-Aid for Scientific Research (11640226) of the Japanese Ministry of Education, Science, Culture, and Sports.


text:
# Superconducting Bulk Magnets: Very High Trapped Fields and Cracking

0.4pt=0pt

\[

## Abstract

Improved trapped fields are reported for bulk melt-textured $`\mathrm{YBa}_2\mathrm{Cu}_3\mathrm{O}_{7\delta }`$ (YBCO) material in the temperature range between 20 K and 50 K. Trapped fields up to 12.2 T were obtained at 22 K on the surface of single YBCO disks (with Ag and Zn additions). In YBCO mini-magnets, maximum trapped fields of 16 T (at 24 K) and of 11.2 T (at 47 K) were achieved using (Zn + Ag) and Zn additions, respectively. In all cases, the YBCO disks were encapsulated in steel tubes in order to reinforce the material against the large tensile stress acting during the magnetizing process and to avoid cracking. We observed cracking not only during the magnetizing process, but also as a consequence of flux jumps due to thermomagnetic instabilities in the temperature range betweeen 20 K and 30 K.

preprint: Gruss et al. APL
\]

Bulk type II superconductors can trap high magnetic fields that are generated by superconducting persistent currents circulating macroscopically within the superconductor. The main features of the resulting field distribution are a maximum trapped field $`B_0`$ in the center of the superconducting domain and a field gradient towards the sample edge which is determined by the critical current density $`j_\mathrm{c}`$ of the supercurrents. Therefore, high trapped fields $`B_0`$ in bulk superconductors require a high critical current density and a large size of the current loops. Large single-domain bulk $`\mathrm{YBa}_2\mathrm{Cu}_3\mathrm{O}_{7\delta }`$ based (YBCO) samples can be produced by melt texture processing, especially by using $`\mathrm{SmBa}_2\mathrm{Cu}_3\mathrm{O}_{7\delta }`$ (Sm-123) as a seed crystal. The pinning effect in melt textured samples can be improved by irradiation methods and alternatively, by zinc-doping. Cracking of the samples was found to limit the trapped field of bulk YBCO at temperatures below 77 K which can be explained by tensile stresses that occur during the magnetization process due to the stored flux density gradient and may exceed the tensile strength of the material. The mechanical properties of bulk YBCO and its tensile strength can be improved considerably by the addition of Ag. Another possibility to avoid cracking during magnetizing is to encapsulate the bulk YBCO disks in steel tubes. A steel tube leads to stress compensation by generating a compressive stress on YBCO after cooling from 300 K to the measuring temperature which is due to the higher thermal expansion coefficient of steel compared to that of YBCO in the $`a,b`$-plane. Maximum trapped fields $`B_0`$ of 11.5 T at 17 K and of 14.4 T at 22 K were reported for single YBCO disks and mini-magnets consisting of two single YBCO disks, respectively. In both cases, Ag was added to the YBCO disks which were placed into austenite Cr-Ni steel tubes. Previous attempts to combine the two beneficial effects of Ag and Zn additions failed due to the solubility of Zn in Ag in oxygen atmosphere. High trapped fields have also been reported for bulk Sm-123 with additions of Ag. At the surface of Sm-123 disks, trapped fields of 2.1 T at 77 K and of 8 T at 40 K have been observed.

In the present paper, the combined influence of Ag additions (improving the mechanical properties of YBCO) and of Zn doping (improving the pinning properties) on the maximum trapped field of single YBCO disks and of YBCO mini-magnets is studied. Furthermore, we report on attempts to obtain high trapped fields in the temperature range around 50 K using Zn-doped YBCO disks with improved pinning behaviour.

Melt textured YBCO bulk samples were prepared with varying content of Zn \[YBCO(+Zn)\] and with 0.12 wt % Zn combined with 10 wt % Ag \[YBCO(+Ag/Zn)\] using Sm-123 seed crystals. Details of the sample preparation are reported elsewhere. The applied melt process results in large-domain monoliths up to 50 mm in diameter and about 12 mm in height. The $`a,b`$-planes of the large grains were oriented parallel to the surface of the cylindrically shaped samples. The microstructure revealed small, homogeneously distributed Y2BaCuO5 precipitates 200-1000 nm in size and, if Ag was added, supplementary inclusions of solidified Ag droplets in the range between 1 and 10 $`\mu `$m in size.

At 77 K, the disks were characterized by field mapping of the trapped field on the surface of the samples using an axial Hall sensor with an active area of 400 $`\mu `$m2. A maximum trapped field of $`B_0`$ = 1.2 T was achieved in zinc-doped YBCO disks 35 mm in diameter, whereas typical values of $`B_0=0.8`$T were found in undoped YBCO disks of the same size. The improved pinning properties of zinc-doped YBCO were confirmed by magnetization measurements on small samples (diameter 3 mm). The volume density of the pinning force $`f=j_\mathrm{c}B`$ of several zinc-doped YBCO samples determined from magnetization-vs.-field curves at T = 75 K is plotted in Fig. 1 versus the reduced magnetic field $`h=H/H_{\mathrm{irr}}`$ with $`H_{\mathrm{irr}}`$ as the irreversibility field. The maximum pinning force density $`f_{\mathrm{max}}`$ is found to increase with increasing zinc concentration $`x`$ in the range of small concentrations $`x<0.12`$ wt % and to decrease for higher concentrations. The improvement of the pinning force by zinc impurities occupying Cu sites in the CuO2 planes can be attributed to pair breaking by locally induced magnetic moments in the neighborhood of the zinc impurities. The size of these disturbed regions within the CuO2 planes was found to coincide approximately with the coherence length of YBCO and, therefore, they can be considered as effective pinning defects at low Zn concentrations. A strong $`T_\mathrm{c}`$ suppression is observed if these disturbed regions around the Zn impurities begin to overlap for higher Zn concentrations. The subsequent decrease of $`f_{\mathrm{max}}`$ for Zn concentrations $`x>0.12`$ wt % (see Fig. 1) is caused by this $`T_\mathrm{c}`$ suppression. Therefore a Zn concentration of 0.12 wt % was chosen to investigate the temperature dependence of the trapped fields for YBCO(+Zn). The same Zn concentration was also used for the investigated YBCO(+Ag/Zn) disks.

The temperature dependence of the maximum trapped field $`B_0`$ was measured on the top of single YBCO disks as well as in the gap between two disks. The measurements were performed in the variable temperature cryostat of a superconducting 18 T magnet. The cylindrical pellets were cut to diameters of 24 mm for YBCO(+Zn) and 22 mm for YBCO(+Ag/Zn) and embedded in austenite Cr-Ni steel tubes with wall thicknesses of 3 mm and 4 mm, respectively (see Fig. 2). The small gap between the steel tube and the YBCO disks was filled with epoxy using a vacuum impregnation technique. For the investigation of single disks, an axial Hall sensor (AREPOC HHP-SA) was positioned on top of the disks in the centre, whereas a transversal Hall sensor (AREPOC HHP-MP) was placed in the gap between the two disks of the investigated mini-magnets. The size of the gap was 2.2 mm (for the YBCO(+Zn) mini-magnet) and 2.6 mm (for the YBCO(+Ag/Zn) mini-magnet). Additionally, a Cernox temperature sensor was fixed on top of one sample in each case. Both, magnetic field and sample temperature, were continuously recorded in steps of 1s during the measurements. After applying an external field at 100 K, the temperature was reduced with a rate of 2.5 K/min to the measuring temperature. When the thermal equilibrium was reached, the external field was ramped down at 0.1 T/min to minimize heating of the samples by dissipative flux motion. Nevertheless, an increase in sample temperature of 0.2 K at 60 K and 1.8 K at 25 K was observed. The trapped fields $`B_0`$ presented in this paper (see Figs. 3 and 4) correspond to the field value recorded immediately after reaching an external field of zero. Afterwards, the disks were heated up to 100 K in order to remove the remanent field. To avoid thermally induced flux jumps the heating rate was limited to 0.5 K/min below 24 K, to 1 K/min between 24 K and 30 K and to 5 K/min above 30 K.

The temperature dependence of the trapped field $`B_0`$ in YBCO single disks and mini-magnets with different combinations of Ag and Zn is shown in Figs. 3 and 4. The $`B_0`$(T) curves for YBCO(+Zn) show a strong increase of the trapped field with decreasing temperature. Compared with the earlier investigated YBCO disks with 10 wt % Ag addition \[YBCO(+Ag)\], the improved pinning properties of this material result in a considerable shift of the $`B_0`$(T) curves to higher temperatures by about 18 K (single disks) and 15 K (mini-magnets). Maximum $`B_0`$ values of 9 T (at 44 K) and 11.2 T (at 47 K) were obtained in YBCO(+Zn) single disks and mini-magnets, respectively, before cracking occured. The addition of Ag to YBCO was found to result in a reduced density of microcracks in the material and, therefore, increases the tensile strength and the achievable trapped field. The highest trapped fields were obtained for YBCO(+Ag/Zn) disks showing improved mechanical and pinning properties. In particular, trapped fields up to 16 T (at 24 K) were achieved in YBCO(+Ag/Zn) mini-magnets. As expected, the corresponding $`B_0`$(T) curves (for single disks and mini-magnets) lie between those for YBCO(+Zn) and YBCO(+Ag). For the YBCO(+Ag/Zn) disks, cracking was observed at higher $`B_0`$ values than for the YBCO(+Ag) disks. This can be explained by the improved reinforcement of the YBCO(+Ag/Zn) disks for which the wall thickness of the steel tubes was enlarged to 4 mm compared to 2 mm in the case of the YBCO(+Ag) disks.

Mostly, cracking of the YBCO disks was observed during magnetizing. However, especially in the temperature range between 20 K and 30 K the disks were damaged by flux jumps that occur during the heating procedure after the measurement for removing the trapped field. Despite the low rate used to warm up the samples, flux jumps due to thermomagnetic instabilities were observed for high trapped fields. The features of a small self-terminating flux jump are shown in the upper panel of Fig. 5:

This partial flux jump released only a small part of the magnetic energy stored in the YBCO disk. The corresponding increase of the temperature by about 10 K caused no degradation of the YBCO disk. On the other hand, the complete flux jump shown in the lower panel of Fig. 5 caused a rapid change of the entire magnetic energy into thermal energy turning the superconductor to the normal state. This flux jump leads to an irreversible damage of the YBCO disk.

In conclusion, very high magnetic fields up to 16 T were trapped in YBCO(+Ag/Zn) disks embedded in a steel tube. These YBCO bulk magnets with improved mechanical and pinning properties resist the strong tensile stress acting during the magnetizing process on the superconductor, because the tensile stress is partly compensated by the compressive stress exerted by the steel tube. In order to remove the trapped field, the YBCO disks have to be warmed up. Thereby, below 30 K, cracking of the disks may occur due to flux jumps.

The authors are indebted to B. Preuß and P. Bartusch for their help in sample characterization. This work was supported by the German Bundesministerium für Bildung und Forschung under contract #13N7677.


text:
# Comment on ‘Two-body random ensembles: from nuclear spectra to random polynomials’

An investigation of the spectroscopy of even-even nuclei with random one- and two-body interactions has led to the surprising result that the ground state has $`J^P=0^+`$ in typically $`6070`$ $`\%`$ of the cases. This holds both for the nuclear shell model  and the interacting boson model . To understand this phenomenon it was suggested to study the distribution of the lowest eigenvalues . To this end, a technique was developed in  for tridiagonal Hamiltonian matrices, such as for the one- and two-body interactions in the vibron model. It relies on the property that for these matrices it is possible to obtain an estimate for the lowest eigenvalue as the minimum of a function .

The purpose of this comment is to point out that for the vibron model the appropriate function is not a quadratic polynomial as used in , but a more general form which in the vibrational basis to leading order in $`1/N`$, is given by

$`f(z,j)=az^2+bz+c+dj^2+e\sqrt{(1z)^2(z^2j^2)}.`$ (1)
Here $`z=n_p/N`$ and $`j=J/N`$ . In the absence of off-diagonal terms ($`e=0`$), the spectrum is that of a three-dimensional oscillator characterized by $`n_p=0,1,\mathrm{},N`$ and $`J=n_p,n_p2,\mathrm{},1`$ ($`n_p`$ odd) or $`0`$ ($`n_p`$ even). In this case, the function $`f(z,j)`$ reduces to a quadratic polynomial, and the results of  can be used to obtain the percentages for ground states with $`J=0`$, $`J=1`$ and $`J=N`$ to be $`64.6`$, $`5.2`$ and $`19.8`$ $`\%`$, respectively, in close agreement with the exact values of $`64.4`$, $`4.1`$ and $`24.2`$ $`\%`$ obtained for $`N=20`$ and 100000 runs. In the case of only off-diagonal interactions ($`a=b=c=d=0`$), the ground state always has $`J=0`$.

In the general case, both diagonal and off-diagonal terms are present. The off-diagonal term introduces a deformation in the system, which gives rise to the occurrence of vibrational bands with $`J=0,1,2,\mathrm{},J_{\mathrm{max}}`$. As a consequence, the $`J=1`$ ground states tend to disappear in favor of an enhancement of the $`J=0`$ ground states. This effect explains qualitatively the change in the $`J=0`$, $`J=1`$ and $`J=N`$ ground states from $`64.4`$, $`4.1`$ and $`24.2`$ to $`72.1`$, $`0.5`$ and $`23.7`$ $`\%`$, respectively.

In order to carry out a more quantitative analysis, we have developed a simple method to analyze arbitrary functions $`f(z,j)`$ by sampling the two-dimensional space ($`z,j`$). As an example, in Figure 1 we show the percentages of ground states with $`J=0`$ as a function of $`x=\sigma _1/(\sigma _0+\sigma _1)`$, where $`\sigma _i`$ denotes the width of the Gaussian distributions for the diagonal ($`i=0`$) and off-diagonal ($`i=1`$) matrix elements. The case of equal widths (discussed in the previous paragraph) corresponds to $`x=0.5`$. Figure 1 shows a good agreement between the values obtained from an analysis of $`f(z,j)`$ (dashed line) and the exact results (solid line). The discrepancies observed for the quadratic polynomials of  (dotted line) are due to a different treatment of the off-diagonal matrix element.

In conclusion, the approach proposed in  provides an interesting method to analyze random interactions in terms of random polynomials which is especially suited for tridiagonal Hamiltonian matrices. The advantage of this method is that it foregoes the diagonalization of thousands of matrices. In this comment we have presented a careful analysis of the dominance of $`J^P=0^+`$ ground states in the vibron model, and found good agreement with the exact calculations.

For cases of interest to nuclear physics, such as the IBM or the shell model, the Hamiltonian matrix has a more complicated form. Although it can be brought to tridiagonal form, the lack of an analytic form for the function $`f`$ leads to complications comparable to those of bringing it to diagonal form right away. For these cases the method becomes difficult to apply (other than by analogy). Many questions remain open to explain the robust features observed in randomly interacting many-body systems.

This work was supported in part by CONACyT under projects 32416-E and 32397-E, and by DPAGA under project IN106400.