Unsigned magnetic flux as a proxy for radial-velocity variations in Sun-like stars

We estimate disc-averaged RV variations of the Sun over the last magnetic cycle, from the single Fe I line observed by SDO/HMI, using a physical model for rotationally modulated magnetic activity that was previously validated against HARPS-N solar observations. We estimate the disc-averaged, unsigned magnetic flux and show that a simple linear fit to it reduces the RMS of RV variations by 62%, i.e. a factor of 2.6. We additionally apply the FF' method, which predicts RV variations based on a star's photometric variations. At cycle maximum, we find that additional physical processes must be at play beyond suppression of convective blueshift and velocity imablances resulting from brightness inhomogeneities, in agreement with recent studies of solar RV variations. By modelling RV variations over the magnetic cycle using a linear fit to the unsigned magnetic flux, we recover injected planets at an orbital period of about 300 days with RV semi-amplitudes down to 0.3 m/s. To reach semi-amplitudes of 0.1 m/s, we will need to identify and model additional physical phenomena that are not well traced by the unsigned magnetic flux or FF'. The unsigned magnetic flux is an excellent proxy for rotationally modulated, activity-induced RV variations, and could become a key tool in confirming and characterising Earth analogs orbiting Sun-like stars. The present study motivates ongoing and future efforts to develop observation and analysis techniques to measure the unsigned magnetic flux at high precision in slowly rotating, relatively inactive stars like the Sun.


INTRODUCTION
The main obstacle we face in detecting, confirming and characterising Neptune-to Earth-mass exoplanets via radial-velocity (RV) monitoring is the intrinsic mass dictates the amount of observing time required to characterise a planet's atmosphere, so it is essential that we know masses reliably to plan atmospheric follow-up observations (Morley et al. 2017;Batalha et al. 2019), e.g. with JWST and ARIEL. To determine accurate and precise planetary masses, we need to develop robust, physically motivated models for stellar variability. We still lack a complete and detailed understanding of how the interplay between magnetic fields and granulation gives rise to RV variations on the Sun and other stars (Blackwood et al. 2020).
The Sun is the only star whose surface we can image directly and at high resolution, making it an ideal test bench to examine the physical phenomena responsible for instrinsic RV variability. It is also the only star whose RV we know independently of spectroscopic measurements (e.g. from HARPS-N).
On timescales of several rotation periods (weeksmonths), RV variability is driven by magnetic activity in the photosphere. The manifestations of magnetic activity relevant to the present analysis are sunspots and faculae. Sunspots are relatively large, dark areas of strong magnetic fields (e.g. Foukal 2004, Chap. 8). Faculae are small, bright magnetic flux tubes (Spruit 1976). They tend to be located in the lanes between supergranular cells and are spread all over the solar surface, thus forming the photospheric magnetic network (e.g. Foukal 2004, Chap. 5, p. 145;Chap. 8). In regions of enhanced magnetic activity, faculae cluster into areas of plage (e.g. Schrijver & Zwaan 2000, Chap. 1, Fig. 1.1). The total surface area covered by network varies throughout the Sun's 11-year cycle, although to a much smaller extent than spots and plage (Meunier 2018(Meunier , 2003. At low activity levels, magnetic elements tend to be spread throughout the solar surface in the form of network. Plage coverage increases with magnetic activity, and at high activity levels, the majority of magnetic elements on the surface are concentrated in plage, rather than in network (e.g. see plage and network filling factors shown in Figure 3). Plage generally decays into network over over timescales of several rotations, which explains the larger network coverage when activity is high.
Magnetic elements inhibit convective motions, thereby suppressing some of the convective blueshift that results from granulation (e.g. Dravins et al. 1981). This suppression of convective blueshift is the dominant contributor to RV variations in the Sun (Saar & Donahue 1997). Solar observations show RV variations, modulated by the Sun's rotation and evolving over timescales of days to weeks, with amplitudes of several m s −1 (Meunier et al. 2010a,b;Dumusque et al. 2014;Haywood et al. 2016;Milbourne et al. 2019). Convective blueshift is suppressed by faculae in concentrated areas of plage, while faculae in the network, being more spatially diffuse, do not perturb convective flows significantly (Milbourne et al. 2019, Sect. 4.4). Sunspots contribute little to observed suppression of convective blueshift, as they are dark (and thus contribute relatively less to observed spectra) and cover very little area in comparison to faculae (e.g. Lagrange et al. 2010;Haywood et al. 2016). At high solar activity levels, sunspots do contribute significantly to RV variations. Because they are much darker than the surrounding photosphere, they produce significant inhomogeneities in surface brightness, which result in RV variations with an rms of 60 cm s −1 and frequent peak-to-peak amplitudes of 2 m s −1 , with variations of up to 5 m s −1 (Lagrange et al. 2010). On the other hand, faculae (both in network and plage) are only 10% brighter than the quiet photosphere at optical wavelengths; moreover, they are distributed much more uniformly longitudinally on the solar disk, so their brightness-induced RV contribution mostly cancels out on rotational timescales. On longer (magnetic cycle) timescales, large-scale changes in the number of faculae can produce long term, bulk RV shifts (e.g. Saar & Fischer 2000;Meunier et al. 2010b, Fig.8).
Magnetic elements enhance chromospheric column density, which strengthens the emission reversals in the Ca ii H&K cores. Thus, the filling factor of magnetic elements correlates strongly with the amount of emission in the cores of the Ca ii H&K lines (e.g. Meunier 2018, Fig.1) as measured by the log R HK index (Vaughan et al. 1978;Noyes et al. 1984). Solar observations show that RV variations and log R HK correlate strongly over long timescales of several years, i.e. over the Sun's 11year magnetic cycle (e.g. Meunier et al. 2010a, Fig.13). However, when we look at shorter timescales of a few weeks to months, i.e. on the solar rotation timescale, the log R HK does not systematically trace RV variations down to sub m s −1 precision. This is the case during both low and high activity phases. At low activity levels (on the rotation timescale), the low correlation between log R HK and RV variations can be explained by the fact that magnetic elements are predominantly found in network, which contributes to Ca ii H&K emission but does not affect RV variations . At high activity levels (on the rotation timescale), RV variations do not correlate well with Ca ii H&K emission, e.g. in Lagrange et al. (2010, Fig.12) and Haywood et al. (2016, Fig.10). Part of this discrepancy may be that the Ca ii H&K emission forms in the chromosphere, and therefore follows a different limb-darkening law and projection effects than RVs, which are measured from photo-spheric absorption lines (see Section A in this paper). To summarise, solar RV variations induced by magnetic activity on timescales of the order of a few rotation periods are not directly correlated and in phase with Ca ii H&K emission, in part because faculae affect RV variations differently depending on whether they are in sparse network or in concentrated regions of plage, and due to changes in spectral line profiles.
We see a similar behaviour in observations of Sunlike stars. As part of the Mt. Wilson HK Project, several dozen slowly rotating, old Sun-like stars were monitored in optical photometry and in Ca II emission over the past several decades (Baliunas et al. 1995;Wilson 1978). Observations showed that as these stars became more magnetically active (as indicated via their Sindex), they also got brighter. Their surfaces are therefore dominated by bright faculae rather than dark spots, just like the Sun (Lockwood et al. 2007;Radick et al. 2018). On timescales of several years, the rms of stellar RV variations increases as the log R HK increases (Saar & Fischer 2000;Lovis et al. 2011). Aigrain et al. (2012) developed a model to estimate stellar activity-induced RV variations based on the star's optical photometric variations. The model accounts for RV variations produced by dark spots and bright plage that are spatially associated with spots, both through rotational flux imbalance and suppression of convective blueshift. Haywood et al. (2014) tested this model on the solar analog CoRoT-7 using simultaneous photometric and RV observations taken at high cadence over a rotation period. While the model of Aigrain et al. (2012) captures a significant part of the rotationally modulated RV signal, it leaves out an equally significant rotationally modulated signal. This additional RV variation likely originates from magnetic regions that have low intensity contrast, and whose brightness-induced RV variations is therefore low. These observations and their interpretation are consistent with the Sun's behaviour (the solar surface is dominated by low-contrast plage).
To confirm and characterise long-period, low-mass exoplanets, we need a proxy that traces RV variations systematically and at sub-m s −1 precision. On both the Sun and other Sun-like stars, Ca ii H&K emission does not systematically correlate as strongly with activityinduced RV variations. For the Sun, a strong correlation has been observed between activity-induced RV variations and the unsigned, full-disc magnetic flux, over the magnetic cycle (Deming & Plymate 1994;Lanza et al. 2016;Meunier 2018) and on the rotation timescale . The disc-averaged RV timeseries estimated from the Michelson Doppler Imager onboard the Solar and Heliospheric Observatory (SoHo/MDI) by Me-unier et al. (2010b) and subsequent papers by Meunier et al. is well-sampled and spans over 4 years of Cycle 23, prior to SDO's launch; however, it cannot be compared to "ground-truth", direct disc-integrated RV observations. To date, systematic RV campaigns of the Sun as a star have been carried out spectroscopically, with HARPS via sunlight reflected from asteroids Lanza et al. 2016) and more recently with dedicated solar feeds at HARPS-N (Dumusque et al. 2015;Collier Cameron et al. 2019) and HARPS (Dumusque 2019, HELIOS).
In the present analysis, we estimate rotationally modulated RV variations of the Sun from SDO/HMI images over 8 years, at daily cadence (Section 2.1), using a technique that has been validated against direct Sun-asa-star HARPS-N observations (Section 3). We present timeseries of RV variations and unsigned magnetic flux in Section 3. We model our RV timeseries using a linear fit in unsigned magnetic flux in Section 4 and with the FF' model of Aigrain et al. (2012) in Section 5. In Section 6, we identify potential additional physical effects giving rise to rotationally modulated RV variations that are not well traced by either models, and identify limitations of the FF' model. We perform simple planet injections to assess the performance of the unsigned magnetic flux for mitigating rotationally modulated RV variations in Section 7. We discuss future prospects, including ways to measure the unsigned magnetic flux in other stars in Section 8, and present our conclusions in Section 9.

SDO/HMI images
We use 720-second HMI exposures of continuum intensity (both uncorrected and corrected for limb darkening by the HMI team), Dopplergrams, and magnetograms, as represented in Figure 1. The HMI instrument takes 6 measurements of intensity across a narrow wavelength range centered on the Fe I line at 6173Å (see Fig. 6 of Schou et al. 2012). These points are fitted with a Gaussian profile to generate the main HMI data products, to generate the main HMI data products including velocity (line shift), intensity (depth), magnetic field strength (width due to Zeeman broadening) and continuum intensity (Schou et al. 2012, Sect.3.3). The line shifts and magnetic field values extracted for each pixel should be independent, physical quantities, obtained from combinations of the different intensities at different points on the measured Fe I line. While the line asymmetry stemming from convection within each pixel is not preserved, we expect this technique to capture asymmetries due to physical processes occuring over scales larger than a Figure 1. SDO/HMI data products for an observation set taken on 2015 November 28 at 20:00:00 UTC. This set of images is representative of the Sun during high activity levels. Faculae and sunspots fill 3.25% and 0.03% of the solar disc, respectively. The areas that suppress convective blueshift, predominantly faculae in regions of plage, shown in the last panel, fill 1.59% of the solar disc. The disc-averaged unsigned magnetic flux is 9.99 G. Notes: The flattened intensity is normalised to the mean intensity. The line-of sight velocity is shown after subtracting the solar rotation profile and the velocity of the spacecraft. pixel. We refer the reader to Schou et al. (2012) for further details on how these images are extracted from the raw filtergrams.
2.1.1. Temporal sampling SDO/HMI has operated almost continuously since the start of the mission, except for spacecraft operations and calibrations, and eclipses that happen due to the geosynchronous orbit of the SDO spacecraft (Hoeksema et al. 2018). There have been very few anomalies requiring interruptions, none of which have been prolonged compared to the seasonal eclipses. We take a set of images every 4 hours (6 times per 24-hour period) from 2010 April 07, 04:00:00 UTC up to 2018 January 12 20:00:00 UTC, amounting to a total of 16855 sets of images spanning 2811 days. We take daily averages to minimise the contribution of short-term processes, namely oscillations and magnetoconvection. We choose not to use the SDO/HMI images at their highest cadence in order to maintain the relevance of this analysis to stellar studies, while still sampling the Sun multiple times a day. Indeed, the cadence of current and planned stellar RV surveys is 1-3 observations per night at most.

Instrument precision and stability
All sources of uncertainty that affect our RV estimates are listed in Table 1. Hoeksema et al. (2018) recently assessed the performance of HMI and reported that the instrument continues to work according to its original specification. The data products are corrected on a regular basis as the calibrations improve (e.g. instrument thermal environment, focus, image distortions, optics alignment, cosmic rays correction, etc.). They report that the quality of the data is very uniform with time. The HMI data products are well calibrated (Hoeksema et al. 2018), with the exception of the long-term stability of the Dopplergrams. The Doppler velocity maps were designed for helioseimology investigations, so they are not calibrated to be stable over timescales longer than a few hours or days (Schou et al. 2011). This is shown in Figure 2, which shows the disc-averaged velocity of the Sun taking irregular jumps of several m s −1 over the course of the SDO mission. To correct for this effect, we perform all of our velocity calculations relative to the disc-averaged velocity of the magnetically inactive, quiet Sun for each Doppler image (as in Meunier et al. 2010b;Haywood et al. 2016;Milbourne et al. 2019). We estimate the velocity of the quiet Sun by summing over all pixels identified as non-magnetic (see Section 3), and excluding pixels considered to magnetically active (even at the peak of the Sun's magnetic cycle, fewer than 5% of all pixels within the solar disc are magnetically active). We are therefore only considering RV variations (∆RV ). Importantly, this subtraction cancels out all velocity flows from the quiet Sun. This means that our RV estimates are free from pressure-mode (p-mode) oscillations, granulation and supergranulation motions, which would otherwise induce uncorrelated noise at the 1 m s −1 level (e.g. Meunier et al. 2015).
Beyond the lack of long-term calibration, Couvidat et al. (2016) report that the most significant source of instrument-related uncertainty that remains in individual HMI Dopplergrams is the orbital velocity of the spacecraft, that is uncertain to 0.01 m s −1 . Indeed, we see a systematic sinusoidal shift with a periodicity of 12 and 24 hours due to the spacecraft's orbit. This systematic should mostly average out through our sampling, but to be on the safe side, we add Couvidat et al. (2016)'s uncertainty of 0.01 m s −1 in quadrature to our RV uncertainties (see Table 1). The total number of pixels inside the solar disc as it appears on the HMI image varies by ±3% over the course of each year, primarily due to Earth's eccentric orbit. To correct for this effect, we normalise all our quantities by the total number of pixels in each image, i.e. we estimate discaveraged quantities. We ran tests to assess the potential effect of hypothetical spurious pixels from the continuum, Doppler and magnetic images. For example, we find that setting the velocity to zero for a large patch of 20 × 20 pixels square incurs changes in RV of 0.04 m s −1 . The RV effect of spurious pixels distributed randomly on the solar disc is less than 0.01 m s −1 . To be conservative, we assign a constant RV uncertainty of 0.1 m s −1 for all RV estimates, to account for uncertainties arising from our pipeline to estimate RVs (detailed in Section 3). Our choice of magnetic and continuum intensity thresholds are based on previous studies, but they remain somewhat arbitrary, and changing them slightly will impact RV estimates at the cm s −1 level. We conservatively remove sets of images with v sin i values more than 3sigma deviant from the mean (6 out of 16855), and those with focal length values more than 3-sigma deviant (2 out of 16855).

SORCE total solar irradiance observations
To apply the FF' method of Aigrain et al. (2012) in Section 5, we use total solar irradiance observations (TSI; akin to a Kepler lightcurve for the Sun) taken by the Total Irradiance Monitor (TIM) onboard the SOlar Radiation and Climate Experiment (SORCE) satellite (Kopp et al. 2001;Lean et al. 2005)   . Disc-averaged radial velocity of the Sun as estimated from SDO/HMI images, before subtracting the discaveraged velocity of the quiet Sun. The jumps in RV are instrument systematics, as HMI is not calibrated for longterm stability. duration of our SDO/HMI timeseries. The TIM takes observations every 50 seconds when the spacecraft faces the Sun, and these observations are then combined to produce daily averages. We concatenate this timeseries with our SDO/HMI timeseries of daily images. Because there is a gap in the SORCE timeseries around 1200-1400 days into the SDO mission, we are left with a combined timeseries of 2535 daily observations, spanning 2811 days. The SORCE timeseries is plotted in the last panel of Figure 3. The TIM achieves a precision of 4-17 ppm per observation (Kopp 2014). For comparison, Kepler achieved the same level of precision on a 7 to 9th V magnitude star in a long-cadence observation 2 (30-minute integration time). TIM has a long-term stability of about 10 ppm per year (Kopp 2014).

S-index from Mt Wilson and HARPS-N
We compare our SDO/HMI-derived quantities (∆RV , |B obs |) against Ca ii H&K emission observations. For this, we use overlapping S-index observations of the Sun seen as a star at the Mount Wilson Observatory as part of the HK Project, fully homogenised and calibrated by Egeland et al. (2017). Their observations run from 1966 until 2015. To cover the 2015-2018 period, we use daily averaged S-index observations taken by the solar telescope that feeds the HARPS-N spectrograph since July 2015 ). The Mt Wilson and HARPS-N datasets overlap for around 45 days in 2015, which we can use to stitch the two S-index timeseries together. We do this by rescaling the overlapping part of the HARPS-N S-index timeseries so it has the same variance as the Mt Wilson S-index timeseries and subtracting the offset between the two datasets. The full timeseries is shown in the second to last panel of

ESTIMATING THE FULL-DISC RV VARIATIONS AND MAGNETIC FLUX OF THE SUN
We estimate disc-averaged active-region filling factors, RV variations and unsigned (unpolarised) magnetic fluxes of the Sun using spatially resolved images from SDO/HMI images according to the same method as Milbourne et al. (2019), adapted from that of Haywood et al. (2016), which builds on the techniques originally developed by Meunier et al. (2010b) and Fligge et al. (2000).

Separating magnetically active regions from quiet Sun
We separate magnetically active regions from quiet-Sun regions by applying a threshold in unsigned radial magnetic field strength for each pixel according to the cutoff found by Yeo et al. (2013): The factor µ accounts for foreshortening, and is equal to cos θ ij , where θ ij is the angle between the outward normal to the feature on the solar surface and the direction of the line-of-sight of the SDO spacecraft. The term σ B obs,ij represents the noise in the observed magnetic field, for each pixel at position i, j on the image. Yeo et al. (2013) estimated σ B obs,ij to be 8 G (photondominated), so the magnetic field threshold is 24 G. We exclude isolated pixels that are above this threshold as they are likely to be false positives.

Filling factors of sunspots & plage
To identify faculae and sunspots, we apply the intensity threshold of Yeo et al. (2013), at 0.89 times the mean flattened intensity over quiet-Sun regions. We further identify faculae in concentrated regions of plage, as opposed to faculae dispersed in the network (cf. Introduction, § 3), according to the area threshold estimated by Milbourne et al. (2019). We identify plage as magnetically active facular regions whose area on the flattened solar disc exceeds 20 microhemispheres (µhem), corresponding to about 60 Mm 2 .
We estimate the disc-averaged filling factors of sunspots and plage as follows: where N pix is the total number of pixels in the solar disc and the weight W ij is set to 1 in sunspot (or plage) pixels, and 0 in quiet-Sun pixels.

RV variations
Milbourne et al. (2019) derived solar RV variations from SDO/HMI images for an 800-day period overlapping disc-integrated RV observations of the Sun with the HARPS-N spectrograph. They reproduced rotation-modulated RV variations in good agreement with the HARPS-N observations, down to an rms level of 1.21 m s −1 , which is consistent with residual motions that are expected from granulation and supergranulation (Meunier et al. 2015). Their model, which we apply here, accounts for the suppression of convective blueshift from magnetic regions, and the velocity imbalances resulting from brightness inhomogeneities. We refer the reader to the Appendix of Milbourne et al. (2019), which fully describes the model. Estimating RV variations according to this technique and model allows us to determine solar RV variations that we can compare directly with spectroscopic measurements of other stars, which are derived from thousands of spectral lines, not just the Fe i line measured by SDO/HMI.

Unsigned magnetic flux |B obs |
We compute the disc-averaged, line-of-sight unsigned (i.e. unpolarised) magnetic flux of the Sun, by summing the intensity-weighted line-of-sight absolute magnetic flux in each pixel according to Haywood et al. (2016): where I ij is the observed, non-flattened continuuum intensity of the Sun. We do not flatten the intensity continuum in order to obtain the observed unsigned magnetic flux.

Correlations between ∆RV , |B obs |, photometry and S-index
The timeseries of RV variations, unsigned magnetic flux, plage and sunspot filling factors are plotted alongside coeval TSI and S-index observations in Figure 3. The SDO/HMI disc-averaged quantities shown in Figure 3 are at the maximum cadence considered in this study (1 observation every 4 hours, i.e. 6 per day). We then average the SDO/HMI quantities over daily bins and concatenate these timeseries with the timeseries of S-index and TSI. We show the RV variations as a function of unsigned magnetic flux, S-index and TSI in the top row of Figure 4. The bottom row shows the three activity indicators plotted as a function of each other. Figure 4 shows that the RV variations correlate much better with the unsigned magnetic flux (R = 0.92) than the S-index (R = 0.75) or optical photometry (R = 0.46). Observations of high sunspot coverage are highlighted in yellow. We see that the Sun is dominated by spots at high activity levels. At the peak of the activity cycle, the Sun's photometric variations are anti-correlated with Ca ii H&K variations (Radick et al. 2018, Fig.1), so the Sun is spot-dominated; at lower activity levels, they are positively correlated, implying that the solar surface is dominated by faculae/plage. This is in agreement with previous studies (Fröhlich & Lean 1998;Krivova et al. 2007;Shapiro et al. 2014).
Following  who reported on a hysteresis between ∆RV and Ca ii H&K in solar Cycle 23, we investigate the hystereses between ∆RV , S-index, and |B obs |in Cycle 24. We observe hystereses between all quantities (shown in Figure A1) and discuss their physical origins in detail in Appendix A.
Time lags of 1-3 days between RV variations and the bisector span and FWHM have been reported previously in spectroscopic HARPS-N observations of the Sun (Collier Cameron et al. 2019, Fig.15). We cross-correlate ∆RV , |B obs |, and the S-index against each other to look for time shifts between them. We do not find any significant time shifts between any of our observables.
3.6. Periodogram analysis Figure 5 shows Generalised Lomb-Scargle periodograms (Zechmeister & Kürster 2009) of the RV variations (panel (a)), the filling factors of plage (panel (b)) and sunspots (panel (c)) and the unsigned magnetic flux (panel (d)). Most of the periodicity below 100 days is confined to periods close to the rotation period and its first harmonic. It is worth noting that these peaks are in fact forests of peaks, in which several peaks are significant above the 0.001% confidence level. This means that depending on when or for how long we might observe the Sun, we may measure rotation periods differing by several days (e.g. Mortier & Collier Cameron 2017;Nava et al. 2019). In the periodogram of the sunspot filling factor, we also detect a significant peak consistent with the 20.8 day peak detected by Lagrange et al. (2010). This peak is possibly related to the lifetime of the spots. Alternatively, it may be associated with global-scale equatorial Rossby waves (rmodes) that produce oscillations on a 19-day recurrence timescale (Lanza et al. 2019). In this periodogram, we see many significant peaks around the rotation period and its harmonics. Additionally, some peaks are significantly different from the rotation period or its harmonics; as Lagrange et al. (2010) previously emphasized, we should be careful when attributing these signals to nonactivity processes. The strong signals in the spot filling factor do not necessarily translate into signals in the total RV variations, because the Sun is faculae-dominated for the majority of its cycle (e.g. Fröhlich & Lean 1998;Krivova et al. 2007;Shapiro et al. 2014). However, in younger, faster rotating Sun-like stars whose behaviour has been observed to be spot dominated (see Lockwood et al. 2007;Radick et al. 2018), we would certainly expect the RV variations to show a more "spot-like" periodogram structure.
4. MODELLING ∆ RV USING |B obs | We model the RV variations estimated in Section 3, ∆RV (t) as a linear model of |B obs |: where α is a constant scaling factor and RV 0 is a constant zero-point offset. |B obs | is the mean of |B obs | over the full timeseries. We optimise the parameters α and RV 0 via a least squares procedure. We model the full ∆RV (t) timeseries of daily averages (plotted in the top panel of Figure 3). The fit over the full magnetic cycle is shown in the top panel of Figure 6. To examine the performance of the linear |B obs | model as a function of magnetic activity levels, we model two separate stretches, at activity maximum and minimum, spanning 600 days each. The low-activity, quiet epoch ranges be-  Table 2. The root mean scatter (rms) of the full dataset is 2.33 m s −1 , and that of the residuals is 0.89 m s −1 . Overall, a simple |B obs | model reduces the rms of the RV variations by 62%, i.e. a factor of 2.6. Although we do see correlated residuals at times of high activity, the residuals over the full cycle are flat.
Rotationally modulated RV variations at times of low activity -The dominant process at play is the suppression of convective blueshift incurred by areas of plage . |B obs | correlates well with their presence, as seen in Figure B2 in Appendix B. There are very few spots; the maximum filling factor of sunspots in this 600-day stretch is 0.03% (compared with 0.14% in the active stretch and 0.09% overall). We therefore expect |B obs | to correlate well with RV variations. Indeed, this model is an excellent fit during activity minimum, as evidenced in Figure 6b. Over the quiet epoch, the RV residuals have an rms of 0.56 m s −1 . Rotationally modulated RV variations at times of high activity -The best-fit estimates of the model parameters (α and RV 0 ) differ significantly from the fit at low magnetic activity (see Table 2). This is because convective blueshift is more suppressed by larger magnetic structures (which are more prevalent in periods of high activity), i.e. the faculae in concentrated areas of plage, as previously found by Meunier et al. (2010b, Fig .7) and Milbourne et al. (2019, Fig.6). The |B obs | model accounts for RV variations down to a residual rms of 1.02 m s −1 . As seen in Figure 6 (a), significant rota-tionally modulated RV variations remain unaccounted for. In Figure 7(a), we zoom-in further on the RV variations over 4 to 5 solar rotations at activity maximum. This stretch spans 120 days from 2014 October 6 (JD = 2456937) to 2015 February 3 (JD = 2457057; days 1620-1740 in the figures). The smooth, rotationally modulated signal that remains in the RV residuals is very clear on this timescale. At activity maximum, we expect suppression of convective blueshift to produce significant RV variations (Meunier et al. 2010a,b;Haywood et al. 2016). Additionally, there are sunspots, which now  Haywood et al. 2016). The relationship between this sunspot flux-blocking RV term and |B obs | is more complex than for the RV due to suppression of convective blueshift. In fact, they do not correlate with each other (R = 0.08). When a spot crosses the central meridian, this RV contribution is zero, while |B obs | would be at its maximum. Therefore, a simple, linear |B obs | model cannot adequately capture RV variations from sunspot flux-blocking.
To investigate this hypothesis, we apply the technique developed by Aigrain et al. (2012). Their FF' term accounts for RV variations incurred by brightness inhomogeneities on a rotating disk.

MODELLING ∆RV USING BOTH |B obs | AND THE FF' MODEL
A simple |B obs | model fits RV variations well overall, but becomes insufficient at the peak of the solar magnetic cycle (see Section 4), where significant, correlated residuals remain, as visible in Figure 6 (a). We apply the method of Aigrain et al. (2012) in an attempt to account for RV variations from flux-blocking from sunspots (and plage).
As derived in Aigrain et al. (2012, Eqn.10), the RV perturbation due to a spot crossing the disc can be expressed as follows: where Ψ(t) is the observed solar flux, Ψ 0 is the solar flux for a non-spotted photosphere andΨ(t) is the first time derivative of Ψ(t). R is the solar radius. The parameter f spot represents the drop in flux produced by a spot at the centre of the solar disc, and corresponds to the sunspot filling factor. Since both dark inhomogeneities (from spots) and bright inhomogeneities (from faculae, mainly in plage) produce RV perturbations (e.g. Meunier et al. 2010a, Fig.7), we use the total magnetic filling factor f all rather than f spot alone. We find that when using f spot , the FF' term has an rms amplitude of 0.02 m s −1 , i.e. 10 times less than when we consider all magnetic elements, including plage. We write the following formulation: where F and F correspond to the TSI (Section 2.2) and its first time derivative, respectively. F and F are the means of the TSI lightcurve and its first time derivative, respectively. To compute F , we interpolate the TSI observations (F ) onto an evenly, oversampled array and then fit them using Gaussian-process regression using a basic square exponential kernel (Rasmussen & Williams 2006). We then compute the derivative using second-order accurate central differences. We multiply the FF' term above by a normalising factor f all F / F so that the full term is of order unity.
Our resulting |B obs | + FF' model is as follows: where β is a scaling factor that we fit for in our leastsquares optimisation, along with α and RV 0 . The resulting best fit when modelling the full timeseries (where the SDO/HMI and SORCE data overlap) is shown in the top panel of Figure 8. We also apply this model to the active and quiet epochs, as shown in panels (a) and (b) of Figure 8. The best-fit estimate of β is non-zero (see Table 2), and the FF' term has an rms of order 0.2 m s −1 . For the full timeseries, the rms of the residuals (0.85 m s −1 ) is slightly lower than that obtained with the |B obs | model (0.89 m s −1 ). However, the overall fits with |B obs | + FF' are very similar to those result- ing from modelling the RVs with |B obs | only (as is done in Section 4). The residuals still display correlated behaviour. Clearly, the FF' model does not fully account for the signals leftover from the |B obs | model. The residuals from both models tested (|B obs | and |B obs |+FF' ) display correlated behaviour during times of high magnetic activity. We propose two explanations.
First, the FF' term does not adequately fit the RV variations resulting from sunspot flux-blocking (Section 5). Several studies have previously found that the FF' method cannot fully account for RV variations (e.g. Oshagh et al. 2017;Bastien et al. 2014;Haywood et al. 2014). We note that the FF' term (or the F 2 term that can be used to account for suppression of convective blueshift) is not expected to match RV variations perfectly, because F includes the derivative of the limb  darkening of F , which should not be a part of the RV model. We demonstrate this in detail in Appendix C. Another possible explanation for the RV residuals is that there are additional processes at play, which are either missing from the RV model of Milbourne et al. (2019) and Haywood et al. (2016) used to estimate RV variations, or that do not correlate directly with |B obs | (or FF' ). This finding is consistent with that of Miklos et al. (2019), who investigated the activitysensitivity of spectral lines observed by HARPS-N and concluded that there must be additional factors, not yet accounted for by current state-of-the-art models such as Meunier et al. (2017). Other types of surface velocity fields not included in our model may give rise to rotationally modulated RV variations, such as: Evershed flows -Sunspots are made of umbral and penumbral regions. Evershed flows, which are contained within penumbral regions and flow radially outward from the central umbra to the outer edge of the sunspot, are tangential to the surface (Evershed 1909). They will be most visible in RV for sunspots located away from disc centre, where their flows are more directed along our line of sight (e.g. Haywood et al. 2016, Fig.6).
Moat flows -Outside of the penumbra, but within the active region, are so-called moat flows (e.g., Solanki (2003) and references therein). These are also tangential to the surface but weaker than penumbral flows (of order 1/4 to 1/2 as strong), but as they are brighter (since they are typically in plage) and normally cover an area ≈ 3× larger than the spot. They may contribute significantly to the total RV signal (e.g. Iampietro et al. 2019). The effect of moat and Evershed flows on the disc-averaged solar RV variations is being investigated by Saar et al. (in prep.).
Active region inflows - Gizon et al. (2001) report the presence of inflows towards active regions on the Sun's surface, with amplitudes up to 50 m s −1 (see Gizon et al. 2010, Sect.7.1.). These are also currently under study (Saar et al. in prep.). Unresolved flows -We could be seeing residual effects arising from unresolved flow motions and magnetic processes taking place within magnetically active SDO pixels. HMI samples the Fe I line profile (at 6173Å) at only six points in wavelength (cf. Section 2.1). The line shift (velocity), depth, width, magnetic field and continuum intensity are then determined by fitting a symmetric Gaussian to these points. This coarse sampling at the pixel level, and the fact that pixels are of a size comparable to that of granules mean that we are missing spectral-line asymmetries from processes taking place below the pixel resolution, e.g. due to convection. Saar (2009) makes some initial attempts (expanding on Saar 2003) at including convective line asymmetries in plage models.
Zeeman broadening -A final item not included in our model is the direct effect of magnetic fields on the line profiles themselves (Reiners et al. 2013). Due to Zeeman broadening in magnetic regions, the line profiles originating there are wider and shaped differently. If the lines are stronger, they can show enhanced equivalent widths as well. These differences lead to subtle RV changes as active regions rotate and change in number and size. Reiners et al. (2013) showed that RV amplitudes resulting from Zeeman broadening were of the order ∆v B ≈ 300 f (B/1kG) 2 (λ/1µm) 2 m s −1 , where f is the magnetic filling factor, B is the local magnetic field strength, and λ is the wavelength. Note that HMI measures fluxes (i.e. B times area) and the actual B in resolved plage flux tubes is ∼ 1.5 kG (e.g. Buehler et al. 2015). Most plage contains a mix of fluxtubes and fieldfree areas. Thus the true solar magnetic f is close to 1-2 % (note that we are ignoring the ubiquitous weak turbulent fields in this estimate). Adopting λ = 0.5µm, we can estimate ∆v B ≈ 1.7 -3.4 m s −1 . Note that Zeeman broadening is only partly and imperfectly removed by the fit to B since the actual RV dependence is ∝ f λ 2 B 2 . Since proper treatment would require a different calculation of the filling factors, and due to the wavelength dependence, a different computation of RVs, we leave exploring this to a future paper. We note, however, that the Zeeman broadening effect would follow B in phase, and could reduce residuals coincident with large B concentrations during active epochs.

USING |B obs | TO CONFIRM AND CHARACTERISE LONG-PERIOD, LOW-MASS PLANETS
Previous studies have shown that we must account for RV variability in order to detect and characterise lowmass, long-period planets (e.g. Saar 2009;Hall et al. 2018;Meunier & Lagrange 2019, and others). Here, we test whether the unsigned magnetic flux could, in principle, be used to mitigate rotationally modulated RV variations to characterise small planets accurately and precisely (e.g. to better than 10% precision in mass; Zeng & Sasselov 2013). The present analysis is not intended to be a comprehensive exploration of parameter space, nor is it meant to reflect realistic ground-based observing conditions for stellar RV surveys. When facing reality, the most important factor to consider will be the precision to which |B obs | can be measured; we discuss this in Section 8.1. The effects of magnetoconvection, i.e. granulation and supergranulation will also need to be accounted for as they produce RV variability at the m s −1 level (e.g. Meunier et al. 2015).

Procedure
Here, we inject synthetic planet signals in the 8-year long, daily averaged SDO/HMI-derived ∆RV timeseries (see Figure 3), down-sampled to 6 months per year, to simulate the visibility pattern of a star observed from the ground (e.g. Hall et al. 2018, Sect.3.3). For simplicity, we do not remove observations to mimic weather losses, as Hall et al. (2018) estimate that targets visible during the summer only experience a 6% loss of nights due to weather at optimal observing sites. We are left with 1456 daily observations spread over 8 consecutive seasons. We consider three RV semi-amplitudes K: 0.5, 0.3 and 0.1 m s −1 . In each case, we consider a circular orbit, with period P orb = 300.38 days and time of transit t 0 = 1501.92 + 2455318 days (close to the mid-time of the datatset). We choose P orb close to but not exactly at 300 days to avoid overlap with any potential aliases arising from SDO's geosynchronous orbit. Assuming an orbital inclination of 90 degrees and a 1 M host star with a typical uncertainty for a bright, solar analogue of 0.03 M , an RV semi-amplitude of 0.5 m s −1 corresponds to a planet with a mass of 5.2 M ⊕ , 0.3 m s −1 corresponds to 3.1 M ⊕ , and 0.1 m s −1 corresponds to 1.04 M ⊕ .
For each scenario, we fit models consisting of a circular Keplerian and zero-point offset, and either the linear function of |B obs | (of Section 4), or the |B obs | + FF' combination (of Section 5). We account for the remaining residuals by adding the residual RMS values (0.89 and 0.87 m s −1 for the |B obs | and |B obs | +FF' models, respectively; see Table 2) in quadrature to the 1-sigma RV uncertainties (0.1 m s −1 , see Table 1). Rounding up, we obtain an effective RV uncertainty of 0.9 m s −1 for both models. In a real-case scenario one would implement a correlated noise framework (e.g. Gaussian process regression) to ensure that the parameter estimates are as accurate and precise as they can be in the presence of correlated noise. This statistically demanding analysis is beyond the scope of the present analysis, whose primary purpose is to determine whether the planets can be recovered. We assume prior indication of a planet in this range of orbits, e.g. through the detection of one or more transits, so we impose broad Gaussian priors of 10 and 20 days on P orb and t 0 respectively. We maximise the likelihood of each model and determine the best-fit parameter values through an MCMC procedure similar to the one described in Haywood et al. (2014), in an affine-invariant framework (Goodman & Weare 2010). In all cases, the MCMCs reveal a parameter space with multiple local maxima in likelihood. As is done routinely in exoplanet analyses, if the majority of the MCMC chains give the same solution, we remove deviant chains before estimating the model parameters and uncertainties. However, if no single area of parameter space is clearly preferred, we deem the MCMC outcome as a non-detection.

Outcomes
The results of all scenarios are presented in Table 3. In all cases, the estimates for α, β and RV 0 match those determined via the least squares optimisation procedures of Sections 4 and 5 within 1 to 2-σ (see Table 2).
Injected planets with K = 0.5 and 0.3 m s −1 -We recover the planet signal for both of these K amplitudes. The RV amplitudes are estimated accurately within 1-σ. The orbital periods are systematically underestimated by up to 3-σ, as was found by Hall et al. (2018, Fig.5) who performed very similar simulations. The orbital phases t 0 , too, are off from their correct value by up to 3-σ, due to P orb being offset. The planet amplitude is recovered with the same significance using either the |B obs | model or the |B obs | and FF' combination. We show the phasefolded orbit of the injected K = 0.3 m s −1 , recovered using the |B obs | model in Figure 9.
Injected planet K = 0.1 m s −1 -Neither activity models are sufficent to recover signals of this amplitude. The parameter space explored by the MCMC chains shows multiple solutions with similar likelihoods at the 300day period of the injected planet, as well as at 330 days and near 400 days. When we inject a planet at 300 days with K = 0.1 m s −1 , we systematically find that the most likely solution is for a signal with a period of 335 days with an RV amplitude of 0.3 m s −1 . The 335-day signal does not appear to be caused by uneven sampling, as we see complete and uniform coverage at all phases. Stellar activity signals in the 300-day range -We see several significant peaks at 200-400 days in the periodogram of the RVs (Figure 5a). |B obs |, too, exhibits significant and similar (but non-identical) periodicities in the same range ( Figure 5d). Panel (e) shows that the RV residuals (after applying the |B obs | model) exhibit comparatively less power in this period range, but we still see two significant peaks in the 350-500 day range. Meunier et al. (2010b, Fig .10) also detect significant peaks at 300-400 days in solar RV variations of Cycle 23. Their RV variations are estimated using an independent method, using catalogues of sunspot and plage records and magnetograms from SoHo/MDI (which is in a different orbit than SDO). The most likely explanation for the nature of the 335-day signal is that it is a long-term signature of magnetic activity.

From this idealised scenario to stellar observations
In this idealised setup, we successfully retrieve 300-day planet orbits with RV amplitudes down to 0.3 m s −1 , while 0.1 m s −1 signals remain out of reach. To break the 0.1 m s −1 barrier, additional RV signals not well traced by either either |B obs | or FF' will need to be modelled adequately. More generally, these planet injection tests show that in order to access planets with periods of a few hundred days, we will need to model all stellar signals that have similar periods and amplitudes much larger than K. The |B obs | and the combined |B obs | + FF' models perform equally well at this orbital period range. This is expected since we obtain very similar RV residuals with both models (see Sections 4 and 5). The planet retrievals carried out here are a best-case scenario that only considers rotationally modulated RV variations from magnetically active regions. In stellar observations, there will be additional intrinsic variability from magnetoconvection. For this dataset, it would have an rms of 1.1 m s −1 given our cadence and sampling strategy (Meunier et al. 2015). Long-baseline stellar observations are expected to feature RV variations from large-scale meridional circulation, which varies with the magnetic cycle (Komm et al. 1993;Meunier 1999;Ulrich 2010). Although meridional flows have not yet been clearly identified in stars other than the Sun, they produce peak-to-peak RV amplitudes of 1-1.4 m s −1 when viewing the Sun edge-on, and 2.3-3.3 m s −1 in a pole-on scenario (Meunier & Lagrange 2020). Meunier & Lagrange (2020) extend their results to other Sun-like stars and predict peak-to-peak amplitudes of up to 4 m s −1 in stars with strong magnetic cycles seen pole-on. The significance of a planet detection will depend strongly on how precisely we can measure the unsigned magnetic flux in Sun-like stars; see further discussion in Section 8.1. Additionally, there will be instrumental systematics of order 0.1-1 m s −1 for current-generation spectrographs (Fischer et al. 2016). In particular, wavelength calibration and long-term stability remain challenging even in current state-of-the-art spectrographs (e.g. Cosentino et al. 2012). Our RV and |B obs | timeseries are sampled simultaneously, and every night for 8 seasons; ground-based surveys will suffer losses from poor weather (e.g. Hall et al. 2018). These caveats will impact the performance of the technique presented here and diminish the significance of the mass determinations. We note that a systematic investigation of the detectability of low-mass, long-period planets is under way (Langellier et al. in prep.). They explore a broad range of parameter space using HARPS-N RV observations of the Sun, so their analysis is highly complementary to the one here. Langellier et al. (in prep.) inject a wide variety of planet signals into solar observations and recover them by treating magnetic activity using Gaussian process regression. For an in-depth investigation of the impact of observation strategies on
the detectability of Earth-mass, long-period planets, we refer the reader to Hall et al. (2018).

Prospects for measuring |B obs | for other stars
The significance of an exoplanet detection will depend strongly on how precisely we can measure the unsigned magnetic flux in stars other than the Sun. Zeeman Doppler imaging (ZDI; Donati & Brown 1997) has long been used to image the large-scale (polarised) magnetic field structures of stars, particularly those that are fast rotating and much more active than the Sun (see Reiners (2012) and references therein). However, RV variations stem from magnetic fields taking place on much smaller spatial scales than those probed by ZDI. In principle, it is possible to measure small-scale, unsigned magnetic flux by examining Zeeman broadening in magnetically sensitive spectral lines of stars (Robinson 1980;Saar 1988). It has been detected, for example, in the younger, moderately active, faster rotating (P rot ≈ 12 d) K2 dwarf Epsilon Eridani, which has an average unsigned magnetic flux in the range 125 -200 G (Valenti et al. 1995;Lehmann et al. 2015). The Sun's average unsigned magnetic field is twenty times smaller (10 G). Several studies have attempted to measure Zeeman broadening in slowly rotating, relatively quiet late-type stars (e.g. Basri & Marcy 1988;Rueedi et al. 1997;Anderson et al. 2010). More recently, Kochukhov et al. (2020) developed a method employing multiple lines with different Zeeman splitting patterns and making use of differential magnetic intensification of line equivalent widths (Basri et al. 1992;Saar et al. 1992). This work looks very promising for accurately pushing |B obs | detections to lower levels: using only eight spectral lines, Kochukhov et al. (2020) detect filling factors as low as f = 7% and unsigned average fields as low as 220 G. Modeling more lines should yield further improvements in determining |B obs | at low levels. Mortier (2016) combine multiprofile least squares deconvolution (originally proposed by Kochukhov et al. 2010) with singular value decomposition to extract the unsigned magnetic flux from thousands of spectral lines. Their preliminary application to high resolution spectra from HARPS and HARPS-N of an inactive K3 dwarf gives encouraging results.
The aforementioned studies identify several avenues to improve our prospects of detecting Zeeman broadening in old, quiet Sun-like stars. These include: improving modelling of spectral lines (better atomic data and better understanding of line broadening), particularly in the (near-)IR where Zeeman broadening is stronger as it has a squared dependence on wavelength; improving our constraints on convection and turbulence in the stellar atmosphere; and improving our understanding of the impact of line blending and telluric contamination. The magnetic-region filling factor and magnetic field strength are, to some extent, degenerate, but their product, B.f is more easily measurable (e.g. Gray 1984;Saar 1988;Reiners 2012, and references therein). For example, the SPIRou spectrograph at CFHT is providing m s −1 precision spectroscopic observations that extend into the near-IR (Artigau et al. 2011) where the Zeeman effect is stronger, and will therefore provide excellent observations to improve our techniques to measure unsigned magnetic flux in slowly rotating, quiet Sun-like stars.

Prospects improving the F F and F 2 methods
In Section 6 and Appendix C, we noted that the F F and F 2 terms of Aigrain et al. (2012) have intrinsic limitations, largely because they introduce extra limbdarkening-like damping to the predicted ∆RV . It may be possible to improve these methods by applying appropriate "anti-limb-darkening" functions to each, timed to the central meridian passage a magnetic feature. We provide further details in Appendix D. There are several difficulties associated with our proposed correction, but we plan to experiment with these ideas in the near future.

Future improvements to SDO/HMI-RV pipeline
The HMI instrument is not stable over timescales longer than a few days, because it was originally designed for helioseismic observations (see Section 2.1.2 and Figure 2). To correct for this, we look at RV variations, which we obtain by subtracting the velocity of the quiet Sun. In the present analysis, we estimate this quiet-Sun velocity by excluding magnetically active pixels. This is a reasonable approximation since active regions occupy only a few % of the solar disc (5% at the peak of the magnetic cycle). Ideally, one should replace each active pixel by a non-active pixel from the same spatial location (e.g., from a few days before or after the active region's presence). This could be achieved by compiling a quiet-Sun "template" image. Such a fix is unlikely to account for the rotationally modulated RV residuals we are seeing, but may improve the accuracy of ∆RV and ultimately enable us to probe deeper into the physical processes at play in RV variations.

CONCLUSIONS
In this paper, we estimate the disc-averaged, rotationally modulated radial-velocity (RV) variations of the Sun as a star over magnetic cycle 24 from spatially resolved images of the Sun taken by the Helioseismic and Magnetic Imager onboard the Solar Dynamics Observatory (SDO/HMI). To do so, we apply a model that was previously validated against overlapping HARPS-N solar observations by Milbourne et al. (2019). We also estimate the disc-averaged, unsigned (i.e. unpolarised) magnetic flux. Our findings and conclusions are summarised here: -The SDO/HMI-derived RV dataset presented here has high cadence and timespan that covers nearly an entire magnetic cycle, and high SNR (Figure 3). It thus provides a testbed to identify and probe the underlying physical processes that are responsible for rotationally modulated RV variations.
-Periodograms of the Sun's RV ( Figure 5) show that the majority of the power is shared between the rotation period (P ) and its first harmonic (P/2). Both peaks are significant, and each are in fact broad forests of significant peaks that lie up to 2-3 days away from P and P/2.
-We fit RV variations with a linear model of the unsigned magnetic flux, and find that it reduces the rms of RV variations by 62% i.e. a factor of 2.6, from 2.33 m s −1 to 0.89 m s −1 (Figure 6, Section 4). The residuals of the fit display rotationally modulated behaviour, particularly at times of high magnetic activity (Figure 7). To try to account for these residuals, we fit RV variations with a combination of a linear |B obs | term and an FF' term from the method of Aigrain et al. (2012) (Figure 8, Section 5). This yields only modest rms improvements, as the combined model gives a residual rms of 0.85 m s −1 for the full timeseries. We show that the FF' model does not adequately account for RV variations from magnetic regions because it over-accounts for limb darkening (Section 6 and Appendix C), and we propose a correction to potentially improve the performance of the FF' method (Section 8.2 and Appendix D).
-Modelling RV variations with |B obs | and the FF' method allows us to identify additional physical processes responsible for rotationally modulated RV variations. These signals are either missing from the RV model of Milbourne et al. (2019) and Haywood et al. (2016) that we use to estimate RV variations from SDO/HMI images, or they are not well traced by |B obs | or FF', or both. Particularly at high magnetic activity levels, the residuals display significant, rotationally modulated variations at the meter-per-second-level. We discuss physical processes that may contribute to these additional RV variations beyond suppression of convective blueshift and brightness inhomogeneities: horizontal flows (such as Evershed flows, moat flows and active region inflows), flows that are not resolved by SDO/HMI's pixels, and Zeeman broadening (Section 6).
-We inject planet signals to test the performance of the unsigned magnetic flux |B obs | for mitigating rotationally modulated RV variations in surveys of low-mass, long-period planets orbiting Sun-like stars (Section 7). We inject planets with orbital periods of 300 days and RV semi-amplitudes of 0.5, 0.3 and 0.1 m s −1 . The |B obs | model and the combined |B obs | + FF' model give very similar results. The parameters of the planets with K = 0.3 and 0.5 m s −1 are detected accurately to within 1σ of the injected parameters. We do not retrieve injected signals with K = 0.1 m s −1 , because of the presence of an activity-induced signal at 330 days.
We conclude that |B obs | could, in principle, enable us to extract planet signals down to 0.3 m s −1 , but we will also need to model additional RV variations to reach 0.1 m s −1 (Section 7.3). The significance of planet detections in stellar observations will depend crucially on how precisely we may be able to measure |B obs |. Stellar RV observations will also be affected by (super)granulation signals at the m s −1 -level (Meunier et al. 2015), instrumental systematics, and ground-based observing schedules. The most promising avenue to measure the unsigned magnetic flux in slowly rotating, relatively inactive stars is by measuring Zeeman broadening of mag-netically sensitive lines in high-precision spectra (e.g. Kochukhov et al. 2020   showed that for a given level of Ca ii H&K emission, RV variations have a comparatively lower amplitude during the descending phase of the magnetic cycle than in its ascending phase. This is likely because the average latitude |φ| of active regions changes over the course of the magnetic cycle, highest in the ascending phase, and lowest in the descending. Signals produced by an active region depend on the region's position on the solar disc; RV and S-index behave differently with line-of-sight angle and follow different limb-darkening laws, because they originate at different heights in the solar atmosphere (further details in Meunier et al. 2019, Sect.6.1). Note that |B obs | is a line-of-sight observable, and only subject to foreshortening.
To examine these long-term effects, we averaged our timeseries with a 300 day boxcar smooth. This is long enough to smooth out both rotational modulation, and the growth-decay timescales of large active regions, thus concentrating on purely cyclic variation. As shown in Figure A1, we observe a hysteresis between ∆RV and S-index (panel (a)), between ∆RV and the unsigned magnetic flux |B obs | (panel (b)), and between the S-index and |B obs | (panel (c)).
The RV-S-index hysteresis looks qualitatively similar to that observed in the previous solar magnetic cycle by (Meunier et al. 2019, Fig .7). In the ascending phase (dark purple line in Figure A1), plages, at higher |φ| , are (more) limb-brightened, and also increasing in total area, yielding a larger S-index per unit projected plage area than later in the cycle. There is also an RV effect, as the difference in RV between non-magnetic pixels and plage pixels peaks at µ ∼ 0.9 (Palumbo et al. 2017). This azimuthal ring lies entirely within φ = ±25 o , the approximate "active latitude zone". Thus, as the cycle progresses and active-region average latitude decreases, the per-pixel average RV difference with the quiet Sun decreases. In the ascending phase, this is more than compensated by the increasing filling factor, but once f plage starts declining in the descending phase, ∆RV drops steadily.
The hysteresis between ∆RV and |B obs | is in some ways simpler to understand. At the 300 day level of smoothing, the twin cycle maxima of cycle 24 at t ∼ 500 days (the weaker Northern hemisphere peak) and t ∼ 1600 d (the stronger Southern hemisphere peak) are both flattened into one slow increase in |B obs |. With an initial φ ∼ 25 o , then mostly decreasing throughout the cycle, the average net RV per plage pixel should also be decreasing. This is counterbalanced, however, by the filling factor, which increases more quickly than RV decreases in both the ascending phase and the peak(s) of the cycle. Thus, ∆RV continues to increase with |B obs | in these phases, more slowly at maximum when |B obs | growth is also reduced, then only finally reversing in the decline phase, when both |B obs | and ∆RV plummet.
Panel (c) of Figure A1 shows a hysteresis between S-index and |B obs |, which is likely due to the different limbdarkening behaviours, with the S-index getting a boost from limb brightening when in the higher |φ| ascending phase.
Further study is needed to better understand the differences in projection effects between these three observables, to correct for them and therefore obtain tighter correlations between ∆RV , S-index and |B obs |.
Attempt to account for the hysteresis between ∆RV and |B obs | -To capture the information in the hysteresis of Figure A1 (b), we fitted the ascending and descending phases of the magnetic cycle separately, using two |B obs | terms and two zero-point offsets. The magnetic cycle has a double-peaked shape, because the active region bands reach maximum activity levels at slightly different times. We identified the peak of the magnetic cycle, i.e. the point separating the ascending and descending phases as the minimum in magnetic flux and active-region coverage between these two peaks, at JD = 2456957.5 (day 1639 of the timeseries shown in Figures 3 and 6). We obtain different model parameters for the ascending (α = 11.15 ± 0.01 m s −1 , RV 0 = -7.98 ± 0.01 m s −1 ) and descending phases (α = 8.64 ± 0.01 m s −1 , RV 0 = -6.26 ± 0.01 m s −1 ). However, the residual rms (over the full cycle) is 0.83 m s −1 (64% reduction in RV variations), which is a only a small improvement compared to fitting the full cycle as one (rms = 0.89 m s −1 , 62% reduction). This is an improvement, but since it is small and adds more parameters and complexity, we leave this avenue open for future investigations.

B. CORRELATIONS BETWEEN FILLING FACTORS AND ACTIVITY INDICATORS
In Figure B2, we show the plage filling factor as a function of |B obs | and S-index, and the spot filling factor as a function of |B obs |. The magnetic filling factors are estimated according to Eqn. 3.2. Figure A1. Hysteresis between RV variations and: (a) Ca ii H&K emission as measured by the S-index, (b) unsigned magnetic flux |B obs | and (c) between |B obs | and the S-index. All timeseries are smoothed over 300-day bins to reveal the long-term hystereses. The points are colour-coded according to time: the ascending phase (increasing activity) is in purple while the descending phase (decreasing activity) is in green and yellow. Here, we show why the FF' method cannot match RV variations perfectly. Consider an equatorial spot. If we assume, for simplicity, linear limb darkening, and solar inclination i = 90 o , the flux from the spot of area A can be written as: where θ is the angle of the surface normal to the line of sight (at disk center), and s is the linear limb darkening coefficient for the spot. The derivative with respect to θ is: To match the RV change due to the rotation of a spot, this simple model (based on Saar & Donahue 1997) yields: ∆v spot = −A sin θ cos θ(1 − s + s cos θ).
The sin θ term captures the RV deficit, and is weighted by the projected area of the spot (cos θ) and its limb darkening, (1 − s + s cos θ).
Thus, the F F method captures ∆v spot but adds an additional limb-darkening-like term. This leads to systematic effects, underestimating ∆v spot progressively more and more as the spot moves away from disk center. Although we do not use the F 2 term proposed by Aigrain et al. (2012) to correct for the convective suppression arising primarily in plage, we note that it is similarly flawed. Following a similar analysis: ∆v plage = A cos 2 θ(1 − p + p cos θ), where p is the linear limb darkening coefficient for plage. But F 2 yields: or, F 2 /A = ∆v plage (1 − p + p cos θ), which contains an extra limb-darkening term.
In the case of a spot, its limb darkening s is unlikely to differ significantly from the quiet Sun value q . This is because limb darkening is wavelength dependent, and the strength-weighted average for HARPS-N RV lines is perhaps λ ∼ 500 nm, not much different from the HMI line (λ = 617.3 nm). FF' captures the RV perturbation due to a spot crossing the disc effectively, but then applies a second,additional, stronger limb darkening.
We thus warn that, while the FF' term is partially successful, it (and the related F 2 term) are also flawed, since they intrinsically add extra, unwanted limb-darkening-like corrections. Therefore, one cannot expect the FF' term (or the F 2 term) to perfectly account for the RV variations of magnetic features.
A similar formula can similarly be applied to plage (detected, e.g. in Ca ii H&K), by using the limb darkening coefficient p . We note that the equivalent correction to F 2 for convective suppression, which occurs predominantly in plage (rather than spots) is: F 2 corrected,plage = F 2 /(1 − p + p cos θ).
In order to apply the above corrections, one has to know when individual active regions cross the solar/stellar disc.
To some extent, one can track the meridian passage of active regions on solar/stellar surfaces via monitoring of disc-averaged photometric and Ca ii H&K emission. However, an active region often contains both spots and plage, resulting in a mixed, degenerate photometric signal. Also, it would be difficult to apply this correction in terms of timing when multiple active regions are present. Additionally, we note that s and p remain poorly known due to lack of realistic models. Better values may be derivable from 3D-MHD simulations (e.g., Cegla et al. 2018). It is also difficult to choose an appropriate value of for spectra that span hundreds of nm. To summarise, applying a correction to the FF' term could potentially improve fits to RV variations incurred by spots and plage, but there will likely remain residual RV variations.