Exoplanet Modulation of Stellar Coronal Radio Emission

Our method employs a numerical model of the corona and wind of a star computed self-consistently with an orbiting planet stationary in the frame rotating with the stellar rotation. To produce a solution for the stellar corona, we use the BATS-R-US MHD model (Powell et al. 1999; Tóth et al. 2012) and its version for the solar corona and solar wind (van der Holst et al. 2014). The model is driven by photospheric magnetic field data, while taking into account the stellar radius, mass, and rotation period. The non-ideal MHD equations are solved on a spherical grid which is stretched in the radial direction, taking into account Alfvén wave coronal heating and wind acceleration that are manifested as additional momentum and heating terms. The model also takes into account thermodynamic heating and cooling effects, Poynting flux that enters the corona, and empirical turbulent length-scales. We refer the reader to Sokolov et al. (2013) and van der Holst et al. (2014) for a complete model description and validation.

The model provides a steady-state, self-consistent, three-dimensional MHD solution for the hot corona and accelerated stellar wind (in the reference frame rotating with the star), assuming the given, data-driven boundary conditions. Thus, it provides the three-dimensional distribution of the plasma density, temperature, velocity and magnetic field (the complete set of MHD plasma properties). With all the plasma parameters defined everywhere, we can also deduce the plasma frequency, ω_p, at each cell, which enables us to track the refraction of the radio waves through the model domain.

In this study, we use a simple, solar-like dipole field of 10 G, and solar values for the radius, mass, and rotation period. As the baseline for our investigation, we choose to study a Sun-like star, but we also study a case with a stellar dipole field of 100 G representing a moderately active star, or non-active M-dwarf star. As we demonstrate in Section 3, this simplified setting might be sufficient to qualitatively capture the main radio modulation effects regardless of the particular star we use.

The strength of the planetary field is chosen to be 0.3 G (Earth-like) and 1 G (Jupiter-like), with semimajor axis, a, of 6, 9, 12, and 15 R_⋆ located along the x = 0 axis. These distances translate to 0.028, 0.042, 0.056, and 0.070 au. The planet is embedded as a second boundary condition as described in Cohen et al. (2011a), where we use planetary boundary number density value of 10⁷ cm⁻³, and planetary boundary temperature value of 10⁴ K. These values produce a thermal outflow from the planet in the range of ${10}^{6}-{10}^{7}\,{\rm{g}}\,{{\rm{s}}}^{-1}$, which is much lower than observed in hot Jupiters (10¹⁰ g s⁻¹, e.g., Vidal-Madjar et al. 2003; Murray-Clay et al. 2009; Linsky et al. 2010), but is sufficiently high to modulate the background coronal density (a stronger outflow will only intensify the modulations). For reference, the escape rate from Saturn is estimated to be between ${10}^{2}\mathrm{and}{10}^{4}\,{\rm{g}}\,{{\rm{s}}}^{-1}$ (Glocer et al. 2007). Future studies that focus on specific planetary systems will require a more detailed planetary and outflow description. Such details could be obtain by coupling the coronal model described here with a model for the planetary magnetosphere (such as code coupling in the Space Weather Modeling Framework, see Tóth et al. 2012).

We require at least 10 grid cells across the planetary body in order to properly resolve it well. Therefore, we use grid refinement with very high resolution around the planet so that the grid size near the planet is Δx ≤ 0.01 R_⋆. In cases where the planet is closer to the star, the initial spherical grid refinement is sufficient. When the planet is further out from the star, we add an additional ring of high resolution along the orbit of the plane. Due to the grid limitations, we set the planet size to be 0.3 R_⋆ which is roughly three Jupiter radii. We performed several tests that have shown that setting a smaller planet size would only require much higher resolution around it, but the results were similar up to R_p = 0.15 R_⋆, and with the radio modulations remain at a similar magnitude. This is because the modulations also occur due to the surrounding plasma around the planet and not only by the planet itself.

We limit ourselves to the case of steady-state solutions, with a stationary planet, that are viewed from different angles within the orbital plane to mimic the orbital phase. Nevertheless, the simulated radio wave modulations from these static cases can only be enhanced by time-dependent effects due to the extra contributions by the dynamic interaction of the planetary magnetosphere with the stellar corona. Thus, here we provide a lower limit for the radio wave modulations. It is possible that the modulations of the radio waves induced by the planet have a similar magnitude to the stellar variation of the background radio flux, just like a dip in the visible flux might be attributed to a starspot instead of a planet transit.

A new tool to create synthetic radio images has been recently added to BATS-R-US (Benkevitch et al. 2010, 2012; Moschou et al. 2018). The tool accounts for the free–free Bremsstrahlung radiation that is created in the corona, and propagates through the non-uniform density of the circumstellar medium (e.g., Kundu 1965; Oberoi et al. 2011; Casini et al. 2017; Mohan & Oberoi 2017). The wave refraction depends on the local plasma density and the wave frequency. Thus, the radio waves of a given frequency propagate along curved rather than straight lines. The new radio image tool performs ray-tracing of the curved (i.e., not straight Line-of-Sight) propagation of the waves for a particular frequency, and calculates the integrated intensity of the radio wave at the end of the ray path, at a given pixel on the observing plane. The collection of all the pixels provides a radio image for the particular frequency.

The intensity of each pixel, I_ν, for a given frequency, ν, is the integral over the emissivity along the ray. Thus, the intensity is given by

$$\begin{eqnarray}&&{I}_{\nu }=\int {B}_{\nu }(T){\kappa }_{\nu }{ds}.\end{eqnarray} \tag{ 1 }$$

For Bremsstrahlung emission, where hν ≪ k_BT, the Planckian spectral blackbody intensity is (Karzas & Latter 1961)

$$\begin{eqnarray}&&{B}_{\nu }(T)=\displaystyle \frac{2{k}_{B}{T}_{e}{\nu }^{2}}{{c}^{2}}\end{eqnarray} \tag{ 2 }$$

with k_B being the Boltzmann constant, T_e is the electron temperature, and c is the speed of light. The absorption coefficient, κ_ν is

$$\begin{eqnarray}&&{\kappa }_{\nu }=\displaystyle \frac{{n}_{e}^{2}{e}^{6}}{{\nu }^{2}{({k}_{B}{T}_{e})}^{3/2}{m}_{e}^{3/2}c}\langle {g}_{{ff}}\rangle \end{eqnarray} \tag{ 3 }$$

Here n_e is the electron number density, m_e is the electron mass, e is the electron charge, and $\langle {g}_{{ff}}\rangle$ is the Gaunt factor, which is assumed to be equal to 10 (Karzas & Latter 1961).

The ray-tracing path for a given angular frequency, ω = 2π ν, is defined by the radio wave refraction from one grid cell to the next one. The index of refraction is related to the dielectric permittivity, , and is given by

$$\begin{eqnarray}&&{n}^{2}=\epsilon =1-\displaystyle \frac{{\omega }_{p}^{2}}{{\omega }^{2}},\end{eqnarray} \tag{ 4 }$$

with ${\omega }_{p}^{2}=4\pi {e}^{2}{n}_{e}/{m}_{e}$ being the plasma frequency (rad s⁻¹). Assuming plasma quasi-neutrality, where the densities of electrons and ions are the same, we can write the mass density as ρ = m_pn_e with m_p being the proton mass. Thus we have

$$\begin{eqnarray}&&\epsilon =1-\displaystyle \frac{\rho }{{\rho }_{c}},\end{eqnarray} \tag{ 5 }$$

where ${\rho }_{c}={m}_{p}{m}_{e}{\omega }^{2}/4\pi {e}^{2}$ is the critical plasma density at which the refraction index equals zero and the wave cannot be transmitted through. Equation (5) demonstrates that higher frequency radio waves can penetrate deeper into the solar atmosphere, where the density is higher. Thus, the synthesized images for higher frequencies capture more detailed structures, such as active regions, in contrast to the lower frequencies synthesized images, that capture the lower density regions at the top of the corona (as demonstrated in Moschou et al. 2018).

When the planet resides at a = 0.028 au, the intensity modulations are driven by the strong star–planet interaction since the planet is located at or close to the Alfvén surface (see top panels in Figure 2). While the majority of the background stellar radio emission comes from the dense, helmet streamer regions that face the observer, in the case of a = 0.028 au, the edge of the helmet streamer is disrupted by the interaction with the planet. This disruption could involve plasma compression, mixing of coronal and magnetospheric plasmas, as well as creating plasma cavities. As a result, the radio intensity in different bands can be modulated by this local interaction region. Since the helmet streamers emissions are blocked during transit, the contribution to the background emission depends on the emissions generated at the interaction region.

Looking closely at the synthetic radio images and Figure 2, we find that for the case of planetary field of 1 G, the intensity contribution of the interaction region in both the 10 and 30 MHz bands is greater than the ambient intensity from the helmet streamers. Thus, there is a significant intensity increase in these bands during transit. Similar trend is found in the 30 MHz intensity for the case with planetary field of 0.3 G. However, the intensity contribution in the 10 MHz with this field strength is found to be negligible comparing to the ambient helmet streamers intensity, which is blocked during transit. Thus, the overall trend we find for this case is an intensity drop in the 10 MHz during transit. In the higher frequency bands, even in the 10 GHz band, we find a small but noticeable drop of about 10% in the intensity due to the blockage of the helmet streamers by the planet.

An interesting feature in the a = 0.028 au and a weak planetary field case is that the emission peaks are slightly shifted from the transit point by about ∼15°. This shift seems to happen due to a late but still strong interaction between the planetary magnetosphere and the stellar corona beyond the transit point, and the fact that the interaction between the planetary and stellar plasma occurs at or within the Alfvénic point. This a-symmetry is not visible in any of the other cases.

For the a = 0.042 au cases, we find similar trends to that of the a = 0.028 au case, but with a significantly reduced magnitude, within the 10% range of intensity increase or decrease. While the planet and the helmet streamers can still interact at this orbit, the interaction is much weaker than the case with a = 0.028 au.

For the a = 0.056 au cases, almost no signs of star–planet interaction is visible, with the exception of a small increase in the 10 MHz intensity. This increase is slightly larger for the case with a stronger planetary field, where the magnetosphere is larger comparing to the 0.3 G planetary field case, and a stronger plasma compression at the magnetopause. Interestingly, the transit shadowing of the helmet streamers emissions in the 100 MHz is still visible, with a decrease of almost 20% in transit. This particular feature, which has significant modulation beyond the very low frequency range, could potentially be observed.

At larger orbital separation of a = 15 R_⋆ and planetary field of 0.3 G, there is a significant enhancement in the 10 MHz band. This enhancement is due to a cavity created in front of the planet (see bottom-left panel of Figure 2), which seems to compress the top of the helmet streamer and increase the emissions in this band. This cavity does not appear in the case with a stronger planetary field, due to the increase in plasma density near its magnetopause. The intensity of the 30 MHz in the case of weaker planetary field is reduce in the form of two "wings" of the light curve. This pattern suggests that the 30 MHz band represents the flanks or edge of the planetary magnetosphere, but it is now shadowing the 30 MHz ambient emissions. The magnetosphere for the case with planetary field of 1 G is larger. As a result, the density of the magnetospheric plasma in this case is slightly lower, leading to a smaller enhancement in the 10 MHz band. The ambient emissions in the 30 MHz are still shadowed by the planet, but the magnetospheric impact on the shape of the light curve is not noticeable in the case of the stronger planetary field. The transit intensity drop in the 100 MHz is still noticeable at this longer orbit.

Figure 6 compares the modulations for the cases where the planetary field is ±1 G (with respect to the stellar dipole polarity). It can be clearly seen that when the planet is further away from the star, the results are not affected at all by the polarity of the planetary field. However, when the planet is close to the star, we find a greater difference between the two cases. The planetary field polarity has a stronger impact on the planetary density profile at closer orbits, while the effect is significantly reduced at further orbits. This happens since at closer orbits, where the planet is located at or below the Alfvén point, the star–planet interaction is more sensitive to the polarity of the planetary field.

Zoom In Zoom Out Reset image size

In particular, the strong enhancements in the 10 and 30 MHz bands for the case where the planetary field polarity is the same as the star are generated by local plasma enhancements through the star–planet interaction. When the field polarity is opposite, magnetospheric plasma is allowed to escape, so the local density enhancements are reduced, resulting in a suppression of the intensity enhancements in the low frequency bands.

Finally, we investigate how a stronger stellar magnetic field affects the modulations of the coronal radio emission induced by the planet. This is an important factor when considering M-dwarf stars, which are known to potentially have extremely strong magnetic fields of up to few kG (Reiners & Basri 2010), and are much more magnetically active than the Sun. These strong stellar fields could potentially produce strong coronal radio emissions (such as in V374 Peg, Llama et al. 2018). We repeat the simulations with a stellar dipole field of 100 G, and the comparison is shown in Figure 7. In general, increasing the stellar dipole strength leads to a larger size of the coronal loops so the helmet streamers extend to a greater distance comparing to the case with a weaker stellar dipole field. The coronal plasma is trapped in these larger closed loops, resulting in an overall enhancement of the coronal density, and a reduction of the density drop with radius. When a planet resides at certain distance from the star, it is surrounded by certain plasma density. By increasing the stellar field strength, we practically move the planet to a higher density region than before. This behavior is clearly seen in Figure 7.

Zoom In Zoom Out Reset image size

The case of a = 0.028 au is initially located (with stellar field of 10 G) at the top of the helmet streamers, and experiences strong interaction with the lower density plasma. This interaction is visible in the 10–30 MHz bands. When we increase the stellar field to 100 G, the planet is surrounded by, and interacts with a much more dense plasma. As a result, the significant modulations in the low frequency range is reduced, and the modulations are more visible in the 100 MHz band. When the planet is at a = 0.070 au, the increase in the stellar field strength practically move the planet inwards, and the intensity modulation trends resemble the case for weaker stellar field, and a planet at a = 0.028 au (top-right panel of Figure 2).

Our simplified approach here uses a idealized, dipolar stellar magnetic field, and a static, steady-state solutions for the structure of the stellar corona with a planet embedded in it. The phase variations are mimicked by viewing the static, three-dimensional solution from different angles. A number of factors, if included, could immediately provide additional variability of the radio intensity.

In reality, the structure of the stellar magnetic field and the stellar corona is more complex than the axisymmetric dipolar geometry we use here, and the corona, which hosts the planet, has sectors with different plasma properties along the planetary orbit. Thus, one should expect variations in the plasma density along the planetary orbit as the planet crosses from one plasma sector to another. For short orbits of few days, the size of the coronal sectors should be of the order of the spatial coverage of a large helmet streamer over the orbital plane. Such variations could be captured by modeling more realistic stellar systems. Simulations using ZDI magnetic maps (Donati et al. 1989) have been done using our code (e.g., Cohen et al. 2010, 2014; Alvarado-Gómez et al. 2016a, 2016b; Garraffo et al. 2016, 2017; Pognan et al. 2018) to simulate the coronae and winds of specific stars. Synthetic radio images could be produced for these more realistic coronal solutions in the same manner of the results presented here.

The orbital motion of the planet with respect to the ambient coronal plasma, not included in our static simulation, could be included in a time-dependent model as presented in Cohen et al. (2011b, 2011c). We expect that the implementation of the planetary orbital motion will enhance the star–planet interaction, and potentially increase the modulations of the radio intensity given the particular geometry of the planetary and stellar magnetic fields (see discussion in Section 4.1).

An additional factor to consider in a realistic case is the temporal variations of the radio signal itself. Variations in the intensity of the radio signal can be due to photospheric convective and diffusive motions, coronal waves, stellar wind, and stellar rotation. All these processes create temporal density variations in the medium through which the radio waves propagate. These variations extend to a wide range of temporal scales—seconds, minutes, hours, and days. Of course, it is important to identify these variations in the radio signal in order to isolate the variations that are associated with the planetary orbital motion, which should be of the order of a few hours or more for a planetary orbit of few days. It is also important to remember that radio observations are not the typical time-series of a source, but it is an observation that is quite diffusive in the spatial manner and it is sometimes hard to identify the exact source location due to refraction and scattering effects.

It is reasonable to assume that the short time variations of the order of less than an hour are associated with plasma variations in the low corona and close to the photosphere. Radio emissions of the dense plasma associated with these regions appear in the higher frequency range of the radio spectrum at 1 GHz or above. Thus, we expect the ambient large-scale variations of the radio signal, originating from higher altitudes and in the range below 1 GHz, to have longer timescales than minutes. The coronal helmet streamers tend to stick around for days and even months. Therefore, signal variations in the 10–100 MHz range, which originates from the helmet streamers, should have a similar, longer timescale comparing to the variation timescales at higher frequencies.

Stellar flares can also dramatically disrupt the coronal density structure and affect the radio signal. There is an active search for such radio signal in an attempt to observe stellar Coronal Mass Ejections (CMEs, e.g., Villadsen 2017; Crosley & Osten 2018). Of course, stellar flares can mask the radio modulation created by an exoplanet. However, stellar flares are typically visible in other wavelengths, and a flare-like different emission mechanism would be easy to distinguish by using, for example, polarization measurements. Therefore, we should have a pretty good idea whether or not a flare occurs during the time of radio exoplanet observation so that period could be excluded to avoid uncertainties.

Our new radio tool could provide predictions for the frequency range at which it is most likely to detect a signal for specific targets. Such predictions can be used by observational radio facilities, such as LOFAR,⁷ MWA,⁸ Effelsberg⁹ and VLA.¹⁰ It is important to note that here we investigate planets with a very short orbital period to maximize the modulation effect, and we find modulation of 50% or more in some cases. While we expect that planets with a larger orbital separation will have a much smaller modulation effect, the modulations can still be of the order of few percent. In addition, as described in Section 4.2, in systems where the stellar field is much stronger, such as M-dwarf systems, the reduction due to the orbital separation can be compensated by the increase in stellar field and coronal density.