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

Variability in the far ultraviolet (FUV) emission produced by stellar activity affects photochemistry and heating in orbiting planetary atmospheres. We present a comprehensive analysis of the FUV variability of GJ 436, a field-age, M2.5V star ($P_\mathrm{rot}\approx44$ d) orbited by a warm, Neptune-size planet ($M \approx 25\ M_\oplus$, $R \approx 4.1\ R_\oplus$, $P_\mathrm{orb}\approx2.6$ d). Observations at three epochs from 2012 to 2018 span nearly a full activity cycle, sample two rotations of the star and two orbital periods of the planet, and reveal a multitude of brief flares. Over 2012-2018, the star's $7.75\pm0.10$ yr activity cycle produced the largest observed variations, $38\pm3$% in the summed flux of major FUV emission lines. In 2018, variability due to rotation was $8\pm2$%. An additional $11\pm1$% scatter at 10 min cadence, treated as white noise in fits, likely has both instrumental and astrophysical origins. Flares increased time-averaged emission by 15% over the 0.88 d of cumulative exposure, peaking as high as 25$\times$ quiescence. We interpret these flare values as lower limits given that flares too weak or too infrequent to have been observed likely exist. GJ 436's flare frequency distribution (FFD) at FUV wavelengths is unusual compared to other field-age M dwarfs, exhibiting a statistically-significant dearth of high energy ($>4\times 10^{28}$ erg) events that we hypothesize to be the result of a magnetic star-planet interaction (SPI) triggering premature flares. If an SPI is present, GJ 436 b's magnetic field strength must be $\lesssim$100 G to explain the statistically insignificant increase in orbit-phased FUV emission. Erratum: Due to an arithmetic error, the published limit on the magnetic field strength is incorrect. The correct limit is $\lesssim$10 G.


INTRODUCTION
Corresponding author: [ astroparke@gmail.comGJ 436 is a nearby M2.5V, 0.45 M ⊙ star (Hawley et al. 2003;Knutson et al. 2011) well known for its lone Neptune-size planet (M = 25.4 +2.1 −2.0 M ⊕ , R = 4.10 ± 0.16 R ⊕ ; Butler et al. 2004;Gillon et al. 2007;Lanotte et al. 2014) that is actively losing atmospheric mass (Kulow et al. 2014;Ehrenreich et al. 2015).This planet is an upcoming target for atmospheric characterization by the James Webb Space Telescope (GTO programs 1177 and 1185, PI Greene), and the system has undergone extensive spectroscopic observation at far ultraviolet (FUV) wavelengths (France et al. 2013(France et al. , 2016;;dos Santos et al. 2019).These FUV observations provide an opportunity to investigate the variability of a fieldage (Torres et al. 2008), partially-convective M dwarf in high-energy emission across a wide range of timescales, from seconds to years.The data can probe variability as short lived as stellar flares and as long lived as stellar activity cycles.This information is relevant to exoplanets, where variations in FUV irradiation affect rates of photochemistry and heating in their atmospheres (e.g., Segura et al. 2007).Because M dwarfs are numerous and favorable for exoplanet observations (e.g., Shields et al. 2016), the constraints on FUV variability available for GJ 436 are broadly applicable to a large number of planets.
The GJ 436 system could be experiencing an interaction between the star and its single close-in planet.Many forms of star-planet interactions have been proposed, such as tidal suppression of activity, tidal spin up, and shrouding by escaped planetary gas (Cuntz et al. 2000;Fossati et al. 2013;Pillitteri et al. 2014;Poppenhaeger & Wolk 2014).Of particular interest is a direct magnetic interaction between the planet and host star that results in energy dissipation at the stellar surface, producing a hot spot that circumnavigates the star at the orbital period of the planet (Saur et al. 2013).For GJ 436, this is 2.64 d (Bourrier et al. 2018).Direct magnetic interactions themselves have many flavors, including particle precipitation from reconnection of stellar and planetary fields, Alfvén wave dissipation, and triggering reconnection of stellar fields (Lanza 2015).For magnetic disturbances to reach the star, the Alfvén Mach number at the planet must be < 1 (Saur et al. 2013).Though difficult to predict, this could well be the case for GJ 436 b (see Section 3.6).
A number of past works have identified evidence of this kind of star-planet interaction (SPI) in hot Jupiter systems (e.g., Shkolnik et al. 2008;Walker et al. 2008;Pagano et al. 2009;Cauley et al. 2019).Because the strength of the magnetic SPI, and hence the amplitude of its periodic signal, is controlled by the planetary magnetic field strength, this possibility provides the prospect of probing the magnetic field of GJ 436 b (Lanza 2015).With this in mind, we conducted targeted observations in 2017 and 2018 to sample the star's FUV emission across the planetary orbital period.
In this paper, we report on a comprehensive variability analysis of GJ 436's FUV emission and a search for a magnetic SPI.This work provides an independent replication of a similar analysis conducted as part of the transit study of dos Santos et al. (2019) (hereafter dS19).Their transit analysis confirmed the stability of the planet's remarkable Lyα transit and demonstrated no significant transit in metal lines.Meanwhile, their variability analysis yielded evidence for rotational modulation in C II, Si III, and N V emission stable over several years as well as magnetic cycle variations leading to optical, Ca II, and FUV variations consistent with changing spot coverage.In comparison to dS19, the present work adds a population study of the star's flares, details of variability fits, and the results of an SPI search.

METHODS
We analyzed all available observations of GJ 436 made with the Cosmic Origins Spectrograph (COS) aboard the Hubble Space Telescope (HST) that utilized the G130M grating.The observations originate from four programs that observed during three separate epochs spanning 5.5 years: HST-GO-12464 (PI France), 13650 (France), 14767 (Sing), and 15174 (Loyd).The data are aggregated under DOI 10.17909/6p65-wg08.Details of the observations are available at the DOI link where they can be downloaded.
The G130M grating covers wavelengths 1150 − 1450 Å at a resolving power of 12,000-16,000.The detector is a photon counter, yielding event lists that can be binned arbitrarily in wavelength and time within the limits of the detector resolution.From these lists, we created flux-calibrated lightcurves following the methodology presented in Loyd & France (2014) and Loyd et al. (2018b;hereafter L18). Figure 1 shows the observed band and Figure 2 indicates the quantity and duration of the data originating from the various observing programs.
We isolated flux within several bands to analyze for variability.The largest of these is the F U V 130 band used in the flare analysis of M dwarf stars by L18. Figure 1 plots the spectrum of GJ 436 within the F U V 130 band, created by coadding all exposures.As in L18, regions masked in gray were excluded from the band because they are contaminated by geocoronal airglow emission and, in the case of the mask near 1300 Å, are also affected by a gap between the two butted detectors used by COS.The coadded spectrum shown in Figure 1 does not exhibit this gap because it was eliminated by dithering.We also analyzed emission from ±100 km s −1 bands covering five strong emission lines and multiplets, C III, Si III, N V, C II, and Si IV, labeled in Figure 1 log 10 (K), within the stellar transition region (Dere et al. 1997;Del Zanna et al. 2021).For multiplets, we summed flux from all components.As a final band, we summed flux between the lines, labeling this "psuedocontinuum" since it is comprised of both continuum sources and weak or undetected emission lines.
We identified flares in the data using the method described in L18, with some modifications.The algorithm identifies flares by isolating runs (series) of points above the mean value of the lightcurve that have an integrated energy at least 5σ larger than the average run.It uses a maximum likelihood estimate for the mean, masking any previously identified flares.Chunks of data with gaps longer than 24 h are processed separately, allowing for, e.g., variations in the mean due to rotation.The process of estimating the mean and identifying flares is iterated until convergence.Newly identified flares are masked with each iteration to mitigate the upward bias in the mean produced by flares.
Rather than using the F U V 130 flux to identify flares as in L18, we used the summed emission line fluxes.These lines are more sensitive to flares than the pseudocontinuum portions of the F U V 130 band.Including the pseudocontinuum added substantial noise without substantially increasing the signal, reducing the sensitivity to flares.We also increased the span of data masked around a flare to start 120 s ahead of and extending to 3 times the length of the data identified as a flare.
After identifying flares in the summed line flux, we measured the equivalent duration, energy, and peak flux of each flare in the F U V 130 band to enable direct comparisons to the results of L18.Equivalent duration is the time the star would have to spend in quiescence to emit the same energy in the same band as was emitted by the flare.
Having cleaned the data of flares, we binned the remaining data to ≈ 10 min and simultaneously fit for variability due to stellar activity cycles, rotation, a possible SPI , and a jitter term.We only used the 2017/2018 epoch of data to constrain rotation and SPI signals because these were the only data that sampled across the planetary orbital period and stellar rotation period without large gaps where phase and amplitude could have shifted.We did not mask or otherwise account for planetary transits captured by the 2015 and 2017/2018 epochs because there is no evidence of a transit in the FUV emission lines that we analyzed (Loyd et al. 2017;Lavie et al. 2017, dS19).We used the MCMC code emcee1 (Foreman- Mackey et al. 2013) to sample the posterior distribution of the fit parameters using 100 walkers to a factor of at least 100× the autocorrelation length for the samples of each parameter, taken as the median of the values for all walkers.
The stellar rotation and activity cycle models were simple sine functions allowed to vary in amplitude and phase but fixed to a rotation period of 44.09 d and a cycle period of 7.75 yr as measured from the optical data (described below).We fit this model to the data integrated to 10 min bins with flares masked.After fitting the summed line fluxes, we used the posterior on phase from that fit as a prior for the other, lower SNR bands to better constrain their variability amplitude.The jitter term was Gaussian white noise added in quadrature to the measurement uncertainty of each data point.This quantity effectively represents how much more scatter was present in the lightcurve than was expected from measurement uncertainties and is analogous to the excess noise metric of Loyd & France (2014).However, we used longer time bins of 10 min instead of the 1 min bins used by Loyd & France (2014) to offset the faintness of GJ 436, which was the fifth faintest star among the 42 analyzed in Loyd & France (2014).We allowed the standard deviation of the added noise to vary between epochs, resulting in three separate white noise parameters.
To fit for a possible magnetic star-planet interaction signal, we assumed a truncated sine model as done for Ca II H & K in Shkolnik et al. (2003Shkolnik et al. ( , 2008)).This model approximates the appearance of a bright, optically-thick spot of constant flux with a viewing angle that shifts as it traverses the visible stellar hemisphere.However, FUV line emission within a hot spot might not be optically thick.If emission originates from an optically thin slab of plasma, and is emitted isotropically, then the flux would be constant once the slab emerged from behind the stellar disk (e.g., Toriumi et al. 2020).To allow for this possibility, we added a shape parameter, α, to the truncated sin model, yielding a model of the form where A is amplitude, P orb is the planetary orbital period of 2.64 d (Bourrier et al. 2018), and Φ is a phase offset.As α → 0, the function approaches a top hat, representative of the case where the hot spot emission is optically thin and the observed flux simply depends on whether the spot is visible or not.This ad hoc parameterization allowed the code to explore sinusoidal, top hat, and intermediate solutions.We fit this model to the data with flares included because they could be triggered by an SPI, simultaneous with rotation and cycle fits but with a separate jitter term.An initial run indicated only upper limits would result.These limits could be unreasonably large if we allowed phase to vary freely since the sampler could place the signal peak in a data gap.To provide a reasonable upper limit, we fixed the signal to be in phase with the orbit.Although phase offsets likely occur for an SPI-induced hot spot, they cannot be easily predicted (Lanza 2013), and we consider a fixed-phase fit reasonable for an order-of-magnitude upper limit.
To provide a consistent comparison to the optical variability of GJ 436, we also analyzed the photometry in the Strömgren b and y band filters described in Lothringer et al. (2018).We used the same techniques outlined above, but included rotation period and the slope and intercept for a linear rise in flux as additional free parameters.In contrast with Lothringer et al. (2018), we fit the (b + y)/2 photometry in flux rather than magnitude space.We set priors for the period of the rotation sinusoid to be < 365 d, the period of the activity cycle sinusoid to be > 365 d, and all amplitudes to be ≥ 0. These priors did not prove restrictive.Figure 3 shows the results of this fit.

Flares
We identified 14 flares in the data, shown in the context of the entire dataset in Figure 2. Table 1 gives the properties of these flares in the F U V 130 band, enabling a comparison to L18, as well as the summed line flux used to identify the flares and shown in the figures.Most flares exhibit a distinct peak.Event 7 (and possibly 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 a manifestation of magneto-acoustic waves or wave interference such as produce minute-timescale variations in transition-region emission on the Sun, though it is unclear if these variations significantly modulate disk-integrated emission (Sangal et al. 2022).Similar variations with a greater amplitude that stands out well above the noise were observed on a young M star by Loyd et al. (2018a).Another possible explanation is that these anomalies are conglomerations of weak flares.
The flares exhibited by GJ 436 are frequent, and 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 normalizes by the star's quiescent emission) in comparison with FFDs made by combining observations of 10 M dwarfs in the same F U V 130 band (L18).Although the L18 FFDs incorporate the 2012 and 2015 epochs of the GJ 436 data, they made up < 10% of the data, so will not significantly bias the comparison.The flattening of the GJ 436 FFD at low energies represents the event detection limit.
The GJ 436 FFD does not appear to follow the "M dwarf average" power law from L18, exhibiting a cliff in the FFD toward large equivalent durations (δ) and energies.Within the cumulative exposure of 0.88 d, the equivalent duration and energy of the strongest expected event would be ≈ 3000 s and ≈ 3 × 10 29 erg (where the L18 FFD lines intersect the 0.88 d −1 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).Observed energies can be underestimated when gaps cut off events (such as event 7), but the bias this introduces is below the statistical uncertainty resulting from a small sample size (see injection-recovery tests with gappy data in Appendix C of L18).
To test the possibility that the FFD cliff is a result of small sample size, we drew flares randomly from the L18 equivalent duration FFD for comparison.In 900,000 trials, 10% yielded no flares with equivalent durations beyond the cliff.However, within this 10%, most trials produced fewer total flares.If 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 no more extreme than those observed from GJ 436.In other words, if GJ 436's FFD were consistent with other M dwarfs, it is highly unlikely that it could yield as many flares as observed and yet produce none with equivalent durations > 860 s.This implies that GJ 436's FFD is steeper than typical.A powerlaw fit to the 8 δ > 500 s events yields an exponent < −2, well below the 0.76 ± 0.1 typical of M dwarfs  b Chronological order of the flare.Corresponds to the numerical labels in the figures.
c Exposure began after the start of the flare as indicated by a lack of any data points below the quiescent level prior to the peak.Energy and equivalent duration are lower limits.
(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 ( §3.6).
To investigate whether GJ 436 can produce more energetic flares over longer timescales, we searched 43 d 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 as 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.To highlight the cliff in the GJ 436 flares, the figure also shows random draws of even equivalent durations and energies assuming the same number of events above 500 s and 2.5 × 10 28 erg and a inverse power-law slope, α, equal to the median observed by L18.The gray horizontal line represents the limit where one event would be expected to occur within the cumulative observing time.

Rotation
Based on our fits to the data sampling 2.5 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.
Taking into account the inclination of the star, the rotational variability can provide some insight on the net contrast of active regions.Bourrier et al. (2022) measured the star's spin axis to be inclined by 35.7 •+5.9 −7.6 relative to the line of sight.2At this inclination, latitudes beyond ±54.3 • are always visible for one pole and invisible for the other.Granzer et al. (2000) found that active regions mostly occur within ±54.3 • for slowly rotating M dwarfs.Therefore, most of GJ 436's active regions likely contributed to the observed variability, albeit with an area foreshortened by the inclination.
dS19 analyzed the same FUV observations for rotational modulations.They found the rotational modulation of individual emission lines only to be significant if they included the 2015 observations under the assump-tion, 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 emission line fluxes resulted in a > 4σ measurement of the rotation signal amplitude without need of the 2015 data.Regardless, we consider rotational modulation of FUV emission as more likely present than absent, given signs of magnetic activity on GJ 436 and rotational modulation of 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 concern that phase and amplitude could have shifted between epochs.
We tested whether the FUV data could recover the rotation period of the star if we included rotation period as a free parameter in the variability fit.In this case, the MCMC sampler identified multiple peaks in the posterior for rotation period, with the two strongest at 21.3 +0.2 −1.1 d and 6.69 ± 0.04 d (Figure 5, second panel) and only a weak peak in the vicinity of 44 d.Amplitudes of the short-period signals were greater than the fixedperiod 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.3.66 +0.17 c Applied as a prior to remaining fits.
In comparison to FUV emission, variability in broadband optical emission is ∼0.1%.Our sinusoidal fit to the optical data (Figure 3) yielded a rotation period of 44.09 ± 1.16 d with zero phase occurring at JD 2455484.47 ± 0.73 (decimal year 2010.7862± 0.0020) and 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 & K and Hα suggest periods nearer to 40 d (Suárez Mascareño et al. 2015;Kumar & Fares 2023).A rotation period within the 40-45 d 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 cospatial with FUV plage, then the two should be in phase.Alternatively if optically dark spots 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 years 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 year 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-year sections of the optical data sampling extrema of the activity cycle (2006-2008, 2010-2012, and 2013-2015).Fits to data near optical minima recovered the same 44 d period, with a larger amplitude near the 2014 minimum (0.27 vs. 0.13%) and large uncertainties in phase allowing for shifts as high as 123 • (1σ) from 2007 to 2014 that could affect an FUV-optical comparison.The fit to data near the 2011 maximum did not clearly recover the 44 d 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, cancellation 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 fitting TESS data from the two available visits, sectors 22 and 49 observed in 2020 March and 2022 March, each lasting only about half of the stellar rotation period.Separate fits to each did not reveal a clear rotation signal or phase constraint.Higher precision optical data are necessary to probe signal shifts over time.Toriumi et al. (2020) provides solar context for how rotational variations, resulting from spots and faculae, manifest at optical versus UV wavelengths.Isolating periods where only a single active region was visible, they found variability in extreme UV (EUV) emission lines can be either in or out-of-phase with optical variability.The determining factor is whether the brightened active region or its dimmed surroundings dominate.The dimming is caused by heating of surrounding plasma to higher temperatures, reducing emission from lower temperature lines.Variations in EUV emission also start earlier and end later than optical variations because EUV emission originates above the photosphere, making it visible when photospheric active regions are hidden just beyond the limb.Finally, 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 we analyzed, the most analogous EUV lines analyzed by Toriumi et al. (2020) include substantial flux from hotter, less dense coronal plasma.

Activity Cycles
The activity cycle of GJ 436 was measured by Lothringer et al. (2018) at optical wavelengths to have a period of roughly 7.4 yr with an amplitude of 5 mmag (0.5%).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 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 year span of observations is also present.The 0.35 yr difference between our cycle period and that of Lothringer et al. ( 2018) is likely a systematic uncertainty due to subjective analysis choices (using flux vs. magnitude, functional form of the long term trend, functional form of the cycle variability, inclusion and form of a jitter term, etc.).An additional search for the star's activity cycle signature in Ca II H & K, Na I, and Hα activity indicators by Kumar & Fares (2023) yielded signals spanning 5.1-6.8yr.All of the aforementioned values fall near the predicted value of 7.6 yr based on the relationship between activity cycle and 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 a 2-5 orbit visit).Fitting a sinusoid to these data, fixed to the period of the optical cycle, yields amplitudes of 30-50%.The amplitude of cycle variability in the pseudocontinuum stands out at 89 +174 −13 %.We interpret this difference with caution because, on a pixel level, flux in the pseudocontinuum was below COS's noise floor 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, optical variability lags the FUV variability by 119 ± 6 • (Figure 6).dS19 found that variations in the Ca II S-index were out of phase with optical variations, though the constraint was not quantified.
The intermediate FUV-optical phase offset could indicate that GJ 436 is undergoing a transition from spotdominated to faculae-dominated activity.Sun like stars exhibit similar intermediate phase offsets between their activity cycles in Ca II H & K and optical emission as they transition from spot to faculae dominance (Reinhold et al. 2019) around a Rossby number of one, 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 −1 to the FUV emission would suffice, as would allowing the FUV cycle to vary in amplitude such that it is about half the value in 2018 as in 2012.Another option is to interpret the 2012 and 2015 observations as sampling a rotational minima, meaning average fluxes were higher at those epochs, which would shift the cycle curve to the 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 activity cycles at optical and EUV emission wavelengths (Fröhlich 2006;Woods et al. 2022;Llama, private comm.).

Excess White Noise
We included an excess white noise parameter ("jitter" term) in our fits so that the MCMC sampler could ex- plore an alternative to forcing 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.
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 (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α (Loyd et al. 2017;Lavie et al. 2017;dS19).Integrated to a 1 h cadence (and assuming a flat power spectral density), the noise level is 5.8 ± 0.8% in emission from the least-ionized ion 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 h cadence.

Erratum: Corrected magnetic field limit for GJ 436 b
This section reflects an erratum relating to the subsequent section submitted on 2024 February 23 and has been copied here for the convenience of readers of the arXiv version.
In the original manuscript, it was reported that the upper limit on variations in far ultraviolet line emission from GJ 436 phased with the orbit of planet b implied an upper limit of ∼100 G.This value is incorrect due to an arithmetic error discovered in the analysis scripts following publication.The revised limit is ∼10 G.In contrast with the previous value, the revised limit is comparable to the magnetic field strength of Jupiter and lower than the magnetic fields reported for 4 hot Jupiter planets using the same methodology with Ca II H & K observations (Cauley et al. 2019).
Later work by Vidotto et al. (2023) has shown that an assumption critical to this estimate is invalid.The estimate is based upon the model of Lanza (2013), which assumes the planet orbits in a potential magnetic field.The Vidotto et al. (2023) model demonstrated that GJ 436 b orbits beyond the range where the stellar field is approximately potential.Although invalid for GJ 436 b specifically, our field constraint remains generally useful for estimating the sensitivity of FUV observations of magnetic SPI to exoplanetary magnetic field strengths.Note that the hypothesis that the unusual flare behavior of GJ 436 might be a result of a magnetic star-planet interaction remains viable, as the Vidotto et al. (2023) model confirmed that GJ 436 b primarily orbits within the Alfvén surface of the star.

Magnetic Star-Planet Interaction
To search for a star-planet interaction, 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 variation 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 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.
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 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 empiricallyscaled estimate of the stellar wind speed at the planet of 235 km s −1 .In contrast, Bourrier et al. (2016) estimate a lower value of 85 km s −1 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 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 star and 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 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 semi-major 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  We speculate that the lack of flares when the planet is near eclipse (phase 0.5) could suggest a magnetic SPI is triggering the star's frequent small flares, but with energies too weak to register as a detection using the truncated sine model.
AD Leo reported by Hawley et al. (2003).These indicate that roughly 3% of the total energy 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 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 d 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 d; Gorrini et al. 2022) and a similar level of FUV emission, yet 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 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 build up of stored energy necessary to give rise to the larger flares exhibited by other M dwarfs.
The comparison sample from L18, which produced stronger flares, are less prone to a magnetic SPI.Within the sample, 8/10 (6/6 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 two of mass or distance.The closest comparison is GJ 581 b, a planet with M sin i = 15.8 ± 0.3 M ⊕ orbiting 0.040 AU from its host (Robertson et al. 2014) versus GJ 436 b's 25.4 +2.1 −2.0 M ⊕ 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 et al. 2012;Kossakowski et al. 2022).
Consistent with an SPI hypothesis is the lack of flares during planetary eclipses, when a sub-planetary 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 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 flares included, the MCMC sampler does favor a slight increase in 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 time-averaged emission with planet phase, and FFD statistics could all lead to confirming or refuting an SPI in this system.The present data hint that flare statistics might have the greatest quantitative power in identifying a magnetic SPI.

Variability Synopsis
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." 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 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 not caught within the limited duration of the observations.
Excess noise, since we assume 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 cycles, 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 amplitude of other variability sources will wax and wane over the course of the stellar activity cycle.
Variability in most of the FUV bands we considered is well correlated.Figure 9 plots the correlation 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.Variations in these bands, particularly at the high end (mostly due to flares), appear subdued relative to other lines, in line with 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.
We expect the extreme UV (EUV) emission of GJ 436 to track the changes in 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 lines fluxes increased by about 30%, while the flux in a 100-920 Å band predicted by the DEMs increased by about half 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 variations in the EUV should exceed those in the FUV for GJ 436.France et al. (2018)  rotation, activity cycles, catch-all "excess noise," and a possible star-planet interaction (SPI), with comparison to the optical variability measured by Lothringer et al. (2018).We focused on the summed emission from the C II, C III, Si III, Si IV, and N V lines in the 1150-1450 Å range, the strongest lines in this range aside from Lyα and O I, which were contaminated by geocoronal airglow.
Figure 8 provides a rapid comparison of the contribution of each variability source to the star's overall variability.The star's 2012-2018 activity cycle was responsible for the greatest degree of observed variability, 38 ± 3% in amplitude.Observed flares increased the time-averaged emission by 15%.Flares not observed, either because they were too weak to register above the noise or too rare to be captured in the limited observing time, would add to this value.Stellar rotation produced variations of 8.1 +2.1 −1.6 % and a mixture of instrumental and astrophysical sources treated with a catch-all white noise term account for an additional 10.8 +1.5  −1.1 % variability in excess of photon noise.We did not detect variability phased with the planet's orbit, with a 2σ upper limit of 9.2%.
Uncertainty in the optical rotational phase at the time of the FUV observations precludes strong conclusions about the relative distribution of spots, faculae, and plage.Some other factors suggest the star could be undergoing a transition from spot-to faculae-dominated activity, including the lack of reddening during optical minima found by Lothringer et al. (2018), the sinusoidal form of the star's optical activity cycle, a possible 119 ± 6 • phase offset between the optical and FUV cycles, and the star's Rossby number near unity.
If the planet is magnetically interacting with its star, its magnetic field must have a strength ≲ 100 G. Our upper limit on the SPI signal assumes variations appear as a truncated sine function, such as would be produced by a planet-induced hot spot traversing the stellar surface.However, the star produced an unusual population of flares that we conjecture could be the result of a magnetic SPI.Although the flares were numerous, totaling 12 (plus 2 dubious) events over the course of about a day of cumulative exposure, their energies were anomalously low in comparison to other M dwarfs.We speculate that this could result from a magnetic SPI triggering a multitude of flares with low equivalent durations (a metric of flare energy normalized by the star's quiescent flux).This would increase the flare rate at low equivalent durations while disrupting the energy buildup necessary for stronger events.The lack of flares during planetary eclipse support this view, but these visits were too brief to have statistical power.Software: emcee (Foreman-Mackey et al. 2013), cal-COS (https://github.com/spacetelescope/calcos),As-tropy (Astropy Collaboration et al. 2013, 2018, 2022), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020)

Figure 1 .
Figure1.Coadded spectrum of GJ 436.The range of the plot matches that of the F U V130 band.The gray regions delineate regions of airglow contamination and a detector gap that falls near 1300 Å that we masked from the F U V130 band.Labels indicate strong emission lines that we also analyzed for variability, and the regions between the labeled lines and gray bands represent the pseudocontinuum.

Figure 3 .
Figure 3.Our fit to the Strömgren (b + y)/2 photometry from Lothringer et al. (2018), mimicking their Figure 1.Left: Data with the sinusoidal fit to the 7.75 ± 0.10 yr activity cycle period and the linear trend overlaid.Line thickness represents uncertainty in the fit.Right: Data with the activity cycle sinusoid and linear trend subtracted, then folded onto the best fit rotation period of 44.09 ± 1.16 d, and binned by groups of 30.The rotation model sinusoid is overlaid with the translucent region representing uncertainty in the fit.Error bars are the sample error on the mean of the binned points.

Figure 4 .
Figure4.Cumulative frequency distribution of GJ 436 flares in the F U V130 band.Gray regions give the uncertainty envelope from fits to a compendium of M dwarf flares in L18.To highlight the cliff in the GJ 436 flares, the figure also shows random draws of even equivalent durations and energies assuming the same number of events above 500 s and 2.5 × 10 28 erg and a inverse power-law slope, α, equal to the median observed by L18.The gray horizontal line represents the limit where one event would be expected to occur within the cumulative observing time.

Figure 5 .
Figure 5. Fit to the rotational modulation of flux in each of the analyzed bands during the 2017/2018 epoch of observations.The data have been cleaned of flares.The fits were made to data binned to a 10 min cadence, but to minimize clutter the plotted data are binned to full exposure length (typically around 45 min).Fits corresponding to alternative periods are shown as dotted and dashed lines in the panel second from top (see text).The slight downward trend in the UV fits is due to the simultaneous fit to the activity cycle.

Figure 6 .
Figure 6.Summed line flux (blue points), binned by HST visit, across all three epochs of HST data, fit with a sinusoidal model for activity cycle variations (orange curve) fixed to the optical cycle period, and shown alongside the (rescaled) fit to optical cycles (green curve) to compare phases.

Figure 7 .
Figure 7. Fluxes plotted against the planetary orbital phase.The planet's optical transit spans a narrow range of ±0.008 in phase.Top: Data from the 2017/2018 epoch binned by exposure with the best-fit rotational modulation subtracted.Blue points exclude and red points include flares.The black solid line shows the 2σ upper limit amplitude for a truncated sine SPI model fit to the data with flares included.Bottom: Data at a 10 s cadence to show individual flares.Blue (quiescent) and red (flare) data are from the 2017/2018 epoch whereas green (quiescent) and orange (flare) data are from the earlier 2012 and 2015 epochs.We speculate that the lack of flares when the planet is near eclipse (phase 0.5) could suggest a magnetic SPI is triggering the star's frequent small flares, but with energies too weak to register as a detection using the truncated sine model.

Figure 9 .
Figure9.Correlations between the variability in each of the analyzed bands.The F U V130 and Summed Line bands are omitted since they are a compendium of all the other bands (thus not independent).Gray lines indicate a 1:1 correlation.The p-value from a Spearman rank-order correlation test is given in each panel.Most band pairs are correlated to very high confidence, with the exception of some correlations with the pseudocontinuum.The largest excursions from the means are driven by flares, indicated by red coloring of the points.Notably, neither N V nor the pseudocontinuum seem to follow a 1:1 relationship with other bands, although the statistical significance of the relationship remains very high despite its shallower slope.

Table 2 .
Fits to variability in FUV emission lines.
also noted an effectively linear correlation between Si IV and N V emission measured by HST and 90-360 Å flux measured (at a different epoch) by the Extreme UltraViolet Explorer (EUVE) 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.Summary of the contributions from various sources to the FUV variability of GJ 436.Arrows indicate upper or lower limits (see text).In parenthesis below each label are the timescale at which we measured each form of variability.