© The Authors 2026

1. Introduction

The Sun provides most of the energy to Earth (e.g. Kren et al. 2017), sustaining life and affecting the climate. The radiative flux or radiative power received per unit area at a distance of 1 AU is termed solar irradiance. The spectrally integrated flux is known as the total solar irradiance (TSI), while the wavelength-dependent quantity is the spectral solar irradiance (SSI). Solar irradiance has been measured from space since 1978, revealing variability on all observable timescales, from minutes to decades and longer (e.g., Willson & Hudson 1988; Fröhlich 2006; Kopp 2016, 2025). However, the instrumental record is too short to be compared directly with climate data and to assess the Sun’s impact on the Earth’s climate. Reconstructions are therefore required to extend irradiance records back in time.

Solar irradiance variability on timescales longer than a few hours is driven by changes in solar surface magnetism, specifically the darkening caused by sunspots and the brightening due to faculae and network regions (e.g., Krivova et al. 2003; Yeo et al. 2017; Shapiro et al. 2017). This means that to model irradiance variations, knowledge of the evolving solar surface magnetic field is required. Direct magnetic field measurements have been available only for the past few decades (see, e.g., Pevtsov et al. 2021). Thus, to go further back in time, proxies of solar magnetic activity are needed, particularly those that distinguish between sunspots and faculae.

Sunspot numbers (SNs) provide the longest direct observational record of solar activity, spanning more than 400 years (see reviews by Arlt & Vaquero 2020; Clette et al. 2023). However, direct information on bright features, such as faculae and the network, is missing. Historical full-disc photographs of the Sun in the Ca II K line, available since the late 19th century (Chatzistergos et al. 2022b), are an excellent source of such information (Chatzistergos et al. 2020) as they capture the chromospheric counterparts of faculae. While significant progress has been made in recovering and processing these archives (Chatzistergos et al. 2018, 2022b), further work is required before they can reliably be used for long-term irradiance reconstructions (Chatzistergos et al. 2024). Consequently, for periods lacking direct facular observations, most existing long-term reconstructions have inferred the facular contribution indirectly from sunspot or cosmogenic-isotope data under various, often uncertain, assumptions. This leads to significant differences in the reconstructed long-term trends of solar irradiance ranging from near-zero secular change to several watts per square metre (for a review of various models and resulting uncertainties, see Chatzistergos et al. 2023). Particularly problematic are extended quiet periods, such as grand minima of solar activity, when only a few to no sunspots emerge over years or even decades. Over such periods, no information on the emergence of the magnetic field on the Sun is available, resulting in an underestimated magnetic flux and overestimated secular irradiance variability (for a discussion, see Krivova et al. 2021).

To address this, Krivova et al. (2021) revised the magnetic flux evolution model originally developed by Solanki et al. (2000, 2002b), to incorporate a more realistic representation of small-scale magnetic features such as faculae and network. The updated model relies on modern observations showing that the emergence rates of magnetic regions follow a single power-law distribution across a wide range of fluxes (Thornton & Parnell 2011). This approach gives a more realistic reconstruction of the solar magnetic flux and an improved agreement with independent observations and reconstructions.

Among existing irradiance models (for a review, see, e.g., Chatzistergos et al. 2023), the Spectral And Total Irradiance REconstructions (SATIRE) framework provides a robust, semi-empirical, physics-based approach (Krivova et al. 2011; Solanki et al. 2013). Different SATIRE versions use different observational data as input. SATIRE-S (for the satellite era; Krivova et al. 2003; Yeo et al. 2014; Chatzistergos et al. 2025a) employs full-disc magnetograms and continuum images, providing the highest accuracy but limited temporal coverage (roughly 50 years). SATIRE-T (for the telescopic era; Krivova et al. 2007, 2010; Wu et al. 2018) extends the reconstruction back to the Maunder minimum by using sunspot numbers to model the evolution of the solar surface magnetic field and subsequently the irradiance variability.

In this work, we employ the SATIRE-T model (Krivova et al. 2007, 2010; Wu et al. 2018) together with the revised magnetic flux model of Krivova et al. (2021) to reconstruct solar surface magnetic field evolution and solar irradiance over the entire telescopic era, extending back to the Maunder minimum. The paper is structured as follows: Sect. 2 describes the data used, Sect. 3 outlines the modelling approach, Sect. 4 presents and discusses the results, and Sect. 5 summarises the main findings.

2. Data

2.1. Input data

To describe the emergence of the solar surface magnetic field (see Sect. 3.1), we used two independent sunspot number (SN) records as inputs: the International Sunspot Number version 2.0 (ISNv2; Clette et al. 2023)¹ and the Group Sunspot Number (GSN) series by Chatzistergos et al. (2017, hereafter CEA17)². ISNv2 provides daily values since 1818, monthly data from 1749, and annual values from 1700. Here we combine daily values from 1818 onwards with monthly ones back to 1749, interpolated to a daily resolution. The CEA17 series begins in 1739, but data prior to 1749 are sparse and irregular and are therefore excluded. This record ends in 2010, since raw group counts in the Vaquero et al. (2016) database are not available afterwards. We linearly regressed the CEA17 record against ISNv2 over their 1976–2010 overlap to match their scales and extended CEA17 to the present using ISNv2.

To extend these records back to the beginning of telescopic observations, we also incorporated the GSN series by Hoyt & Schatten (1998, hereafter HoSc98), which has provided daily data since 1610. However, because the earliest observations contain substantial gaps, we used only the records after 1650. We linearly scaled the HoSc98 series to separately match ISNv2 and CEA17, respectively, and used the scaled versions to extend these records from 1650 to 1749. As an alternative extension, we also used the annual sunspot numbers estimated with the active day fraction method by Carrasco et al. (2022, 2024, 2025). This record covers 1635–1721 and 1725–1726, but we excluded the latter period due to discontinuity with the earlier record. We used this series to extend ISNv2 back to 1635. The gap between 1721 and 1749 in this case was filled by interpolating annual ISNv2 values to a daily cadence. Interpolation was also applied to fill all remaining gaps in the series, since the magnetic flux evolution model requires daily sunspot number inputs.

We note that the HoSc98 GSN has been criticised for potential inaccuracies (Clette et al. 2023; Lockwood et al. 2016; Usoskin et al. 2016; Chatzistergos et al. 2025b), including misclassified zero-group days during the Maunder minimum; that is, days with no available observations were assigned a SNs group count of zero (Vaquero et al. 2016; Clette et al. 2023). For this reason, we used it only when no other, more reliable data exist (i.e. before 1739) and treated it as a lower limit for solar activity during the Maunder minimum. We emphasise, however, that the use of HoSc98 introduces increased uncertainty in our results for those early periods – a limitation shared by all studies extending over and prior to the Maunder minimum that rely on sunspot observations. The series by Carrasco et al. (2022, 2024, 2025) provides a more realistic estimate of solar activity during that time, but its annual resolution limits its use for our purposes. We therefore employed it to obtain a rough upper-limit estimate for our reconstruction during the Maunder minimum. ISNv2 values in the first decades of the 18th century may also be somewhat overestimated (see e.g. Usoskin et al. 2021; Carrasco et al. 2024).

To estimate the fractional area of the solar disc covered by sunspots (Sect. 3.2), we used the sunspot area record by Mandal et al. (2020), updated with recent observations2. We extended this record backwards to 1650 by fitting a second-degree polynomial to the corresponding sunspot number series and filled gaps through interpolation. After 2010, missing values were replaced with areas derived from SDO/HMI (Helioseismic and Magnetic Imager on board the Solar Dynamics Observatory) continuum observations (Chatzistergos et al. 2025a).

2.2. Data used for the optimisation process

As is discussed in Sects. 3 and 4, model parameters were constrained by comparing model outputs to the following independent datasets:

• The total magnetic flux from Carrington-rotation (CR) synoptic magnetograms obtained at the Wilcox Solar Observatory (WSO), Mount Wilson Observatory (MWO), and Kitt Peak National Solar Observatory (NSO/KP) (Arge et al. 2002; Wenzler et al. 2006). For the period 1976–2002, all three datasets are available and have been averaged; for 2002–2009, we averaged data from WSO and MWO.

• Open solar magnetic flux (OSF) from the empirical reconstruction based on the geomagnetic aa-index by Lockwood & Owens (2024), available for the period 1868–2022.8 as both CR and 13-CR averages. For optimisation, we employed the 1-CR dataset.

• The TSI composite of measurements by Montillet et al. (2022)³, covering the period 1980–2025 at daily cadence. For consistency, we scaled it to match the TSIS-1/TIM (Total and Spectral Solar Irradiance Sensor/Total Irradiance Monitor) absolute level at the minimum of cycle 24 (2019.96 ± 0.5 yr).

2.3. Data used for validation and comparisons

To assess the performance of our model, we compared our OSF, TSI, and SSI reconstructions with several datasets not used in the optimization.

2.3.1. Open magnetic flux

• The OSF reconstruction of Wu et al. (2018), derived from cosmogenic isotopes (thus independent of magnetograms and sunspot data).

• In situ OSF measurements of Owens et al. (2017), for 1998–2018 at an annual resolution, based on direct spacecraft measurements of the radial magnetic field.

• The OSF estimates, also from Owens et al. (2017), derived from the OMNI interplanetary magnetic-field data for 1964–1998 at an annual resolution, corrected for solar-wind kinematic effects following Lockwood et al. (2009). These provide an observationally based but proxy-corrected estimate of OSF. For quantitative comparison, we use only the in situ part of the Owens et al. (2017) series, but the full series is still shown as an additional reference.

2.3.2. Total solar irradiance

• The recent high-precision TSI measurements from TSIS-1/TIM⁴ (2018–2025.5; Pilewskie et al. 2018).

• SORCE/TIM (Total Irradiance Monitor on board the SOlar Radiation and Climate Experiment) data⁵ (2003.15–2020.15 Kopp et al. 2005).

• The TSI composite by Dudok de Wit et al. (2017, 1978.9 to 2025.1)⁶.

• The TSI measurements from Variability of Irradiance and Gravity Oscillations (VIRGO; 1996.14–2025.4; Fröhlich et al. 1995; Finsterle et al. 2021)⁷ instrument on board the Solar and Heliospheric Observatory (SoHO; Domingo et al. 1995).

• The TSI reconstruction of Wu et al. (2018), which is based on an earlier version of the SATIRE-T model.

• The TSI reconstruction of Chatzistergos et al. (2025a), which employs the SATIRE-S model.

2.3.3. Spectral solar irradiance

• The Lyman-α composite of Machol et al. (2019), which combines data from multiple instruments and fills data gaps using solar activity proxies: the 10.7 cm radio flux (1947–1955), the 30 cm radio flux (1955–1978), and the Mg II index for the remaining gaps.

• The SSI reconstruction from Chatzistergos et al. (2025a).

3. Model

We used the SATIRE-T model to reconstruct solar irradiance variability. The model requires the distribution and the temporal evolution of magnetic features on the solar surface as input. To obtain these, we first reconstructed the evolution of the solar surface magnetic flux using the model of the magnetic field evolution, originally proposed by Solanki et al. (2000), subsequently refined by Solanki et al. (2002a) and Vieira & Solanki (2010), and most recently revised by Krivova et al. (2021). This latest revision explicitly accounts for the contribution of small-scale magnetic features emerging outside active regions (ARs), as is described in Sect. 3.1. The reconstructed magnetic fluxes of the various components are then used as input to SATIRE-T to compute the TSI and SSI. Although SATIRE-T has been described in detail in earlier papers (e.g., Krivova et al. 2007, 2010; Wu et al. 2018), we summarise its concept and key components in Sect. 3.2 for completeness.

3.1. Magnetic field evolution

To model the evolution of the solar magnetic field over the telescopic era, we used the model by Krivova et al. (2021), which takes the sunspot number as input and describes the temporal evolution of the solar surface magnetic flux by solving a set of coupled ordinary differential equations. The model divides the entire range of the emerging magnetic regions into two classes: (1) ARs, containing sunspots, with magnetic fluxes above 4 × 10²⁰ Mx; and (2) small-scale emergences (SSEs), which do not contain sufficient flux to form sunspots. Active regions include both sunspots and faculae. The SSEs include ephemeral regions (ERs, bipolar magnetic regions with typical fluxes of about 10¹⁸ Mx) and the yet smaller internetwork fields, with fluxes down to 10¹⁶ Mx (see Krivova et al. 2021 for further details, as well as van Driel-Gesztelyi & Green 2015). The number of the smallest internetwork elements (10¹⁶ − 10¹⁷ Mx) remains nearly constant over the solar cycle (Buehler et al. 2013; Lites et al. 2014).

Frozen in the ionized solar wind plasma, part of the surface magnetic flux is carried outwards, forming the OSF, ϕ_open. Following Vieira & Solanki (2010) and Krivova et al. (2021), the model differentiates between the rapidly, ϕ_open^r, and slowly, ϕ_open^s, evolving components of the OSF. The rapid component arises mainly from small coronal holes associated with ARs, which form and disappear quickly on timescales of days to weeks. The slowly evolving component is associated with large-scale polar coronal holes sustained by both ARs and SSEs. Its characteristic timescale range from months to years.

The evolution of the various magnetic flux components in time, t, is described by the set of coupled differential equations:

$$\begin{array}{cc}& \frac{\mathrm{d}{\varphi }_{\mathrm{AR}}}{\mathrm{d}t}={\epsilon }_{\mathrm{AR}}-\frac{{\varphi }_{\mathrm{AR}}}{{\tau }_{\mathrm{AR}}^0}-\frac{{\varphi }_{\mathrm{AR}}}{{\tau }_{\mathrm{AR}}^{\mathrm{s}}}-\frac{{\varphi }_{\mathrm{AR}}}{{\tau }_{\mathrm{AR}}^{\mathrm{r}}},\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\mathrm{d}{\varphi }_{\mathrm{SSE}}}{\mathrm{d}t}={\epsilon }_{\mathrm{SSE}}-\frac{{\varphi }_{\mathrm{SSE}}}{{\tau }_{\mathrm{SSE}}^0}-\frac{{\varphi }_{\mathrm{SSE}}}{{\tau }_{\mathrm{SSE}}^{\mathrm{s}}},\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\mathrm{d}{\varphi }_{\mathrm{open}}^{\mathrm{r}}}{\mathrm{d}t}=\frac{{\varphi }_{\mathrm{AR}}}{{\tau }_{\mathrm{AR}}^{\mathrm{r}}}-\frac{{\varphi }_{\mathrm{open}}^{\mathrm{r}}}{{\tau }_{\mathrm{open}}^{\mathrm{r}}},\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\mathrm{d}{\varphi }_{\mathrm{open}}^{\mathrm{s}}}{\mathrm{d}t}=\frac{{\varphi }_{\mathrm{AR}}}{{\tau }_{\mathrm{AR}}^{\mathrm{s}}}+\frac{{\varphi }_{\mathrm{SSE}}}{{\tau }_{\mathrm{SSE}}^{\mathrm{s}}}-\frac{{\varphi }_{\mathrm{open}}^{\mathrm{s}}}{{\tau }_{\mathrm{open}}^{\mathrm{s}}},\hfill \end{array}$$(4)

$$\begin{array}{cc}& {\varphi }_{\mathrm{open}}={\varphi }_{\mathrm{open}}^{\mathrm{r}}+{\varphi }_{\mathrm{open}}^{\mathrm{s}},\hfill \end{array}$$(5)

$$\begin{array}{cc}& {\varphi }_{\mathrm{total}}={\varphi }_{\mathrm{AR}}+{\varphi }_{\mathrm{SSE}}+{\varphi }_{\mathrm{open}}.\hfill \end{array}$$(6)

Here, ϕ_AR and ϕ_SSE denote the magnetic flux in ARs and SSEs, respectively. The first two equations describe the temporal balance between the emergence (ε_AR and ε_SSE) and decay of magnetic flux in ARs and SSEs, respectively. The following two equations account for the transfer of magnetic flux from closed-field regions to the OSF, subdivided into the rapidly, ϕ_open^r, and slowly, ϕ_open^s, evolving components, whose sum gives the total OSF, ϕ_open. The sum of all components provides the total magnetic flux, ϕ_total, at the solar surface and open into the heliosphere. The parameters τ_AR⁰ and τ_SSE⁰ denote the decay timescales of flux in ARs and SSEs, respectively, while τ_AR^s, τ_AR^r, and τ_SSE^s specify transfer timescales of the AR and ER flux to the open-flux components. Finally, τ_open^r and τ_open^s describe the decay of the rapid and slow components of the OSF, respectively.

Using observational constraints from Parnell et al. (2009) and Thornton & Parnell (2011), Krivova et al. (2021) estimated the typical lifetimes of magnetic flux in ARs and SSEs to be approximately 10 days for τ_AR⁰ and 6 minutes for τ_SSE⁰, respectively. These values were derived by dividing the number of observed magnetic features (with fluxes above or below ϕ_AR for AR and SSE regions, respectively) by their corresponding emergence rates.

The emergence rates of magnetic flux in ARs and SSEs, ε_AR and ε_SSE, respectively, follow the observed flux-emergence distribution derived by Thornton & Parnell (2011):

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}\varphi }=\frac{n_0}{{\varphi }_0}{\left(\frac{\varphi }{{\varphi }_0}\right)}^m,\end{array}$$(7)

where ϕ₀ = 10¹⁶ Mx is the flux in the smallest considered features, n₀ = 3.14 × 10⁻¹⁴ cm⁻²day⁻¹ is the total number density of emergence events, and m = −2.69 is the slope of the distribution.

However, flux emergence is not constant over the solar cycle. Harvey & Zwaan (1993) showed that the emergence rate of large magnetic regions, with areas exceeding 3.5 square degrees (roughly corresponding to AR fluxes of ≥4 × 10²⁰ Mx; Hofer et al. 2024), varied by a factor of ∼8.3 between the solar maximum and minimum of solar cycle 21. In contrast, smaller magnetic regions with areas of 2.5–3.5 square degrees (that is, SSEs with fluxes of ∼3 × 10¹⁸ − 4 × 10²⁰ Mx) varied only by a factor of ∼2.1. This behaviour is consistent with the observed quadratic relationship between sunspot and facular areas (e.g., Shapiro et al. 2014; Yeo et al. 2020a; Chatzistergos et al. 2022a).

To account for this solar cycle–dependent variability, Krivova et al. (2021) allowed the slope, m, of the emergence-rate distribution to vary with activity (quantified by the SN):

$$\begin{array}{c}\hfill m(\mathrm{SN})=m_1-({\mathrm{SN}}_1^{\alpha }-{\mathrm{SN}}^{\alpha })\frac{m_1-m_2}{{\mathrm{SN}}_1^{\alpha }-{\mathrm{SN}}_2^{\alpha }},\end{array}$$(8)

where α is a free parameter constrained through comparison with observations.

Following Krivova et al. (2021), we adopted SN₁  =  217 and SN₂  =  17, corresponding to the maximum and minimum (average over the preceding and following minima) of cycle 21, consistent with the periods analysed by Harvey & Zwaan (1993). To reproduce the observed (Harvey & Zwaan 1993) ratios in emergence rates, the slope is allowed to vary symmetrically about the mean value, m, by an amount, Δm = 0.0946, yielding m₁ = −2.59 and m₂ = −2.78.

The emergence rates of ARs and SSEs were then computed by integrating Eq. (7) over the respective flux ranges:

$$\begin{array}{c}\hfill {\epsilon }_{\mathrm{AR}}={\displaystyle {\int }_{{\varphi }_{\mathrm{AR}}}^{{\varphi }_{\lim }}}\frac{n_0}{{\varphi }_0}{\left(\frac{\varphi }{{\varphi }_0}\right)}^m\varphi \mathrm{d}\varphi =\frac{n_0}{(m+2){\varphi }_0^{(m+1)}}({\varphi }_{\lim }^{(m+2)}-{\varphi }_{\mathrm{AR}}^{(m+2)})\end{array}$$(9)

and

$$\begin{array}{c}\hfill {\epsilon }_{\mathrm{SSE}}={\displaystyle {\int }_{{\varphi }_0}^{{\varphi }_{\mathrm{AR}}}}\frac{n_0}{{\varphi }_0}{\left(\frac{\varphi }{{\varphi }_0}\right)}^m\varphi \mathrm{d}\varphi =\frac{n_0}{(m+2){\varphi }_0^{(m+1)}}({\varphi }_{\mathrm{AR}}^{(m+2)}-{\varphi }_0^{(m+2)}).\end{array}$$(10)

Following Krivova et al. (2021) and Hofer et al. (2024), we used ϕ_AR = 4 × 10²⁰ Mx, corresponding to the smallest regions capable of producing sunspots, and ϕ_lim = 10²³ Mx as the upper limit (cf. Thornton & Parnell 2011). Regions above this limit are extremely rare, so the exact choice of ϕ_lim has a negligible influence (Hofer et al. 2024). We also recall that ϕ₀ = 10¹⁶ Mx corresponds to typical internetwork fluxes.

3.2. Modeling irradiance variability

Variations in solar irradiance on timescales of days to millennia are caused by the evolution of the surface magnetic field through the competing effects of sunspot darkening and facular brightening. The SATIRE model computes irradiance as the sum of contributions from five surface components: sunspot umbrae (u), penumbra (p), faculae (f), network (n), and the quiet Sun (q). With the quiet Sun, we refer to surface areas free of a (detectable) magnetic field.

Thus, the solar irradiance at a given time and wavelength was computed as

$$\begin{array}{cc}\hfill F(\lambda ,t)& ={\alpha }_{\mathrm{q}}(t)F_{\mathrm{q}}(\lambda )+{\alpha }_{\mathrm{u}}(t)F_{\mathrm{u}}(\lambda )\hfill \\ \hfill & +{\alpha }_{\mathrm{p}}(t)F_{\mathrm{p}}(\lambda )+[{\alpha }_{\mathrm{f}}(t)+{\alpha }_{\mathrm{n}}(t)]F_{\mathrm{f}}(\lambda ),\hfill \end{array}$$(11)

where α_i are the fractional disc coverages (filling factors) of component i and F_i the corresponding disc-integrated brightness spectra. The brightness spectra are time-independent but depend on the wavelength and limb angle.

The filling factors were derived from the reconstructed magnetic fluxes (Sect. 3.1). Sunspot filling factors (fractional surface areas) were taken directly from observations (see Sect. 2) and divided into umbra and penumbra with a fixed ratio of α_u/(α_u + α_p) = 0.2, consistent with the mean observed value (Brandt et al. 1990; Solanki 2003). As has been shown by Chatzistergos et al. (2021), the precise umbra–penumbra partition has only a minute effect on the reconstruction. Assuming a mean vertical magnetic field component of ⟨B_u⟩ = 1800 G in umbrae and ⟨B_p⟩ = 550 G in penumbrae (Keppens & Martinez Pillet 1996), the total magnetic flux in sunspots is

$$\begin{array}{c}\hfill {\varphi }_{\mathrm{s}}={\varphi }_{\mathrm{u}}+{\varphi }_{\mathrm{p}}={\alpha }_{\mathrm{u}}⟨B_{\mathrm{u}}⟩+{\alpha }_{\mathrm{p}}⟨B_{\mathrm{p}}⟩.\end{array}$$(12)

Since ARs include sunspots and faculae, the facular flux is the remaining part of the total AR flux: ϕ_f = ϕ_AR − ϕ_s. The magnetic flux contributing to the network includes both the SSE flux and the OSF: ϕ_n = ϕ_SSE + ϕ_open. The filling factors for faculae and network were computed assuming a linear increase with mean flux density until a saturation threshold, above which they remain constant at unity:

$$\begin{array}{c}\hfill {\alpha }_i=min(\frac{B_i}{B_{\mathrm{sat},i}},1),i\in \{\mathrm{f},\mathrm{n}\},\end{array}$$(13)

where B_i = ϕ_i/A_⊙ and A_⊙ denotes the area of the solar hemisphere. For the network we adopted B_sat, n = 800 G, as estimated by Krivova et al. (2007, 2010) from the analyses of magnetograms. The saturation threshold for faculae, B_sat, f, is a free parameter determined from a comparison to observations, in the optimisation described below. For the disc-averaged flux densities used here, the saturation limit is never reached. The quiet-Sun filling factor is then the remaining fraction of the solar surface not occupied by magnetic features, α_q = 1 − α_u − α_p − α_f − α_n.

The brightness spectra of the individual components are unchanged compared to the previous SATIRE implementations. They were computed by Unruh et al. (1999) using the ATLAS9 radiative transfer code (Kurucz 1993, 2005). For the quiet Sun, penumbra, and umbra, we used Kurucz’ model atmospheres with effective temperatures of 5777 K, 5400 K, and 4500 K, respectively. The facular spectrum is based on the FAL-P atmosphere (Fontenla et al. 1993), as modified by Unruh et al. (1999). These spectra represent the emergent intensities, I_i(λ, μ), where μ is the cosine of the heliocentric angle, which we integrated over the visible disc to obtain the disc-integrated component spectra, F_i(λ).

Faculae and sunspots are not evenly distributed over the disc but appear predominantly within the active-region belts: ∼5°–30° for spots and ∼5°–45° for faculae. Therefore, a uniform disc integration would misrepresent their centre-to-limb behaviour. Following Krivova et al. (2007), Wu et al. (2018), we integrated the intensity spectra over the corresponding latitude belts and rescaled the disc-integrated filling factors to the belt areas.

Spectra of the various components were computed under the assumption of local thermodynamic equilibrium (LTE), which becomes increasingly inaccurate below about 300 nm. To remedy this, we applied an empirical correction following the same procedure as in earlier SATIRE versions (Krivova et al. 2006; Yeo et al. 2014; Chatzistergos et al. 2025a). It has been validated against non-LTE calculations by Tagirov et al. (2019). Between 180 and 300 nm, the reconstructed SSI was adjusted to match the whole heliospheric interval (WHI) reference spectrum (Woods et al. 2009); between 115 and 180 nm, it was linearly scaled to agree with SORCE/SOLSTICE (The SOlar Stellar Irradiance Comparison Experiment on board of Solar Radiation and Climate Experiment) observations. These corrections are time-independent, and therefore do not affect the long-term variability.

3.3. Optimisation of model parameters

Several free parameters enter the magnetic-flux evolution (Sect. 3.1) and the irradiance computation (Sect. 3.2). Their values are constrained by comparison with three independent datasets (Sect. 2.2): the total photospheric magnetic flux (Arge et al. 2002; Wenzler et al. 2006), the OSF reconstruction of Lockwood & Owens (2024), and the TSI composite of Montillet et al. (2022). We employed the genetic-algorithm code PIKAIA (Charbonneau 1995) to minimise the combined reduced χ² between the model output and these observations. The resulting best-fit parameter set is listed in the bottom part of Table 1. The top part lists the values of the fixed parameters, which are the same for all runs.

4. Results

4.1. Magnetic flux

The fractional disc coverages used in the irradiance reconstruction (Sect. 3.2) were derived from the modelled evolution of the solar surface magnetic flux. Because several aspects of the modelling differ from the earlier study of Krivova et al. (2021), including updated input and reference datasets as well as the inclusion of TSI in the optimisation, we first compare our reconstructed total and open magnetic flux with the corresponding observations or independent reconstructions, as well as with the Krivova et al. (2021) version. The free parameters of the model, optimised as described in Sect. 3.3, are listed in the bottom part of Table 1.

To evaluate the reconstruction of the total photospheric magnetic flux, we compare the total unsigned flux with the combined measurements from WSO, MWO, and NSO/KP (Sect. 2.2), all at a CR cadence. Our two reconstructions – ISNv2-based (top) and CEA17-based (bottom) – are shown in Fig. 1 as 27-day (roughly one CR) means and are compared with the observational data, also at one-CR cadence. Krivova & Solanki (2004) showed that more than a half of the quiet-Sun flux remains undetected in NSO/KP synoptic charts due to the limited spatial resolution. Therefore, following Krivova et al. (2007, 2021), we scaled the SSE contribution by a factor of 0.4 when comparing with observations. We emphasise that this scaling was only applied for comparison purposes; when reconstructing the irradiance, all magnetic flux was retained. Our tests showed that the exact value of this scaling factor has only a weak effect on the comparison presented in Fig. 1. It mainly introduces a small overall shift in the absolute level but barely affects the variability, consistent with the fact that variations in the magnetic flux on timescales of a solar cycle and shorter on the Sun are dominated by ARs (Harvey & Zwaan 1993; Hofer et al. 2024, see also the discussion of Fig. 2 below).

Fig. 1. Evolution of the total unsigned photospheric magnetic flux over time. Red: Our reconstructions based on the ISNv2 (top) and CEA17 GSN (bottom). Blue: Reconstruction of Krivova et al. (2021) (based on ISNv2 in the top panel, HoSc98 GSN in the bottom panel). Symbols: Total magnetic flux measurements from WSO, MWO, and NSO/KP; black line: their average. All quantities are shown as one-CR averages. The contribution of the SSE is scaled by a factor of 0.4 to compensate for the limited spatial resolution of observations.

Fig. 2. Solar surface magnetic flux reconstruction based on ISNv2. Different components of surface magnetic flux are shown with different colours, both as 11-year moving averages and annual averages shown as thick and thin lines, respectively. Shown are the contributions from: active regions (ARs, blue), open solar flux (OSF, red), small-scale emergence regions (SSE, green), and the sum of OSF and SSE.

Both our reconstructions agree very well with the observations. Comparing the 27-day average of our reconstruction to the averaged observations, we obtain a Pearson correlation coefficient of 0.95 and 0.94, a root mean square error (RMSE) of 8.51 × 10²²Mx and 8.74 × 10²²Mx and a reduced χ² value of 0.07 and 0.065 for ISNv2- and CEA17 GSN-based models, respectively (see Table 2).

Around the maximum of solar cycle 21, the agreement between our reconstructed total magnetic flux and the observation-average becomes somewhat poorer. This period coincides with increased scatter among the magnetogram archives and with documented issues in the NSO/KP record before early 1992 (e.g. zero-point and calibration problems noted by Arge et al. 2002). To account for these issues, Wenzler et al. (2006) applied a linear scaling to NSO/KP data over 1974–1987 to match the level of NSO data corrected by Arge et al. (2002), by multiplying them by 1.242, which might not have been necessary. When we instead use the NSO/KP values prior to that rescaling, the inter-archive scatter decreases, and the agreement with our reconstruction improves slightly (Fig. A.1). Using this ‘restored’ dataset, we obtain Pearson correlation coefficient of 0.95, an RMSE of 7.28 × 10²² Mx, and a reduced χ² of 0.054 for the ISNv2-based reconstruction, while for the CEA-based reconstruction we get a Pearson correlation coefficient of 0.95, RMSE of 7.45 × 10²² Mx, and χ² of 0.06. This suggests that the residual (model–data) differences during cycle 21 are comparable to the archive-to-archive calibration uncertainties.

For comparison, Fig. 1 also shows the total magnetic flux reconstruction by Krivova et al. (2021). Our ISNv2-based model (top panel, red) is in excellent agreement with the reconstruction of Krivova et al. (2021 blue). The blue line is hardly visible, hidden behind the red curve. Our CEA17-based reconstruction has lower amplitudes during solar cycles 21 and 22 than the one by Krivova et al. (2021) using the HoSc98 GSN. This difference is mainly due to the different GSN series used as inputs. Krivova et al. (2021) reported χ² values of 0.058 for both of their reconstructions, ISNv2- and HoSc98-based. The slightly higher values obtained here can be attributed to differences in the input sunspot number records, the datasets used for model optimisation, and stricter constraints introduced by including TSI into optimisation.

Figure 2 shows the ISNv2-based reconstruction of the individual AR (blue) and small-scale (SSE+OSF; black) components of the magnetic flux (annual means and 11-year running means as thin and thick lines, respectively). During the Maunder minimum (∼1650–1700) and Dalton minimum (∼1790–1830), most of the flux originates from SSEs, whereas the AR flux is very low (Dalton minimum) or nearly absent (Maunder minimum). The secular increase in the total flux is dominated by the AR component: the mean flux over 1967–2017 exceeds that over 1650–1700 by ∼3.6 × 10²³ Mx, of which ∼68% originates from ARs and 32% from the combined SSE and OSF components. The CEA17-based reconstruction shows a very similar behaviour (Appendix B, Fig. B.1).

Figure 3 compares the evolution of the computed OSF to that reconstructed by Lockwood & Owens (2024). The top panel presents the ISNv2-based reconstruction, and the bottom panel the CEA17-based reconstruction. For the optimisation of the model parameters (Sect. 3.3), we used the OSF dataset at a cadence of one CR, available from 1868 to 2022.

Fig. 3. Evolution of the open magnetic flux over time. Our reconstruction is shown in red for the ISNv2 input (top panel) and the CEA17 GSN (bottom panel), both at one-CR cadence. In dark and light green, we show the OSF from Lockwood & Owens (2024), reconstructed from the geomagnetic aa-index for 13-CR and 1-CR cadence, respectively. The dashed blue line shows the in situ OSF measurements from Owens et al. (2017) and the dashed black line shows the OSF estimates from Lockwood et al. (2009).

To quantify the performance of our model, we calculated reduced χ² values relative to the one-CR OSF from Lockwood & Owens (2024). Both of our reconstructions (ISNv2- and CEA17-based) show a good agreement with these data. We obtain reduced χ² values of 0.321 and 0.322 for our ISNv2 and CEA17 GSN-based reconstructions, respectively, when compared to the one-CR dataset. Comparing both our reconstructions with the 13-CR series, we get χ² values of 0.228 and 0.225 for the ISNv2 and CEA17 GSN reconstructions, respectively (Table 2). For the ISNv2-based model, the Pearson correlation with the 1-CR OSF is 0.7 and the RMSE of 1.44 × 10²²Mx, while for the CEA17-based model, these values are 0.7 and 1.45 × 10²²Mx, respectively.

Our reconstructions also agree well with the directly measured heliospheric magnetic flux from Owens et al. (2017, dashed blue line in Fig. 3), available at an annual resolution from 1998 to 2018. Using annual means, we obtain for the ISNv2-based model a Pearson correlation coefficient of 0.87, a χ² value of 0.313, and a RMSE of 0.66 × 10²²Mx. For the CEA17 GSN-based model, the corresponding values are a Pearson correlation coefficient of 0.86, an χ-square value of 0.35, and a RMSE of 0.72 × 10²²Mx. For completeness, we also show the earlier part of the Owens et al. (2017) series, providing the OSF estimates derived from the OMNI interplanetary magnetic-field data for 1964–1998, corrected for solar-wind kinematic effects following Lockwood et al. (2009). These provide an observationally based but proxy-corrected estimate of OSF (dashed black line in Fig. 3).

The OSF reconstruction back to 1650 is shown in Fig. 4. For 1750–2020, where the sunspot input is identical, our ISNv2-based OSF agrees well with that of Krivova et al. (2021), with slightly lower cycle amplitudes due to the updated optimisation datasets. Before 1750, differences arise due to the different sunspot records used: Krivova et al. (2021) combined annual ISNv2 and Vaquero et al. (2015), whereas we employed daily HoSc98 GSN before 1749 and extended it further back with annual ISNv2 (to 1722) and the series of Carrasco et al. (2024). The CEA17-based reconstruction likewise diverges from the HoSc98-based one before ∼1870, reflecting differences in the underlying GSN series.

Fig. 4. OSF reconstruction from this study (red) and by Krivova et al. (2021, blue) using ISNv2 (top) and GSN (bottom; CEA17 here and HoSc98 in Krivova et al. 2021) as input. The shaded pink area shows our upper limit estimate for the Maunder minimum, where the SNs from Carrasco et al. (2025) was used as input. Black: OSF reconstructed from the cosmogenic isotope data by Wu et al. (2018).

We also compare our results with the OSF reconstructed from cosmogenic isotope data by Wu et al. (2018, black curve). The OSF from Wu et al. (2018) during the Maunder minimum lies well within our upper and lower estimates. Since the revised model includes SSE emergence, some OSF is maintained even when sunspots are absent. For 1650–1670 we obtain mean OSF values of 0.182 × 10²³ Mx (ISNv2) and 0.194 × 10²³ Mx (CEA17). Given the scarcity and uncertainty of observations in this period, these should be interpreted as lower limits.

4.2. Solar irradiance reconstruction

We now present our irradiance reconstructions. Model optimisation was performed using the TSI composite of Montillet et al. (2022), applying an 81-day moving average. The 81-day smoothing applied to TSI reduces the influence of short-lived irradiance dips caused by individual large sunspot groups. This ensures that the optimisation is driven primarily by variability on timescales of years and longer that the model is designed to reproduce.

Since the TSI composite of Montillet et al. (2022) differs in its absolute level from the TSIS-1/TIM record, we rescaled both the Montillet et al. (2022) composite and our reconstructed TSI to the TSIS-1/TIM level at the minimum of solar cycle 24. The timing of this minimum (2019.96) was taken from the SILSO ISNv2 database⁸. For each series, we computed the mean of the daily TSI values within a 12-month window centred on 2019.96 (±6 months). The ratio of the TSIS-1/TIM mean to the Montillet et al. (2022) mean defines the scaling factor applied to the composite. A similar scaling was applied to the modelled daily TSI.

Figure 5a shows our reconstructed TSI based on the ISNv2 in comparison with the TSI composite of Montillet et al. (2022), both as daily values and their 81-day running means. We find a good agreement between them, with Pearson correlation coefficients of 0.92 for the 81-day smoothed series and 0.80 for the daily values. The corresponding RMSE are 0.17 W m⁻² for the 81-day series and 0.32 W m⁻² for the daily one. The reduced χ² values for the 81-day smoothed comparison is 0.074. Panel (b) presents the residuals between our reconstruction and the Montillet et al. (2022) composite, daily (dots) and smoothed (line). A linear fit to the daily residuals yields a trend of −3.57 ±  0.2 mW m⁻² yr⁻¹. This indicates that our reconstruction reproduces the long-term change of the satellite record without any significant secular bias.

Fig. 5. a) Reconstructed TSI using ISNv2 as input (pink: daily values; red: 81-day smoothing). Blue dots and blue line: TSI composite of Montillet et al. (2022, daily and 81-day smoothing, respectively). b) Difference between our reconstruction and the Montillet et al. (2022) composite (blue dots: daily; blue line: 81-day smoothing). The horizontal dashed black line marks zero, and the orange line is the linear fit to the daily residuals. c) and d) Same as a) and b), but using the CEA17 GSN as an input to the model. e) Our TSI reconstructions based on ISNv2 and CEA17 (red and black). For comparison, the Montillet et al. (blue 2022) composite, the TSI reconstruction of Wu et al. (green 2018), and that of Chatzistergos et al. (orange 2025a) are shown. All datasets in panel e) are shown as 81-day (faint lines) and annual moving means (solid lines).

Figure 5 panel (c) shows the same comparison as panel (a) for our CEA17 GSN-based reconstruction. Similarly to the ISNv2-based case, we find a good agreement with observations, with a Pearson’s correlation coefficient of 0.92 for the 81-day smoothing and 0.8 for the daily values. The RMSE is 0.32 W m⁻² for the daily values and 0.18 W m⁻² for the 81-day smoothed ones. The reduced χ² for the 81-day smoothing is 0.079. A linear regression fit to the daily residuals gives a trend of −2.19 mW m⁻² yr⁻¹, with an uncertainty of ± 0.2 mW m⁻² yr⁻¹.

Further, we have compared the output of our model with the previous SATIRE-T version (Wu et al. 2018) and with SATIRE-S (Chatzistergos et al. 2025a), all shown in panel (e) of Fig. 5. The Wu et al. (2018) model is based on the previous SATIRE-T formulation and uses an earlier magnetic-flux evolution scheme by Vieira & Solanki (2010). The SATIRE-S model of Chatzistergos et al. (2025a), in contrast, is the most accurate version of the SATIRE family, as it employs state-of-the-art full-disc magnetograms and continuum images. It is fully independent of sunspot numbers.

We also evaluated all these versions against the Montillet et al. (2022) composite over their overlap period (1980–2017); see Table 3. Compared to the Wu et al. (2018) reconstruction, our model gives comparable correlation coefficients (marginally lower or higher for a daily or 81-day smoothed cadence), but a considerably lower χ² and a significantly weaker trend in the daily residuals (see Table 3), indicating an improved representation of long-term changes. As was expected, the SATIRE-S model Chatzistergos et al. (2025a) performs best overall. Nevertheless, our model matches SATIRE-S remarkably closely, despite relying solely on the sunspot-number input. This close agreement supports both the present model and the updated magnetic-flux evolution description of Krivova et al. (2021).

To evaluate the performance of our model more broadly, we compared our TSI reconstructions with several independent TSI datasets: the VIRGO TSI measurements (Fröhlich et al. 1995; Finsterle et al. 2021), the composite of Dudok de Wit et al. (2017), TSIS-1/TIM (Pilewskie et al. 2018), and SORCE/TIM measurements (Kopp et al. 2005). Pearson correlation coefficients, computed for both daily and 81-day–smoothed values, are summarised in Table 4. For the 81-day-smoothed series, the correlation with the various TSI datasets ranges between 0.81 and 0.98, while daily correlations range between 0.69 and 0.85. The best agreement is obtained with the TSIS-1/TIM and VIRGO measurements.

Figure 6 presents the reconstructed TSI back to 1650 for both input series. The ISNv2- and CEA17-based reconstructions are nearly identical before 1720 and after 1879, but differ between these periods, with the CEA17-based reconstruction showing weaker cycle amplitudes. Despite these differences, all three reconstructions (ISNv2-based, CEA17-based, and the upper-limit estimate) give a consistent TSI increase between the Maunder minimum and recent activity minima.

Fig. 6. Reconstructed TSI (annual means) since 1700 using ISNv2 (red) and CEA17 GSN (black) as inputs. The shaded pink area shows our upper limit estimate in TSI, where the Carrasco et al. (2025) SNs were used. The TSI composite of Montillet et al. (2022) is shown in blue; the TSI reconstruction of Wu et al. (2018) in green.

Comparing our reconstructed TSI during the Maunder minimum (here defined as the 50-year mean over 1650–1700) with the recent period (1967–2017; 50-year mean), we find a TSI increase of 0.75 W/ m² for the lower-level estimate and 0.67 W/ m² for the upper-level estimate in the ISNv2-based reconstruction. For the reconstruction using the CEA17 GSN series, we obtain an increase of 0.73 W/ m². For comparison, Wu et al. (2018) found an increase of 0.74 W/ m². Overall, our reconstructions lie towards the lower end of published estimates of the TSI change since the Maunder minimum, which reach up to 2.3 W/ m² (Penza et al. 2024) and even 5.3 W/ m² (Egorova et al. 2018). For consistency, these values were computed from the respective published TSI series using the same definition adopted here; namely, the difference between the 50-year mean over 1650–1700 and that over 1967–2017. The inferred secular change remains below the upper limit of 2.0 ± 0.7 W/ m² relative to the 2019 minimum derived by Yeo et al. (2020b), who constrained the minimum irradiance state of the Sun during a grand minimum using state-of-the-art magnetohydrodynamic simulations.

Figure 7 illustrates how the individual magnetic components contribute to the reconstructed TSI in the ISNv2-based model. Both ARs and the small-scale regions outside ARs (here, the sum of SSEs and the OSF) contribute to the secular change in TSI. During the Maunder minimum, when sunspots were essentially absent, the TSI variability originates mostly from the small-scale magnetic flux, while the AR contribution is negligible. During the less deep and shorter Dalton minimum at the beginning of the 19th century, the two components are comparable. Outside grand minimum conditions, ARs dominate nearly all solar cycle and multidecadal variability and provide the largest part of the long-term increase in TSI since the Maunder minimum. Quantitatively, from the Maunder minimum (average over 1650–1700) and the modern period (1967–2017), we obtain a TSI increase of 0.54 W/m² from ARs flux and 0.21 W/m² difference from the small-scale flux. This means that ∼72% of the secular rise originated from ARs and ∼28% from the small-scale flux.

Fig. 7. Total solar irradiance reconstruction based on ISNv2 as an input (red). The green and black curves show the contribution of different components of magnetic flux to TSI. Green shows the contribution of AR magnetic flux to TSI, and black shows the contribution of SSEs. Thin and thick lines show the annual means or 11-year running means, respectively.

The same behaviour was found when using the CEA17 GSN record as an input (see Appendix A, Fig. B.2). In that case, the TSI increase since the Maunder minimum is 0.53 W/m² (∼74%) for the AR component and 0.19 W/m² ∼26% from regions outside ARs. Thus, while the small-scale component determines the baseline irradiance during extended periods of very low activity, the long-term rise in TSI since the Maunder minimum is governed mainly by the evolution of ARs.

This contrasts with the irradiance modelling approach of Penza et al. (2024), in which a substantial fraction of the secular variability is attributed to the network component through an empirically prescribed long-term modulation. In our model, the secular variability arises self-consistently from the emergence, evolution, and decay of magnetic flux, including weak network and open flux contributions, without introducing an additional externally imposed long-term component.

The reconstruction of SSI is crucial for climate studies, particularly in the UV range, where irradiance variations strongly affect the chemistry of the Earth’s atmosphere (Haigh 2007; Solanki et al. 2013). Of special importance is the Lyman-α line at 121.5 nm: it is the brightest emission line in the solar spectrum and varies by about a factor of two over the solar cycle (Woods et al. 2000; Machol et al. 2019).

Figure 8 shows our Lyman-α reconstructions based on ISNv2 and CEA17 GSN, together with the observational composite of Machol et al. (2019). All series are displayed as annual running means. Both reconstructions agree very well with the composite, with correlation coefficients of 0.924 (ISNv2) and 0.921 (CEA17 GSN) at daily cadence. Importantly, no model parameters were explicitly adjusted to match the Lyman-α data.

Fig. 8. Reconstructed Lyman-α irradiance based on ISNv2 (red) and CEA17 (black) inputs. Also shown are the SATIRE-S Lyman-α reconstruction (Chatzistergos et al. 2025a, orange), and the observational composite by Machol et al. (2019, blue curve) with its uncertainties shown in light blue. The shaded grey area marks the period when the F30 data were used by Machol et al. (2019) instead of direct measurements. All series are shown as annual running means.

We also show in Fig. 8 the SATIRE-S Lyman-α reconstruction of Chatzistergos et al. (2025a). Despite being based on different model versions and on entirely independent input data, the two reconstructions agree closely. The Pearson correlation coefficients reach 0.92 for the daily ISNv2- and CEA17GSN-based models, respectively, and 0.99 for the annually smoothed series.

Both our (ISNv2- and CEA17 GSN-based) as well as the SATIRE-S (Chatzistergos et al. 2025a) reconstructions, reproduce the Lyman-α composite well, except over 1960–1975. During this interval (indicated as the shaded interval in the figure), no direct Lyman-α measurements were available, and the composite instead relies on values inferred from the 30 cm radio flux. This suggests that the discrepancy possibly originates from uncertainties in the Machol et al. (2019) composite during this period rather than from our model.

5. Conclusion

We have presented new reconstructions of TSI and SSI over the past four centuries using the updated physics-based SATIRE-T model. In contrast to earlier SATIRE-T implementations, which relied on a more empirical parameterisation of small-scale magnetic flux emergence, the present reconstruction incorporates a revised description of the solar surface magnetic-field evolution, including a more realistic prescription for the emergence of small-scale magnetic flux linked to sunspot activity. This provides a more self-consistent representation of irradiance variability, including its secular component.

Using two independent sunspot-number records (ISNv2 and CEA17 GSN) as inputs, we find that the resulting magnetic-flux and irradiance reconstructions are consistent with each other. Both versions reproduce the measured total photospheric magnetic flux from WSO, MWO, and NSO/KP, the open magnetic flux reconstructed from the geomagnetic aa-index (Lockwood & Owens 2024), and the directly measured heliospheric magnetic flux from Owens et al. (2017). The model also reproduces satellite TSI measurements (VIRGO, SORCE/TIM, TSIS-1/TIM, and the composite of Montillet et al. 2022) and the Lyman-α composite of Machol et al. (2019) without any additional tuning. The model achieves correlation coefficients of 0.81–0.98 with different space-based TSI records at 81-day resolution and 0.69–0.85 at daily cadence. The correlation coefficient with the daily Ly-α composite is 0.92.

On centennial timescales, both sunspot inputs give a relatively small increase in TSI since the Maunder minimum, between 0.67 and 0.75 W m⁻², depending on the sunspot record and in particular on the treatment of the data during the Maunder minimum. In our model, the small-scale magnetic flux dominates the surface magnetic flux and the irradiance level during extended low-activity periods, while the main contribution (about 60–70%) to the long-term rise from the Maunder minimum to the present comes from ARs.

The close agreement between the two versions of the reconstruction, based on independent sunspot-number series, and with independent datasets, together with the more realistic treatment of small-scale magnetic emergence, strengthens confidence in the updated SATIRE-T irradiance series. These new reconstructions provide a physically consistent record of solar irradiance variability over the last 400 years and can be used in studies of the solar influence on the climate.

The SATIRE-T irradiance reconstructions produced in this work supersede earlier SATIRE-T versions and are available to the research community (see data availability). For the period after 1976, we recommend the SATIRE-S version by Chatzistergos et al. (2025a), which is based on spatially resolved solar observations, and therefore provides a more direct and accurate representation of the surface magnetic-field distribution.

Data availability

All datasets used in this study are publicly available. The International Sunspot Number v2.0 can be obtained from the SILSO database (https://www.sidc.be/SILSO/datafiles). The CEA17 group sunspot number series and the updated sunspot-area record of Mandal et al. (2020) are available at https://www2.mps.mpg.de/projects/sun-climate/data.html. The HoSc98 group sunspot number series is accessible at the National Geophysical Data Center. The annual sunspot numbers derived with the active-day fraction method are available from the authors of Carrasco et al. (2022, 2024; 2025).

Magnetogram-based total magnetic flux data from WSO, MWO, and NSO/KP are publicly accessible through the respective observatories. The open-flux reconstruction from the geomagnetic aa-index (Lockwood & Owens 2024) is available from https://www.southampton.ac.uk/~gl/. The in situ heliospheric magnetic-flux data of Owens et al. (2017) and the OMNI-based open-flux estimates corrected following Lockwood et al. (2009) can be downloaded through NASA’s OMNI database (https://omniweb.gsfc.nasa.gov/). The TSI composite of Montillet et al. (2022), VIRGO TSI data, and the TSI composite of Dudok de Wit et al. (2017) are available from PMOD/WRC. SORCE/TIM and TSIS-1/TIM TSI data are available from the LASP data repository. The Lyman-α composite of Machol et al. (2019) is available from LASP. The SATIRE-S irradiance reconstruction of Chatzistergos et al. (2025a) is provided at https://www2.mps.mpg.de/projects/sun-climate/data.html.

The SATIRE-T irradiance reconstructions produced in this work are available under https://www2.mps.mpg.de/projects/sun-climate/data.html.

Acknowledgments

D.T. was supported by the International Max-Planck Research School (IMPRS) for Solar System Science at the Technical University of Braunschweig. S. K. S and T.C. have received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 101097844 – project WINSUN). This research has made use of the Astrophysics Data System (ADS; https://ui.adsabs.harvard.edu/) Bibliographic Services, funded by NASA under Cooperative Agreement 80NSSC21M00561.

References

Appendix A: NSO/KP calibration

Figure A.1 illustrates the comparison between the total magnetic flux reconstructed by the model and measured in magnetograms, discussed in Sect. 4.1. This figure is similar to Fig. 1, except that here NSO/KP data over 1974–1987 were re-scaled back to their original values from Arge et al. (2002) by dividing them by a factor 1.242 introduced by Wenzler et al. (2006), see Sect. 4.1. Restoring the NSO/KP data to the calibration of Arge et al. (2002) slightly improves the agreement with our model for both SNs inputs. Using this ‘restored’ dataset, we obtain Pearson correlation coefficients of 0.95 (ISNv2) and 0.95 (CEA17), with RMSE of 7.28 × 10²² Mx and 7.45 × 10²² Mx and χ² values of 0.054 and 0.06, respectively. This suggests that the additional linear scaling applied by Wenzler et al. (2006) may not be appropriate for all magnetograms.

Fig. A.1. Same as Fig. 1, but here we reduced the NSO KP level in the period 1974–1987 to the level before the application of the Wenzler et al. (2006) calibration by dividing the 1974–1987 data by the factor 1.242. The black line shows the average over the magnetogram archives whenever at least two datasets are available. Using this ‘restored’ NSO/KP dataset, we obtain Pearson correlation coefficients of 0.95 (ISNv2) and 0.95 (CEA17), with RMD values of -7.69% and -7.75%, respectively. This suggests that the additional linear scaling applied by Wenzler et al. (2006) may not be appropriate for all magnetograms.

Appendix B: Contribution of different magnetic flux components to total solar magnetic flux and total solar irradiance

Figure B.1 shows the CEA17-based reconstruction of the individual magnetic-flux components. As for the ISNv2-based case (Fig. 2), during the Maunder and Dalton minima, the solar surface magnetic flux is dominated by SSEs, with only partial (Dalton minimum) or negligible (Maunder minimum) AR contribution. The mean flux over 1967–2017 exceeds that over 1650–1700 by 3.6 × 10²³ Mx, with ∼69% coming from ARs and 31% from the sum of the SSE and OSF components.

Fig. B.1. Same as Fig. 2, but using CEA17 GSN as an input.

Figure B.2 shows the contribution of different magnetic flux components to the TSI reconstruction based on CEA17 GSN. We obtain similar results to those for the ISNv2-based reconstruction, where during periods of low activity (e.g., the Maunder minimum), the main contributor is the SSE flux, while for the rest of the time, the AR flux is the primary contributor to the TSI.

Fig. B.2. Same as Fig. 7, but here using the CEA17 GSN dataset as an input.

All Tables

All Figures

Fig. 1. Evolution of the total unsigned photospheric magnetic flux over time. Red: Our reconstructions based on the ISNv2 (top) and CEA17 GSN (bottom). Blue: Reconstruction of Krivova et al. (2021) (based on ISNv2 in the top panel, HoSc98 GSN in the bottom panel). Symbols: Total magnetic flux measurements from WSO, MWO, and NSO/KP; black line: their average. All quantities are shown as one-CR averages. The contribution of the SSE is scaled by a factor of 0.4 to compensate for the limited spatial resolution of observations.

In the text

Fig. 2. Solar surface magnetic flux reconstruction based on ISNv2. Different components of surface magnetic flux are shown with different colours, both as 11-year moving averages and annual averages shown as thick and thin lines, respectively. Shown are the contributions from: active regions (ARs, blue), open solar flux (OSF, red), small-scale emergence regions (SSE, green), and the sum of OSF and SSE.

In the text

Fig. 3. Evolution of the open magnetic flux over time. Our reconstruction is shown in red for the ISNv2 input (top panel) and the CEA17 GSN (bottom panel), both at one-CR cadence. In dark and light green, we show the OSF from Lockwood & Owens (2024), reconstructed from the geomagnetic aa-index for 13-CR and 1-CR cadence, respectively. The dashed blue line shows the in situ OSF measurements from Owens et al. (2017) and the dashed black line shows the OSF estimates from Lockwood et al. (2009).

In the text

	Fig. 4. OSF reconstruction from this study (red) and by Krivova et al. (2021, blue) using ISNv2 (top) and GSN (bottom; CEA17 here and HoSc98 in Krivova et al. 2021) as input. The shaded pink area shows our upper limit estimate for the Maunder minimum, where the SNs from Carrasco et al. (2025) was used as input. Black: OSF reconstructed from the cosmogenic isotope data by Wu et al. (2018).
In the text

Fig. 5. a) Reconstructed TSI using ISNv2 as input (pink: daily values; red: 81-day smoothing). Blue dots and blue line: TSI composite of Montillet et al. (2022, daily and 81-day smoothing, respectively). b) Difference between our reconstruction and the Montillet et al. (2022) composite (blue dots: daily; blue line: 81-day smoothing). The horizontal dashed black line marks zero, and the orange line is the linear fit to the daily residuals. c) and d) Same as a) and b), but using the CEA17 GSN as an input to the model. e) Our TSI reconstructions based on ISNv2 and CEA17 (red and black). For comparison, the Montillet et al. (blue 2022) composite, the TSI reconstruction of Wu et al. (green 2018), and that of Chatzistergos et al. (orange 2025a) are shown. All datasets in panel e) are shown as 81-day (faint lines) and annual moving means (solid lines).

In the text

	Fig. 6. Reconstructed TSI (annual means) since 1700 using ISNv2 (red) and CEA17 GSN (black) as inputs. The shaded pink area shows our upper limit estimate in TSI, where the Carrasco et al. (2025) SNs were used. The TSI composite of Montillet et al. (2022) is shown in blue; the TSI reconstruction of Wu et al. (2018) in green.
In the text

	Fig. 7. Total solar irradiance reconstruction based on ISNv2 as an input (red). The green and black curves show the contribution of different components of magnetic flux to TSI. Green shows the contribution of AR magnetic flux to TSI, and black shows the contribution of SSEs. Thin and thick lines show the annual means or 11-year running means, respectively.
In the text

Fig. 8. Reconstructed Lyman-α irradiance based on ISNv2 (red) and CEA17 (black) inputs. Also shown are the SATIRE-S Lyman-α reconstruction (Chatzistergos et al. 2025a, orange), and the observational composite by Machol et al. (2019, blue curve) with its uncertainties shown in light blue. The shaded grey area marks the period when the F30 data were used by Machol et al. (2019) instead of direct measurements. All series are shown as annual running means.

In the text

Fig. A.1. Same as Fig. 1, but here we reduced the NSO KP level in the period 1974–1987 to the level before the application of the Wenzler et al. (2006) calibration by dividing the 1974–1987 data by the factor 1.242. The black line shows the average over the magnetogram archives whenever at least two datasets are available. Using this ‘restored’ NSO/KP dataset, we obtain Pearson correlation coefficients of 0.95 (ISNv2) and 0.95 (CEA17), with RMD values of -7.69% and -7.75%, respectively. This suggests that the additional linear scaling applied by Wenzler et al. (2006) may not be appropriate for all magnetograms.

In the text

	Fig. B.1. Same as Fig. 2, but using CEA17 GSN as an input.
In the text

	Fig. B.2. Same as Fig. 7, but here using the CEA17 GSN dataset as an input.
In the text