11institutetext:  Aalto University Metsähovi Radio Observatory (MRO),
Metsähovintie 114, 02540 Kylmälä, Finland 22institutetext: Department of
Computer Science, Aalto University, PO Box 15400, 00076, Finland
22email<EMAIL_ADDRESS>33institutetext: Nordic Institute for
Theoretical Physics, Roslagstullsbacken 23, 106 91 Stockholm, Sweden
44institutetext: Max Planck Institute for Solar System Research, Justus-von-
Liebig-Weg 3, 37077 Göttingen, Germany 55institutetext: Stockholm University,
SE-106 91 Stockholm, Sweden 66institutetext: KTH Royal Institute of
Technology, SE-100 44 Stockholm, Sweden 77institutetext: School of
Mathematics, Statistics and Physics, Newcastle University, NE1 7RU, UK

# Almost fifty years of Metsähovi solar observations on 37 GHz with recovered
digitised historical maps.

Sami Kivistö 0000-0002-7693-6468 1122 Frederick Gent0000-0002-1331-226
2233556677 Merja Tornikoski 0000-0003-1249-6026 11 Joni
Tammi0000-0002-9164-2695 11
Maarit Korpi-Lagg0000-0002-9614-2200 223344

(Received ; accepted TBA)

###### Abstract

Context.  Aalto University Metsähovi Radio Observatory has collected solar
intensity maps for over $45$ years. Most data coverage is on the
$37\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ frequency band, tracking emissions
primarily at the chromosphere and coronal transition region. The data spans
four sunspot cycles or two solar magnetic cycles.

Aims.  We present solar maps, including recently restored data prior to 1989,
spanning 1978 to 2020 after correcting for observational and temporal bias.

Methods.

The solar maps consist of radio intensity sampled along scanlines of the
antenna sweep. We fit a circular disk to the set of intensity samples,
neglecting any exceptional features in the fitting process to improve
accuracy. Applying a simple astronomical model of Sun and Earth, we assign
each radio specimen its heliographic coordinates at the time of observation.
We bin the sample data by time and heliographic latitude to construct a
diagram analoguous to the classic butterfly diagram of sunspot activity.

Results.  Radio butterfly diagram at
$37\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$, spanning solar cycles 21 to 24
and extending near to the poles.

Conclusions.  We have developed a method for compensating for seasonal and
atmospheric bias in the radio data, as well as correcting for the effects of
limb brightening and beamwidth convolution to isolate physical features. Our
observations are consistent with observations in nearby bandwidths and
indicate the possibility of polar cyclic behaviour with a period exceeding the
solar 11 year cycle.

###### Key Words.:

Solar physics, radioastronomy, butterfly diagram, solar cycle

##  1 Introduction

The solar magnetic field is known to switch polarity in a full cycle of
approximately 22 years, which consists of two sunspot cycles (Hale et al.
1919, Babcock 1961). At the beginning of each sunspot cycle, the solar
magnetic field is an axial dipole. Solar differential rotation transforms the
predominantly axisymmetric poloidal magnetic field into an azimuthal field
with opposite polarity between hemispheres. In the presence of convection and
the Coriolis force due to solar rotation, this azimuthal field is further
distorted by the action of turbulent vortical flows stretching, twisting and
looping the field, known as the $\alpha$-effect (Steenbeck et al. 1966, Krause
& Raedler 1980). After approximately 11 years the large scale magnetic field
evolves into a dipolar field with polarity reversed from its original sign and
the next sunspot cycle begins. The full 22 year cycle comprises two such
sunspot cycles.

Starting from a quiet Sun containing a relatively weak axial dipole, each
sunspot cycle emerges as activity on the visible surface of the Sun, typically
near mid-latitudes ($\pm 20^{\circ}$ to $\pm 30^{\circ}$). While the frequency
and coverage of the sunspots increases towards a maximum, the location of the
activity drifts equatorward. This migration continues after the maximum with
decreasing activity closer to the equator, bringing the sunspot cycle to an
end. This was first recorded in the Maunder butterfly diagram by Edward
Maunder in 1904. Since sunspots are a feature of Sun’s magnetic field, and as
the solar magnetic field is mostly axisymmetric, the longitudinal location of
the activity can reasonably be neglected (Pelt et al. 2006). Considering this
two-dimensional time-latitude histogram of sunspot activity or magnetic field
intensity, similar butterfly diagrams are also recovered from solar
magnetograms (Gaizauskas et al. 1983, Vecchio et al. 2012, Leussu et al.
2017). It is also reported for radio brightenings of the 17 and
$34\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ band (Shibasaki 2013, Selhorst et
al. 2014) and of $37\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ (Kallunki et al.
2018).

From sunspots and magnetograms of the photosphere we have relatively rich
well-resolved data about the two-dimensional (2D) evolution of the solar
surface magnetic field, but observation of the three-dimensional (3D)
structure of the field is limited by the weakness of magnetic measurements
from the solar atmosphere relative to the surface. Metsähovi Radio Observatory
(MRO)111 https://www.aalto.fi/en/metsahovi-radio-observatory is observing
emissions, primarily in $37\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ frequency
band, dating back over four decades to the late 1970s. This is in a frequency
range suitable for detecting non-thermal emissions from high energy electrons
associated with active-region magnetic fields, so is a good proxy for the
intensity of the magnetic field. A detailed analysis by Kallunki et al. (2020)
relate brightening at $37\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ emissions
with coronal holes and magnetic bright points in extreme UV emissions. These
are also well suited to observing the thermal emissions from the quiet Sun.
Comparing and contrasting the maps from these emissions with other
observations can help to reveal more about the structure of the magnetic field
in plasma of varying temperature and by inference height. Furthermore, the
data incorporate almost four magnetic cycles, yielding a valuable resource to
investigate long-term variation within and between cycles, alongside and in
contrast to other sources. Similar bands at
$950\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ and
$84\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ have probed the chromesphere
(Loukitcheva et al. 2019), but without the long term historical record
available here.

Previous works (Kallunki et al. 2018) have studied the statistics of active
regions, radio samples with excess of a fixed threshold, sampling individual
MRO solar maps of the cycle 24. Analysis based on manually collected values
has been made, such as Kallunki et al. (2012), Kallunki & Riehokainen (2012)
and Valtaoja et al. (1987). Cycle 21 maps have previously been analysed with
manual methods in Kallunki et al. (2012) and Kallunki & Riehokainen (2012)
using undigitised maps.

For the first time, the whole 45 years of data are available in a form that is
fully accessible to automation, allowing detailed analysis and in accurate
heliographic coordinates (Metsähovi Radio Observatory 2019). The earliest of
these maps feature the solar cycle 21 and have not been available until
appropriate digital processing tools were implemented in Kivistö (2018). The
new analysis includes, along with further development of the calibration
process, calculating temporal averages of radio samples within a given region
on the rotating heliographic surface. Heliographic features can be
distinguished based on their actual separateness.

Across the lifetime of the MRO the equipment has experienced several upgrades
and evolution of the data format. This article shall describe the pipeline of
data processing to transform the raw data into a time series of solar radio
maps, normalised, corrected for instrumental bias, atmospheric noise and
effects, time of day and seasonal effects and solar limb dimming effects. We
present two methods for correctly detecting the position of the solar disk in
each map and three methods for normalising the signal levels, such that the
background _black sky_ has value zero and the Quiet Sun Level (hereafter QSL),
the solar intensity in the absence of surface magnetic activity, has value
unity. The maps are then projected from Right Ascension (RA) and Declination
(Dec) onto their rotating heliographic surface coordinates. This is described
in Sect. 2. We we have also prepared tools for automatic tracking of active
regions, which are suitable for any magnetograms or similar maps, and which we
shall apply in detail in future publication.

We present a time-latitude histogram of the radio intensity, including
previously inaccessible machine-drawn maps from 1978 – 1987 and up to 2024.
This yields the most complete time series thus far available from MRO,
spanning almost completely the last four sunspot cycles. We present our
butterfly diagrams in Sect. 3 and discuss in Sect. 4 how MRO observations
compare with the results of other observations, including Shibasaki (2013).

(a)(b)(c)

Figure 1:  Typical Metsähovi solar observations on
$37\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ (a) 25.07.2011 and (b) 05.07.2016.
Prior to May 2015 the observation beam scans along the equator, with samples
mapped onto a rectangular grid. Since June 2015 a more accurate fit to the
beam path was used, scanning horizontally. For the intensity color map each
pixel is mapped to the nearest observed sample. This is an original
uncalibrated map. (c) Closeup of the 05.07.2016 map which illustrates circle
fitting.

##  2 Methods

Here, we detail the techniques applied to collate Metsähovi observations up to
present, spanning four solar cycles. The data include multiple formats,
cadence and observational conditions, so require normalisation to compare
results over time. To accurately relate observed features to the solar
surface, the maps must be centered and the chromospheric radius of the
observations accurately determined. The telescope has measurement biases and
the measured intensities are affected by the altitude of the observation or
location on the solar disk. We apply corrections to better identify features
near the limb.

###  2.1 The solar maps

During a typical measurement, the antenna beam makes a sweep across the
visible solar disk (Fig. 1). The historical record comprises three major
convention sets applying to the data.

  1. A

Since June 2015 the sweeps have been horizontal, with shifting altitude
between sweeps, as shown in Fig. 1(b). These modern maps respect the true path
of the antenna, typically recorded with $25$ samples per second for two
minutes.

  2. B

Until May 2015 the maps were scanned in the equatorial direction, so that
declination is shifted between adjacent sweeps, as displayed in Fig. 1(a). The
maps assume linear antenna tracks with constant speed and sample rate.

  3. C

Solar maps from 1978 to 1987 were originally recorded on magnetic tapes, from
which the data were rendered onto paper contour maps using a mechanical
plotter. Subsequently the magnetic tapes were lost, so that only the printed
maps were available. Using pattern recognition techniques of Kivistö (2018),
scans of the maps were converted into a digital format, which is compatible
with maps of the categories A and B. The solar disk position as well as the
rotational axis was already marked on these paper maps using a compass and a
ruler.

Each line in a digital solar map file contains the time and coordinates of a
sample as well as the recorded intensity value. The coordinates are relative
to the expected location of the center of the Sun on the sky. The intensity is
measured on $37\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ and contains linear
and logarithmic output from the amplifier as digital values. In this paper, we
focus on the linear scale. However, the logarithmic scale is more sensitive
for very low intensities and also avoids saturation during flare events, thus
offering prospects for future research.

We use the word _calibration_ for describing the combined effort of
normalising and centering a map. The applicable methods are described in Sect.
2.2. The raw intensities of a scan have an arbitrary scale, as in Fig. 1(a-c).
The signal levels must be normalised and also corrected for any minor offsets
in the sample coordinates in order to have a properly centered map. From
normalisation we obtain zero intensity at the background and unit intensity
for the QSL. For the contour maps up to year 1987, we trust the original
compass and ruler markings and no additional calibration is considered.

###  2.2 Solar map centering and normalisation

Solar maps are generated by sweeping the beam of the radio telescope at the
proximity of the solar disk while recording signal intensity. In order to
extract heliographic information from these maps, we first need to superimpose
the visible solar disk over these samples. After correctly locating the disk,
each antenna sample is accurately projected onto the rotating heliographic
surface. The disk position is located by analysing the intensity distribution
on the scans. The size of the visible disk is based on the distance between
the Sun and the observer (Earth) at the time of observation. We assume that
the radius of the Sun is $717.6\text{\,}\mathrm{M}\mathrm{m}$. This value is
based on iterating the limb brightning model using the
$37\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ MRO solar maps between 16th June,
2014 and 14th September, 2022, which are the three most recently completed
quarters of the solar cycle. We define the boundary of the heliographic sphere
to be where its projection on the solar map has intensity $0.5$ QSL.

Each solar map $S$ consists of $N_{S}$ samples $s_{i}$:

$\displaystyle S$ $\displaystyle=$
$\displaystyle\left\\{s_{i}\;:\;i=1,2,...,N_{S}\right\\}\,\text{ where}$
$\displaystyle s_{i}$ $\displaystyle=$
$\displaystyle\left(t_{i},x_{i},y_{i},u_{i},...\right)\text{.}$ (1)

Each sample $s_{i}$ contains the time of observation, $t_{i}$, and location on
the map, $(x_{i},y_{i})$, as well as the recorded intensity $u_{i}$. During
the calibration process described below, additional derived quantities are
attributed to each tuple $s_{i}$. These are denoted as $...$ on Eq. 2.2. The
calibration is iterative, with symbol $c$ used for each step, starting with
$c=0$. Each iterative step $c$ produces a set of coordinates
$(x_{i}^{(c)},y_{i}^{(c)})$ and a calibrated intensity value $v_{i}^{(c)}$,
which equals zero at the background and unity at the average interior disk
brightness, which will be defined below. We denote the radial distance from
the center of the solar disk as:

$\displaystyle r_{i}^{(c)}:=\sqrt{x_{i}^{(c)2}+y_{i}^{(c)2}}\text{.}$ (2)

Location is specified in equatorial coordinates relative to the expected
position, $({\mathsf{RA}}_{\astrosun},{\mathsf{Dec}}_{\astrosun})$, of the
central point of Sun at time $t_{i}$. We do not rely on the targeting accuracy
of the antenna, and instead determine the center of the visual disk from each
map $S$, from which we define the relative location of each $s_{i}$.

With right ascension cosine corrected, such that maps have unit aspect ratio,
we get the final coordinates (for large $c$) as:

$\displaystyle
x_{i}^{(c)}=\frac{{\mathsf{RA}}(s_{i})-{\mathsf{RA}}_{\astrosun}(t_{i})}{\cos\left({\mathsf{Dec}}_{\astrosun}(t_{i})\right)}\text{,}\quad
y_{i}^{(c)}={\mathsf{Dec}}(s_{i})-{\mathsf{Dec}}_{\astrosun}(t_{i})\text{.}$
(3)

The intensity is high at the solar disk, and it drops by orders of magnitude
as we transition to the background sky. A number of factors affect the shape
of this transition:

  1. 1.

Due to the shape and width of the beam used by the MRO $14\text{\,}\mathrm{m}$
disk, near the limb convolution with the black sky reduces the observed solar
disk intensity.

  2. 2.

The optical depth of the chromosphere varies with altitude due to increasing
temperature (Mihalas 1978, Alissandrakis et al. 2017). Nearer the limb, where
the observed line of sight makes a low angle through the corona, we would
expect increased intensity due to the higher temperatures compared to the
center of the disk. We refer to this effect as limb brightening.

  3. 3.

Observed intensity at $37\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$, as a time
average over durations comparable to solar cycle length and butterfly diagram
features (wings), can vary by heliographic latitude.

The latitude dependence, which varies by solar cycle phase, is weak relative
to beam convolution and limb brightening effects. For map calibration
purposes, we ignore these latitude dependency effects. The importance of this
approximation is a subject for future study. We have, however, implemented
tools (see Appendix A) for taking into account the latitude dependency in
order to measure certain chromosphere parameters. For this manuscript we
assume a radially symmetric antenna beam and thus the profile to correct for
limb brightening is also radially symmetric.

For each map, $S$, the raw signal intensity values, $u_{i}$, that come from
the amplifier and analog-to-digital conversion, have an arbitrary scale. These
require a normalisation to extract two representative intensity values,
${u}_{\mathrm{background}}$ and ${u}_{\mathrm{disk}}$. We have trialed three
alternate methods for this initial step:

  1. 1.

S-curve method, described orally by Juha Aatrokoski at MRO and used in MRO
maps such as in Kallunki et al. (2018).

  2. 2.

A variant of a convex hull method described in Kivistö (2018).

  3. 3.

Selecting peak values in intensity histogram, a method adopted from Shibasaki
(2013).

Here we have applied method 1, which we shall describe in detail. The method
for the initial step is not crucial since the final result is derived by
series of converging iterations that improve the centering and normalisation.
We do not rely only on the initial step. When there are bright or dim regions
at the boundary, we discard them as outliers and limit the circle fitting to
the remaining points. While deriving the limb profile, as displayed in
Fig.3(b), we can improve the estimate of a sample’s radial distance based on
its intensity, as illustrated in Fig.2.2(c).

A through analysis and comparison of different calibration methods is a
possible subject of further research. Following method 1, the samples
$s_{i}\in S$ are sorted by absolute intensity. An example of a resulting
intensity vs. sample index is plotted in Fig. 2(a). We observe two vertical
regions in the plot, since the dark background as well as the bright solar
disk cover most of the samples. Between these regions, bounded by the dashed
lines, the curvature of the plot changes sign. Here we place the pivot point
of the normalisation (Fig. 2)(b). The median intensity of samples with
intensity higher than the pivot is taken as representative of the inner disk.
The lowest intensity tertile is taken as representative of the background.
This fixes the signal range from background to QSL as the starting point of
further iterations, to identify the disk centre, radius and intensity
normalisation.

(a)(b)

Figure 2:  Radio samples of a typical solar scan sorted according to
intensity. Panel (a) displays all samples, within the limb transition between
background and disk located with the dashed lines. A zoom in between the
dashed line in (b) is fitted with a polynomial, whose inflection (cross) is a
pivot point for signal level normalisation.

We then perform a linear transformation for the raw intensity values $u_{i}$
and obtain the first normalised intensities $v_{i}^{(0)}$ as

$\displaystyle
v_{i}^{(0)}:=\frac{u_{i}-{u}_{\mathrm{background}}}{{u}_{\mathrm{disk}}-{u}_{\mathrm{background}}}\;\text{with}\;s_{i}\in
S\text{.}$ (4)

The normalised intensities are then used for calculation of the initial center
of the map, $(x^{(0)},y^{(0)})$:

$\displaystyle(x^{(0)},y^{(0)}):=\frac{\sum m_{i}*(x_{i},y_{i})}{\sum
m_{i}}\;\text{.}$ (5)

The samples are weighted by $m_{i}$ as to prevent disk features from affecting
the center of mass.

$\displaystyle m_{i}=\begin{cases}0&\quad\text{when}\quad v_{i}\leq 0.15\\\
\frac{v_{i}-0.15}{0.7}&\quad\text{when}\quad 0.15<v_{i}<0.85\\\
1&\quad\text{when}\quad v_{i}\geq 0.85\end{cases}\text{.}$ (6)

This allows us to refine the equatorial coordinates of the samples as:

$\displaystyle(x_{i}^{(0)},y_{i}^{(0)}):=(x_{i},y_{i})-(x^{(0)},y^{(0)})\;\text{.}$
(7)

From this initial calibration round we iterate, by which we denote the
superscript (c) with $c$ replacing 0 in Eq.s 4 and 7. Each calibration round
produces a new normalised intensity value $v^{(c+1)}$ for each sample. Using
these new values $v^{(c+1)}$, the disk center is further optimised as
$(x^{(c+1)},y^{(c+1)})$. The iterations for producing superscripts $c>0$ are
based on fitting a circle to the boundary points, outliers neglected.

Here we assume a fixed radius for the circle, based on known physical radius
of the Sun as well as known Earth-Sun distance. The physical radius of the
Sun, for $37\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$, is defined by
convergence of iterating the limb profile for the whole MRO dataset in a so
called grand iteration. All maps are first calibrated, after which we fit one
limb profile and temporal latitudinal dependence for each subset of maps
constituting three subsequent quarters of a solar cycle. Then we solve the
relative radius where the average intensity coincides
$0.5\text{\,}\mathrm{Q}\mathrm{S}\mathrm{L}$, and use that radius to fine tune
the calibration for a subsequent grand iteration. In that subsequent grand
iteration, the limb profile at $0.5\text{\,}\mathrm{Q}\mathrm{S}\mathrm{L}$
will be closer to unit radii, and further converges in a few of these grand
iterations. In one grand iteration, any three adjacent quarters of a solar
cycle produce one limb profile and radius correction. For next grand
iteration, we apply them to the middle quarter and for possible boundary
quarters that do not appear in the middle of any triplet of adjacent quarters.
These boundary quarters are the yet incomplete quarter 15th September, 2022 to
15th July, 2025 as well as the first quarter, 3rd March, 1989 to 30th August,
1991, to contain maps in such a format that a limb model is applicable. Thus
every quarter of a solar cycle has its own radius and limbmodel. When
producing a new limb profile and temporal latitudinal dependence as well as a
radius correction, preferably two more quarters, one on both sides, are used
as to prevent ringing artifacts in the middle quarter. The convergent radius
for the most recent three complete quarters (16th June, 2014 to 14th
September, 2022) is $717.6\text{\,}\mathrm{M}\mathrm{m}$, representing the
quarter 16th March, 2017 to 14th December to 16th March, 2019. For other
quarters, the value fluctuates a couple of megameters, and it fluctuates more
for the cycles 22 and 23 which have sparser data. This is a subject of further
study.

(a)(b)

Figure 3: Average quiet Sun intensity as a function of radial distance from
the visual center, for different altitudes (elevations) of observation at MRO.
Fitted from data between 16th March, 2017 and 31st January, 2024. The profile
varies over time, and models fitting each quarter of a Solar cycle have been
derived, which differ only slighly from the displayed example. (a) Limb
brightening is observed when line of sight is almost tangential to the solar
surface. (b) The beam size of $2.4^{\prime}$ mixes the solar limb with the
dark background, resulting in a convolution pattern. On both plots, we observe
more blurring when the altitude of observation is low, due to atmospheric
effects.

###  2.3 Physical radius of the Sun

In this paper, we have assumed that the scale of the relative right ascension
and declination are correct and thus the radius of the disk is based on
geometry. Let $R_{\astrosun}$ be the physical radius of the Sun and $d(t)$ be
the distance between MRO on Earth and the physical center (core) of the Sun at
time $t$. Then:

$\displaystyle
r_{\astrosun}(t)=\arcsin\left(R_{\astrosun}/d(t)\right)\text{.}$ (8)

The apparent radius of the Sun depends on wavelength and method of observation
as well as the exact definition. For radio wavelengths, a quadratic fit of the
apparent radius $\langle r\rangle$ is constructed in Rozelot et al. (2015).
Evaluating this for $37\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$
($\lambda=$8.1\text{\,}\mathrm{m}\mathrm{m}$$ and the avarage distance of
$\langle d\rangle=$1\text{\,}\mathrm{A}\mathrm{U}$$, we obtain
$R_{\astrosun}=$710.06\text{\,}\mathrm{M}\mathrm{m}$$. Selhorst et al. (2019)
examined the variation of the solar radius for
$37\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ from MRO across cycles 22 to 24.
Instrumental changes resulted in differences between upgrades in the measured
radius between 979′′ and 999′′. The values were based on monthly averages of
the radii, which were further averaged over the period of instrumentation
being unchanged. This will normalise for difference of the apparent radius
between perihelion and aphelion, which is ca. $33^{\prime\prime}$. The values
in Selhorst et al. (2019) thus translate to physical diameters
$710\text{\,}\mathrm{M}\mathrm{m}$ and $725\text{\,}\mathrm{M}\mathrm{m}$.

Here, we have analysed the solar radius using two different definitions that
apply to MRO solar maps. We use in this paper the experimental value
$R_{\astrosun}=$717.6\text{\,}\mathrm{M}\mathrm{m}$$, which provides $0.5$ QSL
intensity at the disk boundary when the altitude of the Sun is $25^{\circ}$
above the horizon during observation. This definition is not actually very
sensitive to the exact altitude, as can be seen in Fig. 3(b). However, we fix
the reference altitude in order to establish convergence. The limb profiles
are used for centering the solar maps in the next iteration, from which
subsequent and, eventually final, limb profiles are created. The processes of
correcting for limb brightening and convolution effects are detailed in
Appendix A and accounting for chromospheric effects along different lines of
sight in Appendix B

###  2.4 Heliographic coordinates from radio data

Our solar maps are not sensitive to small systematic pointing errors or
intermap variations in signal levels, since we calibrate each scan
individually. Applying the iterative algorithm outlined in Sect. 2.2 we obtain
a calibrated position and intensity for each sample, with levels scaled such
that zero is for the background and unity for the QSL:

$\displaystyle v_{i}^{(c)}$ $\displaystyle=$
$\displaystyle\frac{u_{i}-{u}_{\mathrm{background}}^{(c)}}{{u}_{\mathrm{QSL}}^{(c)}-{u}_{\mathrm{background}}^{(c)}}\text{.}$
(9)

We use an astronomical model, including the visible solar radius
$r_{\astrosun}$ and the observer’s geolocation (Metsähovi Radio Observatory):

$\displaystyle(\phi^{\earth},\lambda^{\earth})=\left(+60.22^{\circ},+24.39^{\circ}\right)\text{,}$
(10)

to determine the equatorial location of the centre of the solar disk as
observed, denoted by
$(\mathrm{RA}_{\astrosun}(t),\mathrm{dec}_{\astrosun}(t))$. Neglecting the
motion of the Sun across the sky and the rotation of the Sun within the
duration of each scan $S$, for simplicity, we assume the middle time

${t}_{\mathrm{mid}}(S)=\frac{{t}_{\mathrm{min}}(S)+{t}_{\mathrm{max}}(S)}{2}.$

A simple astronomical model is then used to project the cosine-corrected Sun-
relative equatorial coordinates $(x_{i}^{(c)},y_{i}^{(c)})$ into an idealised
heliographic surface set of Carrington coordinates
$(\phi_{\astrosun,i},\lambda_{\astrosun,i})$. This model assumes the Earth in
on a keplerian orbit around the Sun, and the rotation axis of Earth is
precessing at the period of ca. $26\,000\text{\,}\mathrm{a}$. Any inaccuracy
in this simple model will only arise a slight pointing error of the solar
rotation axis on the maps. Since the antenna beam radius is large
($2.4^{\prime}$) this simple astronomical model as adequate.

##  3 Results

(a)(b)(c)

Figure 4: Solar intensity at all heliographic latitudes on
$$\mathrm{3}\mathrm{7}${GHz}$ over cycles 21 to 24. Averaged over samples for
which (a) $r\leq 0.9r_{\astrosun}$ to avoid beam mixing with the background,
and (b) $r\leq 1.0r_{\astrosun}$, such that each sample intensity is corrected
for beam mixing and limb brightening. The latter increases the available
latitude range, revealing the polar magnetic cycle. (c) Count of Metsähovi
solar maps used for analysis after being accepted by automated quality
control.

After disk fitting, each radio sample can be projected onto the Carrington
coordinates of the rotating heliographic surface. Every sample now has time,
intensity relative to QSL and heliographic latitude. We combine these to
construct a time-latitude butterfly diagram. From the series of radio maps
using processing described in Sect. 2, but without applying the limb
correction model (Appendix A), we construct the radio butterfly diagram
displayed in Fig. 4(a).

Fig. 4(a) shows the map radio samples averaged over a grid of time and
heliographic latitude. In the grid, there are $90$ latitude bins, ranging $\pm
90^{\circ}$ of heliographic latitude, and $220$ time bins ranging from March,
1974 to January, 2024. The sample intensities are normalised as described in
2.2 so that for each sample, the intensity is relative to the statistical
Quiet Sun Level of the corresponding map to which the sample belongs.

We present two versions of the butterfly diagram. Fig. 4(a) has no limb
correction, omitting samples that have $r>0.9r_{\astrosun}$ and thus losing
information from the polar regions. Fig. 4(b) applies the limb correction
model identified in Fig. 3, and Appendix A. Comparing the two, it is clear
that without the limb correction model, even for samples at
$0.8r_{\astrosun}<r_{i}<0.9r_{\astrosun}$ the intensities are lower than
physical values due to beam mixing with the background. Data are binned
according to heliographic latitude and time. The value for each bin is an
average of the samples belonging to the bin.

Our time series now extends the previously published data to include the
observations digitally recovered from the historical maps covering cycle 21.
The imbalance in the statistics available across the time span is evident in
more sparse observations during earlier epochs and the seasonal variation in
quantity and quality of the specimens as presented in Fig. 4(c).

##  4 Summary and discussion

Without the limb correction our butterfly diagram Fig. 4(a) during 1993 to
2013 recovers the same activity pattern for latitudes $<|50^{\circ}|$ as
Shibasaki (2013, see their figure 2) at
$17\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$, see also (Kim et al. 2017,
Gopalswamy et al. 2016) and also for
$5.7\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ (Pevtsov 2016). However, only
with the limb correction does our Fig. 4(b) also reproduce the high
intensities approaching polar latitudes. For both observatories the activity
at the Southern pole is more intense and prolonged, and while it appears to
cease in 2012 in the earlier work, from our results we see that this is only a
temporary lull and that this activity has a much longer duration. The results
of Shibasaki (2013) only date from 1993, during which the poles show the onset
of only enhanced activity. Our results including 1989 indicate this was
preceeded by a a period on consistent reduced polar magnetic activity, which
may indicate that the period of the polar magnetic cycle exceeds the 11 year
solar cycle. Extending the observations at
$17\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ to 2018 Gopalswamy et al. (2018)
the activity at the Northern pole reduces after 2013, consistent with our
presentation (see also Fujiki et al. 2019). There is a risk that the signal
near the limb may be corrupted by noise from limb brightening and convolution
within the beam width of limb and background intensity, but the long-term
persistence of weak followed by strong intensities at the poles and
consistency with similar microwave observations supports the conclusion that
these effects are physical rather than instrumental.

Srivastava et al. (2018) explore the properties of the extended solar cycle
(ESC) to around 17 years preceding the 11 year cycle by about 8 years, which
appears to be a residue of earlier cycles evident at higher heliographic
latitude and higher solar atmospheric phenomena. The extended cycle is
identified in the Fe XIV green line emissions in the solar corona by Altrock
(1997) and further analysis using SDO/HMI data (McIntosh et al. 2014) supports
the hypothesis that the magnetic processes leading to the emergence of each
new sunspot cycle have observational signatures apparent a few years in
advance. The extended period of high intensity near the poles apparent in Fig.
4(b) may be relevant to understanding the ESC phenomenon.

North-South asymmetry with a lag as much as two years is identified by Newton
& Milsom (1955) and is evident in asymmetry of our butterfly wings,
particularly for cycles 23 and 24, and in the polar activity after 2010.
Mursula (2023) uses $2.8\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ data to
identify a North-South asymmetry in the overall flux across multiple odd
cycles, as evidence of the presence of a relic field, coexisting with the
dynamo driven magnetic field. As discussed above, there are significant
inbuilt biases in the measurements with respect to North and South poles, so
further work would need to be applied to our data, before we could determine
whether or not our observations support such conclusions.

In summary, we have presented one of the longest continual measurements at
$37\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ of solar chromespheric activity
available, having applied corrections for limb brightening and convolution
effects and atmospheric and orbital biases, providing an accurate historical
catalog of these observations. Our processing of the raw data has permitted us
to include measurements near the solar poles, which indicate a prolonged
period of magnetic activity which may exceed the 11 year solar cycle. The
individual maps contain regions of increased intensity that persist over days
and weeks and indicated enhanced magntic activity. Future work shall present
automatic feature recognition and activity tracking tools and results in order
to relate these to magnetic activity in the chromosphere and elsewhere.

###### Acknowledgements.

This publication makes use of data obtained at the Metsähovi Radio
Observatory, operated by the Aalto University. F.A.G. and M.J.K. acknowledge
support from the Academy of Finland ReSoLVE Centre of Excellence (grant
307411) and the ERC under the EU’s Horizon 2020 research and innovation
programme (Project UniSDyn, grant 818665). S.K.K. acknowledges support from
Niilo Helander Foundation (awarded 2018) as well as Jenny and Antti Wihuri
Foundation (awarded 2019).

## References

  * Alissandrakis et al. (2017) Alissandrakis, C. E., Patsourakos, S., Nindos, A., & Bastian, T. S. 2017, A&A, 605, A78
  * Altrock (1997) Altrock, R. C. 1997, Sol. Phys., 170, 411
  * Babcock (1961) Babcock, H. W. 1961, ApJ, 133, 572
  * Fujiki et al. (2019) Fujiki, K., Shibasaki, K., Yashiro, S., et al. 2019, Sol. Phys., 294, 30
  * Gaizauskas et al. (1983) Gaizauskas, V., Harvey, K. L., Harvey, J. W., & Zwaan, C. 1983, ApJ, 265, 1056
  * Gopalswamy et al. (2018) Gopalswamy, N., Mäkelä, P., Yashiro, S., & Akiyama, S. 2018, Journal of Atmospheric and Solar-Terrestrial Physics, 176, 26
  * Gopalswamy et al. (2016) Gopalswamy, N., Yashiro, S., & Akiyama, S. 2016, ApJ, 823, L15
  * Hale et al. (1919) Hale, G. E., Ellerman, F., Nicholson, S. B., & Joy, A. H. 1919, ApJ, 49, 153
  * Kallunki et al. (2012) Kallunki, J., Lavonen, N., Järvelä, E., & Uunila, M. 2012, Baltic Astronomy, 21, 255
  * Kallunki & Riehokainen (2012) Kallunki, J. & Riehokainen, A. 2012, Astronomische Nachrichten, 333, 20
  * Kallunki et al. (2020) Kallunki, J., Tornikoski, M., & Björklund, I. 2020, Sol. Phys., 295, 105
  * Kallunki et al. (2018) Kallunki, J., Tornikoski, M., Tammi, J., et al. 2018, Astronomische Nachrichten, 339, 204
  * Kim et al. (2017) Kim, S., Park, J.-Y., & Kim, Y.-H. 2017, Journal of Korean Astronomical Society, 50, 125
  * Kivistö (2018) Kivistö, S. 2018, Master’s thesis, Aalto University
  * Krause & Raedler (1980) Krause, F. & Raedler, K. H. 1980, Mean-field magnetohydrodynamics and dynamo theory
  * Leussu et al. (2017) Leussu, R., Usoskin, I. G., Senthamizh Pavai, V., et al. 2017, A&A, 599, A131
  * Loukitcheva et al. (2019) Loukitcheva, M. A., White, S. M., & Solanki, S. K. 2019, ApJ, 877, L26
  * McIntosh et al. (2014) McIntosh, S. W., Wang, X., Leamon, R. J., et al. 2014, ApJ, 792, 12
  * Metsähovi Radio Observatory (2019) Metsähovi Radio Observatory. 2019, Metsähovi Radio Observatory public solar database, http://urn.fi/urn:nbn:fi:att:f371cb6d-f84c-4d76-99e4-c39c639fd0de
  * Mihalas (1978) Mihalas, D. 1978, Stellar atmospheres
  * Mursula (2023) Mursula, K. 2023, A&A, 674, A182
  * Newton & Milsom (1955) Newton, H. W. & Milsom, A. S. 1955, MNRAS, 115, 398
  * Pelt et al. (2006) Pelt, J., Brooke, J. M., Korpi, M. J., & Tuominen, I. 2006, A&A, 460, 875
  * Pevtsov (2016) Pevtsov, A. A. 2016, in Astronomical Society of the Pacific Conference Series, Vol. 504, Coimbra Solar Physics Meeting: Ground-based Solar Observations in the Space Instrumentation Era, ed. I. Dorotovic, C. E. Fischer, & M. Temmer, 71
  * Rozelot et al. (2015) Rozelot, J. P., Kosovichev, A., & Kilcik, A. 2015, ApJ, 812, 91
  * Selhorst et al. (2014) Selhorst, C. L., Costa, J. E. R., Giménez de Castro, C. G., et al. 2014, ApJ, 790, 134
  * Selhorst et al. (2019) Selhorst, C. L., Kallunki, J., Giménez de Castro, C. G., Valio, A., & Costa, J. E. R. 2019, Sol. Phys., 294, 175
  * Shibasaki (2013) Shibasaki, K. 2013, PASJ, 65, S17
  * Srivastava et al. (2018) Srivastava, A. K., McIntosh, S. W., Arge, N., et al. 2018, Frontiers in Astronomy and Space Sciences, 5, 38
  * Steenbeck et al. (1966) Steenbeck, M., Krause, F., & Rädler, K.-H. 1966, Zeitschrift Naturforschung Teil A, 21, 369
  * Valtaoja et al. (1987) Valtaoja, L., Sillanpaa, A., & Valtaoja, E. 1987, A&A, 184, 57
  * Vecchio et al. (2012) Vecchio, A., Laurenza, M., Meduri, D., Carbone, V., & Storini, M. 2012, ApJ, 749, 27

##  Appendix A Correction for limb brightning and antenna convolution

Beam width of the telescope defines the sharpness of the transition from ca.
$1$ QSL to background level (zero) seen in Fig. 3. Diameter of the dish would
suggest a beam diameter $2.4^{\prime}$, but this is further distorted by
atmospheric scattering on Earth. Thus the beam width and shape is also a
function of the altitude angle above horizon as well as local weather.

As discussed in the beginning of Sect. 2.2, signal intensity near the limb
region is affected by beam convolution of the disk and background as well as
increased chromospheric thermal emission along the line of sight. The small
angle of separation between the tangent surface and the chromosphere at the
limb causes the observer to probe higher altitudes of the chromosphere, where
the temperature of the emitting plasma is higher. This produces the limb
brightning, which is seen on Fig. 3 (a) as intensities above $1$ QSL.

We have developed a model to compensate for the effects of limb brightening
and antenna convolution. The model also allows for solar cycle variation and
terrestrial atmospheric scattering. Signals observed from the Sun at very low
altitude travel a longer path through the air than from higher altitudes,
resulting in greater scattering.

The model is iterative, so that we have a set of parameters $M^{(c)}$ for each
step $c$ of the iteration. The parameter set $M^{(c)}$ contains coefficients
of two separate polynomials, $P^{(c)}$ and $Q^{(c)}$:

$\displaystyle M^{(c)}=\left(P^{(c)},Q^{(c)}\right)\text{.}$ (11)

Both polynomials have two-dimensional domain set. The product of polynomials
$P^{(c)}$ and $Q^{(c)}$ is an estimate for the observed intensity based on
four variables:

  * •

$t$, time, phase of solar cycle.

  * •

$l$, heliographic latitude of the intersection of observation line of sight
and the heliographic surface.

  * •

$r$, distance from the center of disk, relative to the apparent radius of the
solar disk.

  * •

$a$, altitude of observation, angle between observation line of sight and the
horizon at MRO.

The model must capture the detail at the solar limb, which is, in practice,
when $r\in[0.8r_{\astrosun},1.2r_{\astrosun}]$. This detail originates from
the convolution by the antenna beam. We map $r$ into adjusted variable $\rho$
using a tuning factor $\kappa$:

$\displaystyle\rho(r)=\arctan\left(\kappa\left(\frac{r}{r_{\astrosun}}-1\right)\right)\;\text{with}\;\kappa=8.8\text{.}$
(12)

The tuning factor $\kappa$ is set by manually searching for a value which
provides the best results.

The estimated intensity, $v$, is then calculated as:

$\displaystyle v=M^{(c)}(t,l,r,a)=P^{(c)}(t,l)\cdot
Q^{(c)}(\rho(r),a)\text{.}$ (13)

The model is optimised using the least squares method in order to minimize the
error between actual observed intensities and the estimates of the model. As a
target set we use the antenna samples of solar cycle 24.

In Fig. 3 (a) radial profiles for limb brightening are plotted for a range of
heliographic latitudes. Fig. 3 (b) for the same range shows the convolution
profiles.

The accuracy of the calibration and centering of the disk is improved by
applying the limb correction to the samples before each iteration round in
Sect. 2.2. It is also of value in determining the radius of the Sun, as
described in Sect. 2.3.

We start with a crude model for the solar disk in order to give each antenna
sample the variables $l$ and $r$ needed in Eq. 13. We then optimize for the
first step of the iteration to get $M^{(1)}$. Thus, the estimated intensity
and the distance from the center of disk are interdependent. When we have the
intensity, we can estimate radial distance as:

$\displaystyle Q^{(c)}(\rho(r),a)=\frac{v}{P^{(c)}(t,l)}\text{.}$ (14)

Given Eq. 14, we can extract $r$ by first substituting $a$ into the polynomial
$Q^{(c)}$ to get one-dimensional polynomial $Q_{a}^{(c)}(\rho)$. We then use
its inverse function to evaluate $\rho$ and $r$:

$\displaystyle r=\left(\rho^{-1}\cdot
Q_{a}^{(c)-1}\right)\left(\frac{v}{P^{(c)}(t,l)}\right)\text{.}$ (15)

The solar map centering is an optimisation problem, which seeks the circle to
best fit the set of limb samples. Using Eq. 15, we can improve the fitting
using the information of how far from the center each limb sample actually
might be. We initially have a circle, centered at $C$, which fits to the limb
samples. We will now translate each sample radially, with respect to $C$, to a
distance set by Eq. 15.

This is illustrated in Fig. 1(c), where the blue dots are true observed limb
samples. The corresponding green dots are translated to a distance given by
Eq. 15. The set of green points provide a less dispersed target set for
fitting a circle. This causes a slight adjustment in the center $C$, so we
will iterate the fitting for a suitable number of rounds.

The accuracy of the calibration and centering of the disk is improved by
applying the limb correction to the samples before each iteration round in
Sect. 2.2. It is also of value in determining the radius of the Sun, as
described in Sect. 2.3.

The limb correction for map centering is illustrated in Fig. 1 (c). While most
observed positions (antenna samples) are denoted with white dots, the subset
at the limb region is shown blue. Each of these blue dots has an arrow
pointing towards a corresponding green dot. Arrows indicate translations
associated with the limb correction. The translation is radial; if the
observed limb sample is more than $0.5\,\textrm{QSL}$, we translate it towards
the center, and if it is less than $0.5\,\textrm{QSL}$, we translate it
outwards. The translations produce the set of green points. The green points
are the target set for fitting a circle, which becomes the disk circumference
and is marked with magenta color. We observe that the green dots appear closer
to the final boundary circle, compared to the blue dots. Since there is less
deviation, we expect a more accurate result for circle fitting when we use the
green dots as a target set.

We note that the interpolation (green) is slightly offset from the
circumference (magenta). The limb correction model compensates variations
associated to altitude of observation. The altitude $a$ relates to line-of-
sight optical thickness $x(a)$ of Earth’s atmosphere as:

$\displaystyle x(a)=\frac{x(90^{\circ})}{\sin a}\text{.}$ (16)

The optical thickness of the air mass along line of sight is lightest when the
Sun is seen high. It is reasonable to assume that the amount of air mass is
significant to the shape of the antenna beam. However, there are other effects
as well, which are not taken into account in this simple model. For example
rain, snow, and turbulence will scatter the beam. The statistical fact of
having winter observations only at low altitudes is taken into account by
normalising the data for altitude.

In coming research, we will test the limb correction on individual maps along
with a more sophisticated observation model. In Fig. 1 (c), we see the antenna
path wiggling along the scanlines. Once the weather effects are compensated,
we will verify the antenna position by analyzing whether the reported wiggling
is consistent with the expected intensity variations due to the beam pointing
slightly towards or away to the disc center at a given sample.

The limb model is based on fitting a polynomial function on the set of all
samples from all accepted solar maps. Since there are various biases in amount
of data, we apply several equalization steps for the weighting of the samples.
This is to avoid ringing artifacts in the limb correction model.

  1. (i)

The number of samples per map varies. Prior to June 2015 samples were more
densely recorded, also taking much longer to scan, as illustrated by Fig. 1
_(a)_. Subsequently the samples are recorded more rapidly with fewer samples
required per map as shown in Fig. 1 _(b)_. The number of samples also varies
annually, since the Sun looks larger when Earth is at perihelion.

  2. (ii)

The heliographic latitude appearing at the centre of each solar map oscillates
annually due to the Earth’s orbital inclination. Thus the number of samples
taken at each heliographic latitude varies between maps, such that
observational altitude is not independent of the heliographic latitude. Within
a map there will be more samples along the visible equator, and on average
along the heliographic equator, where the disc is widest than towards the
poles. During winter, the observation altitude is very low, $\sim 6^{\circ}$,
which adds a bias if favour of one or other hemisphere of the Sun. There is a
risk that observed asymmetry may in fact arise due to this bias.

  3. (iii)

There are around 100 times more maps recorded during summer than during
winter, due to differences in solar altitude and daylight hours.

  4. (iv)

The number of maps, the number of observation days per year, and the number of
observations per day vary greatly over the history of observations (from 1978
to 2024, see also Fig. 4).

To correct for these biases in the data, we apply series of normalisations to
weight the samples equally. This is to prevent ringing artifacts when fitting
the polynomials for the limb model.

  1. 1.

Samples $s_{i}$ for each map are weighted by ca. $r_{i}^{-1}$ as to have all
radial bins an equal contribution.

  2. 2.

The total weight of a map is independent of the amount of samples on the map.

  3. 3.

Each observation day within a year has equal weight.

  4. 4.

Each year from 1978 to 2021 has equal weight. For the incomplete years their
weight is relative to the length of data during that year.

  5. 5.

Each two-dimensional bin of observational sky altitude vs. heliographic
latitude is given equal weight.

In the first two stages, the weights are calculated so as to give each radial
distance an equal representation. Assuming the samples to be uniformly
distributed on the visible solar disk, the number of samples within a radial
bin $r\pm\delta r$ would be $\propto\,\delta r$, yielding a sample weight
$n_{S}^{-1}\left(\operatorname{max}\left\\{r_{i},n_{S}^{-1/2}\right\\}+r_{i}\right)^{-1}$.
Here $n_{S}$ is the amount of samples on the map. The max operator is used to
clip the weight to avoid very large weights for the center samples.

The weighing scheme is not applied to the butterfly diagrams in Fig. 4. We
have, however, tested relevant combinations of different equalisations and the
corresponding butterfly diagrams are available for review. Different
normalisation options do not significantly change the appearance of the
butterfly diagrams.

##  Appendix B Chromosphere effects

In addition to extracting and characterizing features from solar maps, we are
making parallel development with an inversion problem approach. Having
obtained a solar map from MRO, we describe a model for the chromosphere as
well as the MRO telescope which records the map. Each configuration of the
model produces a particular solar map, and we compute this map by ray-tracing
and convolution. We fine-tune the parameters of this model until the map
obtained best fits the actual observation.

As a first approximation, the chromosphere is a spherically symmetric shell
where the physical properties of the plasma depend only on radial distance $r$
from the center of the Sun. When we observe this object using the MRO
$14\text{\,}\mathrm{m}$ dish, we obtain a radially symmetric disk. In the
middle of this disk, we observe an intensity of one Quiet Sun Level. This is
where the line of sight from the observer makes a straight angle with the
chromosphere. When travelling towards the limb, we obtain a profile similar to
what is seen in Fig. 3. Consider angle $\alpha$ between the shell and the line
of sight from the observer. The observed intensity on the map is a function of
this angle, namely, $I(\alpha)$. In the middle of the map,
$\alpha=90^{\circ}$, whereas at the limb, $\alpha\to 0$.

Consider now a point $\bm{w}$ on a line of sight which originates from
observer $\bm{o}$ at Earth and proceeds approximately towards the Sun in unit
direction $\bm{\hat{b}}$:

$\displaystyle\bm{w}(a)=\bm{o}+\left(a+{a}_{\mathrm{offset}}\right)\bm{\hat{b}}\text{.}$
(17)

We denote the travelled length with $a$, starting from offset
${a}_{\mathrm{offset}}$, which is at the effective solar surface of radius
$R_{0}$. This marks a reasonable outer boundary for the chromosphere. We place
the origin at the center of the Sun and assume that the space is transparent
outside radius $R_{0}$ and opaque inside radius $R_{0}-{q}_{\mathrm{max}}$,
thus setting the physical depth of the chromosphere as ${q}_{\mathrm{max}}$.
Function $q(a)$ maps the length travelled along line of sight, measured as
$a$, to the physical depth $q$ as:

$\displaystyle
q(a)=R_{0}-\lVert\bm{o}+\left(a+{a}_{\mathrm{offset}}\right)\bm{\hat{b}}\rVert\text{.}$
(18)

Thus within the chromosphere:

$\displaystyle 0\leq a\leq{a}_{\mathrm{max}}\;\text{and}\;0\leq
q\leq{q}_{\mathrm{max}}\;\text{, with}$ (19) $\displaystyle
q(0)=0\;\text{and}\;q({a}_{\mathrm{max}})\in\\{0,{q}_{\mathrm{max}}\\}\text{.}$
(20)

We emphasize that ${a}_{\mathrm{offset}}$ in Eq. 17 is a function of $\bm{o}$
and $\bm{\hat{b}}$, chosen individually for each line to match the
chromosphere outer boundary. At the central region of the disk,
$q({a}_{\mathrm{max}})={q}_{\mathrm{max}}$ and the line of sight reaches the
bottom shell, while at the limb it is also possible that we see through the
chromosphere and have $q({a}_{\mathrm{max}})=0$. For now, we do not consider
this effect.

We define two functions for describing the chromosphere, $\eta(q)$ and $T(q)$.
These are for opacity and signal temperature, respectively, as measured for
$37\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ radiation as a function of depth,
$q$. Between the two reference spheres, where $0<q<{q}_{\mathrm{max}}$, the
opacity takes finite non-zero values. Otherwise

$\displaystyle\eta(q)$ $\displaystyle=0$ $\displaystyle\;\text{when}\;q\leq
0\;\text{and}$ (21) $\displaystyle\eta(q)$ $\displaystyle=\infty$
$\displaystyle\;\text{when} \;q\geq{q}_{\mathrm{max}}\text{.}$ (22)

The signal temperature $T(q)$ is also measured in
$q\in[0,{q}_{\mathrm{max}}]$, and we denote the bottom temperature as
${T}_{\mathrm{bottom}}=T({q}_{\mathrm{max}})$.

Consider a beam of radiation with intensity $I(a)$ travelling along the line
of sight (Eq. 17) towards the observer at $\bm{o}$. At position measured with
coordinate $a$, the radiation has recently interacted with the plasma along an
infinitesimal path length $\mathrm{d}a$. The radiation at $a$ is thus a sum of
that emitted between positions $a+\mathrm{d}a$ and $a$ and what came at
$a+\mathrm{d}a$ and has not been absorbed until $a$:

$\displaystyle
I(a)=I(a+\mathrm{d}a)*(1-\eta(q(a))\mathrm{d}a)+T(q)\eta(q(a))\mathrm{d}a\text{.}$
(23)

There is an infinitesimal change in intensity,
$\mathrm{d}I=I(a+\mathrm{d}a)-I(a)$, which becomes:

$\displaystyle\mathrm{d}I$ $\displaystyle=$
$\displaystyle\left(I(a+\mathrm{d}a)-T(q)\right)\eta(q(a))\mathrm{d}a$ (24)
$\displaystyle\approx$
$\displaystyle\left(I(a)-T(q)\right)\eta(q(a))\mathrm{d}a\text{.}$ (25)

The opacity $\eta$ thus defines the strength of the coupling between local
plasma and the beam. Since Eq. 24 is linear for $I$, we can obtain the
solution $I(a)$ as sum of solutions for arbitrary temperature profiles, as
long as these profiles add up as $T\left(q(a)\right)$. We use Dirac delta
functions for this and define a particular element $T_{a^{\prime}}(a)$ as:

$\displaystyle
T_{a^{\prime}}(a)=T\left(q(a^{\prime})\right)\delta(a-a^{\prime})\text{.}$
(26)

The opaque background at $a={a}_{\mathrm{max}}$ with temperature
$T({a}_{\mathrm{max}})$ can be treated as an individual contribution to the
observed intensity $I(0)$ as:

$\displaystyle I(0)={I}_{\mathrm{finite}}(0)+{I}_{\mathrm{opaque}}(0)\text{.}$
(27)

Another option is to let the opacity asymptotically reach infinity as:

$\displaystyle\eta({q}_{\mathrm{max}}-\epsilon)\propto\epsilon^{-n}\;\text{with}\;n\geq
1\text{,}$ (28)

so that the spatial integral
$\int\limits_{q=0}^{{q}_{\mathrm{max}}}\\!\eta(q)\,\mathrm{d}q$ becomes
infinite.

The temperature profile $T\left(q(a)\right)$ can be reconstructed as an
integral of contributions from different depth:

$\displaystyle
T\left(q(a)\right)=\int\limits_{a^{\prime}=0}^{{a}_{\mathrm{max}}}\\!T_{a^{\prime}}(a)\,\mathrm{d}a^{\prime}\text{.}$
(29)

We denote the corresponding solution as $I_{a^{\prime}}$, which satisfies:

$\displaystyle\mathrm{d}I_{a^{\prime}}=\left(I_{a^{\prime}}(a)-T_{a^{\prime}}(a)\right)\eta(q(a))\mathrm{d}a\text{.}$
(30)

On the region where $a\neq a^{\prime}$ and thus $T_{a^{\prime}}(a)=0$, we can
integrate Eq. 30 as:

$\displaystyle\int\frac{\mathrm{d}I_{a^{\prime}}}{I_{a^{\prime}}}=\ln\left(I_{a^{\prime}}\right)=\int\eta\left(q(a)\right)\mathrm{d}a\text{.}$
(31)

Thus for small $\epsilon$:

$\displaystyle\ln\left(I_{a^{\prime}}(a^{\prime}-\epsilon)\right)-\ln\left(I_{a^{\prime}}(0)\right)=\int\limits_{a=0}^{a^{\prime}-\epsilon}\\!\eta\left(q(a)\right)\,\mathrm{d}a\text{.}$
(32)

This results in:

$\displaystyle
I_{a^{\prime}}(0)=I_{a^{\prime}}(a^{\prime}-\epsilon)*\exp\left(-\int\limits_{a=0}^{a^{\prime}-\epsilon}\\!\eta\left(q(a)\right)\,\mathrm{d}a\right)\text{.}$
(33)

We set the boundary condition as $I_{a^{\prime}}({a}_{\mathrm{max}})=0$, so
that also:

$\displaystyle I_{a^{\prime}}(a^{\prime}+\epsilon)=0\text{.}$ (34)

To handle the discontinuity of $I_{a^{\prime}}(a)$ at $a=a^{\prime}$, we
integrate over infinitesimal length
$[a^{\prime}-\epsilon,a^{\prime}+\epsilon]$ to get:

$\displaystyle
I_{a^{\prime}}(a^{\prime}-\epsilon)\approx\int\limits_{a=a^{\prime}-\epsilon}^{a^{\prime}+\epsilon}\\!T_{a^{\prime}}(a)\eta(q(a))\,\mathrm{d}a=\eta(q(a^{\prime}))T(a^{\prime})\text{.}$
(35)

We substitute this into Eq. 33 and ignore the infinitesimal gap $\epsilon$:

$\displaystyle
I_{a^{\prime}}(0)=\eta(q(a^{\prime}))T(a^{\prime})*\exp\left(-\int\limits_{a=0}^{a^{\prime}}\\!\eta\left(q(a)\right)\,\mathrm{d}a\right)\text{.}$
(36)

We integrate over all temperature profiles to get the contribution from finite
opacity region $q\in[0,{q}_{\mathrm{max}})$:

$\displaystyle{I}_{\mathrm{finite}}(0)=\int\limits_{a^{\prime}=0}^{{a}_{\mathrm{max}}}\\!I_{a^{\prime}}(0)\,\mathrm{d}a^{\prime}$
(37)
$\displaystyle=\int\limits_{a^{\prime}=0}^{{a}_{\mathrm{max}}}\\!\eta(q(a^{\prime}))T(q(a^{\prime}))*\exp\left(-\int\limits_{a=0}^{a^{\prime}}\\!\eta\left(q(a)\right)\,\mathrm{d}a\right)\,\mathrm{d}a^{\prime}\text{.}$
(38)

For the opaque background contribution in Eq. 27, we have:

$\displaystyle{I}_{\mathrm{opaque}}(0)=T({q}_{\mathrm{max}})*\exp\left(-\int\limits_{a=0}^{{a}_{\mathrm{max}}}\\!\eta\left(q(a)\right)\,\mathrm{d}a\right)\text{.}$
(39)

Instead of using physical distances $a$ and $a^{\prime}$, we’ll switch to
using optical distance $z$ and $z^{\prime}$ as follows:

$\displaystyle\mathrm{d}z=\eta(q(a))\mathrm{d}a\;\text{and}\;\mathrm{d}z^{\prime}=\eta(q(a^{\prime}))\mathrm{d}a^{\prime}\text{.}$
(40)

The scales of $a^{\prime}$ and $z$ will share the same origin, so that:

$\displaystyle
z(a)=\int\limits_{a^{\prime}=0}^{a}\\!\eta(q(a^{\prime}))\,\mathrm{d}a^{\prime}\text{.}$
(41)

We may further approximate that the chromosphere is spatially a thin spherical
shell, so that the observed temperature $I(0)$ is determined only by $\alpha$,
the angle between line of sight and the shell. Thus the dependency between $q$
and $a^{\prime}$ is linear:

$\displaystyle q(a^{\prime})=a^{\prime}*\sin\alpha\text{.}$ (42)

We use symbol $z_{\perp}$ as coordinate for the perpendicular case where
$\alpha=90^{\circ}$. The optical thickness of the chromoshere, seen
perpendically, is ${z}_{\mathrm{max}}$:

$\displaystyle{z}_{\mathrm{max}}$ $\displaystyle=$
$\displaystyle\int\limits_{q=0}^{{q}_{\mathrm{max}}}\\!\eta(q)\,\mathrm{d}q\text{.}$
(43)

We have for general $\alpha$:

$\displaystyle z_{\perp}(a^{\prime})$ $\displaystyle=$ $\displaystyle
z(a^{\prime})*\sin\alpha=\sin\alpha\int\limits_{a=0}^{a^{\prime}}\\!\eta(a)\,\mathrm{d}a\;\text{.}$
(44)

We substitute Eq.s 40 —44 to Eq.s 37 —39:

$\displaystyle{I}_{\mathrm{finite}}(0)=\frac{1}{\sin\alpha}\int\limits_{z_{\perp}=0}^{{z}_{\mathrm{max}}}\\!T(z_{\perp})*\exp\left(-\frac{z_{\perp}}{\sin\alpha}\right)\,\mathrm{d}z_{\perp}\;\text{and}$
(45)
$\displaystyle{I}_{\mathrm{opaque}}(0)=T({z}_{\mathrm{max}})*\exp\left(-\frac{{z}_{\mathrm{max}}}{\sin\alpha}\right)\text{.}$
(46)

When the angle of attack, $\alpha$, is very small, we ultimately observe the
temperature at the very top of chromosphere:

$\displaystyle\alpha\to 0\implies I(0)\to T(0)\text{.}$ (47)

Thus when probing the solar limb, we are measuring the temperature of higher
altitudes of the chromosphere.

Different temperature distributions as $T(z_{\perp})$ produce different limb
brightening patterns when convoluted with the antenna beam We have tried
various polynomial functions, truncated to a variable temperature value at the
opaque background at a variable optical depth ${z_{\perp}}_{\mathrm{max}}$.
The parameters of these functions are optimised using least squares methods,
where we minimize the difference between modellen and observed solar maps.
This requires the antenna beam to be parametrised as well as the physical size
of the Sun. The size of the Sun is optimal when the corresponding temperature
profile is as simple as possible e.g. has no extrema in
$z_{\perp}\in(0,{z_{\perp}}_{\mathrm{max}})$. If we use too large value for
the solar radius, the optimised temperature profile will dip at $z_{\perp}=0$
as it tries to compensate for the background sky, where as too small radius
will result in having a mid-depth temperature minimum when the model
compensates for the oversised rim.

Point sources, such as observed active regions, can be modelled as delta peaks
of intensity, which are smeared out by the antenna convolution. However, since
the emission might be of gyromagnetic origin, it does not radiate
isotropically. This phenomenon has to be treated somehow in future work.