A Possible Mechanism for Driving Oscillations in Hot Giant Planets

To study the stability of pulsations, we developed a code that computes the work function in Equation (1) from the profiles of 1D stellar and planetary models and their associated perturbations. For the 1D models, we employ the Modules for Experiments in Stellar Astrophysics (MESAs) stellar evolution code (Paxton et al. 2011, 2013, 2015). For our purposes, MESA uses the calculations of Ferguson et al. (2005) to implement opacities in the low-temperature regimes necessary for gas giants. These models are then used as input for the GYRE suite of stellar oscillation codes (Townsend & Teitler 2013), which utilizes a multiple shooting scheme to solve the nonadiabatic oscillation differential equations to produce a spectrum of eigenfrequencies, eigenfunctions, and perturbed quantities.

These tools were first tested by replicating the differential work functions of models for two classes of well-studied stars that exhibit pulsations via the κ-mechanism: β Cepheid stars and δ Scuti stars. More specifically, we calibrate by trying to match the work function in Figure 2 of Pamyatnykh (1999), particularly the l = 1, p1 (n = 1) and p4 (n = 4) modes for the β Cep star and the l = 1, p1, and p7 (n = 7) modes for the δ Sct star. l is the angular degree and n is the radial order of the acoustic oscillations.

Each stellar model is created using the parameters listed in Pamyatnykh (1999), and the frequencies of the modes calculated with MESA and GYRE are within a few percent of those published values (β Cep: 61.55 μHz & 106.18 μHz for p1 and p4, respectively; δ Scu: 170.66 μHz & 384.28 μHz for p1 and p7, respectively). The plots of the differential work function given by Equation (1) are shown in Figure 1 for each star and are nearly identical to those found in Pamyatnykh (1999).

Zoom In Zoom Out Reset image size

Oscillations are excited (instability) when the work function is greater than zero. Thus, in Figure 1, if the area under the curve of the differential work is positive, that particular mode is excited, otherwise it is damped. Positive peaks correspond to driving regions, while negative peaks are damping regions. Therefore, in the β Cep star, the p1 (p4) mode is unstable (stable), and likewise, in the δ Sct star, the p1 (p7) mode is unstable (stable), which matches the results in Pamyatnykh (1999). It must be noted that a positive work function does not necessarily mean that a mode is driven. It may be excited, but it must also have a positive growth rate or else it will be quickly damped out. A more detailed discussion of mode growth rates is given in Section 4.3. Finally, the driving regions must also be near the surface, as perturbations deep in the stellar interior must be exceedingly large to overcome the thermal inertia of the surrounding medium.

To study the behavior of any nonadiabatic oscillations excited by the κ-mechanism within Jupiter, we use MESA following this prescription:

• 1.  
  A gas-giant planet is created with a radius of 1.2–2 Jupiter radii with mass fractions of X = 0.74, Y = 0.24, and Z = 0.02. It has a mass of 1 Jupiter mass minus the mass of a currently non-existent core. The model is then evolved for 1000 years to allow for relaxation of the internal equation of state.
• 2.  
  A core with an average density of $10\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and the specified mass needed to bring the planet total to one Jupiter mass, is slowly added over 2000 years to allow the equation of state to relax once more.
• 3.  
  The model is evolved over 4.5 Gyr, allowing the interplay between gravity and internal pressure to evolve the equation of state. During this time, the surface of the planet may be irradiated with some specified bolometric solar flux, penetrating to a specified maximum column depth. At the end of the evolution, the radius of the planet has slightly decreased due to gravitational collapse.
• 4.  
  Each model is then fed into GYRE, where it is linearly scanned for oscillations between 0.1 and 5 mHz for the $l=0,1,2$ angular degrees and radial orders $n=1\mbox{--}10$. The output from both GYRE and MESA includes all the necessary components of Equation (1). It can then be determined which modes, if any, are potentially excited within Jupiter.

Table 1 lists six particular Jupiter models that are explored. It lists whether the planet is irradiated with $\sim {10}^{9}$ erg cm⁻² s⁻¹ (High) or $5.03\times {10}^{4}$ erg cm⁻² s⁻¹ (Solar) of irradiation, its core mass, and its final radius after 4.5 Gyr. Each model is irradiated to a column depth of 300 cm² g⁻¹ ($\approx 0.9967$ R_J). How far the flux penetrates into the Jovian atmosphere is unknown, hence the column depth of irradiation is somewhat arbitrary. We ran an additional 810 models to investigate the dependence of mode excitation on the column depth of irradiation. By varying the column depth throughout these models, we found that even for various flux values (greater than 10⁹ erg cm⁻² s⁻¹), oscillations can be excited as long as the column depth exceeds 50 cm⁻² g⁻¹. For the set of models in Table 1, 300 cm² g⁻¹ corresponds to a pressure depth of ∼1.5 bars (the surface is defined at ∼0.3 bars). Gierasch et al. (2000) cite observations to a depth of ∼3 bars of pressure for infrared observations. If IR light can be observed radiating from the 3 bar level we assume it can penetrate to the 3 bar level (or a column depth greater than 600 cm² g⁻¹). Since shorter wavelength light will be more readily absorbed, we split the difference between our lower and upper limits and select 300 cm² g⁻¹. Therefore, our estimate of ∼1.5 bars may be somewhat conservative. Finally, it is important to note that these are one-dimensional models and thus assume spherical symmetry. Therefore, we must also assume these are rapidly rotating (i.e., on the order of the rotation periods of Jupiter and Saturn) planets such that any consequences of stellar radiation only penetrating half of the planet's surface at a given time can be ignored. However, this may not be a large concern because a tidally locked planet with the driving region always facing the host star could still experience global oscillations just as how Jupiter oscillates due to (probable) excitations of localized moist convective storms.

Table 1.  Parameters of the Jovian Models

Model Name	Irradiation	Core Mass (${M}_{\oplus }$)	Radius (RJ)
A1	High	0	1.1045
A2	Solar	0	1.1045
B1	High	10	0.9691
B2	Solar	10	0.9691
C1	High	15	0.9532
C2	Solar	15	0.9532

Download table as:  ASCIITypeset image

The high amount of irradiation applied to models A1, B1, and C1 is about five orders of magnitude larger than the solar flux at Jupiter. We had initially made a rather trivial mistake implementing an incorrect amount of flux in the models and subsequently found surprising results (explained in Section 4). We thus decided to keep these three models for the analysis (even though they do not represent Jupiter proper), and furthermore, were motivated to consider the more realistic cases of hot Jupiters orbiting very close to their host star.

We generated a set of exoplanet hot-Jupiter models that are created using the same method as described above, with a few differences. Table 2 lists the subset of main-sequence stellar models selected as host stars for these planets, their masses, their effective temperatures, their bolometric luminosities and their ages. The stellar ages are determined by running MESA main-sequence models until 10% of their core hydrogen abundance is depleted, or 3 Gyr, whichever comes first. The planetary models are then evolved to the same age as the corresponding star.

Table 2.  Parameters of the Stellar Models That Are Host Stars to Hot Jupiters

Spectral Type	Mass (${M}_{\odot }$)	$\mathrm{log}{T}_{\mathrm{eff}}$ (K)	Luminosity (${L}_{\odot }$)	Age (Gyr)
M5	0.28	3.565	0.0126	3
M0	0.49	3.584	0.0383	3
K5	0.58	3.598	0.0645	3
K0	0.76	3.670	0.216	3
G5	0.91	3.731	0.523	2.258
G0	1.05	3.767	1.027	1.264
F5	1.33	3.816	3.144	0.535
F0	1.59	3.871	6.828	0.370
A5	1.98	3.944	16.700	0.203

Download table as:  ASCIITypeset image

Representative planetary mass and orbital semimajor axis ranges are taken from NASA's exoplanet archive.¹ We created a parameter space of masses ranging from 1 to 30 Jupiter masses and orbital semimajor axes of 0.01–2 au. The irradiation flux received by each planet is then the bolometric flux at the particular semimajor axis calculated from a stellar blackbody given by the host star's effective temperature. Each planet model is irradiated to a column depth of 300 cm² g⁻¹ over its entire lifetime. The models have a core mass of 10 Earth masses with an average core density of $10\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. We compute planetary models for 9 different host-star spectral types with 11 different planet masses at 12 different semimajor axes for a total of 1188 models. We find that varying the column depth penetration of radiation has negligible effect on the results presented in subsequent sections, as long as it is greater than 100 cm² g⁻¹.

We carry out a nonadiabatic oscillation analysis for the Jupiter models. As one might expect (e.g., Bercovici & Schubert 1987), there is not as large a source of internal energy in Jovian planets as there is in stars, and thus we do not find that the κ-mechanism (or any other mechanism considered here) excites modes within the nominal Jupiter models. All calculated modes in models A2, B2, and C2 are stable and are thus not excited. Frequencies range between 0.129 and 1.74 mHz depending on the l and n values. The differential work function of the models is virtually zero throughout Jupiter except for the outermost ∼1% of its radius. In the outermost layers, there is a strong negative feature, that is, a strong damping region, that contributes to energy lost over an oscillatory cycle. Therefore, these modes are not excited.

The introduction of the significant ($\sim {10}^{9}$ erg cm⁻² s⁻¹) solar irradiation throughout Jupiter's evolution, however, results in excited oscillations. In models A1, B1, and C1, several of the low radial order modes are found to be unstable. For example, model A1 (no core) revealed that modes l = 0, n = 1, and l = 2, n = 1 are excited (recall that angular degrees greater than 2 and radial orders greater than 10 are untested). The reason the oscillations are excited in the presence of extreme stellar radiation is due to a large luminosity perturbation that occurs at about $R=0.993\ {R}_{J}$, a direct consquence of the strong irradiation. A more detailed discussion as to what is happening is given in the following sections for the case of the hot-Jupiter models.

The highly irradiated Jupiter models that demonstrate excited modes motivated us to explore gas-giant planets orbiting very close to their host stars. The parameter space described in Section 3.3 reveals some interesting trends regarding which planets exhibit mode excitation. Figure 2 shows the models for which unstable modes are found in terms of orbital distance and host-star effective temperature. We observe a very well-defined region of excitation favoring short-period planets or hot host stars.

Zoom In Zoom Out Reset image size

Regardless of planetary mass, we find that approximately $\sim {10}^{9}$ erg cm⁻² s⁻¹ of flux are required for oscillations to occur. That flux or higher does not guarantee the existence of oscillations, but it is a requirement for them to be present. Thus, the separation line in Figure 2 is given by

$$\begin{eqnarray}&&d=\sqrt{\displaystyle \frac{\sigma {R}_{* }^{2}{T}_{\mathrm{eff}}^{4}}{F}},\end{eqnarray} \tag{ 2 }$$

where d is the orbital distance, σ is the Stefan–Boltzmann constant, R_* is the stellar radius, and the flux $F={10}^{9}$ erg cm⁻² s⁻¹. We see that this relation does a good job of discriminating between the planets that have excited modes and those that do not.

According to Equation (1), there are three factors that influence the sign of the work function (ignoring perturbations to the energy generation rate). Figure 3 shows the comparison of these three quantities as well as the differential work function between the A5, 30 M_J planets at orbits of 0.01 au (blue) and 0.2 au (red), for the l = 0, n = 2 mode. This particular mode is excited in the 0.01 au model, while in the 0.2 au planet it is not. We choose to display these models, in particular, as they have the most exaggerated profiles so the effects are easiest to see.

Zoom In Zoom Out Reset image size

Throughout the interior ${{\rm{\Gamma }}}_{3}-1$ is positive and similar in magnitude between the models and it is not a critical parameter. We see that for this particular solution the density perturbation for this mode is negative in both models (recall that we are solving the time independent equations). So we are at maximum expansion due to the oscillations. To obtain a postive work function, therefore, a large (positive) luminosity perturbation is thus needed, such that the negative gradient of this term is negative. Indeed this is what is found for the close-in planetary model, where it is observed how much further in radius the light penetrates the outer layers of the planet model at 0.01 au. An oscillation induces a positive luminosity perturbation over a larger fraction of the planet's radius.

To explore this situation in more detail, Figure 4 shows the interior density profiles of the $30{M}_{J}$ planets around the A5 star at increasing orbital distance. As the planet evolves at a further distance from the star, it does not experience significant heating, and therefore becomes a denser planet with a smaller radius. However, the closer the planet is to the star, the more puffed up and diffuse the outermost layers become. This allows light to penetrate deeper into the planet, thereby allowing for a larger region over which significant luminosity perturbations due to modes can occur. Some hot Jupiters have indeed been observed to have larger radii than initially predicted (Bakos et al. 2007; Hartman et al. 2012). The positive luminosity perturbation due to the oscillation heats the gas, causing a negative density perturbation by rarefying it. This newly expanded gas acts like a restoring force, but it is not large enough to win out over the induced rarefaction by the external heat source. This is a general trend seen over all stars and all mass ranges.

Zoom In Zoom Out Reset image size

In addition, for the model with the excited n = 2 radial mode, its frequency is about 0.3 mHz, or a period of an hour. The thermal timescale in the driving region is similar to this, and thus the thermal conditions for the mode are favorable. The driving takes place in a "transition region" (Pamyatnykh 1999). This is also the case for all of the modes we consider that are unstable.

Finally, we must also consider the convective timescale in relation to these oscillations. Typically, convection is expected to stabilize the planet against such pulsations. However, the timescale of convection, estimated in a similar fashion as in the case of the Sun (i.e., the dynamical timescale divided by a measure of superadiabaticity), is quite large. In the driving region for the 0.01 au, 30 Jupiter mass model around the A5 star, for example, the convective timescale is approximately a couple hundred hours, or 2 orders of magnitude larger than thermal timescale and the pulsation period. Therefore, we expect convection to be a passive process regarding these pulsations.

In addition to the work function needing to be positive in order to excite modes, as well as a thermal timescale similar to the mode period, we also examined the η growth rate parameter as discussed in Ando & Osaki (1975). This parameter is given by

$$\begin{eqnarray}&&\eta =\displaystyle \frac{-{\sigma }_{i}}{{\sigma }_{r}}\end{eqnarray} \tag{ 3 }$$

where ${\sigma }_{i}$ is the imaginary component of the (nonadiabatic) mode eigenfrequency and ${\sigma }_{r}$ is the real component. If $\eta \gt 0$ the mode is overstable and will grow, and if $\eta \lt 0$, the mode is understable and will be damped out. Throughout this analysis, if the work function of a particular mode is positive but η is negative, then the mode is considered to be excited but quickly damped. In Figure 2, the modes that are excited and quickly damped are listed as unexcited data points.

The values for η, for both the Jupiter and the hot-Jupiter modes, are mostly within the range calculated by Ando & Osaki (1975). Their values range from 10⁻¹⁰ to 10⁻⁴ whereas our values range from 10⁻¹⁴ to 10⁻³. However, the best indication as to whether these modes are overstable by examining η is actually to compute $\omega /\eta$, where ω is the frequency of the mode (in mHz). The mode is almost certainly overstable for large values of $\omega /\eta$. We calculate a range of 98.8–$2\,\times \,{10}^{13}$ mHz with the vast majority of the modes at the upper end of that range.

The process by which modes are excited via radiative suppression can now be explicitly explained. The combination of the positive (negative) luminosity perturbation and negative (positive) density perturbation is what is leading to unstable modes. The radiative luminosity perturbation can be shown to be (Pamyatnykh 1999)

$$\begin{eqnarray}&&\displaystyle \frac{\delta {L}_{r}}{{L}_{r}}=\displaystyle \frac{{dr}}{d\mathrm{ln}T}\displaystyle \frac{d}{{dr}}\left(\displaystyle \frac{\delta T}{T}\right)-\displaystyle \frac{\delta \kappa }{\kappa }+4\displaystyle \frac{\delta T}{T}+4\displaystyle \frac{\delta r}{r}.\end{eqnarray} \tag{ 4 }$$

In our case, the first term on the RHS of Equation (4) contributes to oscillatory damping and the later two terms are negligible. The dominant term on the RHS is the opacity perturbation, $\delta \kappa /\kappa$, that drives the oscillations. For example, when we do detect an unstable mode, we indeed find at maximum expansion that the opacity perturbation is rather large and negative, such that Equation (4) is positive. Figure 5 shows an example from the same extreme model used earlier.

Zoom In Zoom Out Reset image size

The mechanism at work is as follows. Some weak perturbation in the gas causes it to collapse at a given location. Upon maximum compression, the density, opacity, and temperature all increase. In normal conditions, the temperature increase allows the gas to cool via radiative losses more readily. However, this cooling is suppressed by the high amounts of external radiation heating the gas due to the increased opacity. Fundamentally, it is this absorption of energy that is converted to work to drive the mode. With the increase in energy from the photonic heating, the gas then begins to expand, pushing past the point of equilibrium and eventually to a state of maximum expansion. Here, the opacity perturbation is negative and the gas cools because the luminosity perturbation is positive, increasing outward (Figures 3 and 5). The recompression after maximum expansion overshoots the point of equilibrium and the process repeats.

This process is in contrast with the "classical" kappa mechanism in two fundamental ways. First, the kappa mechanism is intrinsically an internally driven phenomenon with the radiation originating from nuclear fusion whereas the radiative suppression is an externally driven phenomenon. Second, in normal stellar conditions, a compression results in a decrease in opacity rather than an increase. The reason the kappa mechanism works is because in unique regions within the stellar interior, compression raises the temperature of the gas to a level at which hydrogen or helium can be ionized, which thereby increases the opacity. The radiative suppression mechanism does not rely on such ionization as a prerequisite.