© The Authors 2023

1 Introduction

The radial velocity (RV) method has been proven to be one of the most successful methods to date in detecting exoplanets since the discovery of the first exoplanet orbiting a solar-type star (Mayor & Queloz 1995). To detect Earth-like planets orbiting in the habitable zone of its parent star, a precision of a few dozens of cm s⁻¹ must be reached. Although the state-of-the-art spectrographs such as ESPRESSO and EXPRESS are not far from that precision (50 and 58 cm s⁻¹, respectively Pepe et al. 2021; Brewer et al. 2020), the main limitation in detecting Earth-like planets with the RV technique is issues related to stellar activity. Two major physical processes dominating stellar activity on a timescale of the host star’s rotational period are the flux imbalance due to the temperature difference – and the ensuing contrast between active and quiet regions, hereafter referred to as the flux effect. (e.g., Saar & Donahue 1997; Dumusque et al. 2014; Donati et al. 2017) – and the inhibition of convective blueshift, hereafter referred to as the 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 (1D) 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, 2016; Rajpaul et al. 2015; 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 when modeling stellar activity in 1D time series, the community is now shifting toward modeling it in a two-dimensional (2D) 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 changes in shape that are related to stellar activity. A planet signal of amplitude ~40 cm s⁻¹ can be recovered when the algorithm is applied to the HARPS-N solar data (Dumusque et al. 2015, 2021; Collier Cameron et al. 2019). Zhao et al. (2022a) projected CCF time series onto the Fourier basis functions and modeled line variability using different basis. Results on simulated data show a 48% reduction in the 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 HD 10700 (τ Ceti) and HD12861 (α 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 HARPS¹, and (more recently) the solar feed of NEID (Lin et al. 2022) are the best we can get, in terms of the signal-to-noise ratio (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 stars that are 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) that is based on GPU computation that able to efficiently model simplified and realistic stellar activity at the spectral level. In Sect. 2, we revisit the architecture of the SOAP 2.0 code that 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 Zenodo², along with a brief manual and some examples.

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

2 Revisiting SOAP 2.0

In this section, we revisit the Spot Oscillation And Planet (SOAP 2.0, Dumusque et al. 2014) code, which is aimed 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.

2.1 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 2D stellar disk containing N × N cells is initialized (N being the resolution of the disk, with the same parameter referred to as the “grid” in Boisse et al. 2012). The 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(λ), Wallace et al. 1998) as an 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 ~500000, compared to 400 for the CCF – we would need to apply a Doppler-shift each time we wanted to change the velocity of this spectrum, whereas 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).

The next step consists of 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(λ) with the value of λ the same as it is for S_quiet(λ) 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, it is 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, while also including the effect of CB using results from magnetohydrodynamical (MHD) simulations (e.g., with the STARSIM 2 code³ Herrero et al. 2016). We note that injecting observed or synthetic spectra have their advantages and drawbacks. Using observed spectra allows us 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 us to better model the temperature (and therefore the 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 Eqs. (1) and (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 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: $$\text{Δ}{S^{\prime }}_{\text{flux}}\left(X,Y\right)={S^{\prime }}_{\text{quiet}}\left(X,Y\right)-I_{\text{ratio}}\times {S^{\prime }}_{\text{quiet}}\left(X,Y\right),$$(1) $$\text{Δ}{S^{\prime }}_{\text{bconv}}\left(X,Y\right)={S^{\prime }}_{\text{quiet}}\left(X,Y\right)-{S^{\prime }}_{\text{active}}\left(X,Y\right),$$(2) $$\text{Δ}{S^{\prime }}_{\text{tot}}\left(X,Y\right)={S^{\prime }}_{\text{quiet}}\left(X,Y\right)-I_{\text{ratio}}\times {S^{\prime }}_{\text{active}}\left(X,Y\right),$$(3)

where 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: $$S_{\text{integrated,final}}=S_{\text{integrated,quiet}}-\text{Δ}{S}_{\text{integrated,quiet}-\text{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.

2.2 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 ~500 000 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 ~500 000. For example, SOAP 2.0 takes ~800 s to calculate an integrated quiet sun spectrum using a 300 × 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 demonstrate the limitations 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, which is much more efficient in term of computational speed, but also that adds some physical complexity.

3.1 Basic concept of GPU computing

The popularity of artificial intelligence has, in recent years, significantly increased due to the programmability of graphic hardwares. The GPU computing uses graphical card as a coprocessor 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), which is a compiler and toolkit for programming NVIDIA GPUs and is an extension of the C/C++ programming language. It invokes kernel functions by using the syntax of <<< N_blocks, M_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 a 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 follows: $${{\lambda }^{\prime }}_n={\lambda }_n+{\lambda }_{n}f\left(\beta \right),$$(6)

where $f\left(\beta \right)=-\left[1-\sqrt{\frac{\left(1+\beta \right)}{\left(1-\beta \right)}}\right]$ and β = υ/c. The variable υ 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(λ′) → S′(λ), we have to find the index m which satisfies λ′_n < λ_m < λ′_n+1. For the left side, we have: $$\begin{array}{l}{{\lambda }^{\prime }}_n<{\lambda }_m,\hfill \\ {{\lambda }^{\prime }}_n={\lambda }_n+{\lambda }_{n}f\left(\beta \right)={\lambda }_0+nk+{\lambda }_{0}f\left(\beta \right)+nkf\left(\beta \right)<{\lambda }_0+mk,\hfill \end{array}$$

so we have: $$n\left(1+f\left(\beta \right)\right)+\frac{f\left(\beta \right){\lambda }_0}{k}<m.$$(7)

For the right side, we have: $$m<\left(n+1\right)\left(1+f\left(\beta \right)\right)+\frac{f\left(\beta \right){\lambda }_0}{k}.$$(8)

Once the integer m is known, we can estimate the flux for λ_m using the spectrum derivative: $${S^{\prime }}_m=\frac{\text{Δ}{S}_n}{\left({{\lambda }^{\prime }}_{n+1}-{{\lambda }^{\prime }}_n\right)}\times \left({\lambda }_m-{{\lambda }^{\prime }}_n\right)+S_n,$$(9)

where S_n = S(λ_n) and ∆S_n = S_n+1 − S_n.

Equations (7)–(9) can be parallelized 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 × 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.

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 active 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 summarized in Algorithm 5.

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 ω = ω₀ + ω₁ sin²(θ), where ω₀ = 14.371°/day and ω₁ = −2.587°/day for the Sun (Borgniet et al. 2015). To generalize this for other stars, the user can select in the configuration of SOAP-GPU a rotation period and a differential rotation rate. ω0 is then equal to 360/PROT and ω₁ to DIFF_ROT*PROT (PROT = 25.05 and DIFF_ROT = −0.18 for the solar case to reproduce the above equation).

Fig. 2 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 Results

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 analyzed 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 547 840, 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 Fig. 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 disk center with an area of 1% of the entire disk will be contained in a grid of 34 × 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 ~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 modeling thus makes significant improvements in computational speed when considering high-resolution simulations or simulations with large active regions.

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 s to simulate those spectra with SOAP 2.0 for 100 timestamps, whereas it took only 27.9 s 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 fit 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.

Fig. 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 s to simulate those spectra with SOAP 2.0, while only 27.9 s with SOAP-GPU on a Nvidia RTX-3090 card. The computation speed is improved by a factor of 63.

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 Size evolution of active regions

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 10 000 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.

4.2.2 Complex active regions

SOAP-GPU also allows users to simulate complex active region configurations. From the observational point of view, faculae 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 Algorithm5. A simulation of this case is illustrated in Fig. 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 was already demonstrated in several studies (e.g., Meunier et al. 2010; Dumusque et al. 2014; Milbourne et al. 2019).

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 >>>.

The measured RVs are shown in Fig. 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 the same spectrum for the spotted region as we did for 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 show in Sect. 4.3.3, we did 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 contribution from the flux effect. 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, since, after computing the CCF for the quiet and actives regions, we were renormalizing 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 temperature⁴, at the average wavelength of the input spectra 5293 Å. 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 SOAPcontrast <<< N_blocks, N_threads >>>. This new kernel allows us 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.

Fig. 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.

Fig. 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 presented 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.

Fig. 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 547 840 is labeled in black. The blue and red parts of the spectra with length 204 800 are over plotted in blue and red, respectively.

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 diskcenter 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, 2019; Stief et al. 2019).

We note that the varying shape of spectral line with center-to-limb μ angle is also modeled 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 fourth-order polynomial is only valid for line depth as strong as ~0.5 (see Fig. 6 in Herrero et al. 2016)⁵. 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 μ angles provided in Löhner-Böttcher et al. (2019). We first measured the bisectors of all the available iron deep lines at different μ angle in Löhner-Böttcher et al. (2019)⁶ and fit them using polynomial functions. For μ 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 μ angle larger or equal to 0.5, third-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–1.

To obtain the dependency of the spectral line bisector as a function of μ in an active region, we use the observations presented in Cavallini et al. (1985a). We parameterized the bisectors of the FeI at 6301.5008 Å, for the disk center (μ = 1) and different center-to-limb angles (μ = 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 representing the available data, using a fifth-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 the obtained active bisectors from depth 0.0–1.0 in Fig. 8.

Once we have our model for line bisectors at different μ angles, we can use the Python module Convec.model to apply those bisectors to the original spectra, and thus obtain different spectra for different μ angles. Each cell in the stellar disk takes the bisector that has the closest μ 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 selected the same lines as in Löhner-Böttcher et al. (2019) in the input spectra and measured 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 fitted lines are blended, which strongly impacts the bisector measurement. Thus, we rejected bisector points that were significantly off. Then, we modeled 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 μ angles, we removed the original bisectors of the quiet and active Kitt peak solar spectra and then added the bisectors measured for different μ angles. This was done by shifting each point in those spectra depending on their normalized depth.

To inject the proper bisectors for different μ angles in our original spectra, we first removed the original bisector, which changes any CB difference between the quiet and active solar spectra. Thus, we 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. At this point, we make two assumptions: (i) the CB is fully inhibited at μ = 0.2 in the quiet Sun and (ii) it is also fully inhibited for magnetic regions, and this at all μ angles. Using the first assumption, for the quiet Sun we measured the maximum shift between the bisector at μ = 0.85, which is the bisector that is the most blueshifted, and the bisector at μ = 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 rescaled the maximum difference between the μ = 0.85 and μ = 0.2 to be 340 m s⁻¹. We note that this rescaling is negligible in the case of the Sun, however, it will be needed when using PHOENIX spectra as input (as detailed in Sect. 4.3.3). Using the second assumption, we impose that at the same depth of 0.58, the difference in velocity between the quiet bisector at μ = 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 μ 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 μ angles (close to the disk 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 Sect. 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 will be 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 modeled 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 μ are derived from the quiet photosphere (Löhner-Böttcher et al. 2019) and faculae (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.

Fig. 7 Chromatic effect of RVs. SOAP-GPU is initialized with three different spectra: entire wavelength coverage and only red and only 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: measured RVs of the total effect. Middle: measured RVs of the flux effect. An offset of 1 m s−1 is added in the red and blue RVs. Bottom: measured RVs of the CB effect.

Fig. 8 Average bisectors of quiet and active solar regions from the disk center (μ = 1.0) to the limb (μ = 0.2). Continuous lines: fifth-order polynomial fit to the quiet sun bisectors of the FeI 5250.2084, FeI 5250.6453, FeI 5434.5232, FeI 5432.9470, FeI 5576.0881, FeI 6149.2460, FeI 6173.3344, FeI 6301.5008, and FeI 6302.4932 Å 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 Å 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−1 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 μ angles are all shifted by those 340 m s−1 at a depth of 0.58 based on our hypothesis that convection is fully suppressed in magnetic regions (see Sect. 4.3.2 for more information).

Fig. 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 Å bisector because it was significantly off by 2500 m s−1 due to strong contamination by other weak lines.

Fig. 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 μ dependent bisector injection. Dotted line: RVs derived from simulated spectra with μ dependent bisector injection. The input spectra with μ-angle dependency are generated with the Python module Convec.model. The RVs of the flux, CB, and total effect do not change significantly when the μ-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 Å to 5.5 μm with resolutions of 500 000 in the optical. The library covers stellar effective temperature from 2300 to 12000 K. Since the spectra in the PHOENIX library are not normalized, which is critical for performing the injection of CB described in the preceding section, we used the open source package “Rassine” (Cretignier et al. 2020b) to perform the spectral normalization first. The continuum or spectral energy distribution (SED) of input spectra derived from “Rassine” are denoted as SED_quiet and SED_active, respectively.

The inhibition of CB when working with PHOENIX spectra can be rewritten as: $$\text{Δ}{S^{\prime }}_{\text{bconv}}\left(X,Y\right)={S^{\prime }}_{\text{quiet},n}\left(X,Y\right)-{S^{\prime }}_{\text{active},n}\left(X,Y\right),$$(10)

where S′_quiet,n and S′_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: $$\begin{array}{ll}\text{Δ}{S^{\prime }}_{\text{flux}}\left(X,Y\right)=\hfill & {S^{\prime }}_{\text{quiet},n}\left(X,Y\right)\times SE{D}_{\text{quiet}}\hfill \\ \hfill & -LB\left(X,Y\right)\times {S^{\prime }}_{\text{quiet},n}\left(X,Y\right)\times SE{D}_{\text{active}},\hfill \end{array}$$(11)

where LB is the function of limb brightening. For the simulation of spot regions, LB = 1.0 along the disk. For simulation of faculae regions, SOAP 2.0 was using ∆T_f = 250.9−407.7μ + 190.9μ² 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. The 3D MHD simulations used to model faculae on the Sun are able to reproduce the limb-brightening observed for faculae extremely well (Norris et al. 2017). Using similar simulations, Johnson et al. (2021) modeled what would be the limb-brightening on other stars (see Fig. 3 and Table 1 in Johnson et al. 2021). Using the parametrization in Johnson et al. (2021), we derived the limb-brightening curves of faculae for a G2, K0, and M0 dwarfs. We then linearly interpolated between the G2 and K0 to obtain the limb-brightening dependence for a G8 and G9 dwarf and linearly interpolated between the K0 and M0 simulations to obtain the limb-brightening dependence for a K2 dwarf. 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 Fig. 11.

Finally, the total effect from flux and convection can be derived using the following equation: $$\begin{array}{ll}\text{Δ}{S^{\prime }}_{\text{tot}}\left(X,Y\right)=\hfill & {S^{\prime }}_{\text{quiet},n}\left(X,Y\right)\times SE{D}_{\text{quiet}}\hfill \\ \hfill & -LB\left(X,Y\right)\times {S^{\prime }}_{\text{active},n}\left(X,Y\right)\times SE{D}_{\text{active}}.\hfill \end{array}$$(12)

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 the 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–6300 K. The parametrization that we obtain is CB_vel = 95.2388 × ((T_eff − 4400)/1000)³ + 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 μ = 0.85 and μ = 0.2 at depth 0.58 to be equal to the value given by our relation; in addition, we impose the condition 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 properly model 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 used. 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, as we do in Fig. 9 for 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.

In Fig. 13, we show the results of a few simulations considering μ 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 T_eff = 5778 K, the spot temperature set to T_eff = 5115, 5015 and 5215 K and the facula temperature set to T_eff = 5928, 6028, and 6128 K. We also show the result when inputting the Kitt Peak solar quiet and spot or facula spectra (T_eff = 5778 and 5115 K or 6028 K at disk center, respectively) as modified in Sect. 4.3.2 and including (using Eqs. (10)–(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 to 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 ~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 being perfectly modeled in the PHOENIX spectra compared to solar real observations (see the 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 μ = 0.85 and not at disk center, as was shown in Löhner-Böttcher et al. (2019).

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 (F9), 5778 (G2), 5480 (G8), 5380 (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 previously discussed in Sect. 4.3.1) comes from the absorption of molecules, which can 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.

Fig.11 Faculae intensity contrast as a function of μ for different spectral types. The limb brightening curves for the G2, K0 and M0 dwarfs are derived from the parametrization of MHD simulations for 500G faculae, shown in Table 1 in Johnson et al. (2021). 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 obtained the 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.

Fig. 12 CB velocity as a function of spectral type. Top: data from Liebing et al. (2021) and cubic fit to them, giving as relation CBvel = 95.2388 × ((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: bisector of the quiet photosphere measured from Löhner-Böttcher et al. (2019) at μ = 0.85 (thin lines) and the bisector of an active region measured from Cavallini et al. (1985a) at μ = 1.0 (thick lines) for different spectral types. We impose that the maximum difference between the quiet bisector at μ = 0.85 and μ = 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 μ = 0.85 and the active bisectors, independent of their μ angle, to be equal to the same value at a depth of 0.58.

Fig. 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 Teff = 5778 K. The effective temperature of the spot spectra are 5015, 5115, and 5215 K, and for the facula spectra, they are: 5928, 6028, and 6128 K. We also show the results of using the observed Kitt Peak solar spectra, including the PHOENIX SEDs (using Eqs. (10)–(12)). The μ 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 are expected since the same badly normalized 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.

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, such as, 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, such as 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 modeled 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 μ (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 disk-integrated observations (and therefore by analogy at the disk center) varies significantly among luminosity class and spectral type (see Fig. 17.15 in Gray 2008; Baştürk et al. 2011). From those papers, we can see that within the small range from early G’s dwarfs to early K’s dwarfs, the bisector shapes does not change drastically and, thus, although the integrated spectra for those stars will not be realistic in terms of line bisector, the way that SOAP-GPU models the stellar activity should still be quite realistic as in the end – as 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 estimations of the stellar activity. There is perhaps only one exception for dwarfs cooler than early Ks, 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. These could come from 3D MHD simulations (e.g., Johnson et al. 2021; 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 approaches are challenging, however, they could also lead to a much more realistic modeling of stellar activity from other stars tłian our Sun. We note that it is rattier 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.

Fig. 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 ∆Teff for spot and facula. is fixed to 663 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).

5 Conclusions

In this paper, we present a GPU-based improvement to SOAP 2.0, named SOAP-GPU, that allows us 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 оf resolution, while having the same accuracy.

Besides the huge gain in speed, SOAP-GPU also provides a more complex modeling of stellar activity compared to SOAP 2.0. Complex active region scenarios, wish regions overlapping is now handled by the code. This is mainly useful when modeling active phases of a star such as the Sun, with hundreds of active regions. The contrast of the active regions is now wavelength-dependent and it thus changes for each wavelength of the modeled spectra. The dependence of line bisector with center-to-limb angle μ, 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, and for PHOENIX spectra as well. 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 the 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 the discussion in Sect. 4.3.4). The modeling of stellar activity for other stars than the Sun is currently limited by the knowledge we have about disk-resolved bisectors for these kinds of stars. Such information is very challenging to obtain, however, 3D MHD simulations (e.g., Johnson et al. 2021; Dravins et al. 2021) and resolved spectroscopic observation of other stars due to planetary transits (e.g., Dravins et al. 2017) could help significantly.

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_eff = 5100 K) with spots 663 K cooler and give satisfactory results. Modeling later spectral type is challenging mainly due to continuum normalization 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 the spectral content, we still associate the same active bisector with those regions as that 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 μ angles have been published. Cavallini et al. (1985a) did 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 μ dependant spectral line shape in the quiet and active regions is rather small, and while having more solar data on spots and faculae could allow us to better model the physics at play, the 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 anyone working on techniques to separate the activity effect from planetary signals at the spectral level to test their framework on SOAP-GPU simulations, where the photon-noise, instrumental, and telluric systematics are not perturbing the spectral time series.

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.

Appendix A Line bisectors of PHOENIX spectra

As discussed in Sects. 4.3.2 and 4.3.3, before injecting the μ 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 the 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. A.1, 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 (as well as for simulated spot and faculae) to be 663 K cooler or 250 K hotter, respectively. As expected, most of the bisector are close to straight lines. However, we fit the average bisector with a second order polynomial to remove the small curvatures observed before injecting the proper bisectors at different μ 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 these effects are only due to blends, the process does not significantly change the outputs.

Fig. A.1 Bisector of PHOENIX spectra. For each input seed spectrum using PHOENIX spectral library, we use five strong iron lines: FeI 5250.2084Å (green), FeI 5250.6453Å (cyan), FeI 5434.5232Å (purple), FeI 6173.3344Å (orange), and FeI 6301.5008Å (yellow) to measure the average bisector of the input spectra. Bisector outliers outside a window of 0.1Å around each line center are rejected to avoid those points, certainly affected by line blending, to bias our measurement of line bisector. Each line corresponds 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.

References

All Tables

All Figures

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

In the text

Fig. 2 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.

In the text

	Fig. 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 s to simulate those spectra with SOAP 2.0, while only 27.9 s with SOAP-GPU on a Nvidia RTX-3090 card. The computation speed is improved by a factor of 63.
In the text

Fig. 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.

In the text

Fig. 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 presented 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.

In the text

	Fig. 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 547 840 is labeled in black. The blue and red parts of the spectra with length 204 800 are over plotted in blue and red, respectively.
In the text

Fig. 7 Chromatic effect of RVs. SOAP-GPU is initialized with three different spectra: entire wavelength coverage and only red and only 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: measured RVs of the total effect. Middle: measured RVs of the flux effect. An offset of 1 m s−1 is added in the red and blue RVs. Bottom: measured RVs of the CB effect.

In the text

Fig. 8 Average bisectors of quiet and active solar regions from the disk center (μ = 1.0) to the limb (μ = 0.2). Continuous lines: fifth-order polynomial fit to the quiet sun bisectors of the FeI 5250.2084, FeI 5250.6453, FeI 5434.5232, FeI 5432.9470, FeI 5576.0881, FeI 6149.2460, FeI 6173.3344, FeI 6301.5008, and FeI 6302.4932 Å 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 Å 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−1 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 μ angles are all shifted by those 340 m s−1 at a depth of 0.58 based on our hypothesis that convection is fully suppressed in magnetic regions (see Sect. 4.3.2 for more information).

In the text

	Fig. 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 Å bisector because it was significantly off by 2500 m s−1 due to strong contamination by other weak lines.
In the text

Fig. 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 μ dependent bisector injection. Dotted line: RVs derived from simulated spectra with μ dependent bisector injection. The input spectra with μ-angle dependency are generated with the Python module Convec.model. The RVs of the flux, CB, and total effect do not change significantly when the μ-dependent CB is introduced.

In the text

Fig.11 Faculae intensity contrast as a function of μ for different spectral types. The limb brightening curves for the G2, K0 and M0 dwarfs are derived from the parametrization of MHD simulations for 500G faculae, shown in Table 1 in Johnson et al. (2021). 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 obtained the 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.

In the text

Fig. 12 CB velocity as a function of spectral type. Top: data from Liebing et al. (2021) and cubic fit to them, giving as relation CBvel = 95.2388 × ((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: bisector of the quiet photosphere measured from Löhner-Böttcher et al. (2019) at μ = 0.85 (thin lines) and the bisector of an active region measured from Cavallini et al. (1985a) at μ = 1.0 (thick lines) for different spectral types. We impose that the maximum difference between the quiet bisector at μ = 0.85 and μ = 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 μ = 0.85 and the active bisectors, independent of their μ angle, to be equal to the same value at a depth of 0.58.

In the text

Fig. 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 Teff = 5778 K. The effective temperature of the spot spectra are 5015, 5115, and 5215 K, and for the facula spectra, they are: 5928, 6028, and 6128 K. We also show the results of using the observed Kitt Peak solar spectra, including the PHOENIX SEDs (using Eqs. (10)–(12)). The μ 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 are expected since the same badly normalized 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 the text

Fig. 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 ∆Teff for spot and facula. is fixed to 663 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).

In the text

Fig. A.1 Bisector of PHOENIX spectra. For each input seed spectrum using PHOENIX spectral library, we use five strong iron lines: FeI 5250.2084Å (green), FeI 5250.6453Å (cyan), FeI 5434.5232Å (purple), FeI 6173.3344Å (orange), and FeI 6301.5008Å (yellow) to measure the average bisector of the input spectra. Bisector outliers outside a window of 0.1Å around each line center are rejected to avoid those points, certainly affected by line blending, to bias our measurement of line bisector. Each line corresponds 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.

In the text