© I. O. I. Virtanen et al. 2022

1. Introduction

Observations of the photospheric magnetic field were initiated in the early 20th century, when it became possible to measure the polarity as well as (albeit with limited accuracy) the strength of the magnetic field in sunspots (Hale et al. 1919). The first full-disk maps of the magnetic field (magnetograms) were made in the 1950s (Babcock 1953) and continuous daily full-disk observations have been conducted since the 1970s (Howard 1974; Livingston et al. 1976). Stanford University began providing daily full-disk maps in 1976 and observations are still underway at the Wilcox Solar Observatory (WSO, Scherrer et al. 1977). The National Solar Observatory (NSO) at Kitt Peak (KP) and the Mount Wilson observatory (MWO) have also created long continuous series of maps of the photospheric magnetic field at a higher spatial resolution than WSO (for a review of early full-disk magnetographic observations, see Pevtsov et al. 2021). Space-based instruments, such as SOHO/MDI and SDO/HMI, later emerged to provide observations without interference from the Earth’s atmosphere.

Continuous full-disk observations cover only four full sunspot cycles or two full magnetic cycles. Since only one side of the Sun is visible to instruments on the Earth or on a near-Earth orbit at any given time and due to the limited heliolatitude range of the Earth’s vantage point, measuring polar fields accurately is very difficult. For these reasons, our current ability to study the long-term evolution of the photospheric magnetic field is rather limited.

Surface flux transport (SFT) models (for review, see Sheeley 2005; Jiang et al. 2014) can be used to supplement observations. An SFT model takes flux emergence as input and simulates the evolution of the whole photospheric magnetic field based on the transport and diffusion of the emerged flux. Because flux emergence happens mainly in active regions that are relatively easy to observe accurately (since they contain very strong magnetic fields and emerge at relatively low latitudes), SFT simulations are a very helpful tool in studying the evolution of the field, especially the polar fields. Active regions required by the SFT model can be derived from sunspot numbers or areas (Jiang et al. 2010, 2011; Baumann et al. 2004) or assimilated directly from magnetographic observations (Virtanen et al. 2017, 2018; Whitbread et al. 2017, 2018; Yeates et al. 2015). Because sunspot numbers or areas alone contain no information about the structure of the active regions, methods based on them suffer from high uncertainties, while methods based on magnetographic observations are more accurate, but limited to the time period covered by full-disk magnetograms.

Pevtsov et al. (2016) presented a new method to reconstruct maps of the magnetic field of active regions from chromospheric plages and sunspot polarity observations. Daily full-disk observations of the chromospheric Ca II K 393.37 nm spectral line have been conducted since the early 20th century (starting in 1907 at the Kodaikanal Observatory in India and since 1915 at MWO, USA, see Bertello et al. 2016), and plages can be identified from them. Daily sunspots drawings containing (along with other information) the magnetic polarity and location of sunspots have been composed at MWO since 1917. These drawings were later digitized and the parameters tabulated (Pevtsov et al. 2019b; Tlatov et al. 2015). Together with the Ca II K line observations, the two datasets have allowed for the reconstruction of the magnetic fields of active regions since 1917 – and even since 1915, if we assume Hale polarity for active regions between 1915 and 1917. With an SFT model, the reconstruction can be extended to weaker fields outside active regions and near solar poles to cover the entire solar surface.

The reconstruction method has been tested by Virtanen et al. (2019). A reconstruction for the years 1975–1985 was completed and compared with magnetographic observations from the same time period. The reconstruction was found to be relatively accurate at long time scales, although at short time scales, it may deviate from observations considerably due to errors in individual reconstructed active regions. In this article, we extend the reconstruction to cover the time period from 1915 to 1985. We also use the potential field source surface (PFSS) model to reconstruct the coronal magnetic field. We study and discuss the active region reconstruction, the SFT simulation, and the modeled coronal field, then we compare them with observations when available.

The paper is structured as follows. In Sect. 2, we present the data used in this study. Section 3 provides a brief review of the SFT model. In Sect. 4, the active region reconstruction method and the resulting reconstructed active regions are presented. In Sect. 5, the SFT model is used to simulate the entire photospheric magnetic field. In Sect. 6, the coronal magnetic field is modeled and the amount of open field is studied. In Sect. 7, the reconstructed coronal field is compared with eclipse drawings. Section 8 presents our conclusions.

2. Data

Full-disk spectroheliograms of the Ca II K spectral line were taken daily at MWO from August 1915 to July 1985 and based on this data, the synoptic maps of Ca II K line intensity were compiled (Hale 1908; Ellerman 1919). Over the course of this MWO program, about 35 000 spectroheliograms were recorded. The typical number of observations per year is 180 days or better with smaller number of observing days in the first (1915 – about 75 days) and last year (1985 – about 90 days) of the program (see Fig. 1 in Bertello et al. 2020). Except for a small gap in 1939, there are no significant data gaps in the daily observations (see Fig. 2 in Bertello et al. 2020). Total number of observing days in 1939 is about 120. The spectroheliograms were recorded on plates until 1980 and on film henceforth. This change in the recording medium, along with many other poorly documented changes in the observational setup, caused the original dataset to be non-uniform. The data was renormalized by Pevtsov et al. (2016) and Bertello et al. (2020). These updated versions contain less errors and are more homogeneous. We use the renormalized dataset by Bertello & Pevtsov (2020) in our reconstruction process. The synoptic maps cover the time period from August 1915 (Carrington Rotation (CR) 827) to July 1985 (CR 1764). The resolution of the maps is 360 × 180 (longitude – sine of latitude).

Daily recordings of the properties of sunspots started at the MWO in 1917. The maximum strength of the magnetic field inside each sunspot was determined by visually finding the displacement between two Zeeman components of a spectral line. The Fe I 617.3 nm line was used until October 1961 and the Fe I 525.0 nm line since then (Pevtsov et al. 2016, 2019b). A drawing of the visible sunspots on the solar disk was made, and the strength of the magnetic field, its polarity, and other information about each sunspot was written by hand on the drawing. The observational setup has changed many times over the years, which has affected the field strength measurements and caused large uncertainties in the results. The polarity information, however, is very reliable (Pevtsov et al. 2019b). In this study we use the polarity, location, and size of sunspots, along with observation times. This information has been tabulated from the original drawings, as described by Tlatov et al. (2015) and Pevtsov et al. (2019b). MySQL database is accessible via Pevtsov et al. (2019a). Similar to Ca K line spectroheliogram observations, the dataset of sunspot daily recordings is quite uniform in respect to number of observing days per each year. Figure A.5 in Pevtsov et al. (2019b) shows two gaps in 1996–1997 and 2004–2007, with the former gap due to missing time information on the scanned images (daily observations were taken as normal) and the latter gap due to temporary program shut down. Both of these gaps are outside the time period included to our study. Some minor interruptions in daily observations took place as the spectrograph gratings were replaced in course of the project (for reference, see Table 1 in Pevtsov et al. 2019b). With the longest of these was caused by flooding of the spectrograph pit in the fall of 1978. None of these gaps had any significant impact on total number of observing days, as made evident in Fig. A.5 in Pevtsov et al. (2019b).

The photospheric line-of-sight magnetic field has been measured at National Solar Observatory at Kitt Peak since the 1970s. From February 1975 (CR 1625) to March 1992 (CR 1853) observations were made with a 512-channel diode array magnetograph, which measures the Fe I 868.8 nm spectral line (Livingston et al. 1976). The full-disk observations were processed into pseudoradial synoptic maps, which have a resolution of one degree in longitude and 1/90 in sine of latitude and with the field assumed to be purely radial in the photosphere. We used the NSO/KP synoptic maps from February 1975 (CR 1625) until July 1985 (CR 1763) to compare them with the reconstruction. It should be noted that the polar fields of these maps suffer from problems that are mostly related to poor pole filling (Harvey & Munoz-Jaramillo 2015; Virtanen & Mursula 2016). During this period, there is a gap in observations for CR1640–1644 (April 1976 – August 1976). Small gaps in longitudinal coverage (less than 1/3 of the map) can be found in CRs 1625, 1631, 1637, 1650, 1658, 1663, 1665, 1713, 1837, 1860, 1973, and 1981), while larger gaps (≥1/3 of map) are present in CRs 1632, 1635, 1639, 1647, 1648, and 1661).

3. SFT model

The SFT model in this work has been used and validated in our previous studies (Yeates et al. 2015; Virtanen et al. 2017, 2018). Here, we include only a short description of the model. For a more detailed description, testing, and validation, we refer to Yeates et al. (2015), and Virtanen et al. (2017, 2019).

The radial magnetic field, B_r(θ, ϕ, t), is written in the spherical coordinate system in terms of the two-component vector potential [A_θ, A_ϕ]:

$$\begin{array}{c}\hfill B_{\mathrm{r}}(\theta ,\varphi ,t)=\frac{1}{Rsin(\theta )}(\frac{\partial }{\partial \theta }(sin(\theta )A_{\varphi })-\frac{\partial A_{\theta }}{\partial \varphi }),\end{array}$$(1)

where R is the radius of the Sun, θ is colatitude, ϕ is the azimuth angle, and t is time. The vector potential evolves in time as follows:

$$\begin{array}{cc}& \frac{\partial A_{\theta }}{\partial t}=\omega (\theta )Rsin(\theta )B_{\mathrm{r}}(\theta ,\varphi ,t)-\frac{D}{Rsin(\theta )}\frac{\partial B_{\mathrm{r}}(\theta ,\varphi ,t)}{\partial \varphi }+S_{\theta }(\theta ,\varphi ,t),\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial A_{\varphi }}{\partial t}=-u_{\theta }(\theta )B_{\mathrm{r}}(\theta ,\varphi ,t)+\frac{D}{R}\frac{\partial B_{\mathrm{r}}(\theta ,\varphi ,t)}{\partial \theta }+S_{\varphi }(\theta ,\varphi ,t),\hfill \end{array}$$(3)

where ω(θ) is the differential angular velocity of rotation, D is the supergranular diffusivity coefficient, and u_θ(θ) is the velocity of meridional circulation. In addition, S_θ and S_ϕ represent the emergence of new flux in active regions. Each active region is inserted into the simulated radial magnetic field at the moment its center crosses the central meridian, each observed pixel replacing the corresponding pixel in the simulation. The flux imbalance of the replaced area is preserved by calculating the mean imbalance before insertion and adding it back to every pixel after insertion.

The angular velocity of differential rotation in the Carrington frame ω(θ) is:

$$\begin{array}{c}\hfill \omega (\theta )=(0.521-2.396{cos}^2(\theta )-1.787{cos}^4(\theta ))\mathrm{d}egda{y}^{-1},\end{array}$$(4)

and the latitudinal profile of the meridional circulation is:

$$\begin{array}{c}\hfill u_{\theta }(\theta )=u_{0}sin(2(\frac{\pi }{2}-\theta ))exp(\pi -2|\frac{\pi }{2}-\theta |).\end{array}$$(5)

We used a supergranular diffusion coefficient of D = 400 km² s⁻¹ and a peak meridional circulation speed of u₀ = 11 m s⁻¹. The optimization of parameters has been studied in Virtanen et al. (2017) and Whitbread et al. (2017).

It should be noted that while the speeds of meridional circulation and differential rotation are here assumed to be constant, there is evidence that they do change over time (Beck et al. 2002; Ulrich & Boyden 2005; Ulrich 2010; Ulrich & Tran 2013). This may have some effect on the results of long simulations. However, Virtanen et al. (2017) demonstrated that the surface flux transport model is relatively insensitive to such variations.

The model also uses the decay term of Baumann et al. (2006). This term models the slow temporal decay of the radial field with the following equation:

$$\begin{array}{c}\hfill E(\theta ,\varphi ,t)={\displaystyle \sum _{l=1}^{l=9}\frac{1}{{\tau }_l(k_l)}\sum _{m=-l}^{m=l}{c}_{\mathit{lm}}(t)Y_{\mathit{lm}}(\theta ,\varphi ),}\end{array}$$(6)

where E(θ, ϕ, t) is the decay rate, Y_lm(θ, ϕ) are spherical harmonics, c_lm are the harmonic coefficients of the simulated radial field at time t and τ_l are the decay times for each harmonic mode, l. Only low harmonic modes with decay times of more than one year are included, which is why the value of l does not go above l = 9 (see Table 1 in Virtanen et al. 2017).

4. Active region reconstruction

To construct active regions from the calcium K line observations, we followed the method presented in Pevtsov et al. (2016) and Virtanen et al. (2019), albeit with small improvements. The method is here described in full, but for figures and additional information, see Virtanen et al. (2019). We first identified chromospheric plages, which appear above photospheric active regions, by applying a threshold to synoptic maps of the Ca II K line intensity. In our earlier SFT simulations we used the threshold of 50G to define the active regions in NSO/KP synoptic maps of the photospheric magnetic field (Virtanen et al. 2017, 2018). To find the corresponding Ca II K line intensity threshold for chromospheric active regions of a similar size, we first computed the average number of pixels above the threshold of 50G in the NSO/KP synoptic maps of the photospheric magnetic field and then searched for a Ca II K line intensity threshold that would select a similar number of pixels from the normalized Ca II K line maps. In this comparison, we used the latitude range from −50° to 50°, thus excluding the polar regions, which are partly erroneous in the NSO/KP maps and should not contain active regions. We found 1.094 to be the optimal threshold for the identification of plages from the Ca II K line maps (for a more detailed explanation on how this value was obtained, see Virtanen et al. 2019). Using this method, the average number of pixels above the threshold is the same for the Ca II K line maps and the NSO/KP synoptic maps of the magnetic field. All pixels that reside between −50° and +50° latitude and have an intensity above the threshold are then active pixels, meaning that they are inside a chromospheric plage, directly above a photospheric active region.

While the Ca II K line maps give us the locations of plages (active regions), they do not contain information about the polarity of the magnetic field. To assign a positive or negative polarity to the active pixels we used observations of sunspot polarity and field strength from Mount Wilson Observatory. For every pixel above the threshold, we searched for a sunspot that was observed within five days of the time when the pixel crossed the central meridian and was located within a distance of five degrees in both longitude and latitude. This distance was defined to be the distance from the center of the sunspot to the pixel, minus the radius of the sunspot, which was computed from the sunspot area under the assumption that the sunspot is circular. If one or more sunspots fulfilling these criteria were found, the pixel was given the polarity of the closest sunspot. In the resulting synoptic map, the active pixels around a sunspot with a negative polarity have a negative polarity and the active pixels around a sunspot with a positive polarity have a positive polarity, giving the active pixel map a bipolar polarity structure that approximates the polarity structures of the active regions in which the sunspots emerged.

To determine the polarity of active pixels that are not located near an observed sunspot we divided the active pixels into connected regions; in this context, “connected” indicates that the sides or corners of pixels touch each other. If the size of a connected region was at least ten pixels, we used nearest neighbor interpolation on that connected region to define the polarities of the remaining connected pixels and to fill the region with polarity information. This essentially spreads the polarity information from the sunspots, which usually reside at the center of an active region, to a wider area.

If the polarity of a connected active pixel region could not be determined using the above methods, due to either the active region not containing sunspots or the polarity of the sunspots not being known, Hale’s polarity law (Hale et al. 1919; for review see Pevtsov et al. 2014 and references therein) was used to assign polarity. The leading and trailing edges of the active pixel region were assigned positive and negative polarities based on the hemisphere and the Hale polarity of the solar cycle in which the active pixels appeared. The polarity of the rest of the pixels was then determined by running nearest neighbor interpolation until the polarity of every pixel was determined. This essentially fills the region from the leading and trailing edges towards the center, retaining the expected Hale polarity. Using Hale’s polarity law allows us to start the reconstruction from 1915 when the Ca II K observations started, instead of 1917, when the sunspot polarity observations were initiated. However, the accuracy of the first two years suffers from a lack of polarity observations.

To construct active regions from the synoptic maps of active pixels, we combined the synoptic maps into one long continuous map spanning the whole studied time period. We ran a window four pixels wide in both longitude and latitude through the map and connected all active pixels that fit inside the window. We then defined all connected active pixels as belonging to the same active region. This means that not all pixels in the final active region are necessarily next to each other and, thus, small gaps of a few pixels are allowed. It should also be noted that the active regions defined here are not the same as the connected regions used above to determine the polarity of active pixels. An active region may contain multiple connected regions of active pixels.

After the active regions were thus defined, the total magnetic flux inside them was determined. To approximate the total flux within each active region, we compared the area of the active regions identified from the Ca II K line maps to the total unsigned flux contained within the same regions in the synoptic maps of the photospheric magnetic field. The integrated brightness of a plage correlates with the total flux of an active region and is also proportional to the size of the plage, allowing us to make a linear fit between the plage area and its total flux (e.g., Pevtsov et al. 2016; Tähtinen et al. 2022). Because in an even sine-latitude grid, every pixel has the same area, the total unsigned flux Φ of an active region can be written as:

$$\begin{array}{c}\hfill \mathrm{\Phi }=(B\ast N_{\mathrm{p}}-C)A_{\mathrm{p}},\end{array}$$(7)

where N_p is the number of pixels in the active region and A_p is the area of a single pixel. The values of the constants B and C are determined from the linear fit and are B = 111.02G and C = 1546.7G. If the number of pixels in the active region is so small that the total unsigned flux becomes negative, the region is discarded. This effectively sets a lower limit of 14 pixels to the area of an active region.

The distribution of flux within each reconstructed active region is taken to be the same as the intensity distribution of Ca II K line intensity within that region. This may differ from the distribution of magnetic flux in the photosphere but, as was shown in Virtanen et al. (2017), the exact distribution of flux within an active region is not important for the SFT simulation, as long as the total flux and polarity are preserved. Using the active regions with the magnetic fluxes assigned to the pixels, we then normalized the total unsigned flux to match the value given by Eq. (7). To maintain the flux balance, the normalization was performed separately for positive and negative polarities, so that there are equal amounts of positive and negative flux inside each active region. This also causes the flux distributions to slightly differ from the Ca II K line intensity distribution for positive and negative pixels. If the imbalance of positive and negative pixels was at least half of the total area of the region, meaning that there were at least three times as many positive pixels as negative, or vice versa, the whole active region was discarded.

We used the reconstructed active regions for the first simulation run, obtaining a set of simulated synoptic maps of the photospheric magnetic field for the studied time period. If a newly observed plage is above an already existing active region, its magnetic field is likely to have the same polarity as the active region. Thanks to this observation, we can use the simulated synoptic maps to improve the accuracy of polarity information in the reconstructed active regions. We went through all active pixels for which we were unable to determine a polarity from the sunspot observations and assigned them the polarity existing in the simulation – if the simulated field strength was at least 10G during the rotation the active pixel emerged, or the previous rotation. As shown in Virtanen et al. (2019), this treatment slightly improves the accuracy of the reconstructed polarity structure.

Figure 1 shows the active pixel map and the final reconstruction of Carrington rotation 1690 as an example to illustrate the reconstruction process. The reconstructed map is a snapshot of the simulation at the end of the rotation. The active regions on the left side of the map have been recently added to the simulation and have not yet decayed, while the active regions on the right side of the map have already had time to decay and spread to a wider area.

Fig. 1. Active pixels selected from the Ca II K line map of Carrington rotation 1690 (left) and reconstrucrion of the photopheric magnetic field for the same rotation (right).

Figure 2 shows the number and total unsigned magnetic flux of the reconstructed active regions for each Carrington rotation. Total unsigned flux is represented by the blue bars and number of active regions by the red bars directly below. In addition, the monthly sunspot number is shown as a thick green line on top of the total flux in the upper part of the plot. The sunspot number has been multiplied by 10²² to bring it to a similar scale with the total flux for easier comparison.

Fig. 2. Total unsigned magnetic flux in each Carrington rotation (blue bars) and number of active regions in each Carrington rotation (red bars). The monthly sunspot number, multiplied by 1022, is shown as a thick green line on top of the total flux in the upper part of the plot. The numbers above the blue bars are solar cycle numbers. Sunspot data from the World Data Center SILSO, Royal Observatory of Belgium, Brussels.

The total flux, as well as the number of active regions, follows the evolution of the sunspot number fairly closely. The sunspot cycle is clearly visible in the reconstruction and the relative height of the cycles mostly agrees with what one would expect based on the sunspot number. In fact, the total flux often follows the sunspot number even at shorter time scales. Sharp increases or decreases of the sunspot number are often accompanied by a similar change in total flux. The one cycle that clearly deviates from this pattern is cycle 19 in the late 1950s and early 1960s. While cycle 19 is the most active event based on the sunspot number, the amount of flux contained in the reconstructed active regions from the late ascending phase to the late declining phase of the cycle is higher than expected. The sunspot number has increased only moderately compared to cycle 18, but the total flux remains for several years significantly higher than in any other cycle. The number of active regions is also high but does not show an increase that would be as dramatic compared to previous cycles, which implies that the size of active regions has significantly increased. The average total emerging flux per rotation over the whole reconstruction is 1.2 × 10²³ Mx, while the average total emerging flux per rotation for cycle 19 is 2.2 × 10²³ Mx. The corresponding values for the number of emerging active regions per rotation are 8.0 and 9.8. This means that the number of active regions is only 22% above average during cycle 19, while the total emerging flux is 82% above average. Since, as shown in Eq. (7), the reconstruction assumes that the total flux of an active region is linearly dependant on the total area, the increase in total flux must be caused by a similar increase in active region area.

Bertello et al. (2020) reported an increase in plage contrast in the MWO Ca II K line observations during cycle 19, most likely caused by a change in the width of the exit slit of the spectrograph. Even though we have here used recalibrated data that should not suffer from this problem, at least to the same degree, the fact that this known contrast problem perfectly coincides with the unexplained increase in the total flux of reconstructed active regions casts some doubt on whether the unexpected increase is real or simply an artifact related to changes in the observational setup (changing the width of the spectrograph slit), especially since the reconstruction relies on thresholding the Ca II K line maps. We used a constant threshold, which means that changes in contrast can affect the area, number, and total flux of reconstructed active regions. For these reasons, the amount of emerging flux could be too large during cycle 19 compared to other cycles.

On the other hand, a significant number of very large sunspots are present in the MWO sunspot drawings in 1957 and 1959, which supports the idea that the enhanced activity is real and not an artifact. This uncertainty, as well as the issues related to the data gaps discussed below, could be resolved in a future study by using other Ca II K line observations to supplement the MWO maps.

While the data coverage is generally good, there are two gaps of significant length that are visible in Fig. 2. The longest gap is in 1967 during the ascending phase of cycle 20, and a shorter gap can be seen in 1939 in the descending phase of cycle 17. The gap in 1939 can be traced to Ca II K line observations at MWO, as evident from Fig. 2 in Bertello et al. (2020), and the longer one is caused by errors in the sunspot database.

5. Simulated photospheric magnetic field

The upper panel of Fig. 3 shows a simulated supersynoptic map from 1915 to 1985, covering solar cycles 15–21. The simulation uses the reconstructed active regions as input. Because we have no information about the strength of the magnetic field prior to 1917, the simulation starts from an empty map, which means that there are no polar fields at the beginning of the simulation, and it takes a few years for them to form. It was shown in Virtanen et al. (2017) that the initial field does not significantly affect the simulation after about one cycle, so we expect the lack of initial polar fields to cause uncertainties during cycle 15, especially at high latitudes, but afterwards the effect is insignificant.

Fig. 3. Supersynoptic map of the reconstructed photospheric magnetic field from 1915 to 1985 (top). Supersynoptic map of the NSO/KP synoptic maps of the photospheric magnetic field from 1975 to 1985 (bottom).

We can see that the reconstruction creates a realistic solar cycle evolution. The active regions injected to the activity belts form poleward flux surges that create the polar fields and change their polarity once every cycle. The strength of poleward surges and polar fields steadily increases from cycle 15 to cycle 19, and then again decreases. The unusually large flux emergence during cycle 19 leads to extremely strong poleward surges and polar fields. Cycle 20, on the other hand, is very weak. This is probably mainly caused by the long data gap in the ascending phase of the cycle, which causes a large amount of flux to be missing from the simulation.

The lower panel of Fig. 3 shows the NSO/KP observations for comparison. The observations start in 1975, so they overlap with the reconstruction only during cycle 21. The evolution of poleward surges is fairly similar. The polar fields cannot be directly compared due to obvious errors in the observations. For a more detailed comparison between the reconstruction and observations, see Virtanen et al. (2019).

Figure 4 shows the average strength of polar fields above +60° and below −60° latitude. The solid lines are from the reconstruction and the dashed lines are from NSO/KP synoptic maps. As can be seen from unrealistically large variations of the observed values in the late 1970s, the NSO/KP data is partly erroneous especially in the northern pole. As discussed in Sect. 2, the polar fields suffer from known problems related to poor data quality and polar field filling.

Fig. 4. Average strength of the photospheric magnetic field above 60 deg and below −60 deg latitude. Solid red and blue lines are the northern and southern hemispheres of the reconstruction, respectively. Dashed orange and light blue lines are the northern and southern hemispheres of the NSO/KP synoptic maps, respectively.

Starting from unknown polar field, it takes a few years for the simulation to develop realistic polar fields (Virtanen et al. 2017), which is why the field is relatively weak until the early 1920s. Cycle 16 contains a sudden weakening of both polar fields caused by poleward surges of anti-Hale polarity flux around 1930; specifically, “anti-Hale” indicates that the surge consists of opposite polarity to that predicted by Hale’s polarity law. These anti-Hale surges are also visible in Fig. 3. During cycle 17, there is a large hemispheric asymmetry present. The southern hemisphere is activated first and the southern polar field reverses its polarity earlier than the northern polar field. Still, the northern pole obtains a larger maximum strength than the south pole. During cycle 18, the northern polar field reverses polarity first and also grows to be much stronger than the southern polar field due to a few strong poleward surges.

Polar fields in the late declining phase of cycle 19 are extremely strong, especially in the northern hemisphere. These strong polar fields pertain for a very long time, until the maximum of cycle 20, and near the maximum of cycle 20 the polar fields fluctuate around zero for an extended period of time. However, these features are most likely caused by the gap in the active region reconstruction in the ascending phase of the cycle, which halts flux emergence during the gap, delaying the cancelation of old polarity fields and weakening the polar fields even after the polarity reversal.

6. Coronal field

Reconstructed synoptic maps of magnetic field described in Sect. 5 enables the modeling of coronal magnetic fields and the solar wind for historical periods using modern models. Here, we apply the reconstructed magnetic fields to the modeling of solar corona. The potential field source surface (PFSS) model can be used to model the coronal magnetic field from the photospheric magnetic field. It was developed by Schatten et al. (1969) and Altschuler & Newkirk (1969), and later refined by Hoeksema (1984) and Wang & Sheeley (1992). The potential field model includes the assumption that in the corona there are no electric currents affecting the structure of the magnetic field. The lower boundary of the model is the photosphere, and the upper boundary is the so called source surface, where all magnetic field lines are assumed to be radial. Here, we use the reconstructed synoptic maps of the photospheric magnetic field as the lower boundary data and track the path of the field lines to the source surface.

The distance between the photosphere and the source surface is the only physical parameter of the model. The optimal value has been shown to be around 2.5 solar radii (Hoeksema 1984), but it varies over time and depends on the phase of the solar cycle (Koskela et al. 2017). In this work, we apply the most common constant value of 2.5 solar radii.

In Fig. 5, the amount of open field lines at different latitude ranges is visualized. We used the PFSS model to track the magnetic field lines originating from each pixel of reconstructed synoptic maps of the photospheric magnetic field and we then determined whether the field line is open or closed. A field line is open if it goes through the source surface and closed if it turns back before reaching the source surface, forming a closed loop. Each synoptic map was divided into nine 20 deg latitude bands and the number of pixels corresponding to open and closed field lines in each band counted. The nine panels in Fig. 5 show the ratio of pixels with open field lines and all pixels in each synoptic map as a function of time for each latitude band. This is essentially an estimate of the photospheric surface area covered by coronal holes at each latitude band. For comparison, the same ratios were also computed from the NSO/KP synoptic maps and depicted in Fig. 5 as a red line between 1975 and 1985.

Fig. 5. Fraction of photospheric surface area covered by open field lines in nine 20 deg latitude bands. The limits of each band are shown above the panels. Blue lines are from the reconstruction, red lines from NSO/KP maps.

At high latitudes, the reconstruction shows a very consistent solar cycle evolution. At minimum times, both poles are are covered by coronal holes, and at maximum times, when the polarities reverse, there are no coronal holes at the poles. At lower latitudes, there is a lot more variance, but the expected solar cycle behavior, where low-latitude coronal holes are most common in the descending phase of the cycle and disappear in minima, is clearly visible. One interesting feature is the peak in low- and mid-latitude coronal holes in the early 1930s, especially at −10° − +10° in 1930, which is responsible for the abrupt peak of solstice-time geomagnetic activity seen in Fig. 2 of Mursula et al. (2017).

At low latitudes, the reconstruction and the observations follow a similar evolution, but the reconstruction contains more open field lines. Due to the simple nature of the model, the reconstruction lacks small-scale features of the photospheric magnetic field and tends to create a smoother field than what is seen in the observations. Larger areas of a single polarity then lead to somewhat larger areas of open field lines. It is also known that the PFSS model tends to exaggerate coronal hole areas. Some similarity between the observations and the reconstruction can also be seen at high latitudes, although the problems caused by erroneous polar fields in the observations are evident.

At short time scales, the reconstruction often deviates significantly from observations. This was also found by Virtanen et al. (2019). Errors in individual reconstructed active regions have an effect on the whole reconstructed magnetic field for a short time, but over long time scales the errors are smoothed out.

7. Comparison with eclipse drawings

Total solar eclipses offer a unique opportunity to observe the momentary structure of the coronal magnetic field, even with the naked eye. When the bright solar disk is hidden, the differences in the coronal brightness can reveal the locations of coronal holes as well as streamers and the closed field line loops beneath them. Solar eclipses have been observed and recorded for thousands of years, and drawings, and eventually photographs, of the corona exist for numerous historic eclipses.

The accuracy of SFT reconstruction between 1975 and 1985 has previously been studied in Virtanen et al. (2019) and, to a lesser extent, in Sect. 5 of this paper. Here, we compare the reconstructed coronal magnetic field with eclipse drawings to assess the accuracy at times before direct magnetic field observations became possible. We have chosen drawings of two eclipses, the eclipse on September 21, 1922, and the eclipse on September 10, 1923. Early eclipse photographs often do not show the structure of the corona very clearly, which is why drawings were selected for this comparison. The recording of the inner and outer corona required to take separate photographic images with different exposures. The human eye, however, perceives light in a logarithmic manner and, thus, drawings may be better for representing both the inner and outer corona in a single image. Additionally, the orientation of the Sun in the photographs is often not documented. While the selected eclipses occurred soon after the beginning of the reconstruction, meaning that the reconstructed polar fields had not yet reached their full strength, we expect the large-scale structure of the coronal field to be realistic. At this point, the lack of initial field should not affect the reconstruction at low- or mid-latitudes, and at high latitudes the polarity of the field should be correct, even though the strength might not have reached its correct value due to the slowness of poleward transport of flux. In the PFSS model, this would mostly cause errors in the size of the polar coronal holes – and not the locations of streamers.

To carry out these comparisons, we must first determine how the reconstructed coronal magnetic field looks from the Earth during an eclipse. Areas with a high number of closed field lines tend to be bright due to the plasma trapped in magnetic loops, whereas coronal holes are darker. We assume here that coronal brightness depends on the number of closed field lines along the line of sight. We take a cube containing the Sun and divide it into an evenly spaced three-dimensional grid in Cartesian heliocentric coordinates, where the z-axis points toward the north, and the x- and y-axes are in the solar equatorial plane. A cube is used purely because of computational simplicity. In the PFSS model all field lines outside the source surface are open, so the shape of the grid does not matter as long as the whole source surface is contained within it. We rotate the grid so that x-axis points directly towards the Earth, taking into account the location of the central meridian and the angle between the orbit of the Earth and the equator of the Sun at the time of the eclipse. Going through the grid in the direction of the x-axis then gives an evenly spaced grid along the line of sight. We use the PFSS model to determine whether the field line in each point of the grid is open or closed, and sum the number of closed field lines over the x-axis. This results in a two-dimensional (2D) matrix resembling an eclipse photo. The cells of the matrix correspond to pixels in a photo, and the value of each cell is proportional to the number of closed field lines along the line of sight and, consequently, the brightness in the corresponding pixel. It should be noted that this is a highly simplified method that does not take into account Thomson-scattering or other physical processes that might affect coronal brightness as it is seen from the Earth.

The top panel of Fig. 6 shows the reconstructed corona during the 1922 eclipse. The bottom panel shows a drawing of the same eclipse. We can see that the streamers at top right and bottom left are clearly visible in both, the drawing and the reconstruction. The third streamer seen in the top left of the drawing is also present in the reconstruction, although it is much weaker. Coronal holes are seen at both poles, and the southern coronal hole is significantly larger in both the drawing and the reconstruction. For comparison, Fig. 7 shows a photo of the eclipse. The exact orientation of the Sun in the photo is not known, but the photo has been rotated 150 deg counter-clockwise to roughly match the drawing in Fig. 6. The streamers seen in the drawing and the reconstruction can also be seen in the photo, although not as clearly.

Fig. 6. Reconstruction (top) and drawing (bottom) of the solar corona during the eclipse on September 21, 1922. The drawing is from the collection of the Kislovodsk Solar Station (http://old.solarstation.ru).

Fig. 7. Photo of the solar eclipse on September 21, 1922 from the Hamburger Sternwarte Digital Plate Archives (https://plate-archive.hs.uni-hamburg.de).

Figure 8 shows a similar comparison of the reconstruction and a drawing for the eclipse of 1923. The agreement is not as good as in Fig. 6, but many similarities still exist. In the drawing there are two streamers in the right side, but in the reconstruction only one is clearly visible (to the right at about three o’clock using clock bearings). A faint second coronal streamer can be seen in the lower part (between four and five o’clock), but the direction is not exactly the same as in the drawing. It is worth noticing that in the “Solar Eclipse, Lompoc, California, 1923” painting by American artist Howard Russell Butler¹, this streamer is shown farther to the South as compared with drawing shown in Fig. 8, which is in a better agreement with our model. In 1918, H.R. Butler was invited to draw a solar eclipse as part of expedition organized by the US Naval Observatory, and has developed a technique for capturing rapid changes such as during the eclipses. Ultimately, Butler painted a view of the solar corona during the 1918, 1923, 1925, and 1932 eclipses. On the left side, there is one wide streamer in both the drawing and the reconstruction, but in the reconstruction it is directed more towards the north (about 10 o’clock), while in the drawing (and Butler’s painting), it points directly to the East (9 o’clock).

Fig. 8. Reconstruction (top) and drawing (bottom) of the solar corona during the eclipse on September 10, 1923. The drawing is from the collection of the Kislovodsk Solar Station (http://old.solarstation.ru).

In both cases the reconstruction provides a large-scale structure that is fairly similar to the drawing, although the agreement is clearly better in 1922. This provides further evidence that the reconstructed active regions and the SFT simulation based on them creates a realistic magnetic field with large-scale structure that roughly agrees with the observations.

As our final comparison, Fig. 9 shows reconstruction of coronal magnetic fields using PFSS model for the total solar eclipse 11 June 1983. Since NSO/KP magnetograms were available for this eclipse, we can compare the reconstructions using the pseudo and real magnetograms. Our pseudo-magnetograms show a reasonably good agreement with the reconstruction based on real magnetograms. Small differences in location of structures could be contributed to the fact that the reconstruction using pseudo-magnetograms was centered on the day of eclipse, while for real magnetograms, we used full integer solar rotation not centered on the day of eclipse. Still, most major coronal streamers have the magnetic structures in PFSS model with the exception of CS5, which is not reproduced by either pseudo or real magnetograms.

Fig. 9. PFSS model of solar coronal fields for solar eclipse 11 June 1983 derived using (a) pseudo-magnetograms and (b) real magnetograms from NSO/KP, along with (c) their comparison with eclipse observations. Black arrows point to coronal streamers (CS) and the corresponding structures in PFSS model. Red blue lines show open magnetic fields with positive and negative polarity footpoints. Green lines are closed field lines.

8. Conclusions

We reconstructed the photospheric magnetic field from 1915 to 1985 using the Ca II K line and sunspot polarity observations, following the methods presented and tested in Virtanen et al. (2019). The number and total flux of the reconstructed active regions mostly follow the solar cycle evolution of the monthly sunspot number. An exception is cycle 19, which contains an unexpectedly large amount of flux, implying that the size of the reconstructed active regions is not only larger than during other cycles, but also relatively larger than expected based on the sunspot number. This may be caused by enhanced activity during this very strong cycle, but there is also a possibility it is related to the change in the contrast of the MWO Ca II K line maps. This change was found by Bertello et al. (2020) and it occurs at the same time as the rise in the total flux of the reconstruction.

The reconstructed photospheric fields allow the application of modern models to simulate and study coronal magnetic fields and the solar wind for historical periods, which predate direct observations of magnetic fields (magnetograms). Here, we demonstrate the application of this new dataset to represent solar corona and past solar eclipses.

We used the reconstructed active regions and an SFT model to simulate the evolution of the entire photospheric magnetic field. This resulted in a supersynoptic map with a realistic solar cycle and polarity reversals. Due to missing initial polar fields, the first few years of the simulation are inaccurate and data gaps during cycle 20 and, to a lesser extent during cycle 17, lead to some missing flux and weakened polar fields. The unusually large total flux of the reconstructed active regions during cycle 19 cause a very strong cycle also in the simulation, including very strong polar fields. The strong polar fields also persist for an extended period of time due to the data gap during the ascending phase of cycle 20. Between 1975 and 1985 the reconstructed field agrees with the field observed at NSO/KP fairly well at low- and mid-latitudes. At high latitudes, the errors in observations make reliable comparisons practically impossible.

We also calculated the coronal magnetic field from the reconstructed photospheric magnetic field using the PFSS model with a source surface distance of 2.5 solar radii. We found that the reconstruction tends to contain a greater number of open-field regions (coronal holes) than the observations, especially at low latitudes. Despite this difference, the area covered by open field lines evolves with the solar cycle as expected. At solar minimum times the poles are covered with coronal holes, and later in the cycle the number of low-latitude coronal holes increases and peaks in the descending phase of the cycle. The reconstruction also roughly follows the evolution seen in the observations during the overlapping interval from 1975 to 1985.

Finally, we compared the reconstructed coronal field with two eclipse drawings from 1922 and 1923 to assess the accuracy of the early part of the reconstruction. For the eclipse in 1922, the agreement is very good. The reconstruction and the drawing depict a very similar structure of coronal holes and streamers. For the eclipse of 1923 the agreement is not as good, but many similarities can still be seen.

The presented results show that the reconstruction contains quite a realistic evolution of the solar magnetic field from 1915 to 1985 that agrees reasonably well with observations. The large-scale structure of the reconstructed field is similar to observed field from 1975 onwards and to eclipse drawings in the 1920s. A few problems can also be identified, starting with the fact that at short timescales, the reconstruction can deviate from observations considerably due to errors in the reconstruction of a few individual active regions (Virtanen et al. 2019). Cycle 19 clearly differs from the other cycles in terms of flux emergence, which may be caused by the record high activity during the cycle, but the possibility of it being related to changes in the contrast of the Ca II K line maps cannot be ruled out. Data gaps also cause uncertainties in some cycles. These problems could potentially be resolved in the future by using other Ca II K line observations to supplement the MWO Ca II K maps. Overall, we conclude that the large scale photospheric magnetic field derived using our method can be successfully used as input to modern models to represent the 3D distribution of magnetic field in solar corona over the last 100 yr.

Acknowledgments

Authors are members of the international team on Modeling Space Weather and Total Solar Irradiance over the Past Century supported by the International Space Science Institute (ISSI), Bern, Switzerland and ISSI-Beijing, PRC. The National Solar Observatory (NSO) is operated by the Association of Universities for Research in Astronomy (AURA), Inc., under cooperative agreement with the National Science Foundation. We thank staff of Kislovodsk Mountain Astronomical Station of the Pulkovo observatory for providing access to archive of solar eclipse drawings.

References

All Figures

	Fig. 1. Active pixels selected from the Ca II K line map of Carrington rotation 1690 (left) and reconstrucrion of the photopheric magnetic field for the same rotation (right).
In the text

Fig. 2. Total unsigned magnetic flux in each Carrington rotation (blue bars) and number of active regions in each Carrington rotation (red bars). The monthly sunspot number, multiplied by 1022, is shown as a thick green line on top of the total flux in the upper part of the plot. The numbers above the blue bars are solar cycle numbers. Sunspot data from the World Data Center SILSO, Royal Observatory of Belgium, Brussels.

In the text

	Fig. 3. Supersynoptic map of the reconstructed photospheric magnetic field from 1915 to 1985 (top). Supersynoptic map of the NSO/KP synoptic maps of the photospheric magnetic field from 1975 to 1985 (bottom).
In the text

	Fig. 4. Average strength of the photospheric magnetic field above 60 deg and below −60 deg latitude. Solid red and blue lines are the northern and southern hemispheres of the reconstruction, respectively. Dashed orange and light blue lines are the northern and southern hemispheres of the NSO/KP synoptic maps, respectively.
In the text

	Fig. 5. Fraction of photospheric surface area covered by open field lines in nine 20 deg latitude bands. The limits of each band are shown above the panels. Blue lines are from the reconstruction, red lines from NSO/KP maps.
In the text

	Fig. 6. Reconstruction (top) and drawing (bottom) of the solar corona during the eclipse on September 21, 1922. The drawing is from the collection of the Kislovodsk Solar Station (http://old.solarstation.ru).
In the text

	Fig. 7. Photo of the solar eclipse on September 21, 1922 from the Hamburger Sternwarte Digital Plate Archives (https://plate-archive.hs.uni-hamburg.de).
In the text

	Fig. 8. Reconstruction (top) and drawing (bottom) of the solar corona during the eclipse on September 10, 1923. The drawing is from the collection of the Kislovodsk Solar Station (http://old.solarstation.ru).
In the text

Fig. 9. PFSS model of solar coronal fields for solar eclipse 11 June 1983 derived using (a) pseudo-magnetograms and (b) real magnetograms from NSO/KP, along with (c) their comparison with eclipse observations. Black arrows point to coronal streamers (CS) and the corresponding structures in PFSS model. Red blue lines show open magnetic fields with positive and negative polarity footpoints. Green lines are closed field lines.

In the text