Flares, Rotation, Activity Cycles, and a Magnetic Star–Planet Interaction Hypothesis for the Far-ultraviolet Emission of GJ 436

We identified 14 flares in the data, which are shown in the context of the entire data set in Figure 2. Table 1 gives the properties of these flares in the FUV₁₃₀ band, enabling a comparison to L18, as well as the summed line fluxes that are used to identify the flares and are shown in the figures. Most flares exhibit distinct peaks. Event 7 (and possibly also event 10) appears to be the decay phase of a flare that began prior to the start of an exposure. Events 5 and 14 do not exhibit clear peaks and, despite passing the statistical cut, could be false positives. We cannot definitively explain these flux increases. One possibility is that they are manifestations of magnetoacoustic waves or wave interference, which produce minute timescale variations in the transition region emission on the Sun, although it is unclear whether these variations significantly modulate disk-integrated emission (Sangal et al. 2022). Similar variations with greater amplitudes that stand out well above the noise have been observed on a young M star by Loyd et al. (2018a). Another possible explanation is that these anomalies are conglomerations of weak flares.

Table 1. FUV Flares Detected from GJ 436

Summed Lines (Used to Identify Flares and Plotted in Figures 2 and 7)				FUV130 Band (per L18)
δ	E	Fpeak	$\tfrac{{F}_{\mathrm{peak}}}{{F}_{q}}$ a	δ	E	Fpeak	$\tfrac{{F}_{\mathrm{peak}}}{{F}_{q}}$ a	tpeak	No. b
s	1028 erg	10−14 erg s−1 cm−2		s	1028 erg	10−14 erg s−1 cm−2		MJD
1514 ± 130	3.76 ± 0.26	4.85 ± 0.67	23.2 ± 6.0	745 ± 120	3.76 ± 0.49	6.70 ± 0.82	16.1 ± 5.5	58142.88	11
1419 ± 118	3.52 ± 0.28	1.57 ± 0.38	8.2 ± 2.7	863 ± 129	4.35 ± 0.61	1.93 ± 0.54	5.4 ± 2.3	58143.15	12
1337 ± 112 c	3.68 ± 0.27	3.05 ± 0.53	13.6 ± 4.0	660 ± 90	3.82 ± 0.51	3.63 ± 0.68	8.1 ± 1.8	58110.89	7
1018 ± 117	2.29 ± 0.24	4.77 ± 0.66	25.1 ± 5.2	594 ± 111	2.81 ± 0.52	5.99 ± 0.81	15.4 ± 2.7	58108.72	6
1001 ± 95	2.92 ± 0.25	2.23 ± 0.46	9.7 ± 2.9	591 ± 97	3.16 ± 0.50	2.66 ± 0.65	6.7 ± 2.1	58177.33	13
884 ± 99	2.43 ± 0.26	1.44 ± 0.37	7.0 ± 2.4	562 ± 106	3.25 ± 0.61	2.55 ± 0.82	6.0 ± 2.0	58111.42	8
776 ± 87	2.14 ± 0.22	1.49 ± 0.37	7.2 ± 2.5	357 ± 86	2.07 ± 0.49	2.26 ± 0.59	5.4 ± 1.5	58111.70	9
663 ± 98	1.15 ± 0.17	1.88 ± 0.42	13.3 ± 3.0	732 ± 173	1.84 ± 0.39	2.26 ± 0.58	11.2 ± 5.2	56101.31	1
646 ± 92	2.14 ± 0.28	1.28 ± 0.34	5.4 ± 2.3	434 ± 89	2.48 ± 0.51	2.21 ± 0.54	5.4 ± 1.3	58137.65	10
627 ± 74	2.51 ± 0.28	1.38 ± 0.36	4.9 ± 1.9	527 ± 75	3.54 ± 0.50	2.23 ± 0.54	4.8 ± 1.1	57199.12	5
626 ± 85	1.82 ± 0.24	1.13 ± 0.33	5.4 ± 2.0	460 ± 104	2.47 ± 0.54	2.16 ± 0.60	5.6 ± 1.9	58177.46	14
476 ± 71	0.83 ± 0.12	2.08 ± 0.43	14.6 ± 3.0	280 ± 131	0.70 ± 0.28	2.69 ± 0.61	13.2 ± 5.9	56101.37	2
395 ± 64	1.58 ± 0.24	1.68 ± 0.40	5.8 ± 2.1	78 ± 67	0.53 ± 0.45	1.97 ± 0.59	4.3 ± 1.3	57199.10	4
317 ± 55	1.27 ± 0.19	2.22 ± 0.44	7.3 ± 2.5	198 ± 47	1.33 ± 0.31	3.00 ± 0.57	6.1 ± 1.2	57198.99	3

Notes.

^a The ratio of the peak flux to the quiescent flux. ^b The chronological order of the flares, corresponding to the numerical labels in the figures. ^c The exposure began after the start of the flare, as indicated by the lack of any data points below the quiescent level prior to the peak. The energy and equivalent duration are lower limits.

Download table as:  ASCIITypeset image

The flares exhibited by GJ 436 are frequent, yet they are of unusually low energy. Figure 4 shows the cumulative flare frequency distribution (FFD) of the observed flare energies and equivalent durations (a metric of flare energy that is normalized by the star's quiescent emission) in comparison with FFDs that are made by combining observations of 10 M dwarfs in the same FUV₁₃₀ band (L18). Although the L18 FFDs incorporate the 2012 and 2015 epochs of the GJ 436 data, they make up <10% of the data, so they will not significantly bias the comparison. The flattening of the GJ 436 FFD at low energies represents the event detection limit.

Zoom In Zoom Out Reset image size

The GJ 436 FFD does not appear to follow the "M-dwarf average" power law from L18, since a cliff in the FFD is exhibited toward large equivalent durations (δ) and energies. Within the cumulative exposure of 0.88 days, the equivalent duration and energy of the strongest expected event would be ≈3000 s and ≈3 × 10²⁹ erg (where the L18 FFD lines intersect the 0.88 day⁻¹ occurrence rate in Figure 4), but the strongest observed event is several times weaker. The observed flare rate is not affected by exposure gaps, since flares are effectively random in time (Wheatland 2000). The observed energies can be underestimated when gaps cut off events (such as event 7), but the bias that this introduces is below the statistical uncertainty resulting from a small sample size (see the injection recovery tests with gappy data in Appendix C of L18).

To test the possibility that the FFD cliff is the result of a small sample size, we randomly drew flares from the L18 equivalent duration FFD for comparison. In 900,000 trials, 10% yielded no flares with equivalent durations beyond the cliff. However, of this 10%, most trials produced fewer total flares. When we considered only trials that yielded at least as many flares with δ > 500 s as GJ 436 (180,000 trials yielding ≥8 flares), only 0.02% produced events that were no more extreme than those observed from GJ 436. In other words, if GJ 436's FFD were to be consistent with the FFDs of other M dwarfs, it is highly unlikely that it could yield as many flares as observed, yet produce none with equivalent durations >860 s. This implies that GJ 436's FFD is steeper than is typical. A power-law fit to the 8δ > 500 s events yields an exponent < −2, well below the 0.76 ± 0.1 that is typical for M dwarfs (L18), although the fit is poor due to the low number of events. We hypothesize that this unusual FFD is the result of a magnetic SPI (Section 3.5).

To investigate whether GJ 436 is able produce more energetic flares over longer timescales, we searched 43 days of optical observations by the Transiting Exoplanet Survey Satellite (TESS). The presence of a flare in those data would have been evidence against the cliff in UV flare energies being a real feature. However, we found no flares in the TESS data. Because flares produce very low contrasts at optical wavelengths (e.g., MacGregor et al. 2021), the possibility remains that flares with large UV energies occurred, but their optical signals were buried in the noise of the TESS time series.

Based on our fits to the data sampling two and a half stellar rotation periods during the 2017/2018 epoch, we detected stellar rotational modulation at >2σ in each band and 4.3σ in the summed line flux, with typical amplitudes spanning 5%–13%. Figure 5 shows these fits and Table 2 gives the fit parameters. MCMC chains of the fits and corner plots are available upon request to the corresponding author.

Zoom In Zoom Out Reset image size

Table 2. Fits to the Variability in the FUV Emission Lines

Band	Mean	Cycle	Cycle	Rotation	Rotation	SPI	Excess
	Flux	Amplitude	t0 a	Amplitude	t0 a	Amplitude	Noise b
	10−16 erg s−1 cm−2	%	yr	%	d	%	%
Summed Lines	${31.92}_{-0.74}^{+0.59}$	${37.5}_{-3.0}^{+3.3}$	${3.905}_{-0.075}^{+0.079}$	${8.1}_{-1.6}^{+2.1}$	−2.1 ± 1.5 c	<9.2	${10.8}_{-1.1}^{+1.5}$
C ii	${10.01}_{-0.29}^{+0.35}$	${41.9}_{-3.8}^{+6.2}$	3.86 ± 0.10	${9.5}_{-2.2}^{+2.6}$	$-{3.7}_{-1.1}^{+1.2}$	<7.9	${13.6}_{-1.6}^{+2.4}$
C iii	${5.25}_{-0.27}^{+0.32}$	${30.4}_{-9.5}^{+9.2}$	${4.04}_{-0.31}^{+0.18}$	${4.8}_{-2.3}^{+3.1}$	$-{1.7}_{-1.5}^{+1.4}$	<13	<12
Si iii	${3.55}_{-0.16}^{+0.22}$	${50.8}_{-9.5}^{+8.0}$	${4.02}_{-0.15}^{+0.13}$	${13.0}_{-3.7}^{+4.1}$	$-{1.2}_{-1.3}^{+1.2}$	<16	${21.5}_{-3.0}^{+4.1}$
Si iv	${4.74}_{-0.21}^{+0.31}$	${29.1}_{-7.1}^{+8.7}$	${3.66}_{-0.42}^{+0.17}$	${12.5}_{-3.1}^{+3.0}$	${0.1}_{-1.2}^{+1.1}$	<11	${8.1}_{-5.3}^{+3.6}$
N v	${7.89}_{-0.27}^{+0.23}$	${30.3}_{-4.0}^{+6.8}$	${3.95}_{-0.15}^{+0.11}$	${6.7}_{-1.6}^{+1.8}$	$-{3.0}_{-1.2}^{+1.1}$	<4.8	${6.9}_{-3.6}^{+1.9}$
Pseudocontinuum	${14.4}_{-5.1}^{+1.5}$	${89}_{-13}^{+174}$	${4.538}_{-0.065}^{+0.052}$	${11.7}_{-4.5}^{+7.8}$	$-{4.3}_{-1.4}^{+1.8}$	<10	${14}_{-8}^{+12}$
FUV130	${48.1}_{-3.3}^{+1.8}$	${47}_{-8}^{+16}$	${4.24}_{-0.11}^{+0.06}$	${8.3}_{-2.0}^{+2.4}$	−3.5 ± 1.2	<7.8	${10.7}_{-1.9}^{+2.1}$

Notes.

^a The time at which the phase of the sin function reaches zero after 2018 January 1 00:00:00 UT (58119.0 MJD). ^b Additional white noise hyperparameter for the 2018 epoch, excluding flares. ^c Applied as a prior to the remaining fits.

Download table as:  ASCIITypeset image

Taking into account the inclination of the star, the rotational variability can provide some insight into the net contrast of active regions. Bourrier et al. (2022) measured the star's spin axis as being inclined by $35.7{^\circ }_{-7.6}^{+5.9}$ relative to the line of sight. ¹¹  At this inclination, latitudes beyond ± 543 are always visible for one pole and invisible for the other. Granzer et al. (2000) found that active regions mostly occur within ±543 for slowly rotating M dwarfs. Therefore, most of GJ 436's active regions likely contributed to the observed variability, albeit with an apparent area foreshortened by the inclination.

dS19 analyzed the same FUV observations for rotational modulations. They only found the rotational modulation of the individual emission lines to be significant if they included the 2015 observations under the assumption, justified by optical observations (Bourrier et al. 2018), that the signal phase did not shift between the 2015 and 2017/2018 epochs. In our analysis, summing the emission line fluxes resulted in a >4σ measurement of the rotation signal amplitude, without needing the 2015 data. Regardless, we consider the rotational modulation of the FUV emission as more likely present than absent, given the signs of magnetic activity on GJ 436 and the rotational modulation of the UV emission due to magnetic activity on the Sun (Toriumi et al. 2020). The ds19 fits yielded rotational amplitudes of roughly 15% in C ii and Si iii and 10% in N v. Our fits yielded lower levels of rotational variability (Table 2), likely because, in contrast with dS19, we did not include the higher fluxes of the 2015 epoch when fitting for rotation, due to concerns that the phase and amplitude could have shifted between the epochs.

We tested whether the FUV data could recover the rotation period of the star if we included the rotation period as a free parameter in the variability fit. In this case, the MCMC sampler identified multiple peaks in the posterior for the rotation period, with the two strongest at ${21.3}_{-1.1}^{+0.2}$ days and 6.69 ± 0.04 days (Figure 5, second panel), and only a weak peak in the vicinity of 44 days. The amplitudes of the short-period signals were greater than the fixed-period fit, but within 2σ. The multiple modes could be harmonics of the longer optical period, perhaps indicating that multiple longitudes of enhanced activity are contributing to the FUV variability.

In comparison to the FUV emission, the variability in the broadband optical emission is ∼0.1%. Our sinusoidal fit to the optical data (Figure 3) yielded a rotation period of 44.09 ± 1.16 days, with the zero phase occurring at JD 2455484.47 ± 0.73 (decimal year 2010.7862 ± 0.0020), with an amplitude of 0.127% ± 0.013%, consistent with the results of Bourrier et al. (2018) and Lothringer et al. (2018). In contrast, rotation measurements made using Ca ii H and K and Hα suggest periods nearer to 40 days (Suárez Mascareño et al. 2015; Kumar & Fares 2023). A rotation period within the 40–45 days range is typical of a field-age M dwarf with a mass ≈0.45 M_⊙ (Bourrier et al. 2018; Newton et al. 2018). If the rotational variability is due to dark spots, the star should appear redder during optical lows, but Lothringer et al. (2018) could not detect this, suggesting that the variability could be dominated by faculae.

The phase difference between the UV and optical rotation curves encodes information about the surface structure of the stellar activity. If the rotation signals result from optically bright faculae that are cospatial with FUV plage, then the two should be in phase. Alternatively, if optically dark spots that are cospatial with FUV plage are responsible, then the two curves should be 180° out of phase.

In this case, however, a clear comparison of the UV and optical phases is precluded by the possibility of phase shifts in the optical rotation signal. The optical variability fit assumes a static phase over the 14 yr of data, a span that covers two activity cycles. Although the data, folded onto the best-fit period, show convincing modulation (Figure 3), hidden phase shifts over the 14 yr span could be present, particularly if the signal is dominated by only one or a few short periods of high variability.

To investigate possible shifts in the optical phase, we isolated 3 yr sections of the optical data, sampling the extrema of the activity cycle (2006–2008, 2010–2012, and 2013–2015). The fits to data near optical minima recovered the same 44 day period, with a larger amplitude near the 2014 minimum (0.27% versus 0.13%) and large uncertainties in phase allowing for shifts as high as 123° (1σ) from 2007 to 2014, which could affect an FUV–optical comparison. The fit to the data near the 2011 maximum did not clearly recover the 44 day period, yielding multiple modes in period and phase, all with amplitudes near 0.9%. The unclear period and lower-variability amplitude near the 2011 maximum could result from a lower activity level (fewer or weaker spots/faculae), a more homogeneous distribution of spots/faculae across the star, cancellations of spots by faculae, or any combination of the above. Since the 2017/2018 FUV observations occurred near the subsequent optical maximum, they could also be affected by reduced or more spatially homogeneous activity.

We also attempted to fit TESS data from the two available visits, sectors 22 and 49, observed in 2020 March and 2022 March, each only lasting for about half of the stellar rotation period. Separate fits to each did not reveal clear rotation signals or phase constraints. Higher-precision optical data are necessary to probe signal shifts over time.

Toriumi et al. (2020) provide solar context for how rotational variations, resulting from spots and faculae, manifest at optical versus UV wavelengths. By isolating periods where only a single active region was visible, they found that the variability in extreme UV (EUV) emission lines could either be in or out of phase with the optical variability. The determining factor was whether the brightened active region or its dimmed surroundings dominated. The dimming is caused by the heating of the surrounding plasma to higher temperatures, reducing the emission from lower-temperature lines. The variations in EUV emission also start earlier and end later than the optical variations, because the EUV emission originates above the photosphere, with the result that it is visible when photospheric active regions are hidden just beyond the limb. Finally, the EUV variations exhibit top hat–like shapes, due to low optical depths. GJ 436's FUV variability could exhibit similar effects, though, unlike the transition region FUV lines that we analyzed, the most analogous EUV lines that were analyzed by Toriumi et al. (2020) included substantial flux from hotter, less dense coronal plasma.

The activity cycle of GJ 436 was measured by Lothringer et al. (2018) as having a period of roughly 7.4 yr, with an amplitude of 5 mmag (0.5%), at optical wavelengths. Our analysis of the same data yielded a period of 7.75 ± 0.10 yr and an amplitude of 0.376% ± 0.015% (4 mmag), with a zero phase occurring at JD 2454826 ± 18 (decimal year 2008.982 ± 0.049). A long-term increase of 0.050% ± 0.003% over the 14 yr span of observations is also present. The difference of 0.35 yr between our cycle period and that of Lothringer et al. (2018) is likely a systematic uncertainty, due to subjective analysis choices (using flux versus magnitude, the functional form of the long-term trend, the functional form of the cycle variability, the inclusion and form of a jitter term, etc.). An additional search for the star's activity cycle signature in Ca ii H and K, Na i, and Hα activity indicators by Kumar & Fares (2023) yielded signals spanning 5.1–6.8 yr. All of the aforementioned values fall near the predicted value of 7.6 yr, based on the relationship between the activity cycle and the rotation period established by Suárez Mascareño et al. (2016). The simple sinusoidal shape of the activity cycle echoes that of field-age Sun-like stars, with weak faculae-dominated activity.

The three epochs of HST data span nearly a full activity cycle and show significant variability beyond what can be explained by stellar rotation, visit-to-visit systematic flux errors (<2%; James 2022), or excess noise (2%–4% when binned over visits of two to five orbits). Fitting a sinusoid to these data, fixed to the period of the optical cycle, yields amplitudes of 30%–50%. The amplitude of the cycle variability in the pseudocontinuum stands out at ${89}_{-13}^{+174}$%. We interpret this difference with caution because, on a pixel level, the flux in the pseudocontinuum was below the noise floor of COS and vulnerable to errors in the background subtraction and time-varying flat-field correction.

Similar to rotational variability, if optically dark, FUV-bright active regions explain this variability, then the optical and FUV activity cycles should be 180° out of phase (Reinhold et al. 2019). For GJ 436, the optical variability lags the FUV variability by 119° ± 6° (Figure 6). dS19 found that the variations in the Ca ii S-index were out of phase with the optical variations, though the constraint was not quantified.

Zoom In Zoom Out Reset image size

The intermediate FUV–optical phase offset could indicate that GJ 436 is undergoing a transition from spot-dominated to faculae-dominated activity. Sun-like stars exhibit similar intermediate-phase offsets between their activity cycles in Ca ii H and K and optical emission, as they transition from spot to faculae dominance (Reinhold et al. 2019) around a Rossby number of 1, similar to GJ 436's Rossby number of 0.74 (Loyd et al. 2021). A transitional state for GJ 436 could explain its mixed indications of spot- and faculae-driven variability.

Systematic unknowns could make a 180° FUV–optical phase difference possible. A global trend of about 7% yr⁻¹ to the FUV emission would suffice, as would allowing the FUV cycle to vary in amplitude, with it being about half the value in 2018 as it was in 2012. Another option would be to interpret the 2012 and 2015 observations as sampling rotational minima, meaning that the average fluxes were higher at those epochs, which would shift the cycle curve to earlier times. The star's increase in optical variability near 2015 (Section 3.2) could indicate larger FUV variations at that time. On the Sun, the amplitude of rotational variability waxes and wanes over the activity cycles at optical and EUV emission wavelengths (Fröhlich 2006; Woods et al. 2022; Llama 2023, private communication).

We included an excess white noise parameter ("jitter" term) in our fits, so that the MCMC sampler could explore an alternative to forcing the rotational, cyclic, or SPI signals to produce the observed point scatter. The sampler identified excess white noise, in addition to other variability, at levels of 10%–20% in all bands, except C iii. The amplitudes of excess noise and other variability were uncorrelated, i.e., increasing the noise did not require decreasing other variability, to obtain a similar goodness of fit.

Excess noise could come from a mixture of instrumental and astrophysical sources. A known instrumental source is the movement of fixed pattern noise when the position of the spectrum on the detector is shifted (by as much as 36 Å, or 10% of the wavelength range) between observations. However, instrumental noise of this kind is unlikely to account for the entirety of the measured excess noise, given the 2% visit-to-visit accuracy of COS (James 2022). Astrophysical sources could include undetected flares, "transition region bombs," shocks initiated by convective motions, and others (Loyd & France 2014).

Excess noise, instrumental or astrophysical, limits the detectability of planet b's transit in FUV emission lines, aside from Lyα (Lavie et al. 2017; Loyd et al. 2017; dS19). Integrated to a 1 hr cadence (and assuming a flat power spectral density), the noise level is 5.8% ± 0.8% in the emission from the least-ionized ion that we analyzed, C ii. Loyd & France (2014) estimated excess noise with an analytical method and a 60 s cadence, using only the 2012 data, and found values equating to 2.5% ± 1.4% in C ii at a 1 hr cadence.

To search for an SPI, our HST program (HST-GO-15174) was designed to add 8 points, sampling across the orbital phase of GJ 436 b. The data augmented the transit investigations of program HST-GO-14767. This sampling was motivated by a peak in the N v emission of GJ 436, in phase with the planet transit during the 2015 epoch.

We found no detectable evidence of an SPI, in the form of variations in the emission from GJ 436, in N v or any other line phased with the planetary orbit. Our simultaneous fits limit SPI variability in the summed line flux to <9.4% (2σ). The MCMC sampler strongly preferred sinusoidal solutions to a top hat. In Figure 7, we plot a model set at the 2σ upper limit on the SPI amplitude against the data, after subtracting the best-fit rotational modulation signal from those data.

Zoom In Zoom Out Reset image size

The upper limit on the SPI signal enabled us to place a corresponding limit on the strength of the planetary magnetic field, under certain assumptions.

An initial assumption is that the magnetic disturbances from the planet propagate back to the star. For this to occur, the Alfvén Mach number at the planet must be <1 (Lanza 2015). Saur et al. (2013) estimated the Alfvén Mach number of GJ 436's wind near GJ 436 b to be 1.01, very near the limit. This value relied on an empirically scaled estimate of the stellar wind speed at the planet of 235 km s⁻¹. In contrast, Bourrier (2016) estimate a lower value of 85 km s⁻¹, based on matching a model of the planetary outflow's dynamical interaction with the stellar wind to Lyα transit data. This would imply an Alfvén Mach number ≈0.3, making a magnetic SPI capable of propagating to the star. However, this still relies on an estimate of the stellar magnetic field that is based on a scaling relationship for Sun-like stars (Saur et al. 2013). A direct measurement of GJ 436's magnetic field is necessary to better gauge whether GJ 436 b orbits within a region where the Alfvén Mach number is <1.

A second assumption is that the interaction takes the form of the "flux tube dragging" scenario set forth by Lanza (2013). In this scenario, a persistent magnetic flux tube links the star and the planet. The orbital motion of the planet drags this tube through the ambient stellar magnetic field, triggering the stellar field to relax to lower-energy states, releasing energy in the process. This is the only SPI configuration that predicts energy release at a level that is consistent with past observations of stellar activity–related emission that is modulated at the orbital period of a close-in planet (Cauley et al. 2019). Note that for other SPI models, the relative orientation of the planetary and stellar fields can play a critical role (e.g., Strugarek 2016).

We followed the methodology of Cauley et al. (2019) to estimate an upper limit on GJ 436 b's magnetic field. The method requires estimates for a number of quantities. For the stellar field at the planet, we took the value of 667 nT, as estimated by Saur et al. (2013). We set the relative velocity between the planet and the star's magnetic field to be the same as the mean orbital velocity, estimated from the period and semimajor axis found by Lanotte et al. (2014). This simplification is justified, because the planet's orbital velocity is over an order of magnitude larger than the rotational velocity of the stellar magnetic field at the orbital distance of the planet. We used the planetary radius value of 4.10 ± 0.16 R_⊕ from Lanotte et al. (2014). Finally, to estimate the fraction of SPI power radiated in the summed FUV line flux, we used the flare energy budgets from panchromatic observations of stellar flares on the active M dwarf AD Leo, as reported by Hawley et al. (2003). These indicate that roughly 3% of the total energy that is emitted in a flare is emitted in the five lines that we summed. From the above values, we estimated an upper limit on the planetary magnetic field of ≲100 G. This limit is well above Jupiter's magnetic field strength and within the range of values that Cauley et al. (2019) estimated for several hot Jupiters.

Besides inducing a hot spot, a close-in planet could affect the activity level of its host in several ways. It could suppress activity by tidally disrupting the stellar convective dynamo (Pillitteri et al. 2014; Fossati et al. 2018). It could obscure activity by enshrouding the star in absorbing gas from a radiation-powered planetary outflow (Fossati et al. 2013; Staab et al. 2017). Alternatively (or even concurrently), it could enhance activity, by tidally inhibiting the spin-down of the host star (Poppenhaeger & Wolk 2014; Tejada Arevalo et al. 2021). GJ 436's activity, in the form of FUV emission lines, falls a factor of a few below the predicted level for an M star with a 44 day rotation period (Loyd et al. 2021; Pineda et al. 2021), but is not an outlier. For comparison, the early M star GJ 832 has a slightly shorter rotation period (37.5 ± 1.5 days; Gorrini et al. 2022) and a similar level of FUV emission, yet radial velocity (RV) surveys have not revealed a close-in massive planet (Bailey et al. 2009; Gorrini et al. 2022). We conclude that GJ 436's level of FUV emission is not strongly influenced by the presence of its warm Neptune, but that a magnetic SPI could nonetheless be present.

Although we did not detect an SPI in the form of orbit-modulated emission, we speculate that GJ 436's odd FFD could indicate a magnetic SPI. The enhancement of relatively low energy (equivalent duration ≈500 s) flares and the statistically significant lack of more energetic events could indicate that a magnetic SPI is triggering frequent releases of magnetic energy at low levels, thereby preventing the buildup of stored energy that is necessary to give rise to the larger flares exhibited by other M dwarfs.

The stars in the comparison sample from L18, which produced stronger flares, are less prone to magnetic SPIs. Within the sample, eight of 10 (six of six in the inactive sample) are known to host planets (including GJ 436; see Section 3.1). None are as massive and as close-in as GJ 436 b, though several come to within about a factor of 2 in terms of mass or distance. The closest comparison is GJ 581 b, a planet with $M\sin i=15.8\pm 0.3$ M_⊕ orbiting 0.040 au from its host (Robertson et al. 2014), versus GJ 436 b's ${25.4}_{-2.0}^{+2.1}\,{M}_{\oplus }$ mass and 0.031 au distance from its host (Lanotte et al. 2014). For the two stars without known planets, RV limits do not exclude planets like GJ 436 b (Bailey & White 2012; Kossakowski et al. 2022).

Consistent with an SPI hypothesis is the lack of flares during planetary eclipses, when a subplanetary hot spot on the star would be most likely to be invisible. The eclipse data are too brief for this difference to be statistically significant.

Theoretically, the increase in the time-averaged emission caused by SPI-triggered flares should also produce a detectable orbit-modulated signal. For this reason, we included the flux added by flares when fitting for a possible SPI signal. With the flares included, the MCMC sampler favors a slight increase in the emission within the −0.25 to 0.25 phase range, but it is marginal, with a likelihood ratio near unity. Within this phase range, the average flux (after removing the best-fit rotation and cycle variability) is 1.7σ higher than the two visits near the planetary occultation. With more data, changes in the flare rate with planet phase, variations in the time-averaged emission with planet phase, and FFD statistics could all lead to the confirmation or refutation of an SPI in this system. The present data hint that flare statistics might have the greatest quantitative power in identifying a magnetic SPI.

To place each form of variability in context, we compare their amplitudes side by side in Figure 8. Each type of variability has its own time structure, so each type has its own definition of "amplitude."

Zoom In Zoom Out Reset image size

Flares do not have a well-defined amplitude, since they occur stochastically, with large variations in peak flux (Figure 4). As a metric for comparison, we adopted the time-averaged contribution of the observed flares to the star's emission. This number represents a lower limit on the true contribution of flares to the star's time-averaged emission, because it omits both more frequent, weaker flares below the detection limit and larger, rarer flares that are not caught within the limited duration of the observations.

Excess noise, since we assumed it to be white noise, does not have time structure. However, better data could reveal structure on timescales of seconds to minutes that the present data do not resolve. As a metric for comparison, we adopted the 1σ level of the noise.

Our models for variability due to rotation, activity cycle, and SPI all have well-defined amplitudes. Since we did not detect SPI variations, Figure 8 shows only the 2σ upper limits on the amplitudes of the fits.

Activity cycles appear to be responsible for the greatest share of the variability observed across the observations. If the Sun is any guide, activity cycles will not only affect the average level of FUV emission, but also the frequency of flares and the amplitude of rotational variability. In other words, the amplitudes of other variability sources will wax and wane over the course of the stellar activity cycle.

The variability in most of the FUV bands that we considered is well correlated. Figure 9 plots the correlations between the fluxes of each pair of bands. For most line pairs, the points are consistent with a 1:1 correlation. The notable exceptions are N v and the pseudocontinuum. The variations in these bands, particularly at the high end (mostly due to flares), appear to be subdued, relative to the other lines, in line with the observations of M-dwarf flares (France et al. 2016; L18; France et al. 2020). Despite having a shallower slope for N v and the pseudocontinuum, correlations are still present at high significance in all cases, except the correlation between the pseudocontinuum and Si iv.

Zoom In Zoom Out Reset image size

We expect the EUV emission of GJ 436 to track the changes in the FUV line fluxes. Bourrier et al. (2021) conducted differential emission measure modeling of the similar early M1.5 dwarf GJ 3470 across several observation epochs. From 2018 to 2019, GJ 3470's FUV line fluxes increased by about 30%, while the flux in the 100–920 Å band predicted by differential emission measure distributions (DEMs) increased by about half of that, or 15%. In contrast, the Sun's emission variability across activity cycles increases toward shorter wavelengths (Woods et al. 2022), suggesting the opposite trend: that the variations in the EUV should exceed those in the FUV for GJ 436. France et al. (2018) also noted an effectively linear correlation between the Si iv and N v emission measured by HST and the 90–360 Å flux measured (at a different epoch) by the Extreme UltraViolet Explorer for a sample of 11 F–M dwarfs. Ultimately, simultaneous FUV–EUV observations of GJ 436 or similar stars are needed to quantify the relationship between FUV and EUV variability, but, at least qualitatively, EUV variations will track FUV variations. We leave DEM or stellar modeling reconstructions of GJ 436's EUV emission to future work.