11institutetext: Department of Astronomy of the University of Geneva, 51
chemin de Pegasi, 1290 Versoix, Switzerland
11email<EMAIL_ADDRESS>
# SOAP-GPU: Efficient Spectral Modelling of Stellar Activity Using Graphical
Processing Units

Yinan Zhao 11 Xavier Dumusque 11

###### Abstract

Context. Stellar activity mitigation is one of the major challenges for the
detection of earth-like exoplanets in radial velocity measurements. Several
promising techniques are now investigating the use of spectral time-series, to
differentiate between stellar and planetary perturbations. In this context,
developing a software that can efficiently explore the parameter space of
stellar activity at the spectral level is of great importance.

Aims. The goal of this paper is to present a new version of the Spot
Oscillation And Planet (SOAP) 2.0 code that can model stellar activity at the
spectral level using graphical processing units (GPUs).

Methods. We take advantage of the computational power of GPUs to optimise the
computationally expensive algorithms behind the original SOAP 2.0 code. For
that purpose, we developed GPU kernels that allow to model stellar activity on
any given wavelength range. In addition to the treatment of stellar activity
at the spectral level, SOAP-GPU also includes the change of spectral line
bisectors from center to limb, and can take as input PHOENIX spectra to model
the quiet photosphere, spots and faculae, which allow to simulate stellar
activity for a wider space in stellar properties.

Results. Benchmark calculations show that for the same accuracy, this new code
improves the computational speed by a factor of 60 compared with a modified
version of SOAP 2.0 that generates spectra, when modeling stellar activity on
the full visible spectral range with a resolution of R=115’000. Although the
code now includes the variation of spectral line bisector with center-to-limb
angle, the effect on the derived RVs is small. We also show that it is not
possible to fully separate the flux from the convective blueshift effect when
modeling spots, due to their lower temperature and thus the appearance of
molecular absorption in their spectra. Rather negligible for the Sun, this
degeneracy between the flux and convective blueshift effect become more
important when we move to cooler stars, however, this issue does not impact
the estimation of the total effect (flux plus convection), and therefore users
can trust this output.

Conclusions. The publicly available SOAP-GPU code allows to efficiently model
stellar activity at the spectral level, which is essential to test further
stellar activity mitigation techniques working at the level of spectral
timeseries not affected by other sources of noise. Besides a huge gain in
performance, SOAP-GPU also includes more physics and is able to model
different stars than the Sun, from F to K dwarfs, thanks to the use of the
PHOENIX spectral library. We however note that due to the limited
understanding of stellar convection and activity on other stars than the Sun,
the more we go away from the solar case, the more the output of the code
should be taken with care.

###### Key Words.:

Methods: data analysis – Techniques: radial velocities – Techniques:
spectroscopic - Stars: activity

##  1 Introduction

The radial velocity (RV) method has been proven to be one of the most
successful method to detect exoplanets since the discovery of the first
exoplanet orbiting a solar-type star (Mayor & Queloz, 1995). In order to
detect earth-like planets orbiting in the habitable zone of its parent star, a
precision of a few dozens of $\rm{cms^{-1}}$ must be reached. Although the
state-of-the-art spectrographs such as ESPRESSO, and EXPRESS are not far from
that precision (50 and 58 $\rm{cms^{-1}}$, respectively Pepe et al., 2021;
Brewer et al., 2020), the main limitation to detect Earth-like planets with
the RV technique is stellar activity. Two major physical processes dominating
stellar activity on a time scale of the host star’s rotational period are the
flux imbalance due to the temperature difference and therefore contrast
between active and quiet regions (hereafter flux effect. (e.g. Saar & Donahue,
1997; Dumusque et al., 2014; Donati et al., 2017)) and the inhibition of
convective blueshift (hereafter CB effect). The CB effect is due to the
presence of strong local magnetic fields inside active regions, which suppress
the CB inside those regions and leads to positive RV variations (e.g.
Cavallini et al., 1985a; Meunier et al., 2010).

Many methods have been proposed to mitigate activity-induced variations using
photometric and spectroscopic time series. In the one-dimensional time series
space, many parametric models based on analytic forms or different Gaussian
process (GP) frameworks have been developed to model stellar activity using
photometry or spectroscopic activity indicators (e.g. Aigrain et al., 2012;
Rajpaul et al., 2015; Aigrain et al., 2016; Gilbertson et al., 2020b; Barragán
et al., 2022). Jointly modeling the data with Keplerians to model planets in
addition to a GP to model stellar activity may significantly reduce the
stellar activity but may also lead to overfitting when the GP kernel or priors
are not wisely set. This is particularly dangerous when the planetary
properties are not constrained from transit observations.

Due to inherent problems in modeling stellar activity in one-dimensional time
series, the community is now shifting toward modeling it in a two-dimensional
space. Collier Cameron et al. (2021) calculated the autocorrelation function
(ACF) of cross-correlation function (hereafter CCF Baranne et al., 1996), to
isolate Doppler shift from shape shift variations and applied principle
component analysis (PCA) on the obtained ACFs to model shape changes related
to stellar activity. A planet signal of amplitude $\sim 40\,\rm{cm/s}$ can be
recovered when the algorithm is applied to the HARPS-N solar data (Dumusque et
al., 2015; Collier Cameron et al., 2019; Dumusque et al., 2021). Zhao et al.
(2022a) projected CCFs time series onto the Fourier basis functions and
modelled line variability using different basis. Results on simulated data
show a 48$\%$ reduction in RV rms. de Beurs et al. (2022) trained a
convolutional neural network (CNN) on both simulated CCFs and HARPS-N solar
CCFs and were able to significantly reduce stellar activity effects.

The idea behind building the CCF is to extract with the best precision the RV
information contained in a spectrum. However, key variations at the spectral
level related to stellar activity may be lost when performing the
dimensionality reduction imposed by the CCF. Therefore, several methods have
been proposed to disentangle stellar activities at the spectral level. Davis
et al. (2017) applied PCA to simulated spectral time series and demonstrated
that eigen-vectors are spectral line dependent. Rajpaul et al. (2020) used GP
to directly derive RV information from spectral time series. Jones et al.
(2017) also applied multivariate GP to model stellar activity on PCA-reduced
spectral dataset. Cretignier et al. (2022), based on the knowledge that the
impact of stellar activity is line-depth dependant (Cretignier et al., 2020a),
used PCA to model stellar activity in the flux-flux gradient space (named the
“shell” space) and results on HD10700 ($\tau$ Ceti) and HD12861 ($\alpha$ Cen
B) indicates the method can successfully remove variations from non-Doppler
origin. Last by not least, Binnenfeld et al. (2020, 2022) are developing the
unit-sphere representation periodogram (USuRPER), to seperate Doppler from
other RV variations. This technique is based on representing spectra as unit
vectors in a multidimensional hyperspace.

The spectral time series used to evaluate the performance of the algorithms
developed to mitigate stellar activity at the spectral level are either
obtained from simulated data or real observations. The major issue with
simulations, is that most of them only model the RV activity effect at the CCF
level (e.g. Dumusque et al., 2014; Herrero et al., 2016) due to computational
inefficiency. A few other libraries of simulated spectra affected by stellar
activity exist, but generating them takes hours to run, which is not
convenient when exploring the parameter space in stellar activity and
properties (e.g. Gilbertson et al., 2020a; Dumusque, 2016). Regarding real
observations, solar data obtained by the HARPS-N solar telescope (Collier
Cameron et al., 2019; Dumusque et al., 2021), HELIOS on
HARPS111https://www.eso.org/public/announcements/ann18033/ and more recently
the solar feed of NEID (Lin et al., 2022) are the best we can get, in terms of
S/N and sampling. However, those spectra corresponds for the most part to
quiet activity phases of the Sun (end of cycle 24 end beginning of cycle 25)
and can only used to mitigate stellar activity for star very similar to the
Sun. When moving to stellar observations, the recent Extreme precision
Spectrograph (EXPRES) Stellar Signals Project (ESSP) shared some valuable
data. However, due to the small number of stars and the rather small number of
spectra available, it was rather difficult to compare different activity
mitigation techniques together (Zhao et al., 2022b). As a conclusion of this
discussion, it is essential for the community to have access to a code that
can simulate efficiently stellar activity at the spectral level, and for a
wide range of stellar properties.

In this paper, we present a new code, Spot Oscillation And Planet Graphical
Process Unit (SOAP-GPU) based on GPU computation that can efficiently model
simplified and realistic stellar activity at the spectral level. In Sect 2, we
revisit the architecture of the SOAP 2.0 code it is based on (Dumusque et al.,
2014) and discuss about its limitations. The algorithms behind SOAP-GPU are
presented in Sect 3. In Sect 4, we explore the physical parameters of stellar
activity and simulation of different cases are presented. Finally, we draw our
conclusion in Sect 5. The SOAP-GPU code is publicly available on Github and
Zenodo222code available here https://github.com/YinanZhao21/SOAP_GPU and
https://doi.org/10.5281/zenodo.7499461 along with a brief manual and some
examples.

##  2 Revisiting SOAP 2.0

In this section, we revisit the code Spot Oscillation And Planet (SOAP 2.0
Dumusque et al., 2014). This code aims at modeling both the flux effect and
the CB effect of active regions affecting RV measurements. Although the public
version of the SOAP 2.0 code can only simulate stellar activity at the level
of the CCFs, modeling the effect at the spectral level follow the same ideas.
In this section, we first discuss the basic algorithms behind SOAP 2.0 and
demonstrate the limit of the code, in terms of computational efficiency, when
we want to model stellar activity at the spectral level.

Figure 1: The stellar disk is initialized with velocity and intensity fields.
_Left:_ The intensity in each cell is computed depending on a limb darkening
law. _Right:_ The velocity in each cell is computed considering rotational
period, stellar inclination and radius of the star. As we can see, iso-
velocity lines are not vertical as we implemented differential rotation in
SOAP-GPU, which was not the case in SOAP 2.0.

###  2.1 The structure of SOAP 2.0

SOAP 2.0 first computes the “quiet” (without any active region) emission
spectrum of the star. To do so, a 2-dimension stellar disk containing N
$\times$ N cells is initialized (N being the resolution of the disk, the same
parameter called “grid” in Boisse et al., 2012). Velocity and intensity of all
disk cells are computed based on the physical configuration of the star
(rotational period, stellar inclination and radius of the star) and a limb
darkening law (as shown in Fig. 1). In each cell, the quiet photosphere
spectrum is injected, weighted by the cell intensity (limb-darkening), and
Doppler-shifted to the projected velocity of that cell (rotation). Linear
interpolation is applied at this step to project the Doppler-shifted spectrum
into the original wavelength grid, to make sure that spectra in different
cells are on a common wavelength grid. We note that the public version of SOAP
2.0 was using the CCF of the high-resolution Kitt Peak Observatory Fourier
Transform Spectrograph (FTS) quiet photosphere spectrum ($S_{quiet}(\lambda)$,
Wallace et al., 1998) as approximation of the quiet Sun to increase
computational speed. However, injecting the original spectrum is possible,
with the only difference that the dimension of the input is $\sim$500000,
compared to 400 for the CCF, and that we will need to apply a Doppler-shift
each time we want to change the velocity of this spectrum, while a simple
translation was sufficient in the case of the CCF. After injecting the quiet
photosphere spectrum in each cells, the integrated quiet solar spectrum is
obtained by summing the content of all the cells together. All those processes
are summarized in the pseudo code below (Algorithm 1).

Algorithm 1 Quiet spectrum integration

1:for $X_{location}=1,2,\ldots N$ do

2: for $Y_{location}=1,2,\ldots,N$ do

3: Shift $S_{quiet}(\lambda)$ with velocity $\mathbf{vel}_{X,Y}$.
$S_{quiet}(\lambda)\rightarrow\ S_{quiet}(\lambda^{{}^{\prime}})$.

4: Do linear interpolation to project the spectrum back to the original
wavelength grid. $S_{quiet}(\lambda^{{}^{\prime}})\rightarrow\
S^{{}^{\prime}}_{quiet}(\lambda)$

5: Weight $S^{{}^{\prime}}_{quiet}(\lambda)$ by limb-darkening intensity
$\mathbf{I}_{X,Y}$ and integrate spectra along disk surface.
$\mathbf{S}_{quiet}+=\mathbf{I}_{X,Y}\times S^{{}^{\prime}}_{quiet}(\lambda)$

6: end for

7:end for

The next step consists in initializing the active regions using the following
parameters: the number of active regions, their size, their corresponding
latitudes and longitudes, their types (either spot or faculae) and the
resolution of the active region contour. An active region spectrum is also
needed at this step to model the CB effect. The original SOAP 2.0 code uses
the CCF of the observed spot spectrum in the visible obtained from the Kitt
Peak Observatory FTS ($S_{active}(\lambda)$ with $\lambda$ the same as for
$S_{quiet}(\lambda)$ Wallace et al., 2005). The spectrum used for faculae
regions is the same, with the difference that the contrast of such a region
follow what is observed in the Sun (e.g. Fig. 3 in Meunier et al., 2010), thus
brighter than the quiet sun and with a center-to-limb brightening. Other
groups use synthetic spectra at different temperature to model the quiet
photospehere, spots and faculae, and include the effect of CB using results
from magneto-hydrodynamical simulations (e.g. the STARSIM 2 code333code
available here https://github.com/rosich/starsim-2 Herrero et al., 2016). We
note that injecting observed or synthetic spectra have their advantages and
drawbacks. Using observed spectra allows to better model the inhibition of
convection inside active regions, but we note that if we use the observed
spectra of a spot to model a facula (because an observed spectrum of a facula
across the entire visible spectral range does not seem to exist), molecular
features will be present in the facula spectrum despite the temperature being
higher. Using synthetic spectra on the contrary allows to better model the
temperature, and therefore spectral features, of the injected spectra. The
choice of the spectra will be addressed in the later sections.

In order to simulate a spectral time series, we need to calculate the disk
location of active regions at each timestamp. As shown in Equation 1 and
Equation 2 of (Boisse et al., 2012), active regions are first put in the
center of the disk and their initial configuration is obtained using a
rotation matrix. Next, to get the position of those active regions as a
function of time, another rotation matrix is used. At each timestamp, the code
evaluates which active regions are visible, and which ones are hidden behind
the star. This is performed by the function
$\textbf{Localize}(lat,long,i,ph)$, where $lat$ and $long$ are the latitude
and longitude of the active region center, $i$ is the inclination angle of the
stellar disk and $ph$ is the rotational phase. The output of this function is
a binary; one if visible, zero otherwise. If a region is visible, the code
proceed to estimate the difference between the quiet solar spectrum and active
spectrum at the location of the active regions.

The difference for the flux effect, the CB effect and the combination of the
two (total) in each cell can be calculated using the following equations:

$\Delta\mathbf{S}^{{}^{\prime}}_{flux}(X,Y)=\mathbf{S}^{{}^{\prime}}_{quiet}(X,Y)-\mathbf{I}_{ratio}\times\mathbf{S}^{{}^{\prime}}_{quiet}(X,Y),$
(1)
$\Delta\mathbf{S}^{{}^{\prime}}_{bconv}(X,Y)=\mathbf{S}^{{}^{\prime}}_{quiet}(X,Y)-\mathbf{S}^{{}^{\prime}}_{active}(X,Y),$
(2)
$\Delta\mathbf{S}^{{}^{\prime}}_{tot}(X,Y)=\mathbf{S}^{{}^{\prime}}_{quiet}(X,Y)-\mathbf{I}_{ratio}\times\mathbf{S}^{{}^{\prime}}_{active}(X,Y),$
(3)

where $\mathbf{I}_{ratio}$ is the contrast of the spot or faculae region. The
code then integrates over all the cells covered by active regions to get final
difference between the quiet spectrum and active spectrum. The final spectrum
of each effect at each timestamp can be calculated by:

$\mathbf{S}_{integrated,final}=\mathbf{S}_{integrated,quiet}-\Delta\mathbf{S}_{integrated,quiet-
active}.$ (4)

Once a spectrum for each timestamp is obtained, the code then lowers the
resolution of the integrated spectrum to match the resolution provided in the
configuration file. The pseudo code that describes how active regions are
included, and how the final integrated spectra is obtained is summarized
below.

Algorithm 2 Active region updates

1:for $n_{region}=1,2,\ldots N$ do

2: for $t_{timestep}=1,2,\ldots,T$ do

3: $\textbf{Localize}(lat,long,i,ph)$.

4: if $\textbf{Localize}=1$ then

5: Shift $S_{quiet}(\lambda)$ with velocity $\mathbf{vel}_{X,Y}$.
$S_{quiet}(\lambda)\rightarrow\ S_{quiet}(\lambda^{{}^{\prime}})$.

6: Do linear interpolation to project the spectrum back to the original
wavelength space. $S_{quiet}(\lambda^{{}^{\prime}})\rightarrow\
S^{{}^{\prime}}_{quiet}(\lambda)$

7: Weight $S^{{}^{\prime}}_{quiet}(\lambda)$ by limb-darkening intensity
$\mathbf{I}_{X,Y}$

8: Shift $S_{active}(\lambda)$ with velocity $\mathbf{vel}_{X,Y}$.
$S_{active}(\lambda)\rightarrow\ S_{active}(\lambda^{{}^{\prime}})$.

9: Do linear interpolation to project the spectrum back to the original
wavelength space. $S_{active}(\lambda^{{}^{\prime}})\rightarrow\
S^{{}^{\prime}}_{active}(\lambda)$

10: Weight $S^{{}^{\prime}}_{active}(\lambda)$ by limb-darkening intensity
$\mathbf{I}_{X,Y}$

11: Use Equations 1 to 3 to calculate the difference of each effect

12: end if

13: Compute summation of $\Delta\mathbf{S}^{{}^{\prime}}(X,Y)$ for each
effect.

14: end for

15:end for

16:for $t_{timestep}=1,2,\ldots,T$ do

17: for $n_{region}=1,2,\ldots N$ do

18: Use Equation 4 to update final spectrum at $t_{T}$.

19: Lower the resolution of final spectrum at $t_{T}$ to match the HARPS-N
observation.

20: end for

21:end for

###  2.2 The limitation of SOAP 2.0

The structure of SOAP 2.0 provides an efficient way to estimate stellar
activities on spectroscopic measurement by simulating CCFs at different
timestamps. The major drawback when changing the input from CCFs to spectra is
the dimension of the data. The dimension of the input CCFs in SOAP 2.0 was
$400$ in velocity space while the input high-resolution spectra we want to use
have a dimension of $\sim 500000$ in the wavelength domain. From Algorithm 1
and Algorithm 2, we clearly see that the linear interpolation is repeatedly
called when injecting the spectrum in each cell, which is computationally
expensive for an array with dimension of $\sim 500000$. For example, SOAP 2.0
takes $\sim 800$ seconds to calculate an integrated quiet sun spectrum using a
$300\times 300$ disk-grid. Another issue is how the code handles multiple
active regions. Each active region is modeled independently, without
information from other regions. This algorithm cannot handle the case in which
some active regions overlap with each other. From real observations, we know
that some active regions have complicated configurations. For example, most of
active regions are a combination of a large faculae presenting a small spot in
its center. In this context, a more computationally efficient and generalized
algorithm is needed.

##  3 Description of SOAP-GPU

In the previous section we’ve demonstrated the limitation of SOAP 2.0 when
modeling stellar activity at the spectral level. Here, we present a new
version of SOAP, based on Graphical Processing Unit (GPU) computing, that is
much more efficient in term of computational speed, but also that adds some
physical complexity.

###  3.1 The basic concept of GPU computing

The popularity of artificial intelligence has in recent year significantly
increased due to the programmability of graphic hardwares. GPU computing uses
graphical card as a co-processor for parallel computing. Compared with CPU,
GPU solves problems by breaking them into separate tasks and processing them
simultaneously. The basic computational unit that can independently perform
simple calculation in a graphic card is called a $thread$. A group of threads
that communicate and share memory with each other is called a $block$.

The new version of SOAP presented here, SOAP-GPU, is written using the Compute
Unified Device Architecture (CUDA). CUDA is a compiler and toolkit for
programming NVIDIA GPUs, and is an extension of the C/C++ programming
language. CUDA invokes kernel functions by using the syntax of
$<<<N_{blocks},N_{threads}>>>$. This syntax allows the user to define the
thread hierarchy before launching in parallel the same program function called
$kernel$ to many threads. In order to launch the computation at the level of
the GPU, a $host$ function defined in CPU controls the data transfer between
CPU and GPU and can execute the kernel function inside the GPU.

Threads in the same block can be accessed as 1D, 2D or 3D structures. In order
to perform the thread level calculation, the index of individual thread and
block need to be accessed. The index of each thread in the same block can be
expressed as $threadIdx$. If the block is launched as the 1D structure, each
thread in the same block can be accessed as $threadIdx.x$. The number of the
treads used in each 1D block can be obtained as $blockDim.x$. $Grid$ is a
group of blocks. It can be either 1D, 2D or 3D. For the 1D grid, the index of
each block in the grid can be expressed as $blockIdx.x$. Since the input
spectra of quiet sun and active region are both 1D, we used the configuration
of 1D grid with 1D block and the global index is
$index=blockIdx.x*blockDim.x+threadIdx.x$.

###  3.2 Fast linear interpolation with GPU

As mentioned in previous sections, the major limitation in SOAP 2.0 is the way
it handles linear interpolation in each disk cell. A GPU provides thousands of
cores which can be implement for linear interpolation for large data array.
Considering that both quiet sun and spot spectra are evenly sampled in the
wavelength domain, then the input wavelength can be described as:

$\lambda_{n}=\lambda_{0}+nk,$ (5)

where $n$ is the pixel number and $k$ is the step size. When a Doppler shift
is applied, the wavelength array is modified as follow:

$\lambda^{{}^{\prime}}_{n}=\lambda_{n}+\lambda_{n}f(\beta),$ (6)

where $f(\beta)=-\left[1-\sqrt{\frac{(1+\beta)}{(1-\beta)}}\right]$ and
$\beta=v/c$. The variable $v$ is the velocity for each cell and $c$ is the
speed of light. Since we need to project the shifted spectrum back to the
original wavelength space $S(\lambda^{{}^{\prime}})\rightarrow\
S^{{}^{\prime}}(\lambda)$, we have to find the index $m$ which satisfies
$\lambda^{{}^{\prime}}_{n}<\lambda_{m}<\lambda^{{}^{\prime}}_{n+1}$. For the
left side, we have:

$\lambda^{{}^{\prime}}_{n}<\lambda_{m},$
$\lambda^{{}^{\prime}}_{n}=\lambda_{n}+\lambda_{n}f(\beta)=\lambda_{0}+nk+\lambda_{0}f(\beta)+nkf(\beta)<\lambda_{0}+mk,$

so we have:

$n(1+f(\beta))+\frac{f(\beta)\lambda_{0}}{k}<m.$ (7)

For the right side, we have:

$m<(n+1)(1+f(\beta))+\frac{f(\beta)\lambda_{0}}{k}.$ (8)

Once the integer $m$ is known, we can estimate the flux for $\lambda_{m}$
using the spectrum derivative:

$S^{{}^{\prime}}_{m}=\frac{\Delta
S_{n}}{(\lambda^{{}^{\prime}}_{n+1}-\lambda^{{}^{\prime}}_{n})}\times(\lambda_{m}-\lambda^{{}^{\prime}}_{n})+S_{n},$
(9)

where $S_{n}=S(\lambda_{n})$ and $\Delta S_{n}=S_{n+1}-S_{n}$.

Equations 7 to 9 can be parallelised using GPU. We launch 1D grid of 1D blocks
with $<<<N_{blocks},N_{threads}>>>$ to perform the linear interpolation
mentioned above and the number of blocks and threads satisfies
$Dim_{input\\_spectrum}=N_{blocks}\times N_{threads}$. The pseudo code for
this part is summarised in Algorithm 3 and the quiet sun spectra integration
can be rewritten as Algorithm 4.

Algorithm 3 Fast interpolation with GPU

1:$index=blockIdx.x*blockDim.x+threadIdx.x$

2:$index_{target}=ceil(index*(1+f(\beta))+f(\beta)*\lambda_{0}/k)$

3:$S^{{}^{\prime}}_{index_{target}}=\frac{\Delta
S_{index}}{(\lambda^{{}^{\prime}}_{index+1}-\lambda^{{}^{\prime}}_{index})}\times(\lambda_{index_{target}}-\lambda^{{}^{\prime}}_{index})+S_{index}$.

Algorithm 4 Quiet spectrum integration with GPU

1:for $X_{location}=1,2,\ldots N$ do

2: for $Y_{location}=1,2,\ldots,N$ do

3: Apply Doppler shift with velocity $\mathbf{vel}_{X,Y}$ and derive
$S^{{}^{\prime}}_{quiet}(\lambda)$ using GPU fast interpolation

4: Weight $S^{{}^{\prime}}_{quiet}(\lambda)$ by limb-darkening intensity
$\mathbf{I}_{X,Y}$ and integrate spectra along disk surface.
$\mathbf{S}_{quiet}+=\mathbf{I}_{X,Y}\times S^{{}^{\prime}}_{quiet}(\lambda)$

5: end for

6:end for

###  3.3 Active region updates

As addressed in the previous section, one of the disadvantage of SOAP 2.0 is
that each active region is modeled independently, which makes the code unable
to handle complicated active region configurations: some active regions may
overlap with each other; spots may be surrounded by facualae regions. Here we
propose a revised algorithm to update active regions: an empty disk map called
$InfoMap$ is allocated in the GPU first. At each timestamp, A list of active
regions with their properties is uploaded and the code calculates the location
of active regions projected on the disk map. If some regions are visible, we
update the corresponding pixels with their active region types in the
information map. For example, if a faculae region is visible at
($x_{n},y_{n}$), $InfoMap(x_{n},y_{n})=1.0$. If there are multiple regions
with the same type overlapping with each other, the overlapping region in the
information map will remain the same. This will avoid the over-calculation for
the overlapping region issue in the SOAP 2.0 since each acitve region is
calculated independently. This algorithm can also simulate complicated active
region configurations. For example, a spot surrounded by a large faculae can
be simulated by updating the information map with a faculae first. If the spot
region is embedded inside the faculae, the overlapping region in the
information map will be updated with the type of the spot. The pseudo code of
this part is summarised in Algorithm 5.

Algorithm 5 Active region updates with GPU

1:for $t_{timestep}=1,2,\ldots,T$ do

2: for $n_{region}=1,2,\ldots N$ do

3: $\textbf{Localize}(lat,long,i,ph)$.

4: if $\textbf{Localize}=1$ then

5: Updating $InfoMap(x_{n},y_{n})=$ the type of active region.

6: end if

7: end for

8: Inject velocity $\mathbf{vel}_{X,Y}$ with GPU fast interpolation for the
active regions in the information map. and derive
$S^{{}^{\prime}}_{active}(\lambda)$.

9: Weight $S^{{}^{\prime}}_{active}(\lambda)$ by limb-darkening intensity
$\mathbf{I}_{X,Y}$

10: Use Equations 1 to 3 to calculate the difference of each effect

11: Compute summation of $\Delta\mathbf{S}^{{}^{\prime}}(X,Y)$ for each
effect.

12: Use Equation 4 to derive the final spectrum at $t_{T}$.

13: Lower the resolution of final spectrum at $t_{T}$ to match the HARPS-N
observation.

14:end for

###  3.4 Differential rotation

In the original SOAP 2.0 code, there is no differential rotation implemented.
In order to better model stellar activity, differential rotation is included
when the stellar disk is initialized according to the equation
$\omega=\omega_{0}+\omega_{1}\sin^{2}(\theta)$, where
$\omega_{0}=14.371^{\circ}/\rm{day}$ and $\omega_{1}=-2.587^{\circ}/\rm{day}$
for the Sun (Borgniet et al., 2015). To generalise this for other stars, the
user can select in the configuration of SOAP-GPU a rotation period and a
differential rotation rate. $\omega_{0}$ is then equal to 360/PROT and
$\omega_{1}$ to DIFF_ROT*PROT (PROT=25.05 and DIFF_ROT=-0.18 for the solar
case to reproduce the above equation).

##  4 Results

Figure 2: The computation speed comparison between SOAP 2.0 and SOAP GPU: The
integrated quiet sun spectrum is calculated with different disk resolutions.
When disk resolution is below 10, SOAP 2.0 is faster than SOAP-GPU since the
communication between CPU and GPU in SOAP-GPU is time-consuming. For
resolutions above 10, SOAP-GPU is significantly faster than SOAP 2.0. With a
typical resolution value of 300, the quiet disk spectrum integration in SOAP-
GPU is 100 times faster than in SOAP 2.0

###  4.1 Performance and precision comparison with SOAP 2.0

We examined the performance of SOAP-GPU in two aspects: computational speed
and accuracy. SOAP-GPU code is executed on a Nvidia RTX-3090 card while we run
the modified SOAP 2.0 that generates spectra in a MacBook Pro with 2.6 GHz
6-Core Intel Core i7. We analysed the speed performance of SOAP-GPU by
calculating the time it takes to obtain an integrated quiet sun spectrum. The
input quiet sun spectrum has a dimension of 547840, thus the kernel function
fast interpolation is launched with
$<<<N_{blocks},N_{threads}>>>=<<<1070,512>>>$. We note that $N_{threads}$ is
fixed to 512 and $N_{blocks}$ is an adaptive number based on the dimension of
the input. SOAP 2.0 is executed with the same simulation configuration on a
single CPU. We computed the integrated quiet sun spectrum with different disk
resolution and their computational time is show in Figure 2. When the disk
resolution is very low, smaller than 10, SOAP 2.0 is faster than SOAP-GPU.
This is not surprising since the data transfer between GPU and CPU in SOAP-GPU
is the dominating factor. When the disk resolution increases, SOAP-GPU is
significantly faster than SOAP 2.0. When the resolution is above 100, the
quiet sun spectrum integration of SOAP-GPU is 100 times faster than SOAP 2.0
and both computational curves linearly increase in log-log space.

Boisse et al. (2012) found no significant change in their results with
resolution beyond 300, therefore, we used this disk resolution for the
following of the paper. For the typical disk resolution of 300, a spot at disc
center with an area of $1\%$ of the entire disk will be contained in a grid of
$34\times 34$ cells. Due to the small size of the grid for such a
configuration, the fast interpolation algorithm (see Algorithm 3) is only able
to gain a factor of $\sim$10 in computation time. If the spot size increases
to $9\%$ of the entire disk, the simulation can then gain almost the full
speed boost from fast interpolation (100 time faster). Fast interpolation at
the level of the active region modelisation makes therefore significant
improvements in computational speed when considering high-resolution
simulations or simulations with large active regions.

Figure 3: Comparison of the RVs derived from the simulated spectra modeled by
SOAP2.0 and SOAP-GPU. A single equatorial spot with 1% area of the entire disk
surface is simulated. It took 1749.3 seconds to simulate those spectra with
SOAP 2.0 while only 27.9 seconds with SOAP-GPU on a Nvidia RTX-3090 card. The
computation speed is improved by a factor of 63.

We also examined the accuracy of SOAP-GPU. We simulated a single equatorial
spot with a $1\%$ area of the entire disk surface using a disk resolution of
300. It took 1749.3 seconds to simulate those spectra with SOAP 2.0 for 100
timestamps while only 27.9 seconds using SOAP-GPU, which corresponds to a gain
of a factor 63. The modeled RVs relative to the flux effect, the CB effect and
the total effect are derived from the simulated spectra by cross-correlating
them with the same mask originally used in SOAP 2.0, and measuring the RV as
the mean of a Gaussian profile fitted to the obtained CCFs. Figure 3
illustrates that the simulated spectra from SOAP-GPU provides the same RVs as
the spectra from SOAP 2.0.

###  4.2 Exploration of active region properties

The dynamics of active regions plays an important role for understanding the
stellar activity-induced RVs. Most of previous study aimed at investigating
these effects with real observations. For example, Meunier et al. (2010)
derived the stellar activity induced RVs by using Michelson Doppler
Imager/Solar and Heliospheric Observatory (MDI/SOHO) magnetograms images. At
the simulation level, Gilbertson et al. (2020a) investigated the effect of
spot evolution on the long-term and at the spectral level, using a modified
version of SOAP 2.0. However, they only considered spots, and only their
decaying phase. In order to illustrate the effects of active region dynamics,
we discuss in this section the photometric and RV variations observed when an
active region changes in size, when different number of active regions are
present and when the active region configuration changes.

####  4.2.1 The size evolution of active region

To explore different active region evolution scenarios, we developed and
included an evolution module in SOAP-GPU. This module can model evolution in
three different ways: i) a linear growing phase, ii) a linear decaying phase
or iii) a growing and decaying phase modeled by an an asymmetric Gaussian
function (Muraközy et al., 2014). Other user-defined functions can be added to
this module if desired. We show in Fig. 4 the impact of active region
evolution on the light-curve and on the different RVs derived (flux, CB and
total effects). For the asymmetric Gaussian evolution phase, the maximum size
is set to 10000 millionths of solar hemisphere (MSH) equivalent to 1% of the
visible hemisphere, the FWHM to 10 days and an asymmetry factor of 0.09. For
the growing only, or decaying only evolution phases, the initial size is set
to 10000 MSH and the growth or decay rate is set to 400 MSH/day. We found that
both flux and CB effects are sensitive to the evolution of active regions.

Figure 4: SOAP-GPU simulation of different active region evolution curves: A
single spot with a latitude of 30°and longitude of 180°is simulated. Three
spot size evolution types are demonstrated: i) a fast growth and slow decay
evolution is shown in the first column, ii) a linear decay evolution curve in
the second column and iii) a linear growth curve in the third column. The
evolution curves and the simulated light curves are shown in the first two
rows. The RVs of the total effect, the flux effect and the CB effect are
present in the rest of the rows. The simulation of an non-evolving spot (red
dashed line ) is also shown in each figure for comparison.

####  4.2.2 Complex active regions

SOAP-GPU also allows users to simulate complex active region configurations.
From the observational point of view, facluae and spots are not independent
from each other. The facula distribution is based on the spot distribution.
This leads to a complex configuration in which spots may overlap faculae
(Borgniet et al., 2015; Chapman et al., 2001). In order to model such a
configuration, the SOAP-GPU config file allows users to define the
distribution of active regions, as a sequence of spots and faculae with given
properties (size, initial longitude, initial latitude). For example, in order
to simulate a spot surrounded by a facula, the user can define the location
and the size of the large facula first and then define a smaller spot at the
same location. The region of overlap will be replaced by the spot as mentioned
in Algorithm 5. A simulation of this case is illustrated in Figure 5, with a
$1\%$ spot surrounded by a $9\%$ facula. Since the spot has a higher contrast
than the faculae, the light curve and RVs of the flux effect is dominated by
the spot while the RVs of the CB effect is dominated by facula. Overall, the
CB RV effect induced by the facula dominates all the other contributions, and
thus the total RVs is affected mainly by the facula, as it was already
demonstrated in several studies (e.g. Meunier et al., 2010; Dumusque et al.,
2014; Milbourne et al., 2019).

Figure 5: Three different active region configurations are simulated. A single
spot located at the latitude of 30°and the longitude of 180°with a fixed size
of $1\%$ of the entire solar disk is present in black. A single facula with
the same coordinates and a fixed size of $9\%$ of the entire solar disk is
shown in blue line. The $1\%$ spot region surrounded with $9\%$ faclua region
is labeled in red dashed line. The top panel demonstrates the light curves of
different configurations and the rest of the panels shows the RVs of the total
effect, the flux effect and the CB effect.

###  4.3 Exploration of spectral properties

In this section, we explore the input spectra properties and demonstrate how
the derived RV behaves depending on the wavelength domain.

####  4.3.1 Chromatic effects of different wavelength coverage

To explore the effect induced by different wavelength coverage, we injected
into SOAP-GPU only the red or only the blue part of the quiet sun and spot
spectra (see Fig. 6). The red and blue parts have the same dimension of
204800, which is different from the full spectra. As the fast interpolation
kernel function depends on the dimensions of the input spectra, the code
automatically configures the kernel with the option $<<<400,512>>>$.

Figure 6: Injecting different size of spectra in SOAP-GPU input. Three
different sets of quiet sun and spot spectra, with different wavelength ranges
are used as input to SOAP-GPU: The entire spectra with length 547840 is
labeled in black. The blue and red parts of the spectra with length 204800 are
over plotted in blue and red, respectively. Figure 7: The chromatic effect of
RVs. SOAP-GPU is initialized with three different spectra: entire wavelength
coverage, and only red and blue spectral parts (see Sect. 6). A single spot
with a size of 1%, a latitude of 30°and a longitude of 180°is modeled by the
code. The measured RVs with different input spectra are labeled with black,
red and blue, respectively. Top: The measured RVs of the total effect. Middle:
The measured RVs of the flux effect. An offset of 1 $\rm{ms}^{-1}$ is added in
the red and blue RVs. Bottom: The measured RVs of the CB effect.

The measured RVs are shown in Figure 7. The RVs of the CB effect are different
between the blue and red parts. One notable thing is that the RVs of the CB
effect simulated from the red inputs goes below zero, while we expect the CB
effect to only be positive, as it corresponds to an inhibition of CB. To
confirm that nothing was wrong at the level of the code, we injected for the
spotted region the same spectrum as the quiet Sun, but we red-shifted it by
300 m/s to model at first order the inhibition of CB. In this idealist case,
the CB effect does not provide negative values. After further investigation,
those negative values comes from the fact that the spot temperature is lower
than the quiet photosphere. Thus, spectral lines will change in depth, which
will induce a flux effect even when not considering the contrast of the active
regions when estimating the CB effect (see Eq. 2). In the case of the Sun,
this flux effect seen in the CB derived RVs is mainly coming from the red part
of the spectrum due to molecular absorption that can be seen in the spot
spectrum, but not in the quiet photosphere spectrum. As we will see in Sect.
4.3.3, we do not obtain negative values when injecting PHOENIX solar
equivalent spectra for the quiet and active Sun instead of the Kitt Peak solar
quiet and active atlases, however we still see a slight asymmetry in the
derived CB RV effect, pointing toward a small flux effect contribution. This
is likely because PHOENIX spectra are not able to model all the absorptions
coming from molecular bands, and thus the flux effect seen in the CB
estimation only is stronger for the real solar spectra than for the synthetic
spectra. This issue prevent us of fully separating the flux from the CB
effect, however, we note that the total effect (flux + CB) should be modeled
properly. We note that this feature was not visible in SOAP 2.0, as after
computing the CCF for the quiet and actives regions, we were renormalising
them.

We note that in SOAP 2.0, we used a fixed contrast to model the flux effect of
active regions. This contrast was derived by comparing the Planck function of
the quiet Sun effective temperature and of the spot or facula
temperature444the config file in SOAP 2.0 allows to give an effective
temperature for the quiet photosphere, 5778 $K$ for the Sun, and a temperature
difference with respect to this former value for the spot spectrum (663 $K$ as
the default value in SOAP 2.0). For a facula, the temperature is dependent on
the center-to-limb angle and was following what is observed on the Sun (e.g.
Fig. 3 in Meunier et al., 2010)., at the average wavelength of the input
spectra 5293 $\AA$. Now that we use spectra as input, and not CCFs, we
implemented a contrast that is wavelength dependant to model the chromatic
effects of stellar activity. To do so, we introduced a new GPU kernel function
called SOAP$contrast\,<<<N_{blocks},N_{threads}>>>$. This new kernel allows to
perform the wavelength dependent contrast calculation: each wavelength pixel
is first accessed by the global index of the kernel. Next on each thread, it
derives the contrast by calculating the ratio of two Planck functions at two
different effective temperatures. The absolute value of the contrast in the
blue part of the spectrum is higher than in the red part. This implies that
the flux effect for the blue part is stronger than in the red part, which can
be seen in the middle panel of Fig. 7.

####  4.3.2 Convection as a function of center-to-limb angle

Solar spectral line profiles become asymmetric due to convective motions
varying with physical depth inside the solar photosphere (e.g. Dravins et al.,
1981; Gray, 2009). This effect also leads to a change in shape of the bisector
of spectral lines from disk-center to the limb, as photons are coming from
different physical depths (e.g. Cavallini et al., 1985b). In order to better
model the effect of convection in SOAP-GPU, we derived this effect from very-
high spatial and spectral resolution observations of the Sun (Löhner-Böttcher
et al., 2018; Stief et al., 2019; Löhner-Böttcher et al., 2019).

We note that the varying shape of spectral line with center-to-limb $\mu$
angle is also modelled in the STARSIM 2 code Herrero et al. (2016) by fitting
a fourth-order polynomial function on magneto-hydrodynamic CIFIST 3D models
(Ludwig et al., 2009). However, this fifth-order polynomial is only valid for
line depth as strong as $\sim 0.5$ (see Fig. 6 in Herrero et al., 2016)555The
parameters mentioned here are derived from the published code (
https://github.com/rosich/starsim-2). This was enough to model the shape
change of the CCFs in STARSIM 2, which does not go deeper. However, when
working at the spectral line level, this polynomial will give completely wrong
estimate for the core of deep lines, due to extrapolation. We therefore used
the quiet sun observations at different $\mu$ angles provided in Löhner-
Böttcher et al. (2019). We first measured the bisectors of all the available
iron deep lines at different $\mu$ angle in Löhner-Böttcher et al. (2019)666We
use the following lines: FeI $5250.2084\AA$, FeI $5250.6453\AA$, FeI
$5434.5232\AA$, FeI $5432.9470\AA$, FeI $5576.0881\AA$, FeII $6149.2460\AA$,
FeI $6173.3344\AA$, FeI $6301.5008\AA$ and FeI $6302.4932\AA$ , and fitted
them using polynomial functions. For $\mu$ angle smaller then 0.5, we used a
straight line to fit the bisectors of the selected deep lines, to prevent
strong divergence when extrapolating the fit towards very large depths. For
$\mu$ angle larger or equal to 0.5, 3rd order polynomial functions are used to
capture the curvature of the bisectors around disk center. The final
bisectors, shown in Fig. 8 are obtained by interpolating and extrapolating
those bisectors from depth 0 to 1.

To obtain the dependency of the spectral line bisector as a function of $\mu$
in an active region, we use the observations presented in Cavallini et al.
(1985a). We parameterised the bisectors of the FeI at $6301.5008\AA$, for the
disk center ($\mu=1$) and different center-to-limb angles ($\mu=0.82$, 0.66
and 0.44). The bottom part of the bisectors, below 0.5, is fitted using a
straight line, the upper part for which we have data, using a 5th-order
polynomial. Rather than extrapolating the fitted polynomial towards very
shallow depths, which can give unrealistic redshifted values, we decided to
use the more redshifted data value of the top bisector for extrapolation. We
show in Fig. 8 the obtained active bisectors from depth 0.0 to 1.0.

Once we have our model for line bisectors at different $\mu$ angles, we can
use the Python module $Convec.model$ to apply those bisectors to the original
spectra, and thus obtain different spectra for different $\mu$ angles. Each
cell in the stellar disk takes the bisector that has the closest $\mu$ angle.
However, the code first has to remove the original bisector from the spectral
lines of the input quiet and active Kitt Peak solar spectra. To do so, we
select the same lines as in Löhner-Böttcher et al. (2019) in the input spectra
and measure their individual bisectors. To model the average bisector of the
lines selected in the quiet spectrum, we use a second-order polynomial. For
the active spectrum, due to the lower effective temperature, the wings of
certain lines fitted are blended, which strongly impact the bisector
measurement. We therefore rejected bisector points that are significantly off.
Then, we model the average active bisector of the lines by fitting the regions
below and above a depth of 0.5 with two different linear models. Fitting a
higher-order polynomial for those active bisectors was giving unrealistic
values when extrapolating to very small or very large depths. The measured
individual bisectors with our models are shown in Fig. 9. To finally obtain
quiet and active spectra with proper bisector shape as a function of $\mu$
angles, we remove the original bisectors of the quiet and active Kitt peak
solar spectra, and then add the bisectors measured for different $\mu$ angles.
This is done by shifting each point in those spectra depending on their
normalised depth.

Figure 8: Average bisectors of quiet and active solar regions from the disk
center ($\mu=1.0$) to the limb ($\mu=0.2$). _Continuous lines:_ Fifth-order
polynomial fit to the quiet sun bisectors of the FeI $5250.2084\AA$, FeI
$5250.6453\AA$, FeI $5434.5232\AA$, FeI $5432.9470\AA$, FeI $5576.0881\AA$,
FeI $6149.2460\AA$, FeI $6173.3344\AA$, FeI $6301.5008\AA$ and FeI
$6302.4932\AA$ lines as measured by the Laser Absolute Reference Spectrograph
(LARS) at the German Vacuum Tower Telescope (Löhner-Böttcher et al., 2019).
_Dashed lines:_ Fit of the bisectors of the FeI $6301.5008\AA$ spectral line
inside a faculae region, as measured by the Fabry-Perot interferometer at the
Donati Solar Tower (Cavallini et al., 1985a). Below a depth of 0.5, a linear
fit is performed, while a fifht-order polynomial is used to model the top part
of the bisector. To prevent unrealistic value when interpolating the
polynomial above a normalised flux of 0.9 where no measurement exists, we
selected the most redshifted part of the top bisector, explaining the vertical
values for very shallow depths. The two vertical lines are shifted by 340 m/s
which corresponds to the solar convective blueshift value derived from a fit
to the data of Liebing et al. (2021) (see Sect. 4.3.3). The active bisectors
at different $\mu$ angles are all shifted by those 340 m/s at a depth of 0.58
as we make the hypothesis that convection is fully suppressed in magnetic
regions (see Sect 4.3.2 for more information). Figure 9: Bisectors of the FeI
lines used in Löhner-Böttcher et al. (2019) and fitted model to account for
the CB. _Left:_ Bisectors from the quiet Kitt Peak solar spectrum. _Right:_
Bisectors from the active Kitt Peak solar spectrum. We rejected the bottom
part of the $6301.5008\AA$ bisector because it was significantly off by 2500
m/s due to strong contamination by other weak lines.

To inject the proper bisectors for different $\mu$ angles in our original
spectra, we first remove the original bisector, which changes any CB
difference between the quiet and active solar spectra. We therefore need to
impose a shift between the bisectors of quiet and active regions in order to
properly model the inhibition of CB inside active regions. We here make two
assumptions: i) the CB is fully inhibited at $\mu$=0.2 in the quiet Sun, and
ii) it is also fully inhibited for magnetic regions, and this at all $\mu$
angles. Using the first assumption, we measure for the quiet Sun the maximum
shift between the bisector at $\mu$=0.85, which is the bisector that is the
most blueshifted, and the bisector at $\mu$=0.2. This maximum happens at a
depth of 0.58 and equals to 375 m/s. This value is extremely similar to the
340 m/s CB value derived from a fit to the data of Liebing et al. (2021) (see
Sect. 4.3.3). To match the CB relation derived in Sect. 4.3.3, we rescale the
maximum difference between the $\mu=0.85$ and $\mu=0.2$ to be 340 m/s. We note
that this rescaling is negligible in the case of the Sun, however, it will be
really needed in Sect. 4.3.3 when using PHOENIX spectra as input. Using the
second assumption, we impose that at the same depth of 0.58, the difference in
velocity between the quiet bisector at $\mu$=0.85 and all active bisectors is
also 340 m/s. We show the proper shift between the quiet and active bisectors
in Fig. 8.

We show the RV impact of considering the $\mu$ angle dependency on the
observed solar spectra in Fig. 10. As we can see, the impact is not
significant when looking at the shape of the signal as a function of phase.
This come from the fact that due to limb-darkening and the projection of
active regions on the limb, most of the signal comes from larger $\mu$ angles
(close to disc center). The only significant difference is for the amplitude
of the CB effect. This is because we forced the CB difference between the
quiet and active sun to be 340 m/s, while the CB difference between the quiet
and active Kitt peak solar spectra is less than 300 m/s when measuring the
average difference between the quiet and active CCF bissectors (see Fig. 2 in
Dumusque et al., 2014). We note that the complexity of modifying the bisectors
depending on the center-to-limb angle is not strongly justified when using
real solar spectra as input due to the small difference observed in the
estimated RVs, however, this step is critical when working with synthetic
spectra that does not include the proper bisectors, as described in Setc.
4.3.3.

We are conscious that depending on the magnetic field of an active region, the
inhibition of the CB will be different and therefore the bisectors more or
less redshifted compared to the quiet Sun, as seen in Fig. 1 in Cavallini et
al. (1985a). Also, faculae tend to have weaker magnetic fields than spots and
in our case, we model those two active regions with the same bisectors and the
same CB inhibition. It is therefore likely that the CB effect for faculae is
slightly overestimated, and this will translate in larger RV amplitudes when
modeling the CB effect for faculae.

In summary, in this subsection we present a framework to model the Sun but
also other stars (see also next subsection). Different bisectors at different
$\mu$ are derived from the quiet phototsphere (Löhner-Böttcher et al., 2019)
and facuale (Cavallini et al., 1985b) and are injected into the input spectra
for which we have removed any variation in line bisector from a vertical line.

Figure 10: Comparison of the RVs derived using two different configurations
for the input spectra. A single equatorial spot with 1% area of the entire
disk surface is simulated in both cases. Red line: RVs derived from simulated
spectra with observed quiet sun and spot spectra without $\mu$ dependent
bisector injection. Dotted line: RVs derived from simulated spectra with $\mu$
dependent bisector injection. The input spectra with $\mu$-angle dependency
are generated with the Python module $Convec.model$. The RVs of the flux, CB
and total effect don’t change significantly when the $\mu$-dependent CB is
introduced.

####  4.3.3 Simulation based on the PHOENIX spectral database

The implementation of convective motions described in the preceding section
allow us to use synthetic spectra as input, since the effect of convection can
be injected using the $Convec.model$ module. In order to study stellar
activity affecting the data used in RV, a high resolution spectral library is
needed. For SOAP-GPU, we decided to make it easy for the user to use as input
PHOENIX high-resolution spectra (Husser et al., 2013). We note however that
SOAP-GPU can accept other spectral libraries, but it might be a little more
difficult for the users to properly setup the inputs since the parameters to
remove bisectors of input spectra are only optimized for the PHOENIX spectra
and the solar atlas from the Kitt Peak Observatory FTS (Wallace et al., 2005).
The PHOENIX library propose a collection of spectra with the wavelength
coverage from $500\AA$ to $5.5\mu m$ with resolutions of 500,000 in the
optical. The library covers stellar effective temperature from 2300K to
12000K. Since the spectra in the PHOENIX library are not normalized, which is
critical to perform the injection of CB described in the preceding section, we
used the open source package “Rassine” (Cretignier et al. (2020b)) to perform
spectral normalization first. The continuum or spectral energy distribution
(SED) of input spectra derived from “Rassine” are denoted as
$\mathbf{SED}_{quiet}$ and $\mathbf{SED}_{active}$ respectively.

The inhibition of CB when working with PHOENIX spectra can be rewritten as:

$\Delta\mathbf{S}^{{}^{\prime}}_{bconv}(X,Y)=\mathbf{S}^{{}^{\prime}}_{quiet,n}(X,Y)-\mathbf{S}^{{}^{\prime}}_{active,n}(X,Y),$
(10)

where $\mathbf{S}^{{}^{\prime}}_{quiet,n}$ and
$\mathbf{S}^{{}^{\prime}}_{active,n}$ are the normalized quiet and active
spectra. Since the contrast between the quiet and active region is naturally
included in the continuum of input PHOENIX spectra, the flux effect can be
rewritten as:

$\displaystyle\Delta\mathbf{S}^{{}^{\prime}}_{flux}(X,Y)$
$\displaystyle=\mathbf{S}^{{}^{\prime}}_{quiet,n}(X,Y)\times\mathbf{SED}_{quiet}$
(11)
$\displaystyle-\mathbf{LB}(X,Y)\times\mathbf{S}^{{}^{\prime}}_{quiet,n}(X,Y)\times\mathbf{SED}_{active},$

where $\mathbf{LB}$ is the function of limb brightening. For simulation of
spot regions, $\mathbf{LB}=1.0$ along the disk. For simulation of faculae
regions, SOAP 2.0 was using $\Delta T_{f}=250.9-407.7\mu+190.9\mu^{2}$ to
model the limb brightening in temperature domain (Meunier et al., 2010). Since
this equation is only valid for the Sun, we use the empirical equation derived
from 3D MHD simulations to model other spectral types. 3D MHD simulations
modelling faculae on the Sun can reproduce extremely well the limb-brightening
observed for faculae (Norris et al., 2017). Using similar simulations, Johnson
et al. (2021a) model what would be the limb-brightening on other stars (see
Figure 3 and Table 1 in Johnson et al. (2021a)). Using the parametrisation in
Johnson et al. (2021a), we derived the limb brightening curves of faculae for
a G2, K0 and M0 dwarfs. We then linearly interpolated between the G2 and K0
and K0 and M0 simulations to obtain the limb-brightening dependence for a G8
and G9 dwarf, and a K2 dwarf, respectively. To obtain the dependence for a F9
dwarf, we linearly extrapolated from the G2 and K0 models. The derived limb
brightening curves from F9 to K2 are shown in Figure 11.

Finally, the total effect from flux and convection can be derived using the
following equation:

$\displaystyle\Delta\mathbf{S}^{{}^{\prime}}_{tot}(X,Y)$
$\displaystyle=\mathbf{S}^{{}^{\prime}}_{quiet,n}(X,Y)\times\mathbf{SED}_{quiet}$
(12)
$\displaystyle-\mathbf{LB}(X,Y)\times\mathbf{S}^{{}^{\prime}}_{active,n}(X,Y)\times\mathbf{SED}_{active}.$

Figure 11: Faculae intensity contrast as a function of $\mu$ for different
spectral types. The limb brightening curves for the G2, K0 and M0 dwarfs are
derived from the parametrisation of MHD simulations for 500G faculae, shown in
Table 1 in Johnson et al. (2021a). The limb brightening curves for other
spectral types are linearly interpolated by using the two closet limb
brightening curves. For spectral types hotter than G2, we obtain limb
brightening by performing linear extrapolation using the G2 and K0 curves. The
limb brightening derived from Meunier et al. (2010) is labeled with orange.

Convection velocity changes with spectral type, and decreases towards cooler
stars than the Sun. This is a well known effect, that comes out from
observations (e.g. Meunier et al., 2017; Liebing et al., 2021) and magneto
hydrodynamic models of stellar phostophere (e.g. Allende Prieto et al., 2013).
As in Liebing et al. (2021), we show in the top panel of Fig. 12 that the CB
is a cubic function of effective temperature in the range 4800 to 6300 $K$.
The parametrisation that we obtain is
$\mathrm{CB}_{vel}=95.2388\times((\mathrm{Teff}-4400)/1000)^{3}+91.2791$. Once
we have this relation to measure the velocity of CB as a function of effective
temperature, we simply have to rescale the difference between the quiet Sun
bisectors for $\mu=0.85$ and $\mu=0.2$ at depth 0.58 to be equal to the value
given by our relation, and we also impose that the active bisectors are
shifted by the same value, at the same depth (see Sect. 4.3.2 for
justification). This allows us to model properly the change of convective
velocity as a function of stellar effective temperature. We note that as in
the case of the Sun, before injecting the proper bisector for the quiet and
active regions, we first have to remove the bisector present in the PHOENIX
spectra that we use. This is a necessary step as the PHOENIX spectral library
is obtained from 1D atmospheric models and cannot properly reproduce line
bisector shape. We show in the appendix, like for Fig. 9 in the case of the
Sun, how the bisector of the original PHOENIX spectra for the quiet stellar
region, a faculae and a spot, are fitted before being removed.

Figure 12: CB velocity as a function of spectral type. _Top:_ Data from
Liebing et al. (2021) and cubic fit to them, giving as relation
$\mathrm{CB}_{vel}=95.2388\times((\mathrm{Teff}-4400)/1000)^{3}+91.2791$. We
show with coloured stars the CB velocity value for different spectral types
that we model in the paper. _Bottom:_ The bisector of the quiet photosphere
measured from Löhner-Böttcher et al. (2019) at $\mu=0.85$ (thin lines) and the
bisector of an active region measured from Cavallini et al. (1985a) at
$\mu=1.0$ (thick lines) for different spectral types. We impose that the
maximum difference between the quiet bisector at $\mu=0.85$ and $\mu=0.2$ (not
shown here), happening at a depth of 0.58 is equal to the CB velocity given by
the relation found in the top panel. We also force the difference between the
quiet bisector at $\mu=0.85$ and the active bisectors, independent of their
$\mu$ angle, to be equal to the same value at a depth of 0.58.

In Fig. 13, we show the results of a few simulations considering $\mu$
dependent input spectra and a single 1% equatorial spot (left panel) or a
single 1% facula (right panel). We highlight the RV outputs of SOAP-GPU when
using PHOENIX spectra, with the quiet Sun temperature set to
$\rm{T_{eff}}=5778K$, the spot temperature set to
$\rm{T_{eff}}=5115K,5015K\,\mathrm{and}\,5215K$ and the facula temperature set
to $\rm{T_{eff}}=5928K,6028K\,\mathrm{and}\,6128K$. We also show the result
when inputting the Kitt Peak solar quiet and spot or facula spectra
($\rm{T_{eff}}=5778K$ and $5115K$ or $6028K$ at disk center, respectively) as
modified in Sect. 4.3.2, and including (using Eq. 11, 10 and 12) the SED from
corresponding PHOENIX spectra. As we can see, the amplitude of the RV flux
effect increases with an increase in temperature difference between the quiet
photosphere and the spot or facula. The CB effect does not change in amplitude
with temperature difference and thus the larger amplitude observed for larger
difference in temperature between quiet and active regions is solely driven by
the change in contrast of the active regions with temperature. While for the
flux effect, the simulation using PHOENIX spectra or observed solar spectra as
input gives the same results, which is not surprising as we use the same SED,
this is not the case for the CB effect. The amplitude derived show a small
discrepancy of $\sim$20%, and the maximum has a slight phase shift. This
likely comes from a different flux effect contribution seen in the derived CB
RV, due to molecular absorption not perfectly modeled in the PHOENIX spectra
compared to solar real observations (see discussion in Sect. 4.3.1). We note
that this asymmetry was already something seen in the original SOAP 2.0 paper
(see Fig. 6 in Dumusque et al., 2014). Something also interesting to note,
that we see when using as input both the PHOENIX and solar spectra is the bump
in the CB RV effect when the active region crosses the center of the disk
(phase=0.5). This is induced by the fact that CB is maximum at $\mu=0.85$ and
not at disk center, as was shown in Löhner-Böttcher et al. (2019).

Figure 13: Simulation of an equatorial 1% spot (left) and 1% facula (right)
with different temperatures using PHOENIX spectral library. The quiet Sun
spectrum is extracted from PHOENIX spectral library with log(g) of 4.5 and
$\rm{T_{eff}}=5778K$. The effective temperature of the spot spectra are 5015K,
5115K and 5215K, and for the facula spectra 5928K, 6028K and 6128K. We also
show the results of using the observed Kitt Peak solar spectra including the
PHOENIX SEDs (using Eq. 11, 10 and 12). The $\mu$ dependent CB is included in
the simulations. _Left:_ For the spot we see that when the difference in
temperature between the spot and the quiet photosphere increases, the RV flux
effect becomes larger. We note a rather large discrepancy in the CB effect
between the observed solar and PHOENIX input spectra. This is likely due to
the fact that the observed solar spectra are not well normalized (see Sect
4.3.3). _Right_ : For the facula, we note the same problem of negative values
for the CB effect, which expected as the same badly normalised solar active
spectrum is used. We also see that considering the corresponding PHOENIX SED
when using the observed solar spectra significantly change the flux
contribution (see Sect 4.3.3). Considering the PHOENIX SED gives results much
closer to the PHOENIX simulations.

In Fig.14 we show the result of the estimated CB RV effect for an equatorial
1% spot or facula for stars of different temperature (i.e. spectral type):
6050 $K$ (F9), 5778 $K$ (G2), 5480 $K$ (G8), 5380 $K$ (G9) and 5100 $K$ (K2).
For the spot, we see a positive only effect for the F9 and G2 simulations,
which is expected, but for later spectral type, we start to see the emergence
of a flux effect. This effect, as already discussed in Sect.4.3.1 comes from
the absorption of molecules, that change significantly over a few hundreds of
Kelvin for the quiet photosphere at 5480 $K$ and a spot at 4817 $K$ for the G8
simulation for example. Also, we see that the more we go towards cooler stars,
the more the CB effect show a flux contribution, and this comes from the fact
that molecular absorption is not linear with effective temperature. Although
molecular absorption prevent us of clearly separating the flux effect of spots
from the CB effect like in the case of the F9 or G2 star, the total RV effect
including both contributions is still properly estimated. Therefore, users
should be careful when interpreting the estimated RV induced by the inhibition
of convection for spots on stars with spectral type later than the Sun,
however, they can trust the total RV effect estimated. Regarding the facula,
we observe a positive only effect for all spectral type, and thus the CB
effect is properly modeled for this type of active regions.

Figure 14: Simulation of the CB RV effect of an equatorial 1% spot (top) and
1% facula (bottom) for different temperature of the quiet photosphere (i.e.
different spectral type) using the PHOENIX spectral library. The temperature
difference between the quiet photosphere $\Delta\rm{T_{eff}}$ for spot and
facula is fixed to 663K and 250 $K$, respectively. Due to the strong
absorption effect of molecules in spots on stars cooler than the Sun, we see
the appearance of a flux effect in the derivation for the CB effect only (see
Sect 4.3.3).

####  4.3.4 Limitation of SOAP-GPU when moving away from solar twins

On a careful note, we want to warn the user that a lot of the physics included
in SOAP-GPU is based on solar observations, like for example the variation of
the bisector of the quiet photosphere with respect to the center-to-limb angle
(Löhner-Böttcher et al., 2019), or the bisector of an active region extracted
from solar observations (Cavallini et al., 1985a). Therefore, although we try
to correct for some effects, like the variation of CB velocity as a function
of spectral type, the more we go away from the solar case, the more we should
be careful about interpreting the results coming out of SOAP-GPU.

We note that when modifying the bisector shape of solar or PHEONIX spectra for
the disk center, to account for different convection velocity across spectral
type, we model the effect of the “third signature” of granulation (Gray,
2009). The inherent shape of the bisector (known as the “second signature” of
granulation) and how it varies with $\mu$ (Löhner-Böttcher et al., 2019) is
still the one observed for the Sun and it is well known in the literature that
the bisector shape of disc-integrated observations, and therefore by analogy
at disk center, varies significantly among luminosity class and spectral type
(see Figure 17.15 in Gray, 2008; Baştürk et al., 2011). In those papers, we
see that inside the small range from early G’s dwarfs to early K’s dwarfs, the
bisector shapes does not change drastically and therefore, although the
integrated spectra for those stars won’t be realist in terms of line bisector,
the way SOAP-GPU models stellar activity should still be quite realistic as in
the end, what counts is the differential between the activity and quiet
phases. As we can see in the preceding subsection, besides SOAP-GPU not being
able anymore to separate the inhibition of the CB effect from the flux effect
for early K’s, the tests we performed show that the code seems to behaves
quite well in estimating the total RV effect induced by spots and faculae.
However, outside of this small range from G to K dwarfs, bisectors are
completely different, and we warn the users that the present version of SOAP-
GPU might give very unrealistic stellar integrated spectra and estimation of
stellar activity. There is perhaps only one exception for dwarfs cooler than
early K’s, for which the velocity of convection decreases to level that are
very difficult to measure. For those stars, stellar activity is dominated by
the flux effect from spot and faculae. Therefore, for such stars, users should
ignore the output from the CB effect and only consider the flux effect.

To better model stellar activity for stars different than the Sun, we would
require disk-resolved spectra for other stars. Those could come from 3D MHD
simulations (e.g. Johnson et al., 2021b; Dravins et al., 2021), or thanks to
spatially resolved spectroscopic observation across stellar surfaces thanks to
transiting planets (e.g. Dravins et al., 2017). Both approach are challenging,
however could led to a much more realistic modelisation of stellar activity
from other stars than our Sun. We note that it is rather easy to input
different spectra in SOAP-GPU than the ones provided (the Sun and PHEONIX
spectra), and therefore if users have access to disk-resolved spectra with
more realistic line-bisectors, SOAP-GPU could still be used to obtain
efficiently disk-integrated spectra and model the corresponding stellar
activity effect. We note that the CB injection part is an individual module in
SOAP-GPU that users can easily modify.

##  5 Conclusions

In this paper, we present a GPU-based improvement to SOAP 2.0, named SOAP-GPU,
that allows to efficiently model stellar activity at the spectral level. With
the implementation of GPU interpolation and summation, benchmark calculations
demonstrate that SOAP-GPU improves the computational speed by a factor of 60
when modeling stellar activity on a full visible spectral range at R=115’000
of resolution, while having the same accuracy.

Beside the huge gain in speed, SOAP-GPU also provides a more complex
modelisation of stellar activity compared to SOAP 2.0. Complex active region
scenarios, with regions overlapping is now handled by the code. This is mainly
useful when modeling active phases of a star like the Sun, with hundreds of
active regions. The contrast of the active regions is now wavelength
dependant, and therefore change for each wavelength of the modeled spectra.
The dependence of line bisector with center-to-limb angle $\mu$, following the
work of Löhner-Böttcher et al. (2019) and Cavallini et al. (1985a), is also
now accounted for for the Kitt Peak observed quiet and active atlases, but
also for PHOENIX spectra. Although the induced RV effect is rather negligible
when using the observed solar spectra as input, including the framework to
change line bisector is crucial to properly model convection when injecting
PHOENIX spectra as input. The use of PHOENIX spectra allows us now to model a
wide variety of stars with different stellar and active region properties, and
allows us as well to better model faculae, as the corresponding spectrum now
has the proper effective temperature (SOAP 2.0 was using the Kitt Peak sunspot
atlas to model faculae).

When modeling the inhibition of CB effect using as input the solar Kitt Peak
quiet and active spectral atlases, we noticed that the derived RVs go
negative, which is not expected. This comes from molecular absorption that can
be seen in the spot spectrum due to lower temperature compared to the quiet
Sun. Even though we do not include the contrast of the active region when
modeling the RV CB effect only (see Eqs. 2 and 10), the difference in flux at
the level of molecular absorption bands will show up as a flux effect in the
estimated CB RVs. Positive values will be added to the CB RV effect before the
spot crosses the stellar center, and negative values after, therefore creating
an asymmetry.

When modeling other stars than the Sun using PHOENIX spectral library, users
should be aware that a lot of physics included in SOAP-GPU are based on solar
observations, and although the code tries to correct for known effects like
the variation of CB velocity as a function of effective temperature (i.e
spectral type, see Sect. 4.3.3), the more we go away from the Sun, the more
the results should be interpreted with caution (see discussion in Sect.
4.3.4). The modelisation of stellar activity for other stars than the Sun is
currently limited by the knowledge we have about disk-resolved bisectors for
such stars. Such information is very challenging to obtain, however, 3D MHD
simulations (e.g. Johnson et al., 2021b; Dravins et al., 2021) and resolved
spectroscopic observation of other stars due to planetary transits (e.g.
Dravins et al., 2017) could significantly help.

Also, when modeling stars of later spectral type than the Sun, we are not able
anymore to separate clearly the inhibition of CB effect from the flux effect
due to the strong absorption of molecules. However, the output for the total
RV effect (flux plus inhibition of CB effects) should be modeled properly.
SOAP-GPU have been tested up to a K2 star (T${}_{\mathrm{eff}}=5100\,K$) with
spots 663 $K$ cooler and give satisfactory results. Modeling later spectral
type is challenging mainly due to continuum normalisation of the PHOENIX
spectra and injection of spectral line bisector due to line blending, and
users should be very careful about the interpretation of the results for such
stars with the present code.

There are still some improvements that could be made to better model the
physics at play. Although the spot and facula spectrum used when considering
PHOENIX spectra as input are of different temperature, and therefore in
spectral content, we still associate to those regions the same active bisector
as measured for the Sun on a facula (Cavallini et al., 1985a). Spots are
induced by stronger magnetic fields than facula, and thus it is likely that
the bisector of spectral lines will be slightly different. Although it is
possible to know what is the bisector of a few spectral lines inside a spot at
disk center, to our knowledge, no measurement of spot line bisectors for
different $\mu$ angles are published. Cavallini et al. (1985a) also show in
their Fig. 1 that depending on the facula observed, the bisector shape changes
due likely to different magnetic field strength and therefore different level
of CB inhibition. As can be seen in Fig. 10, the effect of inducing $\mu$
dependant spectral line shape in the quiet and active regions is rather small,
and although with more solar data about spots and faculae we could better
model the physics at play, results in terms of RV derivation would be rather
similar. This likely comes from the fact due to limb-darkening, most of the
weight is put on the disk center, where spectral lines does not change
significantly in shape.

With the performance of SOAP-GPU, it is now possible to model activity at the
spectral level for complex stellar surfaces with many active regions and for a
long period of time. A solar activity simulator, either based on statistical
properties of solar active regions (similar to Borgniet et al., 2015) or on
the observed distribution of those (similar to e.g. Meunier et al., 2010) will
be published in a forthcoming paper. We encourage any person working on
techniques to separate the activity effect from planetary signals at the
spectral level, to test their framework on SOAP-GPU simulations, where photon-
noise, instrumental and telluric systematics are not perturbing the spectral
timeseries.

###### Acknowledgements.

We thank the anonymous referee for the insightful and constructive comments on
this paper. We thank Michael Crerignier for his help in normalizing PHOENIX
spectra with RASSINE. We also thank Xiang Gao for the constructive comments on
GPU computing. This project has received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement SCORE No 851555). This work has been carried out
within the framework of the National Centre of Competence in Research PlanetS
supported by the Swiss National Science Foundation. The authors acknowledge
the financial support of the SNSF.

## References

  * Aigrain et al. (2016) Aigrain, S., Parviainen, H., & Pope, B. J. S. 2016, MNRAS, 459, 2408
  * Aigrain et al. (2012) Aigrain, S., Pont, F., & Zucker, S. 2012, MNRAS, 419, 3147
  * Allende Prieto et al. (2013) Allende Prieto, C., Koesterke, L., Ludwig, H.-G., Freytag, B., & Caffau, E. 2013, A&A, 550, A103
  * Baştürk et al. (2011) Baştürk, Ö., Dall, T. H., Collet, R., Lo Curto, G., & Selam, S. O. 2011, A&A, 535, A17
  * Baranne et al. (1996) Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373
  * Barragán et al. (2022) Barragán, O., Aigrain, S., Rajpaul, V. M., & Zicher, N. 2022, MNRAS, 509, 866
  * Binnenfeld et al. (2022) Binnenfeld, A., Shahaf, S., Anderson, R. I., & Zucker, S. 2022, A&A, 659, A189
  * Binnenfeld et al. (2020) Binnenfeld, A., Shahaf, S., & Zucker, S. 2020, A&A, 642, A146
  * Boisse et al. (2012) Boisse, I., Bonfils, X., & Santos, N. C. 2012, A&A, 545, 109
  * Borgniet et al. (2015) Borgniet, S., Meunier, N., & Lagrange, A. M. 2015, A&A, 581, A133
  * Brewer et al. (2020) Brewer, J. M., Fischer, D. A., Blackman, R. T., et al. 2020, AJ, 160, 67
  * Cavallini et al. (1985a) Cavallini, F., Ceppatelli, G., & Righini, A. 1985a, A&A, 143, 116
  * Cavallini et al. (1985b) Cavallini, F., Ceppatelli, G., & Righini, A. 1985b, A&A, 150, 256
  * Chapman et al. (2001) Chapman, G. A., Cookson, A. M., Dobias, J. J., & Walton, S. R. 2001, ApJ, 555, 462
  * Collier Cameron et al. (2021) Collier Cameron, A., Ford, E. B., Shahaf, S., et al. 2021, MNRAS, 505, 1699
  * Collier Cameron et al. (2019) Collier Cameron, A., Mortier, A., Phillips, D., et al. 2019, MNRAS, 487, 1082
  * Cretignier et al. (2020a) Cretignier, M., Dumusque, X., Allart, R., Pepe, F., & Lovis, C. 2020a, A&A, 633, A76
  * Cretignier et al. (2022) Cretignier, M., Dumusque, X., & Pepe, F. 2022, A&A, 659, A68
  * Cretignier et al. (2020b) Cretignier, M., Francfort, J., Dumusque, X., Allart, R., & Pepe, F. 2020b, A&A, 640, A42
  * Davis et al. (2017) Davis, A. B., Cisewski, J., Dumusque, X., Fischer, D. A., & Ford, E. B. 2017, ApJ, 846, 59
  * de Beurs et al. (2022) de Beurs, Z. L., Vanderburg, A., Shallue, C. J., et al. 2022, AJ, 164, 49
  * Donati et al. (2017) Donati, J. F., Yu, L., Moutou, C., et al. 2017, MNRAS, 465, 3343
  * Dravins et al. (1981) Dravins, D., Lindegren, L., & Nordlund, A. 1981, A&A, 96, 345
  * Dravins et al. (2017) Dravins, D., Ludwig, H.-G., Dahlén, E., & Pazira, H. 2017, A&A, 605, A91
  * Dravins et al. (2021) Dravins, D., Ludwig, H.-G., & Freytag, B. 2021, A&A, 649, A16
  * Dumusque (2016) Dumusque, X. 2016, A&A, 593, A5
  * Dumusque et al. (2014) Dumusque, X., Boisse, I., & Santos, N. C. 2014, ApJ, 796, 132
  * Dumusque et al. (2021) Dumusque, X., Cretignier, M., Sosnowska, D., et al. 2021, A&A, 648, A103
  * Dumusque et al. (2015) Dumusque, X., Glenday, A., Phillips, D. F., et al. 2015, ApJ, 814, L21
  * Gilbertson et al. (2020a) Gilbertson, C., Ford, E. B., & Dumusque, X. 2020a, Research Notes of the American Astronomical Society, 4, 59
  * Gilbertson et al. (2020b) Gilbertson, C., Ford, E. B., Jones, D. E., & Stenning, D. C. 2020b, ApJ, 905, 155
  * Gray (2008) Gray, D. F. 2008, The Observation and Analysis of Stellar Photospheres, ed. U. C. U. P. Cambridge
  * Gray (2009) Gray, D. F. 2009, ApJ, 697, 1032
  * Herrero et al. (2016) Herrero, E., Ribas, I., Jordi, C., et al. 2016, A&A, 586, A131
  * Husser et al. (2013) Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6
  * Johnson et al. (2021a) Johnson, L. J., Norris, C. M., Unruh, Y. C., et al. 2021a, MNRAS, 504, 4751
  * Johnson et al. (2021b) Johnson, L. J., Norris, C. M., Unruh, Y. C., et al. 2021b, MNRAS, 504, 4751
  * Jones et al. (2017) Jones, D. E., Stenning, D. C., Ford, E. B., et al. 2017, arXiv e-prints, arXiv:1711.01318
  * Liebing et al. (2021) Liebing, F., Jeffers, S. V., Reiners, A., & Zechmeister, M. 2021, A&A, 654, A168
  * Lin et al. (2022) Lin, A. S. J., Monson, A., Mahadevan, S., et al. 2022, AJ, 163, 184
  * Löhner-Böttcher et al. (2019) Löhner-Böttcher, J., Schmidt, W., Schlichenmaier, R., Steinmetz, T., & Holzwarth, R. 2019, A&A, 624, A57
  * Löhner-Böttcher et al. (2018) Löhner-Böttcher, J., Schmidt, W., Stief, F., Steinmetz, T., & Holzwarth, R. 2018, A&A, 611, A4
  * Ludwig et al. (2009) Ludwig, H. G., Caffau, E., Steffen, M., et al. 2009, Mem. Soc. Astron. Italiana, 80, 711
  * Mayor & Queloz (1995) Mayor, M. & Queloz, D. 1995, Nature, 378, 355
  * Meunier et al. (2010) Meunier, N., Desort, M., & Lagrange, A.-M. 2010, A&A, 512, A39
  * Meunier et al. (2017) Meunier, N., Lagrange, A.-M., Mbemba Kabuiku, L., et al. 2017, A&A, 597, A52
  * Milbourne et al. (2019) Milbourne, T. W., Haywood, R. D., Phillips, D. F., et al. 2019, ApJ, 874, 107
  * Muraközy et al. (2014) Muraközy, J., Baranyi, T., & Ludmány, A. 2014, Sol. Phys., 289, 563
  * Norris et al. (2017) Norris, C. M., Beeck, B., Unruh, Y. C., et al. 2017, A&A, 605, A45
  * Pepe et al. (2021) Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96
  * Rajpaul et al. (2015) Rajpaul, V., Aigrain, S., Osborne, M. A., Reece, S., & Roberts, S. 2015, MNRAS, 452, 2269
  * Rajpaul et al. (2020) Rajpaul, V. M., Aigrain, S., & Buchhave, L. A. 2020, MNRAS, 492, 3960
  * Saar & Donahue (1997) Saar, S. H. & Donahue, R. A. 1997, ApJ, 485, 319
  * Stief et al. (2019) Stief, F., Löhner-Böttcher, J., Schmidt, W., Steinmetz, T., & Holzwarth, R. 2019, A&A, 622, A34
  * Wallace et al. (1998) Wallace, L., Hinkle, K., & Livingston, W. 1998, An atlas of the spectrum of the solar photosphere from 13,500 to 28,000 cm-1 (3570 to 7405 A)
  * Wallace et al. (2005) Wallace, L., Hinkle, K., & Livingston, W. C. 2005, An atlas of sunspot umbral spectra in the visible from 15,000 to 25,500 cm-1 (3920 to 6664 Å)
  * Zhao et al. (2022a) Zhao, J., Ford, E. B., & Tinney, C. G. 2022a, ApJ, 935, 75
  * Zhao et al. (2022b) Zhao, L. L., Fischer, D. A., Ford, E. B., et al. 2022b, AJ, 163, 171

##  Appendix A Line bisectors of PHOENIX spectra

As discussed in Sects. 4.3.2 and 4.3.3, before injecting the $\mu$ dependant
bisector for solar or PHOENIX spectra to properly model CB and its inhibition
close to the limb and in active regions, we need to remove any bisector shape
already present in the input spectra. As PHOENIX spectral library is generated
from 1D spectral synthesis, the line bisectors cannot include properly the CB
effect and therefore should be close to straight. In Fig. 15, we show for each
simulation of different spectral types the bisector of a few iron lines that
are used in Löhner-Böttcher et al. (2019). For each spectral type simulated,
we show the bisectors for the quiet photosphere, but also for simulated spot
and faculae, 663 $K$ cooler or 250 $K$ hotter, respectively. As expected, most
of the bisector are close to straight lines. We however fitted the average
bisector with a second order polynomial to remove the small curvatures
observed before injecting the proper bisectors at different $\mu$ angles (see
Sect. 4.3.2). It is not clear if those curvatures are real effect in the
spectral synthesis, or simply due to blends. The correction performed is small
compared to the bisectors that we inject afterward, therefore if only due to
blends, this process does not significantly change the outputs.

Figure 15: Bisector of PHOENIX spectra. For each input seed spectrum using
PHOENIX spectral library, we use five strong iron lines: FeI $5250.2084\AA$
(green), FeI $5250.6453\AA$ (cyan), FeI $5434.5232\AA$ (purple), FeI
$6173.3344\AA$ (orange) and FeI $6301.5008\AA$ (yellow) to measure the average
bisector of the input spectra. Bisector outliers outside a window of $0.1\AA$
around each line center are rejected to avoid those points, certainly affected
by line blending, to bias our measurement of line bisector. Each line
correspond to a different spectral type, and from left to right, we can see
the bisector of the spectrum used for the quiet photosphere, a spot region
(663 $K$ cooler) and a facula region (250 $K$ hotter). We average those line
bisectors at certain depth (as shown by the red dots) and fit the obtained
data with a second order polynomial. The fitted bisector is used to remove the
bisector of input seed spectrum.