Shadowing the Sun I: Testing Synthetic Granulation Spectra with NEID Solar Eclipse Observations

Introduction

The first exoplanet orbiting a main-sequence star was discovered in October 1995 using the radial velocity (RV) method ([62]). Since this discovery, the RV method remains a crucial technique for confirming and characterizing candidate exoplanets, and is the only current method that can potentially measure the masses of terrestrial exoplanets orbiting within the Habitable Zone of nearby stars, some of which will be observed by the planned Habitable Worlds Observatory (HWO, Alice et al. [4]). Advancements in present-day spectrographs, such as optical fiber double scramblers (S. Halverson et al. [39]), stable environment control systems (P. Robertson et al. [88]), and laser frequency comb-based wavelength calibration (C.-H. Li et al. 2008), allow for sub-$\mathrm{m\ s^{-1}}$ RV precision. However, an Earth-mass exoplanet orbiting within the Habitable Zone of a solar-mass star produces a Doppler shift of only $8.9\ \mathrm{cm\ s^{-1}}$, which currently remains beyond our reach, not primarily due to instrumentation, but also due to astrophysical stellar variability. Several strategies have been designed to model and mitigate the various forms of stellar variability, which occur on many different timescales. For example, pressure modes (acoustic modes excited by near-surface convection) can be reduced by binning over the pressure mode period (W. J. Chaplin et al. [18]). Rotationally-modulated magnetic activity, such as dark spots, bright plages, and faculae, can be efficiently modeled via Gaussian processes (e.g., V. M. Rajpaul [82]; R. Haywood [42]; M. R. Kosiarek & I. J. M. Crossfield [50]; C. Beard et al. [6]). In addition, several statistical methods, such as StarSim (E. Herrero et al. [44]), starry (R. Luger et al. [58]), and SOAPv4 (E. Cristo et al. [22]), simulate the photometric and spectroscopic impacts of magnetically active regions. However, should such sources of stellar activity be successfully removed, granulation would still hinder the detection of Earth-mass exoplanets (N. Meunier et al. [69]; J. K. Luhn et al. [59]).

Granulation, a convective phenomenon, is the interplay between rising hot plasma called granules and cool intergranular lanes of falling plasma. Granulation can be further decomposed based on magnetic field strength, with magnetohydrodynamic (MHD) studies finding distinct contributions from granules, magnetic and non-magnetic intergranular lanes, and magnetic bright points (H. M. Cegla et al. [14]). Granules on solar-type stars have lifetimes of approximately 10 minutes and range in size from 1300 to 10000 km (D. F. Gray [36]). Granulation RV variability depends on line formation temperature (K. Al Moulla et al. [3, 2]), wavelength region (M. Ellwarth et al. [30]), and ionization and excitation potential (K. Sowmya et al. [99]). The collective motion of granules and intergranular lanes integrated over the solar disk results in 0.1-1.0 $\mathrm{m\ s^{-1}}$ RV shifts (F. Tincon & M. Reuitford 2018). Therefore if left uncorrected, granulation can easily overshadow the Doppler signatures of Earth-mass planets, leaving these exoplanets undetectable.

The effects of granulation in perturbing the position and shape of spectral lines have been widely discussed in the literature. Following the nomenclature of D. F. Gray [36], these perturbations are often referred to as first, second, and third signatures. The first signature refers to the supra-thermal broadening of spectral lines that makes the lines wider and shallower. The second signature refers to the asymmetries of spectral lines, and the third signature refers to the net convective blueshift of spectral lines. The latter two signatures are a direct result of the granules being hotter, larger and therefore contributing more light than the intergranular lanes (D. Dravins [24]). Center-to-limb variations (CLV) have been measured in disk-resolved solar observations (e.g., J. Löhner-Böttcher et al. [57]; F. Stief et al. [100]; J. Löhner-Böttcher et al. [56]). At disk center, spectral line bisectors manifest in a “C” shape, whereas towards the solar limb they resemble more of a “\” shape. In addition, the net convective blueshift declines towards the limb, with some lines demonstrating slight redshifts near-limb (J. Halm [38]; M. Ellwarth et al. [30]). These CLVs are a direct result of different viewing angles across the solar disk. At disk center, granules are viewed deeper into the photosphere, whereas towards the limb the optical path length increases so granules are viewed higher in the photosphere. As a result, different temperature contrasts, velocity flows, and therefore line profiles are observed across the solar disk.

It has been claimed that the RV effects of granulation can be reduced via observation strategies such as averaging multiple observations spaced by hours in a single night (X. Dumusque et al. [29]; N. Meunier et al. [68]). In addition, 3D MHD simulations combined with radiative transfer have been used to quantify granulation-driven line perturbations (H. M. Cegla et al. [15, 17, 16]; G. Frame et al. [33]). A strong correlation has been found between spectral line bisector statistics and the granulation-induced RV shift (H. M. Cegla et al. [16]), suggesting that precise measurements of line bisectors could trace, and therefore help mitigate, granulation-driven RV shifts. However, such MHD simulations are based on a single simulated line, whereas real stellar spectra have thousands of lines. To address this limitation, M. L. Palumbo et al. [75, 72] developed the GRanulation And Spectrum Simulator (GRASS, available on Zenodo, M. L. Palumbo III et al. [76]) that creates synthetic spectra with time resolved granulation-driven variability. Unlike the MHD simulations, GRASS uses high spectral resolution disk-resolved Laser Absolute Reference Spectrograph (LARS, J. Löhner-Böttcher et al. 2017) observations of the Sun to create synthetic disk-integrated spectra. M. L. Palumbo et al. [72] found that use of bisector statistics can only remove 25–35% of granulation noise. Hence, detection of Earth-mass exoplanets requires accurately accounting for granulation-driven RV variability. GRASS successfully reproduces the disk-integrated, time-averaged line profiles of the Institut Für Astrophysik, Göttingen (IAG) Solar Atlas ([84]), but it has not, prior to our work, been tested against highly precise RV observations.

While granulation-induced CLV of line profiles can be directly observed on the Sun due to its spatially resolved surface, such information can only be indirectly probed for other stars using planetary transits ([26, 25]). Using archival data of HD 209458 and HD 189733 from the UVES ([23]) and HARPS ([78]) spectrometers, D. Dravins et al. [27, 25] extract center-to-limb observations of Fe I lines with remarkable accuracy to predictions from hydrodynamic 3D models. Yet, the Sun, with its spatially resolved surface and precisely known dynamics, provides an ideal opportunity to test the accuracy of RV models that account for the granulation signatures. During a solar eclipse, solar RVs show the Rossiter-McLaughlin (RM, [91, 64]) effect, which is characterized by blue/red-shifted deviations in the RV curve that scale with the fraction of surface occulted, the star’s rotational velocity, and the spin-orbit alignment between the stellar spin axis and the planet’s orbital plane ([35]). Estimates of spin-orbit alignment angles inform planet migration theories, as the degree of misalignment serves as a proxy for the degree of violent interactions during migration. However, ignoring granulation-driven convective blueshift in RM observations leads to RV deviations, which impact estimates of spin-orbit alignment angles ([96, 12 b]). A solar eclipse provides an unique opportunity to investigate the influence of the granulation signatures on RVs, as the large occultation from the Moon results in a large RM effect and drastically altered spectral line shapes.

The first study of this manner was performed by A. Reiners et al. [83] for the 2015 March 20 solar eclipse, which was observed with the IAG Fourier Transform Spectrograph ([52]). With these observations, A. Reiners et al. [83] tested various RV models with distinct treatments of convective blueshift to highlight the influence of convective blueshift on the RM curve. They found that inclusion of convective blueshift in their RV models reduced the peak velocity residual by roughly 50% regardless of how the RVs are computed. However, even with convective blueshift accounted for, there still remained a significant data-model discrepancy with peak residuals up to ~25 m s$^{-1}$. They suspected their models were limited by lack of understanding line profile CLV, for which they adopted a sigmoid function with arbitrarily prescribed coefficients. Here we follow up on the [83] study by modeling the 2023 October 14 solar eclipse using GRASS to determine whether the empirically based treatment of granulation therein is able to better replicate solar eclipse RVs.

Spectral observations of the Sun as a star were obtained during the 2023 October 14 solar eclipse using the NEID solar feed ([54]) and spectrograph ([95] (S. Halverson et al. 2016)) at Kitt Peak National Observatory. We model the solar RVs using parameterizations for limb darkening, differential extinction, differential rotation, and the granulation information from GRASS. The data used are described in § 2. The modifications to GRASS and the RV models used are discussed in § 3, followed by a discussion of our results in § 4.

Eclipse Observations

Solar eclipse RVs are sensitive to granulation signatures and so allow us to test the accuracy of granulation-informed RV models. However, such RV observations require high spectral resolution in order to best capture the granulation signatures. We designed an experiment to obtain such data with the NEID spectrometer, and describe the obtained observations below.

Instrument

The 2023 October 14 solar eclipse was observed with the NEID spectrograph ([95]), which is attached to the WIYN 3.5 m telescope³ at Kitt Peak National Observatory. NEID is a state-of-the-art high-resolution spectrometer ($R \sim115{,}000$, 380–930 nm) that obtains spectroscopic data of nearby stars. NEID is highly stabilized for an instrumental RV precision of better than 27 cm s$^{-1}$ (S. Halverson et al. 2016) and is wavelength-calibrated using a laser-frequency comb and Fabry-Pérot etalon. This analysis specifically makes use of the NEID solar feed ([54]) for Sun-as-a-star observations. The NEID solar feed assembly includes a solar tracker to maintain pointing precision and a pyrheliometer as a cloud sensor. More details of the NEID solar feed can be found in [54]. Sun-as-a-star NEID observations over the course of 3.5 years when modeled with SCALPELS [21] reveal an on-sky RV precision of 37 cm s$^{-1}$ — suggesting that the capability of detecting Earth-mass exoplanets orbiting solar-mass stars may be coming into reach [32].

Observations

At its maximum, the 2023 October 14 solar eclipse covered a total of 76% of the solar disk when viewed from Kitt Peak. The eclipse began at 15:11 UTC (8:11 AM local time) at an airmass of 2.86 and ended at 18:03 UTC (11:03 AM local time) at an airmass of 1.38. The maximum occurred at 16:31 UTC (9:31 AM local time). The eclipse configurations are shown in Figure 1. NEID solar observing automatically triggers at 9:30 AM local time, which is after the start time of the eclipse. Therefore, observations were manually triggered starting at 8:04 AM local time immediately after a shortened (high-resolution mode only) calibration sequence following the on-sky observing from the previous night. Due to clouds at the beginning and end of the observations, our analysis data begin at 15:26 UTC (8:26 AM local time) and ends at 18:24 UTC (11:24 AM local time) for a total of 130 observations. The exposure time of each observation followed the standard NEID solar observation sequence of 55 seconds exposure and a read out time of 27 seconds for a total observing cadence of 82 seconds. We note that daily at 10 AM local time a liquid nitrogen fill occurs that does induce a change in the intrinsic velocity drift of the NEID spectrometer [54]. However, this drift is tracked and removed by a simultaneous calibration source observed with the solar data. The measured signal-to-noise ratio (SNR) per resolving element of the observed spectra for this work’s case study line, 5434.5 Å (air wavelength, vacuum wavelength of 5436.0 Å), is ∼378 outside of the eclipse and drops to ∼161 during the time of maximum eclipse. During normal operations an automated routine periodically ensures that the output image of the fiber carrying the solar light is well centered on the NEID input fiber at the telescope port. However, during this experiment we opted not to engage this routine to prevent issues due to the rapidly changing intensity of the Sun from the eclipse.

Measuring Radial Velocities

The NEID data reduction pipeline⁴ (DRP) RVs are calculated by computing a cross-correlation function (CCF) between each observed spectrum and publicly available ESPRESSO [79] masks. Since the NEID DRP RVs are essentially an average RV computed from a large array of lines simultaneously, they do not allow for investigation of line-by-line variability. Therefore, we separately compute line-by-line RVs from the NEID observed spectra for a selected set of lines between 525–630 nm. We use RvSpecTM [31], which provides multiple algorithms for estimating RVs from a CCF, to extract a specific spectral line from the NEID spectra and compute line-by-line RVs following the CCF algorithm from F. Pepe et al. [77]. The CCF is computed by shifting a template mask along the spectra at varying velocities. Our template masks are top-hat functions centered at the rest-frame wavelength of each individual spectral line with the height determined by the spectral line depth. Following F. Pepe et al. [77], we determine the RV from the CCF by fitting a Gaussian around the CCF minimum (defined using 75% of the CCF full width at half maximum) using the Levenberg-Marquardt algorithm. The mean of the Gaussian fit represents the apparent Doppler velocity of the spectral line and therefore is our line-by-line RV. The same procedure is used to compute the simulated RVs described in Section § 3.2.

Eclipse Simulations

The spectral signatures of granulation are modeled using a modified version of GRASS v2.0. GRASS creates synthetic spectra with granulation-driven variability from spatially and temporally resolved observations of 22 solar absorption lines [75, 72]. Spectroscopic parameters for these spectral lines can be found in Table 1 in M. L. Palumbo et al. [72]. Each input spectral line used by GRASS v2.0 was observed by LARS ($R \sim700{,}000$, J. Löhner-Böttcher et al. [55]) at up to eleven positions along the North, South, East and West axes from disk center with a field of view up to 40x30 arcseconds and a pointing accuracy of $\sim1$ arcsecond [57, 100, 56]. Such spatially resolved observations are critical when attempting to capture CLV. The observations were done in $\sim20$ min sequences with each exposure having a cadence of 1.5 s then binned into ten exposures (15 s). By using these data, GRASS captures temporal variability in line shape and convective blueshift. Furthermore, the statistical behavior of granulation used in GRASS has been shown to closely match state-of-the-art hydrodynamic simulations [33]. We refer the reader to M. L. Palumbo et al. [75] and M. L. Palumbo et al. [72] for specific details on the spatial interpolation and integration procedure used by GRASS.

Spectral Line (Å)	$a$	$b$
Fe I 5250.2	1.625	-2.294
Fe I 5250.6	1.421	-2.193
Fe I 5379	1.519	-2.180
Ti II 5381	0.497	-0.938
Fe I 5382	2.075	-3.088
Fe I 5383	0.978	-1.646
Mn I 5432	0.642	-0.945
Fe I 5432	1.616	-2.355
Fe I 5434	0.637	-1.221
Ni I 5435	1.431	-1.997
Fe I 5436.3	1.893	-2.697
Fe I 5436.6	1.379	-2.047
Fe I 5576	0.783	-1.356
Ni I 5578	2.003	-2.737
Fe II 6149	1.442	-2.139
Fe I 6151	0.919	-1.390
Ca I 6169.0	0.638	-1.089
Ca I 6169.5	0.721	-1.180
Fe I 6170	0.653	-1.067
Fe I 6173	1.410	-2.065
Fe I 6301	0.660	-1.230
Fe I 6302	1.548	-2.336

GRASS Repurposed for Solar Observations

GRASS was intentionally created for stars and makes extensive use of large-distance approximations that assume the observing geometry does not drastically change with time. Therefore, as written, GRASS cannot be used to model the resolved Sun from the perspective of the Earth orbiting it. To use GRASS to interpret solar observations, we modify the geometry (§ 3.1.1), explore various new intensity models (§ 3.1.2, § 3.1.3) and incorporate sunspots (§ 3.1.4). These modifications are available on the GRASS GitHub⁵ under the Julia submodule GRASS.Eclipse, referred to as GRASS-E herein.

Geometry Modifications

The geometry used in GRASS-E generates synthetic granulation-driven spectra of the eclipsed Sun. The state vectors of the Earth, Moon and Sun during the eclipse are acquired using NASA’s Navigation and Ancillary Information Facility SPICE toolkit (C. H. Acton [1]), which provides highly accurate ephemerides of Solar System objects with an accuracy of a few cm s$^{-1}$. Such state vectors are corrected for light travel time within the SPICE toolkit by iteratively solving the light time equation until a solution converges. Both the Moon and Sun are assumed to be perfect spheres; whereas, Earth’s oblateness, defined as the fractional difference between the equatorial and polar radii, is accounted for in calculating the position of the observatory. We maintain the same solar surface gridding algorithm described in Appendix A.§ 2 in M. L. Palumbo et al. [72]. For each cell on the solar grid, the differential rotation velocity is determined using the following functional form:

where the angular rotational velocity $\omega$ is a function of solar latitude $l$ [45]. We use the spectroscopic coefficients from H. B. Snodgrass & R. K. Ulrich [97] ($A = 14.113^\circ\ \mathrm{day}^{-1}$, $B = -1.698^\circ\ \mathrm{day}^{-1}$, $C = -2.346^\circ\ \mathrm{day}^{-1}$). The differential rotation velocity is projected onto the line of sight from each individual solar cell to the observer, which is specified by user-given latitude, longitude and altitude⁶. In order for our models to best mimic the involved dynamics under which the observations were obtained, we include the velocity from the Earth’s rotation and orbital motion towards each individual solar cell. Hence, the final velocity of each solar cell is the sum of the projected differential rotation velocity and the projected rotation and orbital velocity of the Earth towards the solar cell. By including the barycentric motion of the observer on Earth in our models, we can apply a consistent barycentric correction to both the NEID-observed and simulated RVs using barycentry (S. K. Kandoia & J. J. Wright 2018).

The key solar-related modification in GRASS-E is that for each epoch we must update the SPICE query of the Earth, Moon, and Sun state vectors, which are used to determine (a) whether a solar grid cell is eclipsed or not, and (b) the time-evolving projected cell velocity and area grid. Based on this new geometry, GRASS-E generates time-series absorption line profiles that are shifted by each solar cell’s projected velocity, then integrated (weighted by projected area and intensity) to produce synthetic flux spectra.

Limb Darkening Modifications

In addition to geometry modifications, we also modify the intensity model in GRASS-E. The default limb darkening model used by GRASS is the wavelength-independent, quadratic limb darkening law (Equation (2)) from Z. Kopal (1950) with fixed coefficients of $u_1 = 0.4$ and $u_2 = 0.26$ [98].

However, commonly used limb darkening laws, such as the quadratic law, fail to reproduce true center-to-limb intensity variations [71]. N. M. Kostogryz et al. [51] show that magnetic fields have a direct impact on limb darkening and can help explain the discrepancies between stellar atmosphere models and observations. N. M. Kostogryz et al. [51] use MURaM [104], a three-dimensional radiative MHD code, to simulate limb darkening values at $\mu= 0.1$–$1.0$ (with a $\Delta\mu$ of $0.1$) for the following magnetic field strengths: hydrodynamic (i.e. no field), small-scale turbulent dynamo (SSD), 100 G, 200 G, and 300 G. The SSD generates a baseline level of magnetic flux (70 G) across the entire stellar surface as a result of convective motion distorting magnetic field lines [85].

The four-parameter limb darkening law is considered the most accurate [81] and is as follows:

We fit the quadratic Equation (2) and four-parameter Equation (3) limb darkening laws to the SSD and 300 G MURaM simulated limb darkening. Figure 2 shows the simulated limb darkening values as a function of $\mu$ at SSD and 300 G taken from N. M. Kostogryz et al. [51] for 5428.5 Å, which is the closest to our case study line of Fe I 5434 Å (air wavelength). The best fit quadratic and four-parameter laws are shown in the dashed and dot-dashed lines, respectively, in the corresponding magnetic field strength color. In addition, we include in the solid black line the V-band four-parameter law as determined by A. Claret [19]. The A. Claret [19] limb darkening prescription is based on plane-parallel stellar atmosphere models, which fail to capture the impact of magnetic fields on limb darkening. As shown in Figure 2, the A. Claret [19] limb darkening prescription aligns with the SSD simulated limb darkening, which differs from the 300 G simulated limb darkening. This is consistent with the N. M. Kostogryz et al. [51] finding that stellar atmosphere models are only able to reproduce limb darkening when little to no magnetic fields are physically involved. N. M. Kostogryz et al. [51] found that for the quiet Sun the SSD limb darkening is the most appropriate. Since GRASS–E is informed by quiet-Sun LARS observations, we decide to use the simulated SSD limb darkening values to determine our limb darkening model. However, we do caution that the solar eclipse occurred when the Sun was near a maximum activity phase, so using SSD limb darkening and quiet Sun center-to-limb observations to inform CLV may not capture all magnetic effects.

Given its flexibility near-limb, we use the four-parameter law fit of the simulated SSD limb darkening over the quadratic law. We find that although the A. Claret [19] V-band limb darkening prescription aligns well with the Fe I 5434 Å simulated SSD limb darkening values (Figure 2), towards the redward end of the V-band the prescription starts to fail. For example, the reddest template line from GRASS–E, Fe I 6302 Å, has a growing deviation towards smaller $\mu$ between the A. Claret [19] V-band limb darkening prescription and the N. M. Kostogryz et al. [51] simulated SSD limb darkening. This suggests that line-by-line limb darkening is preferred over the V-band limb darkening prescription. Hence, our intensity models use the four-parameter limb darkening law (Equation (3)) with the line-by-line fitted coefficients from the simulated SSD limb darkening.

Modeling Extinction

In addition to updating the limb darkening model used in GRASS–E, we also include differential extinction in our intensity model. Since the eclipse occurred early in the morning while the Sun was at a high airmass, we expect extinction to be strong. In addition, due to the Sun’s large angular size (0.5 degrees), extinction will differ significantly across the solar disk. Under cloud-free conditions, extinction follows an exponential law as a function of airmass ([20]). The corresponding extinction coefficients are seeing-dependent, and so are generally unknown. To determine the wavelength-dependent extinction coefficients ($\alpha(\lambda)$), we perform a non-linear least-squares fit (via the Levenberg-Marquardt algorithm) of the following:

where $I_{\mathrm{NEID}}$ is the observed NEID intensity, $I_{\mathrm{LD}}$ is the four-parameter limb darkening law (Equation (3)) as a function of wavelength and time, and $z$ is the zenith angle. Figure 3 highlights the impact of including extinction in our intensity model: including the best-fit extinction reduced the intensity residuals by nearly 10x. We confirm that our estimated wavelength-dependent extinction coefficients are reasonable and similar in value to the computed extinction coefficients by E. B. Ford et al. [32] for solar NEID observations. We also allow the extinction coefficients to slightly vary with time to reflect extinction changes during the eclipse and find no significant improvements in our models. Hence, we proceed with constant-in-time wavelength-dependent extinction coefficients fitted from the extinction law.

Modeling Sunspots

Magnetically active regions, such as spots, degrade RV precision in two distinct ways: (1) the flux effect and (2) the suppression of convective blueshift effect. The flux effect arises from the temperature contrast between an active region and the surrounding quiet photosphere (X. Dumusque et al. [28]). These active regions emit different amounts of light than the surrounding quiet photosphere, breaking the symmetry between the rotationally red- and blue-shifted hemispheres and introduce asymmetries in disk-integrated spectral line profiles and quasi-sinusoidal RV signals. The suppression of convective blueshift effect results from magnetic fields within an active region suppressing convective motion (F. Cavallini et al. [11]). In the Sun, this suppression can reach $\sim300\ \mathrm{m\ s^{-1}}$ in spots, making magnetically active regions appear redshifted relative to the average photosphere (N. Meunier et al. [65]). Images taken by the Helioseismic and Magnetic Imager (HMI, J. Schou et al. [94]) onboard the Solar Dynamics Observatory (SDO, W. D. Pesnell et al. (2012)) reveal the presence of sunspots during the 2023 October 14 solar eclipse. Figure 4 shows the downloaded HMI continuum image (exposure time of 720 s) at 16:35 UTC, near maximum eclipse. All sunspots are shown in Figure 4 with the largest three sunspots highlighted in red. The HMI image is used to extract the coordinates, effective radius, and mean intensity contrast for each of the largest three sunspots. Given this information, we are able to include both the flux and convective blueshift effects in our RV model. We assume the spots were circular with no intensity gradient and not changing in size or position. Based on the spots’ coordinates and angular sizes, cells in our solar grid that are found inside a spot have their intensity reduced by the determined intensity contrast and their convective blueshift suppressed to zero. We find that accounting for the flux and convective blueshift effects only improves our RV residuals RMS by $< 1.0\ \mathrm{m\ s^{-1}}$. We acknowledge that other sources of magnetically active regions, such as plages, can impact RV precision. However, given the minuscule improvement in RV RMS when including spots, we do not explore other sources of magnetically active regions.

Radial Velocity Models

We simulate the 2023 October 14 solar eclipse RVs at the mid-exposure time of the NEID observations using three distinct models. We find that during the 55 second NEID exposure, solar flux from the NEID exposure meter varies linearly throughout the eclipse. In addition, we find that simulating the RVs at fine samples (down to one second) within the 55 second exposures and then computing a weighted average RV for an individual 55 second exposure does not significantly improve our models. Each model follows the same geometry and intensity computations described within § 3.1.

• Model I computes RVs as the weighted average projected velocity by weighting each solar cell’s projected velocity with its corresponding intensity and projected area. This model does not account for any granulation signatures.

• Model II computes RVs using the CCF-algorithm described in § 2.3 on simulated GRASS-E spectra that exclude all time-varying granulation signatures (convective blueshift jitter, line asymmetry).

• Model III also computes CCF RVs but from GRASS-E simulated granulation-driven spectra.

Model III is the only model that runs the original GRASS granulation-driven spectra simulator. While model II is run through GRASS-E for consistency in surface gridding, spectra simulation, and CCF RV computation, the generated spectra from model II do not account for any granulation signatures. In other words, model II grids

the Sun as model III does, but rather than assigning a line with granulation signatures to each cell, model II instead assigns a symmetric non-blueshifted line that is static with time. From there, model II performs a spatial integration to get a disk-integrated spectrum (no granulation signatures) from where a CCF RV is computed. Since the only difference between models II and III is the inclusion of the granulation signatures, we can use these models to quantify the impact of granulation and validate GRASS-E’s granulation modeling. We note that the observations informing GRASS-E’s granulation model are not contemporaneous with the solar eclipse. Therefore, while model III captures the statistical behavior of granulation, the asynchronous nature of the data acts as uncorrelated noise at short timescales. However, the intrinsic granulation-induced RV variability is on the order of $\sim$45–75 cm s$^{-1}$ ([72]), which is small compared to our photon noise limit of $\sim$2–3 m s$^{-1}$ for individual lines. Hence, the temporal mismatch contributes minimally to the overall RMS and does not drive the larger structures seen in the residuals.

In order to directly compare the RVs from models II and III to the NEID line-by-line RVs, the GRASS-E simulated spectra must be degraded from the resolving power set by the LARS input data ($R \sim700{,}000$) down to NEID’s resolving power of $R \sim115{,}000$. We test three line spread functions (LSF): a Gaussian, a symmetric NEID LSF⁷, and an asymmetric NEID LSF (Joshi in prep.). The NEID LSF is a non-Gaussian kernel and is modeled as a convolution of: the pixel response function, a 1D integrated circular optical fiber image, and a Gaussian with a standard deviation of 0.6 pixels. However, significant residuals are observed when comparing this symmetric LSF model with the laser frequency comb (LFC) spectrum. Joshi in prep. account for both the symmetric and asymmetric deviations from the observed LFC spectrum by perturbing the model LSF with orthogonal polynomials up to degree 8 in order to model the true NEID LSF. We find that the two NEID LSF models perform comparably and tend to outperform the Gaussian LSF. Therefore, we use the asymmetric NEID LSF to degrade the GRASS-E spectra to NEID’s resolving power.

All three models are computed separately for each of the 22 GRASS-E template lines in order to incorporate wavelength-dependent limb darkening, extinction, and granulation signatures. Line-by-line models allow us to evaluate the variability of each line (see Section § 4.4). We compare the observed and simulated line profiles of the 22 GRASS-E template lines and confirm that, within 2%, GRASS-E successfully reproduces the line profiles from the NEID spectra.

Results

The RM curve of the 2023 October 14 solar eclipse was simulated using the three models described in Section § 3.2. We here use the Fe I 5434 Å spectral line as an initial case study.

Figure 5 shows the NEID line-by-line RVs for Fe I 5434 Å along with the simulated RVs from each model. Models I and II show similar RM curves with an offset towards the blue (negative velocities) during the eclipse. The simulated RM curves of models I and II are significantly similar based on their residuals following the same shape (peaking at $\sim75$ m s$^{-1}$) and their binned RMS both at $\sim30$ m s$^{-1}$. To average over the $\sim5.4$ min solar p-mode period, we bin RVs into sets of 4 exposures (328 s, 5.4 min). Neither model I nor II accounts for any granulation signatures. The key difference between models I and II is how the RVs are computed: model I computes an average weighted projected velocity across the solar disk, while model II computes a CCF RV from simulated spectra. Given the similarities between the RM curves of models I and II, how the RVs are computed does not seem to significantly influence the shape of the simulated RM curve. Model III, which computes a CCF RV from GRASS-E simulated granulation-driven spectra, performs better than models I and II by a factor of more than 2 in RV RMS. Figure 5 highlights that model III does experience a similar offset towards the blue like the other two models. The key difference between models II and III (both CCF RVs from simulated disk-integrated spectra) is that model III accounts for convective blueshift and line asymmetries, whereas model II does not account for any granulation signatures. As a result, model III’s better performance can be attributed to its inclusion of granulation signatures. The results from these models emphasize the known need to account for granulation signatures in solar eclipse RV models ([83]).

Following F. Bouchy et al. [10], we estimate the photon noise for Fe I 5434 Å to be: 2.75 m s$^{-1}$ for all considered observations, 3.01 m s$^{-1}$ during the eclipse, and 2.05 m s$^{-1}$ outside of the eclipse. Outside of the eclipse is defined as 18:05 UTC up to 18:24 UTC, where clouds start to affect the RVs. We obtained a total of 15 out-of-eclipse observations, each with a cadence of 82 s. Our full data binned RMS (depicted in Figure 5) do not come close to the photon noise limit, suggesting that there may be missing or incorrect astrophysics in our RV models. We explore some potential sources of errors in the following sections.

Pressure Modes

Pressure modes, or p-modes, are oscillations on a stellar surface that arise due to traveling sound waves that are generated by pressure fluctuations from convection in the interior of stars. Although p-modes introduce $\sim1$ m s$^{-1}$ RV variability, they can be averaged down to sub-10 cm s$^{-1}$ level by binning RVs according to the principal p-mode period ([18]), which is $\sim$5.4 minutes for the Sun. Hence, we bin the entire eclipse RVs into sets of 4 exposures ($\sim$5.4 min). We find that binning reduces the RMS very slightly, suggesting p-modes are involved. However, the improvement is very small compared to the residuals. We test how RMS changes as a function of different bin sizes and find that additional binning does not contribute to a significant reduction in RMS.

Hence, p-modes alone cannot explain the data-model discrepancy during the eclipse. We do acknowledge that during the eclipse partial blocking of the solar disk occurs and results in the symmetry breaking of the p-mode averaging; therefore, it is unclear whether averaging over a 5.4 min period continues to successfully remove p-mode variability almost entirely, as expected in precision RV surveys [18, 37]. For consistency with how a normal star would be observed in such surveys, we present all our RV RMS with this 5.4 minute binning, though our plots show the individual 55 second RV observations.

Limitations of GRASS-E

The GRASS-E treatment of granulation signatures is informed by the center-to-limb observations of 22 solar absorption lines taken by LARS [57, 100, 56]. Each absorption line was observed at up to eleven heliocentric positions along four radial axes (North, South, East and West) with sampling of $\mu= [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.85, 0.9, 0.95, 1.0]$, but for some lines only up to $\mu$ of 0.3. For each LARS observation, GRASS-E computes and stores the line convective blueshift, along with the line bisector and width, as a function of line depth. This information is then used to generate synthetic disk-resolved line profiles. Although the LARS observations used by GRASS-E have a pointing accuracy of $\sim$1 arcsecond and a high resolving power of $\sim$700,000, the field of view of the observations covers up to 40x30 arcseconds in order to average over the approximate size of a supergranule. Such large field of view introduces uncertainties of the observed heliocentric position ($\mu$), which worsen towards the limb. As a result, the measured CLV of bisector asymmetry and convective blueshift used by GRASS-E are for a superposition of $\mu$ regions, especially towards the limb. The LARS-measured center-to-limb convective blueshift for Fe I 5434 Å is highlighted in Figure 6.

GRASS-E takes a nearest-neighbor approach to determine the convective blueshift and line asymmetry at a given solar position. In other words, to create the synthetic disk-resolved line profiles, GRASS-E determines the granulation signatures at a given solar position based on the nearest solar axis and $\mu$ value from the input data sampling. As a result, the line bisector and convective blueshift CLV used by GRASS-E resembles a piecewise constant function, which is not reflective of observed CLV. Observed convective blueshift CLVs tend to resemble more of a polynomial shape as a function of $\mu$, with maximum convective blueshift occurring off-center at $\mu$ values between 0.5 and 0.9 [56, 30]. During an eclipse, where partial blocking of the solar disk occurs, problems may arise from using a piecewise constant function for CLV. As the Moon crosses the solar disk during an eclipse, specific $\mu$ regions will be probed where the use of nearest-neighbor line information may be inadequate. In addition, near maximum eclipse, light is only collected from near-limb regions of which GRASS-E only has input observations up to $\mu$ of 0.2 and for some lines only up to $\mu$ of 0.3 (as is the case for Fe I 5434 Å). Therefore, near maximum eclipse GRASS-E’s nearest-neighbor approach is likely unable to accurately reproduce near-limb (past the smallest observed $\mu$) line changes. This may explain the residual feature at maximum eclipse shown in Figure 5. Therefore, when modeling a solar eclipse, it may be crucial to have center-to-limb observations at more finely sampled $\mu$ values. This is especially true near the limb, where a nearest-neighbor approach would fail to capture any extremes in the CLV, such as near-limb redshifts.

Optimizing Convective Blueshift

As an alternative to GRASS-E’s granulation model, we fit a second-order polynomial for convective blueshift as a function of $\mu$. The polynomial form allows for continuous sampling of $\mu$ from 1.0 to 0.0, which is currently not captured by GRASS-E. The second-order polynomial coefficients are determined by optimizing the residual RMS of model I to the NEID observations using an ordinary least squares regression. The regression fits for the shape-dependent coefficients ($a$ and $b$) along with a constant velocity offset. Since we normalize both the NEID-observed and simulated RVs before comparing them, the constant offset cancels out and does not affect the simulated RVs. Similarly, we omit the $+c$ term in the second-order polynomial and only apply the shape-dependent coefficients, which capture the CLV shape. For each cell on the solar grid, the convective blueshift polynomial value is added to the projected velocity of the cell as a function of the cell’s $\mu$ value. This convective blueshift is static with time, unlike the convective blueshift in GRASS-E, model III.

Similar to model I, model II also does not account for any granulation signatures; therefore, model II, which computes CCF RVs unlike model I, can also be used to determine the second-order polynomial coefficients. Using the Nelder-Mead algorithm (an iterative downhill optimization method, [34]), we optimize the residual RMS of model II to the NEID observations to determine the second-order polynomial coefficients ($a$ and $b$). We find that the optimized coefficients using model II agree (within error bars) with the optimized coefficients found using model I. In addition, the residual RMS of models I-CB and II-CB remain consistent within $< 0.05$ m s$^{-1}$ when using either set of optimized coefficients. This is not surprising since the only difference between models I and II is how they compute the RVs, and this does not impact the results as discussed in § 4. We proceed with the optimized coefficients found using model I since ordinary least squares regression requires less runtime for optimizing second-order polynomials for multiple lines as done in § 4.4.

Models I and II, now models I-CB and II-CB, were rerun with the optimized polynomial for convective blueshift. Figure 7 shows the NEID line-by-line RVs for Fe I 5434 Å along with the simulated RVs from models I-CB and II-CB with the optimized convective blueshift polynomial for Fe I 5434 Å (polynomial shown in Figure 6). Figure 7 highlights that models I-CB (RMS of 3.51 m s$^{-1}$) and II-CB (RMS of 2.9 m s$^{-1}$) have significant improvements in RMS compared to the original models I (RMS of 29.88 m s$^{-1}$) and II (RMS of 29.95 m s$^{-1}$) that exclude convective blueshift. In addition, models I-CB and II-CB outperform model III (RMS of 13.47 m s$^{-1}$), which uses the GRASS-E granulation model, by a factor of more than 3 in RV RMS. While model III utilizes an empirical granulation model, it is unable to reach the level of precision achieved by models I-CB and II-CB. This discrepancy highlights that the current GRASS-E

granulation model does not fully capture the complex CLV required to accurately reproduce the solar eclipse observations. Furthermore, the inclusion of an optimized convective blueshift polynomial allows models I-CB and II-CB to reach the estimated $\sim3$ m s$^{-1}$ photon noise of the entire eclipse for the Fe I 5434 Å line.

Line-by-Line Variability

The Fe I 5434 Å line provides a good case study to investigate the different RV models; however, line-by-line variability is known to exist due to line varying excitation and ionization potentials, which influence the line formation location in the atmosphere. Since individual line formation probe unique ranges of heights in the solar atmosphere, different line lists will generate different measured RVs [3, 2]. To investigate potential line-by-line variability, we run our models for all 22 GRASS-E template lines. We proceed with model I-CB that includes the optimized convective blueshift polynomial, model II-CB that also includes the optimized convective blueshift polynomial, and model III that keeps the GRASS-E’s convective blueshift treatment. Table 1 shows the optimized coefficients of the line-by-line convective blueshift polynomials. We only report the shape-dependent coefficients ($a$ and $b$) since the simulated RVs are insensitive to the $+c$ term due to the normalization of the RVs. We confirm that each of these optimized convective blueshift polynomials follow the expected trend of reduced convective blueshift towards the limb. Figure 8 highlights the four-exposure binned RMS of each model and each line. Across all 22 lines, model III performs the weakest — emphasizing that an RV model with a simple convective blueshift polynomial outperforms GRASS-E’s advanced empirical granulation treatment that accounts for convective blueshift and line asymmetry. For the majority of the spectral lines, model II-CB performs better than model I-CB, suggesting that once convective blueshift is accounted for, the method of RV computation can now impact the results. Yet, the majority of the models fail to reach the average photon noise limit for the considered observations (depicted in black x’s). In addition, the RV RMS varies across the different spectral lines. We explore correlations between the RV RMS and various line parameters, but find no obvious correlations with any single parameter, which is consistent with the findings in M. L. Palumbo et al. [72]. The variability of model III could be driven instead by differences in the LARS data used in GRASS-E, such as varying $\mu$ coverage from center to limb. The combined model II-CB weighted RMS of all 22 GRASS-E template lines (weighted by inverse squared RMS) is $\sim5.2$ m s$^{-1}$ over the RM signal.

The Impact of Charge-Transfer Inefficiency

During readout, charge-coupled device (CCD) detectors may experience errors when transferring recorded charge between pixels, an effect known as charge transfer inefficiency (CTI). Since CTI is dependent on the total collected flux, precision RV surveys set exposure SNR triggers for all observations of a given star [9, 37]. Therefore, our solar eclipse observations, with 76% reduction in total collected flux at maximum eclipse and with fixed 55 second exposures, may be sensitive to CTI uncertainties. We consider how the data-model I-CB RV deviations vary as a function of line-by-line SNR and find no trend with SNR. Even the drastic 76% drop in total flux during the eclipse does not induce significant CTI-related uncertainties, which is not surprising given NEID’s 5 cm s$^{-1}$ CTI induced RV error budget (S. Halverson et al. 2016) and the consistency with on-sky tests during NEID commissioning (NEID Team, private communication).

Optimizing A Non-Radial Velocity Field

The fact that the majority of the spectral lines we consider do not reach their photon noise limit leads us to believe that additional sources of astrophysical variability are missing from our models. We suspect longer length-scale velocity fields, such as meridional flows (e.g., [60, 67]) and supergranulation (e.g., [87, 66]).

Meridional Flows

We begin exploring additional velocity fields by incorporating to our RV models a ground-based version of the meridional flows velocity field model used in M. L. Palumbo et al. [73] (code available on GitHub, M. L. Palumbo et al. [74]). Meridional flows are a direct result of turbulent-driven transports of angular momentum on stars’ surfaces. Related plasma flows occur symmetrically poleward at velocities up to 40 m s$^{-1}$, antisymmetric from the equator [48]. Such velocity fields are largely ignored in RV models, despite being shown to introduce 1-2 m s$^{-1}$ RV variations [67]. Our meridional flows model is based on S. G. Kashyap & S. M. Hanasoge [48], which utilize even-degree ($l = 2$, $l = 4$) poloidal components of vector spherical harmonic expansion. We refer the reader to § 2 and Appendix A of S. G. Kashyap & S. M. Hanasoge [48] for specific details on the parameterization of meridional flows. The meridional flows model is a linear combination of the $l = 2$, $l = 4$ poloidal components. We determine the optimized meridional flows velocity field by optimizing the coefficients of the $l = 2$, $l = 4$ linear combination. Such optimization is completed similarly to how we determine the optimized convective blueshift polynomial in Section § 4.3, by optimizing the residual RMS of model I-CB, with the $l = 2$, $l = 4$ Legendre polynomials for meridional flows, to the NEID observations using an ordinary least squares regression. Much like the optimized convective blueshift, the meridional flows velocity is added to the projected velocity of the cell but as a function of cell solar latitude and longitude. The poloidal axis is fixed to the Sun’s rotation axis.

The optimized non-radial velocity field is consistent with the expected meridional flows velocity field found using the model from M. L. Palumbo et al. [73]. When including the optimized non-radial velocity field, model I-CB’s binned RMS drops from 3.51 to 3.06 m s$^{-1}$ for Fe I 5434 Å. However, this improvement manifests across all 22 GRASS-E template lines inconsistently, with some lines having little to no improvement and other lines having a drop in RMS up to 2 m s$^{-1}$ when including the optimized non-radial velocity field. This incoherence across the spectral lines suggests that rather than capturing the true velocity field, the optimized non-radial velocity field may provide the model with additional degrees of freedom that could absorb line-by-line variability.

Supergranulation

Supergranulation refers to the convection process associated with supergranules, which consist of clusters of individual granules and have diameters ranging from 30-35 Mm [87]. Unlike granulation that has predominantly vertical plasma flows, supergranulation occurs mostly horizontally with 300-400 $m\ \mathrm{s}^{-1}$ horizontal flows and weaker 20-30 $m\ \mathrm{s}^{-1}$ vertical flows ([87]).

Suspecting supergranulation, we turn to using solar velocity fields from the Solar Dynamics Observatory’s (SDO) Helioseismic Magnetic Imager (HMI; [93]). Following the methods from R. D. Haywood et al. [43] and M. L. Palumbo et al. [73], we correct the HMI Dopplergrams for satellite motion and differential rotation in order to estimate a solar velocity field that accounts for all solar processes, such as convective blueshift, supergranulation and meridional flows. We input the SDO/HMI velocity field for 6173 \AA{} into our model I, which computes weighted average projected RVs excluding the optimized convective blueshift polynomial, and find that the resulting RM curve is very similar to the original results of model I for 6173 \AA{} without any convective blueshift. In other words, the HMI Dopplergrams do not provide any additional information past our initial weighted average projected RVs. A. Reiners et al. [83] also found that HMI Dopplergrams do not include information of convective blueshift

Summary

The discovery of Earth-mass exoplanets requires accurate mitigation of granulation. Although the impact of granulation on spectral lines is well understood, we currently lack the ability to accurately account for granulation within our RV models. Convective blueshift, specifically, is well-known to influence the shape of a RM curve. Therefore, we use RVs of the 2023 October 14 solar eclipse from the NEID solar feed and NEID spectrograph to test different RV models and probe our ability to model the influence of convective blueshift. We present three models: model I computes weighted average projected RVs, model II computes CCF RVs from simulated spectra with no granulation signatures, and model III computes CCF RVs from simulated GRASS-E granulation-driven spectra. We find that GRASS-E performs better than the models with no knowledge of granulation signatures. However, when optimizing and including a convective blueshift polynomial in the models that previously did not account for granulation signatures, we are able to outperform GRASS-E’s ability to recover the solar eclipse RM curve by more than a factor of 3 in RV RMS. In addition, we find that including a non-radial velocity field in our model allows us to further reduce our RV RMS for the majority of the considered spectral lines.

We cannot fully dismiss sources of uncertainty such as: potential inaccuracies in the adopted limb darkening and differential extinction corrections, incomplete p-mode mitigation, and the imperfect granulation-induced CLV. Regardless, our efforts show improvements from previous efforts to model convective blueshift during a solar eclipse. For example, the 2015 March 20 solar eclipse observed by the IAG Fourier Transform Spectrograph was modeled similarly to our model I and I-CB but using a sigmoid parameterization of convective blueshift (A. Reiners et al. [83]). Their version of our model I-CB (using 5000 – 6200 and 6500 – 6700 \AA{}) has structured residuals up to $\sim40\ \mathrm{m}\ \mathrm{s}^{-1}$. The residual structure they find is similar in shape to the residuals from our models I and II, which peak at $\sim75\ \mathrm{m}\ \mathrm{s}^{-1}$. When using interpolated synthetic line profiles instead to compute RVs, A. Reiners et al. [83] are able to reduce their peak residual to $\sim25\ \mathrm{m}\ \mathrm{s}^{-1}$, but the residual structure persists. On the other hand, our models using an optimized convective blueshift polynomial have unstructured residuals, suggesting that the residual structure is dependent on the correction of convective blueshift.

We find that when attempting to model the solar eclipse RM curve, GRASS-E is empirically limited by the sampling of existing center-to-limb observations as GRASS-E currently takes a piecewise constant function approach for its CLV models. To overcome this limitation, inter-/extra-polations of the GRASS-E v2.0 template line asymmetries and convective blueshift could be conducted in order to better resemble continuous CLV. Such modifications to GRASS-E are beyond the scope of this paper, which solely aims to evaluate GRASS-E’s current granulation model. In addition, more finely sampled spectral observations at well-defined $\mu$ values would aid our understanding and modeling of CLV. Such observations may be obtained with the 4 m Daniel K. Inouye Solar Telescope (DKIST, T. R. Rimmele et al. [86]), upcoming Paranal solar ESPRESSO Telescope (PoET, N. C. Santos et al. [92]), and upcoming Micromirrors for Resolved Adaptable Solar Observations Limb-to-limb (MiraSOL, C. Robles & S. Mahadevan [89]; C. Robles et al. [90]). The advent of the 30 m class ELTs and their associated high-resolution spectrometers (e.g., The GMT-Consortium Large Earth Finder (A. Szentgyorgyi et al. [101]) and ANDES (A. Marconi [61])) will allow us to perform similar tests on stars with transiting planets, probing granulation across cool, Sun-like stars. However, since granulation-induced CLV remains extremely challenging to measure spectroscopically, the most promising way to parameterize granulation-induced CLV may just be through highly precise RV observations of a solar eclipse with instruments such as NEID. Our work is an important step in this direction.

Acknowledgments

We thank the anonymous referees for their thorough reviews, which have improved the clarity of this work.

AI tools were used to assist in debugging the source code presented in this work; the tool used was Claude Opus 5.0.

Based on observations at Kitt Peak National Observatory, NSF NOIRLab, managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. Data presented were obtained by the NEID spectrograph built by Penn State University and operated at the WIYN Observatory by NOIRLab, under the NN-EXPLORE partnership of the National Aeronautics and Space Administration and the National Science Foundation. Development of NEID was performed for the Jet Propulsion Laboratory, California Institute of Technology, sponsored by the United States Government under the Prime Contract 80NM0018D0004 between Caltech and NASA. We acknowledge support from NASA Grant 80NSSC26K0554. WIYN is a joint facility of NSF NOIRLab, Indiana University, the University of Wisconsin–Madison, the Pennsylvania State University, and Princeton University. Kitt Peak National Observatory sits atop I’oligam Du’ag (Manzanita Shrub Mountain). Astronomers are honored to be permitted to conduct scientific research on the sacred mountain located in the homelands of the Schuk Toak District within the Tohono O’odham Nation. We honor their past, present, and future generations, who have lived here for time immemorial and will forever call this place home.

The Center for Exoplanets and Habitable Worlds is supported by Penn State and its Eberly College of Science. The authors recognize the Penn State Institute for Computational and Data Sciences (ICDS) (RRID:SCR_025154) for providing access to computational research infrastructure (RRID:SCR_026424). This content is solely the responsibility of the authors and does not necessarily represent the views of the Institute for Computational and Data Sciences.

EG was supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE1255832. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

MLP was supported by the Flatiron Research Fellowship at the Flatiron Institute, a division of the Simons Foundation.

This research has made use of NASA’s Astrophysics Data System Bibliographic Services and NASA’s Navigation and Ancillary Information Facility data services.

Facilities: NEID/WIYN (NEID solar feed)

Software: Julia ([8]), JuliaGPU ([7]), GRBSS ([76]), RecteSpecML ([31]), Matplotlib ([46]), barycorrpy ([47]), OptIm ([70]), pandas ([63]), SPICE ([1]), astropy ([5]), sunpy ([102]), numpy ([103])