Challenges in 2D Stellar Modeling

Lovekin, Catherine C.

doi:10.3389/fspas.2019.00077

REVIEW article

Front. Astron. Space Sci., 09 January 2020

Sec. Stellar and Solar Physics

Volume 6 - 2019 | https://doi.org/10.3389/fspas.2019.00077

This article is part of the Research Topic The Future of Asteroseismology View all 19 articles

Challenges in 2D Stellar Modeling

$\nCatherine C. Lovekin$ Catherine C. Lovekin^*

Physics Department, Mount Allison University, Sackville, NB, Canada

Traditionally, stellar structure and evolution have been modeled with a series of concentric spherical shells. This description allows the star to be modeled in 1 dimension, greatly simplifying the calculations. However, as our understanding of stars becomes more advanced, the effects of non-symmetric effects must be included, which necessitates 2 or even 3 dimensional simulations. In this work, I discuss how 2D stellar models can help understand stars, improving our models of their pulsation frequencies, and allow us to place better constraints on their internal convection.

1. Introduction

The standard technique for modeling stars takes a 1-dimensional approach, and assumes the stars are relatively close to spherical. However, a large number of stars, especially high mass stars, are known to be rapid rotators (Huang and Gies, 2006; Huang et al., 2010), which breaks the spherical symmetry. In addition, a large fraction of stars are members of binary or multiple star systems, which can introduce tidal distortion. This distortion in surface shape and resulting variation in effective temperature has now been measured interferometrically for a number of stars (van Belle et al., 2001, 2006; Domiciano de Souza et al., 2003; Aufdenberg et al., 2006; Peterson et al., 2006; Monnier et al., 2007; Zhao et al., 2009; Che et al., 2011). To reproduce the observations requires 2D models of stars.

The effects of distortion are also becoming important as recent space missions like Kepler, CoRoT, MOST, BRITE, and TESS allow us to measure pulsation frequencies with very high precision over long times. These more precise frequencies make it easier to measure the effects of rotation and other distorting effects on the pulsation frequencies of the stars. Advances in spectropolarimetry have also given us a wealth of information on the magnetic fields of stars other than the Sun (Grunhut et al., 2017), which has triggered a number of new studies on the effects of magnetic fields on stellar evolution (Keszthelyi et al., 2018). As a result of these advances, it is becoming ever more important to improve our stellar models to be able to account for these multi-dimensional effects.

In the most general case, a static stellar structure model must obey four basic equations: Poisson's equation:

\begin{array}{l} \nabla^{2} ϕ = 4 π G ρ & (1) \end{array}

where ϕ is the gravitational potential and ρ is the density; mass conservation:

\begin{array}{l} \nabla \cdot (ρ v) = 0 & (2) \end{array}

where v is the fluid velocity; the momentum equation:

\begin{array}{l} ρ v \cdot \nabla v = - \nabla P - ρ \nabla ϕ + F_{v} & (3) \end{array}

where P is the pressure and F_v is the viscous force; and an entropy (or energy) equation:

\begin{array}{l} ρ T v \cdot \nabla S = - \nabla F + ϵ_{*} & (4) \end{array}

where S is the entropy, T is the temperature, F is the heat flux, and ϵ_* is the energy generation rate. This set of equations, combined with an equation of state and prescriptions for F and F_v fully define the behavior of the gas in the stellar interior, and are applicable for any rotation law.

Rather than use the full equations of fluid dynamics above, another approach uses the Roche approximation to get the pressure gradient:

\begin{array}{l} \frac{1}{ρ} \nabla P = - \nabla ϕ + \frac{1}{2} Ω^{2} {(r sin θ)}^{2} & (5) \end{array}

where ϕ is the gravitational potential, unmodified by rotation. This gives a purely radial gravitational acceleration g = ∇ϕ. Very little is known about the internal rotation profiles of stars, but arguments are often in favor of either cylindrical or shellular rotation laws. For conservative rotation laws like these, where Ω depends only on the distance from the rotation axis (cylinderical) or center of the star (shellular), the centrifugal acceleration can be derived from an effective potential V. Indeed, if the rotation is not conservative, it is impossible to define an effective potential, and it is not clear how the total potential should be determined. For conservative rotation laws, the centrifugal acceleration, potential, and rotation rate are all functions of the distance from the rotation axis, ϖ = r sin θ and the total potential is given by Φ = ϕ + V. This potential still follows Poisson's equation, so ∇²Φ = ∇²ϕ + ∇²V. For conservative rotation laws, ∇²V = −2Ω², so the total potential is given by ∇²Φ = 4πGρ − 2Ω², and the equation of hydrostatic equilibrium becomes:

\begin{array}{l} \frac{1}{ρ} \nabla P = - \nabla Φ = g_{eff} . & (6) \end{array}

The other equations of stellar structure must be written in 2D in the presence of rotation or other distortion, but the equations themselves take the same form as in the 1D case.

In 1D models, the equations can easily be written in terms of mass instead of radius. In 2D models, this is not generally possible. In the case of shellular rotation (Zahn, 1992), the radius is associated with isobars, and since most quantities are constant on isobars, the problem becomes a pseudo-1D calculation, and the equations can be written in terms of the mass as for the non-rotating case. For a fully 2D calculation, the mass inside a spherical shell of radius r is difficult to define, and is less useful as a coordinate. In general, 2D calculations either use the radius r and colatitude θ as coordinates, or map these coordinates to spheroidal coordinates r(ζ, θ).

As in the 1D problem, 2D models of stellar structure must apply boundary conditions at the surface, but in 2D, it is more difficult to determine what the appropriate conditions are. In 1D, the mass is used as the coordinate, and the surface is given by the radius where the enclosed mass is equal to the total mass of the star. As discussed above, in 2D the radius and colatitude (or some mapping of these variables) is used, and it is difficult to define the mass inside a radius r since the density is a function of colatitude θ as well. Some other method must be used to determine the radius at the surface of the star. Generally it is assumed that the surface of the star is an equipotential (Deupree, 1990; Maeder and Meynet, 2012; Espinosa Lara and Rieutord, 2013), but this still leaves the problem of how the location of the equipotential is determined.

When stars are not spherically symmetric the surface becomes distorted, the temperature varies with colatitude (discussed below) and this introduces large-scale meridional currents. The best-known of these is the meridional circulation introduced by rotation, but tidal or magnetic distortions can also introduce flows (Tassoul, 1978). These currents are expected to be slow, with theoretical velocities on the order of 3 × 10⁻¹⁰ cm s⁻¹ (Sweet, 1950). Despite the slow velocities, they are a significant source of mixing within the star. For shellular rotation laws, the effective chemical diffusion depends on the square of the velocity of the meridional circulation, so the contribution to chemical mixing is expected to be small in most cases (Zahn, 1992). The angular momentum transport is governed by an advection diffusion equation which can be horizontally averaged in the case of shellular rotation to give:

\begin{array}{l} ρ \frac{d}{d t} (ρ r^{2} \bar{Ω}) = \frac{1}{5 r^{2}} \frac{\partial}{\partial r} (ρ r^{4} U (r) \frac{\partial \bar{Ω}}{\partial r}) + \frac{1}{r^{2}} \frac{\partial}{\partial r} (ρ D r^{4} \frac{\partial \bar{Ω}}{\partial r}) & (7) \end{array}

where $\bar{Ω}$ is the horizontally averaged angular rotation rate, U(r) is the radial velocity of the meridional circulation, and D is the total vertical diffusion due to all transport mechanisms (Zahn, 1992; Ekström et al., 2012). However, because the velocity of the meridional circulation is so slow, it is very difficult to model these flows self-consistently on an evolutionary time scale.

Indeed, one of the biggest challenges to 2D modeling is the variety of timescales involved. To model the meridional flows or the details of internal convection requires high resolution simulations with short time steps, but these requirements are not possible with the computational power currently available. Modeling of multidimensional stars has so far taken two approaches: either to model the detailed structure and ignore the longer evolution timescales; or to neglect many of the internal processes and model the evolutionary timescales. However, unless we can find good models for the internal processes, the evolution problem is not completely separable from the short timescale simulations. To fully understand the evolution of distorted stars, we will need to understand how both angular momentum and internal composition gradients are transported through convection and turbulence, how the distortion affects radiatively driven winds and subsequent angular momentum loss, and the influence of magnetic fields.

In this review, I discuss some of the recent progress in modeling the structure, evolution, and pulsation frequencies of stars in 2D in section 2. In section 3 I summarize the effects of rotation on stellar evolution, highlighting the importance of 2D modeling. In section 4 I briefly discuss 2D modeling of binary stars, and in section 5 I discuss the challenges that 2D modeling still faces. A brief summary is presented in section 6.

2. A Brief History of 2D Modeling

2.1. Stellar Modeling

The first attempts to model stellar structure in 2D were performed by Sweet and Roy (1953), using first order perturbations to Cowling models. They assume that the stars were uniformly rotating, so that the total effective potential is given by ∇²Φ = −4πGρ + 2Ω². They then introduce the dimensionless variable $λ = Ω^{2} / 2 π G ρ_{c} (0)$ , where ρ_c(0) is the central density of the non-rotating star. All of the structure equations can then be expanded as a power series in λ. This method only worked for slowly rotating stars, or the perturbations become too large for the method to apply.

More realistic models were developed by Roxburgh et al. (1965), which divided the star into two regions. The innermost region, which contained most of the mass, was treated as spherical, while the outer, low density region is allowed to deform. This works since the innermost region is relatively unaffected by rotation, while the gravitational potential in the low density region depends primarily on the potential due to the inner region, which is approximately a Roche potential. Once the potential is known, the other structure equations can be integrated as for 1D stars. This method was found to work best for uniformly rotating stars.

Two dimensional static models have also been produced using the self-consistent field (SCF) method (Ostriker and Mark, 1968), later extended by Jackson (1970). This method uses an iterative approach to solve the equations of hydrostatic equilibrium and the gravitational potential to determine the stellar structure. Early implementations encountered convergence problems for stars below 9 M_⊙, limiting the applicability. Later implementations reformulated the method, and were able to avoid this problem (Jackson et al., 2005; MacGregor et al., 2007).

Eriguchi (1978) calculated the structure of rotating polytropes by transforming the equations of stellar structure to a complex plane. This work was expanded to include relativistic polytropes (Eriguchi, 1980), non-axisymmetric polytropes (Hachisu et al., 1982) and binary systems (Hachisu and Eriguchi, 1984a,b). However, this method still faced many limitations. In particular, the method breaks down if there are any discontinuities in the physical properties inside the star. Such discontinuities can be caused by a change in the composition or ionization state, which means this method cannot be used to study evolution.

However, Eriguchi and Sugimoto (1981) came up with another method, this time mapping the equation of hydrostatic equilibrium and Poisson's equation onto a 2D grid such that r_ij = r_i(θ_j). This allows the radial grid to adjust to the shape of the star at each angle. Using this method, it is possible to calculate the structure for a number of rapidly rotating, axisymmetric polytropes (Eriguchi and Sugimoto, 1981; Eriguchi and Mueller, 1985). The method can also be extended to piece-wise polytropes, which can then be used to approximate the structure of real stars.

A different approach to solving the equations of stellar structure used a 2D finite difference technique (Clement, 1974, 1978, 1979, 1994) on a 2D grid in colatitude θ and total potential. Rather than radius, the equations are rewritten to relate the mass and volume of the star, with the constraint that in the outermost zone, the enclosed mass must equal the total mass of the star. This technique was still limited to relatively slowly rotating stars, and is restricted to equilibrium models.

In recent years, 2D hydrostatic models have been calculated for use in asteroseismology (Roxburgh, 2004, 2006). These models use spectral methods to model the angular dependence of the physical quantities, and depend on knowing the density along one particular angle as a reference point. In this way, the density profile produced by 1D models that assumes shellular rotation can be expanded to provide a 2D hydrostatic model as input to stellar pulsation codes. This technique is important since several 2D codes have been developed in recent years to calculate stellar pulsation frequencies (Lignières et al., 2006; Reese et al., 2006; Ouazzani et al., 2012).

The first 2D evolution calculations with no theoretical rotation limit were performed by Deupree using ROTORC (Deupree, 1990, 1995, 1998). This code uses the fractional radius and the colatitude θ as the independent variables. One colatitude (typically the equator) is chosen to be a reference colatitude, and the surface location there has a fractional radius of 1. This defines the potential at the surface, and at each angle, the surface is defined to be the radial zone with the total potential closest to the value of the reference potential. Once the structure of the star is determined, the radius at the equator is determined by integrating over the density and requiring the total to equal the total mass of the star. In this way, the surface deformation caused by rotation can be modeled without including assumptions like von Zeipel's law (von Zeipel, 1924).

The ESTER code (Espinosa Lara and Rieutord, 2013; Rieutord et al., 2016) has been developed to model the internal processes of a rapidly rotating star in two dimensions. They use a modified coordinate system which can adapt to the shape of the star. These spheroidal coordinates (ζ, θ) are scaled such that ζ = 1 at the stellar surface, and then divide the star into several discrete radial domains. The transformation from spherical to spheroidal coordinates is given by:

\begin{array}{l} r & = & a_{i} ζ + A_{i} (ζ) [R_{i + 1} (θ) - a_{i} η_{i + 1}] + B_{i} (ζ) [R_{i} (θ) - a_{i} η_{i}] \\ η_{i} \leq ζ \leq η_{i + 1} & (8) \end{array}

where B_i(ζ) = 1 − A_i(ζ), and R_i+1 and R_i are the inner and outer boundaries of the domain. The functions A_i(ζ) and a_i are chosen to satisfy properties at the boundaries between the domains. Each of the radial domains is discretized into a spectral grid, so the horizontal dependence is expanded in terms of spherical harmonics. However, this code focuses on the shorter time-scale processes like meridonial circulation, rather than attempting to model the long-term evolution of a star.

A new code under development, 2DStars (Halabi et al., 2017) is intended to solve the stellar structure equations in 2D, assuming an axisymmetric geometry. Development of this code is still in the very preliminary stages, but if successful, it will be able to model the evolution of stars in a variety of axisymmetric environments, including the effects of binarity, disks, and rotation.

2.2. Pulsation Modeling

Fewer attempts have been made at modeling pulsation frequencies in 2 dimensions and these attempts have focused on the effects of rotation. Normally, eigenfrequencies of rotating stars are calculated using perturbation theory (Dziembowski and Goode, 1992; Soufi et al., 1998), but this is generally limited to slow rotation rates (see, for example Soufi et al., 1998, and references therein). Exactly how slow rotation is defined depends on the application, but it has been found that rotation rates as small as 50–75 km/s require more realistic calculations (Reese et al., 2006). Few people would argue that perturbation theory remains valid above rotation rates of 200–300 km/s.

The earliest attempts at 2D calculations were done by Clement (1981, 1989, 1998), and involved solving the linear pulsation equations on a 2D grid. This program, Non-Radial Oscillations (NRO) uses the degree ℓ to determine the parity of the mode, and odd modes are scaled by a factor of cos θ to make them symmetric. This allows the same method of solution to be used for both even and odd modes. The pressure perturbation is then described as expansion of k spherical harmonics, where k is the number of angular points in the grid. All other components of the eigenfunction can be expressed the same way, and the final solution is scaled to some fixed value at the surface, typically δr/r = 1 at the equator. While the initial attempts were limited to slow rotation, NRO (Clement, 1998) was able to solve the pulsation equations for very rapidly rotating stars. As initially developed, this method was restricted to uniform rotation but was later extended to include the effects of differential rotation (Lovekin et al., 2009).

Another approach uses spectral methods to calculate the eigenfunctions in 2D, including the effects of the Coriolis force and centrifugal distortion (Lignières et al., 2006; Reese et al., 2006). Using acoustic ray tracing of the frequencies calculated using this method, the resulting modes can be classified into chaotic modes, island modes, and whispering gallery modes based on their behavior in phase space (Lignières and Georgeot, 2009). Initially, this code was only applicable to polytropic models, but has since been expanded to more realistic stellar models (Reese et al., 2009). ACOR (Adiabatic Code of Oscillation including Rotation) (Ouazzani et al., 2012) has been developed using similar techniques, but has been designed to work with realistic stellar models.

For post-main sequence models, a number of attempts have been made to use full hydrodynamics to model the convection-pulsation interactions. This interaction is particularly important for classical pulsators, as the red edge of the instability strip is determined by the interaction between convection and pulsation (Tuggle and Iben, 1973). Early attempts focused on the interaction in one dimension, as was done in DYNSTAR (Ostlie, 1990). DYNSTAR has been used to model outbursts in Luminous Blue Variables (Guzik and Lovekin, 2014; Lovekin and Guzik, 2014). However, convection is inherently a multi-dimensional process, and the applicability of 1D models is limited. For this reason, a number of groups have developed 2 and 3D codes to model convective processes in pulsating stars. Cepheids have been modeled in 2D using the ANTARES (A Numerical Tool for Astrophysical RESearch) code (Mundprecht et al., 2013, 2015), focusing on the energetics of the convective zone and the possibility of developing descriptions that can be implemented in 1D models. Similarly, RR Lyrae have been modeled using SPHERLS (Stellar Pulsation with a Horizontal Eulerian Radial Lagrangian Scheme) (Geroux and Deupree, 2011, 2013, 2014, 2015). SPHERLS is able to model stars in either 2 or 3 dimensions, and has been used to calculate full pulsation amplitude models. The results agree well with observed RR Lyrae light curves. Although both ANTARES and SPHERLS are promising, calculations are time consuming, and neither has been widely used for modeling stellar pulsation.

3. Effects of Rotation

3.1. Effects on Evolution

The most pronounced effect of rapid rotation is the deformation of the stellar surface. The centrifugal force causes the star to become oblate in shape, and both temperature and emitted flux vary from pole to equator as a result. The first calculations of this effect found that the effective temperature is related to the effective gravity by:

\begin{array}{l} T_{eff} \propto g_{eff}^{β} & (9) \end{array}

where β is theoretically predicted to be 0.25 for purely radiative stars (von Zeipel, 1924). As a result of the rapid rotation, this gravity darkening will cause the equatorial regions to be cooler, and hence appear dimmer than the polar regions of the star. This effect has been directly measured using interferometry in several stars, including Achernar, Alderamin, Altair, Vega, Regulus, Rasselhague, α Cephei, α Leonis, α Ophiuchi, and β Cassiopeia (van Belle et al., 2001, 2006; Domiciano de Souza et al., 2003; Aufdenberg et al., 2006; Peterson et al., 2006; Monnier et al., 2007; Zhao et al., 2009; Che et al., 2011). All of these studies have found that gravity darkening does occur, but the observed temperature variation is inconsistent with β = 0.25. For example, a best fit value of β = 0.22 was found for α Cephei (Zhao et al., 2009). Values derived for other stars are similarly lower than the traditional β = 0.25 value given in von Zeipel's law (von Zeipel, 1924). One explanation for this difference is that these stars have thin surface convection zones, where the value of β is expected to be much lower, β ≈ 0.08 (Lucy, 1967).

The original gravity darkening model required the assumption of barotropicity, and is actually incompatible with solid body rotation in a radiative zone. To address these problems, a more realistic model for gravity darkening was derived (Espinosa Lara and Rieutord, 2011). This model relaxes von Zeipel's assumption of barotopicity, and is thought to be more applicable to real stars, particularly in cases of extremely rapid rotation. This new model was compared to calculations using 2D models constructed with ESTER (Espinosa Lara and Rieutord, 2011) which calculates the gravity darkening naturally as a result of the stellar structure equations, and the results were in good agreement with the new theoretical model. These results were later extended to include differential rotation (Zorec et al., 2017).

As a consequence of the flux variation, the flux seen by an observer varies with the inclination angle i, between the rotation axis and the line of sight to the observer. Inclination is very difficult to measure directly. Spectroscopic measurements give the Doppler broadening of spectral lines, which gives the component of the rotational velocity along the line of sight, v sin i. It is possible to decouple v and i (Stoeckley, 1968; Stoeckley and Buscombe, 1987; Reiners and Schmitt, 2003), but these methods depend on high-precision spectroscopy and very accurate synthetic line profiles. However, these techniques do not work well for slowly rotating stars (v_eq ≤ 200 km/s). As a result of the large quantity of high resolution photometric data from space missions including Kepler, CoRoT, and TESS, techniques have been developed to measure stellar rotation periods using asteroseismology, as discussed in more detail below.

The variation in observed flux as a function of inclination makes it difficult to assign an exact effective temperature and luminosity to a rapidly rotating star. The observed magnitude is a function of inclination, and the star is observed as a single point on a curve in the HR diagram, parameterized by the inclination. This effect is well known (Collins, 1963, 1966; Hardorp and Strittmatter, 1968; Maeder and Peytremann, 1970), but has recently been improved using more realistic stellar structure and atmosphere models (Lovekin et al., 2006; Gillich et al., 2008).

Rotation increases mixing through shear-layer instabilities, likely driven by differential rotation (e.g., Mathis et al., 2004). This mixing is likely to be latitude dependent, an effect that is not modeled by 1D parameterizations. The extra mixing introduced by these instabilities results in changes in the evolutionary tracks. Close to the Zero Age Main Sequence (ZAMS), rotation makes a star appear cooler and less luminous than a non-rotating star of the same mass. On the later main sequence, the behavior depends on how angular momentum transport is treated in the star, with some models finding that the rotating star remains cooler and dimmer (Deupree, 1995) while others find the rotating models become brighter (Ekström et al., 2012). In stars with M > 2M_⊙, rotation has also been found to extend the main sequence lifetime by about 25%, which in turn increases the observed width of the main sequence. This effect is illustrated in Figure 1 for 2.0 M_⊙ 1D-models calculated with MESA. Models with moderate amounts of rotation show good agreement with the observed width of the main sequence (Ekström et al., 2012). The rotation provides a source of extra mixing inside the star, increasing the amount of hydrogen available, and also potentially changing the surface abundances of other elements. In extreme cases, models suggest that high levels of rotation can produce so much mixing that the star evolves in a chemically homogeneous fashion (Szécsi et al., 2015; Song et al., 2016).

FIGURE 1

Figure 1. Evolution tracks along the main sequence for 2.0 M_⊙ models calculated with MESA. Shown are non-rotating models (blue) and models rotating at 0.3 v_crit (orange) and 0.5 v_crit (green). Initially, rotating models appear cooler and dimmer than non-rotating counterparts. As the star evolves, the rotating models become brighter, and spend a longer time on the main sequence.

Although stars generally spin down as they evolve, the effects of rotation can persist into later phases of evolution (Anderson et al., 2014). This is particularly evident for Cepheids, which have main sequence progenitors that typically rotate at v/v_crit = 0.3 − 0.4 (Huang et al., 2010). During the post main sequence evolution, the effects of rotation give these stars a luminosity that corresponds to that of a non-rotating model of higher mass. As a result, rotating stellar models have been used to resolve the Cepheid mass discrepancy (Anderson et al., 2014, 2016). Models indicate that rotating Cepheids are also older than non-rotating Cepheids, as rotation increases the main sequence lifetime of a star. This affects the period-age relationship, which in turn has implications for the use of Cepheids in Galactic archaeology. Finally, rotation extends the blue loop, and so rotating Cepheids traverse the instability strip more slowly than their non-rotating counterparts (Anderson et al., 2014).

The effects of rotation can be detected observationally via surface abundances. Models predict that the surface abundances of certain isotopes of nitrogen (¹⁴N) and helium (⁴He) should be enhanced, while light elements like lithium and boron should be depleted (Heger and Langer, 2000). Rotationally-induced mixing brings the lighter elements down to layers that are hot enough for these elements to be destroyed. These layers are also hot enough for the first steps in the CNO cycle to occur, converting the ¹²C and ¹⁶O to be converted into ¹⁴N, which is then brought to the surface. The theoretical predictions are generally in good agreement with the observations. For example, in this model, stars that have enhanced nitrogen abundances should also show depleted boron abundances, as was found by Venn et al. (1996). This model was also tested by comparing the surface abundance of ¹⁴N to the projected rotational velocity (v sin i) of stars observed as part of the VLT-FLAMES survey (Hunter et al., 2008). These results (see Figure 2) show that the nitrogen enhancement does correspond well with rotation velocity for most stars. Two populations, one with low rotation velocity and high nitrogen enhancement, and one with high rotation velocity and low nitrogen enhancement are thought to be the result of binary evolution, or possibly the presence of magnetic fields (Hunter et al., 2008).

FIGURE 2

Figure 2. The observed nitrogen enhancement plotted as a function of v sin i. The background density contours are predicted enhancements from population synthesis models, with individual stars plotted as either circles or triangles. The evolution tracks are for 13 M_⊙ stars at various initial velocities. Taken from Brott et al. (2011). Reproduced with permission.

The nucleosynthetic yields of a rotating star are also expected to differ from that of a non-rotating star, increasing the amount of metals produced. This effect is especially pronounced for carbon and oxygen, which increase by a factor of 1.5–2.5 for models with M ≲ 60M_⊙ (Maeder and Meynet, 2012). As discussed above, rotation provides a source of extra mixing, which increases the size of the convective core, so a rapidly rotating model will have yields similar to a more massive star. Above 60M_⊙, the total mass lost is the same in both rotating and non-rotating models, and this effect disappears (Maeder and Meynet, 2012).

3.2. Effects on Mass Loss

Stars with initial masses above ~ 10M_⊙ lose significant amounts of mass through radiatively driven winds. However, mass loss rates are highly uncertain, even when rotation is not considered. Mass loss rates can be derived from spectroscopic observations of wind speeds, can as high as $1 0^{- 5} M_{⊙} / y r$ in hot OB stars. These derivations often assume the winds are smooth, and detailed observations and simulations now suggest that the mass loss may be clumpy, which lowers the mass loss rates significantly. Mass loss rates calculated using Monte Carlo simulations (Vink et al., 2001) have been found to be in good agreement with observations (Mokiem et al., 2006). However, there are many indications that radiatively driven winds are clumpy rather than smooth, and this can reduce mass loss rates by a factor of 2 or 3 (Maeder and Meynet, 2012).

Many massive stars are rapid rotators, and so rotation is expected to have a significant impact on the way mass and angular momentum are lost from a star. Theoretical mass loss rates are sensitive functions of luminosity and effective temperature (Castor et al., 1975; Vink et al., 2001; Kudritzki, 2002), and so the variation in surface effective temperature and flux discussed above can be expected to cause variations in mass loss rates as a function of colatitude. Investigations into rotating winds have shown that there are two solutions: a classic Castor, Abbott & Klein (CAK) model (Castor et al., 1975) or fast wind solution (FWS) and a slow wind solution (SWS) which exist when the rotation speed is a significant fraction of the critical velocity (Curé et al., 2011). However, these models still assume that the mass loss from the star can be effectively described in 1D. Much work has been done on the 2D hydrodynamics of rotating stellar winds (e. g., Müller and Vink, 2014), but here I will discuss the effects on the star rather than the circumstellar environment.

The effect of gravity darkening on wind-driven mass loss has been investigated using stellar evolution codes in both 1D (Maeder, 2002; Georgy et al., 2011) and 2D (Lovekin, 2011). Initial investigations showed that rotating stars could be expected to produce a bipolar nebula, as the increased temperature near the poles enhances the mass loss rate in this region (Maeder, 2002). These models use a mass loss prescription that is based on the standard radiatively driven wind model of Castor et al. (1975) and Kudritzki (2002). An updated version of this model found an increase in the global mass loss rate as the rotation rate increased, and as expected, more mass was lost from the polar regions than the equatorial regions. This effect increased with the rotation rate (Georgy et al., 2011). Lovekin (2011) used fully 2D stellar structure models and mass loss rates calculated using Castor et al. (1975); Vink et al. (2001); Kudritzki (2002). In these 2D models, the surface is assumed to be an equipotential surface rather than relying on the Roche approximation, as is done in the Geneva models (Georgy et al., 2011). In the 2D models, the global mass loss rate actually decreased slightly as the rotation rate increased, as the increase in mass loss at the pole was offset by the decrease in mass loss from the equator. This ran counter to the conventional scaling relations for mass loss as a function of rotation rate (Bjorkman and Cassinelli, 1993). However these scaling relations are for stars at constant temperature and luminosity, while Lovekin (2011) used a fixed mass and allowed the temperature and luminosity to vary.

The differential mass loss can also affect the evolution of the angular momentum of the star. Georgy et al. (2011) and Lovekin (2011) found that the majority of the mass is lost from the polar regions of the star, however, the majority of the angular momentum is in the equatorial regions. The mass loss rates at the equator will have a significant effect on the angular momentum evolution. According to the 1D models, the majority of the angular momentum is lost from regions close to the equator, ranging from θ = 0 in the non-rotating models to θ ~ 70° at critical rotation (Georgy et al., 2011). As a result, they find that rotating models with anisotropic mass loss lose less angular momentum than those with uniform mass loss. They were able to evolve their models until the surface reached critical rotation, and found small differences in the evolution, with rotating models generally slightly cooler than their non-rotating counterparts. Lovekin (2011) also found that the angular momentum loss peaks near ~ 70° in the most rapidly rotating models, and that the rate of angular momentum loss increases overall as the rotation rate increases. In general, the 2D models lose less angular momentum than those with isotropic mass loss, as was found by Georgy et al. (2011).

3.3. Effects on Pulsation Frequencies

In non-rotating models, pulsation frequencies are described by two quantum numbers, n and ℓ. The first identifies the number of radial nodes in a particular eigenfunction, and ℓ determines the number of horizontal nodes in the latitudinal direction. The modes can be modeled as spherical harmonics, which have a third (azimuthal) quantum number, m, where m ≤ ℓ. For non-rotating models, the pulsation frequencies are independent of m for a given (ℓ, n) combination. In rotating stars, the degeneracy between modes of different m is lifted, and the modes are visible as multiplets in the pulsation spectrum. In the simplest case, the splitting is expected to follow:

\begin{array}{l} ν_{n ℓ m} = ν_{n ℓ} + m Δ ν & (10) \end{array}

where Δν is given by:

\begin{array}{l} Δ ν = \frac{1}{2 π} \int_{0}^{R} K_{n ℓ} (r) Ω (r) d r & (11) \end{array}

where Ω is the angular rotation velocity of the star and K_nℓ(r) is a rotation kernel, which depends on the real and imaginary components of the eigenfunction. Since the azimuthal number m can range from ℓ to −ℓ, the splitting should be uniformly spaced to first order. However, rapid rotation can complicate the spectrum by introducing non-uniform mode splitting (Reese et al., 2006). In extreme cases, the difference in frequency between members of a rotational multiplet can become larger than the difference in frequency between successive degrees ℓ, further complicating the interpretation of the spectrum. In cases where a rotational multiplet can be identified, the rotation frequency of the star can be determined, independent of the inclination. This type of analysis has been done for stars across the HR diagram (Charpinet et al., 2018; Paparó et al., 2018). If multiple triplets are detected, they can be expected to probe different radii inside the star, and it may be possible to probe the rotation profile (Aerts et al., 2003; Suárez et al., 2010).

Rotation is also known to decrease pulsation frequencies overall (Lignières et al., 2006; Reese et al., 2006; Lovekin and Deupree, 2008; Lovekin et al., 2009). As this difference depends on the order and harmonic degree (ℓ) of the mode, the rotation also changes the large and small frequency separation. However, one of the biggest challenges lies in identifying the correspondence between observed and theoretical modes for use in asteroseismology (Deupree, 2011).

A further challenge to interpretation of observed frequencies in rotating stars is the destruction of the regular frequency patterns expected for both p-mode and g-mode pulsators. For p modes, the frequencies are predicted to be regularly spaced, with a large separation that is approximately constant (Tassoul, 1980):

\begin{array}{l} Δ ν = ν_{n + 1, ℓ, m} - ν_{n, ℓ, m} . & (12) \end{array}

According to the ray-tracing 2D models developed by Lignières et al. (2006), the regular patterns seen in the p-modes of non-rotating stars disappear, and modes are instead divided into three categories: island modes, chaotic modes, and whispering gallery modes. Regular patterns do seem to exist for the island modes (Lignières et al., 2006), and these patterns have been successfully tested for more realistic models (Ouazzani et al., 2015).

Similarly, for g−modes, it is predicted that the periods will show regular spacings (Miglio et al., 2008). The g-mode patterns are less affected by rotation, but some new patterns do appear (Ballot et al., 2012). The regular period spacings in the g-modes do show the effects of rotation quite clearly, as shown in Figure 3. While the m = 0 modes are expected to have completely regular period spacings, the slope of the pattern is affected by the rotation rate. Retrograde modes (m = −1) show an upward slope, as shown in Figure 3, while prograde modes (m = +1) show a downward slope. In both cases, the effects are increased by rotation (Van Reeth et al., 2016; Lovekin and Guzik, 2017; Ouazzani et al., 2017). These patterns can be matched to observed g-mode spacings to determine stellar rotation rates, as has been done for many of the Kepler γ Doradus stars (Li et al., 2019).

FIGURE 3

Figure 3. Period difference as a function of period for ℓ = 1, m = −1 g−modes in a 1.6 M_⊙ model. As the rotation rate increases, the upward slope in the period spacing becomes more pronounced. Taken from Lovekin and Guzik (2017) ©AAS. Reproduced with permission.

4. Binary Systems

In close binary systems, the stars are distorted as a result of the gravitational pull of the companion, and this is not currently accounted for in most work on the evolution of binary stars. For example, de Vries et al. (2014) use 1D stellar models up to the point of Roche Lobe Overflow, and then convert to a 3D hydrodynamics simulation. Using a 2D stellar structure code, Deupree and Karakas (2005) have shown that the shape of both the surface and the convective core are distorted in a binary system, although the total mass of the convective core is essentially unchanged. For the closest binaries, they found that the shape of the surface differs slightly from the Roche potential usually used in 1D models. However, this assumes that the binary star system is coplanar and tidally locked. To fully simulate the distortion in a binary star system is a 3D problem.

Indeed, observations of heartbeat stars with the Kepler telescope have shown that even in close eccentric binaries, the rotation period can be different from the orbital period (Hambleton et al., 2018). In close systems like this, the gravitational field can even produce tidally excited oscillations, and pulsations have also recently been detected in η Carinae, which have been proposed to be tidally excited (Richardson et al., 2018). Tidally excited oscillations were first suggested in binary stars where the rotation period of the components is different from the orbital period of the system (Cowling, 1941), so these systems will necessarily be at least spheroidal, and are highly likely to be triaxial.

5. The Future: Challenges and Needs

Multi-dimensional models are common in simulations of short-timescale events inside stars. Several groups have modeled the internal convection of stars (Woodward et al., 2018; Arnett et al., 2019), but even these simulations often don't include rotation. Typically, convection simulations model a part of a star, either a wedge or a shell. In either case, these simulations are much shorter in timescale than evolution calculations, and the base structure is effectively static. Rotation also breaks the symmetry between the latitudinal an azimuthal (θ and ϕ) directions, so the simulations must be 3D to capture all the details (Viallet et al., 2013). The rotation also changes the way the composition of the star changes, introducing non-local effects (Arnett and Meakin, 2010). Although great progress is being made in these 3D convection simulations, the computational power required to model both short and long timescales in a single simulation does not exist. Much work is being done on determining algorithms that can be used to replace the standard Mixing Length Theory (MLT) (Böhm-Vitense, 1958), in 1D stellar models, but this will still only approximate some of the effects.

In most 1D stellar evolution codes, rotation is included using a shellular approximation, which essentially reduces the problem to 1D (Zahn, 1992). This type of rotation law is thought to arise from strong horizontal turbulence, which forces the rotation rate to be constant along isobars. Transport of material and angular momentum via the meridional circulation can be modeled as an effective diffusion in this model. Models suggest that shellular rotation is probably a good approximation at low rotation rates, while at higher rotation rates, models produce rotation profiles that do depend on latitude (Kitchatinov and Rüdiger, 1999; Espinosa Lara and Rieutord, 2013; Rieutord et al., 2016). Nevertheless, this is the approximation made by the vast majority of stellar evolution codes in use today. The exceptions are ESTER (Espinosa Lara and Rieutord, 2013; Rieutord et al., 2016), which fixes a surface velocity V_eq and allows the angular momentum to evolve, and ROTORC (Deupree, 1990, 1995, 1998), which accepts an initial rotation profile, and then evolves the star with local conservation of angular momentum.

It is reasonable to expect that stars are uniformly rotating, at least initially, but as they evolve, angular momentum will be transported inside the star. Models that are created with a uniform rotation profile on the ZAMS will quickly evolve to a differential rotation profile with the core rotating more rapidly than the surface (Denissenkov et al., 1999). This profile is in good agreement with the few stars that have had core to surface rotation rates measured using asteroseismology (Aerts et al., 2003; Lovekin and Goupil, 2010). Internal transport mechanisms like meridional circulation tend to bring angular momentum from the core to the surface, where it can be carried away with the stellar winds (Ekström et al., 2012). Since the mass loss is so much more efficient in more massive stars, these models spin down more rapidly, which in turn affects the size of the core, effective temperature, and luminosity of the star. In order to accurately model these properties, the details of internal transport must be correctly included.

In massive stars, many observed magnetic fields are approximately dipolar, the magnetic axis is frequently offset from the rotation axis of the star, making this a 3D problem (Grunhut et al., 2017). In lower mass stars, magnetic fields are produced by a dynamo, which is also a 3D problem, and needs to be parameterized to include the effects in 1 or 2D models. This remains challenging to implement, and at present can only be done for simplified systems (Jennings et al., 1990). Nevertheless, magnetic fields have been included in a few 1D codes (see, for example Feiden, 2016; Keszthelyi et al., 2018).

The effects of magnetic fields on stellar evolution can be significant. Magnetic fields couple with the rotation of a star, and so affect the internal transport and mass loss of a star. The effects on internal energy transport are particularly pronounced in the convective regions, and this has been suggested to be the origin of the radius inflation problem in M dwarfs (Jackson et al., 2018). One promising line of inquiry uses post-common-envelope binaries as probes of the magnetic field (Navarrete et al., 2019). In this case, it appears that the magnetic field can have a measurable effect on the stellar structure. Clearly, magnetic fields are an important part of stellar evolution and need to be accounted for in the models. This remains a challenging problem though, as magnetic fields are intrinsically 3D and act and evolve on timescales much shorter than typical stellar evolution.

Magnetic fields are also responsible for the roAp stars, where the oscillations arise from magneto-acoustic waves (Elkin et al., 2008) which will require both 2D structure and pulsation calculations to understand fully. In addition, a number of chemically peculiar stars show horizontal inhomogeneities in the surface abundance patterns (Khokhlova, 1976). Some work has been done on the effect of magnetic fields on the pulsation frequencies in roAp stars (Saio, 2014), but the 2D effects of magnetic fields on rotation have not been investigated in detail.

Of course, many groups are working on the 2D and 3D properties of stellar convection and turbulence (Davis et al., 2018; Arnett et al., 2019), but these simulations usually cover only a small fraction of the star, and cover a timescale that is short compared to the evolutionary timescale. The ESTER code described above (Espinosa Lara and Rieutord, 2013; Rieutord et al., 2016) models the stellar structure in 2D, and includes the effects of meridional circulation and rotation, but is not capable of modeling the longer timescales required to describe the evolution in 2D.

To date, there is no code capable of modeling the evolution of a star in 2D, while also taking into account the effects of short term processes like meridional circulation and convection. This is a serious gap in our abilities, and limits our ability to model stars with strong magnetic fields, in tidally distorted systems, or with rapid rotation. In order to progress, models for the short term processes that can be included in the evolutionary calculations need to be developed. Since 2D calculations will intrinsically take longer than their 1D counterparts, it would be appropriate to develop scaling laws that can be used in 1D, as has been done for pulsations in rotating stars (Castañeda and Deupree, 2016).

6. Conclusions

The equations of stellar structure in 1D have been well understood for nearly 100 years, but almost immediately, it was realized that the effects of rotation, binarity, and magnetic fields would introduce 2D or 3D effects that needed to be taken into account. The first attempts were made in the 1950s (Sweet and Roy, 1953), yet 70 years later, we still do not have a way of realistically modeling the evolution of a star in 2D. This remains a challenging problem, and trade-offs must be made between the ability to model detailed internal processes and the long-term evolution of the star. This multi-dimensional structure becomes ever more important as improved observations give us more details about the surface properties, binary interactions, magnetic field strengths, and pulsation frequencies. In order to accurately reproduce these observations, our models must be able to take into account the multi-D effects in a self-consistent fashion.

From a variety of box-in-star simulations, as well as structure simulations with ESTER (Espinosa Lara and Rieutord, 2013; Rieutord et al., 2016) we know that convection and rotation produce 2D effects, and that these will in turn change the evolution of the rotating stars when compared to non-rotating stars. From the nucleosynthesis in the core to the mass loss from the surface, these multidimensional processes influence the evolution of the star. A good understanding the process of stellar nucleosynthesis and the nature of supernova progenitors are the building blocks for Galactic archaeology, as well as the history of the universe. Previous studies have shown that the details of stellar evolution can influence stellar population synthesis and the derived masses of galaxies, affecting the way we understand the history of our universe (Conroy, 2013).

The evolution of a star in 2D is a very challenging problem, and requires improvements in both theoretical underpinnings and computational resources. At the same time, the disparity in size and time scales required to calculate the internal processes and the evolution suggests that no improvement in computational resources is going to make it possible to include all aspects of the multi-D structure in a single calculation. Further complicating the problem, convection, rotation, and magnetic fields are all coupled, and will interact in ways that are likely non-linear. Rather than hoping to create an ideal multi-D calculation, we must focus on finding better models for these internal processes that can be incorporated into models that include their 2D effects on evolutionary timescales. In particular, a better model for convection is needed to replace the current MLT theory, as well as a better understanding of angular momentum transport inside stars. This is an active area of research, and several groups are working on an updated convection model to replace MLT in 1D models (Garaud et al., 2010; Canuto, 2011; Arnett et al., 2019). Magnetic fields have been included in 1D models, but the details of their 2D influence on chemical mixing and angular momentum transport cannot be ignored. Many groups are working on large-scale simulations of these processes, with a goal of better understanding how they work, and how they can be modeled in evolution simulations. We are starting to progress in all of these areas, and computational resources are becoming faster, and more widely available (more cores), which will make it feasible to run 2D or 3D evolution simulations in the near future.

Author Contributions

The author confirms being the sole contributor of this work and has approved it for publication.

Funding

This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1125915. This work was supported in part by NASA ATP grant 12-ATP12-0130 and in part by the Natural Sciences and Engineering Research Council of Canada (NSERC).

Conflict of Interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The author gratefully acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC). This research has made use of NASA's Astrophysics Data System.

References

Aerts, C., Thoul, A., Daszyńska, J., Scuflaire, R., Waelkens, C., Dupret, M. A., et al. (2003). Asteroseismology of HD 129929: core overshooting and nonrigid rotation. Science 300, 1926–1928. doi: 10.1126/science.1084993

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, R. I., Ekström, S., Georgy, C., Meynet, G., Mowlavi, N., and Eyer, L. (2014). On the effect of rotation on populations of classical Cepheids. I. Predictions at solar metallicity. Astron. Astrophys. 564:A100. doi: 10.1051/0004-6361/201322988