The magnetically quiet solar surface dominates HARPS-N solar RVs during low activity

Using images from the Helioseismic and Magnetic Imager aboard the Solar Dynamics Observatory (SDO/HMI), we extract the radial-velocity (RV) signal arising from the suppression of convective blue-shift and from bright faculae and dark sunspots transiting the rotating solar disc. We remove these rotationally modulated magnetic-activity contributions from simultaneous radial velocities observed by the HARPS-N solar feed to produce a radial-velocity time series arising from the magnetically quiet solar surface (the ‘inactive-region radial velocities’). We find that the level of variability in the inactive-region radial velocities remains constant over the almost 7 year baseline and shows no correlation with well-known activity indicators. With an RMS of roughly 1 m s − 1 , the inactive-region radial-velocity time series dominates the total RV variability budget during the decline of solar cycle 24. Finally, we compare the variability amplitude and timescale of the inactive-region radial velocities with simulations of supergranulation. We find consistency between the inactive-region radial-velocity and simulated time series, indicating that supergranulation is a significant contribution to the overall solar radial velocity variability, and may be the main source of variability towards solar minimum. This work highlights supergranulation as a key barrier to detecting Earth twins.


INTRODUCTION
The radial-velocity (RV) method is one of the most valuable tools in the planet hunters' arsenal, as both a detection method, with over 1000 confirmed exoplanets to date, and a vital follow-up tool for transit missions such as the Kepler (Borucki et al. 2010) , Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014), and upcoming PLAnetary Transits and Oscillations of stars (PLATO; Rauer et al. 2014) missions.Precise RV measurements are often necessary to independently confirm planetary candidates and to determine planetary masses.As instrumental precision improves, the intrinsic RV variability of the host star becomes the primary obstacle to detecting and characterising low-mass planets (see Crass et al. 2021, and references therein).Stellar variability can mimic non-existent planets (Rajpaul et al. 2016), lead to inaccurate mass measurements (e.g., Blunt et al. 2023;Meunier et al. 2023), or completely obscure the signal of a low-mass planet.
Photospheric magnetic activity primarily causes RV variations via two processes.Firstly, strong magnetic fields act to inhibit convective blueshift of rising stellar material, imparting a net redshift onto the overall stellar RVs (Meunier et al. 2010a;Dumusque et al. 2014).Secondly, bright and dark active regions break the axial symmetry of the rotating stellar disc.This results in an imbalance of flux from the approaching and retreating stellar limbs, causing a net Doppler shift (Saar & Donahue 1997).As these phenomena occur in isolated, magnetically active regions which can persist for multiple rotation cycles, the resulting RV variability typically displays quasi-periodic behaviour with amplitudes of several m s −1 .
Traditional activity-mitigation techniques include the use of magnetically sensitive activity indicators which correlate with activity-induced RVs (e.g., Boisse et al. 2011;Haywood et al. 2022), modelling RV variations with Gaussian processes to leverage the quasi-periodicity of isolated active regions rotating with the stellar surface (e.g., Haywood et al. 2014;Rajpaul et al. 2015), and using simultaneous photometry and ancilliary time series to estimate the effect of active regions (e.g., Aigrain et al. 2012;Barragán et al. 2022).
Whilst these techniques have proven effective in mitigating the RV variability originating in isolated, magnetically active stellar regions, they do not take into account variability originating from the magnetically quiet stellar surface.Nevertheless, there are a number of phenomena on the quiet stellar surface which can cause RV variability.In particular, observations and simulations of solar convective flows (e.g., granulation and supergranulation) have shown RV variability on the order of 0.5-1 m s −1 (e.g., Meunier et al. 2015;Meunier & Lagrange 2019;Meunier & Lagrange 2020;Al Moulla et al. 2023).In this paper, we use spatially resolved solar images to calculate the RV impact of isolated magnetically active regions (following Meunier et al. 2010b;Milbourne et al. 2019).We remove these contributions from Sun-as-a-star RVs from the HARPS-N spectrograph to isolate and directly characterise the RV impact of the magnetically quiet solar surface (hereafter the 'inactive-region RVs').This paper is structured as follows.In Section 2 we describe the data used in this work and we isolate the inactive-region RVs in Section 3. We analyse the statistical properties of the inactive-region RVs in Section 4 and demonstrate that the inactive-region RVs are not simply an instrumental artefact (Section 5).It is worth noting that, up until this point, our analysis is purely statistical and therefore independent of the physical cause of inactive-region RVs.In Section 6, we show that the inactive-region RVs match simulations of oscillation, granulation, and supergranulation and propose supergranulation as the dominant source of variability.We conclude in Section 7.

HARPS-N Solar telescope
The solar telescope at the Telescopio Nazionale Galileo is a 7.6-cm achromatic lens which feeds the sunlight to an integrating sphere and then into the High Accuracy Radial velocity Planet Searcher for the Northern hemisphere (HARPS-N) spectrograph (Dumusque et al. 2015;Phillips et al. 2016;Collier Cameron et al. 2019;Dumusque et al. 2021).Observations are taken almost continuously during the day time, with 5-min integration times in order to average over the solar p-mode variations.RVs are calculated from the extracted HARPS-N spectra using the HARPS-N Data Reduction Software (DRS, Dumusque et al. 2021).By using the same instrument and DRS used for night-time exoplanet searches, the solar feed provides a realistic Sun-as-a-star RV time series and has yielded information about sub-m s −1 level instrumental systematics present in the HARPS-N spectrograph (Collier Cameron et al. 2019).To better emulate night-time stellar observations, we average three successive observations to give 15-minute effective integration times.
In this paper, we use a data set spanning from 2015-Jul-29 to 2021-Nov-12.The RVs have been transformed into the heliocentric rest frame, allowing the data to be free of planetary signals.Additionally, a daily differential extinction correction was done following the model described in Collier Cameron et al. (2019).The HARPS-N Solar Telescope operates continuously, observing from under a plexiglass dome.Data affected by clouds or other bad weather thus need to be accounted for afterwards.To remove potential bad data points, we use two metrics for each observation.Firstly, data quality factor  is calculated from a mixture model approach described in Collier Cameron et al. (2019).This probability value can be between 0 and 1 where 1 represents the reliable data.Secondly, we extract the maximum and mean count of the exposure meter from the headers of the fits files.The ratio of max and mean, , should be close to 1 for uniform exposures (and is always higher than 1 by definition).A first quality cut is done using the conditions  ≥ 0.99 and  < 1.5.From the remaining data points, we fit the histogram of  with a Gaussian function.A slight delay in the shutter opening causes this distribution to resemble a Gaussian, rather than the more expected pile-up around 1. In a second cut we remove data where  is higher than 3 standard deviations above the mean.Finally, as a last cut, we perform a 5- clipping on the remaining RVs.The cuts made are very strict and may indeed remove reliable data too, but it is done to ensure the very high quality of data.The remaining data-set has 20 435 observations.

SDO/HMI
The Helioseismic and Magnetic Imager aboard the Solar Dynamics Observatory (SDO/HMI) uses six narrow bands around the magnetically-sensitive 6173 Å Fe I line.A Gaussian fit of the narrow-band fluxes as a function of wavelength is used to derive continuum intensitygrams, magnetograms, and line-of-sight velocity maps (Dopplergrams), each spatially resolved with approximately 1" resolution.In this work we use 720-second exposures, every four hours from 2015-Jul-29 to 2021-Nov-12, totalling 12 434 sets of images.We briefly describe how we calculate RVs from these images.For more details, we refer the reader to Haywood et al. (2016) and Milbourne et al. (2019) which describe in detail the process of extracting disc-averaged physical observables from SDO/HMI images.
Owing to poor long-term stability of the SDO/HMI instrument, we are not able to directly calculate the solar RVs by simply intensity-weighting the Dopplergram pixels.Fig. 2 of Haywood et al. (2022) shows the effects of long-term systematics when calculating the RVs in this way.There are clear jumps in the RVs caused by instrumental effects.This means that we cannot use the Dopplergram to produce long-baseline Sun-as-a-star RV measurements.To overcome this instrumental limitation, we use a physically-motivated model to estimate the RV variations arising from two processes on the Sun.In doing so, we calculate the RV variability relative to the quiet Sun (Meunier et al. 2010b;Haywood et al. 2016;Milbourne et al. 2019;Haywood et al. 2022).
The first of these is a photometric shift, RV phot , caused by a bright plage or dark sunspot breaking the rotational symmetry of the solar disc.In the absence of inhomogeneities on the solar disc, the contribution of the blue-shifted approaching limb and red-shifted retreating limb are in balance, with no overall effect on the solar RVs.
If a dark or bright region is present on one of the limbs, the Doppler contribution from that limb is diminished or enhanced, respectively, resulting in an overall RV shift.The second component, RV conv , arises from regions with large magnetic flux suppressing convection on the surface of the Sun.The reduction in the convective blue-shift in these regions imposes a net red-shift onto the overall solar RV.By applying continuum intensity and magnetic flux thresholds described by Haywood et al. (2016), respectively, we isolate the RV contribution of active regions on the Sun.
Since SDO/HMI derives its RVs from a single iron line, as opposed to a whole spectrum, Haywood et al. (2016) and Milbourne et al. (2019) fit a linear combination of these components to solar RVs from HARPS and HARPS-N, respectively, to account for the different scaling of each component.They found that RV conv is the dominant of these two effects in the Sun.Milbourne et al. (2019) only consider the RV contributions from active regions larger than 20 µHem.Those authors find that this cutoff obviates the need to vary the relative contributions of the photometric and convective effects over time.Owing to the fact that the two components forming this RV series are driven by magnetic processes, and that RVs are calculated relative to the quiet Sun, we will refer to the RV series derived from the SDO/HMI images as the 'magnetic activity' RVs.Since the formal uncertainties on the magnetic activity RVs are significantly lower than the uncertainties for the HARPS-N RVs (Haywood et al. 2022), we will refrain from extensive error analysis of the magnetic activity RVs.

ISOLATING RV CONTRIBUTIONS FROM MAGNETICALLY-INACTIVE REGIONS
In this paper, we are interested in the RV variability originating from magnetically quiet regions of the solar surface.To isolate these contributions, we subtract the SDO/HMI RVs (which are only sensitive to large magnetically active regions) from the HARPS-N RVs (which probe the entire solar disc).The resultant time series of residual RVs therefore probes the RV signals not arising from the large, magnetically active regions identified by Milbourne et al. (2019).Throughout this paper, we will refer to this time series as the inactive-region RVs.To account for the difference in cadence of SDO/HMI and HARPS-N, we interpolate the magnetic activity RVs onto the HARPS-N timestamps.In principle, this could introduce a signal into the inactive-region RVs, since we are smoothing over all variability in the magnetic activity RVs over timescales shorter than 4 hours.Recalling, however, that the magnetic activity RVs encode variability occurring on rotational timescales, and noting that shorter-timescale processes (e.g., oscillations/granular flows) do not contribute to the magnetic activity RVs, we are justified in this approach.We provide a quantitative justification in Appendix A. It is worth highlighting that the model of Milbourne et al. (2019), which we use in this paper, does not include the effect of the smallest active regions.However, those authors show that by only considering the largest active regions, they reproduce the majority of the rotationally-modulated, activity-induced RVs, without having to introduce an arbitrary trend.The models of Milbourne et al. (2019) reduce the HARPS-N RV scatter from 1.65 m s −1 to 1.21 m s −1 and 1.18 m s −1 with and without the 20 µHem active-region cutoff, respectively 1 .We are therefore justified in describing the RV time series produced by removing this contribution from the disc-integrated HARPS-N RVs the 'inactive-region' RVs.We do note however, that there will inevitably remain a small component in the inactive-region RVs that corresponds to the smallest-scale magnetic regions.
To avoid any effects of long-term instrumental effects, we de-trend the inactive-region RVs with a 100-day rolling mean, following the treatment of Al Moulla et al. (2023).We note that such a de-trending preserves both the amplitude and timescale of any short-term variability present in the RV time series, and so will not affect the results of this work.We opt to smooth the inactive-region RVs, as opposed to the HARPS-N RVs, since we are interested in the difference between the RVs measured by HARPS-N and those derived from SDO/HMI images.Directly smoothing the HARPS-N RVs would risk smoothing over the rotationally modulated activity signal present in both series, thereby imposing an artificial signal in the residuals when we subtract the un-smoothed magnetic activity RVs.
Fig. 1 shows the correlation between the daily mean RV and both log R ′  and unsigned magnetic flux (Haywood et al. 2022) values for the HARPS-N RVs and the inactive-region RVs, as well as the respective Pearson correlation coefficients.We see an almost total removal of any correlation between the RVs and the activity indices when the magnetic activity RVs are subtracted from the HARPS-N RVs, indicating that we are successfully removing the majority of the RV signal caused by magnetic activity.
In addition, following Al Moulla et al. (2023), we calculate the velocity power spectrum (PS) of the HARPS-N RVs and the residual RVs as follows.Each RV time series can be represented as a linear combination of sinusoidal components, where () and () are frequency-dependent coefficients, and  is a constant offset.The velocity power spectrum is defined as and is converted to velocity power spectral densities (PSDs) by multiplying by the effective length of the observing window (Kjeldsen et al. 2005).
Fig. 2 shows the two PSDs as well as the ratio of power in the HARPS-N RVs to the inactive-region RVs, calculated for logarithmically-spaced bins in frequency.We see evidence of decreased power in the inactive-region RVs at timescales comparable to and longer than the solar rotation period, whereas the two PSDs are similar for shorter timescales.This further suggests that we are successful at removing most of the rotationally modulated variability, whilst retaining the short-timescale variability characteristics of the HARPS-N RVs.

INACTIVE-REGION RVS EXHIBIT A CONSTANT LEVEL OF VARIABILITY
Having established that we are successful in removing the RV signature of large magnetically active regions (Section 3), we now analyse the inactive-region RVs.To investigate how the inactive regions vary throughout the solar cycle, we compare the RMS of the RVs with the mean sunspot number, taken from SILSO World Data Center (2023), for 45 50-day intervals, as a proxy for overall activity level.In Fig. 3 we show that there is no correlation.Additionally, despite the sunspot number changing by more than two orders of magnitude, the standard deviations of the RVs vary only slightly.This shows that the inactive-region RV variability remains relatively constant at a wide range of magnetic-activity levels.In Appendix  The ratio of HARPS-N RV PSD to inactive-region RV PSD.Although there is still a noticeable feature in the residual RV PSD at the rotation period, and at one day and its alias, we see lower power at timescales longer than ∼ 10 days when compared to the full HARPS-N RVs.By contrast, the two PSDs are very similar at shorter timescales.This indicates we are successful in removing much of the variability at timescales around or above the solar rotation period, but we leave the shorter timescale variability largely unaffected.B, we show that the 17 per cent scatter in the RMS shown here is consistent with being dominated by sampling noise.
To investigate how the inactive-region RVs behave at different timescales, we use structure functions (SFs; see e.g,.Simonetti et al. 1985;Sergison et al. 2020;Lakeland & Naylor 2022).Structure functions quantify how the variability within a time series changes with timescale.In Appendix C1 we show that, for a continuous, uncorrelated signal, the structure function is equivalent to where  is the RMS of the signal.We therefore choose to use √︃ 1 2 SF to quantify the RV variability at a given timescale , emphasising the connection with the RMS.We provide a more in-depth description of structure functions and their key properties required for this analysis in Appendix C.
We now segment the RVs by calendar year and show the corresponding structure functions in Fig. 4. We opt to use longer sub-series than the 50-day samples in Fig. 3 to ensure that each timescale is well sampled, with many pairs of data points contributing to the calculation for each timescale.We see that despite covering a wide range of magnetic activity levels, the inactive-region RVs have similar variability at all timescales.We also note that, for each 1-year segment, there is no significant increase in variability at timescales longer than ∼ 3 days.The typical year-to-year variation is 5-10 per cent.This is comparable to the typical lower-bound sampling error we estimate in Appendix C2, indicating that the true year-to-year variation is very low.We contrast this with the behaviour of the magnetic-activity RVs over the same baseline.We calculate structure functions of the magnetic-activity RVs (again segmented by calendar year).We plot structure functions calculated from both the inactive-region RVs and magnetic activity RVs on the same axis scales in Fig. 5 to highlight the contrast.Whilst the magnetic-activity RVs vary by more than an order of magnitude, the inactive-region RVs exhibit an almost unchanging level of variability.Fig. 5 also demonstrates that, for the majority of the 7-year baseline, the variability in the residual RVs is comparable to or larger than the variability shown in the magnetic activity RVs.This shows that, for the quieter period of solar cycle 24, the majority of the RV variability originated from magnetically-inactive regions, as opposed to magnetically-active regions.The two contributions are approximately equal during the 2016-2017 period, where the average sunspot number was 39.8.Fig. 6 shows the distribution of yearly-averaged sunspot numbers, dating back to 1700 from the World Data Center SILSO, Royal Observatory of Belgium, Brussels (SILSO World Data Center 2023).The vertical line corresponds to the 2016-2017 average.Roughly 35 per cent of recorded years have a lower average sunspot number.This implies that the inactiveregion RVs may dominate solar RVs for roughly one third of the time.

INACTIVE-REGION RVS ARE NOT DOMINATED BY INSTRUMENTAL EFFECTS
It is in principle plausible that the variability shown in Fig. 4 at timescales of around two days is due to the instrumental characteristics of the HARPS-N spectrograph, rather than genuine solar effects.To address this concern, we repeat our analysis with RVs from the NEID solar feed (Lin et al. 2022) taken between 2021-Jan-01 and 2022-Jun-13.We use RVs reported by version 1.1 of the NEID Data Reduction Pipeline2 and select observations that are free from known issues due to the daily wavelegnth calibration, a known cabling issue or evidence of poor observing conditions based on the pyrheliometer and exposure meter data.We select only observations taken before 17:30 and with airmass less than 2.25 to reduce calibration and atmospheric effects.Due to the different wavelength ranges and sensitivities of the NEID and HARPS-N spectrographs, the photometric shift, RV phot , and suppression of convective blueshift, RV conv , (see Section 2.2) may have different contributions to the overall RVs.To account for this, following Milbourne et al. (2019), we fit a linear combination of these components to the NEID RVs.
It is worth noting that significant intra-day systematics appear in the NEID solar RVs due to, amongst other effects, unaccounted-for differential extinction.We model the sub-day variations with a bestfit cubic trend to the 1-day phase-folded RVs, which we then subtract from the NEID solar RVs.Whilst this is not as sophisticated as the treatment of Collier Cameron et al. ( 2019), we wish simply to recover the general variability properties seen in Fig. 4 and so a full treatment of NEID instrumental effects is both beyond the scope of this paper and not necessary for our purposes.As with the HARPS-N RVs, we use the recalculated SDO/HMI RVs and NEID RVs to produce a time series of inactive-region RVs.The HARPS-N and NEID RVs temporally overlap between 2021-Apr-01 and 2021-Nov-12 and we compare the HARPS-N-and NEID-derived inactive-region RV time series during this overlap to investigate the effects of instrumental systematics.Fig. 7 shows structure functions of the two inactiveregion RV series.We find that the variability of these two time series at all timescales typically differs by less than 10 per cent, comparable to the lower-bound sampling uncertainty we estimate in Appendix C2.This indicates that the results of section 4 are genuine, and are not simply instrumental artefacts, and that any effects of instrumental differences are smaller than the intrinsic variability of the solar RVs; as such, we will focus the remainder of our analysis on the longerbaseline time series of inactive-region RVs derived from HARPS-N.We note that this analysis qualitatively agrees with the results of Zhao et al. (2023), who compare solar RVs from state-of-the-art spectrographs, including HARPS-N and NEID.Those authors find remarkable agreement between the different instruments.

SUPERGRANULATION AS THE DRIVER OF INACTIVE-REGION RVS
In Section 4, we demonstrated that the inactive-region RVs have a roughly constant level of variability over the decline of solar cycle 24.To determine the cause of this variability, we compare the inactive-region RVs to a simulation of the oscillation, granulation, and supergranulation of a G2 star from Meunier & Lagrange (2019), based on the technique of Harvey (1984) and the results of Meunier et al. (2015).The simulated RV time series is averaged over 15 minutes to best match the HARPS-N observations.2019) (pink), the combined simulated time series (red), and the inactive-region RVs (black).The granulation (Gra) and supergranulation (Sgr) components are re-scaled by a constant factor to match the observations; the oscillation (Osc) component is assumed to be negligible (e.g., Chaplin et al. 2019) and so is not re-scaled.However, we opt to include the oscillation component for completeness.Since granulation and supergranulation dominate at different timescales, we are able to re-scale the simulated RVs to constrain the contribution of each process by matching the power in the inactive-region RVs at short and long timescales, respectively.By combining and re-scaling the simulated RVs in this way, we are able to accurately reproduce both the overall variability amplitude of the inactive-region RVs, and the timescale of the inactive-region RVs.Meunier et al. (2019) found that the amplitude for the RV signals from granulation is 0.8 m s −1 , though there is evidence that an amplitude of 0.4 m s −1 is a more appropriate level (e.g., Sulis et al. 2020Sulis et al. , 2023)).We see in Fig. 8 that the choice to use the lower level of granulation variability is justified; a granulation RV time series with RMS of 0.8 m s −1 would have more power at short timescales than is seen in the true RVs.We find granulation and supergranulation variability amplitudes of 0.37 m s −1 and 0.86 m s −1 , respectively.We therefore find a similar level of granulation to the lower level of Meunier & Lagrange (2020), and to the 0.33 m s −1 prediction of Dalal et al. (2023).We note that the 0.86 m s −1 level of supergranulation we find is around 25 per cent higher than the 0.68 m s −1 level found by Al Moulla et al. (2023), though the granulation amplitudes are similar.Given the difference in methodology between this work and that of Al Moulla et al. (2023), who use asteroseismology techniques to analyse Sun-as-a-star RVs from the HARPS and HARPS-N spectrographs, it is unsurprising that there is a difference in the exact amplitude of supergranulation derived.Both analyses demonstrate that supergranulation contributes a significant fraction of the total solar RV variability.
Our interpretation that the inactive-region RVs are primarily caused by granulation and supergranulation is supported by the fact that the level of variability is constant at various levels of the solar cycle (Figs. 3 & 5).Muller et al. (2018) used HINDODE/SOT images of the Sun from 2006-Nov to 2016-Feb to investigate photometric properties of solar granules.Over this period, they did not find changes in either the granulation contrast or granulation scale of the 3 per cent level at which they were sensitive, indicating no significant variation of granulation properties over the solar cycle.
We note that Fig. 8 does not include any instrumental systematics.Whilst this makes it difficult to evaluate numerical uncertainties on the level of supergranulation quoted above, it is unlikely to significantly change the amplitude.Firstly, a white noise profile (such as from photon noise) would inject power at all timescales and appear flat as a structure function.Since we constrain the granulation amplitude by using the short-timescale power, we can be confident that we have not underestimated the uncorrelated noise.Thus the re-scaling we apply to the supergranulation time series to match the long-timescale power is reasonably robust against the level of photon noise, which is estimated to be as low as 0.24 m s −1 (Al Moulla et al. 2023).Secondly, we show in Section 5 that the level of variability of inactive-region RVs derived from HARPS-N and from NEID observations are consistent within typical sampling errors at all timescales, indicating that the variability observed is not dominated by instrument-specific artefacts.This shows that the supergranulation amplitude derived from Fig. 8 is representative.
Fig. 8 shows that the combined simulated RVs show slightly reduced power on intermediate timescales (between a few hours and a few days) as compared to the inactive-region RVs.Unaccountedfor correlated instrumental noise (e.g., an imperfect correction for differential extinction) could inject variability at intermediate timescales to account for this discrepancy.Alternatively, the intrinsic variability spectrum of supergranulation may differ slightly to that of granulation.The simulations of Meunier et al. (2015) assume that supergranules evolve similarly to granules.The resulting simulated time series therefore share a common variability power law.A relaxation of this assumption may account for the mismatched power at intermediate timescales.Additionally, the smallest-scale active regions, unaccounted for in the model of Milbourne et al. (2019), may inject RV variability at these intermediate timescales.
Regardless of the cause of this discrepancy, the agreement of the inactive-region RVs with both the amplitude and timescale of the simulations of Meunier et al. (2019) and analysis of Al Moulla et al. (2023) offers good evidence that the variability seen in the inactive-region RVs is predominantly caused by supergranulation at timescales of a few days and longer.Given that these RVs exhibit a relatively constant level of variability over the solar cycle (Figs. 4 & 5), we provide evidence that supergranulation dominates the solar RVs on the approach to cycle 24 minimum and therefore poses a significant barrier to the detection of Earth twins whose Doppler shifts are ∼ 0.1 m s −1 .

CONCLUSIONS
We have generated solar RVs from SDO/HMI images following the technique of Milbourne et al. (2019).The two components accounted for in this time series are the RV signatures arising from the suppression of convective blueshift in magnetically active regions and the effect of photometrically imbalanced regions transiting the rotating solar disc.We subtract these 'magnetic activity' RVs from the disc-integrated RVs obtained from the HARPS-N solar feed to isolate the RV contribution from the magnetically quiet solar surface.
We find that the resulting 'inactive-region' RV time series shows no correlation with the well known activity indicators log R ′  and unsigned magnetic flux.This, combined with the reduced Fourier power at timescales longer than roughly half a solar rotation period, show we are successful in removing the majority of the RV signal arising from the active regions on the Sun.
We show that the inactive-region RV variability is relatively stable, showing no significant change in level with either time or sunspot number.This contrasts starkly with the magnetic-activity RVs, which show a change in variability of more than an order of magnitude over the same baseline.This relative constancy in the inactive-region RV variability means that, when the Sun is relatively quiescent, the solar RVs are dominated by the RV signal originating in inactive regions, not active regions.We find that the two contributions are roughly even when the average sunspot number is around 40.This level of magnetic activity represents the 35 th percentile recorded in archival sunspot data dating from 1700.
Finally, we compare the inactive-region RV time series to recent models of oscillation, granulation, and supergranulation (Meunier & Lagrange 2019).We find good agreement with both variability amplitude and timescale between the inactive-region and simulated RV time series.We find that granulation and supergranulation induce RV variability with amplitudes of 0.37 m s −1 and 0.86 m s −1 , respectively.In doing so, we provide empirical evidence that supergranulation can dominate solar RVs and therefore pose a significant barrier to detection of Earth twins.This is in agreement with the findings of Meunier et al. (2023).
To match the high cadence of the HARPS-N solar RVs, we calculate RVs from the SDO/HMI images (see section 2.2) every four hours, and interpolate these RVs onto the timestamps of the HARPS-N observations.
To assess the impact that interpolation has on the signal, we produce a month-long time series of RVs, generated from every available 720s exposure set of SDO/HMI images from 2017-Jan-01 to 2017-Feb-01.From this high-cadence RV series we produce a sub-sampled RV series with six RV measurements per day, to match the long-baseline series we use in the bulk of this paper.We then interpolate this sub-sampled RV series back onto the timestamps of the original high-cadence time series.
We find that the median absolute deviation between the highcadence and interpolated RV series is 2.1 cm s −1 .This is significantly lower than any other level of variability we discus in this paper, justifying that our choice to interpolate the RVs calculated from the SDO/HMI images does not impact the results of this paper.

APPENDIX B: THE EFFECT OF SAMPLING NOISE ON THE RMS OF THE INACTIVE-REGION RVS.
In Fig. 3, we see that a variation in the RMS of inactive-region RV sub-series.Whilst each RMS is calculated from a relatively large number of observations, we show in Fig. 4 that the inactive-region RVs become uncorrelated at timescales longer than a few days.This relatively large timescale, compared to the 50-day baseline of each sub-series, means that the sampling error on the RMS is larger than would be expected by simply considering the number of observations.
To test this, we take the simulated supergranulation time series of Meunier et al. (2019) and calculate the RMS for the same 50-day sub-series as in Fig. 3.This time series is constructed in a stationary manner, so that any variation seen between values of the RMS can with good confidence be attributed to sampling error.Fig. B1 shows the distribution of sub-series RMS values.As we are only interested in the fractional difference between the individual RMS values, we normalise by the mean RMS.We see a relative scatter of the RMS values,   < > , of 14 per cent, where ) This value of   < > is comparable to the inactive-region   < > of 17 per cent.This indicates that the majority of the variation seen in Fig. 3 is due to statistical noise.

APPENDIX C: STRUCTURE FUNCTIONS
The structure function (SF) is an analysis tool which shows how the variability within a time series changes as a function of timescale.It has had success in analysing photometric time series data (e.g., Hughes et al. 1992;de Vries et al. 2003;Sergison et al. 2020).SFs use a pairwise comparison of points in a time series  () to calculate the level of variability as a function of timescale.To calculate the variability of points with a given separation  = | 1 −  2 |, we define the structure function as where the angular brackets indicate an average over all pairs of observations with separation .In practice, SF is calculated for logarithmic bins in .That is to say, in the case of discrete data, the structure function is calculated as where the summation is taken over all pairs of data points which are separated by  1 < |  −  | <  2 and  p is the number of such pairs.below the noise floor, so the variability is dominated by uncorrelated instrumental noise.In Region 2, there is a power-law increase in the structure function, the gradient of which is determined by the spectrum of variability being investigated.Region 3 shows timescales longer than the dominant variability timescale in the data.Therefore, after this timescale has been reached, the value of the structure function stops increasing, again becoming uncorrelated.The transition between Regions 2 and 3, referred to as the 'knee', tells us the characteristic variability timescale.Fig. 5 of Lakeland & Naylor (2022) demonstrates that, for a sinusoidal signal, the knee is at roughly a quarter of the period.This highlights that the timescale identified by a structure function, whilst related to any period present in a signal, is not the same as a period.Another feature of Fig. C1 is that the structure function becomes larger as we move towards longer timescales.Whilst not strictly a cumulative measure, this behaviour is typical for most types of variability. 3This property highlights how short-timescale variations still impact the overall level of variability at longer timescales.
An advantage of the structure function over other analysis tools is the direct relationship between the structure function and the level of variability at a given timescale.We show in Appendix C1 that, for an uncorrelated signal, the structure function will tend to the value 2 2 , where  is the standard deviation of the data.This clear relationship between structure function and standard deviation allows a much more direct interpretation of the level of variability within a signal for a given timescale than, say, a Fourier power  Sergison et al. (2020).Region 1 shows white noise present at short timescales, often interpreted as instrumental noise.Region 2 shows a typical power-law relationship between timescale and level of variability.Region 3 shows timescales at which there is no further increase in variability.The transition between Regions 2 and 3 shows the characteristic timescale of the variability.spectrum.
For a more detailed description of some properties of structure functions, see section 5 of Lakeland & Naylor (2022) and references therein.To avoid spurious results arising from bins with too few points, we require at least 100 pairs of observations in each  bin throughout this paper.

C1 The relationship between the structure function and standard deviation
Consider an infinite, uncorrelated time series  ().This time series will have mean This is equivalent to eq.A4 of Simonetti et al. (1985).

Figure 1 .
Figure 1.Scatter plots showing the correlation between daily mean RVs (orange for the HARPS-N RVs and blue for the inactive-region RVs) and log  ′ HK (left) and unsigned magnetic flux (right).Both scatter plots demonstrate an almost complete removal of any correlation between the RVs and the respective activity indicator after subtracting the magnetic activity RVs.For each data set, we also numerically show the Pearson correlation coefficient.RV data are shown with an offset for clarity.

Figure 2 .
Figure 2. (Top)The power spectral density (PSD) of HARPS-N RVs (orange) and of the inactive-region RVs (blue), calculated in 150 logarithmically-spaced bins in frequency.(Bottom) The ratio of HARPS-N RV PSD to inactive-region RV PSD.Although there is still a noticeable feature in the residual RV PSD at the rotation period, and at one day and its alias, we see lower power at timescales longer than ∼ 10 days when compared to the full HARPS-N RVs.By contrast, the two PSDs are very similar at shorter timescales.This indicates we are successful in removing much of the variability at timescales around or above the solar rotation period, but we leave the shorter timescale variability largely unaffected.

Figure 3 .Figure 4 .
Figure 3. Scatter plot of the mean sunspot number and the standard deviation of inactive-region RV, calculated for non-overlapping 50-day periods from 2015 to 2021.We see no correlation between RV RMS and mean sunspot number.

Figure 5 .FractionFigure 6 .
Figure 5. (Left) The structure functions shown in Fig. 4 of inactive-region RVs.(Right) The structure functions of the corresponding year-long run of magneticactivity RVs.To illustrate the relative level of variability in the inactive-region and magnetic activity RVs, both sets of structure functions are shown on the same scale.The inactive-region RVs consistently show a level of variability around 1 m s −1 .The relatively constant level of variability in the inactive-region RVs contrasts drastically to the large range of variability levels demonstrated in the magnetic activity RVs.Additionally, the level of variability shown by the inactive-region RVs is often much larger than that in the magnetic activity RVs.As in Fig. 4, both sets of structure functions are coloured by the average sunspot number during the calendar year in which the RVs were obtained.

Figure 7 .Figure 8 .
Figure 7. Structure function comparison of inactive-region RVs calculated by subtracting the RV signatures of magnetically-active regions (Section 2.2) from HARPS-N (black) and NEID (red) RVs between 2021-Apr-01 and 2021-Nov-12.Both structure functions demonstrate similar variability characteristics despite the different spectrographs, indicating that the variation seen in the inactive-region RVs is predominately astrophysical in nature, not instrumental.

Fig. 8
Fig.8shows structure functions of each component in the simulation ofMeunier et al. (2019) (pink), the combined simulated time series (red), and the inactive-region RVs (black).The granulation (Gra) and supergranulation (Sgr) components are re-scaled by a constant factor to match the observations; the oscillation (Osc) component is assumed to be negligible (e.g.,Chaplin et al. 2019) and so is not re-scaled.However, we opt to include the oscillation component for completeness.Since granulation and supergranulation dominate at different timescales, we are able to re-scale the simulated RVs to constrain the contribution of each process by matching the power in the inactive-region RVs at short and long timescales, respectively.By combining and re-scaling the simulated RVs in this way, we are able to accurately reproduce both the overall variability amplitude of the inactive-region RVs, and the timescale of the inactive-region RVs.

Fig. 7 of
Fig.7ofSergison et al. (2020), reproduced here as Fig.C1, demonstrates the typical form of a SF.The short timescale plateau shown in Region 1 indicates any intrinsic variability in the signal is

Figure C1 .
Figure C1.An illustration of the typical properties of structure functions, reproduced fromSergison et al. (2020).Region 1 shows white noise present at short timescales, often interpreted as instrumental noise.Region 2 shows a typical power-law relationship between timescale and level of variability.Region 3 shows timescales at which there is no further increase in variability.The transition between Regions 2 and 3 shows the characteristic timescale of the variability.