$\texttt{GRASS}$: Distinguishing Planet-induced Doppler Signatures from Granulation with a Synthetic Spectra Generator

Owing to recent advances in radial-velocity instrumentation and observation techniques, the detection of Earth-mass planets around Sun-like stars may soon be primarily limited by intrinsic stellar variability. Several processes contribute to this variability, including starspots, pulsations, and granulation. Although many previous studies have focused on techniques to mitigate signals from pulsations and other types of magnetic activity, granulation noise has to date only been partially addressed by empirically-motivated observation strategies and magnetohydrodynamic simulations. To address this deficit, we present the GRanulation And Spectrum Simulator ($\texttt{GRASS}$), a new tool designed to create time-series synthetic spectra with granulation-driven variability from spatially- and temporally-resolved observations of solar absorption lines. In this work, we present $\texttt{GRASS}$, detail its methodology, and validate its model against disk-integrated solar observations. As a first-of-its-kind empirical model for spectral variability due to granulation in a star with perfectly known center-of-mass radial-velocity behavior, $\texttt{GRASS}$ is an important tool for testing new methods of disentangling granular line-shape changes from true Doppler shifts.


INTRODUCTION
Since the detection of the first exoplanet orbiting a main sequence star (Mayor & Queloz 1995), the radial velocity (RV) method has enabled the discovery of several hundred planets to date (see the NASA Exoplanet Archive; Akeson et al. 2017). Despite these successes, the RV method will need additional development to address stellar variability and characterize Earth-analog planets around solar-type stars. In principle, the RV method uses the planet-induced Doppler reflex motion of host stars to infer the existence of planetary companions. In practice, the host stars exhibit intrinsic astrophysical variability that complicate or preclude the measurement of minute reflex velocities due to smaller, Earth-mass planets. Although advances in instrumentation have dramatically improved spectrograph stability and spectral sensitivity, the challenges posed by intrinsic stellar variability remain a significant barrier to the discovery of lower-mass, rocky planets in or near the habitable zone of other stellar systems.
Several physical mechanisms, which operate on many different timescales, drive stellar variability, including pulsations, magnetic activity, and granulation (for a detailed review of these and other phenomena see Meunier 2021). Pulsations refer to the rich spectrum of acoustic arXiv:2110.11839v1 [astro-ph.SR] 22 Oct 2021 modes excited by instabilities in near-surface convection. For dwarf F-, G-, and K-type (FGK) stars, the principal pulsations occur on timescales of minutes and lead to residual Doppler velocity amplitudes large enough to be problematic for the search for Earth-like exoplanets (i.e., multiple tens of centimeters per second from pulsations alone, compared to the ∼10 cm s −1 signal expected for an Earth-twin; e.g., O'Toole et al. 2008 andMichel et al. 2008). Recently, Chaplin et al. (2019) have shown that carefully-planned exposure times equal to or greater than the period of the principal p-mode can effectively mitigate the observed pulsation amplitude down to the 0.1 m s −1 level necessary for Earth-twin detection.
In comparison to pulsations alone, stellar magnetic activity has posed a relatively complex barrier. Starspots, cooler photospheric regions where penetrating magnetic fields inhibit surface convection, produce a rotationallymodulated, wavelength-dependent flux deficit that alters the shape of spectral lines and can therefore mimic the Doppler wobble induced by a planetary companion. In the chromosphere, concentrated magnetic field lines can produce bright regions known as plages, which tend to spatially coincide with smaller-area photospheric faculae, introducing further perturbations. The behavior of these chromospheric regions is complicated: not only do they exhibit more complex center-to-limb brightness variations, but they can be upwards of an order-of-magnitude more extended in area than spots (e.g., Chapman et al. 2001;Meunier et al. 2010). Indeed, posited planets have been shown to actually be stellar activity masquerading as Keplerian signals (e.g., Robertson et al. 2014 andpotentially Jeffers et al. 2014). In one especially insidious case, stellar activity coupled with the observing window function led to the false detection of an Earth-mass planet in the α Cen B system (Dumusque et al. 2012;Rajpaul et al. 2016).
Approaches to mitigating spurious velocities introduced by stellar activity have been multifaceted. In the observational realm, suspected planets with orbital periods near integer multiples of the stellar rotation period are generally disregarded since spot groups can persist for multiple stellar rotations and inject power into harmonics of the rotation period (Robertson et al. 2014). Recent astrophysical insights into radial velocity "jitter" have also produced guidance for identifying and prioritizing stars with minimal intrinsic radial velocity variability (Luhn et al. 2020). To combat activityinduced noise in observed stars, astronomers have used different spectral lines, such as Hα (e.g., Kürster et al. 2003;Hatzes et al. 2015) and the Ca H&K doublet (e.g., Wilson 1968;Wright et al. 2004;Isaacson & Fischer 2010), to simultaneously diagnose activity in stellar photospheres and chromospheres. However, these activity indicators are not universally exploitable: diagnostic lines will show variability among stars of differing spectral types, and even different tracers may differ on the same star. Often, different indicators must be calibrated to suit specific instrumental spectral ranges, especially in the case of M dwarfs (Robertson et al. 2016).
Recently, astronomers have begun applying advanced statistical models to RV data series in an attempt to disentangle activity from true Doppler shifts. Gaussian Process (GP) models for stellar variability have shown particular promise in separating out these two effects (e.g., Haywood et al. 2014;Rajpaul et al. 2015;Jones et al. 2017). Indeed, a large blind test conducted in Dumusque et al. (2017) suggested that GPbased analyses provided a promising way to model correlated signals due to stellar variability and could assist with extracting short-period planetary signals with low RV semi-amplitudes in the presence of stellar variability (Dumusque 2016). To aid in the further development of these novel statistical methods, the community has developed algorithms to simulate the photometric and spectroscopic impacts of magnetically-active regions present on the stellar photosphere (e.g., Herrero et al. 2016;Boisse et al. 2012;Dumusque et al. 2014). The respective tools produced in these publications, StarSim and SOAP 2.0, have seen widespread use, since they are able to efficiently produce simulated data sets such as those used in the development of GP models (e.g., Rajpaul et al. 2015;Jones et al. 2017;Gilbertson et al. 2020).
Astronomers have also used our nearest star, the Sun, as a tool to understand the RV impact of individual magnetic features. For example, Haywood et al. (2016) use simultaneous observations of sunlight reflected off the asteroid 4/Vesta and full-disk solar observations from NASA's Solar Dynamics Observatory (SDO) to show that observational proxies for stellar plage coverage may enable better inference of activity-induced RV variations. Follow-up direct observations of the Sun with a dedicated solar telescope feed into the HARPS-N spectrograph (Dumusque et al. 2015;Phillips et al. 2016;Collier Cameron et al. 2019) have reinforced this finding (Milbourne et al. 2019).
Although astronomers have made strides toward understanding activity from observational, theoretical, and statistical standpoints, novel approaches will be required to overcome another principal obstacle to centimeterper-second precision RVs: granulation. Granules, individual magnetoconvective cells which blanket the surfaces of cool stars, are comprised of hot and bright upwelling plasma, surrounded by rings of cooler, denser downwelling plasma known as intergranular lanes. The collective motion of granules and larger convective cells known as "supergranules" can introduce RV perturbations on the order of meters per second on timescales ranging from minutes to days (Rieutord & Rincon 2010;Rincon & Rieutord 2018). Unlike spots and plages resulting from activity, these convective phenomenon are not rotationally modulated, and so do not exhibit the same quasi-periodic time dependence. Therefore, distinct statistical methods and observational techniques may be required to distinguish spectral variability due to granulation from true Doppler shifts.
Although an effective treatment for granular noise has not been devised, the physical processes by which granules perturb spectral lines is well-understood. In disk-integrated flux, typical absorption lines for solartype stars exhibit a characteristic "C"-shaped asymmetry that is observable in the line bisector curve (Stathopoulou & Alissandrakis 1993). This characteristic shape results from the distinct motions within granules: more light is emitted by blueshifted upwelling plasma in granules compared to the redshifted downwelling plasma in lanes. An asymmetric absorption line results when the light from these distinct components is combined (Gray 2008).
In disk-resolved observations, solar line bisectors vary in shape with limb position (see Figure 1 and §2). At disk-center, typical bisectors exhibit the characteristic "C"-shape seen in disk-integrated light. Toward the limb, however, bisectors become distorted, more resembling a "\"-shape. This transformation is an effect of viewing angle: at disk-center, the line of sight passes down vertically through the convective cell, whereas at the limb, the granules are viewed in profile, sampling different line-of-sight velocities. Due to these effects and their short lifetime (of order minutes), granules induce spectral perturbations that vary with time in both diskresolved and disk-integrated light.
Thus far, attempts to quantify the precise RV impact of granulation-induced line perturbations have been conducted via complex magnetohydrodynamic (MHD) simulations of small patches of the stellar surface, as in Cegla et al. (2018Cegla et al. ( , 2019. Using their MHD (from Vögler et al. 2005) and radiative transfer (from Socas-Navarro et al. 2015) simulations, these authors synthesized realistic disk-resolved Fe I 6302Å absorption line profiles. In Cegla et al. (2019), the authors then tiled a sphere with simulated patches to model a star and reconstruct a disk-integrated line profile. Using this approach, they were able to show that diagnostics of line bisector shape strongly correlated with the size of the measured spurious RV shift. However, these results are based on the synthesis of a single line. Real stellar spectra contain a myriad of lines with varying formation heights, magnetic sensitivity, etc. Given the computational time required to first conduct the magnetoconvective simulations and then run the radiative transfer, this approach is not currently feasible for generating full synthetic spectra with dense time sampling. Moreover, the amplitude of granular variability naturally depends on the assumptions and physical simplifications made in solving the MHD and radiative transfer equations. Likewise, the presence, strength, and configuration of magnetic fields can greatly impact the inferred RV RMS. For example, Sulis et al. (2020) find an RV RMS of ∼0.507 m s −1 , similar to the value observed by Pallé et al. (1999). Yet with the inclusion of strong magnetic fields, Cegla et al. (2019) find a much smaller amplitude of only ∼10 cm s −1 .
Although these sorts of MHD simulations have yielded great physical insight, developing practical solutions to address convectively driven RV noise will require new, complementary approaches. As we have seen from other astronomers' solar observations and data-driven photospheric models (e.g., Dumusque et al. 2014;Haywood et al. 2016), we can glean a great deal of practical knowledge from codes designed to produce observationallyinformed synthetic data sets. Indeed, Sulis et al. (2017) indicate that simulated data sets could be used to train detection processes that are normally strongly affected by granular noise. Although simulation tools such as StarSim and SOAP 2.0 have had great success modeling line perturbations from stellar active regions, neither includes a comprehensive treatment of time-variable perturbations from granulation. Other past works have used the observed and simulated properties of granules (Meunier et al. 2015) and supergranules (Meunier & Lagrange 2019, 2020a to assess the impact of these convective motions on exoplanet detection prospects, though they only provide velocities and not spectra. To fill in these gaps, we have developed an observationallyinformed simulation tool to create synthetic spectra with time-resolved, granulation-induced variability. We dub this tool GRASS, the GRanulation And Spectrum Simulator. We envision that GRASS may be useful for training and validating various strategies for reducing the effects of stellar variability on precision Doppler measurements. In particular, mitigation strategies that make use of lineshape information could be tested on data sets generated by GRASS.
In the following sections, we describe GRASS in detail. In §2, we summarize the input observations used to create our model spectra, and briefly comment on the output of the simulation tool. In §3, we discuss the processes by which GRASS creates synthetic spectra from the input solar observations. In §4, we validate the GRASS modeling approach by comparison to disk-integrated reference observations. In §5, we present some first results obtained from GRASS single-line simulations. In §6, we comment on some potential uses, discuss future avenues, and also highlight the limitations and caveats of GRASS.

DATA
Rather than using MHD and radiative transfer simulations to produce synthetic spectra (as in previous studies; e.g., Cegla et al. 2018Cegla et al. , 2019, we instead use disk-resolved, high-resolution observations of the Sun to create our synthetic disk-integrated spectra. This modeling approach offers a number advantages. Principally, we know that our model captures the most important physics precisely because it is based on real data. Although the current version of GRASS explicitly excludes phenomena such as spots, this is in principle a solvable problem with additional observations or modeling. Our input observations are described below, in §2.1. In §2.2, we briefly describe the output of GRASS, reserving a more comprehensive description of our methodology for §3.

Input Data Description
Since line-shape variations driven by granulation vary with time and location on the disk, synthesizing our model disk-integrated spectra requires temporally and spatially resolved input data. To account for these variations in line shape, we use a set of solar observations originally published in Löhner-Böttcher et al. 2019 (hereafter L19) as a template for line synthesis. We refer to these data as the "input data." These input data consist of reduced spectra covering the Fe I 5434.5Å line, from which we measure bisectors and widths as a function of depth. Ideally, these data would be derived from high-resolution observations conducted simultaneously across the solar disk. To our knowledge, however, no such data set exists and obtaining such data would require a special-purpose instrument. Relevant spectroscopic parameters and quantities for the Fe I 5434.5Å line are summarized in Table 1.

Summary of Input Data Observations
The observations of this Fe line were carried out as described in Löhner-Böttcher et al. 2018 (hereafter L18). The observations, consisting of 20-minute sequences, were performed using the Laser Absolute Reference Spectrograph (LARS, Doerr 2015;Löhner-Böttcher et al. 2017) at the German Vacuum Tower Telescope Table 1. Parameters for spectroscopic lines in this work. Quantities denoted with a dagger ( † ) are from L19. All other quantities are from the NIST Atomic Spectra Database (Kramida et al. 2020 (VTT). LARS employs VTT's high-resolution echelle spectrograph in combination with a laser frequency comb (LFC) to achieve superb spectral resolution R ∼ λ/∆λ ∼ 7×10 5 with an absolute wavelength calibration accuracy of ∼0.02 mÅ (about 1 m s −1 ) in the wavelength region of interest. To sample line bisectors along different lines of sight, observations were carried out at eleven positions along each of four directions from disk center, oriented along and perpendicular to the sky-projected solar rotation axis (see Figure 1). The positions are radially parameterized as µ = cos(θ), where the heliocentric angle θ is the angle subtended by the line of sight and the solar surface vertical. The parameter µ is therefore 1 at diskcenter and 0 at the limb edge. Individual observations had a ∼10 field of view set by the fiber-coupling unit and were performed with a pointing accuracy of ∼1 .
To average over supergranular motions which would hinder the precise convective blueshift measurements performed in L18 and L19, observations at each limb position were quickly scanned over ellipses with area corresponding to the typical angular area of supergranules. The sizes and orientations of the ellipses were varied with position on the disk in order to account for the effects of foreshortening. At disk center, a circular scan pattern was used, with a total diameter of ∼30 accounting for the size of the fiber-coupling unit. Slightly away from disk center (from around µ = 0.95 to 0.7), where the amplitudes of p-modes dominate over motions from supergranulation, the scanning ellipses had total axes dimensions of about 40 ×10 (i.e., including the width of the fiber-coupling unit). Toward the limb (µ ≤ 0.5), where horizontal supergranular motions dominate over p-modes, the total ellipse axes dimensions were about 40 ×30 .
Observations were carried out in sequences of ∼20 minutes in length, with individual exposures conducted at a cadence of 1.5 seconds and then binned into tenexposure (15-second) bins. To avoid the systematic effect of changes in convective blueshift introduced by active regions, simultaneous G-band images from the LARS context camera and recent magnetograms from the Solar Dynamics Observatory (SDO) Helioseismic Magnetic Imager (HMI, Pesnell et al. 2012;Scherrer et al. 2012) were used to ensure that only quiet Sun regions were observed. In the case of the HMI magnetograms, the automatically identified active region patches, which are flagged by the data reduction pipeline as pixels with an absolute line-of-sight magnetic field strength greater than 100 G, were used to reject active regions (Hoeksema et al. 2014). Subsequent data reduction, including absolute wavelength calibration and corrections for barycentric motion and solar rotation, were performed as described in Löhner-Böttcher et al. (2017) and L18. For lines observed by L18 and L19, the observed wavelengths of solar lines were empirically measured with respect to the lab-frame wavelengths of Fe I lines to an accuracy of below 0.1 mÅ, or about 2 − 4 m s −1 , using the LARS LFC in combination with a hollow cathode lamp.

Output Data
As we describe in detail in §3, we spatially interpolate and integrate the input data to create synthetic time series spectra like those for an unresolved point source. We refer to these spectra, and the associated cross-correlation functions and velocities, as the "output data" of the simulation.

Program Language, Installation, and Performance
The latest version of GRASS and its associated documentation are available from GitHub 1 . The frozen version of GRASS accompanying this work (v1.0) is archived on Zenodo (Palumbo et al. 2021). GRASS is written entirely in Julia 2 , a high-level, performance-oriented language. To generate synthetic spectra, the user must specify some input parameters in a brief Julia script (see the associated documentation and tutorials on GitHub). In addition to the output spectra, users can choose to store derived data products such as crosscorrelation function (CCF) profiles and/or apparent ra-dial velocities measured from the spectra using one of several algorithms.
GRASS was designed to produce synthetic diskintegrated spectra more quickly than can be done through MHD and radiative transfer simulations. The compute-time and memory requirements of GRASS are primarily dictated by three parameters: the grid size of the simulated disk (N 2 , see §4.2), the number of timesteps simulated (N t , with one step corresponding to one 15-second input observation), and the number of spectral resolution elements synthesized (N λ ). We find that the computation time scales quadratically with N , and linearly with both N t and N λ . For this work, we simulate single-line spectra to validate GRASS against reference disk-integrated observations. For these runs of GRASS, we use N 2 = 132 2 , N t = 50, and N λ = 193 (corresponding to a spectrum 1.5Å wide with spectral resolution R = 7 × 10 5 at λ = 5434.5Å). Depending on the exact hardware configuration, we find that we can generate 50 instances of output spectra in less than thirty seconds on desktop hardware.

Creating Synthetic Spectra
GRASS is designed to create realistic time series of synthetic spectra with granulation-driven line-shape variation. The creation of these synthetic spectra from the input data is carried out in three broad steps, which we describe in the subsequent sections: • initial"pre-processing" of the input data, in which the reduced spectra are used to measure line bisectors and widths as functions of depth within the line ( §3.2.1), • simulation of a spatially-resolved stellar disk, in which the pre-processed data are used to re-create synthetic disk-resolved line profiles on a model stellar disk ( §3.2.3), • and spatial integration of the model stellar disk, in which the simulated disk-resolved spectra are coadded into a disk-integrated spectrum like those observed for point sources ( §3.2.4).
The output data consisting of disk-integrated, timeseries synthetic spectra can optionally be binned and used to compute cross-correlation functions and apparent radial velocities, as is typical for current RV surveys. We detail our procedure for these computations in §3.2.5 and §3.3.

Input Data Pre-processing
As described in §2.1, the input solar data are reduced spectra binned to a 15-second cadence. To uniquely specify line shape from the core to the wings, we directly measure line bisectors and widths as functions of depth from these reduced spectra. To do so, we isolate the section of the reduced spectrum containing only the observed line of interest (i.e., the Fe I 5434.5Å line), which we then interpolate onto an evenly-sampled intensity grid. The width at each depth is then calculated as the difference in wavelength between the red and blue sides of the line, and the bisector as the average.
In implementation, we perform a number of additional steps to produce smooth and accurate bisector and width measurements. We first average the flux across adjacent pixels in areas where the line profile is not smooth. To circumvent the impact of shallow line blends in the line wings, we model the far line wings above 93% continuum level as Voigt functions. We use separate non-linear least squares fits of the red-ward and blue-ward wings (excluding pixels containing significant blends) to determine the best model parameters. Since the uncertainty in the bisector measurement increases where the slope in the line is small (Gray 2008), we find that the top-most bisector measurements are quite inaccurate. To avoid the impact of these spurious measurements, we model the bisectors above 80% continuum flux as the linear extrapolation of the bisector data measured between 70% and 80% continuum flux.
For many limb positions, the full sequence of input data were not observed continuously. Often, several minutes of spectra were observed in one session, with the other session observed months later. In these cases, we adjust the mean velocity of the latter bisector series to match that of the first. For all results presented in this paper, we have performed additional simulations using only the longest set of contiguous data for each limb position. We find no significant difference in our results between these two sets of simulations. However, caution should nevertheless be taken when producing or interpreting series of spectra longer than ∼20 minutes.
In some cases, we exclude highly spurious instances of the input data. We discard any epochs for which the bisector is more than 3σ removed from the mean bisector measurement at a given intensity and limb position. Due to apparent errors in the raw spectra reduction, we also entirely discard the data collected on 2017/10/17 for µ = 0.7 on the northern solar axis. We instead only use the data for this position collected in 2016. Although this reduces the total time-baseline of available data, this section of the Sun only contributes ∼4% to the total disk-integrated flux in our simulations.
Example pre-processed bisectors and width measurements are shown in Figure 2. Variability from pulsations and granulation are clearly visible in the bisectors. Whereas pulsations largely only change the horizontal position (i.e., velocity) of the bisector curves, granulation can change the curvature, shape, and position of individual bisectors. The line width measurements do exhibit a small degree of variability, though the level of variability is very small compared to that of the bisectors.
The initial pre-processing described above is performed only once, and the resulting bisectors and width measurements are written to output files that are included in the GRASS GitHub repository. Upon execution of GRASS, the files containing these pre-processed data are read into memory and sorted into arrays that are assigned µ and axis (i.e., North, East, South, West with respect to the solar rotation axis) identifiers. For each unique (µ, axis) pair, the corresponding array contains the full time-series of pre-processed bisectors and width measurements for that disk position. As detailed in the subsequent sections, these data are then used to generate synthetic absorption line profiles.

Adjusting Input Line Depths
Following from the detailed explanation of stellar line formation presented in Gray (2008), we can use the information encoded in the deep Fe I 5434.5Å line to also model shallower lines. As explained in Chapter 17 of Gray (2008) and shown in their Figure 17.14, the bisectors of shallower lines resemble the upper portions of the bisectors of deeper lines. This fact is driven by the nature of deep line formation in stellar atmo-spheres. Because the intensity in a given portion of a line is controlled by the source function at a corresponding physical depth in the stellar atmosphere, the wings of deep lines, which form deep in the photosphere, may be thought of as shallow lines themselves. Lines from species with identical lower energy-level abundance distributions should therefore share line profiles. Given this fact, we have written GRASS to allow the user to generate synthetic lines of arbitrary continuum-normalized depths, assuming that the information encoded in the deep Fe I 5434.5Å line samples velocity structure along deep slant depths in the solar atmosphere.
To model shallower lines, GRASS must transform the pre-processed input data. In the case of the width measurements, GRASS simply scales the data onto the new range of intensities for the shallower line. For the bisector measurements, GRASS discards measurements for intensities below the depth of the line and interpolates the remaining points onto a same-length grid. Since the input line bisectors are extrapolated above 80% of the continuum level, synthetic lines of depths less than 20% are based solely on extrapolated data. The consequences of this approach with respect to Doppler velocity measurements are analyze and discussed in §5.1.

Simulating a Stellar Disk
Following the initial pre-processing of the input data, GRASS models a stellar disk by creating a uniform Cartesian xy spatial grid of N × N pixels. The model stellar disk is described by a unit circle inscribed within the square grid. At each point on this model disk, GRASS creates a synthetic, spatially-resolved spectrum consisting of any number of synthetic lines specified by the user. To create the synthetic line profiles, GRASS first finds the nearest solar axis and value of µ to the point of interest on the grid. Then, it identifies the corresponding input solar observations and chooses a starting epoch from the input data. This epoch is chosen randomly for each spatial grid cell, which ensures that synthetic lines in separate spatial cells do not move in concert.
On a given spatial grid cell, GRASS synthesizes each absorption line by using the input data to assign the wavelength at each depth for the left-and right-hand sides of the synthetic line. The resulting absorption line profile, which is initially specified on a uniform grid of intensities, is then interpolated onto a uniform wavelength grid of the spectral resolution specified by the user. This process is repeated for any additional synthetic lines. When lines are sufficiently close (i.e., in the case of blends), the separate absorption profiles are multiplied to create a blended line profile.
After all lines have been synthesized within a given spatial cell, GRASS applies a rotational Doppler shift based on the position of the cell on the disk. We use the differential rotation model from Snodgrass & Ulrich (1990), where the angular rotational velocity ω varies as a function of a solar latitude φ: We use the best-fit values reported in Snodgrass & Ulrich (1990) for the coefficients A, B, and C. The true rotational velocity v rot at some latitude φ is then calculated as: where c is the speed of light, and the stellar radius R * is set to unity. This velocity is then projected onto the line of sight for the given sky-projected inclination of the rotation axis. The default inclination is given by the 3D unit vector (0, 1, 0) corresponding to an equator-on star.
We note that this calculation, as well as the spatial xy grid, assume a spherical star (observed as a circle in projection). In the case of rapidly rotating stars, this assumption may not be valid.
In addition to the effects of differential rotation, GRASS must also account for limb darkening. To obtain I µ , the limb-darkened intensity at a given µ, we use a simple wavelength-independent, quadratic limb darkening law from Kopal (1950): where the default values for the coefficients are u 1 = 0.4 and u 2 = 0.26, typical solar values for this wavelength region (as verified via Southworth 2015). We take I 0 , the intensity at disk center, to be 1. For simulations of other wavelength regions or stars with different limb darkening behavior, the limb darkening coefficients may be easily changed by the user.

Integrating Over the Stellar Disk
As described in the previous sections, GRASS essentially constructs the intensity as a function of time t, wavelength λ, and position on the disk (x, y): I(t, λ, x, y). To create a spatially-integrated spectrum for a point source I(t, λ), we sum intensities over the spatial dimensions: where C is a normalization constant such that the diskintegrated continuum level is normalized to unity when the individual disk-resolved spectra are summed over spatial dimensions. Positions outside the stellar disk contribute zero intensity. We compute C = π/(2N 2 ), where N is the length of the N × N spatial grid. In implementation, we allocate only a N λ × N t array (where N λ is the number of spectral resolution elements and N t is the number of simulated observation epochs) and perform the spatial summation without storing the intensity at each location. This process circumvents the need to allocate a much larger N × N × N λ × N t array, significantly reducing memory overhead.
We note that our method of disk simulation and integration destroys any coherency of p-modes that is naturally present in the input data. Since each spatial grid cell is initialized with a random starting epoch selected from the input data, the contribution of the p-modes to the RV RMS is effectively averaged over. Thus the RV RMS measurements presented in §5 may appear lower than those that pulsations would be expected to produce. Because we are attempting to probe only the granulation-induced noise, this is actually desirable. However, we do note that an analysis of the power spectrum of synthetic RV time series produced by GRASS (not shown) suggests that the peak residual pulsation amplitude persists at a level of 5 cm s −1 .

Simulating Realistic Observations
While GRASS can synthesize spectra with arbitrary spectral resolution, the process described in the previous sections creates synthetic spectra at the same spectral resolution as the input data (R ∼ 7 × 10 5 ). Though these sorts of spectra may be useful for some purposes, they are not typical for current Doppler planet survey observations, which are conducted on instruments with R ∼ 50, 000 − 150, 000.
In order to degrade the spectral resolution to more typical values, we perform a convolution with a truncated Gaussian line spread function (LSF) with characteristic width σ LSF computed from the desired lower spectral resolution. By default, the output spectra are degraded to R ∼ 1.17 × 10 5 , comparable to the average spectral resolutions of the HARPS-N and NEID spectrographs. As in real RV spectrographs, we oversample these spectra with about ∼6 pixels per spectral resolution element in the region of 5434.5Å. To simulate photon shot noise, we then add Gaussian white noise to each pixel in order to achieve a specified signal-to-noise ratio (SNR). In doing this, we assume that the Poisson noise is in the Gaussian regime.
By default, the synthetic spectra are generated at a 15-second cadence (the same as the input solar data described in §2.1). However, real observations of stellar spectra (i.e., for stars other than the Sun) use larger integration times and have practical constraints on the observing window function. For the analysis presented in §5.3, we compare synthetic observation patterns by defining parameters such as number of observations, exposure times, and any time between exposures (including readout and off-target time). We achieve this by appropriately weighting and binning the consecutive 15second "snapshots" produced by GRASS. To save computation time and resources, spectra are not synthesized for simulated read-out or off-target times.

Measuring Velocities from Synthetic Spectra
The spectral imprints of true planet-induced stellar velocities and granular velocities are distinct. Whereas Doppler reflex motions shift the stellar spectrum, granular motions perturb individual line shapes, creating skewness and kurtosis in their profiles. Although these perturbations should not necessarily be characterized or interpreted as wholesale velocities, these small effects are captured as an apparent source of stellar variability noise in cross correlation function (CCF) analyses. In order to explore the impact of granulation on these analyses, we have implemented a CCF-based algorithm from Ford et al. (2021) for measuring apparent RVs in the GRASS package.
The implementation of our CCF algorithm follows that of Pepe et al. (2002) and additionally accounts for the appropriate treatment of Doppler shifts discussed in Wright (2019). We calculate the CCF by shifting a template "mask" in velocity along each spectrum: where I(λ) is the observed or simulated spectrum and M (λ v R ) the mask at some velocity-shifted wavelength λ v R given by the classical Doppler formula where λ 0 is the rest-frame wavelength of the line. The mask M (λ v R ) is the weighted sum of individual top-hat functions M i (λ v R ) centered on rest-frame line centers: We take the weights w i to be the continuum-normalized depth of line i. For actual observations, the weights may also take into to account the instrument efficiency at the location of each line. Keeping with common practice, we shift the mask in velocity increments smaller than the spectral resolution of the synthetic spectrum to compute an over-sampled CCF. As in Pepe et al. (2002), we compute a velocity from the CCF by fitting with a Gaussian model. We first estimate a full-width at half-maximum of the CCF profile, and then perform the Gaussian fit via non-linear least squares on the data around the CCF minimum. We find that fitting within ±0.69σ of the minimum reliably retrieves injected velocities for lines with widths typical of the Sun with errors well below 1 cm s −1 . We report the mean of this Gaussian fit as the apparent Doppler velocity of the spectrum, which is expressed relative to the velocity zero point defined by the rest wavelengths of the lines in the CCF mask.

Validating GRASS against real solar observations
In order to validate and assess our line synthesis model, we compare our synthetic disk-integrated absorption line profiles and bisectors to those extracted from the IAG solar atlas (Reiners et al. 2016). To compare line profiles, we first use GRASS to synthesize a disk-integrated absorption line profile for a model Fe I 5Å absorption lines in disk integrated light are overplotted. We additionally plot the IAG spectrum cleaned of shallow line blends in green, which nearly completely covers the black synthetic curve. Residuals between the IAG spectrum and the synthetic spectrum are shown in the bottom panel. Our synthetic spectrum deviates from the observed IAG spectrum at most by ∼5% (blue squares), but only by ∼0.5% from the cleaned IAG spectrum (green circles). Right: The observed IAG (blue) and synthetic (black) Fe I 5434.5Å bisectors, which have been arbitrarily shifted horizontally (i.e., in velocity) to achieve maximal alignment, generally agree in shape up to 60% of the continuum level. The cleaned IAG (green) and synthetic bisectors are in much greater agreement. The residuals shown in the right panel further show that the blend cleaning greatly improves the agreement between model and observation, especially in the upper portion of the bisector.
5434.5Å line. We set the continuum-normalized depth to match that of the line in the IAG atlas. We then shift the synthetic and observed line profiles to the same rest-frame velocity. The two profiles are shown at left in Figure 3 as the black and blue curves. The synthetic and IAG spectra are in general agreement, with deviations in flux between the two only as large as ∼5%. We attribute these larger deviations to blends with shallow lines that were modeled out of the input data in the pre-processing stage (see §3.2.1). These blends are more clearly visible in the residuals shown in the bottom panel of Figure 3. Any additional deviations are likely introduced by the presence of active regions in the IAG spectrum, which was observed from a different site during a period of heightened solar activity in 2014 (Hathaway 2015;Reiners et al. 2016).
To more cleanly compare and evaluate our line synthesis model, we model line blends in the IAG spectrum by iteratively fitting Gaussian profiles to the largest deviations between the IAG and synthetic spectra. Although real stellar spectra will contain many blends between both stellar and telluric lines, we nonetheless choose to model a "clean" line profile in this work for three primary reasons. Firstly, our input data were observed at lower SNR than the IAG spectrum, meaning that line blends in these spectra are not strongly detected. We find that smoothing over the line wings in the preprocessing stage ( §3.2.1) allows greater numerical accuracy in bisector and width calculations. Secondly, telluric line blends will naturally vary with observing time and location, so we should not expect those tellurics seen in our input data to match those observed by another instrument. Thirdly, we prefer the approach of using a single, "clean" template line to create synthetic spectra. Given this template line, GRASS can create simulated blended profiles of arbitrary overlap.
The cleaned IAG spectrum is shown at left in Figure 3 as the green dotted curve, with residuals plotted as green points in the bottom panel. We find that cleaning the IAG spectrum by dividing out the best Gaussian fit for each line reduces the largest flux deviation between the spectra to the ∼0.5% level. Indeed, in the top panel of Figure 3, the curves for the cleaned IAG and synthetic spectra essentially overlap.
To further validate our model, we compare CCF bisectors for the synthetic and IAG line profiles. In each case, the CCF is computed using a single top-hat mask entry at the rest wavelength of 5434.5232Å. We then measure a bisector from each CCF using the same method described in §3.2.1. The bisectors are arbitrarily shifted in velocity to achieve maximum alignment to aid in visual comparison. These bisectors are shown at right in Figure 3. The synthetic (black) and IAG (blue) bisectors show remarkable agreement from ∼20% to ∼60% of the continuum-normalized flux. Above ∼60%, the presence of the aforementioned line blends move the ob-  Figure 4. Effect of simulation spatial grid size on the generated synthetic spectra. Points indicate the mean RMS calculated from many RV time series, and error bars the standard error of the mean. Two measurements of the RV variability from solar granulation are shown as the orange forwardhatched (Elsworth et al. 1994) and green back-hatched (Pallé et al. 1999) regions. We find that our ideal grid size of N 2 = 132 2 (vertical dashed line), chosen to produce grid cells with angular area matching the weighted-average spatial resolution of the input observations, intersects the power law relation within the Pallé et al. (1999) region. At spatial resolutions greater than N ∼ 400, the angular size of the grid cells becomes much smaller than that of the smallestfootprint input observations and consequently loses physical validity.
served solar bisector artificially to the right. Comparing the synthetic (black) and cleaned IAG (green) curves, the improvement in agreement is notable. In all cases, we only compute bisectors up to 90% flux, since the measurement uncertainty in the bisector computation inflates drastically above this level.

Effects of Spatial Resolution
To create synthetic spectra, GRASS creates diskresolved line profiles on a uniform spatial grid (as described in §3.2.3). The ideal resolution of this spatial grid over the stellar surface would match that of the input observations. However, because the spatial resolution of the input observations was designed to approximate the typical angular area of supergranules and consequently varies with position on the disk (as described in §2.1), simulating an ideal grid is impractical. Therefore, we choose a rectangular grid and set N (the number of grid cells along each Cartesian xy direction in the sky plane) such that the area of the spatial grid cell matches the weighted-average area of the input data observations. Given the sizes of the elliptical oscillations described in §2.1 and L18, we compute N 2 ∼ 132 2 as our optimized number of grid cells.
Given the complexity of the spatial scanning patterns in the input data, we further validate our choice of spatial resolution by comparing the RV RMS of synthetic time series to direct measurements of the granulation RMS in the literature. To perform this comparison, we repeatedly generate time series of synthetic spectra consisting of a single absorption line at several spatial grid sizes. These noiseless spectra were generated at the full 15-second cadence of the input solar data and with the native spectral resolution of R = 7×10 5 . The single synthetic line was centered on 5434.5Å and was assigned a continuum-normalized depth of 0.8, comparable to the depth of the template line in the input data. We then measured velocities from these spectra using the CCF method described in §3.3. This process was performed many times to measure a robust mean and sample variance of the mean-adjusted RMS of the RV time series.
The results of this test are shown in Figure 4. We find that the RMS of the measured velocities falls off as a power law with a best-fit index of α ∼ 1 (dashed gray curve). Ideally, our chosen grid size should produce a synthetic RV RMS similar to that observed in the real Sun. Unfortunately, as summarized by Cegla et al. (2019), this value is not precisely known and varies depending on both the line measured and time of observation. In Figure 4, we show the 2σ confidence bands for two measurements of the granulation RMS performed by Elsworth et al. 1994 (31.9 ± 9 cm s −1 , orange forwardhatched) and Pallé et al. 1999 (46.1 ± 10 cm s −1 , green back-hatched). These horizontal curves intersect our relation in the approximate region of N = 2 7 to N = 2 8 . Moreover, we find that our chosen grid size of N 2 = 132 2 (vertical dashed line) intersects the power law relation within the 2σ confidence interval for the solar granulation RMS measured by Pallé et al. (1999), suggesting that our chosen N produces physically valid levels of variability. Given this validation, we choose N = 132 (i.e., N 2 = 132 2 ) as the default length of our spatial grid for the other simulations and analyses presented in this work. As a utility, we allow users to change the spatial resolution of simulations in GRASS. However, the default value is set to N 2 = 132 2 , and GRASS reports a warning if another value is chosen.

APPLICATIONS AND FIRST RESULTS
GRASS is a powerful tool for creating synthetic stellar spectra. By design, these spectra contain minimal magnetic activity-induced velocities and are free from planet-induced Doppler reflex motions. Ostensibly, the only sources of variability encoded in each epoch origi- . Doppler velocity variability weakly scales with synthetic line depth for lines deeper than 20% of continuum level. As in Figure 4, points indicate the mean RMS calculated from many RV time series, and error bars the standard error of the mean. The shaded blue region encloses 68% of all RMS values across all runs. From line depths of 0.5% to 20%, the RMS of an RV time series obtained from GRASS output spectra appears to decrease. From 20% to 95% depth, the RMS then appears to steadily increase. However, the shallowest lines (gray hatched region) are synthesized exclusively from extrapolated input data (see Figure 2 and associated text). The trend in this region may therefore be unreliable.
nate with granulation and pulsation. Therefore, by creating these spectra, we can assess in detail those perturbations created by granulation. In the following subsections, we analyze single-line spectra generated by GRASS and offer preliminary observation advice for the mitigation of granulation noise.

Synthetic Line Depth and RV Variability
Given the asymmetric nature of stellar absorption lines and the trends in bisector shape seen with line depth (Gray 2008), it is conceivable that lines of different depths will show varying levels of apparent Doppler variability. To test this possibility, we again generate many synthetic spectra consisting of single lines at different depths. We again add no noise to these spectra, and preserve the native spectral and temporal resolution of the input data. We plot the mean RMS of the radial velocities measured from these spectra against line depth in Figure 5. We find that the RV RMS scales with depth above a threshold of 20% continuum flux. Lines shallower than 20% continuum flux also show increased variability. However, we caution that these lines of lower depths are modeled solely from extrapolated input data (see Figure 2 and associated text), potentially compromising the reliability of these RV RMS measurements. Pole-on Equator-on Figure 6. RV time-series variability and stellar inclination. As in Figure 4, points indicate the mean RMS calculated from many RV time series, and error bars the standard error of the mean. As in Figure 5, the shaded blue region encloses 68% of all RMS values across all runs. Equator-on stars do not appear to show any significant difference in granulationdriven RV variability compared to pole-on stars. Additionally, neither the width of the RV RMS distributions (blue shaded region) nor the error on the mean (error bars) appear to scale with inclination, suggesting that the impact of additional uncertainty introduced by a broader CCF profile for the Sun when viewed equator-on is small in comparison to the variability introduced by granulation signals.
Considering only lines deeper than 20%, the total amplitude of the trend is only about several cm s −1 . In principle, measuring radial velocities separately for shallow and deep lines may therefore be advantageous. However, even if one were able to achieve the desired signal to noise from lines of similar, shallow depths there would still be an RMS at the level of 45 cm s −1 .

Effects of Inclination
As an additional test of GRASS, we examine the relationship between sky-plane stellar inclination and synthetic RV variability. As in the previous sections, we use GRASS to generate many synthetic spectra at the native spectral and temporal resolutions of the input data. These spectra again consist of a single line at 5434.5 A with a continuum-normalized depth of 0.8, similar to the template line in the input data. As seen in Figure 6, the mean RMS of the measured RV time series shows no strong trend with sky-projected inclination. Additionally, we note no strong trend between the widths of the RMS distributions and the stellar inclination. This suggests that the effects of the rotationally-broadened CCF profile for equator-on stars does not significantly affect the amplitude of granulation-induced RV signals. Each distribution is built up from these averaged velocity measurements. The black histograms correspond to simulated observing plans with short, 30-second lags between simulated exposures. The blue histograms correspond to long, 60-minute lags between simulated exposures. The relative values of the measured velocities are themselves arbitrary, and correspond to the difference in rest-frame and CCF-template wavelengths. As expected, the distributions visibly narrow as the number of exposures is increased. However, there is no significant difference in the width of the short-wait and long-wait distributions in all panels (see Table 2), in apparent disagreement with past studies which have shown that long waits between observations produced lower scatter in measured RVs (e.g., Dumusque et al. 2011;Meunier et al. 2015;Collier Cameron et al. 2019). In the future, input data spanning longer than 20 minutes would allow a more robust comparison of observing strategies, since our short baseline input data do not preserve longer timescale variability (e.g., as in supergranulation).

Informing Observing Strategy
In the previous sections, we generated evenly sampled, noiseless time series at the full 15-second cadence of the input solar observations. For real observations of stars other than the Sun, such data are impractical to obtain. Therefore, we use GRASS to simulate more realistic observations to evaluate the effectiveness of various observation strategies at mitigating the impact of granulation.
To generate these spectra, we apply simulations of various observational limitations detailed in §3.2.5 for a variety of synthetic observation patterns. In each observation plan, we simulate a series of 300-second observations (each corresponding to about one full period of the principal solar p-mode). Given practical limitations as well as the ∼20-minute baseline of our input data, we only simulate up to four 300-second observations per night (i.e., 20 total minutes of exposure time) to avoid excessively cycling through the input data. In addition to varying the number of observations, we also vary the time between observations. We use either a short wait time of 30 seconds (corresponding to typical readout times for current-generation optical RV spectrographs), or a long wait of 60 minutes.
As in the previous sections, we generate our singleline synthetic spectra at the native resolution of R = 7×10 5 , but we now convolve the observed spectrum with a truncated Gaussian to degrade the spectral resolution to R ∼ 117, 000 (see §3.2.5). This spectral resolution is comparable to those seen in the current generation of optical RV spectrographs.
We additionally add simulated photon noise to achieve an SNR of ∼5000 per pixel. In realistic observations, the SNR per pixel is often closer to ∼150 (Gupta et al. 2021). Following Murphy et al. (2007), the error in velocity for a single line at this noise level is of order meters per second. In real spectra, this is overcome by the use of many thousands of lines in the CCF computation. For this work, rather than generating much larger spectra, we instead only explicitly synthesize one line. Thus we choose to scale the SNR to the level such that the single-measurement velocity precision is 30 cm s −1 . To verify our choice of SNR, we repeated our analysis with several higher SNR values (not shown) and found no discernible differences in the distributions of measured velocities, suggesting that the impact of granulation noise dominates over photon noise in this high-SNR regime.
For each simulated observation plan, we measure apparent Doppler velocities by fitting to computed CCF profiles. We repeat this procedure a few thousand times in order to sample a distribution of measured Doppler velocities for each synthetic observation plan. The results of this experiment are plotted in Figure 7. As expected, the distributions becoming increasingly narrow with increasing number of observations. However, we note no significant difference in the widths of the short-wait and long-wait distributions for all observation plans. Indeed, simple two-sample Kolmogorov-Smirnov (KS) tests fail to reject the null hypothesis that the short-and long-wait distributions are drawn from the same parent distribution (see Table 2 for a summary of the distribution statistics). This result is not necessarily expected. Previous investigations into empiricallymotivated observations strategies have suggested that returning to a star multiple times over the course of a night (in a manner akin to our long-wait observation simulations) leads to a lower scatter in the measured velocities (e.g., Dumusque et al. 2011;Meunier et al. 2015;Collier Cameron et al. 2019), whereas these distributions suggest no such result.
The short baseline of our input data is likely responsible for this apparent discrepancy. Because our simulations are based on 20-minute observation sequences, GRASS is completely insensitive to longer-timescale phenomena, principally supergranulation. Additionally, to produce the long-wait simulated observations, GRASS must cycle through the same input data some number of times, possibly leading to spurious correlations in the measured velocities. Potential avenues for addressing the limited time-baseline are discussed in §6.2.3.

DISCUSSION AND FUTURE PROSPECTS
This paper presents GRASS, a computational tool for generating observationally-informed stellar spectra with granulation-driven perturbations. We include these perturbations by using spectra of the quiet Sun observed at various positions on the disk to synthesize line profiles. Table 2. Comparison of observing strategies with a short (30-second) or long (one-hour) wait between observations within a night (see Figure 7). Two-sample KS tests fail to reject the null hypothesis that the observed velocity distributions are drawn from the same parent distributions. These perturbations, when interpreted as bulk Doppler velocities (like in CCF analyses), are a significant source of RV noise in planet-search surveys.
GRASS is an important new tool in the suite of alreadyexisting RV-focused star simulators. Whereas tools like SOAP 2.0 (Dumusque et al. 2014) are concerned with magnetic activity and the associated suppression of convective blueshift, GRASS is solely concerned with granulation noise on much shorter timescales. Unlike previous simulations of magnetoconvectively-driven RV noise (e.g., Cegla et al. 2019), GRASS is not based on underlying MHD and radiative transfer simulations. Instead its use of solar line observations allow it to simulate spectra much more quickly.
GRASS will benefit from additional development. In the following sections, we briefly discuss some of these opportunities, as well as some potential applications and limitations of GRASS.

Implications for Observing Strategy
Our initial analyses of GRASS-generated synthetic spectra indicate a correlation between RV RMS and line depth, suggesting that deeper lines contain more information about the range of convective motions present in the solar atmosphere. Our observation simulations, as expected, showed that a greater number of exposures produces a narrower distribution of measured velocities (i.e., a more robust mean). Somewhat unexpectedly, these simulations suggested that the wait-time between observations does not impact the observed scatter in measured velocities. multiple closely-spaced (∼30-second off-target time) exposures produce lower RV scatter than taking multiple distantly-spaced (60-minute off-target time) observations. This result is in apparent disagreement with past studies, which have suggested observations spaced by hours are favorable for reducing the RMS of RVs binned across nights (e.g., Dumusque et al. 2011;Meunier et al. 2015;Collier Cameron et al. 2019). We suggest that the short time baseline of our input data is responsible for this discrepancy. Because we only have 20-minute sequences of input data, GRASS cycles repeatedly through the input data in order to create the long-wait simulations, potentially introducing spurious correlations.
As we mention in §6.2.3, longer baseline input observations could alleviate this limitation and also enable analyses of longer-timescale phenomena, such as supergranulation and meridional flows (see Cegla 2019 for a review). Although Makarov (2010) and Meunier & Lagrange (2020b) find that meridional flows may have large impacts on exoplanet detectability, these motions are variable on very long timescales (i.e., from months to years). It is therefore unlikely that these additional unmodeled motions contribute significantly to the apparent inconsistency of this result with past findings, which consider observations spaced only by hours. Future investigations should also examine the relationship and any correlation between measured RVs and traditionally-used line shape diagnostics (e.g., bisector span, bisector inverse slope, etc.).

Modeling Stellar Spectral Types
GRASS specifically uses solar observations to inform granular variability in spectra. Thus the applicability of this tool for spectral types beyond G2V may be limited. Studies of line bisectors for dwarf stars with convective envelopes have shown that the typical bisector shape is extremely sensitive to spectral type ( Figure 17.15 of Gray 2008). This is especially true for stars hotter than the Sun; as effective temperature increases, the velocity span of the bisectors increases, and the bisectors themselves reverse in curvature ( Figure 17.16 of Gray 2008). This change in bisector velocity span and curvature indicates a change in granulation that may be caused by the thinning of the surface convective zone toward earlier spectral types (Kippenhahn et al. 2012). Toward cooler effective temperatures, typical bisector velocity spans and shapes change less rapidly. Therefore GRASS may be conservatively used to model stars later than G2V, but this should only be done with extreme caution and very strong caveat. In the future, GRASS could use line bisectors produced by MHD simulations of other stellar types as input (Dravins et al. 2021). In this way, we believe GRASS is highly complementary to the suite of Sun-as-a-star MHD simulations.

Modeling Other Lines
In this work, we explicitly validate the efficacy of GRASS for only the Fe I 5434.5Å line. Although the bisectors of shallower lines generally resemble the upper portion of bisectors of deeper lines, there is still some variability between lines. Currently, caution must be taken when simulating synthetic lines in spectral regions far beyond the spectral neighborhood of 5434.5Å. Future iterations of GRASS may include observations of other solar lines. Input data derived from observations of solar lines with very different excitation potentials and/or magnetic sensitivies are two especially important additions that should be made. Following Gray (2008), species of different excitation potentials may have very different lower energy-level abundance distributions and therefore different line profiles. Including some new input data could therefore enable us to model a wider swath of new lines.
In addition to the above variability, the strength of limb darkening also has a wavelength dependence. At shorter wavelengths, limb darkening is stronger, increasing the importance of the bisector shape near disk center. We have also not included the effects of any potential limb brightening, since we do not model any lines for which this is a concern in this work. For future iterations of GRASS that include plages and chromospheric lines, limb brightening would then become an important modeling consideration (Meunier et al. 2010;Dumusque et al. 2014).
We also emphasize that our template line is a magnetically-insensitive line (i.e., it has Landé g-factor g eff = 0.0). Magnetically-sensitive lines, which can exhibit Zeeman splitting and Zeeman enhancement in the presence of a sufficiently strong magnetic field, are known to have their bisectors influenced by magnetic fields (e.g., Cavallini et al. 1988;Brandt & Solanki 1990). These effects are unmodeled by the current version of GRASS.

Input data timescale
The current version of GRASS uses 20-minute baseline input observations. Sampled at a 15-second cadence, these observations are well-suited for studying the short timescales relevant for solar-type granulation, but fail to capture longer-timescale magnetoconvective phenomena. Future versions of GRASS could take advantage of longer baseline input data to include the contribution of longer timescale phenomena, such as supergranulation. To the authors' knowledge, such long-scale data sets do not exist. If such observations are impractical for direct observations, MHD simulations could provide an alternative path to the requisite input data. Although MHD simulations are computationally expensive compared to GRASS, these two methods could be valuable complements to each other. Given a single realization of the appropriate MHD-generated spectra, GRASS could then produce many other realizations at a fraction of the computational cost.

CONCLUDING REMARKS
GRASS is a new tool which uses the first empirical model of convective line-shape perturbations to generate stellar spectra with granular variability, as opposed to previous models which have relied upon more expensive magnetohydrodynamic and radiative transfer simulations. In this work, we have shown that GRASS is a valuable and useful new tool for disentangling stellar RV variability from true Doppler shifts: • GRASS accurately reproduces disk-integrated solar line profiles within ∼0.5% and bisector profiles within a few meters per second ( §4.1 and Figure  3).
• GRASS can probe the effects of line depth, stellar inclination, and other variables on apparent RV variability ( §5.1 and §5.2).
• GRASS in principle is useful for designing and evaluating observation strategies for minimizing noise from granulation ( §5.3), • but in practice, GRASS is currently limited to timescales shorter than 20 minutes, set by the length of the input data which excludes signal from effects such as supergranulation and meridional flows.
As a new entry in the suite of programs designed to model stellar spectra with Doppler variability, we hope that GRASS will serve as a valuable tool in the road toward the detection and characterization of potentially Earth-like planets.  (Hunter 2007)