# Using Doppler imaging to model stellar activity and search for planets around Sun-like stars

Baptiste Klein, Suzanne Aigrain, Michael Cretignier, Xavier Dumusque, Khaled Al Moulla, Jean-François Donati, Niamh K. O'Sullivan, Haochuan Yu, Andrew Collier Cameron, Oscar Barragán, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Baptiste Klein, Suzanne Aigrain, Michael Cretignier, Xavier Dumusque, Khaled Al Moulla, et al. Using Doppler imaging to model stellar activity and search for planets around Sun-like stars. Monthly Notices of the Royal Astronomical Society, 2025, 542 (4), pp.2714-2732. 10.1093/mnras/staf1337. hal-05376965

HAL Id: hal-05376965

https://cnrs.hal.science/hal-05376965v1

Submitted on 22 Nov 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Using Doppler imaging to model stellar activity and search for planets around Sun-like stars

Baptiste Klein<sup>1,★</sup> Suzanne Aigrain,<sup>1</sup> Michael Cretignier<sup>1,†</sup> Xavier Dumusque,<sup>2</sup> Khaled Al Moulla<sup>1,‡</sup>,<sup>3</sup> Jean-François Donati<sup>3,‡</sup>,<sup>4</sup> Niamh K. O'Sullivan,<sup>1</sup> Haochuan Yu<sup>1,‡</sup>,<sup>5</sup> Andrew Collier Cameron<sup>1,§</sup>,<sup>5</sup> Oscar Barragán<sup>1,6</sup>,<sup>7</sup> Annelies Mortier<sup>1,7</sup> and Alessandro Sozzetti<sup>1,8</sup>

Accepted 2025 August 7. Received 2025 August 6; in original form 2025 May 13

## ABSTRACT

Doppler imaging (DI) is a well-established technique to map a physical field at a stellar surface from a time series of high-resolution spectra. In this proof-of-concept study, we aim to show that traditional DI algorithms, originally designed for rapidly rotating stars, also have the ability to model the activity of Sun-like stars, when observed with new-generation highly stable spectrographs, and search for low-mass planets around them. We used DI to retrieve the relative brightness distribution at the surface of the Sun from radial velocity (RV) observations collected by HARPS-N (High Accuracy Radial velocity Planet Searcher for the Northern hemisphere) between 2022 and 2024. The brightness maps obtained with DI have a typical angular resolution of ~36° and are a good match to low-resolution disc-resolved Dopplergrams of the Sun at epochs when the absolute, disc-integrated RV exceeds ~2 m s⁻<fnref n="1" />. The RV residuals after DI correction exhibit a dispersion of about 0.6 m s⁻<fnref n="1" />, comparable with existing state-of-the-art activity correction techniques. Using planet injection-recovery tests, we also show that DI can be a powerful tool for blind planet searches, so long as the orbital period is larger than ~100 d (i.e. 3–4 yr stellar rotation periods), and that it yields planetary mass estimates with an accuracy comparable to, for example, multidimensional Gaussian process regression. Finally, we highlight some limitations of traditional DI algorithms, which should be addressed to make DI a reliable alternative to state-of-the-art RV-based planet search techniques.

**Key words:** line: profiles – techniques: radial velocities – Sun: activity – planets and satellites: detection – stars: activity.

## 1 INTRODUCTION

The detection and characterization of temperate Earth-mass planets around Sun-like stars is one of the most exciting and challenging prospects in contemporary astrophysics. New-generation extreme-precision spectrographs like ESPRESSO (Echelle SPECTograph for Rocky Exoplanets and Stable Spectroscopic Observations; Pepe et al. 2021; Faria et al. 2022), EXPRES (EXtreme PREcision Spectrograph; Blackman et al. 2020; Brewer et al. 2020; Petersburg et al. 2020), and NEID (NN-explore Exoplanet Investigations with Doppler spectroscopy; Schwab et al. 2018; Gupta et al. 2025) routinely deliver radial velocities (RVs) with sub-metre per second precision. With ongoing and forthcoming long-baseline planet search programs, such as the PLATO (PLAnetary Transits and Oscillations of stars) mission (including its RV follow-up; Rauer et al. 2014),

the NEID Earth Twin Survey (Gupta et al. 2021, 2025), or the Terra Hunting Experiment (Thompson et al. 2016; Hall et al. 2018), the detection of an Earth twin in the solar neighbourhood seems to be within reach.

Stellar variability is the main obstacle to such a detection (Fischer et al. 2016; Crass et al. 2021; Meunier et al. 2023). Photospheric flows (e.g. p-mode oscillations, supergranulation, and meridional circulation; see Kjeldsen & Bedding 1995; Dumusque et al. 2011a; Meunier et al. 2015; Chaplin et al. 2019; Meunier & Lagrange 2019; O'Sullivan & Aigrain 2024; O'Sullivan et al. 2025) and magnetic activity (e.g. magnetically active regions and magnetic cycles; see Saar & Donahue 1997; Desort et al. 2007; Meunier, Desort & Lagrange 2010; Dumusque et al. 2011b; Meunier, Mignon & Lagrange 2017) distort the stellar absorption line profiles, giving rise to RV signals that mask planet signatures.

The modelling of stellar activity RV signals has become a very active field of research in the last decade (see Zhao et al. 2022).

<footnote n="1"><sup>1</sup>Department of Physics, University of Oxford, Oxford OX1 3RH, UK</footnote>

<footnote n="2"><sup>2</sup>Observatoire Astronomique de l'Université de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland</footnote>

<footnote n="3"><sup>3</sup>Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, P-4150-762 Porto, Portugal</footnote>

<footnote n="4"><sup>4</sup>CNRS, IRAP, Université de Toulouse, 14 avenue Belin, F-31400 Toulouse, France</footnote>

<footnote n="5"><sup>5</sup>Centre for Exoplanet Science/SUPA, School of Physics & Astronomy, University of St Andrews, North Haugh, St Andrews, Fife KY16 9SS, UK</footnote>

<footnote n="6"><sup>6</sup>Department of Physics, University of Warwick, Coventry CV4 7AL, UK</footnote>

<footnote n="7"><sup>7</sup>School of Physics & Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK</footnote>

<footnote n="8"><sup>8</sup>INAF – Osservatorio Astrofisico di Torino, via Osservatorio 20, I-10025 Pino Torinese, Italy</footnote>

<footnote>* E-mail: baptiste.klein@physics.ox.ac.uk</footnote>
These signals are often modelled in the time domain with Gaussian processes (GPs; see Aigrain & Foreman-Mackey 2023, for a review). In particular, multidimensional GPs, where RV signals are modelled jointly with activity proxies, provide a trackable, robust, and flexible stellar activity modelling framework, often considered as state of the art in the community (e.g. Rajpaul et al. 2015; Barragán et al. 2022b, 2023; Delisle et al. 2022; Hara & Delisle 2025). Despite their popularity, GPs suffer from key limitations in the search for low-mass long-period planets. In general, current GP models lack the flexibility to capture the RV signals of stars with rapidly evolving activity (i.e. where active regions evolve on time-scales similar to the star’s rotation period), which can lead to inaccuracies in the estimates of the planet parameters (e.g. Blunt et al. 2023). Moreover, magnetic cycles induce significant variations in the properties of stellar activity signals (e.g. rotation period and evolution time-scale, due to changes in the filling factor of active regions; Klein et al. 2024), which are generally not accounted for in GP covariance functions.

To overcome these limitations, the community is increasingly focusing on wavelength-domain activity filtering techniques. Approaches exploiting activity-induced changes in the shape of spectral lines have been shown to robustly filter cycle-induced activity variations, increasing the RV sensitivity to long-period planet signatures (e.g. Jones et al. 2022; Collier Cameron et al. 2021; de Beurs et al. 2022; John, Collier Cameron & Wilson 2022; Wilson et al. 2022; Cretignier et al. 2023; John et al. 2023; Klein et al. 2024; Yu et al. 2024; Zhao et al. 2024). On the other hand, by leveraging the wealth of information delivered by new-generation échelle spectrographs, line-by-line techniques enable a more detailed physical characterization of stellar activity, which is known to affect different spectral lines in distinct ways (e.g. Thompson et al. 2017; Dumusque 2018; Cretignier et al. 2020a, 2022; Al Moulla et al. 2022; Bellotti et al. 2022; Cretignier, Dumusque & Pepe 2022; Lienhard et al. 2022, 2023; Al Moulla, Dumusque & Cretignier 2024; Rescigno & Al Moulla 2025).

In this study, we investigate how Doppler imaging (DI; see Vogt, Penrod & Hatzes 1987; Donati & Landstreet 2009; Kochukhov 2016) can be used to mitigate stellar activity and search for planet signatures in the spectra of Sun-like stars. This physically motivated technique aims at inverting a time series of high-resolution spectra into a distribution of a physical field at the stellar surface (e.g. relative brightness, temperature, magnetic field, and chemical composition). Historically, the mapping of brightness inhomogeneities from intensity spectra, which this study focuses on, has been limited to rapidly rotating stars such as giants and subgiants (e.g. Donati 1999; Strassmeier 1999; Roettenbacher et al. 2017), T Tauri stars (e.g. Collier Cameron & Unruh 1994; Hatzes 1995; Donati et al. 2000, 2012, 2014; Yu et al. 2019; Finociety et al. 2021), stars in close-in binary systems (e.g. Zaire, Donati & Klein 2021, 2022), and even brown dwarfs (Crossfield et al. 2014; Luger et al. 2021a). New-generation highly stable high-resolution spectrographs have the ability to resolve the profile distortions induced by active regions at the surface of more slowly rotating stars (as already showcased in Donati, Henry & Hall 1995; Hébrard et al. 2016), which suggests that DI could be applied to Sun-like stars, provided that the signal-to-noise ratio (SNR) of the spectra is high enough, and that the temporal sampling is sufficiently dense.

This proof-of-concept study aims to show that traditional DI cannot only be applied to Sun-like stars, but also provide a reliable alternative to time-domain activity modelling techniques in the search for low-mass planetary signatures. In Section 2, we give a short review of DI techniques in stellar physics and describe our framework to apply it to Sun-like stars. In Section 3, we describe

the High Accuracy Radial velocity Planet Searcher for the Northern hemisphere (HARPS-N) Sun-as-a-star observations used in this case study, before modelling them with DI and comparing the results to resolved images of the Sun in Section 4. Finally, we assess the sensitivity of DI to planet signatures in Section 5, before discussing the current limitations of our framework and avenues to address them in Section 6.

# 2 METHOD

## 2.1 Doppler imaging in stellar physics

DI is a long-standing method to map inhomogeneous structures on stellar surfaces from high-resolution spectra. DI has been historically used for the chemical mapping of Ap and Bp stars, which exhibit strong chemically peculiar absorption lines with periodic variations (e.g. Goncharskij et al. 1983; Vogt & Penrod 1983; Khokhlova, Rice & Wehlau 1986; Vogt et al. 1987; Hatzes 1991; Kochukhov et al. 2004, 2007). Zeeman-Doppler imaging (ZDI), an extension of DI to map stellar magnetic field topologies from polarized spectra, is one of the most successful stellar imaging techniques (see Donati & Landstreet 2009; Kochukhov 2016, for a review). ZDI has been applied to a very wide range of stars, including massive early-type stars (e.g. Donati et al. 2002, 2006a), chemically peculiar stars (e.g. Landstreet & Mathys 2000; Bagnulo et al. 2002; Kochukhov & Wade 2010), RS CVn stars (e.g. Donati 1999; Donati et al. 2003; Kochukhov et al. 2013), classical and weak-line T Tauri stars (e.g. Donati et al. 2008, 2012; Folsom et al. 2016; Yu et al. 2017, 2019), Sun-like stars (Petit et al. 2008, 2021; Marsden et al. 2014; Folsom et al. 2020), and M dwarfs (e.g. Morin et al. 2008, 2010; Klein et al. 2021a; Lehmann et al. 2024). The output magnetic topologies are excellent showcases of the dynamo processes in stellar interiors (e.g. Donati et al. 2003; Gastine et al. 2013; Brun & Browning 2017), and the evolution of the field's properties with time is one of the most direct ways to unambiguously unveil magnetic cycles (e.g. Fares et al. 2009; Boro Saikia et al. 2016, 2018; Lavail, Kochukhov & Wade 2018; Lehmann et al. 2021; Bellotti et al. 2025). Moreover, magnetic topologies are essential ingredients for modelling magnetized winds and their interaction with close-in planets (e.g. Kavanagh et al. 2021; Callingham et al. 2024; Strugarek & Shkolnik 2025).

In this study, we focus on the mapping of brightness inhomogeneities from intensity spectra. This method has historically been effective with rapidly rotating active low-mass stars, such as classical and weak-line T Tauri stars, members of close-in binary systems (e.g. RS CVn stars), and rapidly rotating giant stars (e.g. K Com stars). In contrast with temperature mapping, which uses a small number of spectral lines to constrain the spot properties (see Berdyugina 2005; Strassmeier 2009; Afram & Berdyugina 2015), brightness reconstructions focus on retrieving a fractional spot coverage, losing the absolute information on the spot temperature, but allowing one to combine thousands of spectral lines at the same time (e.g. through least-squares deconvolution; Donati et al. 1997). Mapping the star's relative brightness is a way of modelling stellar activity, bypassing the RV computation process. In fact, DI has been used to jointly model the activity of pre-main-sequence stars whilst searching for planetary signatures in the recent literature (e.g. Petit et al. 2015; Donati et al. 2017; Yu et al. 2017; Klein et al. 2021b, 2022).

DI is an ill-posed inverse problem, which means that multiple solutions can fit the same data set equally well. To lift these degeneracies, the problem is generally regularized using, for example, Tikhonov regularization (e.g. Khokhlova et al. 1986; Piskunov, Tuominen & Vilhu 1990), or the maximum entropy method (Vogt & Penrod 1983;
Skilling & Bryan 1984; Vogt et al. 1987). As pointed out by Luger et al. (2021a), these methods, albeit well understood (e.g. Donati & Brown 1997; Kochukhov & Piskunov 2002; Piskunov & Kochukhov 2002; Rice 2002), are both known to alter the reconstructed maps (e.g. with spurious artefacts on smaller spatial scales) and are more sensitive to longitudinal than latitudinal variations, especially for near-equator-on stars. Bayesian implementations of the DI algorithm have been proposed recently and offer a promising route to account for the framework's inherent degeneracies and their effect on the inferred stellar maps (e.g. Luger et al. 2021a; Asensio Ramos, Díaz Baso & Kochukhov 2022).

## 2.2 Description of the framework

In this study, we use the DI code described in Semel (1989), Donati, Semel & Praderie (1989), Brown et al. (1991), Semel, Donati & Rees (1993), Donati & Brown (1997), Donati et al. (2006b, 2014), which is designed to perform Doppler and Zeeman-Doppler analyses of unpolarized and polarized spectra, using maximum-entropy regularization. This DI framework has been validated on various types of low-mass stars, such as classical and weak-line T Tauri stars (e.g. Donati et al. 2000, 2012, 2014; Yu et al. 2017, 2019; Finociety et al. 2021), and K/M dwarfs (e.g. Hébrard et al. 2016; Klein et al. 2021b; Zaire et al. 2021, 2022).

Our DI framework describes the star's surface as a grid of cells (typically about 10 000 in this study). Each cell $i$ is assigned a relative brightness coefficient $b_i$ (in [0,1] for a dark spot and >1 in the case of a bright plage) and a local RV. In the direct approach, the code computes, at all epochs of observation, a local line profile for each grid cell using the Unno-Rachkovsky analytical solution to the radiative transfer equations in a Milne-Eddington atmosphere<fnref n="1" /> (Unno 1956; Landi Degl'Innocenti & Landolfi 2004). At a given time stamp (i.e. rotational phase), the algorithm computes the sky-projected area $\omega_i$ and the RV of each cell of the star. Disc-integrated line profiles are computed by averaging the local profiles, shifted at the local RV, and weighted by the local limb-darkening coefficient, projected cell area, and relative brightness. The limb-darkening coefficients are computed following the linear law $1 - \epsilon(1 - \mu)$, where $\epsilon$ is the limb-darkening coefficient (typically taken from Claret, Hauschildt & Witte 2013) and $\mu$ is the cosine of the limb angle.

In the inverse approach, our DI framework starts from an unspotted star (i.e. $b_i = 1$ for all the cells) and uses a conjugate-gradient process to iteratively compute relative brightness maps until the synthetic line profiles (computed from the model stellar surface) match the observed profiles down to a given level of $\chi^2$. The degeneracy of this inversion problem is lifted by imposing a maximum-entropy regularization condition on our estimator (Skilling & Bryan 1984). In our case, the entropy $S$ is defined by the Shannon entropy of the relative brightness of the cells:

$$S = - \sum_{i}^{N_c} \omega_i b_i \left[ \log \left( \frac{b_i}{b_0} \right) - 1 \right], \quad (1)$$

where $N_c$ is the total number of cells at the surface of the star. Conceptually, this can be seen as choosing the brightness distribution with the filling factor $f_{DI}$ of active regions that can reproduce the observed profiles at a given level of $\chi^2$. This filling factor is defined as the sum of the filling factors in dark spots $(\sum_i \omega_i(1 - b_i)/N_c)$ and in bright plages $(\sum_i \omega_i(b_i - 1)/N_c)$.

### 2.2.1 Applying Doppler imaging to slow rotators

Applying DI to slowly rotating stars (which we define as having a rotational velocity $v \sin i \lesssim 10$ km s$^{-1}$) is more complicated for two main reasons. The first, and most intuitive, is that the angular size of the resolution element at the stellar surface scales as $(Rv \sin i)^{-1}$, where $R$ is the resolving power of the instrument (see section 9.2.2. of Kochukhov 2016). Therefore, in principle, a single snapshot can only probe the largest spatial scales of the brightness variations. This limitation can none the less be overcome with a dense temporal sampling of the star with highly stable high-precision spectrographs, capturing the evolution of the activity-induced distortions as the star rotates and the active regions evolve. In other words, as long as (i) the SNR of the line profiles is large enough to resolve the activity-induced distortions, and (ii) enough observations are collected on the time-scales over which the star's activity evolves, DI can be, in principle, applied to very slow rotators (e.g. down to $v \sin i$ as low as 1 km s$^{-1}$ in Hébrard et al. 2016).

The other and the most challenging limitation comes from the fact that the intrinsic profile (i.e. the local profile within each DI cell) is unknown. Unlike fast rotators, whose line profiles are generally dominated by the star's rotational broadening, the shape of the disc-integrated line profiles of slow rotators encompasses the effects of formation temperature, microturbulence, granulation, magnetic field, and rotation, which makes it impossible to predict the intrinsic profile to the level of the photon noise. We must thus assume a local line profile (in the Unno-Rachkovsky's framework), which will likely exhibit systematic differences with the actual local line profile of slow rotators like the Sun (e.g. because the 3D effects of convection are ignored). In this case, our DI code will naturally focus on reconstructing these systematic differences rather than the activity-induced distortions, leading to unrealistic brightness distributions. To overcome this problem, we use the iterative process introduced in Klein et al. (2021b, 2022). We use our DI code to perform a maximum-entropy fit of all the observed cross-correlation functions (CCFs), $I_{obs}$. We subtract the median difference between $I_{obs}$ and $I_{syn}$, the best-fitting synthetic profiles, from $I_{obs}$ and repeat the procedure until the median difference between $I_{obs}$ and $I_{syn}$ is flat (i.e. with a dispersion significantly lower than the photon noise). This process should preserve most of the activity-induced CCF distortions (as confirmed in Section 4), but could affect axisymmetric structures generating non-modulated profile distortions though.

### 2.2.2 The specific case of Sun-like stars

Since active regions at the surface of the Sun evolve on a time-scale similar to the rotation period (e.g. Leighton 1964; Foukal 1998; Klein et al. 2024), the solar disc will likely look dramatically different from one solar rotation to the next. Our DI framework assumes that active regions do not evolve intrinsically, to avoid adding further

<footnote n="1">¹Note that the Unno-Rachkovsky framework is used by default to compute local line profiles in our DI code. Equally good results would be obtained with more simplistic models like Voigt profiles since (i) the Sun is slowly rotating and, thus, the intrinsic line profile is unknown (see Section 2.2.1), (ii) we are working with CCFs, whose link with the physical properties of the stellar atmosphere remains obscure, and (iii) we are only interested in the relative variations of line profiles.</footnote>
degeneracies to the model.<fnref n="2" /> In this study, we therefore apply DI on data covering no more than a 27-d rotation cycle, as described in Section 3.2.

Furthermore, our DI code assumes that the profile distortions are only due to brightness inhomogeneities at the stellar surface. This means that the inhibition of the convective blueshift within faculae, which dominates the solar activity RV budget (e.g. Meunier et al. 2010), is approximated by brightness contrast effects in our model. The flexibility of DI will naturally reproduce the effects of faculae on the line profiles, even if it means that the resulting brightness distribution loses some of its realism. This will not be a major issue for this study, which focuses on the mitigation of stellar activity signals for RV planet searches. None the less, the realism of the DI maps will be assessed in Section 4.2. Additionally, note that our code does not include the effects of granulation and supergranulation (e.g. Meunier et al. 2015).

# 3 INPUT DATA AND PRE-PROCESSING

## 3.1 Observations and data processing

The cross-dispersed échelle spectrograph HARPS-N, at the 3.58-m Telescopio Nazionale Galileo (TNG) at Roque de Los Muchachos Observatory (La Palma, Spain), has been monitoring the Sun since 2015 (Cosentino et al. 2012; Dumusque et al. 2015, 2021; Phillips et al. 2016; Collier Cameron et al. 2019; Dumusque et al., in preparation). Disc-integrated solar spectra spanning most of the optical domain (383–690 nm) at a resolving power of $R = 115\,000$ are collected at a 5-min cadence during day time. Our input data set contains the observations acquired between 2021 December and 2024 January (i.e. in the start of solar cycle 25), reduced with version 3.0.1 of the ESPRESSO data reduction software (DRS; Dumusque et al. 2021; Pepe et al. 2021). Working with these 2 yr of data is motivated by two main reasons. First, the Sun is particularly active and the RV signals exhibit clear quasi-periodic modulations, well suited for DI. Secondly, the data were collected after the refurbishment of the instrument’s CCD camera in 2021, and no major maintenance operations were carried out on the instrument since then.

We use the same approach as Klein et al. (2024) to build the input data set. Following Collier Cameron et al. (2019), we only select epochs for which the probability that the Sun is partially masked by clouds in the Earth’s atmosphere is less than 10 per cent. Similarly, observations with velocity corrections larger than 0.1 m s⁻¹ were discarded to ensure that differential extinction is minimum (i.e. that the Sun is close to the zenith). Daily-stacked S1D spectra are continuum normalized with the open-source software RASSINE (Cretignier et al. 2020b) and corrected from known instrumental contamination (i.e. ghosts, stitchings, interference patterns, ThAr bleeding, defocus of the point spread function, and tellurics) using the software YARARA (Cretignier et al. 2021).

CCFs are computed using the line list tailored for the Sun described in Cretignier et al. (2020a, 2022). The RVs, full width at half-maximum (FWHM) and bisector velocity span ($V_s$; Queloz et al. 2001) are extracted from the derived CCFs. Our final input data set contains 511 daily-binned CCFs, collected between 2021 December and 2024 January (2.1 yr) and sampled at the resolution of the

instrument (0.82 km s⁻¹; Dumusque et al. 2021). The typical photon noise in the CCF continuum is about $2 \times 10^{-5}$, corresponding to a typical RV uncertainty of 0.11 m s⁻¹.

## 3.2 Data selection

As described in Section 2.2, we split our input data set into 27-d subsets (hereafter called chunks) that we analyse independently in the rest of the paper. We discard all chunks with either fewer than 10 observations or with a gap larger than 25 per cent in the Sun’s rotational phase. These assumptions, discussed in more detail in Section 4.3, are motivated by the need for a dense temporal sampling, in order to capture the activity-induced profile distortions as they are modulated by the Sun’s rotation. Our input data set, shown in Fig. 1, contains 439 observations spread into 23 chunks, with an average of 19 points per chunk. Information on the different chunks is given in Table 1. Finally, the time $t$ of each observation is also expressed in units of solar rotation phase $\phi_{rot}$ using

$$ \phi_{rot} = \frac{t - T_0}{P_{rot}}, \qquad (2) $$

where $P_{rot}$ refers to the Carrington solar rotational period, set to 27.2753 d, and $T_0 = 2457\,220.55$ is a reference time, arbitrarily chosen as the first epoch of the HARPS-N data set. Note also that this value of $P_{rot}$ is a good match with the Sun’s rotation period measured from the same data set in Klein et al. (2024).

# 4 RESULTS

## 4.1 Doppler imaging maps

We apply DI to the continuum-normalized HARPS-N solar CCFs within each 27-d chunk defined in Table 1. The inversion process is performed assuming a solar inclination of 80° to avoid including additional north–south degeneracies in our fit. This assumption, which is discussed in Section 4.4, should not have a major effect on the reconstructed maps since DI is generally not very sensitive to variations of stellar inclinations of 10°–20°, especially for slowly rotating stars. As a conservative approach, we use CCF continuum dispersion as formal CCF uncertainties, in order to account for sources of noise, such as the instrument instability and stellar variability (e.g. supergranulation), not accounted for in the formal CCF continuum uncertainty. Our final uncertainties are typically around 0.01 per cent, which is about three times smaller than the rms of the observed activity-induced distortions in the core of the CCFs. After computing the best-fitting intrinsic profile using the method described in Section 2.2, the DI fit systematically converges to a reduced $\chi^2$ of 1.

The best-fitting brightness maps, shown in orthographic projection in Fig. 2, exhibit filling factors ranging from 0.09 per cent to 0.57 per cent (see Fig. A1 for the relative brightness distribution on the full solar surface). We also show the best fit to the line profiles in Figs A2 and A3. To estimate the typical variation of the filling factor in active regions induced by noise, we create a series of mock data sets by adding Gaussian white noise, drawn from the uncertainties adopted on the CCFs (see Table 1), to the best-fitting profiles of our model. We then apply DI to the mock CCF time series. For each chunk, we repeat the process for 10 white noise realizations and found a typical relative variation of 20 per cent on the filling factor. The robustness of the reconstructed maps is discussed in more detail in Section 4.3.

<footnote n="2">²Note that recent efforts to account for the evolution of stellar magnetic fields using GPs in ZDI reconstructions have shown promising results (Finociety & Donati 2022).</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Solar rest-frame RV time series extracted with HARPS-N DRS. The vertical colour bands indicate the different subsets of data used in the DI analysis (see Section 3.2). Grey points were not used in the analysis. The number of each chunk as listed in Table 1 is indicated in dark blue at the bottom of the figure.</figcaption></figure>

<table><caption>**Table 1.** List of relevant quantities for each of the 23 chunks defined in Section 3.2 and shown in Fig. 1. Columns 2 and 3 give the day and BJD of the epoch corresponding to a rotation phase of zero in equation (2), which is also the phase of display of the brightness maps in Fig. 2. Columns 4–7 give the number of points $N_{pt}$, the RV dispersion $\sigma_{RV}$, the average continuum dispersion in the CCFs $\sigma_{CCF}$, and the recovered active region filling factor $f_{DI}$ for each chunk, respectively.</caption><thead><tr><th>Chunk</th><th>Date</th><th>BJD</th><th>N<sub>pt</sub></th><th>σ<sub>RV</sub><br>(m s<sup>-1</sup>)</th><th>σ<sub>CCF</sub><br>(per cent)</th><th>f<sub>DI</sub><br>(per cent)</th></tr></thead><tbody><tr><td>1</td><td>2021-12-17</td><td>59566.23</td><td>17</td><td>1.39</td><td>0.008</td><td>0.20</td></tr><tr><td>2</td><td>2022-01-14</td><td>59593.50</td><td>13</td><td>0.88</td><td>0.007</td><td>0.57</td></tr><tr><td>3</td><td>2022-02-10</td><td>59620.78</td><td>14</td><td>1.13</td><td>0.007</td><td>0.36</td></tr><tr><td>4</td><td>2022-03-09</td><td>59648.05</td><td>17</td><td>1.56</td><td>0.008</td><td>0.35</td></tr><tr><td>5</td><td>2022-04-05</td><td>59675.33</td><td>22</td><td>1.34</td><td>0.010</td><td>0.34</td></tr><tr><td>6</td><td>2022-05-03</td><td>59702.60</td><td>23</td><td>1.80</td><td>0.010</td><td>0.26</td></tr><tr><td>7</td><td>2022-05-30</td><td>59729.88</td><td>19</td><td>2.02</td><td>0.012</td><td>0.34</td></tr><tr><td>8</td><td>2022-06-26</td><td>59757.15</td><td>25</td><td>1.40</td><td>0.010</td><td>0.43</td></tr><tr><td>9</td><td>2022-07-23</td><td>59784.43</td><td>22</td><td>1.93</td><td>0.012</td><td>0.37</td></tr><tr><td>10</td><td>2022-08-20</td><td>59811.70</td><td>20</td><td>1.02</td><td>0.007</td><td>0.35</td></tr><tr><td>11</td><td>2022-10-13</td><td>59866.25</td><td>16</td><td>1.66</td><td>0.011</td><td>0.39</td></tr><tr><td>12</td><td>2022-11-10</td><td>59893.53</td><td>19</td><td>1.47</td><td>0.007</td><td>0.29</td></tr><tr><td>13</td><td>2023-02-27</td><td>60002.63</td><td>11</td><td>1.91</td><td>0.009</td><td>0.40</td></tr><tr><td>14</td><td>2023-03-26</td><td>60029.91</td><td>20</td><td>1.51</td><td>0.010</td><td>0.24</td></tr><tr><td>15</td><td>2023-04-22</td><td>60057.18</td><td>25</td><td>1.19</td><td>0.006</td><td>0.29</td></tr><tr><td>16</td><td>2023-05-19</td><td>60084.46</td><td>23</td><td>1.74</td><td>0.009</td><td>0.20</td></tr><tr><td>17</td><td>2023-06-16</td><td>60111.73</td><td>23</td><td>2.76</td><td>0.014</td><td>0.32</td></tr><tr><td>18</td><td>2023-07-13</td><td>60139.01</td><td>27</td><td>1.98</td><td>0.009</td><td>0.22</td></tr><tr><td>19</td><td>2023-08-09</td><td>60166.28</td><td>18</td><td>0.73</td><td>0.005</td><td>0.09</td></tr><tr><td>20</td><td>2023-09-06</td><td>60193.56</td><td>19</td><td>1.09</td><td>0.006</td><td>0.16</td></tr><tr><td>21</td><td>2023-10-03</td><td>60220.83</td><td>17</td><td>1.68</td><td>0.008</td><td>0.31</td></tr><tr><td>22</td><td>2023-10-30</td><td>60248.11</td><td>16</td><td>1.13</td><td>0.007</td><td>0.29</td></tr><tr><td>23</td><td>2023-11-26</td><td>60275.38</td><td>13</td><td>2.88</td><td>0.016</td><td>0.56</td></tr></tbody></table>

As a sanity check, we extract the RVs from the best-fitting CCFs and compare them to the solar RVs extracted from HARPS-N DRS in Fig. 3. We find that both RV time series match well, with a residuals root mean square (rms) of 0.58 m s⁻¹, compared to a dispersion of 2.2 m s⁻¹ in the input RVs. Using generalized Lomb–Scargle (GLS) periodograms (Zechmeister & Kürster 2009), we found no periodicity in the RV residuals (see Fig. 4). In particular, most power at the Sun’s rotation period and first harmonic has been removed in the residuals, which indicates that the rotationally modulated activity component is well captured by our model. The rms of the RV residuals remains significantly larger than the formal RV uncertainties (~0.1 m s⁻¹) but is consistent with that obtained with a GP regression on the RV time series (see Klein et al. 2024). Since

oscillation- and granulation-induced RV variations have been, for the most part, averaged out in our daily-binned data (Dumusque et al. 2011a; Chaplin et al. 2019), the dispersion budget in the RV residuals is most likely a mix of instrument stability (~0.5 m s⁻¹; Dumusque et al. 2021) and supergranulation signals (estimated at around 0.7 m s⁻¹, but likely reduced by the daily binning process; Meunier et al. 2015; Al Moulla et al. 2023; Lakeland et al. 2024).

## 4.2 Comparison with SDO observations

A clear advantage of using Sun-as-a-star observations is that our best-fitting brightness maps of Fig. 2 can be compared with resolved images of the Sun. The Helioseismic and Magnetic Imager (HMI;
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Best-fitting relative brightness distribution of the Sun for each of the 23 chunks listed in Table 1 and shown in Fig. 1. The colour scale depicts the logarithm of the relative surface brightness. All maps are shown at a solar rotation phase of 0, computed using equation (2), and the chunk number is indicated on top of each map.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Time series of RVs extracted from HARPS-N DRS (black points) and from the best-fitting DI line profiles (red solid lines, top panel), and residuals (bottom panel).</figcaption></figure>

Scherrer et al. 2012; Schou et al. 2012) of the *Solar Dynamics Observatory* (SDO; Pesnell, Thompson & Chamberlin 2012) has been collecting high-cadence high-resolution observations of the solar surface since 2010. We used the open-source pipeline SOLASTER (Ervin et al. 2022) to download the HMI intensitygrams, Dopplergrams, and magnetograms at the epochs corresponding to a rotational phase of zero as defined in equation (2) (i.e. the same phase as the maps shown in Fig. 2), for each of the 23 chunks. The Dopplergrams are corrected for the Sun's rotational velocity and spacecraft motion, and intensitygrams are corrected for limb darkening. To make sure that the SDO data are not affected by short-term variability (e.g. p-mode oscillations and granulation), we average 10 SDO maps

collected on a 10-h window centred on the Sun's rotation phase of interest.

We model the SDO Dopplergrams, magnetograms, and mean-subtracted intensitygrams using spherical harmonics (SH) of maximum degree $l_{max} = 5$, corresponding to the typical spatial resolution probed by our DI algorithm. Hereafter, we refer to the SH model of the SDO maps as large-scale SDO maps. As a starting point, we compare the orthographic projection of our DI brightness maps at rotational phase 0 (as shown in Fig. 2) to the large-scale SDO outputs at the same phase. We find that the Pearson correlation coefficient between the large-scale SDO Dopplergrams and the DI brightness maps varies significantly from one chunk to the next, with
<figure><img src="image_6.png" /></figure>

**Figure 4.** GLS periodogram of the HARPS-N RVs (top, black line) and of the RV residuals (bottom, yellow line). The false alarm probability of the most prominent peak in the periodogram of the RV residuals, computed using the method of Baluev (2008), is about 0.5. The vertical dotted lines indicate the Sun's rotation period (27.2753 d here) and first two harmonics. The periodograms were computed using the ASTROPY PYTHON module (Astropy Collaboration 2013, 2018, 2022).

<figure><img src="image_7.png" /></figure>

**Figure 5.** Correlation coefficients between SDO large-scale Dopplergrams (DSDO) and DI brightness maps (IDI), against the solar RVs at rotational phase 0, as defined in equation (2). The two time series exhibit a Pearson correlation coefficient ρRV of 0.85.

values as high as 0.76, but with a mean absolute value of 0.34 and standard deviation of 0.22. Correlation coefficients with the large-scale SDO intensitygrams and magnetograms do not reach absolute values greater than 0.5.

As shown in Fig. 5, the correlation coefficients ρ between the DI maps and SDO Dopplergrams are themselves correlated with the RV (respectively FWHM) at phase 0, with a Pearson correlation coefficient of 0.85 (respectively 0.75). We interpret this as follows. Line profile distortions inducing RV signals of absolute value smaller than ~2 m s⁻¹ have an amplitude comparable to the dispersion of the CCF continuum and, therefore, are not well separated from the noise in the DI inversion. On the other hand, the positions of active regions inducing RV signatures with absolute values larger than ~2 m s⁻¹ are accurately recovered within ~36° (i.e. the typical angular resolution

<figure><img src="image_8.png" /></figure>

**Figure 6.** Qualitative comparison between SDO Dopplergrams and DI brightness maps at the rotational phase that maximizes the absolute RV value within chunks 1, 7, and 11, as defined in Table 1. From left to right, we show the SDO Dopplergram corrected from rotational flows, the SH model of the SDO Dopplergram with $ℓ_{max}$ = 5, and the best-fitting relative brightness maps obtained with DI. The colour bars are defined between ±0.5 km s⁻¹, ±0.15 km s⁻¹, and ±0.5 per cent, for the left-hand, middle, and right-hand column, respectively.

of SH with $ℓ_{max}$ = 5). To bring further evidence to this suggestion, we compared the large-scale SDO Dopplergrams to the DI maps, this time at the stellar phase which maximizes the absolute value of the RV in each chunk (see Fig. 6). This time, the median absolute value of ρ rises to 0.63 ± 0.1, with a maximum of 0.83 (see Fig. A4). We also note that positive RVs (respectively negative RVs) lead to positive (respectively negative) correlation coefficients between large-scale SDO Dopplergrams and DI maps. Since the input line profiles of each chunk are centred (i.e. the RVs are mean subtracted), negative RVs correspond to negative correlations between the DI maps and SDO Dopplergrams.

The fact that ρ does not reach values greater than 0.83 and varies significantly from one rotation to the next, even at the phase of maximum distortion, can be explained by several factors. First, the accuracy of the DI maps is intrinsically limited by the fact that the mean absolute value of the maximum RV of each chunk is only 1.8 m s⁻¹. Secondly, our assumed stellar inclination of 80° will slightly affect the latitude of the active regions. In addition, choosing the maximum entropy maps will likely favour maps with more equatorial active regions. These two factors can induce mismatches with the solar surface, especially at the beginning of the magnetic cycle, where active regions appear at higher latitudes. Finally, the presence of more complex signals on the SDO Dopplergrams (e.g. supergranulation), not included in our DI framework, could also affect the comparison with the brightness maps.

We also find a significantly better match between the large-scale SDO intensitygrams and the DI maps at the phase of maximum profile distortion. The absolute value of the correlation coefficient between these maps increases from 0.11 ± 0.02 to 0.32 ± 0.10, with a maximum of 0.7. Since the activity-induced distortions in the line profiles are dominated by the inhibition of convective blueshift in faculae (e.g. Meunier et al. 2010, 2015, 2022; Milbourne et al. 2019), we naturally expect the DI maps to be better matched by the
<figure><img src="image_9.png" /><figcaption>Figure 7. Average Pearson correlation coefficient between the input brightness distribution and the recovered DI maps as a function of the noise level in the continuum-normalized line profiles (left panel), and the number of epochs (right panels). In both panels, the dashed lines indicate the best-fitting straight line (note that, in the right panel, the first point is not included in the linear fit).</figcaption></figure>

SDO Dopplergrams than the intensitygrams, which are primarily sensitive to sunspots. At the stellar phase of maximum profile distortion, the locations of the most prominent regions of convective blueshift inhibition are relatively well recovered in our brightness maps. Since faculae tend to cluster around sunspots, the brightness maps will naturally roughly recover the position of sunspots, hence the increased correlation with the large-scale SDO intensitygrams. Conversely, we do not observe good matches between the DI maps and the large-scale SDO magnetograms (or their absolute values), probably because of the higher level of complexity of the small-scale magnetic field network.

## 4.3 Robustness tests

As shown in the previous section, DI is able, in most cases, to accurately reproduce large-scale SDO Dopplergrams from densely sampled HARPS-N solar spectra. This naturally raises the question of its application to other stars. In this section, we perform simulations to estimate the accuracy of DI maps as a function of the number of points, noise level, and RV value. From one representative DI maps (chunk 7 in this case), we generate a time series of synthetic profiles, to which we add normally distributed noise. The synthetic profiles are then modelled with DI, and the best-fitting map is compared to the input brightness distribution at different stellar phases (i.e. RV values). In each case, we select a number $N_{pt}$ of epochs (between 5 and 20) and randomly assign their stellar phases between 0 and 1.3 Relative noise levels between $10^{-5}$ and $3 \times 10^{-4}$ are considered, the latter being of similar amplitude to the largest distortions in the input synthetic profiles. We repeat the simulation 10 times for each number of epochs and noise level, each time randomly assigning different observational phases and noise realization to the synthetic profiles.

Fig. 7 shows the evolution of the Pearson correlation coefficient between the input and recovered brightness maps, as a function of the CCF noise per observation and number of observations. We find that the correlation coefficient between the input and the recovered brightness distributions linearly increases with the number of epochs

for $N_{pt} \gtrsim 7$. In most cases, a smaller number of observations per rotational phase does not ensure a sufficient phase coverage to accurately recover the input brightness distribution. The correlation coefficient decreases roughly linearly with the noise level, with values consistently larger than 0.5 for relative levels larger or equal to ~$10^{-4}$ with respect to the continuum, corresponding to the typical amplitude of the activity-induced distortions in the input data. It is worth noting that, for more than seven observations, it is the CCF noise rather than the number of points that controls the accuracy of the retrieved brightness map. As expected, we also find that the input map is more accurately reproduced at phases with larger RV values. A two-dimensional colour map of the correlation between the input and retrieved brightness maps as a function of the number of observations and CCF noise is provided in Fig. A5.

## 4.4 Impact of stellar parameters

### 4.4.1 Differential rotation

DI is a powerful tool to constrain the latitudinal differential rotation (DR) of rapidly rotating stars (e.g. Donati et al. 2000; Petit, Donati & Collier Cameron 2002; Yu et al. 2019; Zaire et al. 2021). Assuming a solar-like DR, the rotation rate $\Omega$ as a function of the colatitude $\theta$ is given by

$$ \Omega(\theta) = \Omega_{eq} - (\cos \theta)^2 \, d\Omega, \qquad (3) $$

where $\Omega_{eq}$ is the rotation rate at the stellar equator and $d\Omega$ is the difference in rotation rate between the equator and the pole. Traditionally, DI inversions are performed for a grid of ($\Omega_{eq}$, $d\Omega$) to a fixed level of information (i.e. to a fixed level of $f_{DI}$), and the best-fitting DR parameters are estimated from the resulting $\chi^2$ map. This method will most likely fail in the case of the Sun, since active regions evolve on time-scales similar to the rotation period. Yet, our input data set covers several consecutive solar rotation cycles. Thus, by studying how the position of the largest features evolves between consecutive maps, one may still be able to place constraints on the solar rotation and DR.

The optimal brightness distribution of each chunk is projected into an rectangular view, and normalized by its standard deviation. For each pair of consecutive rectangular projections, we differentially rotate the rows of the first map for a grid of ($\Omega_{eq}$, $d\Omega$) values using equation (3), and compute the correlation between the rotated map and the second map. The inferred DR parameters are obtained by modelling the correlation map with a second-order polynomial near the maximum correlation. Uncertainties on the DR parameters are obtained independently for each pair of consecutive chunks using a bootstrap. This is done by adding normally distributed Gaussian white noise to the best-fitting line profiles of Section 4.1, and, then, inverting these line profiles into brightness distributions with DI. The best-fitting DR parameters are then computed between the pair of synthetic brightness maps using the technique described above. We repeat this process 200 times with different white noise realizations, and take the standard deviation on the distribution of $\Omega_{eq}$ and $d\Omega$ as $1\sigma$ uncertainties for the parameters. For each chunk, we consider that the Sun's rotation has been detected if the best-fitting value of $\Omega_{eq}$ lies between 0.224 and 0.256 rad $d^{-1}$ (i.e. rotation periods of 28 and 24.5 d, respectively). We could not draw a direct link between the variation of $\Omega_{eq}$ and the Sun's magnetic cycle, which is not surprising as the time-scale of our observations is only 2 yr.

The top panel of Fig. 8 shows the average correlation map for all pairs of chunks with detected solar rotation, in the ($\Omega_{eq}$, $d\Omega$) space. The Sun's equatorial rotation period is accurately recovered

<footnote n="3">³Note that our assumption to limit our observations to one stellar rotation phase is purely driven by the fact that, for the Sun, active regions evolve on the time-scale of one rotation period. For other stars, chunks should be defined on the time-scale on which stellar activity varies rather than on one rotation phase.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Correlation between consecutive DI brightness maps averaged over all pairs of chunks with detected solar rotation (top panel) and between chunks 8 and 9 in Table 1 (bottom panel), in the ($\Omega_{eq}$, $d\Omega$) space. In both panels the white dashed lines indicate the solar DR values adopted in the study.</figcaption></figure>

for 14 of the 20 pairs of brightness maps analysed, with an average equatorial rotation rate of $0.227 \pm 0.018 \text{ rad d}^{-1}$ (i.e. a rotation period of $27.7 \pm 2.1 \text{ d}$ at the equator). This value is fully consistent with the solar rotation period measured by Klein et al. (2024) from the same data set with GPs ($27.4 \pm 0.4 \text{ d}$). In contrast, the solar rotation rate is not well retrieved when at least one of the chunks exhibits a low activity level (i.e. a peak-to-peak RV variation smaller than ~3 m s$^{-1}$, like in chunks 2 and 19). We also report a marginal $2.5\sigma$ detection of DR between chunks 8 and 9, with $d\Omega = 0.05 \pm 0.02 \text{ rad d}^{-1}$, consistent with the solar value of ~0.07 rad d$^{-1}$ (see the bottom panel of Fig. 8). However, the Sun's DR remains undetected in all other pairs of chunks, which is most likely due to the poor latitudinal resolution of our brightness maps and the fast-evolving solar activity.

### 4.4.2 Stellar inclination

It is challenging to constrain the inclination of the star's rotation axis from the activity-induced distortions of unpolarized line profiles, especially for slowly rotating stars. We found that, at the level of one chunk, near-pole-on configurations (i.e. inclinations smaller than ~30°) can be confidently excluded, but that no precise estimate of the star's inclination could be recovered, given the relatively small length of the chunks and angular resolution of the recovered brightness maps.

In an attempt to refine the constraints on the star's inclination, we performed a DI inversion of all the 439 line profiles simultaneously, whilst varying the assumed stellar inclination in the fit. The resulting grid of reduced $\chi^2$ reaches a minimum of 2.2, significantly larger than 1, due to the fact that activity evolution is not accounted for in the model. However, we find a minimum in $\chi^2$ grid, yielding a stellar inclination estimate of $69^\circ \pm 16^\circ$. Inclinations smaller than ~40° are rejected at a $3\sigma$ level. We conclude that DI can help setting lower limits on stellar inclinations, but cannot provide precise

<table><caption>**Table 2.** Orbital periods $P_{\text{orb}}$, RV semi-amplitudes $K_p$, masses $M_p$, and orbital phases $\phi_p$ of the synthetic planets considered in this study.</caption><thead><tr><th>Case</th><th><i>P</i><sub>orb</sub><br>(d)</th><th><i>K</i><sub>p</sub><br>(m s<sup>-1</sup>)</th><th><i>M</i><sub>p</sub><br>(M<sub>⊕</sub>)</th><th><i>φ</i><sub>p</sub></th></tr></thead><tbody><tr><td>1</td><td>6.0</td><td>1.0</td><td>2.8</td><td>{0.0, 0.33, 0.7}</td></tr><tr><td>2</td><td>6.0</td><td>0.4</td><td>1.1</td><td>{0.0, 0.33, 0.7}</td></tr><tr><td>3</td><td>100.0</td><td>0.5</td><td>3.7</td><td>{0.0, 0.33, 0.7}</td></tr></tbody></table>

constraints on this parameter for near-equator-on stars. In particular, our assumption of an inclination of 80° in the DI process does not have a major impact on the results of this paper.

## 5 PLANET INJECTION AND RECOVERY

We perform a series of planet injection-recovery tests in order to estimate the sensitivity of the DI framework to low-mass planet signatures. We consider the three different planet signatures listed in Table 2. Two short-period planets with an orbital period $P_{\text{orb}}$ of 6 d and RV semi-amplitudes of 1 m s$^{-1}$ (our fiducial case) and 0.4 m s$^{-1}$ (Earth-mass planet), and one super-Earth at an orbital period of 100 d. In each case, we linearly interpolate the observed CCFs and shift them according to the RV signature of a single of the three planets. The planet's orbit is assumed circular, and, thus, the RV signature $v_p$ as a function of time $t$ is given by

$$v_p = K_p \sin 2\pi \left[ \frac{T_0 - t}{P_{\text{orb}}} + \phi_p \right], \quad (4)$$

where $K_p$, $P_{\text{orb}}$, and $\phi_p$ are the planet RV semi-amplitude, orbital period, and orbital phase, respectively, and where $T_0$ is the reference time of equation (2). For each planet, we consider three different orbital phases to ensure that the results are not biased. Note that we also extract the CCF RVs, $v_{\text{obs}}$, using a Gaussian fit. Multiplanetary systems are not considered in this proof-of-concept study.

### 5.1 Doppler imaging as a planet detection tool

We first perform a simple test to assess the ability of DI to be used as an activity filtering tool. As in Section 4.1, we blindly model the CCFs using DI, regardless of the presence of planet signatures in the data. From the best-fitting brightness distributions, we generate synthetic CCFs at the epochs of observations, compute the RVs of the synthetic CCFs, and subtract these RVs from $v_{\text{obs}}$. Ideally, the DI correction should remove most of the stellar activity signals whilst preserving the injected planet signature.

The GLS periodograms of the DI-corrected RV time series are shown in Fig. 9. Similarly to Section 4.1, most of the power at $P_{\text{rot}}$ and its harmonics has disappeared in the DI-corrected RVs, which indicates that DI has efficiently filtered the quasi-periodic activity signals from the data. We also note that the most prominent peak in the periodogram of the DI-corrected RVs corresponds, in all cases, to the orbital period of the injected planet. We then model the DI-corrected RVs using equation (4), fitting for $K_p$, $P_{\text{orb}}$, and $\phi_p$, as well as for an uncorrelated jitter term, quadratically added to the formal RV uncertainties. The parameter space is sampled using the Bayesian Markov chain Monte Carlo (MCMC) process implemented in the open-source software PYANETI (Barragán, Gandolfi & Antoni-ciello 2019; Barragán et al. 2022a). Non-informative uniform prior distributions are adopted for all the model parameters.

The optimal planet parameters are shown in Table 3. As intuited from Fig. 9, the orbital period and phase of the injected planets
<figure><img src="image_11.png" /><figcaption>**Figure 9.** GLS periodograms of the HARPS-N solar RVs after the injection of the planet RV signatures listed in Table 2. In each case, the top black and bottom yellow lines represent the periodogram of the RV time series before and after applying DI to filter stellar activity signals. In each panel, the blue vertical dashed line indicates the orbital period of the injected planet, and the three red dotted lines mark the Sun's rotation period and first two harmonics.</figcaption></figure>

<table><caption>**Table 3.** Best estimates of the orbital period, RV semi-amplitude, and orbital phase of the planet signatures described in Table 2, when retrieved from the DI-corrected RVs.</caption><thead><tr><th>Case</th><th><i>P</i><sub>orb</sub><br>(d)</th><th><i>K</i><sub>p</sub><br>(m s<sup>-1</sup>)</th><th><i>φ</i><sub>p</sub></th></tr></thead><tbody><tr><td rowspan="3">Planet 1</td><td>6.002 ± 0.002</td><td>0.58 ± 0.04</td><td>0.02 ± 0.02</td></tr><tr><td>6.000 ± 0.002</td><td>0.57<sup>+0.04</sup><sub>-0.05</sub></td><td>0.33 ± 0.02</td></tr><tr><td>6.000 ± 0.002</td><td>0.62 ± 0.04</td><td>0.68 ± 0.02</td></tr><tr><td rowspan="3">Planet 2</td><td>6.006 ± 0.005</td><td>0.23<sup>+0.04</sup><sub>-0.05</sub></td><td>0.04<sup>+0.07</sup><sub>-0.06</sub></td></tr><tr><td>5.997 ± 0.005</td><td>0.24<sup>+0.05</sup><sub>-0.06</sub></td><td>0.31<sup>+0.06</sup><sub>-0.04</sub></td></tr><tr><td>5.998 ± 0.004</td><td>0.29<sup>+0.04</sup><sub>-0.05</sub></td><td>0.66 ± 0.05</td></tr><tr><td rowspan="3">Planet 3</td><td>99.5<sup>+0.8</sup><sub>-1.0</sub></td><td>0.35<sup>+0.05</sup><sub>-0.06</sub></td><td>-0.01 ± 0.04</td></tr><tr><td>100.2<sup>+1.0</sup><sub>-1.2</sub></td><td>0.29<sup>+0.04</sup><sub>-0.05</sub></td><td>0.33 ± 0.04</td></tr><tr><td>99.5<sup>+0.7</sup><sub>-0.8</sub></td><td>0.29<sup>+0.04</sup><sub>-0.05</sub></td><td>0.68<sup>+0.04</sup><sub>-0.03</sub></td></tr></tbody></table>

are fully consistent with their injected counterparts. However, the recovered RV semi-amplitudes are systematically underestimated by about 40 per cent, which suggests that the DI process has partly filtered out the planet signatures. In other words, planet signatures in the CCF can be, to an extent, modelled with DI alone. The planet CCF signature is a simple Doppler shift and, as such, affects only the first derivative of the CCF with respect to the wavelength. On the other hand, stellar activity affects mostly higher order derivatives, which is the basis of data-driven methods to filter stellar activity in the wavelength space (e.g. Jones et al. 2022; Collier Cameron et al. 2021; Klein et al. 2024). However, there is no reason for the first derivative of the CCF to be orthogonal to higher odd-order derivatives (e.g. due to cross-talks between these derivates and the first-order derivative). This cross-talk could explain why our activity model is able to filter out part of the planet signature. Therefore, similarly to other wavelength-based activity modelling frameworks (e.g. John et al. 2022, 2023; Wilson et al. 2022), we need to simultaneously model stellar activity signals and planet signatures. However, it is

worth noting that the blind application of DI was able to provide accurate estimates for the orbital period and phase of the planets, which can be used as informative priors to other models.

### 5.1.1 Planet detection limit

Building on the results obtained above, we now assess the ability of DI to be used as a blind planet search tool by generating 1000 data sets from the solar HARPS-N observations. For each data set, we inject a single planetary signal directly into the same 439 CCFs selected in Section 3.2 (corresponding to the 23 chunks listed in Table 1). The planet signatures are computed using equation (4), assuming circular planetary orbits. The planet orbital period, $P_{\text{orb}}$, is randomly drawn from a log-uniform law between 5 and 300 d, whereas $K_p$ and $\phi_p$ are randomly drawn from uniform laws between 0.2 and 1.0 m s$^{-1}$ and 0.0 and 1.0, respectively (baseline of 2.1 yr).

For each data set, we use three independent methods to filter stellar activity signals. First, we simply model the stellar activity RV signal by a sine wave at the rotation period of the Sun and its first harmonic. Secondly, we apply the Doppler-constrained principal component analysis (DCPCA) framework introduced in Jones et al. (2022), as implemented in Klein et al. (2024). This framework, conceptually similar to the SCALPELS algorithm (e.g. Collier Cameron et al. 2021; John et al. 2022), has been shown to robustly filter long-term stellar activity variations whilst preserving planet RV signatures in most cases. Finally, we blindly apply DI to the input CCFs and subtract the RVs derived from the best-fitting line profiles from the input RVs. For each of the three cases, we compute a GLS periodogram from the activity-filtered RVs and estimate the period of the most prominent peak and its uncertainty using $\chi^2$ statistics. In each case, the planet is said to be detected if (1) the injected orbital period lies within $3\sigma$ of the recovered one, and (2) the false alarm probability of the recovered period, computed using Baluev (2008)'s criterion, is smaller than 0.1 per cent. It is worth noting that, without activity filtering, none of the injected planets would be retrieved with this criterion.

The completeness maps of the planet injection-recovery tests for the three activity-filtering techniques are shown in Fig. 10. The sine wave is too simplistic and only filters part of the stellar activity signal. Therefore, if planets with RV semi-amplitudes larger than ~0.7 m s$^{-1}$ are relatively well recovered, the completeness decreases quickly with decreasing RV semi-amplitude. The DCPCA framework is significantly more sensitive to smaller amplitude planet signatures, with about 87 per cent of completeness for planet semi-amplitudes larger than 0.4 m s$^{-1}$. The completeness completely drops for smaller semi-amplitudes, which we attribute to the fact that DCPCA only provides a partial filtering of quasi-periodic stellar activity signals (see fig. 7 of Klein et al. 2024). DI performs exceptionally well with long-period planets (i.e. $P_{\text{orb}} > 100$ d), with the totality of injected planets with semi-amplitude larger than 0.4 m s$^{-1}$ and a completeness of 67 per cent for semi-amplitudes smaller than 0.4 m s$^{-1}$. On the other hand, almost all planets with orbital periods near the Sun's rotation period and first harmonics are lost in the process, which is not surprising as DI is driven by the modulation of CCF distortions at the star's rotation period. In the DI case, we also note that most of the discarded signals with $P_{\text{orb}} \lesssim 10$ d and $K_p \gtrsim 0.4$ m s$^{-1}$ have an orbital period between 8.5 and 9.5 d, which corresponds to the third harmonic of the solar rotation period (see also Fig. 4). None the less, this test demonstrates that DI can, in principle, be used as a planet search tool, and complements well other methods in the planet parameter space to which it is sensitive.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Completeness estimates for the recovery of the 1000 planet RV signatures injected in the HARPS-N solar CCFs for different activity filtering techniques (see Section 5.1.1). The number and the colour of each cell both indicate the completeness of the planet injection recovery (in per cent).</figcaption></figure>

## 5.2 Joint activity and planet fit

The previous section has demonstrated the need to incorporate the planet search directly into the DI framework. This framework is first described in Section 5.2.1, and then applied to recover the three planetary signals listed in Table 2.

### 5.2.1 Method

We use the method introduced by Petit et al. (2015), and validated on various cases (e.g. Donati et al. 2017; Yu et al. 2017; Heitzmann et al. 2021; Klein et al. 2021b, 2022). We start from a grid of parameters ($K_p$, $P_{orb}$, and $\phi_p$ in this case). For each set of parameters in the grid, we compute the planet RV signature using equation (4), and use this signature to shift the wavelength of each observed CCF. We then apply DI to the shifted CCFs to a given level of entropy, which gives a value $\chi^2$ quantifying the goodness of the fit. We repeat this process for all the grid parameters and estimate the best-fitting parameters and uncertainties from the resulting $\chi^2$ distribution (Press et al. 1992). The planet search is performed independently on each of the 23 chunks of Table 1. The best estimate of the planet parameters is simply obtained by averaging the best-fitting parameters over all chunks and propagating uncertainties.

For comparison, we use one- and two-dimensional GPs to model stellar activity signals while fitting for the planet RV signature. We use the GP framework of Rajpaul et al. (2015) as implemented in the PYANETI software of Barragán et al. (2022a). This framework uses a quasi-periodic GP to represent the activity signal, with covariance kernel between epochs $t_i$ and $t_j$ given by

$$k(t_i, t_j) = A^2 \exp \left[ - \frac{(t_j - t_i)^2}{2\lambda_e^2} - \frac{\sin^2 \pi(t_j - t_i)/P_{GP}}{2\lambda_p^2} \right], \quad (5)$$

where A is the GP amplitude, and where the GP period $P_{GP}$, evolution time-scale $\lambda_e$, and inverse harmonic complexity $\lambda_p$ are the three hyperparameters of our model. In the one-dimensional case, we model the RVs only, and our model is composed of a planet component computed with equation (4) (i.e. three free parameters, $K_p$, $P_{orb}$, and $\phi_p$), a stellar activity signal modelled by a quasi-periodic GP (following equation 5), and a white noise term, defined as the quadratic sum of the photon noise and an additional uncorrelated

jitter term $\sigma_j$, to absorb variations not accounted for by the GP. In the two-dimensional case, the RVs are modelled jointly with the CCF FWHM, known to be a reliable proxy of the solar activity RV signals (Klein et al. 2024). As described in Rajpaul et al. (2015), the RV and FWHM activity signals are modelled as linear combinations of a latent variable (typically the square of the photometric flux) and its first temporal derivative, following the FF' framework of Aigrain, Pont & Zucker (2012). This latent variable is modelled as a GP with the quasi-periodic covariance kernel of equation (5).

### 5.2.2 Results

We proceed to recover the three planetary signals listed in Table 2. Recovering the planet RV semi-amplitude, orbital period, and orbital phase is challenging, even for multidimensional GPs. However, as shown in Section 5.1, the blind CCF modelling with DI has enabled us to accurately estimate the planet's orbital period and phase. We thus adopt Gaussian priors for these two parameters, using the best-fitting values and uncertainties reported in Table 3. For DI, the grids on $P_{orb}$ and $\phi_p$ are centred on the mean prior values and extend up to 10 standard deviations. A uniform prior between 0 and 5 m s⁻¹ is adopted for $K_p$. Non-informative priors are adopted for the other parameters of the model.

The planet RV semi-amplitudes recovered by our different models are given in Table 4 and shown in Fig. 11. We find that the recovered planet RV semi-amplitudes lie systematically within ~2σ of the injected value, for each injected planet and activity modelling technique. Additionally, we do not identify any systematic trend between the recovered semi-amplitude and the injected planet orbital phase. As already shown in Klein et al. (2024), the one-dimensional GP is only able to provide a marginal ~1.8σ detection of the 100-d period planet, due to long-term activity evolution. The two-dimensional GP provides significantly more precise semi-amplitude estimates for Planet 3. However, as we change the planet orbital phase, we do observe significant variations of ±2σ of the recovered value with respect to the injected one. In contrast, the semi-amplitudes estimates from DI are systematically consistent with the injected values at 1σ, with error bars comparable to those obtained with the two-dimensional GP.
<table><caption>**Table 4.** Best-fitting estimates of the RV semi-amplitude of the planets described in Table 2 and recovered using one- and two-dimensional GPs and DI. All semi-amplitudes are given in m s⁻¹ with 1σ uncertainties.</caption><thead><tr><th>Orbital phase</th><th>Planet 1 (P<sub>orb</sub> = 6 d, K<sub>p</sub> = 1 m s<sup>-1</sup>)</th><th>Planet 2 (P<sub>orb</sub> = 6 d, K<sub>p</sub> = 0.4 m s<sup>-1</sup>)</th><th>Planet 3 (P<sub>orb</sub> = 100 d, K<sub>p</sub> = 0.5 m s<sup>-1</sup>)</th></tr></thead><tbody><tr><td></td><td>0.0<br>0.33<br>0.7</td><td>0.0<br>0.33<br>0.7</td><td>0.0<br>0.33<br>0.7</td></tr><tr><td>1D GP</td><td>0.95<sup>+0.05</sup><sub>-0.07</sub></td><td>0.37 ± 0.05</td><td>0.46 ± 0.25</td></tr><tr><td>2D GP</td><td>0.96<sup>+0.08</sup><sub>-0.07</sub></td><td>0.95<sup>+0.06</sup><sub>-0.07</sub></td><td>0.43<sup>+0.26</sup><sub>-0.25</sub></td></tr><tr><td>DI</td><td>1.02 ± 0.05</td><td>0.94 ± 0.06</td><td>0.54 ± 0.06</td></tr><tr><td></td><td></td><td>1.05 ± 0.06</td><td>0.36<sup>+0.27</sup><sub>-0.23</sub></td></tr><tr><td></td><td></td><td>0.95 ± 0.05</td><td>0.47 ± 0.05</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Best estimates of the planet RV semi-amplitude for the three different planets listed in Table 2 and retrieved using three different stellar activity models. In each panel, the injected RV semi-amplitude is indicated by the horizontal dashed line.</figcaption></figure>

# 6 DISCUSSION AND CONCLUSIONS

In this proof-of-concept study, we assessed the ability of traditional DI to model the magnetic activity and search for low-mass planets around of Sun-like stars. We divided 2 yr of Sun-as-a-star observations, collected with the high-precision spectrograph HARPS-N between 2022 and 2024, into 27-d long chunks of data, and modelled the time series of CCFs within each chunk with DI. The output relative brightness topologies, shown in orthographic projection in Fig. 2 and in rectangular view in Fig. A1, exhibit typical filling factor of 0.3 per cent, and angular resolution of about 36°, and vary significantly from one rotation cycle to the next, as expected from the fast-evolving solar activity. We found that the RVs derived from the best-fitting line profiles match well the HARPS-N solar RVs, with residuals of 0.58 m s⁻¹ rms (see Fig. 3), similar to the level obtained with GP modelling (Klein et al. 2024), and no sign of modulation at the solar rotation period and harmonics.

## 6.1 Accuracy of the brightness maps

We found that the brightness maps obtained with DI match well the large-scale structures of contemporaneous SDO Dopplergrams (when modelled with SH of degree $ℓ_{max}$ = 5), provided that the absolute value of the RV at the rotational phase of comparison is larger than ~2 m s⁻¹. Smaller RVs reflect CCF distortions that are barely larger than the noise and are generally less well recovered with DI. The DI brightness maps match significantly less well the SDO/HMI intensitygrams (probably because the main cause of CCF distortions is the inhibition of convective blueshift rather than brightness contrast effects), and the magnetograms, which exhibit structures with a higher degree of complexity.

The match between Doppler images and SDO Dopplergrams could bring constraints on the location of stellar active regions for the characterization of exoplanet atmospheres with space-based
transmission spectroscopy. Stellar activity is currently one of the main limitations in the characterization of the atmosphere of rocky planets with the James Webb Space Telescope (e.g. Rackham, Apai & Giampapa 2018, 2019; Fournier-Tondreau et al. 2024; Radica et al. 2025). Knowing the position of the largest active regions beforehand would help refine stellar activity models in exoplanet atmosphere retrievals. For the time being, we have shown that DI is able to pick up spatially resolved information that correlates moderately well with SDO/HMI Dopplergrams. However, more work is needed to understand what physical quantities are extracted with DI and how to convert DI maps into physical constraints for planet atmosphere retrievals. Realistic simulated spectra (e.g. Sowmya et al. 2023; Smitha et al. 2025) or disc-resolved high-resolution solar spectra, like those that the PoET (Paranal solar ESPRESSO Telescope) instrument will collect (Santos et al. 2025), will help improve the accuracy of the DI modelling.

## 6.2 Towards an operational version of Doppler imaging for Sun-like stars

As shown in Section 5, DI is a reliable algorithm for searching for long-period low-mass planet signatures, with accurate retrieval of the planet's orbital period and phase. However, blindly applying DI to a data set without jointly fitting for planet signature will result in largely underestimated planet masses. In contrast, when the DI activity modelling is jointly performed with a planet fit, the planet mass estimates are as accurate as those obtained with multidimensional GP regression, which is arguably the current state-of-the-art activity filtering algorithm. These promising results suggest that DI could become a reliable alternative to traditional RV-based activity modelling techniques in the search for long-period low-mass planets, provided that a number of limitations, detailed below, are addressed.

Performing the DI inference and planet search in a Bayesian framework is a natural next step, motivated by the results of this study. However, this requires the DI inverse process to be made significantly faster and more scalable than most of the DI codes used in the literature, including the one used in this study. Recently, Luger et al. (2021a) made significant progress in this direction by developing a linear, scalable, and open-source version of a DI code.<fnref n="4" /> Though the linearization of the DI process may make it difficult to incorporate some physical processes like DR or convective blueshift inhibition, this code is currently our best hope for a Bayesian mapping of stellar surfaces from high-resolution spectra. In fact, information field theory, which describes a physical field using Bayesian statistics (see Enßlin, Frommert & Kitaura 2009; Enßlin & Frommert 2011), provides a flexible framework fully adapted to DI, bypassing the constraints of regularization (e.g. maximum entropy, Tikhonov).

Another fundamental limitation of current DI algorithms is that they generally do not account for the evolution of stellar activity with time. This is a problem for the search for long-period planets, which requires DI to be performed on long baselines, and for stars with fast-evolving activity, like the Sun. Describing the temporal evolution of the coefficients of the stellar surface map with GPs is one of the most promising ways of addressing these issues. This, in fact, has been recently incorporated in light-curve modelling (Luger, Foreman-Mackey & Hedges 2021b) and for ZDI (Finociety &

Donati 2022). However, this implementation carries the risk of adding further degeneracies or assumptions to the inversion process (e.g. accounting for the different evolution time-scales of active regions of different sizes, latitudes, and depths). For Sun-like stars in particular, the rapidly evolving activity might be difficult to capture with GPs.

Finally, incorporating more physically accurate descriptions of the stellar surfaces is one of the most exciting long-term perspectives of DI, taking advantage of the unrivalled precision of new-generation spectrographs. Inverting the Sun's CCF distortions into a velocity field rather than into a relative brightness distribution will likely improve the fit quality, especially in the line cores where the largest deviations are observed (see Figs A2 and A3). This should also arguably improve the match with SDO Dopplergrams. To go further, modelling the full spectra of Sun-like stars, or, in the first instance, CCFs computed using tailored line masks, will allow us to separate different spectral contributions of stellar activity. For example, whereas contrast-induced distortions of stellar line profiles, dominated by spots, have a well-known chromatic dependence, the inhibition of the convective blueshift in faculae depends on the line formation temperature (Al Moulla et al. 2022; Rescigno & Al Moulla 2025). These spectral differences could, in principle, enable one to produce an intensitygram and a Dopplergram from the same time series of spectra, provided that the SNR of the spectra is high enough. Similarly, new-generation near-infrared (NIR) spectra of M dwarfs and pre-main-sequence stars are primarily sensitive to two effects of stellar activity: the spot-induced contrast and the Zeeman splitting in magnetically sensitive lines (Reiners et al. 2013; Klein et al. 2021b; Donati et al. 2023). Accounting for the different magnetic sensitivities (i.e. Landé factors) of the spectral lines used in the DI inversion could, in principle, enable one to separate the Zeeman-induced distortions from contrast effects, provided that enough spectral lines with accurate Landé factors are available in the NIR. Lastly, other causes of stellar variability like supergranulation (Rincon & Rieutord 2018; Lakeland et al. 2024; O'Sullivan & Aigrain 2024; O'Sullivan et al. 2025) remain poorly understood and might, in fact, represent an intrinsic limitation to the accuracy of DI maps.

## ACKNOWLEDGEMENTS

We thank an anonymous referee for valuable comments that improved the paper. This publication is part of a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Framework Programme (grant agreement no. 865624). KAM acknowledges support from the Swiss National Science Foundation (SNSF) under the Postdoc Mobility grant P500PT_230225. ACC acknowledges support from STFC consolidated grant number ST/V000861/1 and UKRI/ERC Synergy Grant EP/Z000181/1 (REVEAL). AM acknowledges funding from a UKRI Future Leader Fellowship, grant number MR/X033244/1 and a UK Science and Technology Facilities Council (STFC) small grant ST/Y002334/1.

This study is based on observations made with the Italian Teleoscopio Nazionale Galileo (TNG) operated on the island of La Palma by the Fundación Galileo Galilei of the INAF (Istituto Nazionale di Astrofisica) at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. We thank the HARPS-N solar team and TNG staff for processing the solar data and maintaining the solar telescope.

<footnote n="4">⁴Available through the STARRY PYTHON package: https://starry.readthedocs.io/en/latest/notebooks/DopplerImaging_Intro/.</footnote>
DATA AVAILABILITY

The HARPS-N solar data used in this study can be accessed from
the Data Analysis Center for Exoplanet (DACE), either through the
web platform: //dace.unige.ch/sun/ or through the DACE PYTHON
API. The data will be extensively described in Dumusque et al. (in
preparation). The SDO/HMI images are publicly available at: //sdo.
gsfc.nasa.gov/data/.

REFERENCES

Afram N., Berdyugina S. V., 2015, A&A, 576, A34

Aigrain S., Foreman-Mackey D., 2023, ARA&A, 61, 329

Aigrain S., Pont F., Zucker S., 2012, MNRAS, 419, 3147

Al Moulla K., Dumusque X., Cretignier M., Zhao Y., Valenti J. A., 2022, A&A, 664, A34

Al Moulla K., Dumusque X., Figueira P., Lo Curto G., Santos N. C., Wildi F., 2023, A&A, 669, A39

Al Moulla K., Dumusque X., Cretignier M., 2024, A&A, 683, A106

Asensio Ramos A., Díaz Baso C. J., Kochukhov O., 2022, A&A, 658, A162

Astropy Collaboration, 2013, A&A, 558, A33

Astropy Collaboration, 2018, AJ, 156, 123

Astropy Collaboration, 2022, ApJ, 935, 167

Bagnulo S., Szeifert T., Wade G. A., Landstreet J. D., Mathys G., 2002, A&A, 389, 191

Baluev R. V., 2008, MNRAS, 385, 1279

Barragán O., Gandolfi D., Antoniciello G., 2019, MNRAS, 482, 1017

Barragán O., Aigrain S., Rajpaul V. M., Zicher N., 2022a, MNRAS, 509, 866

Barragán O. et al., 2022b, MNRAS, 514, 1606

Barragán O. et al., 2023, MNRAS, 522, 3458

Bellotti S., Petit P., Morin J., Hussain G. A. J., Folsom C. P., Carmona A., Delfosse X., Moutou C., 2022, A&A, 657, A107

Bellotti S. et al., 2025, A&A, 693, A269

Berdyugina S. V., 2005, Living Rev. Sol. Phys., 2, 8

Blackman R. T. et al., 2020, AJ, 159, 238

Blunt S. et al., 2023, AJ, 166, 62

Boro Saikia S. et al., 2016, A&A, 594, A29

Boro Saikia S. et al., 2018, A&A, 620, L11

Brewer J. M. et al., 2020, AJ, 160, 67

Brown S. F., Donati J. F., Rees D. E., Semel M., 1991, A&A, 250, 463

Brun A. S., Browning M. K., 2017, Living Rev. Sol. Phys., 14, 4

Callingham J. R. et al., 2024, Nat. Astron., 8, 1359

Chaplin W. J., Cegla H. M., Watson C. A., Davies G. R., Ball W. H., 2019, AJ, 157, 163

Claret A., Hauschildt P. H., Witte S., 2013, A&A, 552, A16

Collier Cameron A., Unruh Y. C., 1994, MNRAS, 269, 814

Collier Cameron A. et al., 2019, MNRAS, 487, 1082

Collier Cameron A. et al., 2021, MNRAS, 505, 1699

Cosentino R. et al., 2012, in McLean I. S., Ramsay S. K., Takami H., eds, Proc. SPIE Conf. Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV. SPIE, Bellingham, p. 84461V

Crass J. et al., 2021, preprint (arXiv:2107.14291)

Cretignier M., Dumusque X., Allart R., Pepe F., Lovis C., 2020a, A&A, 633, A76

Cretignier M., Francfort J., Dumusque X., Allart R., Pepe F., 2020b, A&A, 640, A42

Cretignier M., Dumusque X., Hara N. C., Pepe F., 2021, A&A, 653, A43

Cretignier M., Dumusque X., Pepe F., 2022, A&A, 659, A68

Cretignier M., Dumusque X., Aigrain S., Pepe F., 2023, A&A, 678, A2

Crossfield I. J. M. et al., 2014, Nature, 505, 654

de Beurs Z. L. et al., 2022, AJ, 164, 49

Delisle J. B., Unger N., Hara N. C., Ségransan D., 2022, A&A, 659, A182

Desort M., Lagrange A. M., Galland F., Udry S., Mayor M., 2007, A&A, 473, 983

Donati J. F., 1999, MNRAS, 302, 457

Donati J. F., Brown S. F., 1997, A&A, 326, 1135

Donati J. F., Landstreet J. D., 2009, ARA&A, 47, 333

Donati J. F., Semel M., Praderie F., 1989, A&A, 225, 467

Donati J. F., Henry G. W., Hall D. S., 1995, A&A, 293, 107

Donati J. F., Semel M., Carter B. D., Rees D. E., Collier Cameron A., 1997, MNRAS, 291, 658

Donati J. F., Mengel M., Carter B. D., Marsden S., Collier Cameron A., Wichmann R., 2000, MNRAS, 316, 699

Donati J. F., Babel J., Harries T. J., Howarth I. D., Petit P., Semel M., 2002, MNRAS, 333, 55

Donati J. F. et al., 2003, MNRAS, 345, 1145

Donati J. F., Howarth I. D., Bouret J. C., Petit P., Catala C., Landstreet J., 2006a, MNRAS, 365, L6

Donati J. F. et al., 2006b, MNRAS, 370, 629

Donati J. F. et al., 2008, MNRAS, 386, 1234

Donati J. F. et al., 2012, MNRAS, 425, 2948

Donati J. F. et al., 2014, MNRAS, 444, 3220

Donati J. F. et al., 2017, MNRAS, 465, 3343

Donati J. F. et al., 2023, MNRAS, 525, 455

Dumusque X., 2018, A&A, 620, A47

Dumusque X., Udry S., Lovis C., Santos N. C., Monteiro M. J. P. F. G., 2011a, A&A, 525, A140

Dumusque X. et al., 2011b, A&A, 535, A55

Dumusque X. et al., 2015, ApJ, 814, L21

Dumusque X. et al., 2021, A&A, 648, A103

Enßlin T. A., Frommert M., 2011, Phys. Rev. D, 83, 105014

Enßlin T. A., Frommert M., Kitaura F. S., 2009, Phys. Rev. D, 80, 105005

Ervin T. et al., 2022, AJ, 163, 272

Fares R. et al., 2009, MNRAS, 398, 1383

Faria J. P. et al., 2022, A&A, 658, A115

Finociety B., Donati J. F., 2022, MNRAS, 516, 5887

Finociety B. et al., 2021, MNRAS, 508, 3427

Fischer D. A. et al., 2016, PASP, 128, 066001

Folsom C. P. et al., 2016, MNRAS, 457, 580

Folsom C. P., Ó Fionnagáin D., Fossati L., Vidotto A. A., Moutou C., Petit P., Dragomir D., Donati J. F., 2020, A&A, 633, A48

Foukal P., 1998, ApJ, 500, 958

Fournier-Tondreau M. et al., 2024, MNRAS, 528, 3354

Gastine T., Morin J., Duarte L., Reiners A., Christensen U. R., Wicht J., 2013, A&A, 549, L5

Goncharskij A. V., Ryabchikova T. A., Stepanov V. V., Khokhlova V. L., Yagola A. G., 1983, Astron. Zh., 60, 83

Gupta A. F. et al., 2021, AJ, 161, 130

Gupta A. F. et al., 2025, AJ, 169, 1

Hall R. D., Thompson S. J., Handley W., Queloz D., 2018, MNRAS, 479, 2968

Hara N. C., Delisle J.-B., 2025, A&A, 696, A141

Hatzes A. P., 1991, MNRAS, 248, 487

Hatzes A. P., 1995, ApJ, 451, 784

Hébrard É. M., Donati J. F., Delfosse X., Morin J., Moutou C., Boisse I., 2016, MNRAS, 461, 1465

Heitzmann A. et al., 2021, MNRAS, 505, 4989

John A. A., Collier Cameron A., Wilson T. G., 2022, MNRAS, 515, 3975

John A. A. et al., 2023, MNRAS, 525, 1687

Jones D. E., Stenning D. C., Ford E. B., Wolpert R. L., Loredo T. J., Gilbertson C., Dumusque X., 2022, Ann. Appl. Stat., 16, 652

Kavanagh R. D., Vidotto A. A., Klein B., Jardine M. M., Donati J.-F., Ó Fionnagáin D., 2021, MNRAS, 504, 1511

Khokhlova V. L., Rice J. B., Wehlau W. H., 1986, ApJ, 307, 768

Kjeldsen H., Bedding T. R., 1995, A&A, 293, 87

Klein B. et al., 2021a, MNRAS, 500, 1844

Klein B. et al., 2021b, MNRAS, 502, 188

Klein B. et al., 2022, MNRAS, 512, 5067

Klein B. et al., 2024, MNRAS, 531, 4238

Kochukhov O., 2016, in Rozelot J.-P., Neiner C., eds, Cartography of the Sun and the Stars. Springer, Cham, Switzerland, p. 177

Kochukhov O., Piskunov N., 2002, A&A, 388, 868

Kochukhov O., Wade G. A., 2010, A&A, 513, A13

Kochukhov O., Drake N. A., Piskunov N., de la Reza R., 2004, A&A, 424, 935
Kochukhov O., Ryabchikova T., Weiss W. W., Landstreet J. D., Lyashko D., 2007, MNRAS, 376, 651

Kochukhov O., Mantere M. J., Hackman T., Ilyin I., 2013, A&A, 550, A84

Lakeland B. S. et al., 2024, MNRAS, 527, 7681

Landi Degl'Innocenti E., Landolfi M., 2004, Polarization in Spectral Lines. Springer, Dordrecht

Landstreet J. D., Mathys G., 2000, A&A, 359, 213

Lavail A., Kochukhov O., Wade G. A., 2018, MNRAS, 479, 4836

Lehmann L. T., Hussain G. A. J., Vidotto A. A., Jardine M. M., Mackay D. H., 2021, MNRAS, 500, 1243

Lehmann L. T. et al., 2024, MNRAS, 527, 4330

Leighton R. B., 1964, ApJ, 140, 1547

Lienhard F., Mortier A., Buchhave L., Collier Cameron A., López-Morales M., Sozzetti A., Watson C. A., Cosentino R., 2022, MNRAS, 513, 5328

Lienhard F., Mortier A., Cegla H. M., Cameron A. C., Klein B., Watson C. A., 2023, MNRAS, 522, 5862

Luger R., Bedell M., Foreman-Mackey D., Crossfield I. J. M., Zhao L. L., Hogg D. W., 2021a, preprint (arXiv:2110.06271)

Luger R., Foreman-Mackey D., Hedges C., 2021b, AJ, 162, 124

Marsden S. C. et al., 2014, MNRAS, 444, 3517

Meunier N., Lagrange A. M., 2019, A&A, 625, L6

Meunier N., Desort M., Lagrange A. M., 2010, A&A, 512, A39

Meunier N., Lagrange A. M., Borgniet S., Rieutord M., 2015, A&A, 583, A118

Meunier N., Mignon L., Lagrange A. M., 2017, A&A, 607, A124

Meunier N., Kretzschmar M., Gravet R., Mignon L., Delfosse X., 2022, A&A, 658, A57

Meunier N., Pous R., Sulis S., Mary D., Lagrange A. M., 2023, A&A, 676, A82

Milbourne T. W. et al., 2019, ApJ, 874, 107

Morin J. et al., 2008, MNRAS, 390, 567

Morin J., Donati J. F., Petit P., Delfosse X., Forveille T., Jardine M. M., 2010, MNRAS, 407, 2269

O'Sullivan N. K., Aigrain S., 2024, MNRAS, 531, 4181

O'Sullivan N. K. et al., 2025, MNRAS, 541, 3942

Pepe F. et al., 2021, A&A, 645, A96

Pesnell W. D., Thompson B. J., Chamberlin P. C., 2012, Sol. Phys., 275, 3

Petersburg R. R. et al., 2020, AJ, 159, 187

Petit P., Donati J. F., Collier Cameron A., 2002, MNRAS, 334, 374

Petit P. et al., 2008, MNRAS, 388, 80

Petit P. et al., 2015, A&A, 584, A84

Petit P. et al., 2021, A&A, 648, A55

Phillips D. F. et al., 2016, in Navarro R., Burge J. H., eds, Proc. SPIE Conf. Vol. 9912, Advances in Optical and Mechanical Technologies for Telescopes and Instrumentation II. SPIE, Bellingham, p. 99126Z

Piskunov N., Kochukhov O., 2002, A&A, 381, 736

Piskunov N. E., Tuominen I., Vilhu O., 1990, A&A, 230, 363

Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical Recipes in FORTRAN: The Art of Scientific Computing, 2nd edn. Cambridge Univ. Press, Cambridge

Queloz D. et al., 2001, A&A, 379, 279

Rackham B. V., Apai D., Giampapa M. S., 2018, ApJ, 853, 122

Rackham B. V., Apai D., Giampapa M. S., 2019, AJ, 157, 96

Radica M. et al., 2025, ApJ, 979, L5

Rajpaul V., Aigrain S., Osborne M. A., Reece S., Roberts S., 2015, MNRAS, 452, 2269

Rauer H. et al., 2014, Exp. Astron., 38, 249

Reiners A., Shulyak D., Anglada-Escudé G., Jeffers S. V., Morin J., Zechmeister M., Kochukhov O., Piskunov N., 2013, A&A, 552, A103

Rescigno F., Al Moulla K., 2025, MNRAS, 536, 3601

Rice J. B., 2002, Astron. Nachr., 323, 220

Rincon F., Rieutord M., 2018, Living Rev. Sol. Phys., 15, 6

Roettenbacher R. M. et al., 2017, ApJ, 849, 120

Saar S. H., Donahue R. A., 1997, ApJ, 485, 319

Santos N. C. et al., 2025, The Messenger, 194, 21

Scherrer P. H. et al., 2012, Sol. Phys., 275, 207

Schou J. et al., 2012, Sol. Phys., 275, 229

Schwab C. et al., 2018, in Evans C. J., Simard L., Takami H., eds, Proc. SPIE Conf. Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII. SPIE, Bellingham, p. 1070271

Semel M., 1989, A&A, 225, 456

Semel M., Donati J. F., Rees D. E., 1993, A&A, 278, 231

Skilling J., Bryan R. K., 1984, MNRAS, 211, 111

Smitha H. N. et al., 2025, ApJ, 978, L13

Sowmya K., Shapiro A. I., Rouppe van der Voort L. H. M., Krivova N. A., Solanki S. K., 2023, ApJ, 956, L10

Strassmeier K. G., 1999, A&A, 347, 225

Strassmeier K. G., 2009, A&AR, 17, 251

Strugarek A., Shkolnik E., 2025, preprint (arXiv:2502.13262)

Thompson S. J. et al., 2016, in Evans C. J., Simard L., Takami H., eds, Proc. SPIE Conf. Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI. SPIE, Bellingham, p. 99086F

Thompson A. P. G., Watson C. A., de Mooij E. J. W., Jess D. B., 2017, MNRAS, 468, L16

Unno W., 1956, PASJ, 8, 108

Vogt S. S., Penrod G. D., 1983, PASP, 95, 565

Vogt S. S., Penrod G. D., Hatzes A. P., 1987, ApJ, 321, 496

Wilson T. G. et al., 2022, MNRAS, 511, 1043

Yu L. et al., 2017, MNRAS, 467, 1342

Yu L. et al., 2019, MNRAS, 489, 5556

Yu H., Aigrain S., Klein B., Cretignier M., Lienhard F., Roberts S. J., 2024, MNRAS, 535, 634

Zaire B., Donati J. F. Klein B., 2021, MNRAS, 504, 1969

Zaire B., Donati J. F., Klein B., 2022, MNRAS, 513, 2893

Zechmeister M., Kürster M., 2009, A&A, 496, 577

Zhao L. L. et al., 2022, AJ, 163, 171

Zhao L. L., Bedell M., Hogg D. W., Luger R., 2024, ApJ, 977, 140

## APPENDIX: DOPPLER IMAGING MODELLING

In this appendix, we give additional details on the DI reconstruction of the surface distribution of the Sun's brightness from HARPS-N spectra. In Fig. A1, we show a flattened polar view of the best-fitting brightness distributions for all the chunks. We also show the best-fitting line profiles in Figs A2 and A3. In Fig. A4, we show the correlation coefficients between SDO large-scale Dopplergrams and the DI relative brightness maps for each chunk at rotation phase 0 and at the phase with maximum absolute RV value. Finally, as a complement to Section 4.3, we show in Fig. A5 the evolution Pearson correlation coefficient between an input brightness distribution and the DI map recovered, as a function of the number of points and noise level.
<figure><img src="image_14.png" /><figcaption>**Figure A1.** Best-fitting relative brightness maps for the 23 chunks listed in Table 1. As in Fig. 2, the colour scale depicts the logarithm of the relative surface brightness (in per cent). Low latitudes are truncated given our assumed inclination of 80°.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure A2.** Time series of all observed (thin black lines) and best-fitting (thick red lines) mean-subtracted line profiles in the first 15 chunks of Table 1. The ±1σ uncertainties on the line continuum is indicated on the right-hand side of each profile. For each chunk, the line profiles are stacked vertically according to the Sun's rotational phase, computed with equation (2), and indicated on the left-hand side of each observation. Note that, in the rare event that two line profiles from the same data chunk are obtained at the same rotational phase (e.g. following a complete solar rotation), only one profile is shown in the figure for clarity (which happened in chunk 8 at phase 0.36, and in chunk 15 at phase 0.35).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure A3.** Same as Fig. A2, but for chunks 16–23. As a reference, we also superimposed all continuum-normalized CCFs in the bottom right panel. To make the figure clearer, one observation (chunk 18 at phase 0.35) is not shown in this figure.</figcaption></figure>
<figure><img src="image_17.png" /></figure>

**Figure A4.** Pearson correlation coefficients between SDO large-scale Dopplergrams ($D_{SDO}$) and DI brightness maps ($I_{DI}$) for each chunk at rotational phase 0 (blue dots) and the phase with maximum absolute RV value (orange stars).

<figure><img src="image_18.png" /></figure>

**Figure A5.** Pearson correlation coefficients between an input brightness distribution and the Doppler maps recovered for different noise levels and number of observations, as described in Section 4.3.

This paper has been typeset from a TeX/Latex file prepared by the author.