text:
# Simulations of core convection in rotating A-type stars: Differential rotation and overshooting

## 1 INTRODUCTION AND MOTIVATION

Convection within the cores of massive stars has major impact on their structure and evolution, yet little is known about the detailed properties of such convection. These convective motions, driven by the steep temperature gradient arising from the vigorous burning of the CNO cycle, carry outwards a large fraction of the stars’ luminosities. In standard 1–D stellar models (e.g., Maeder & Meynet 2000), the effects of such convection are usually calculated using simple mixing-length prescriptions, but such approaches involve considerable uncertainties. Mixing-length theory can provide only rough estimates of the energy flux carried by the convection, and no effective estimates of the differential rotation or meridional circulation generated by the convective flows, which may have important consequences for mixing and for redistribution of angular momentum within massive stars.

Overshooting from convection zones raises other problems. It has long been realized that convective motions are unlikely to come to a halt at the boundary between the convective core and the stable radiative zone above it. Indeed, upward moving fluid parcels will penetrate into the stable zone, decelerating and mixing with their surroundings (e.g., Roxburgh 1965). Such overshooting motions might bring fresh fuel into the core, thereby prolonging a star’s lifetime on the main sequence, and have noticeable effect on its evolution (e.g., Woo & Demarque 2001). Yet estimating the extent of this overshooting is challenging.

Likewise uncertain is the differential rotation achieved within convective cores. In the solar convection zone, Reynolds stresses, meridional circulations, and viscous forces give rise to a prominent differential rotation, now being probed extensively by helioseismology (Gough & Toomre 1991; Schou et al. 1998). Explaining that differential rotation within the sun has been a major challenge for theory and simulation in recent years, with 3–D spherical shell simulations of turbulent convection now beginning to make contact with the helioseismic results (e.g., Miesch et al. 2000; Elliott et al. 2000; Brun & Toomre 2002). The presence of such non-uniform rotation deep within more massive stars would have major consequences for the properties of such stars: differential rotation can give rise to shear instabilities that stir and mix material, and may serve to build strong magnetic fields through sustained dynamo action.

Indeed, the generation of magnetic fields within stars by dynamo action must result from large-scale convection in the stellar plasma interacting with rotation. In some stars with convective envelopes, such as the sun, the building of orderly magnetic fields and cyclic activity is thought to depend sensitively on highly turbulent convection yielding strong differential rotation, including a tachocline of shear at the interface between the bottom of the convection zone and the radiative interior (e.g., Thompson et al. 2003). More generally, correlations between magnetic fields and rotation have been inferred in main-sequence stars ranging from F5 to M9 (Noyes et al. 1984; Mohanty et al. 2002), with these fields also thought to be a result of dynamo action driven by convection, whether occupying an envelope or the full interior. Yet a comprehensive understanding of how this occurs remains elusive. More massive stars with convective cores are similarly likely to admit magnetic dynamo action, for they too possess the necessary ingredients: intensely turbulent convection, a highly conductive medium, and global effects of rotation, all thought to be crucial for the building of magnetic fields (Charbonneau & MacGregor 2001). Whether the resulting fields are chaotic, or global and orderly, may be a sensitive matter.

### 1.1 Estimates of Overshooting From Convective Cores

Some theoretical work has provided constraints on the size of the convective core and the extent of overshooting. Roxburgh (1978, 1989, 1998; see also Zahn 1991) showed that an upper limit to the total size of the convective region can be deduced by considering the integral

$$_0^{r_c}\left(L_{rad}L\right)\frac{1}{T^2}\frac{dT}{dr}𝑑r=_V\frac{\mathrm{\Phi }}{T}𝑑V>0,$$
(1)
with $`L`$ the total nuclear luminosity, $`L_{rad}`$ the radiative luminosity, $`\mathrm{\Phi }`$ the viscous dissipation per unit volume, and $`r_c`$ the radius of the convective core including overshooting. Since within the convective core some of the nuclear energy must be carried by convection, $`L_{rad}<L`$ there. Thus if viscous dissipation is neglected, there must be an overshooting region where $`L_{rad}>L`$, whose extent can be estimated from this equation. For small convective cores, these considerations yield an overshoot region of size $`d0.18r_0`$, where r0 is the radius of the convectively unstable region (Roxburgh 1992). Rosvick & VandenBerg (1998) found that observations of the stellar cluster NGC 6819 were best fit by overshooting that is about half the upper limit provided by Roxburgh’s constraint.

Various attempts have been made to estimate the extent of overshooting from observations, typically expressed in terms of pressure scale heights. Meynet, Mermilliod, & Maeder (1993), for example, found that best-fit isochrones for a large number of clusters required overshooting that extends for about $`0.20.3`$ times the pressure scale height. Perryman et al. (1998) found evidence for a similar degree of overshooting in Hyades stars. In the open cluster M67, whose main-sequence turnoff stars are thought to be near the low mass end of stars possessing convective cores, some authors have suggested overshooting that is of order $`0.1`$ pressure scale heights (e.g., Maeder & Meynet 1991; Carraro et al. 1994).

### 1.2 Challenges Raised by A-type Stars

Convective cores are realized for main-sequence stars more massive than about $`1.2`$ solar masses, thereby providing many stars in which the effects of overshooting could be assessed. However, A-type stars exhibit a variety of peculiarities that have made them the subject of particularly detailed study (see Wolff 1983, hereafter W83, for a broad review). We recall that A-type stars possess both a convective core as well as multiple shallow convection zones near the surface. Some of these stars display strong abundance anomalies, with greatly enhanced rare-earth element abundances relative to normal stars (e.g., Kurtz 1990). Surface abundance anomalies have also been observed in a variety of other stars (Preston 1974; Gehren 1988; Pinsonneault 1997; Abt 2000). Radiative diffusive separation, wherein radiation pressure drives outward some elements while others sink, is a favored explanation for these abundance features, but requires a very stable radiative zone in which this may be occurring. Latour, Toomre, & Zahn (1981) showed that the H and He envelope convection zones in A-type stars are likely linked by penetrative motions that would upset the delicate quiescence needed for radiative diffusion to work effectively. This problem can be avoided if helium gravitationally settles out to greater depths, in which case the deeper convection zone driven primarily by that element’s second ionization would not exist (e.g., Vauclair, Vauclair, & Pamjatnikh 1974). Such settling is impeded by meridional circulations within the radiative zone (Michaud et al. 1983), which are thought to increase in amplitude with increasing stellar rotational velocity. It has recently been shown (e.g., Richard, Michaud, & Richer 2001) that in some stars a further iron convection zone is established below the other two surface convective regions, which may lead to deeper mixing within such stars (e.g., Vauclair 2003; Richer, Micaud, & Turcotte 2000). While this basic picture explains many of the observed abundance features, many of the details are not precisely known.

A subset of the A-type stars also possesses strong surface magnetic features that appear to persist for many rotation periods (W83). Such fields may well be oblique rotators of primordial origin (e.g., Mestel 1999). However, if dynamo-generated fields within the convective core were able to rise through the radiative zone, possibly by means of magnetic buoyancy instabilities, they might also influence the surface fields. Recent numerical studies (MacGregor & Cassinelli 2003) have provided tantalizing indications that it may indeed be possible for strong magnetic fields generated in the core to rise to the surface.

The A-type stars also hold interest because of the rich set of observational constraints beginning to be provided by asteroseismological probing of such stars. The pulsating Ap stars exhibit high-order nonradial p-modes, which allow some deductions about stellar radii, temperature, and magnetic fields (e.g., Matthews 1991; Cunha 2002).

### 1.3 Modeling Convection in Full Spherical Domains

The extensive observations of A-type stars, and the challenges raised, have encouraged us to consider the core convection influenced by rotation that is occurring deep within their interiors. The surface pathologies of A stars clearly raise many puzzles about their interior dynamics, thereby lending vibrancy to their study. We hope that our modeling will further serve to reveal general dynamical properties of the core convection that is also occuring in other massive stars.

The major uncertainties associated with core convection – the extent of its overshooting into the surrounding radiative envelope, and the differential rotation and circulations it establishes – now lead us to undertake simulations of such convection in full spherical geometries that permit global connectivity. We aim to capture much of the essential physics, admitting highly nonlinear flows that can extend over multiple scale heights as they mold the dynamical structure of the convective core. We realize that there may be few observable consequences at the surfaces of these stars of what may be proceeding dynamically at their centers. However, the circulations and gravity waves that may be induced by the core convection may well have implications at the surface, as might magnetic fields being produced by dynamo action deep within these stars. All may depend somewhat sensitively on the rotation rates of the stars.

Advances in supercomputing are now beginning to allow us to examine the properties of core convection in some detail. We begin here by considering 3–D hydrodynamic simulations of the inner regions of a 2-solar mass A-type star. In subsequent papers, we plan to explore the magnetic dynamo action that may be realized in such convective cores, and to examine possible instabilities that might allow the resulting magnetic fields to rise to the surface.

In Section $`2`$ we describe our formulation of the problem, and briefly summarize the computational tools being used. Section $`3`$ discusses the general properties of the nonlinear convective flows realized in the core, including the transport of energy achieved. The redistribution of angular momentum by the convection yields prominent differential rotation profiles that are presented in Section $`4`$. Section $`5`$ considers the temporal evolution of convective patterns within the core, and examines the gravity waves that are excited within the radiative envelope. Analysis of the penetration and overshooting by the convection is discussed in Section $`6`$. The meridional circulations induced by the convection, within both the core and the radiative envelope, are considered in Section $`7`$. Section $`8`$ evaluates the manner in which the convection yields strong differential rotation. A summary of our principal findings and their implications is presented in Section $`9`$.

## 2 POSING THE PROBLEM

### 2.1 Anelastic Equations

The simulations described here were performed using our Anelastic Spherical Harmonic (ASH) code. ASH solves the three-dimensional anelastic equations of motion in a rotating spherical geometry using a pseudospectral semi-implicit approach (Clune et al. 1999; Miesch et al. 2000). These equations are fully nonlinear in velocity variables and linearized in thermodynamic variables with respect to a spherically symmetric mean state. This mean state is taken to have density $`\overline{\rho }`$, pressure $`\overline{P}`$, temperature $`\overline{T}`$, specific entropy $`\overline{S}`$; perturbations about this mean state are written as $`\rho `$, $`P`$, $`T`$, and $`S`$. Conservation of mass, momentum, and energy in this rotating reference frame are therefore written as

$$\mathbf{}(\overline{\rho }𝐯)=0,$$
(2)
$`\overline{\rho }\left({\displaystyle \frac{𝐯}{t}}+(𝐯\mathbf{})𝐯+2𝛀_𝐨\times 𝐯\right)=\mathbf{}P`$
$`+\rho 𝐠\mathbf{}𝓓[\mathbf{}\overline{P}\overline{\rho }𝐠],`$ (3)
$`\overline{\rho }\overline{T}{\displaystyle \frac{S}{t}}`$ $`=`$ $`\mathbf{}[\kappa _r\overline{\rho }c_p\mathbf{}(\overline{T}+T)+\kappa \overline{\rho }\overline{T}\mathbf{}(\overline{S}+S)]`$
$``$ $`\overline{\rho }\overline{T}𝐯\mathbf{}(\overline{S}+S)+2\overline{\rho }\nu \left[e_{ij}e_{ij}1/3(\mathbf{}𝐯)^2\right]+\overline{\rho }ϵ,`$
where $`c_p`$ is the specific heat at constant pressure, $`𝐯=(v_r,v_\theta ,v_\varphi )`$ is the local velocity in spherical geometry in the rotating frame of constant angular velocity $`𝛀_𝐨`$, $`𝐠`$ is the gravitational acceleration, $`\kappa _r`$ is the radiative diffusivity, $`ϵ`$ the heating rate per unit mass due to nuclear energy generation, and $`𝒟`$ is the viscous stress tensor, with components

$`𝒟_{ij}=2\overline{\rho }\nu [e_{ij}1/3(\mathbf{}𝐯)\delta _{ij}],`$ (5)
where $`e_{ij}`$ is the strain rate tensor. Here $`\nu `$ and $`\kappa `$ are effective eddy diffusivities for vorticity and entropy. To close the set of equations, linearized relations for the thermodynamic fluctuations are taken as

$$\frac{\rho }{\overline{\rho }}=\frac{P}{\overline{P}}\frac{T}{\overline{T}}=\frac{P}{\gamma \overline{P}}\frac{S}{c_p},$$
(6)
assuming the ideal gas law

$$\overline{P}=\overline{\rho }\overline{T},$$
(7)
where $``$ is the gas constant. The effects of compressibility on the convection are taken into account by means of the anelastic approximation, which filters out sound waves that would otherwise severely limit the time steps allowed by the simulation.

Convection in stellar environments occurs over a large range of scales. Numerical simulations cannot, with present computing technology, consider all these scales simultaneously. We therefore seek to resolve the largest scales of the nonlinear flow, which we think are likely to be the dominant players in establishing differential rotation and other mean properties of the core convection zone. We do so within a large-eddy simulation (LES) formulation, which explicitly follows larger scale flows while employing subgrid-scale (SGS) descriptions for the effects of the unresolved motions. Here, those unresolved motions are treated as enhancements to the viscocity and thermal diffusivity ($`\nu `$ and $`\kappa `$), which are thus effective eddy viscosities and diffusivities. For simplicity, we have taken these to be functions of radius alone, and to scale as the inverse of the square root of the mean density. We emphasize that currently tractable simulations are still many orders of magnitude away in parameter space from the highly turbulent conditions likely to be found in real stellar convection zones. These large-eddy simulations should therefore be viewed only as indicators of the properties of the real flows. We are encouraged, however, by the success that similar simulations (e.g., Miesch et al. 2000; Elliott et al. 2000; Brun & Toomre 2002) have enjoyed in matching the detailed observational constraints for the differential rotation within the solar convection zone provided by helioseismology.

### 2.2 Computational Approach

Thermodynamic variables within ASH are expanded in spherical harmonics $`Y_{\mathrm{}}^m(\theta ,\varphi )`$ in the horizontal directions and in Chebyshev polynomials $`T_n(r)`$ in the radial. Spatial resolution is thus uniform everywhere on a sphere when a complete set of spherical harmonics of degree $`\mathrm{}`$ is used, retaining all azimuthal orders $`m`$. We truncate our expansion at degree $`\mathrm{}=\mathrm{}_{max}`$, which is related to the number of latitudinal mesh points $`N_\theta `$ (here $`\mathrm{}_{max}=(2N_\theta 1)/3`$), take $`N_\varphi =2N_\theta `$ latitudinal mesh points, and utilize $`N_r`$ collocation points for the projection onto the Chebyshev polynomials. We employ a stacked Chebyshev representation, wherein the computational domain is split into two regions and separate Chebyshev expansions performed for each. The interface between these two regions, here taken to be the approximate boundary between the convective and radiative zones, is thus treated with greater resolution in order to capture elements of the penetrative convection occurring there. We have taken $`N_r=49+33=82`$ and $`\mathrm{}_{max}=85`$ in this study for the cases considered here, and have used $`\mathrm{}_{max}=170`$ for some companion simulations to verify that the spectral resolution was adequate. An implicit, second-order Crank-Nicholson scheme is used in determining the time evolution of the linear terms, whereas an explicit second-order Adams-Bashforth scheme is employed for the advective and Coriolis terms. The ASH code has been optimized to run efficiently on massively parallel supercomputers such as the IBM SP-3 and the Compaq TCS-1, and has demonstrated excellent scalability on such machines.

### 2.3 Formulation of the Model

The model considered here is intended to be a simplified description of the inner 30% by radius of a main sequence A-type star of 2 solar masses, consisting of the convective core (approximately the inner 15% of the star) and a portion of the overlying radiative zone. Contact is made with a 1–D stellar model (at an age of 500 Myr) for the initial conditions, adopting realistic values for the radiative opacity, density, and temperature. This 1–D model was obtained using the CESAM stellar evolution code (Morel 1997, Brun et al. 2002), which employs the OPAL equation of state and opacities (Iglesias & Rogers 1996; Rogers, Swenson, & Iglesias 1996) and the nuclear cross sections of Adelberger et al. (1998) in describing the microscopic properties of the stellar plasma. Convection is computed within CESAM using a classical mixing length treatment calibrated on solar models, with some convective overshooting (typically $`r=0.25H_p`$) taken into account. Our simulations were initialized using the radial profiles of gravity $`g`$, radiative diffusivity $`\kappa _{rad}`$, mean density $`\overline{\rho }`$, and mean entropy gradient $`d\overline{S}/dr`$ as the starting points for an iterative Newton-Raphson solution for the hydrostatic balance and for the gradients of the thermodynamic variables. The mean temperature $`\overline{T}`$ is then deduced from equation (7). This technique yields mean profiles in reasonable agreement with the 1–D stellar model (Fig. $`1`$). The nuclear source of energy $`ϵ`$ in our simulations is also deduced from the stellar model. For computational convenience, we take $`ϵ=ϵ_0\overline{T}^8`$, determining the constant $`ϵ_0`$ by requiring that the integrated luminosity match the 1–D model’s surface luminosity.

We have softened the steep entropy gradient contrast encountered in going from the convective core to the surrounding radiative zone, which would otherwise favor the driving of small-scale, high-frequency internal gravity waves that are presently unresolvable given reasonable computational resources. This lessening of the ‘stiffness’ of the system has impact on the extent to which convective motions may overshoot into the radiative region. We have excluded the inner 2% of the star from our computational domain for numerical reasons, namely that the coordinate systems employed in ASH possess both coordinate singularities at $`r=0`$ and decreasing meshsize (and hence highly limited time steps) with decreasing radius. It is difficult to gauge with certainty the effects of excluding this central region, which might in principle project both a Taylor column aligned with the rotation axis (e.g., Pedlosky 1987) into the surrounding fluid, as well as give rise to boundary layers. However we have seen no evidence that spurious physical responses have been generated by its omission. In trial simulations with both smaller and larger excluded central regions (ranging from 0.3% to 4% of the stellar radius), the properties of the developed mean flows were very similar to those described here.

As boundary conditions on the spherical domain, we utilize impenetrable and stress-free conditions for the velocity field and constant entropy gradient (i.e constant emergent flux) both at the inner and outer boundaries. As initial conditions, some simulations have been started from quiescent conditions of uniform rotation whose radial stratification is given by the 1–D stellar model; others used evolved states from prior solutions to initiate simulations with reduced diffusivities.

The main parameters of our simulations are summarized in Table 1. In brief, we have modeled the inner regions of 2-solar mass A-type stars at a range of rotation rates, with nonlinear flows of varying complexity achieved by modifying the effective eddy viscosities and diffusivities $`\nu `$ and $`\kappa `$. Several values of $`\nu `$ and $`\kappa `$ are utilized to achieve both laminar and more complex convective flows in some companion cases, all at a Prandtl number $`P_r=\nu /\kappa =0.25`$. We consider rotation rates from one-tenth to four times the solar rate of $`\mathrm{\Omega }_o`$$`=2.6\times 10^6`$ $`s^1`$ (or $`414`$ nHz), though concentrating primarily on models at the solar rotation rate. Three different values of the maximum entropy gradient $`d\overline{S}/dr`$ in the radiative region were taken, thus effectively varying the stiffness of the boundary between that region and the convective core, primarily as a way to study the effects on overshooting motions induced by softening that boundary. We have chosen only representative simulations for summary in this paper (and in Table 1) out of a larger total set of computations.

Figure $`1`$ shows the mean radial profiles for temperature, density, and pressure attained in a representative simulation (case C), together with profiles of the same variables in a typical 1–D stellar structure model. The evolved simulation possesses somewhat less steep gradients for these variables than does the 1–D stellar model, but the functional form of these quantities is largely unchanged by our detailed treatment of the convection.

## 3 NATURE OF CORE CONVECTION

The flows in all the cases studied are highly time dependent, with complex and intermittent features emerging as the simulations evolve. The convection is characterized by large-scale flows that sweep through much of the unstable core, coupling widely separated sites. Thus influences of the convective motions are decidely nonlocal. Such global motions can result because of the topologically connected nature of full spheres. The addition of rotation to these spherical domains admits Coriolis forces and resulting velocity correlations (Reynolds stresses) that yield some surprising effects. In particular, we find that the cores exhibit prominent differential rotation both in radius and latitude, involving in most of our simulations a central column of slow rotation aligned with the rotation axis. Further, the extent of penetration of the convection into the surrounding stable radiative envelope varies from equator to pole, yielding a prolate shape to the central region of nearly adiabatic stratification.

The rich time dependence typical of core convection implies that extended simulation runs must be conducted in order to sample the dynamics. Attaining a statistically equilibrated state in the dynamics typically requires evolving the flows for about $`1500`$ days of physical time, or roughly 60 rotational periods at the nominal solar rate. Deductions about the mean flows accompanying the convection, such as differential rotation and meridional circulation, or of the mean heat transport achieved by the convection, require averaging the flows over lengthy periods of time, typically greater than 100 days in the late stages of a simulation. We begin by describing some of the general attributes of the convection realized in our simulations.

### 3.1 Topology and Scales of Convective Flows

Core convection clearly involves global-scale overturning motions that extend over much of the unstable domain. This is illustrated in Figure $`2`$ which presents both a 3–D rendering of the instantaneous radial velocity near the outer surface of the convective core studied in case E, and an equatorial cut at the same instant through that core domain. Specific regions of upflow and downflow occupy large fractions of the irregular surface (Fig. $`2a`$), and extend in radius over much of that unstable volume (Fig. $`2b`$). Evidently there is global connectivity in the flows over well separated locations, spanning wide ranges in longitude and latitude in both hemispheres. The equatorial cut reveals that the flow patterns involve a prominent azimuthal wavenumber $`m=2`$ component, though many other wavenumbers are certainly also involved in describing these overturning motions. Little asymmetry is apparent between upflows and downflows: radially inward and outward motions occur on roughly the same spatial scales, and possess comparable flow amplitudes (see also Table 2). This appears to be a consequence of the relatively modest radial density contrast (about 2.5) present across the unstable core. We do not here obtain the fast narrow downflows and much broader upflows realized in far more compressible configurations representative of the solar convection zone (Brummell et al. 1998, Brun & Toomre 2002).

The convectively unstable region that results is mildly prolate in shape (Fig. $`2a`$), having greater spatial extent near the poles than near the equator. This non-spheroidal shape (cf. §6.1) for the extent of the vigorous convection is a striking feature shared by all our simulations. The rumpled appearance of $`v_r`$ in this 3–D volume rendering comes from emphasizing the large-scale penetrative structure of such velocity fields. In studying time sequences of such renderings, the interface between the well-mixed interior and the radiative exterior looks much like a throbbing heart, as penetrative motions wax and wane and propagate laterally, launching gravity waves in the process (cf. §5.2). Analysis of the time evolution of equatorial cuts (Fig. $`2b`$) reveals that some flows dive through the central regions of the core to emerge on the other side, with such patterns exhibiting distinctive retrograde motion.

### 3.2 Velocity and Thermal Patterns

We present in Figure $`3`$ global mappings for cases B and C of both the radial velocities $`v_r`$ and temperature fluctuations $`T`$ (relative to the average of temperature on a spherical surface at a given radius). These Mollweide projections for the fields are shown both well within the convective interior (at $`r=0.10R`$, lower panels) and at the transition into the stable exterior ($`r=0.16R`$, upper panels). These mappings reveal that the radial velocity patterns in any given snapshot have a noticeable alignment with the rotation axis. (Such projections exhibit meridian lines as increasingly curved arcs away from the central meridian, which appears as a vertical line; here the lines of constant latitude are indeed parallel.) The more laminar case B (Figs. $`3ad`$, left panels) possesses structures that are largely symmetric about the equator, whereas the lower diffusivities ($`\nu `$ and $`\kappa `$) in case C (Figs. $`3eh`$, right) lead to flows with broken symmetries in the two hemispheres as more intricate flows are realized. In both cases the flow structures at high latitudes are more isotropic and of smaller scale, and somewhat faster flows are realized for case C (see also Table $`2`$).

Radial motions in these simulations are intimately linked to fluctuations in temperature. These buoyantly driven flows should typically involve positive $`T`$ in regions of ascending motion within the unstable core, and negative $`T`$ in descending flows there. Such a positive correlation between $`v_r`$ and $`T`$ is largely realized in the lower panels in Figure $`3`$ for both cases showing the $`r=0.10R`$ depth where the convection is vigorous. The less complex flows in case B exhibit such correlation readily, whereas in case C there are also sites where the sense of the correlation is reversed. The nature of the correlations changes as motions penetrate into the radiative surroundings. The buoyancy braking that is expected within this region of stable stratification, as at the $`r=0.16R`$ depth (upper panels, Fig. $`3`$), arises from the radial velocities and temperature fluctuations becoming anti-correlated. The sense of $`T`$ in both regions is partly controlled by the source term in equation (4) that involves the product of $`v_r`$ with the radial mean entropy gradient. When the latter gradient changes sign in going from an unstable stratification (superadiabatic in the convective core) to stable (subadiabatic) stratification as flows extend upward into the radiative envelope, so does the sense of $`T`$. This correspondence between radial velocities and temperature fluctuations is generally true only in the large: on finer scales, secondary instabilities can drive motions that locally violate these principles. Similarly, the nonlinear advection of temperature fluctuations by motions that twist and turn can lead to unusual combinations of $`v_r`$ and $`T`$ at some sites. The anti-correlation of $`v_r`$ and $`T`$ in the region of penetration is clearest in case $`B`$ involving the simpler flow fields. Because the thermal diffusivity is larger than the viscosity in these simulations (i.e. the Prandtl number $`P_r`$ is less than unity, namely 0.25), thermal structures have a somewhat larger physical scale than do associated features in the radial velocity.

The temperature fluctuations at both depths in Figure $`3`$ have amplitudes of only a few K (compared to the mean temperatures of millions of K found at such depths), which combined with radial velocities of order 50 m s-1 in the bulk of the convective interior and the high heat capacity of the gas, are sufficient to carry outward a good fraction of the emerging flux (cf. §3.3). The radial velocities have plummeted to much smaller values just beyond the edge of the convection zone where the motions are beginning to penetrate upward into the radiative zone.

The convective fields in these simulations can be assessed using probability distribution functions (pdfs). Figure $`4`$ shows pdfs for the radial velocity and temperature fluctuations on a spherical surface within the convective core ($`r=0.10R`$) and the region of overshooting ($`r=0.16R`$) of case C, exhibiting both instantaneous pdfs and time-averaged ones. The distributions of radial velocities are largely symmetric Gaussian functions at both depths. The detailed features in the instantaneous pdfs change as the flows evolve. The instantaneous temperature field distributions (Fig. $`4d,h`$) possess spikes, indicating the presence of large-scale coherent temperature structure; at other instants, the spikes may be replaced by broader plateaus. The time-averaged pdfs for T are nonsymmetric, with broad tails extending to positive values in the region of overshooting and to negative ones in the core. These tails serve to compensate for the off-center most probable values that are positive in the core and negative in the overshooting region. Such complexities in the pdfs are realized in all of our simulations, and thus are at some variance with simpler relations assumed in mixing-length treatments of convection.

### 3.3 Energy Transport

The release of energy by CNO-cyle nuclear burning in the core leads to a progressively increasing local luminosity with radius within the inner parts of our domain. The energy produced in such a distributed fashion gets carried outward primarily by radiative diffusion and mechanical transport of heat (or enthalpy) by the convection. To assess the role convection plays in the overall energy balance, we recall that the rate at which convection transports energy in the radial direction is given by the enthalpy flux

$$F_e(r)=\overline{\rho }c_p\overline{v_rT},$$
(8)
involving a product of the radial velocity with fluctuations in the heat content of the fluid, proportional to $`T`$, and averaged over a spherical surface and in time (indicated by the ‘overbar’). We may compare $`F_e`$ with the flux of energy carried by radiation,

$$F_r(r)=\kappa _r\overline{\rho }c_p(\frac{d\overline{T}}{dr}+\overline{\frac{dT}{dr}}),$$
(9)
with the latter term arising from the mean of radial gradients of temperature fluctuations, which are here negligible. The enthalpy and radiative fluxes largely suffice to carry outwards the energy produced by nuclear reactions within the core, yielding a total flux

$$F_t=\frac{L(r)}{4\pi r^2},$$
(10)
where $`L`$ is the luminosity of the star that increases with radius until the surface luminosity $`L_{}`$ is attained. (The kinetic energy flux, the viscous flux, and the flux carried by unresolved subgrid motions are all small compared to $`F_e`$ and $`F_r`$.)

Figure $`5`$ displays the variation with radius of $`F_r`$ and $`F_e`$, expressed as luminosities, for case C. The enthalpy flux is maximized near the middle of the convective core (at $`r=0.08R`$), serving to carry about 57% of the local stellar luminosity $`L`$ there, with the remainder transported by radiation. Since the core convection establishes a nearly adiabatic stratification (with $`_{ad}10^7`$), the associated temperature gradient that arises serves to specify a radiative flux $`F_r`$ that increases steadily with radius. Such a growing $`F_r`$ forces $`F_e`$ to decrease in the upper half of the unstable core. Beyond the convectively unstable region, at $`\overline{r}_c0.14R`$ (cf. §6.2), the enthalpy flux becomes negative, resulting from the anti-correlation of radial velocity and temperature fluctuations as the motions are braked. The radiative flux there must compensate for that inward-directed enthalpy flux. Our simulations have not been evolved for a sufficient amount of time to allow this compensation to occur fully, since the relevant thermal relaxation time is long compared to other timescales in the problem. Thus the total luminosity in Figure $`5`$ displays a small dip at the interface between the convective and radiative zones, which in real stars (or a fully evolved simulation) would be absent.

The correlations of radial velocity and fluctuating temperature that yield an outward-directed enthalpy flux are achieved in a complicated manner within the bulk of the convective interior. In mapping these variables at one instant in time, it is evident that a number of sites yield an unfavorable correlation, whereas in the main the requisite positive correlation is realized.

The nearly adiabatic stratification achieved by the convection in the bulk of the core is much as in 1–D stellar models employing standard mixing-length treatments. However, our self-consistent 3–D analysis of convective dynamics and its transports also provides estimates of the stratification in the region of overshooting, which are far more difficult to assess using mixing-length approaches.

## 4 RESULTING DIFFERENTIAL ROTATION

Convection occurring within deep spherical domains should be strongly influenced by the rotation of the star when the convective overturning time is of the order or greater than the rotation period. This implies a convective Rossby number $`R_{oc}`$ (cf. Table 1) which is of order unity or smaller. The global-scale convection in our models has an overturning period of about one month, comparable to the rotation period of a star rotating at the solar rate. Thus we expect the convection influenced by Coriolis forces to possess correlations between the components of its velocity field. The resulting Reynolds stresses can serve to continuously redistribute angular momentum, leading to angular velocities within the core that are far from uniform rotation. Further, the connected topology of the convective core appears to play an important role in the prominent differential rotation profiles that are realized in our simulations.

### 4.1 Central Column of Slow or Rapid Rotation

We find that strong angular velocity contrasts are established within the convective cores in all our simulations. For model stars that rotate at least as rapidly as the sun, these involve roughly cylindrical central columns of decidedly slow rotation (retrograde relative to the reference frame), accompanied by equatorial regions of somewhat enhanced rotational velocity (prograde). As shown in Figure $`6`$ (left panels) for the three cases B, C and C4, the time-averaged mean longitudinal velocity $`\widehat{v}_\varphi `$ exhibits some variations along that central column, with those fluctuations becoming more pronounced as the complexity of the convection increases in going from case B to C (cf. Fig. $`3`$). There is likewise greater asymmetry in $`\widehat{v}_\varphi `$ between the two hemispheres (delineated by the equator) in these cases. Such symmetry breaking may be anticipated as the convection becomes increasingly nonlinear. The convection itself is not symmetric about the equator, and thus mean flows associated with it, including the differential rotation, are likely to exhibit differences in the two hemispheres. The columns of slow rotation extend slightly into the radiative envelopes at high latitudes, in keeping with the prolate shape of the convective core. The presence of the central slow cylinder of rotation does not appear to result from the small inner sphere omitted from our computational domain, for we have found similar columns of slowness for several different sizes of central spheres. Figure $`6`$ (right panels) shows the variation of angular velocity $`\widehat{\mathrm{\Omega }}`$ with radius in three latitudinal cuts, providing another assessment of the slowness of the central column and of the speeding up of the equatorial region.

These central columns of slow rotation are largely unexpected, and constitute a striking finding of these simulations. A general result from previous studies of convection in deep shells under strong rotational constraints (with $`R_{oc}`$ small) is that equatorial regions of fast rotation are obtained (e.g., Gilman 1979; Miesch et al. 2000; Brun & Toomre 2002). The conservation of angular momentum requires that other regions must be slowed down, and this often appears as higher latitude regions of slowness. The analog in our spherical domain may be the central column of particularly retrograde $`\widehat{v}_\varphi `$, for this could compensate for the equatorial speeding up (prograde $`\widehat{v}_\varphi `$) of regions with a large moment arm. The continuous redistribution of angular momentum by the convection, involving variously effects of Reynolds stresses, meridional circulations and viscous stresses, will be examined in §8 as we seek to understand the balances needed to account for such central columns of slowness.

We have also considered a more slowly rotating case F, with a rotation rate of $`\mathrm{\Omega }_o/10`$, to examine whether the central column of slowness is still achieved under a weaker rotational influence (with $`R_{oc}`$ large). The prior studies of convection in shells (Gilman 1978) indicated that under such conditions the equatorial regions may exhibit relatively slow rotation due to the modified transport of angular momentum. If this property of shell convection were to carry over to full spherical domains, then it might follow that the overall conservation of angular momentum could yield central regions of fast rotation to compensate for the equatorial slowness. Figure $`7`$ shows that for case F this is indeed realized. The central region of relatively fast rotation is decidedly less columnar than those of slow rotation in Figure $`6`$, probably because the Taylor-Proudman constraint is weakened. We thus suspect that central columns of slowness may be ubiquitous for stars that rotate sufficiently rapidly that $`R_{oc}`$ is less than unity. In the remainder of this paper, we concentrate on our simulations that satisfy this criterion (all but case F), since real A-type stars are generally quite rapidly rotating.

### 4.2 Contrasts in Angular Velocity

We now turn to the time-averaged angular velocity $`\widehat{\mathrm{\Omega }}`$ profiles associated with the differential rotation in our convective cores with $`R_{oc}`$ small, some examples of which are shown as radial cuts at constant latitudes in the right-hand panels of Figure $`6`$ that accompany the views of $`\widehat{v}_\varphi `$. The $`\widehat{\mathrm{\Omega }}`$ profiles in Figure $`6`$ emphasize that most of the radiative envelope rotates at the nominal rate of the reference frame. The variations in $`\widehat{\mathrm{\Omega }}`$ are largely confined to the convective core, though within the overshoot region there are still some weak gradients in $`\widehat{\mathrm{\Omega }}`$ evident. In the simulations of semi-turbulent convection rotating at the solar rate (cases A, C and E, with identical viscosities but increasing ‘stiffness’ in the stratification of the radiative exterior), we obtain angular velocity contrasts $`\mathrm{\Delta }\widehat{\mathrm{\Omega }}/\mathrm{\Omega }_o`$ (from the equator to 60°) of order 60% (see also Table 2). These contrasts appear to be insensitive to the stratification of the radiative zone considered here; likewise the character of the differential rotation throughout the convective core is similar in these three cases. In more laminar simulations at the same rotation rate (cases B and D, involving greater viscosities than their companion cases C and E), though the angular velocity contrast is lessened to about 25%, a central cylinder of slow rotation is still evident. At a four-fold higher rotation rate considered as case C4, $`\mathrm{\Delta }\mathrm{\Omega }`$ increases by about a factor of two, yielding a $`\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega }_o`$ that is reduced to about 33%; a central column of slowness is retained. As we increased the rotation rate in going to case C4 from case C, we sought to maintain roughly the same degree of supercriticality in the convection, which required decreasing the effective viscosities to compensate for the stabilizing effects of rotation (e.g., Chandrasekhar 1961). Without such decreases, the angular velocity contrast was lessened to about 22% in the resulting more laminar flow.

Another assessment of the resulting differential rotation in these simulations is provided by examining the kinetic energies associated variously with the convection (CKE), meridional circulation (MCKE), and differential rotation (DRKE), together with their sum (KE), all relative to the rotating frame. As detailed in Table 2, the CKE is reasonably invariant across all our simulations rotating at the solar rate. In going from the more laminar flows of cases B and D to the semi-turbulent flows in cases C and E, the stronger differential rotation that results is reflected in the greater DRKE, and thus also in the total kinetic energy. The MCKE is likewise increased, but remains a small fraction of KE. Thus the more complex flows more effectively transfer energy – from the vast reservoirs of energy associated with the rotation of the frame and with the internal energy of the star – into KE, which nonetheless remains a small fraction of the total energy in these simulations.

From our limited sampling of rotation rates, we cannot reliably estimate how $`\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega }`$ is likely to vary with $`\mathrm{\Omega }`$ at the much higher rates typical of many A-type stars. However, we suspect that the strong differential rotation obtained in our simulations is likely to be realized even for core convection within more rapidly rotating stars with much stronger rotational constraints.

## 5 EVOLVING CONVECTION AND WAVE EXCITATION

The inherently time-dependent convection involves substantial pattern evolution and propagation within the bulk of the convective core, and the excitation of gravity waves in the surrounding radiative envelope. Some large features in the convection appear, persist, and are sheared and advected by the mean flows that they drive. The radiative zone contains the rippling signature of unremitting gravity waves excited by the convection below, accompanied by weak circulations. We now consider such time variability.

### 5.1 Propagation and Shearing of Convective Structures

The flows in these simulations evolve in a complicated manner. This may be seen in Figure $`8`$, which examines the radial velocity patterns within the convective core (at $`r=0.10R`$) in case C at three instants each separated by seven days in time. Large scale features are recognizable in all three images, but show evolution and propagation. Some orderly downflow structures within the convection persist for long periods of time. Yet they also clearly evolve over relatively short time intervals: major downflow lanes present in Figure $`8a`$, for instance, are recognizable in Figure $`8c`$, but have changed in shape and intensity. They do not evolve in isolation: patterns merge, are sheared, and cleave into smaller patterns, as they interact with surrounding flows. Such features are not simply advected by the mean motion of the region in which they are embedded, but can additionally propagate relative to the differential rotation that they establish. The rate and sense of pattern propagation is most readily seen by turning to time-longitude mappings, such as in sampling the radial velocities at 65° north (Fig. $`8d`$) and at the equator (Fig. $`8e`$) over nearly 200 days. In these mappings, persistent strong downflow lanes are revealed as dark bands tilted to the right in Figure $`8e`$, and sheared to the left in Figure $`8d`$. Such tilts reveal that patterns at this radius and near the equator propagate in a prograde sense, taking about 300 days to complete one revolution, and at the higher latitude their propagation is retrograde with a period of about 60 days. These mappings also emphasize that the convection possesses larger scales at the equator than at high latitudes, with a distinctive change in the sense of pattern propagation relative to the frame at mid-latitudes.

### 5.2 Internal Gravity Wave Responses

The vigorous core convection extends upward into the stable stratification that bounds it, with such pummeling serving to excite internal gravity waves within the radiative envelope. Such wave responses are shown vividly in the upper panels of Figure $`9`$, which displays the radial velocities at multiple depths in cases E and C4. These can appear both as localized ripples that are one signature of the waves (Fig. $`9a`$), arising either from distinctive events or from the interference of many waves, or as striking large-scale global resonances (Fig. $`9d`$) that persist (once established) for as many rotation periods as we continued the simulations. The latter response may be more easily realized in our simulations, with their artificially imposed outer boundary conditions that reflect gravity waves, than in real stars where the radiative envelope is extensive and probably has no simple reflecting surfaces. However, they may also mimic gravity waves of moderate radial order $`n`$ (such as $`n=8`$ or $`10`$) in a real star, though the stable stratifications in such stars are much stiffer than the ones we consider here.

The generation of gravity waves by the upward-directed plumes is fully expected. Waves with frequencies less than or equal to the Brunt-Vaisala frequency $`N`$ (proportional to the square root of the local entropy gradient $`d\overline{S}/dr`$ of the stable stratification) are able to propagate in the radiative outer envelope and may be excited in abundance by the time-varying convection, as in case E at $`r=0.26R`$ (Fig. $`9a`$). This also implies that the spectrum will shift as the stiffness of the surrounding envelope is varied. More surprising is the large-scale global resonance of an internal gravity mode (of spherical harmonic degree $`\mathrm{}=3`$ and azimuthal order $`m=2`$) that is readily achieved within the more rapidly rotating case C4, and likewise within some of its more laminar relatives. This may come about because the rotating convection is likely to possess a range of inertial oscillations with frequencies less than 2$`\mathrm{\Omega }_o`$ (Greenspan 1969), which accompany the broad band of lower frequencies present in the time-dependent convection. When there is a suitable overlap between the temporal spectrum of the rotating convection and the admissible internal gravity wave frequencies that depend on $`\mathrm{}`$ and radial order $`n`$ (e.g., Unno et al. 1979; Gough 1993; Dintrans, Rieutord, & Valdettaro 1999), a prominent wave response may result. In our models $`N`$ varies somewhat with radius; for case C4 its representative value is $`1.1\times 10^5`$ Hz, whereas $`\mathrm{\Omega }_o`$ is $`1.04\times 10^5`$ Hz. It appears that a global gravity wave (with low $`\mathrm{}`$) has been selected that propagates retrograde relative to the frame (with a period of about $`40`$ days for a complete revolution), within a region that exhibits almost no differential rotation. In addition to the global mode, other gravity waves appear as ripples near the equator in Figure $`9d`$.

## 6 CONVECTIVE PENETRATION AND RADIATIVE ENVELOPE

Convective motions do not stop abruptly at the interface between the convective core and the stable radiative envelope. Rather, outward moving parcels of fluid possess momentum, and therefore must travel some distance into the stable radiative zone before they come to a halt. There are two effects. Firstly, if convective motions are sufficiently vigorous, and possess a large enough filling factor, they can create a nearly-adiabatic stratification within a portion of the overlying radiative zone that in a 1–D stellar structure model would be stably stratified. This may be termed a region of *penetration*, as discussed in Zahn (1991), Hurlburt et al. (1994), and Brummell, Clune, & Toomre (2002). Secondly, there will be a further region of stable stratification into which convective motions extend, but whose entropy gradient they fail to modify significantly. Within such a region of *overshooting*, outward moving parcels of fluid are braked until they come to a halt. Throughout the domain of penetrative convection, and the extended region of overshooting, the fluid motions should be effective in mixing the chemical composition.

### 6.1 Shape of Convective Core

In the 1–D stellar models used for our initial conditions, the region of unstable stratification is of course spherical. In the evolved simulations, penetrative convective motions have yielded a nearly adiabatically stratified core region (weakly superadibatic) that is now prolate in shape and aligned with the rotation axis. The radius $`r_c`$ of this core region thus varies with latitude, as sketched in Figure $`10`$. The simulations further reveal that the overshooting motions extend outward to a radius $`r_o`$ that is largely independent of latitude, and thus roughly spherical in shape, as also sketched.

Such geometries for the penetrative core and the further region of overshooting are shown in Figure $`11`$ for case C. These regions are delineated by turning to the time- and longitude-averaged enthalpy flux, given by

$$\widehat{F_e}(r,\theta )=\overline{\rho }c_p<v_rT>_{t,\varphi },$$
(11)
where the brackets and their subscripts indicate the averaging. Within the convective core that includes the region of penetration, the enthalpy flux is typically positive, owing to the positive correlations of radial velocity and temperature fluctuations that are largely realized there (cf. §3.2). As the motions extend outward into the subadiabatic stratification of the region of overshoot, the enthalpy flux becomes negative since the temperature fluctuations typically flip their sign almost immediately, becoming in the main anti-correlated with the radial velocity. Thus the inner boundary of the region of overshoot is usefully also revealed by the radius where $`\widehat{F_e}`$ changes sign.

The outer boundary $`r_o`$ of that overshooting region is taken to be the radius at which the enthalpy flux has attained an amplitude of 10% of its most negative value. We introduce this definition because small negative values of the enthalpy flux persist wherever any radial motions anti-correlated with temperature fluctuations exist, even if these radial motions are very weak. Our intent is to identify regions where overshooting motions are reasonably vigorous, recognizing that such a constraint on the amplitude of the enthalpy flux is somewhat arbitrary. For case C, $`r_o`$ has a value of about $`0.158R`$, independent of latitude, whereas $`r_c`$ ranges from about $`0.131R`$ at the equator to $`0.148R`$ at the poles (Fig. $`11`$). The values for other cases are presented in Table 2. Since in all our cases the inner boundary $`r_c`$ of the overshoot region is a function of latitude and the outer boundary $`r_o`$ is largely not, the extent of the overshooting $`d_o`$ varies with latitude, possessing the greatest spatial extent near the equator. Figure $`11`$ also shows that the enthalpy flux $`\widehat{F_e}`$ has prominent variations with latitude and radius within the convective core, with a greater flux at mid and high latitudes than at the equator. As the convection evolves, the time-averaged $`\widehat{F_e}`$ as determined over different epochs shows some changes in its structure and its symmetries with respect to the equator.

The prolate geometry of the convective core may result from the effects of rotation yielding convective flows that typically involve greater radial velocities at high latitudes than at the equator. Figure $`12`$ shows the rms radial velocity in case C within both the convective core ($`r=0.10R`$) and the region of overshooting ($`r=0.15R`$). It is evident that radial velocities are generally lowest near the equator and highest near the poles. Outward-directed motions at low latitudes are turned by the Coriolis forces, thus decreasing their radial velocity component; radial motions that are closely aligned with the rotation axis, such as at high latitudes, are deflected less, and thus retain most of their radial velocity.

### 6.2 Stable Stratification and Overshooting

As convective motions enter the stable stratification of the radiative envelope, they are buoyantly braked (cf. §3.2). The amount of braking they experience is controlled by the magnitude of the stable entropy gradient. By varying that gradient in our simulations, we have thus effectively altered the ability of flows to overshoot into the radiative envelope: a larger $`d\overline{S}/dr`$ implies stronger buoyancy braking, and thus a smaller extent to the overshooting region. A larger $`d\overline{S}/dr`$ also leads to a smaller extent to the penetration as indicated by $`r_c`$. Both these effects are apparent in Figure $`13`$, which displays $`F_e(r)`$ at several latitudes, for three simulations of varying $`d\overline{S}/dr`$. With increasing stiffness comes decreasing $`r_o`$ and $`r_c`$.

For our stiffest and most complex case E, the latitude-averaged extent of overshoot, $`\overline{d_o}`$, is $`1.906\times 10^9`$ cm, or $`0.208H_p`$, where $`H_p`$ is the local pressure scale height averaged over the region of overshoot. In the less stiff case C, a larger $`\overline{d_o}`$ is realized; such measures are quoted for all cases in Table 2. Increasingly complex flows, as achieved in our simulations by lowering the viscosity in going from case B to case C (or from case D to case E, at a different stiffness), yield smaller $`\overline{d_o}`$ than do their more laminar relatives. This seems at first paradoxical, given the larger rms velocities associated with the increased complexity of the convection. However, this trend is largely in keeping with nonlinear modeling of 3–D penetrative compressible convection in localized planar domains (cf. Brummell et al. 2002), which showed that a smaller filling factor for the plumes as realized in the more turbulent convection translates into a smaller extent to the overshooting.

Increasing $`\mathrm{\Omega }_o`$ without altering the viscosity leads to more laminar convection, due to the stabilizing effects of rotation. This in turn yields a greater overshooting. In devising case C4 with its fourfold greater rotation rate, we reduced the viscosity relative to case C in order to achieve comparable complexity in the convection as measured by the Reynolds numbers. In case C4, $`\overline{d_o}`$ is greater than that in the slower rotating case C, suggesting that increasing the rotation rate enhances somewhat the overshooting.

## 7 INDUCED MERIDIONAL CIRCULATIONS

Prominent meridional circulations are induced by the convection within the unstable core, accompanied by much weaker circulations in the radiative surroundings. Within the convective core, such mean meridional flows arise from a combination of buoyancy forces, Reynolds stresses, latitudinal pressure gradients, and of Coriolis forces acting on the mean zonal flow (differential rotation). In the radiative exterior beyond the zone of overshooting, slow circulations can be driven both by viscous stresses that couple that region to the flows in the convective interior, and by latitudinal temperature gradients arising from the nonuniform heating from below. Our simulations reveal that the resulting meridional circulations are complex in structure, possess multiple cells in both radius and latitude, and exhibit noticeable variations in time. Such multi-cellular patterns are in sharp contrast to the single monolithic cells often assumed in mean-field models of differential rotation (e.g., Rekowski & Rüdiger 1998; Durney 2000). On the other hand, our models include viscous stresses that may be largely absent in real stars, and thus our estimates of the meridional circulations in the immediate surroundings of the convective core are likely to be somewhat unrealistic.

### 7.1 Patterns Within Convective Core

The time-averaged meridional circulations achieved in three of our simulations are shown in Figure $`14`$ as streamfunctions of mass flux (cf. equation  of Miesch et al. 2000). Within the convective core (left panels), these flows have amplitudes of order $`5`$ m s-1; they are thus appreciably smaller than the instantaneous radial velocities in the same simulations (cf. Table 2). Indeed, the kinetic energy of these axisymmetric meridional flows (MCKE) is in all cases at least an order of magnitude smaller than that of the non-axisymmetric convection (CKE) and of the differential rotation (DRKE). More rapid flows, and more intricate cellular patterns, are realized in the less viscous cases C and C4 than in case B.

The patterns of meridional flow are quite intricate, with multi-celled structures achieved in radius and latitude. These structures exhibit some alignment with the rotation axis, particularly in the more rapidly rotating case C4. As the convection itself evolves (Fig. $`8`$), so do the mean meridional circulations associated with it. Thus Figure $`14`$ represents typical examples of these time-dependent circulations in a particular epoch. These circulations occupy most of the volume delineated by $`r_c`$, and weaken in the region of overshooting. Thus their extent varies with the stiffness of the exterior stratification.

### 7.2 Weak External Circulations

The meridional circulations achieved in the radiative envelope within our models are much weaker than those in the core. As shown in the right panels of Figure $`14`$, these flows are of large scale and possess amplitudes that in cases B and C are about $`0.06`$ m s-1. The more rapidly rotating case C4 exhibits faster exterior circulations. Such flows evolve more slowly than do the meridional circulations within the core. They are driven by a combination of viscous stresses communicated across these domains (with a characteristic diffusion timescale across the radiative envelope of about $`1.8`$y for case C) and by weak latitudinal thermal gradients resulting from the energy deposited by the latitudinally-varying $`\widehat{F_e}`$ (Fig. $`11`$).

These mean flows within the radiative envelopes in our models are considerably greater in amplitude than estimates of Eddington-Sweet circulations. The latter, which arise from baroclinic effects due to the flattened shape of a rotating star, would involve flows of order $`10^4`$ m s-1 in our A-type stars rotating at $`\mathrm{\Omega }_o`$ (see review in Pinsonneault 1997). Our modeling ignores centrifugal acceleration and thus such distortions, and further imposes a spherical outer boundary upon which a uniform $`dS/dr`$ is specified. Therefore our model does not accurately address the very slow circulations of the Eddington-Sweet variety.

## 8 INTERPRETING THE ROTATION PROFILES

Our simplified model of rotating core convection coupled by penetration to a radiative envelope yields a nonlinear dynamical system with complex feedbacks. The rich variety of responses and marked time variability that is exhibited cannot now be predicted or explained from first principles. However, the numerical simulations have the great advantage that they provide full snapshots that can be assessed for the dynamical balances that are achieved. This can yield insights into the operation of these complicated dynamical systems. Such an approach is particularly useful in identifying the processes involved in continuously redistributing the angular momentum to establish the differential rotation. Because we have adopted stress-free boundary conditions at the top and bottom boundaries of our simulations, no net torque is applied to these rotating spheres of convection. Thus total angular momentum within our simulations is conserved, and we now examine the manner in which it is redistributed.

### 8.1 Redistributing the Angular Momentum

We may assess the transport mechanisms of angular momentum in our simulations, which must combine to give rise to the observed differential rotation, in the manner of Elliott et al. (2000) and Brun & Toomre (2002). We consider the mean radial ($`_r`$) and latitudinal ($`_\theta `$) angular momentum fluxes, and write the $`\varphi `$-component of the momentum equation, expressed in conservative form and averaged in time and longitude, as

$$\frac{1}{r^2}\frac{(r^2_r)}{r}+\frac{1}{r\mathrm{sin}\theta }\frac{(\mathrm{sin}\theta _\theta )}{\theta }=0,$$
(12)
involving the mean radial angular momentum flux

$$_r=\widehat{\rho }r\mathrm{sin}\theta [\nu r\frac{}{r}\left(\frac{\widehat{v}_\varphi }{r}\right)+\widehat{v_r^{^{}}v_\varphi ^{^{}}}+\widehat{v}_r(\widehat{v}_\varphi +\mathrm{\Omega }r\mathrm{sin}\theta )]$$
(13)
and the mean latitudinal angular momentum flux

$$_\theta =\widehat{\rho }r\mathrm{sin}\theta [\nu \frac{\mathrm{sin}\theta }{r}\frac{}{\theta }\left(\frac{\widehat{v}_\varphi }{\mathrm{sin}\theta }\right)+\widehat{v_\theta ^{^{}}v_\varphi ^{^{}}}+\widehat{v}_\theta (\widehat{v}_\varphi +\mathrm{\Omega }r\mathrm{sin}\theta )].$$
(14)
In the above expressions for both fluxes, the first terms in each bracket are related to the angular momentum flux due to viscous transport (which we denote as $`_{r,V}`$ and $`_{\theta ,V}`$), the second term to the transport due to Reynolds stresses ($`_{r,R}`$ and $`_{\theta ,R}`$) and the third term to the transport by the meridional circulation ($`_{r,M}`$ and $`_{\theta ,M}`$). The Reynolds stresses in the above expressions arise through correlations of the components of velocity.

In Figure $`15`$, we show the components of $`_r`$ and $`_\theta `$ for cases $`D`$ and $`E`$ with a strong rotational influence, integrated along co-latitude and radius respectively to deduce the net fluxes through shells at various radii and through cones at various latitudes, such that

$$I__r(r)=_0^\pi _r(r,\theta )r^2\mathrm{sin}\theta d\theta \text{ , }I__\theta (\theta )=_{r_{bot}}^{r_{top}}_\theta (r,\theta )r\mathrm{sin}\theta dr,$$
(15)
We have identified the contributions from viscous (V), Reynolds stresses (R) and meridional circulation (M) terms. This representation is helpful in considering the sense and amplitude of the transport of angular momentum within the convective shells by each component of $`_r`$ and $`_\theta `$.

Turning first to the radial fluxes of angular momentum (Fig. $`15a`$, $`b`$), we see that viscous forces act to transport angular momentum radially inward in both cases. In this they are opposed by a combination of the Reynolds stress and meridional circulation fluxes, to yield a total radial angular momentum flux that is nearly zero, as noted in Fig. $`15`$ by the solid line. While the systems here are highly variable in time, by allowing the system to evolve for extended periods of time (typically thousands of days) and performing long time averages, we appear to be sensing the equilibrated states reasonably well. In going from the more laminar flows of case $`D`$ to the complex ones of case $`E`$, we see that the viscous flux has dropped and that the Reynolds stresses and meridional circulations have changed accordingly to maintain equilibrium. This includes the development of an appreciable integrated $`_{r,M}`$ that is directed radially inward, near the boundary of the radiative and convective zones; this is presumably a result of the stronger circulations and shear generated near that boundary by the more complex flows of case $`E`$.

Examining the latitudinal transport of angular momentum (Fig. 15$`c`$, $`d`$), we note that the effect of the Reynolds stresses in both cases is primarily to speed up the equator, since $`_{\theta ,R}`$ is positive in the northern hemisphere and negative in the southern. It is opposed by the meridional circulation and viscous fluxes, which act to speed up the poles. The manner in which each of the different components of $`_\theta `$ acts does not appear to vary appreciably in going from case $`D`$ to $`E`$. As the level of complexity is increased, however, we see a decrease in the magnitude of all components of $`_\theta `$.

We conclude that for these cases involving a strong rotational constraint, the Reynolds stresses act latitudinally to speed up the equator, and radially to slow down the convective core. These latidudinal fluxes are opposed by the meridional circulations and viscous effects, whereas the radial fluxes are aided somewhat by the meridional circulations. Examining these fluxes gives us a picture of how the angular momentum is continuously transported, and which are the significant processes serving to speed up or slow down certain regions. However, it is difficult to infer from the fluxes alone what will be the detailed $`\mathrm{\Omega }`$ profile that results.

### 8.2 Baroclinicity and Thermal Winds

Rotating convection involves both radial and latitudinal heat transport, with the likelihood that latitudinal gradients in temperature and entropy may result within the convective core. This implies that surfaces of mean pressure and density will not coincide, thereby yielding baroclinic terms in the vorticity equations (Pedlosky 1987; Zahn 1992). Under sufficiently strong rotational constraints, a ‘thermal wind balance’ might be achieved in which departures of the angular velocity from being constant on cylinders (aligned with the rotation axis) are controlled by those baroclinic terms. Indeed, some mean-field approaches have invoked such a balance to obtain differential rotation profiles with bearing on the solar convection zone (e.g., Kichatinov & Rüdiger 1995). As discussed in Brun & Toomre (2002), such a balance effectively implies that

$$\frac{\widehat{v}_\varphi }{z}=\frac{1}{2\mathrm{\Omega }_o\widehat{\rho }c_P}\widehat{S}\mathbf{\times }\widehat{p}|_\varphi =\frac{g}{2\mathrm{\Omega }_orc_P}\frac{\widehat{S}}{\theta }\text{ ,}$$
(16)
where $`z`$ is parallel to the rotation axis. Thus latitudinal entropy gradients could serve to break the Taylor-Proudman constraint, which would otherwise require the rotation to be constant on cylinders with $`\widehat{v}_\varphi /z`$ equal to zero. That constraint may also be broken by Reynolds and viscous stresses. We next show that these terms are at least as important as the baroclinic terms in establishing the differential rotation in the bulk of the convective core.

Figure $`16`$ assesses for case E the extent to which latitudinal entropy gradients serve to drive the temporal mean zonal flows $`\widehat{v}_\varphi `$ seen as the differential rotation. Figure $`16b`$ displays the gradient $`\widehat{v}_\varphi /z`$ of those zonal flows, and Figure $`16c`$ the corresponding latitudinal entropy term on the right hand side of equation (16), which in an exact thermal wind balance would be identical. Figure $`16d`$ shows the difference between this baroclinic term and the actual $`\widehat{v}_\varphi /z`$, thus providing a measure of departures from such a balance. Within the bulk of the convective core, Reynolds and viscous stresses lead to zonal flows that are not so balanced. At the interface between the core and the radiative exterior, however, baroclinicity accounts for most of the differential rotation.

Thus examining the nature of this thermal wind balance, combined with the assessment of fluxes of angular momentum, reveals that Reynolds and viscous stresses have a major role in establishing the differential rotation in the bulk of the core, with latitudinal thermal gradients influencing the profile near the interface between the convective region and the radiative envelope.

## 9 SUMMARY AND PERSPECTIVES

Our highly simplified treatment of core convection within A-type stars has begun to capture some of the complex dynamics that must be occurring when turbulent convection and rotation meet in a full spherical domain. The most severe simplifications in our models are that we are able to simulate only the central portions of these stars, that compositional gradients that must be present are ignored, and that the levels of turbulence studied here are only modest compared to what may be realized in real stars. The dynamics may well be influenced by each in turn, and in ways that are difficult to predict at this stage. Though our work is thus quite preliminary, it has revealed several dynamical properties that may well turn out to be robust.

Firstly, these convective cores rotate differentially. Prominent angular velocity contrasts are established in all cases studied here, with a central column of slow rotation realized in stars that are strongly influenced by rotation (with $`R_{oc}`$ small). Within the bulk of the convection zone, this differential rotation is driven primarily by Reynolds and viscous stresses, whereas near the interface between the core and the radiative envelope, baroclinicity plays an important role.

Secondly, convective motions penetrate into the radiative envelope, yielding a prolate shape to the region of nearly adiabatic stratification. The surrounding region of overshooting, in which convective motions persist but do not appreciably modify the prevailing entropy gradient, is found to have an outer boundary that is largely independent of latitude, yielding a spherical shape to the domain experiencing convective stirring.

Thirdly, small-amplitude circulations and gravity waves are excited within the radiative envelope. The latter appear in some simulations as a rippling patchwork of waves, and in our more rapidly rotating case C4 as a striking global-scale resonance that once established, appears to persist.

The simplifications we have made in our modeling are likely to impact the penetrative properties to some degree. In particular, the $`\mu `$-gradients of chemical composition present in real stars should diminish the extent of penetrative and overshooting motions realized there (Zahn 1991). They also probably reduce the vigor of circulations within the radiative envelope (Lebreton & Maeder 1987). More turbulent flows are likely to yield smaller spatial filling factors for the convective plumes, and hence probably also lessen the extent of overshooting (Brummell et al. 2002). Further, real stars possess stiffer entropy stratifications in their radiative envelopes than those we have been able to consider. This too would lead to less overshooting. Thus our estimates here provide upper limits for the extent of penetration and overshooting that may be achieved in real stars.

The strong regions of shear realized in all our simulations may contribute to the building of orderly magnetic fields by dynamo action. Poloidal fields can be converted into toroidal ones by such differential rotation (the $`\mathrm{\Omega }`$-effect), with regeneration of the poloidal fields perhaps accomplished by the helical turbulence (the $`\alpha `$-effect; cf. Mestel 1999). Acting together, these processes may build well-organized fields that could, if subject to magnetic buoyancy instabilities, rise to the surfaces of these stars where they might be observed (MacGregor & Cassinelli 2003). We are in the process of beginning to explore the generation of magnetic fields within these cores (cf. Browning, Brun, & Toomre 2003).

The penetrative convection and overshooting in our simulation may likewise have some observable consequences. Although we have not investigated in any detail the bringing of fresh fuel into the convective core by these motions, some replenishment of the hydrogen stockpile available for nuclear reactions must occur, and thereby lead to a prolonged lifetime on the main-sequence.

We thank Douglas Gough and Jean-Paul Zahn for very helpful discussions. J.T. thanks the Observatory of Paris, Meudon and the CEA-Saclay for their hospitality during the interpretation of these simulations. This work was partly supported by NASA through SEC Theory Program grant NAG5-11879, and through the Graduate Student Researchers Program (NGT5-50416). Various phases of the simulations here were carried out with NSF PACI support of the San Diego Supercomputing Center (SDSC), the National Center for Supercomputing Applications (NCSA), and the Pittsburgh Supercomputing Center (PSC). Much of the analyses of the extensive data sets were conducted in the Laboratory for Computational Dynamics (LCD) within JILA.


format: text

text:
# OGLE small amplitude red giant variables in the Galactic Bar

## 1 Introduction

While the existence of the Galactic Bar was pointed out decades ago (de Vaucouleurs, 1964), the Bar existence received a broader support only following the work of Blitz & Spergel (1991), based on 2.4 micron light distribution in the inner Galaxy as measured by Matsumoto et al. (1982), and the work of Nakada et al. (1991), based on the asymmetry in the distribution of Galactic Bulge stars as measured by IRAS. The presence of the Galactic Bar has been observed in many different ways; recent reviews are provided by Garzón (1999), Gerhard (2002), Dehnen (2002), and Merrifield (2003).

The Optical Gravitational Lensing Experiment (OGLE) is a long term observing project carried out at the Las Campanas Observatory in Chile, operated by the Carnegie Institution of Washington. The instrumentation of the three stages of the project was described by Udalski et al. (1992, OGLE–I), Udalski et al. (1997, OGLE–II), and Udalski et al. (2002a, OGLE–III). A catalogue of over 200,000 variable stars discovered within approximately 11 square degrees in the direction of the Galactic Bulge was published by Woźniak et al. (2001). The public domain data set contains $`200400`$ I-band photometric data points obtained in the observing seasons 1997, 1998, 1999 with the OGLE-II instrumentation: a dedicated 1.3-meter telescope with a $`2K\times 2K`$ pixel CCD camera, built and operated by the Warsaw University Observatory.

The presence of the Galactic Bar was already detected in the apparent distribution of red clump giants in the OGLE-I data (Stanek et al. 1994, Stanek et al. 1997). In this paper, we set out to determine whether the evidence for the Galactic Bar may also be found among OGLE variables. We calculated periodograms (Lomb 1976, Press et al. 1992) for each of the 200,000 objects in the OGLE-II variable star catalogue (Woźniak et al., 2001) that had at least 100 I-band measurements, and accepted as real those for which the probability of a false signal (Press et al., 1992) was smaller than 0.001. Because fields 45 and 46 have only two seasons of observations (fewer than 100 measurements), periodograms were not calculated for any stars in these fields. We were struck by the two clear sequences apparent in the period – amplitude diagram, presented in Fig. 1. For many stars, several significant periods were found. The plotted period has the largest amplitude. Our convention is that amplitudes are defined as half of the peak-to-peak amplitude for a given mode.

Originally, we had intended to use Mira-type long period variables as possible tracers of the Galactic Bar. However, it is clear that the variables defining the two sequences in Fig. 1 are more numerous and, having small amplitudes, they appear to be more useful than Miras as tracers of the Galactic structure. Therefore, this paper is devoted to these low amplitude variables only.

In addition to the two sequences present in Fig. 1 there is a clump of variables around $`\mathrm{log}P2.1`$ and $`\mathrm{log}A<2`$; this is mostly caused by an artefact, and should be ignored.

We divided the variables into types A and B, with the division line shown in Fig. 1. Type A is above and to the left, and type B below and to the right, of the division line. Variables within the parallelograms:

$$\begin{array}{ccccc}2.31& <& \mathrm{log}(A)& <& 1.31,\hfill \\ 6.56& <& \mathrm{log}(A)3.223\mathrm{log}(P)& <& 5.66,\hfill \end{array}$$
$`(1a)`$
$$\begin{array}{ccccc}1.82& <& \mathrm{log}(A)& <& 0.90,\hfill \\ 7.62& <& \mathrm{log}(A)3.223\mathrm{log}(P)& <& 6.56,\hfill \end{array}$$
$`(1b)`$
were selected in 33 OGLE-II fields spanning a large range of Galactic longitudes, but only a small range of Galactic latitudes. The colour – magnitude diagram, corrected for interstellar reddening and extinction with the maps provided by Sumi (2003), is shown in Fig. 2 for the field BUL\_SC1. Type A variables are marked with large blue dots; B variables are marked with green open circles. A relatively small range of apparent I-band magnitudes and a large number of variables indicates that these stars might be useful ‘standard candles,’ and that it might be possible to use them as tracers of the Galactic Bar.

Most variables shown in Fig. 2 are bright red giants. It is very likely that these stars are similar to the variables found by Kiss & Bedding (2003) and Ita et al. (2003) in the OGLE catalogue for the LMC and SMC, to those found by Cook et al. (1997), Wood et al. (1999), and Wood (2000) in the MACHO database for the LMC, to the M giant variables in NGC 6522 (Glass & Schultheis, 2003), and to those found by Cook et al. (1997) in the MACHO database for the Galactic Bulge. Minimum variability levels have been detected for red giants from the ESO project “long term photometry on variables” (Jorissen et al. 1997), and also from Hipparcos photometry (Eyer & Grenon 1997), the cooler stars having larger minimum variability thresholds and larger mean amplitudes. Furthermore, a period-amplitude relation has been noticed for a few red giants stars by Jorissen et al. (1997), and also in a large sample of stars from ASAS data (Eyer & Blake 2002). A possible theoretical mechanism for the oscillations was studied by Dziembowski et al. (2001) and Christensen-Dalsgaard et al. (2001).

The aim of our paper is to select from the OGLE catalogue red giant variables populating the two period – amplitude sequences in Fig. 1, to provide many parameters for these variables in electronically accessible tables, and to check if they trace the Galactic Bar.

## 2 Selection of variables

We used reddening maps published by Sumi (2003) for 48 OGLE-II Galactic Bulge fields. All variables of type A in a subset of these fields located at the Galactic latitude $`b4^{}`$ are shown in Fig. 3. It is clear that they form three distinct groups: luminous red giants, less luminous red giants located near the Bulge red clump region, and relatively blue Galactic disc stars. These three groups are well separated. We selected stars located above and to the right of the line in Fig. 3 for our study, and we refer to them as OGLE Small Amplitude Red Giant variables (OSARGs).

Some type B variables are close to the saturation limit of OGLE-II, and therefore we restrict our study of the Galactic structure to OSARG type A variables only. The OGLE saturation limit causes a problem for the study of the brightest OSARGs. This technical difficulty will be overturned when the ASAS I-band data (Pojmański, 2002) become available, which have a limiting I-mag of about 14, and a considerable overlap with OGLE-II photometry. Three years of Galactic Bulge data have been acquired by ASAS already. For our analysis we omit fields located far from $`b4^{}`$, and those fields for which a reliable extinction map is not available. This leaves us with a total of 33 OGLE-II fields.

## 3 Distribution of variables

The OSARG type A variables form a nearly horizontal sequence in the colour – magnitude diagram shown in Fig. 3. We fitted their distribution with a straight line using the least square method, and looking for the best global solution of the form:

$$I_{k,i,0}=I_{k,0}+a\left[(VI)_{k,i,0}2.881\right],$$
$`(2)`$
where $`I_{k,i,0}`$ and $`(VI)_{k,i,0}`$ are the reddening corrected average magnitude and colour for variable number $`i`$ in the field $`k`$, $`I_{k,0}`$ is a constant, different for every field, representing the ‘average’ magnitude of OSARG type A variables in field $`k`$, $`a`$ is the same constant for all fields, and $`2.881`$ is the mean $`(VI)_0`$ colour for all OSARG type A variables in all 33 fields represented in Fig. 3. We computed the value of the parameter $`a`$ as 0.169. The values of the parameter $`I_{k,0}`$ and the standard deviations from the relation (2) are given in Table 1 for all 33 fields. This table also lists the total number of type A and type B variables in each field and the Galactic coordinates of the field centers.

Fig. 4 presents the variation of the average magnitude of OSARG type A variables with the Galactic longitude. A clear trend is present, practically identical to that shown by the red clump giants (Stanek et al. 1994, Stanek et al. 1997, Sumi 2003), demonstrating that OSARG variables are located in the Galactic Bar.

We repeated the analysis for OSARG type B variables and it was qualitatively similar to that shown in Fig. 4, but disturbed by the OGLE-II saturation limit.

## 4 Multiple periods

Multiple periods were detected for most OSARGs. The ratios of those periods to the dominant period (i.e., the one with the largest amplitude) are shown in Fig. 5 as a function of the dominant period. Several bands corresponding to the period ratios close to 2, 1.4, 0.9, 0.75, and 0.5 are clearly seen. The interpretation of this rich spectrum of oscillations is beyond the scope of our paper. Our intent is to present this richness with a hope that it may attract specialists to investigate the nature of the oscillations.

Virtually all OSARGs were identified with objects in the 2MASS catalogue, and the values of J, H, and K magnitude are included in the tables available on Internet/Simbad. A sample of these tables is shown in our Table 2. These magnitudes were corrected for interstellar extinction using maps developed by Sumi (2003) and using extinction ratios between the infrared bands and V-band provided by Schlegel et al. (1998). Note that the extinction in K-band is small, and any errors in the estimate are of relatively little consequence. The average K-band magnitude, corrected for extinction, varies with Galactic longitude in a way similar to that seen in Fig. 4.

Next, extinction corrected $`K_0`$ magnitudes were also corrected for the Galactic Bar effect (as shown in Fig. 4 of this paper) to the distance corresponding to Galactic longitude $`l=0`$; we label the corresponding magnitudes $`K_{0,0}`$. A relation between the dominant period of oscillations and $`K_{0,0}`$ is shown in Fig. 6. The two groups of OSARGs that are visible in Fig. 1 are also clearly visible in Fig. 6.

The two groups of OSARGs are also plotted in the infrared color – magnitude diagram in Fig. 7. Only 20% of all variables are plotted to avoid crowding. There is almost full overlap between the two groups, distinctly different from the clear separation in Fig. 1 and Fig. 6.

Fig. 8 shows the dependence of the dominant period on the reddening corrected $`(VI)_0`$ color, and clear period – color relations are apparent. Note: K-band magnitude is close to the bolometric magnitude, and hence the period – luminosity relations are clearly visible in Fig. 6. However, the bright giants form a nearly horizontal line in the $`I_0(VI)_0`$ diagram shown in Fig. 2. Hence, the $`(VI)_0`$ color is for them a better absolute luminosity indicator than the $`I_0`$ magnitude, and this is clearly demonstrated in Fig. 8.

Finally, it is time for us to present a random sample of type A and type B OSARGs in Fig. 9. These are 12 variables of each type in the field $`\mathrm{BUL}\mathrm{\_}\mathrm{SC3}`$, sorted with increasing dominant period. These appear as semi-regular variables because of their rich period structure. A much longer time baseline of photometric measurements will become available soon when the current OGLE-III data are processed and published (cf. Udalski et al. 2002a). With a longer time baseline, it will be much easier to decide to what extent the multiple periods we found are stable, and to what extent OSARGs are multiply-periodic rather than semi-regular variables.

## 5 The Catalogue

We provide a catalogue of 8970 OSARG type A variables and 6399 OSARG type B variables in electronic form at:

http://www.astro.princeton.edu/$``$jwray/OSARG/

in the files A.txt and B.txt, respectively. Separate files for each field are available at the same URL, and are called A1.txt, A2.txt, …, A49.txt and B1.txt, B2.txt, …, B49.txt. The catalogue is a subset of the general OGLE-II catalogue of variable stars in the Galactic Bulge (Woźniak et al., 2001).

A sample of the tables, listing data for the first 25 type A stars from the field $`\mathrm{BUL}\mathrm{\_}\mathrm{SC1}`$, is shown in Table 2, which is split into two parts because of its considerable width. The consecutive columns are as follows: the variable star number according to Woźniak et al. (2001); the star number according to OGLE maps (Udalski et al. 2002b); 2000 RA and DEC; X and Y pixel coordinates; average V-band magnitude; V-band extinction; average (V-I) color; (V-I) reddening; average I-band magnitude; I-band extinction; J, H, K-band magnitudes; the dominant period and its amplitude; the next three significant periods and amplitudes. In cases where fewer than three significant periods were found, the extra columns contain the value -99.99999. Similarly, for a few stars there was no match found in the 2MASS catalogue, so the J, H, K columns for those stars contain the value -9.999.

The coordinates and I and V-band magnitudes are given following OGLE maps (Udalski et al. 2002b), the reddening and extinction is given following Sumi (2003), and the J, H, K magnitudes are given following 2MASS (All-Sky Point Source Catalogue, Release 2003 March 25). We provide all these data in electronically accessible files to facilitate analysis of OSARG variables to all interested users.

## 6 Discussion

There can be no doubt that OSARGs have been found in the MACHO database for the Galactic Bulge (Cook et al., 1997), and that they are related to red giant variables found in the LMC and SMC by Cook et al. (1997), by Kiss & Bedding (2003) and by Ita et al. (2003). However, there is also a striking difference between the Galactic Bulge and the LMC and SMC: the multiple narrow bands seen in the $`\mathrm{log}PK`$ diagrams for the LMC and SMC are not present in the Galactic Bulge, where we see only two groups of stars. It appears that our type A and B variables correspond to type $`\mathrm{A}^{}`$ and $`\mathrm{B}^{}`$ variables, respectively, of Ita et al. (2003, fig. 6). Note: the distance modulus to the Galactic Bulge is 4 mag smaller than that to the LMC, so the red giant tip, which is at $`\mathrm{K}12`$ in the LMC, is expected to be at $`\mathrm{K}8`$ in the Galactic Bulge. Unfortunately, this is close to the I-band saturation limit of OGLE, and this complicates the comparison and biases our diagrams.

We have an impression that our two sequences of variables are broader in the apparent magnitude than the corresponding sequences seen in the LMC and SMC. This may be caused by at least two effects: the Galactic Bar has a relatively large radial depth, which remains an effect even after correction is made for the inclination (cf. Fig. 4), and there is likely also a broad range of metallicities in our Galaxy as compared to LMC and SMC.

Another striking difference between OSARGs and similar variables in the LMC is the presence of a large number of LMC stars with color $`(\mathrm{J}\mathrm{K})>1.4`$, as shown in fig. 3 of Kiss & Bedding (2003); these are carbon stars. We found very few such variables, as is apparent in our Fig. 7, where there are less than a handful of stars off the right hand limit of the figure. This difference may be due to the higher metallicity of the Galactic Bulge as compared with the LMC.

The period – amplitude relations as presented in our Fig. 1 have not been reported in print to the best of our knowledge. However, similar relations were found for the LMC and SMC variables by Soszyński (2003).

We restrict this paper to the announcement of the OSARG variables in the Galactic Bulge/Bar, to the catalogue of 15,369 objects, and to the demonstration that they are located in a bar inclined to the line of sight. We make no attempt at modelling the Bar, as considerably more complete data should be available as soon as ASAS photometry (cf. Pojmański 2002) becomes available, and the problems with the OGLE saturation limit are alleviated.

## Acknowledgments

This work was supported with the NSF grant AST-0204908, and NASA grant NAG5-12212. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation