# GeV $\gamma$-ray Emission of Compact Steep-Spectrum Source 4C +39.23B

Ying Gu  Guangxi Key Laboratory for Relativistic Astrophysics, School of
Physical Science and Technology,
Guangxi University, Nanning 530004, People’s Republic of China  Hai-Ming Zhang
School of Astronomy and Space Science, Nanjing University,
Nanjing 210023, People’s Republic of China  Ying-Ying Gan  School of Physics,
Beijing Institute of Technology, Beijing 100081,
People’s Republic of China<EMAIL_ADDRESS>Jin Zhang  School of Physics,
Beijing Institute of Technology, Beijing 100081,
People’s Republic of China<EMAIL_ADDRESS>Xiao-Na Sun  Guangxi Key
Laboratory for Relativistic Astrophysics, School of Physical Science and
Technology,
Guangxi University, Nanning 530004, People’s Republic of China  En-Wei Liang
Guangxi Key Laboratory for Relativistic Astrophysics, School of Physical
Science and Technology,
Guangxi University, Nanning 530004, People’s Republic of China

###### Abstract

Thirteen yr observation data of 4FGL J0824.9+3915 with the Large Area
Telescope on board the _Fermi_ Gamma Ray Space Telescope (_Fermi_ /LAT) are
analyzed for revisiting whether 4C +39.23B, a compact steep-spectrum (CSS)
source closed to a flat-spectrum radio quasar (FSRQ) 4C +39.23A in the
$\gamma$-ray emitting region of 4FGL J0824.9+3915, is a $\gamma$-ray emitter.
We find that the time-integrated $\gamma$-ray emission of 4FGL J0824.9+3915 is
overwhelmingly dominated by 4C +39.23A. It shows significant variability at a
6.7$\sigma$ confidence level and the average $\gamma$-ray flux in the 0.1–300
GeV energy band is $(1.60\pm 0.15)\times 10^{-8}$ ph cm-2 s-1 with a power-law
photon spectral index of $2.48\pm 0.05$. During MJD 57500–58500, 4FGL
J0824.9+3915 is in a low state with a steady $\gamma$-ray flux. Analyzed the
_Fermi_ /LAT observation data in this time interval, it is found that the TS
values of the $\gamma$-ray emission from 4C +39.23A and 4C +39.23B are $\sim
5$ and $\sim 31$, respectively, indicating that the $\gamma$-ray emission in
this time interval is dominated by the CSS 4C +39.23B. The derived average
flux in this time interval for 4C +39.23B is $(9.40\pm 4.10)\times 10^{-9}$ ph
cm -2 s-1 with $\Gamma_{\gamma}=2.45\pm 0.17$. Attributing the spectral energy
distribution (SED) of 4C +39.23B to the radiations from its core and extended
region, we show that the SED can be represented with a two-zone leptonic
model. Its $\gamma$-ray emission is contributed by the core region. The
derived magnetic field strength and Doppler boosting factor of the core are
0.13 G and 6.5. Comparing 4C +39.23B with other $\gamma$-emitting CSSs and
compact symmetric objects (CSOs) in the $\Gamma_{\gamma}-L_{\gamma}$ plane, it
resembles CSSs.

galaxies: active—galaxies: jets—radio continuum: galaxies—gamma rays: galaxies

##  1 Introduction

As a young active galactic nucleus (AGN) population, compact steep-spectrum
sources (CSSs) make up significant fraction ($\sim 30\%$) of centimeter-
wavelength-selected radio sources (O’Dea & Saikia 2021 for a review). They are
characterized by a steep radio spectrum ($\alpha\geq 0.5$,
$F_{\nu}\propto\nu^{-\alpha}$) and a compact dimension ($\leq 15$ kpc, angular
size $\leq 1\arcsec\sim 2\arcsec$) without considering the projection effect
(O’Dea 1998 for a review). More than 70% CSSs may be intrinsically small
rather than a projection effect (Fanti et al. 1990). It was proposed that CSSs
may be the parent population of radio-loud narrow-line Seyfert 1 galaxies (RL-
NLS1s, Caccianiga et al. 2014; Berton et al. 2016; Liao & Gu 2020) because
both CSSs and RL-NLS1s have similar radio properties (Komossa et al. 2006;
Caccianiga et al. 2014; Gu et al. 2015), black hole (BH) masses, and Eddington
ratios (Wu 2009). However, the link between RL-NLS1s and the young radio
sources is also debated, for example, large linear sizes or the lack of
absorption in the X-ray spectrum are observed for some $\gamma$-ray emitting
RL-NLS1s (e.g., Orienti et al. 2015). Superluminal motion was observed in a
few CSSs (Cotton et al. 1997; Taylor et al. 1995; Gawroński & Kus 2006). Their
kinematic ages estimated with the separation velocity of the hotspots and
their radiative ages derived from the synchrotron radiation spectral break are
roughly consistent within the error bars, indicating that CSSs are likely of
young objects. A fraction of CSSs may eventually evolve into typical
Fanaroff–Riley (FR) I or II radio galaxies (RGs; O’Dea 1998; Polatidis &
Conway 2003; Randall et al. 2011).

The $\gamma$-ray emitting CSSs are good candidates for investigating the jet
properties of young AGNs. Only several CSSs so far are detected in the
$\gamma$-ray band with the Large Area Telescope on board the _Fermi_ Gamma Ray
Space Telescope (_Fermi_ /LAT, Abdollahi et al. 2020; Zhang et al. 2020). They
share similar jet radiation properties to the other kinds of $\gamma$-ray
emitting AGNs (Zhang et al. 2020). Recently, Principe et al. (2021) performed
an analysis with the $\sim$11 yr _Fermi_ /LAT observation data to investigate
the $\gamma$-ray emission of 162 young radio sources, and 11 sources among
them have the significant $\gamma$-ray emission.

4C +39.23B is a CSS. It has a double morphology (component-E and component-W)
with a total extent of $\sim$60 mas in the very long baseline array (VLBA)
images at 5.0 and 8.4 GHz (Orienti et al. 2004). The two components have a
steep-spectrum $\alpha\geq 1$ in the range from 1.7 to 5.0 GHz. Component-E is
brighter than component-W, and its core is not detected. The association
between 4C +39.23B and the $\gamma$-ray source 3FGL J0824.9+3916 is reported
in the Third _Fermi_ LAT AGN catalog (3LAC, Ackermann et al. 2015, see also
Massaro et al. 2015). However, it is controversial whether 4C +39.23B is a
$\gamma$-ray emitter since a flat-spectrum radio quasar (FSRQ) 4C +39.23A
($z=1.21568$, Pâris et al. 2018) is also reported to be the counterpart of
3FGL J0824.9+3916 in the 3LAC (Ackermann et al. 2015). Massaro et al. (2015)
also argued that 4C +39.23A may be responsible for the $\gamma$-ray detections
of 3FGL J0824.9+3916. In the updated Fourth _Fermi_ LAT Source Catalog (4FGL,
Abdollahi et al. 2020), 4FGL J0824.9+3915 is suggested to be associated with
4C +39.23A. We dedicate to investigate whether 4C +39.23B is a $\gamma$-ray
emitter by analyzing the $\sim$13 yr _Fermi_ /LAT observational data of 4FGL
J0824.9+3915. The photometric redshift of 4C +39.23B is $z\sim 1.0-1.4$ (Lilly
1989) and $z=1.21$ is used in our analysis (see also Law-Green et al. 1995;
Roche et al. 1998). Throughout, $H_{0}$=71 km s-1 Mpc-1, $\Omega_{m}$=0.27,
and $\Omega_{\Lambda}$=0.73 are used.

##  2 _Fermi_ /LAT Data Analysis

The 13 yr (covering from 2008 August 4 to 2021 October 31) _Fermi_ /LAT
observational data are adopted for our analysis. The publicly available
software Fermitools (ver. 1.2.23) with the binned likelihood analysis method
are used to analyze the data. The region of interest (ROI) is selected as a
$14\arcdeg\times 14\arcdeg$ square region centered at 4FGL J0824.9+3915
(R.A.=126.24$\arcdeg$, Decl.=39.26$\arcdeg$). We select only the $\gamma$-ray
events in the $0.1-300$ GeV energy range with the standard data quality
selection criteria “$(DATA\\_QUAL>0)\&\&(LAT\\_CONFIG==1)$”. Considering the
contamination of the $\gamma$-ray Earth limb, we set the maximum zenith angle
to be 90$\arcdeg$. The instrument response function of $P8R3\\_SOURCE\\_V2$ is
used.

We first examine the possible new sources in the background. We set a model of
known $\gamma$-ray sources in the ROI. The model include the isotropic
emission (“iso_P8R3_V3_v1.txt”), the diffuse Galactic interstellar emission
(IEM, “$gll\\_iem\\_v07.fits$”), and all $\gamma$-ray point sources and
extended sources listed in the _Fermi_ LAT Fourth Source Catalog Data Release
2 (4FGL-DR2, Abdollahi et al. 2020; Ballet et al. 2020). The parameters of the
spectral models of the $\gamma$-ray point sources within ROI and the IEM and
isotropic emission models are kept free. The maximum likelihood test statistic
(TS) is used to estimate the significance of $\gamma$-ray signals, which is
defined by TS$=2(\ln\mathcal{L}_{1}-\ln\mathcal{L}_{0})$, where
$\mathcal{L}_{1}$ is the likelihood value for the background including the
point source and $\mathcal{L}_{0}$ is the likelihood value for the background
without the point source (null hypothesis). Using the package _gtfindsrc_ , we
find a new $\gamma$-ray source located at (R.A.=127.36$\arcdeg$$\pm
0.04\arcdeg$, Decl.=40.33$\arcdeg$$\pm 0.04\arcdeg$) with TS$\sim$45.8, which
is 1.55$\arcdeg$ away from 4FGL J0824.9+3915. The average flux of this new
source is $(2.60\pm 1.00)\times 10^{-9}$ ph cm-2 s-1 with a photon spectral
index of 2.28$\pm 0.15$ in the 0.1–300 GeV energy band.

We add the new point source as a background source to our source model and
perform the analysis again. We use TS${}_{\rm curv}=2\log(\mathcal{L}_{\rm
log-normal}/\mathcal{L}_{\rm power-law})$ to evaluate the curvature
significance of a spectrum. A spectrum is regarded as significantly curved if
its $\rm TS_{curv}>9$ ($\rm TS_{curv}=9$ corresponding to 3$\sigma$, Abdollahi
et al. 2020). We find TS${}_{\rm curv}\sim 4$ for 4FGL J0824.9+3915,
indicating that the spectral shape of 4FGL J0824.9+3915 is well described by a
power-law function. We thus adopt a power-law spectral function to describe
the spectrum of 4FGL J0824.9+3915 (Figure 1), i.e.,
$dN(E)/dE=N_{0}(E/E_{0})^{-\Gamma_{\gamma}}$, where $\Gamma_{\gamma}$ is the
photon spectral index and $N(E)$ is the photon distribution as a function of
energy. The TS map of 4FGL J0824.9+3915 derived with the $\sim$13 yr _Fermi_
/LAT observation data is shown in Figure 2. Using the package _gtfindsrc_ , we
obtain the best-fit position of 4FGL J0824.9+3915, i.e., (R.A.=126.26°$\pm
0.02\arcdeg$, Decl.=39.27°$\pm 0.02\arcdeg$). The position of the $\gamma$-ray
source reported in the 4FGL-DR2 (Abdollahi et al. 2020; Ballet et al. 2020) is
within the 95% error circle of our analysis result. The TS value of the source
is 492.7. The average flux of 4FGL J0824.9+3915 is $(1.60\pm 0.15)\times
10^{-8}$ ph cm-2 s-1 with $\Gamma_{\gamma}=2.48\pm 0.05$ in the 0.1–300 GeV
energy band.

Figure 1 shows the the light curve of 4FGL J0824.9+3915 in a time bin of 180
days with a criterion of TS$\geq$9 (approximately corresponds to 3$\sigma$
detection, Mattox et al. 1996). Significant variability is observed. The
variability index TSvar is commonly used to estimate the significance of
variability for a source, which is defined as

$\rm TS_{\rm
var}=2\sum_{i=1}^{N}[log(\mathcal{L}_{i}(F_{i}))-log(\mathcal{L}_{i}(F_{\rm
glob}))],$ (1)

where $F_{i}$ is the fitting flux for the bin $i$, $\mathcal{L}_{i}(F_{i})$ is
the likelihood corresponding to the bin $i$, and $F_{\rm glob}$ is the best
fit flux for the glob time by assuming a constant flux (Nolan et al. 2012).
$\rm TS_{\rm var}=51.7$ indicates that the variability is at $3\sigma$
confidence level in a $\chi^{2}_{N-1}$(TSvar) distribution with $N-1=25$
degrees of freedom, $N$ is the number of time bins. We obtain TSvar=104.3,
which corresponds to a 6.7$\sigma$ confidence level, indicating that 4FGL
J0824.9+3915 is a variable source111The average flux of 4FGL J0824.9+3915
reported in the 4FGL-DR2 is $F_{1-100\ \rm GeV}=(6.13\pm 0.45)\times 10^{-10}$
ph cm-2 s-1 with $\Gamma_{\gamma}=2.43\pm 0.05$ and the variability index is
TS${}_{\rm var}=35.7$ (the bin size of light curve is one year, Abdollahi et
al. 2020; Ballet et al. 2020) corresponding to a 3.7$\sigma$ confidence level.
With the 13 yr _Fermi_ /LAT data, we obtain $F_{1-100\ \rm GeV}=(5.40\pm
0.50)\times 10^{-10}$ ph cm-2 s-1 with $\Gamma_{\gamma}=2.48\pm 0.05$ and
TS${}_{\rm var}=104.3$ corresponding to a 6.7$\sigma$ confidence level for
4FGL J0824.9+3915. Our results are consistent with that reported in the 4FGL-
DR2..

##  3 Association between 4FGL J0824.9+3915 and 4C +39.23A/B

We show the radio positions of 4C +39.23A (R.A.=126.23$\arcdeg$,
Decl.=39.28$\arcdeg$) and 4C +39.23B (R.A.=126.35$\arcdeg$,
Decl.=39.33$\arcdeg$) together with the 68% and 95% error circles of the best-
fit position for 4FGL J0824.9+3915 in Figure 2. The radio position of 4C
+39.23A falls into the 95% error circle of the best-fit position, indicating
that 4FGL J0824.9+3915 may be spatially associated with 4C +39.23A. 4C +39.23B
is located at 0.12$\arcdeg$ away from the best-fit position and outside of the
95% error circle. Note that the angular resolution of the _Fermi_ /LAT
increases with the energy band from $<3^{\circ}$.5 at 100 MeV to
$<0^{\circ}$.15 at $>$10 GeV (Atwood et al. 2009). Using the gtsrcprob tool,
we find that the maximum energy of the detection photons for 4FGL J0824.9+3915
is only 16.9 GeV. Considering the low spacial resolution of the _Fermi_ /LAT,
one still cannot convincingly exclude the contribution from 4C +39.23B to the
$\gamma$-ray flux of 4FGL J0824.9+3915. It is also possible that there is
another unidentified source nearby the position of 4FGL J0824.9+3915. Thus,
the $\gamma$-ray emission of 4FGL J0824.9+3915 may be contributed from 4C
+39.23A, 4C +39.23B, and/or the unidentified source. We further check this
issue in the following scenarios (see also in Table 1).

  * •

Scenario I: The $\gamma$-ray detections are fully attributed to an
unidentified source at the position of 4FGL J0824.9+3915. In this scenario, we
obtain TS=492.7, $F_{0.1-300\ \rm GeV}=(1.60\pm 0.15)\times 10^{-8}$ ph cm-2
s-1, and $\Gamma_{\gamma}=2.48\pm 0.05$.

  * •

Scenario II: The $\gamma$-ray detections of 4FGL J0824.9+3915 are assumed to
be fully attributed to 4C +39.23A by setting the $\gamma$-ray position at 4C
+39.23A. In this scenario, we obtain TS=491.2, $F_{0.1-300\ \rm GeV}=(1.60\pm
0.16)\times 10^{-8}$ ph cm-2 s-1, and $\Gamma_{\gamma}=2.48\pm 0.05$.

  * •

Scenario III: The $\gamma$-ray detections of 4FGL J0824.9+3915 are assumed to
be fully attributed to 4C +39.23B by setting the $\gamma$-ray position at 4C
+39.23B. In this scenario, we obtain TS=458.6, $F_{0.1-300\ \rm GeV}=(1.60\pm
0.16)\times 10^{-8}$ ph cm-2 s-1, and $\Gamma_{\gamma}=2.50\pm 0.05$.

  * •

Scenario IV: The $\gamma$-ray detections of 4FGL J0824.9+3915 are attributed
to both 4C +39.23A and 4C +39.23B. In this scenario, we get TS=316.1,
$F_{0.1-300\ \rm GeV}=(1.30\pm 0.34)\times 10^{-8}$ ph cm-2 s-1, and
$\Gamma_{\gamma}=2.48\pm 0.07$ for 4C +39.23A, and TS=7.8, $F_{0.1-300\ \rm
GeV}<0.59\times 10^{-8}$ ph cm-2 s-1 (an upper-limit is presented if TS$<$9),
and $\Gamma_{\gamma}=2.52\pm 0.28$ for 4C +39.23B, respectively.

  * •

Scenario V: The $\gamma$-ray detections are attributed to 4C +39.23A, 4C
+39.23B, and an unidentified source at the position of 4FGL J0824.9+3915. In
this scenario, the fits are not convergent.

One can find that the results from the scenarios I–IV are consistent,
indicating that the $\gamma$-rays of 4FGL J0824.9+3915 are dominated by 4C
+39.23A and the hypothetically unidentified source in the position of 4FGL
J0824.9+3915 should be 4C +39.23A. In the scenario V, if the $\gamma$-ray
detections are attributed to 4C +39.23A, 4C +39.23B, and an unidentified
source in the position of 4FGL J0824.9+3915, the fits are not convergent,
indicating that this scenario is unreasonable. As shown in scenarios I–III,
the derived fluxes and photon spectral indices are almost the same, and thus
it is also possible that 4C +39.23B contributes to the $\gamma$-rays. Taking
the contributions from both 4C +39.23A and 4C +39.23B into account, we obtain
TS=316.1 for 4C +39.23A and TS=7.8 for 4C +39.23B, respectively. In this case,
4C +39.23A approximatively contributes to 81% of the $\gamma$-ray flux. One
can also find that the contribution from 4C +39.23B to the $\gamma$-rays is
not negligible.

Note that 4C +39.23A is a FSRQ, a sub-class of blazar, which is characterized
by violent variability. The above analysis shows that 4FGL J0824.9+3915 is a
variable source with a 6.7$\sigma$ confidence level. These results
convincingly indicate that the detected $\gamma$-rays are overwhelmingly
dominated by 4C +39.23A. The weak $\gamma$-ray emission from 4C +39.23B is of
our interest. As shown in Figure 1, 4FGL J0824.9+3915 is in the low state from
MJD 57500 to MJD 58500, which would be offer an opportunity for estimating the
contribution of 4C +39.23B to the $\gamma$-rays of 4FGL J0824.9+3915. We thus
reanalyze the _Fermi_ /LAT observation data in this time interval. The new
source is removed in our analysis. The derived TS map is given in Figure 3.
The gtfindsrc best-fit position is located at (R.A.=126.32$\arcdeg$$\pm
0.07\arcdeg$, Decl.=39.34$\arcdeg$$\pm 0.07\arcdeg$) with TS=52.3 and an
average flux of $(1.10\pm 0.30)\times 10^{-8}$ ph cm-2 s-1. 4C +39.23B is
closer to the best-fit position of the $\gamma$-ray source than 4C +39.23A and
falls into the 68% error circle of the best-fit position. The position of 4C
+39.23A is out of the 68% error circle, but is still within the 95% error
circle. Again, assuming that both 4C +39.23A and 4C +39.23B are responsible
for the $\gamma$-ray detections in this time interval, we obtain TS=31.0
($>5\sigma$) with an average flux of $(9.40\pm 4.10)\times 10^{-9}$ ph cm-2
s-1 for 4C +39.23B while only TS=4.8 is derived for 4C +39.23A. We estimate
the probability of association between the $\gamma$-ray source and 4C +39.23B
in this time interval with a VLBA catalog222To further check our result, we
also estimate the probability of association between the $\gamma$-ray source
and 4C +39.23A using the $\sim$13 yr _Fermi_ /LAT data and obtain 99.1% while
99.85% was reported in the 4FGL-DR2 (Abdollahi et al. 2020; Ballet et al.
2020). , where a Bayesian method with the gtsrcid tool (Abdo et al. 2010) is
used, and the derived association probability is 98.6%. Hence, the
$\gamma$-rays of 4FGL J0824.9+3915 in this time interval are dominated by 4C
+39.23B. The derived average energy spectrum covering MJD 57500–58500 for 4C
+39.23B is shown in Figure 4, which is fitted by a power-law model with
$\Gamma_{\gamma}=2.45\pm 0.17$.

##  4 Jet Properties derived from Broadband SED Modeling

As discussed above, 4C +39.23B would be one more valuable member of the
$\gamma$-ray emitting CSSs. Using the derived average energy spectrum of 4C
+39.23B at the GeV energy band (Figure 4) and other multi-wavelength data from
the NASA/IPAC Extragalactic Database (NED) and ASI Science Data Center
(ASDC)333https://tools.ssdc.asi.it/SED/, we establish the broadband spectral
energy distribution (SED) of 4C +39.23B to investigate its multi-wavelength
radiation properties, as illustrated in Figure 5. We also use the two-zone
leptonic model as reported in Zhang et al. (2020) to reproduce the broadband
SED of 4C +39.23B. The synchrotron (syn), synchrotron-self-Compton (SSC), and
external Compton (EC) scattering of the relativistic electrons in both the
core and extended region are considered. The radio emission below 10 GHz in
the SED should be dominated by the radiation of the extended region on account
of the significant synchrotron-self-absorption effect on the radio emission of
the core.

The electron distribution in both regions is assumed as a broken power-law,
characterized by a density parameter $N_{0}$, a break energy $\gamma_{\rm b}$,
and indices $p_{1}$ and $p_{2}$ in the range of
[$\gamma_{\min},\gamma_{\max}$]. The radiation region is taken as a
homogeneous sphere with radius $R$. The predicted spectrum by the leptonic
model is decided by the following parameters: the size of the radiating region
$R$, the strength of magnetic field $B$, and the Doppler boosting factor
$\delta$ together with the six electron spectrum parameters, where
$\delta=1/(\Gamma-\sqrt{\Gamma^{2}-1}\cos\theta)$, $\Gamma$ and $\theta$ are
the bulk Lorenz factor and viewing angle of the radiation region. The Klein-
Nishina effect and the absorption of high energy $\gamma$-ray photons by
extragalactic background light (Franceschini et al. 2008) are also taken into
account during SED modeling.

For the core region, we take $R=\delta c\Delta t/(1+z)\sim 2.1\times
10^{17}\delta$, where $\Delta t=180$ days, and thus we consider the photons
from the torus provide the seed photons for the inverse Compton process
(IC/torus) of the relativistic electrons in the core region. No viewing angle
is available for 4C +39.23B, we thus assume $\delta$=$\Gamma$ during SED
modeling. For the extended region, the size of the emitting region $R$ is
estimated by the angular radius of 60 mas at the radio band (Orienti et al.
2004). We do not consider the relativistic effect of the extended region and
assume $\Gamma=\delta=1$ during SED modeling. The IC process of the cosmic
microwave background (IC/CMB) by the relativistic electrons in the extended
region is considered. The syn+SSC+IC/CMB model under the equipartition
condition ($U_{B}=U_{\rm e}$) is used to reproduce the radiation of the
extended region, where $U_{B}$ and $U_{\rm e}$ are the energy densities of
magnetic fields and electrons, respectively. The fitting result of the
broadband SED for 4C +39.23B is given in Figure 5 and the derived modeling
parameters are listed in Table 2.

As shown in Figure 5, the SED of 4C +39.23B can be well explained by the model
and the $\gamma$-rays are produced by the IC/tours process. Note that we
cannot give a constraint on the fitting parameters due to the limited
observation data and only obtain a set of model parameters by artificially
adjusting the model parameters to make an acceptable fit. For comparison, the
SED fitting parameters of other $\gamma$-ray emitting CSSs (taken from Zhang
et al. 2020) are also given in Table 2, the derived fitting parameters of 4C
+39.23B are consistent with them.

##  5 Conclusions and Discussion

We have reanalyzed the $\sim$13 yr _Fermi_ /LAT data of 4FGL J0824.9+3915 for
evaluating the contributions from radio sources 4C +39.23A and 4C +39.23B to
the $\gamma$-ray detections. We showed that 4FGL J0824.9+3915 is detected with
a TS value of 492.7. The derived average $\gamma$-ray flux of 4FGL
J0824.9+3915 is $(1.60\pm 0.15)\times 10^{-8}$ ph cm-2 s-1 with a photon
spectral index of $2.48\pm 0.05$ in the 0.1–300 GeV energy band. In the
derived TS map, the radio position of 4C +39.23A falls into the $95\%$ error
circle of the best-fit position, while the radio position of 4C +39.23B is
away from the best-fit position and outside of the $95\%$ error circle.
Considering the contributions of the two sources, we found that the TS values
are 316.1 and $\sim 7.8$ for 4C +39.23A and 4C +39.23B, respectively,
indicating that 4C +39.23A is dominant to the total $\gamma$-rays of 4FGL
J0824.9+3915. The $\gamma$-ray emission of 4FGL J0824.9+3915 is significantly
variable at a 6.7$\sigma$ confidence level, which is overwhelmingly dominated
by the radiations of 4C +39.23A. During MJD 57500–58500, the $\gamma$-ray
emission of 4FGL J0824.9+3915 is in a low state with a steady flux. To
reanalyze the _Fermi_ /LAT observation data in this time interval, it was
found that the TS values of the $\gamma$-ray emission from 4C +39.23A and 4C
+39.23B are 4.8 and 31.0, respectively, indicating that the $\gamma$-ray
emission in this time interval is dominated by the CSS 4C +39.23B. The derived
probability of association between the $\gamma$-ray source and 4C +39.23B in
this time interval with a Bayesian method is 98.6%. The average flux of 4C
+39.23B covering MJD 57500–58500 is $(9.40\pm 4.10)\times 10^{-9}$ ph cm -2
s-1 with a photon spectral index of $2.45\pm 0.17$. Attributing the broadband
SED of 4C +39.23B to the radiations from the core and extended region, we
showed that the SED can be represented by a two-zone leptonic model. Its
$\gamma$-ray emission is contributed by the core region, which has a magnetic
field strength of 0.13 G and a Doppler boosting factor of 6.5.

Our results suggested that 4C +39.23B is a weak $\gamma$-ray emitter. The VLBA
images of 4C +39.23B show a double morphology without detected core (Orienti
et al. 2004) and the total extent of $\sim 60$ mas corresponds to $\sim 0.5$
kpc ($<$1 kpc). 4C +39.23B would be a compact symmetric object (CSO)
considering these radio properties. We compare it with other $\gamma$-ray
emitting young sources in the $\Gamma_{\gamma}-L_{\gamma}$ plane in Figure
6(a), where the five CSOs are CTD 135 (Gan et al. 2021), PMN J1603–4904
(Müller et al. 2014), PKS 1718–649 (Migliori et al. 2016), NGC 3894 (Principe
et al. 2020), and TXS 0128+554 (Lister et al. 2020), and the five CSSs include
3C 138, 3C 216, 3C 286, 3C 309.1, and 3C 380 (Zhang et al. 2020 and references
therein)444We do not take 4C 15.05 into account here since its type is still
debated.. One can observe that 4C +39.23B almost has the highest $\gamma$-ray
luminosity (lower than CTD 135, Gan et al. 2021). In the
$\Gamma_{\gamma}-L_{\gamma}$ plane, the location of 4C +39.23B is consistent
with other five $\gamma$-ray emitting CSSs, and almost overlaps with that of
3C 380. Comparing with the five $\gamma$-ray emitting CSOs, except for CTD
135, 4C +39.23B has the higher $\gamma$-ray luminosity than those CSOs, but
has a steeper spectrum. As reported in Zhang et al. (2020), except for 3C 216,
the other $\gamma$-ray emitting CSSs appear as variable sources and display
obvious variability. For the five $\gamma$-ray emitting CSOs, both CTD 135 and
PMN J1603–4904 display the obvious variability in $\gamma$-rays while the
other three CSOs are not ( Principe et al. 2021; Lister et al. 2020). As
illustrated in Figure 6(b), these young radio sources have the lower ratios of
the $\gamma$-ray luminosity to the radio luminosity ($L_{\rm 5GHz}$) than
blazars. The $\gamma$-rays of blazars are totally dominated by the radiations
of their core jets, and thus the $\gamma$-rays are strongly amplified due to
the Doppler boosting effect while the radio radiations are mainly from the
more extended regions. Comparing with blazars, these young radio sources have
the weakly relativistic effect and thus the lower ratios of $L_{\gamma}$ to
$L_{\rm 5GHz}$.

On the basis of the above analysis results, the $\gamma$-rays of 4FGL
J0824.9+3915 should be attributed to 4C +39.23A and 4C +39.23B. We examined
whether other counterparts make contributions to the $\gamma$-rays of 4FGL
J0824.9+3915. We searched for the low-energy counterparts with the SIMBAD
Astronomical Database and found 19 sources within the 95% error circle of 4FGL
J0824.9+3915 (covering MJD 57500–58500). We further investigated these sources
with the NED. Besides 4C +39.23B and 4C +39.23A, only other two sources (FIRST
J082519.5+391646 and SDSS J082507.90+391859.0) are detected by the Wide-field
Infrared Survey Explorer (WISE). However, the radio flux at 1.4 GHz of FIRST
J082519.5+391646 is much lower than that of 4C +39.23B (more than two orders
of magnitude), and no observational data at higher radio frequencies are
available, while SDSS J082507.90+391859.0 has no radio observations. Among the
four sources, only 4C +39.23A has the X-ray detection. Therefore, the most
likely counterparts of the $\gamma$-ray source 4FGL J0824.9+3915 would be 4C
+39.23A and 4C +39.23B, although we still cannot conclusively rule out the
other sources within the error circle of its position.

This work is supported by the National Natural Science Foundation of China
(grants 12022305, 11973050, U1738136, U1731239, 11851304, and 12133003), and
Guangxi Science Foundation (grants 2017AD22006, 2019AC20334, and
2018GXNSFGA281007).

Figure 1: _Left panel_ : the $\sim$13 yr average energy spectrum along the
power-law function fit for 4FGL J0824.9+3915. _Right panel_ : the long-term
light curve for 4FGL J0824.9+3915 with time bins of 180 days. The inverted
triangle indicates TS$<$9 for this time bin and the horizontal dashed line
represents the $\sim$13 yr average flux of 4FGL J0824.9+3915 observed by the
_Fermi_ /LAT. Figure 2: _Left panel_ : 2.0$\arcdeg\times 2.0\arcdeg$ TS map in
the 0.1–300 GeV energy band observed by the _Fermi_ /LAT covering MJD
54682–59518 for 4FGL J0824.9+3915 (green cross). The black symbols indicate
the best-fit position (black cross) of the $\sim$13 yr _Fermi_ /LAT
observation data together with the 68% and 95% error circles (black circles).
The blue and red crosses show the radio positions of 4C +39.23A (A) and 4C
+39.23B (B), respectively. _Right panel_ : zoom in the central region of the
left panel. The TS map is created with a pixel size of 0.05$\arcdeg$ and
smoothed with a Gaussian kernel of 0.3$\arcdeg$. Figure 3: _Left panel_ :
2.0$\arcdeg$ $\times$ 2.0$\arcdeg$ TS map in the 0.1–300 GeV energy band
observed by the _Fermi_ /LAT covering MJD 57500–58500 for 4FGL J0824.9+3915
(green cross). The black symbols indicate the best-fit position (black cross)
of the _Fermi_ /LAT observation data (covering MJD 57500–58500) together with
the 68% and 95% error circles (black circles). The blue and red crosses show
the radio positions of 4C +39.23A (A) and 4C +39.23B (B), respectively. _Right
panel_ : zoom in the central region of the left panel. The TS map is created
with a pixel size of 0.05$\arcdeg$ and smoothed with a Gaussian kernel of
0.3$\arcdeg$. Figure 4: The average energy spectrum observed by the _Fermi_
/LAT covering MJD 57500–58500 for 4C +39.23B. The black solid line represents
the fitting result with a power-law function and the red bow shows the
1$\sigma$ error. If TS$<$9, an upper-limit is presented for this energy bin.
Figure 5: Observed SED along with the two-zone leptonic model fitting curves.
The data marked as black solid symbols are taken from the NED and ASDC. The
red solid symbols in the $\gamma$-ray band indicate the average energy
spectrum of the _Fermi_ /LAT observations for 4C +39.23B (same as in Figure
4), where the down-triangles represent upper-limits. The thick black solid
line is the sum of emission from each component: synchrotron radiation (red
lines), SSC process (blue lines), and EC process (magenta lines), and among
them the solid lines and dashed lines represent the radiations from the core
and extended region, respectively.

Figure 6: _Panel (a)_ : Comparison of 4C +39.23B with other CSSs and CSOs in
the $\Gamma_{\gamma}-L_{\gamma}$ plane. $\Gamma_{\gamma}$ and $L_{\gamma}$ of
4C +39.23B are derived with the _Fermi_ /LAT observation data covering MJD
57500–58500. $\Gamma_{\gamma}$ and $L_{\gamma}$ of the $\gamma$-ray emitting
CSSs and CSOs are taken from Zhang et al. (2020), Abdollahi et al. (2020), and
Gan et al. (2021), respectively. _Panel (b)_ : Comparison of 4C +39.23B with
other CSSs and COSs as well as blazars in the $L_{\rm 5GHz}-L_{\gamma}$ plane.
The dashed line is the equality line while the solid line is the linear fit in
log scale for all the blazar data, i.e., $\log L_{\rm 5GHz}=(0.47\pm
1.16)+(0.92\pm 0.02)\log L_{\gamma}$. $L_{\rm 5GHz}$ of these $\gamma$-ray
emitting CSSs and CSOs are taken from NED and Lister et al. (2020). $L_{\rm
5GHz}$ of 4C +39.23B are taken from the ASDC. The blazar data are taken from
Fan & Wu et al. (2018).

Table 1: Analysis Results of the $\sim$13 yr _Fermi_ /LAT Observation Data with the Different Templates Templates | $R.A.$ | $Decl.$ | TS | $F_{0.1-300~{}\rm GeV}$ | $\Gamma_{\rm\gamma}$
---|---|---|---|---|---
| [deg] | [deg] |  | [10-8 photon cm-2 s-1] |
4FGL J0824.9+3915 | 126.236 | 39.257 | 492.7 | 1.60$\pm$0.15 | 2.48$\pm$0.05
4C +39.23A | 126.231 | 39.278 | 491.2 | 1.60$\pm$0.16 | 2.48$\pm$0.05
4C +39.23B | 126.349 | 39.329 | 458.6 | 1.60$\pm$0.16 | 2.50$\pm$0.05
4C +39.23A | 126.231 | 39.278 | 316.1 | 1.30$\pm$0.34 | 2.48$\pm$0.07
+ |  |  |  |  |
4C +39.23B | 126.349 | 39.329 | 7.8 |  $<$0.59aaAn upper-limit of flux is given when TS$<$9\.  | 2.52$\pm$0.28
4FGL J0824.9+3915 | 126.236 | 39.257 |  $\cdots$ |  $\cdots$ |  $\cdots$
+ |  |  |  |  |
4C +39.23A | 126.231 | 39.278 |  $\cdots$ |  $\cdots$ |  $\cdots$
+ |  |  |  |  |
4C +39.23B | 126.349 | 39.329 |  $\cdots$ |  $\cdots$ |  $\cdots$
Table 2: SED Fitting Parameters of 4C +39.23B together with that of the other fiveaaWe do not take 4C 15.05 into account here since its type is still debated. CSSs in Zhang et al. (2020) |  | Compact Core | Extended Region |
---|---|---|---|---
Source | z | $R$ | $B$ | $\delta$ | $\Gamma$ | $\theta$ | $\gamma_{\min}$ | $\gamma_{\rm b}$ | $\gamma_{\max}$ | $N_{0}$ | $p_{1}$ | $p_{2}$ |  | $R$ | $B$ | $\gamma_{\min}$ | $\gamma_{\rm b}$ | $\gamma_{\max}$ | $N_{0}$ | $p_{1}$ | $p_{2}$
|  | [cm] | [G] |  |  | [deg] |  |  |  | [cm-3] |  |  |  | [kpc] | [$\mu$G] |  |  |  | [cm-3] |  |
4C +39.23B | 1.21 | 1.37E18 | 0.13 | 6.5 | 6.5 | 8.8 | 1 | 2192 | 1E5 | 0.13 | 1.56 | 4.36 |  | 0.5 | 860 | 500 | 1500 | 7.5E3 | 9.0E-02 | 2.2 | 4
3C 138 | 0.759 | 7.44E17 | 0.6 | 2.8 | 5.5 | 18 | 1 | 531 | 5E5 | 139.8 | 1.8 | 3.46 |  | 2.95 | 277 | 5E2 | 5E3 | 2.5E5 | 6E-03 | 2.2 | 4
3C 216 | 0.67 | 1.01E18 | 0.4 | 3.6 | 3.6 | 15.9 | 1 | 1427 | 1E5 | 25.3 | 1.8 | 4.04 |  | 17.54 | 66 | 1E2 | 2E4 | 1E6 | 8.8E-04 | 2.46 | 4
3C 286 | 0.849 | 1.13E18 | 0.2 | 4.5 | 4.5 | 12.7 | 1 | 1666 | 3E5 | 5.1 | 1.6 | 4.5 |  | 19.18 | 72 | 1E3 | 1.9E4 | 9.5E5 | 4.7E-04 | 2.22 | 4
3C 309.1 | 0.904 | 8.75E17 | 0.22 | 3.6 | 2.8 | 15.5 | 50 | 1914 | 1E5 | 15 | 1.5 | 4.2 |  | 15.65 | 77 | 5E2 | 1.5E4 | 7.5E5 | 4.3E-04 | 2.22 | 4
3C 380 | 0.692 | 1.11E18 | 0.15 | 4.0 | 2.4 | 9.5 | 1 | 1845 | 1E5 | 40.1 | 1.68 | 4.0 |  | 7.11 | 190 | 1E2 | 1.5E4 | 7.5E5 | 5E-03 | 2.4 | 4

## References

  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJS, 188, 405.
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33.
  * Ackermann et al. (2015) Ackermann, M., Ajello, M., Atwood, W. B., et al. 2015, ApJ, 810, 14.
  * Atwood et al. (2009) Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071.
  * Ballet et al. (2020) Ballet, J., Burnett, T. H., Digel, S. W., et al. 2020, arXiv:2005.11208
  * Berton et al. (2016) Berton, M., Caccianiga, A., Foschini, L., et al. 2016, A&A, 591, A98.
  * Caccianiga et al. (2014) Caccianiga, A., Antón, S., Ballo, L., et al. 2014, MNRAS, 441, 172.
  * Cotton et al. (1997) Cotton, W. D., Dallacasa, D., Fanti, C., et al. 1997, A&A, 325, 493
  * Fan & Wu (2018) Fan, X.-L. & Wu, Q. 2018, ApJ, 869, 133.
  * Fanti et al. (1990) Fanti, C., Fanti, R., & O’Dea, C. P. 1990, Compact Steep-Spectrum & GHZ-peaked Spectrum Radio Sources
  * Franceschini et al. (2008) Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837.
  * Gan et al. (2021) Gan, Y.-Y., Zhang, H.-M., Zhang, J., et al. 2021, arXiv:2104.04436
  * Gawroński & Kus (2006) Gawroński, M. P. & Kus, A. J. 2006, Proceedings of the 8th European VLBI Network Symposium, 20
  * Gu et al. (2015) Gu, M., Chen, Y., Komossa, S., et al. 2015, ApJS, 221, 3.
  * Komossa et al. (2006) Komossa, S., Voges, W., Xu, D., et al. 2006, AJ, 132, 531.
  * Law-Green et al. (1995) Law-Green, J. D. B., Eales, S. A., Leahy, J. P., et al. 1995, MNRAS, 277, 995.
  * Liao & Gu (2020) Liao, M. & Gu, M. 2020, MNRAS, 491, 92.
  * Lilly (1989) Lilly, S. J. 1989, ApJ, 340, 77.
  * Lister et al. (2020) Lister, M. L., Homan, D. C., Kovalev, Y. Y., et al. 2020, ApJ, 899, 141.
  * Massaro et al. (2015) Massaro, F., Thompson, D. J., & Ferrara, E. C. 2015, A&A Rev., 24, 2.
  * Mattox et al. (1996) Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396.
  * Migliori et al. (2016) Migliori, G., Siemiginowska, A., Sobolewska, M., et al. 2016, ApJ, 821, L31.
  * Müller et al. (2014) Müller, C., Kadler, M., Ojha, R., et al. 2014, A&A, 562, A4.
  * Nolan et al. (2012) Nolan, P. L., Abdo, A. A., Ackermann, M., et al. 2012, ApJS, 199, 31.
  * O’Dea (1998) O’Dea, C. P. 1998, PASP, 110, 493.
  * O’Dea & Saikia (2021) O’Dea, C. P. & Saikia, D. J. 2021, A&A Rev., 29, 3.
  * Orienti et al. (2004) Orienti, M., Dallacasa, D., Fanti, C., et al. 2004, A&A, 426, 463.
  * Orienti et al. (2015) Orienti, M., D’Ammando, F., Larsson, J., et al. 2015, MNRAS, 453, 4037. doi:10.1093/mnras/stv1845
  * Polatidis & Conway (2003) Polatidis, A. G. & Conway, J. E. 2003, PASA, 20, 69.
  * Principe et al. (2020) Principe, G., Migliori, G., Johnson, T. J., et al. 2020, A&A, 635, A185.
  * Principe et al. (2021) Principe, G., Di Venere, L., Orienti, M., et al. 2021, MNRAS, 507, 4564.
  * Pâris et al. (2018) Pâris, I., Petitjean, P., Aubourg, É., et al. 2018, A&A, 613, A51.
  * Randall et al. (2011) Randall, K. E., Hopkins, A. M., Norris, R. P., et al. 2011, MNRAS, 416, 1135.
  * Roche et al. (1998) Roche, N., Eales, S., & Rawlings, S. 1998, MNRAS, 297, 405.
  * Taylor et al. (1995) Taylor, G. B., Ge, J., & O’Dea, C. P. 1995, AJ, 110, 522.
  * Wu (2009) Wu, Q. 2009, MNRAS, 398, 1905.
  * Zhang et al. (2020) Zhang, J., Zhang, H.-M., Gan, Y.-Y., et al. 2020, ApJ, 899, 2.


11institutetext: Brunnenstr. 58, 61231 Bad Nauheim, Germany 22institutetext:
Max-Planck-Institut für Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077
Göttingen, Germany
22email<EMAIL_ADDRESS>
# No evidence for synchronization of the solar cycle by a “clock”

E. Weisshaar 11 R. H. Cameron 22 M. Schüssler 22

(Received ; accepted)

The length of the solar activity cycle fluctuates considerably. The temporal
evolution of the corresponding cycle phase, that is, the deviation of the
epochs of activity minima or maxima from strict periodicity, provides relevant
information concerning the physical mechanism underlying the cyclic magnetic
activity. An underlying strictly periodic process (akin to a perfect “clock”),
with the observer seeing a superposition of the perfect clock and a small
random phase perturbation, leads to long-term phase stability in the
observations. Such behavior would be expected if cycles were synchronized by
tides caused by orbiting planets or by a hypothetical torsional oscillation in
the solar radiative interior. Alternatively, in the absence of such
synchronization, phase fluctuations accumulate and a random walk of the phase
ensues, which is a typical property of randomly perturbed dynamo models. Based
on the sunspot record and the reconstruction of solar cycles from cosmogenic
14C, we carried out rigorous statistical tests in order to decipher whether
there exists phase synchronization or random walk. Synchronization is rejected
at significance levels of between 95% (28 cycles from sunspot data) and beyond
99% (84 cycles reconstructed from 14C), while the existence of random walk in
the phases is consistent with all data sets. This result strongly supports
randomly perturbed dynamo models with little inter-cycle memory.

###### Key Words.:

Sun: activity – Sun: dynamo

##  1 Introduction

The historical sunspot record reveals significant variability in the amplitude
and length of the individual activity cycles. In order to understand the
physical mechanism underlying the variable cyclic activity, it is important to
answer the question of whether the phase of the activity cycle shows long-term
stability. In this case, although the length of the individual cycles
fluctuates, the cycle appears to remain synchronized to some strictly periodic
process (akin to a “clock”). This means that the fluctuations of the cycle
length do not correspond to real phase changes but reflect, for instance,
observational errors or varying transit times of magnetic flux rising through
the convection zone. Such a situation would occur if the activity cycles were
triggered, for example, by a hypothetical high-Q torsional oscillation in the
solar interior (Dicke, 1970) or by minuscule tidal forces exerted by the
planets (e.g., Stefani et al., 2021). Alternatively, in the absence of such
synchronization in the system, phase fluctuations can lead to a random walk of
the phase. This would happen in the case of a cyclic dynamo with a period that
varies because of random fluctuations of dynamo excitation or meridional flow.
We note that the clock case could also correspond to a dynamo, albeit one with
an externally enforced fixed period.

(Yule, 1927) took an ideal pendulum as a simple mechanical example to
illustrate both cases: if the motion of the pendulum is undisturbed, but the
observation of its position is subject to measurement error, we have a case of
phase stability and synchronization. Alternatively, if the pendulum were to be
randomly bombarded with peas, the resulting phase fluctuations would
accumulate and the phase would show a random walk.

There have been previous attempts in the literature to explore the phase
stability of the solar cycle on the basis of empirical data. The sunspot
record since 1700 provides a relatively consistent dataset, but has generally
been considered to be too short for a reliable statistical phase analysis
(e.g., Hoyng, 1996). Further extension of the cycle into the past has been
attempted on the basis of reports of pre-telescopic sightings of sunspots and
aurorae. The problem there is that the sporadic and often ambiguous reports
are much too sparse to faithfully determine the times of maxima and minima of
individual cycles with yearly precision (e.g., Stephenson, 1988, 1990).
Nevertheless, Schove (1955, 1983) gave years of sunspot minima and maxima back
to the year 649 B.C., mainly relying on aurora reports. However, as he imposed
a fixed number of nine cycles per century, he assumed phase stability
beforehand (Nataf, 2022), meaning that these data are useless for our
purposes. Catalogs of naked-eye observations of sunspots in the pre-telescopic
era (Wittmann, 1978; Wittmann & Xu, 1987) are even more sparse and show no
clustering around the maxima of a hypothetical phase-stable cycle (Stix, 1983,
1984). Furthermore, Carrasco et al. (2020) compared naked-eye observations
with telescopic records in the 19th century and found that the naked-eye
observations do not necessarily correspond to high solar activity or large
sunspot groups, meaning that such observations cannot be used offhandedly to
infer the timing of individual cycle maxima (see also Usoskin et al., 2015).
Similar arguments apply to historical reports considered to represent aurorae,
although long-term activity variations and grand minima may tentatively be
identified in these data. A comprehensive and critical summary of the use of
historical documents to infer solar phenomena in the past was given by Vaquero
& Vázquez (2009).

Although claims of phase stability of the solar cycle can be found in various
papers (e.g., Lomb, 2013; Russell et al., 2019; Stefani et al., 2020), a
conclusive answer to the question of phase stability or otherwise requires
both a clear definition of phase stability and proper statistical analysis.
For instance, nonmonotonic evolution of phase in the sense that sequences of
positive excursions of phase over a number of cycles are seemingly
“compensated for” by later cycles (e.g., Charbonneau & Dikpati, 2000) does not
necessarily entail phase stability, synchronization, or memory (see examples
in Fig. 1). This refers also, for instance, to the well-known weak anti-
correlation between cycle length and amplitude (Hoyng, 1996; Lomb, 2013).

Figure 1: Examples of Monte Carlo simulations of the phase evolution for case
R (random walk, left panels) and for case C (synchronization, right panels)
and different numbers of cycles.

The first systematic statistical studies of the phase stability of the solar
cycle were presented by Dicke (1978) and Gough (1978), who carried out
conceptually similar analyses based on the telescopic sunspot record. Both
authors concluded that the available data were not sufficient to decisively
discriminate between random walk and phase synchronization. Subsequently,
Gough (1981, 1983, 1988) corrected and modified his earlier analysis, still
concluding that the results do not permit the rejection of one of the
alternatives.

In this paper, we revisit the problem using the additional cycle maximum and
minimum epochs until 2019 (28 cycles compared to 24 cycles available in 1981)
as well as the reconstruction of 84 solar cycles between the years 976 and
1888 based on a new analysis of 14C in tree rings (Brehm et al., 2021; Usoskin
et al., 2021). Extending the analyses of Gough (1981, 1983), we find that
these data are statistically consistent with random walk of the cycle phase
and exclude clock synchronization with at high levels of significance.

The paper is organized as follows. In Sect. 2 we consider the disparate phase
evolution of a randomly perturbed, memory-less oscillator (random walk of
phase) and of a perturbed clock-synchronized system and reconsider the
statistical method of Gough (1981, 1983), thereby correcting an error in his
original analysis. We apply the method to the up-to-date sunspot record and to
the 14C-based reconstruction in Sect. 3. Further tests concerning the
significance of the results are presented in Sect. 4. Section 5 contains our
conclusions.

##  2 Random walk of phase versus synchronization: statistical analysis

###  2.1 The method of Gough(1981)

Following Gough (1981), we consider two extreme cases as the basis of our
statistical analysis. The first case corresponds to the absence of inter-cycle
memory: stochastic fluctuations of cycle length lead to random walk of phase.
The alternative case is a perturbed cyclic system whose period is fixed
through synchronization by a perfect clock: any fluctuation in the timing of
cycle number $n$ does not affect the timing of cycle $n+1$. For convenience,
we call the first model random walk (in short: case R) and the second model
clock (in short: case C).

We consider epochs of the cycle minima (or maxima), $t_{n}$ ($n=0,\dots N$).
In both cases (R and C), we have random fluctuations of the cycle length, for
example due to its dependence on fluctuating parameters (case R) or to time
delays of the transfer of magnetic flux through the turbulent convection zone
(case C). In case C, the resulting phase fluctuations do not accumulate,
because the governing periodic process runs with a fixed period and phase,
meaning that previous fluctuations appear to be “corrected” by
synchronization, and we have

$t_{n,\mathrm{C}}=n\cdot P_{\mathrm{C}}+\tau_{n}\,,$ (1)

where $P_{\mathrm{C}}$ is the clock period and $\tau_{n}$ is the (random)
phase fluctuation in cycle $n$. In contrast, for case R, the phase
fluctuations, $\psi_{i}$, accumulate, so that

$t_{n,\mathrm{R}}=n\cdot P_{\mathrm{R}}+\sum_{i=0}^{n}\psi_{i}\,,$ (2)

where $P_{\mathrm{R}}$ is the period at which the undisturbed system would
operate in the absence of fluctations. In practice, both $P_{\mathrm{C}}$ and
$P_{\mathrm{R}}$ correspond to the average of the cycle length over an
infinite set of cycles, so that, in principle, they could have been denoted by
the same symbol. For clarity, we continue to use separate symbols in what
follows. These two simple models confer an advantage in that they are quite
general and independent of the detailed nature of the phase changes (e.g.,
discrete random jumps or more gradual variations during the cycle). Of course,
more explicit models could be analyzed more deeply and therefore possibly
yield more detailed results.

We illustrate the disparate phase evolution for both cases using a simple
example. We assume $P_{\mathrm{R}}=P_{\mathrm{C}}=11\,$yr and normally
distributed phase fluctuations with zero mean and a standard deviation of
$2\,$yr. Figure 1 shows the evolution of the phases for both cases ($\tau_{n}$
and $\sum\psi_{i}$, respectively) and for different total numbers of cycles
(note the different scaling of the ordinate axes). The two cases are virtually
indistinguishable if only ten periods are considered. With 84 cycles, which is
the length of the reconstructed sunspot record by Usoskin et al. (2021),
differences are visible and finally become obvious for 1024 cycles. We note
that random walk of phase does not imply that the phase evolves monotonically;
indeed in our case of a probability distribution for phase fluctuations that
is symmetric with respect to zero, a random walk corresponds to a recurrent
Markov chain that returns to zero infinitely often. Such returns have
therefore nothing to do with any kind of “clock correction of phase
deviations”. The random walk with 1024 cycles clearly shows such returns as
well as intervals with approximately linear phase drift, which, if viewed in
isolation, could easily be mistaken as evidence for phase stability with
respect to a matching period.

In the case of the solar cycle, available datasets are not sufficiently long
to decide the case simply by visual inspection of the phase evolution: a
proper statistical analysis needs to be carried out. Here, we follow the
approach of Gough (1981, 1983), which is sketched below (an example of the
detailed calculations is given in the Appendix).

We have a series of epochs of cycle minima (or maxima), $t_{n}$ $(n=0,\dots
N)$, covering $N$ cycles. The aim of the analysis is to define statistical
criteria that allow us to decide whether these data are consistent with (or
exclude) the clock case (C, synchronization) or random phase walk (case R) as
given in Eqs. (1) and (2). The observed period of cycle $n$,
$P_{n}=t_{n}-t_{n-1}$, is given by

$P_{n,{\mathrm{C}}}=P_{\mathrm{C}}+\tau_{n}-\tau_{n-1}$ (3)

for case C and by

$P_{n,{\mathrm{R}}}=P_{\mathrm{R}}+\psi_{n}$ (4)

for case R. As the “true” periods, $P_{\mathrm{C}}$ and $P_{\mathrm{R}}$, are
not known from the data, they are estimated by the mean period over the $N$
empirical cycles, $\langle P\rangle=(t_{N}-t_{0})/N$, meaning that we have

$\langle P\rangle{{}_{\mathrm{C}}}=P_{\mathrm{C}}+{1\over
N}(\tau_{N}-\tau_{0})$ (5)

for case C and

$\langle P\rangle{{}_{\mathrm{R}}}=P_{\mathrm{R}}+{1\over
N}\sum_{i=1}^{N}\psi_{i}$ (6)

for case R. The phase deviations with respect to the mean period are generally
defined as

$\phi_{n}=t_{n}-t_{0}-n\langle P\rangle\,,$ (7)

meaning that we obtain

$\phi_{n,{\mathrm{C}}}=t_{n,\mathrm{C}}-t_{0,\mathrm{C}}-n\cdot\langle
P\rangle_{\mathrm{C}}=\tau_{n}-\tau_{0}-{n\over N}(\tau_{N}-\tau_{0})$ (8)

for case C, and

$\phi_{n,{\mathrm{R}}}=t_{n,\mathrm{R}}-t_{0,\mathrm{R}}-n\cdot\langle
P\rangle_{\mathrm{R}}=\sum_{i=1}^{n}\psi_{i}-{n\over N}\sum_{i=1}^{N}\psi_{i}$
(9)

for case R. In both cases, we therefore start with zero phase deviation for
$n=0$.

The basis for the statistical analysis are the expectation values of the
variances of cycle period and phase deviation, namely

$\sigma{{}_{\mathrm{P}}}^{2}={1\over N}\sum_{i=1}^{N}(P_{i}-\langle
P\rangle)^{2}$ (10)

for the period and

$\sigma{{}_{\phi}}^{2}={1\over
N+1}\sum_{i=0}^{N}\left(\phi_{i}-\langle\phi\rangle\right)^{2}={1\over
N+1}\sum_{i=0}^{N}\phi_{i}^{2}-\left({1\over
N+1}\sum_{i=0}^{N}\phi_{i}\right)^{2}$ (11)

for the phase deviations, where the angular brackets indicate the arithmetic
mean. For both cases, we assume uncorrelated phase fluctuations ($\tau_{i}$
and $\psi_{i}$, respectively) with zero mean. A somewhat lengthy, but
straightforward calculation (see Appendix) yields the expectation values of
the variances as functions of the number of cycles ($N$) and the amplitudes of
the fluctuations ($\mathcal{T}^{2}$ and $\Psi^{2}$, respectively). For the
expectation value of the variance of period, we obtain

$E(\sigma_{\mathrm{P}}^{2})_{\mathrm{C}}=2\left(1-{1\over
N^{2}}\right)\mathcal{T}^{2}$ (12)

for case C and

$E(\sigma_{\mathrm{P}}^{2})_{\mathrm{R}}={N-1\over N}\Psi^{2}$ (13)

for case R, both in agreement with Gough (1981). The expectation value of the
variance of phase for case C is given by

$\displaystyle E(\sigma_{\phi}^{2})_{\mathrm{C}}$ $\displaystyle=$
$\displaystyle E\left({1\over
N+1}\sum_{i=0}^{N}\phi_{i}^{2}\right)-E\left(\left[{1\over
N+1}\sum_{i=0}^{N}\phi_{i}\right]^{2}\right)$ (14) $\displaystyle=$
$\displaystyle{5N^{2}-6N+1\over 3N(N+1)}\mathcal{T}^{2}-{(N-1)\over
2(N+1)}\mathcal{T}^{2}$ $\displaystyle=$ $\displaystyle{(7N-2)(N-1)\over
6N(N+1)}\mathcal{T}^{2}\,.$

Here, Gough (1981) apparently omitted the second term on the r.h.s. of Eq.
(11), meaning that our result differs from his. Likewise, for case R, we
obtain

$E(\sigma_{\phi}^{2})_{\mathrm{R}}={N-1\over 6}\Psi^{2}-{N(N-1)\over
12(N+1)}\Psi^{2}={(N+2)(N-1)\over 12(N+1)}\Psi^{2}\,.$ (15)

The final statistic is defined as the ratio of the variances, meaning that the
unknown fluctuation amplitudes cancel out. The result is

$S_{\mathrm{C}}={E(\sigma_{\phi}^{2})_{\mathrm{C}}\over
E(\sigma_{\mathrm{P}}^{2})_{\mathrm{C}}}={N(7N-2)\over 12(N+1)^{2}}$ (16)

for case C and

$S_{\mathrm{R}}={E(\sigma_{\phi}^{2})_{\mathrm{R}}\over
E(\sigma_{\mathrm{P}}^{2})_{\mathrm{R}}}={N(N+2)\over 12(N+1)}$ (17)

for case R. For large values of $N$, $S_{\mathrm{C}}$ approaches a constant
value of $7/12,$ while $S_{\mathrm{R}}$ approaches linear growth corresponding
to $N/12$. This is half of the slope obtained by Gough (1981). We cross-
checked the correctness of our analytical calculations using Monte Carlo
simulations.

###  2.2 Modified method (Gough, 1983)

In follow-up papers, Gough (1983, 1988) modified his analysis, referring to
Dicke (1978). Instead of taking the arithmetic mean of the cycle periods
provided by the data as in Eqs. 5 and 6, an estimate of the “true” periods,
$P_{\mathrm{C}}$ or $P_{\mathrm{R}}$, is determined by minimizing the variance
of phase deviation, $\sigma_{\phi}^{2}$. The phases are generally given by
$\phi_{n}=t_{n}-t_{0}-nP$ and $P_{\rm min}$ is determined by the
minimalization process, as

$P_{\rm min}={6\over N(N+1)(N+2)}\sum_{i=0}^{N}(2i-N)t_{i}\,.$ (18)

The resulting periods for the two cases are then given by

$P_{\rm min}^{\mathrm{C}}=P_{\mathrm{C}}+{6\over
N(N+1)(N+2)}\sum_{i=0}^{N}(2i-N)\tau_{i}$ (19)

for case C, and

$P_{\rm min}^{\mathrm{R}}=P_{\mathrm{R}}+{6\over
N(N+1)(N+2)}\sum_{i=0}^{N}(2i-N)\sum_{j=1}^{N}\psi_{j}$ (20)

for case R. Clearly, for growing values of $N$, these expressions converge
more rapidly to the “true” periods than the estimates derived from the average
period given by Eqs. (5) and (6).

Based on the these period estimates, the corresponding ratios of expectation
values for the variances of phase deviation and period are given by

$S_{\mathrm{C}}={E(\sigma{{}_{\phi}}^{2})_{\mathrm{C}}\over
E(\sigma{{}_{\mathrm{P}}}^{2})_{\mathrm{C}}}={N^{2}\over 2(N+1)^{2}},$ (21)

for case C and

$S_{\mathrm{R}}={E(\sigma{{}_{\phi}}^{2})_{\mathrm{R}}\over
E(\sigma{{}_{\mathrm{P}}}^{2})_{\mathrm{R}}}={N(N+3)\over 15(N+1)}$ (22)

for case R, both in agreement with the expressions given by Gough (1983).

Figure 2: Application of the statistical methods of Gough (1981, 1983) to the
epochs of sunspot minima and maxima of 24 activity cycles between 1712 and
1976. Upper panel: (Corrected) Original method of Gough (1981). Lower panel:
Modified method (Gough, 1983). In both panels, the symbols show the ratio
$S=\langle\sigma_{\phi}^{2}\rangle/\langle\sigma_{\mathrm{P}}^{2}\rangle$ as a
function of $N$ (number of cycles in the data segments) for cycle minima
(squares) and maxima (rhombs). The lines indicate the ratio of the expectation
values, $E(\sigma_{\phi}^{2})/E(\sigma_{\mathrm{P}}^{2})$ for case R (random
walk of phase; blue curve) and case C (clock synchronization; red curve),
respectively. The shaded areas show one-sided 68% $(p_{68})$ and 95%
$(p_{95})$ percentiles determined from 10 000 Monte-Carlo simulations for each
case. In both panels, the data points align well with the curve corresponding
to case R (random walk) while case C (synchronization) is rejected at a
confidence level of at least 95%.

Figure 3: Same as Fig. 2, now on the basis of the sunspot data for 28 solar
cycles between 1712 and 2019. Compared to Fig. 2, the confidence intervals
have become somewhat narrower, owing to the longer data set, thus
strengthening the conclusion: clock synchronization is rejected at a
confidence level of at least 95%, while random walk of phase is consistent
with the data.

##  3 Application to observed and reconstructed sunspot cycles

We consider a set of empirical epochs of solar cycle minima (or maxima),
$t_{n}\,(n=0\dots M)$, corresponding to $M$ cycles. The individual periods are
then $P_{i}=t_{i}-t_{i-1}$ $(i=1\dots M)$ and the phase deviations with
respect to the estimate, $P_{\rm est}$, of the “true” period are determined as
$\phi_{n}=t_{n}-t_{0}-n\cdot P_{\rm est}\,(n=0\dots M)$. $P_{\rm est}$ is
calculated either as the arithmetic mean of the individual periods (original
Gough method) or by minimalization of $\sigma_{\phi}^{2}$ (modified method).

The full data set of $M$ cycles is divided into contiguous segments of
$N(q)=M/q$ cycles each, where the values of $q$ are the divisors of $M$. For
each $q$, the mean variances, $\langle\sigma_{\mathrm{P}}^{2}\rangle$ and
$\langle\sigma_{\phi}^{2}\rangle$, are determined according to Eqs. (10) and
(11) as averages over the corresponding segments of length $N(q)$. The ratio
$S=\langle\sigma_{\phi}^{2}\rangle/\langle\sigma_{\mathrm{P}}^{2}\rangle$ is
then plotted as a function of $N(q)$, to be compared with the curves
corresponding to the ratio of the analytically calculated expectation values,
$S_{\mathrm{C}}$ and $S_{\mathrm{R}}$.

Confidence intervals with regard to $S_{\mathrm{C}}$ and $S_{\mathrm{R}}$ are
determined by Monte Carlo simulations. We assume normally distributed phase
fluctuations $\tau_{n}$ and $\psi_{n}$, respectively, with zero mean and a
standard deviation of 2 yr to compute 10 000 realizations of $M$ perturbed 11
yr cycles for each of the two cases according to Eqs. (1) and (2). The assumed
standard deviation is consistent with the value inferred for 84 solar cycles
reconstructed from 14C by Usoskin et al. (2021), but the precise choice is
irrelevant, as the statistic $S$ neither depends on the amplitude nor on the
standard deviation of the distribution of the phase fluctuations. Following
the procedures described in the previous section, we then calculate the values
of $S$ for each value of $N(q)$ and for all realizations. From the resulting
distributions of $S_{\mathrm{C}}$ and $S_{\mathrm{R}}$, we determine the one-
sided 68%, 95%, and 99% percentiles: as the distributions are generally not
symmetric, we start from the average value and separately count the number of
realizations until we reach the given percentile on each side. For a normal
distribution, these percentiles would (approximately) correspond to the $\pm
1\sigma$, $\pm 2\sigma$, and $\pm 3\sigma$ ranges, where $\sigma$ denotes the
standard deviation.

As a first application, we consider the telescopic sunspot record. Gough
(1981, 1983) considered the maximum and minimum epochs of 24 cycles between
the years 1712 and 1976. Using the epochs provided by the SILSO database
(www.sidc.be/silso), we redid his analysis using both the (corrected) original
method and the modified method as presented in Sect. 2. The results are shown
in Fig. 2. Already with the 24 cycles available to Gough (1981), the data
points for both variants of the method align well with the blue line
representing the expectation values for case R. Clock synchronization (case C)
can be rejected at the 95% confidence level. Owing to its better convergence
property, the modified method provides somewhat narrower confidence intervals
than the original method. The differences between the diagram in the left
panel of Fig. 2 and Fig. 7 of Gough (1981) result from the correction of the
calculation of the variance of phase (see Sect. 2). This leads to a downward
shift of the theoretical curves for both cases. The result is that even the
limited data available in the early 1980s clearly favor case R,that is, a
random walk of the phase of the solar cycle. This is in accordance with the
results presented in Gough (1983) and Gough (1988).

In the meantime, data for 28 cycles have become available. Figure 3 shows the
results of the statistical analysis. The longer data series confirms and
strengthens the results obtained on the basis of 24 cycles. In particular, the
confidence intervals have become narrower compared to Fig. 2, meaning that the
rejection of the synchronization hypothesis becomes even more significant.

The recent reconstruction of 84 cycles (between the years 976 and 1888) on the
basis of the cosmogenic isotope 14C by Usoskin et al. (2021) provides a
significant extension of the database available for the analysis. Figure 4
shows the results of applying both variants of the Gough method to these data;
they confirm and significantly reinforce the results drawn from the sunspot
data: clock synchronization is now rejected at a confidence level of at least
99%, while random walk remains consistent with the data.

Figure 4: Application of the statistical analysis to the series of 84 cycles
covering the period between the years 976 and 1888 as reconstructed from 14C
data by (Usoskin et al., 2021). Upper panel: (Corrected) Original method of
Gough (1981). Lower panel: Modified method (Gough, 1983). The graphical
representation of the results is the same as in the previous figures. The data
are consistent with random walk of phase (case R) while, owing to the much
longer data set, clock synchronization is rejected at levels exceeding 99%.

The longer set of reconstructed cycles allows us to also consider the
distribution of the phase deviations, $\phi_{n}$, as a means to qualitatively
discriminate between the cases C and R. The upper row of Fig. 5 shows the
distribution of phase deviations for the 84 reconstructed cycles from the 14C
record (relative to the mean period of 10.9 yr) in comparison to the
corresponding distributions from one realization of the Monte Carlo
simulations of 84 artificial cycles for case R and the same for case C. Both
simulations assume a base period of 11 years with a standard deviation of 2.2
years for the phase fluctuations. As the phase deviations accumulate in case R
(upper middle panel), the corresponding distribution is significantly broader
($\sigma=6.3\,$yr) than that of case C (upper right panel, $\sigma=2.2\,$yr),
the latter reflecting the input Gaussian distribution. Likewise, the
distribution for the 84 reconstructed cycles (upper left panel,
$\sigma=6.2\,$yr) is also markedly broader than that of the simulated case C.
For case R, we can retrieve the individual phase fluctuation, $\psi_{n}$, of
cycle $n$ as the difference $\phi_{n}-\phi_{n-1}$. This is shown in the lower
row of Fig. 5. For the simulated case R, the result (based on the mean period)
therefore closely reflects the input distribution (lower middle panel).
However, if we apply the same procedure to the simulated case C, this
necessarily leads to a broader distribution than the input, because it results
from the difference of two Gaussian distributions, $\tau_{n}-\tau_{n-1}$
(lower right panel). We can therefore use this procedure as a qualitative
method to distinguish whether a given dataset corresponds to case R or to case
C: if the distribution of $\phi_{n}-\phi_{n-1}$ is narrower than the original
distribution, this indicates case R, while a broader distribution favors case
C. The two panels corresponding to the 14C reconstructions on the left side of
Fig. 5 clearly show that the reconstructed solar cycles behave as expected for
case R, that is, the distribution of the inferred phase fluctuations for the
individual cycles is significantly narrower ($\sigma=2.2\,$yr) than that of
the original phase deviations ($\sigma=6.2\,$yr).

Figure 5: Distributions of phase deviations (upper panels) and inferred phase
fluctuations (differences between subsequent phase deviations, lower panels)
for the 84 reconstructed cycles from the 14C data (left panels) and for
simulated cycles with Gaussian fluctuations corresponding to case R (middle
panels) and case C (right panels). The dotted curves indicate Gaussian fits.

##  4 Further tests

###  4.1 Systematic effects

We consider two systematic effects that can influence the epochs of activity
maxima and minima. Firstly, the “Waldmeier effect”, which is the correlation
between the rise time from activity minimum to maximum and the amplitude of a
cycle (Waldmeier, 1935), affects the times of the maxima. This effect enhances
the variance of the period determined from the maxima and possibly explains
why the empirical $S$ values corresponding to the set of cycle maxima tend to
be somewhat lower than those for the minima (see Figs. 2–4). Following Dicke
(1978), Gough (1981, 1983) corrected the sunspot data for the Waldmeier effect
and found that the results were not significantly affected.

Secondly, activity cycles overlap, that is, flux emergence of a new cycle
already starts at mid solar latitudes while the previous cycle is still
present at low latitudes. In combination with the Waldmeier effect, this leads
to an amplitude-dependent shift of the minimum epochs (Cameron & Schüssler,
2007). Incidentally, this shift also explains the empirical statistical
relationship between cycle length and amplitude. Hathaway et al. (1994) and
Hathaway (2011) used curve fitting to an empirical standard cycle shape to
obtain starting dates of sunspot cycles 1 to 23 (1755 until 1996) that are not
affected by cycle overlap. The result of applying the modified Gough test
(Sect. 2.2) to these data is given in Fig. 6. Comparison with the square
symbols presented in Figs. 2b and 3b shows that the shift of the minimum
epochs due to cycle overlap and the Waldmeier effect does not significantly
affect the statistical results of the Gough test.

Figure 6: Application of the modified method of Gough (1983) to the unbiased
starting dates of solar cycles 1 to 23 determined by Hathaway et al. (1994)
and Hathaway (2011).

###  4.2 Dependence on the length of the data set

The length of the data set, that is, the number of observed cycles, is crucial
for the capability of a statistical test to distinguish between case C and
case R. Our results show that the set of 84 cycles reconstructed from the 14C
record permits a much more stringent answer than the 28 cycles provided by the
sunspot record. In order to systematically assess the reliability of the
original and modified Gough methods, we used Monte Carlo simulations to
determine the false-alarm probability (FAP) with respect to the 99% confidence
interval as a function of the total number of cycles, $M$. For each given $M$,
we computed $10^{6}$ realizations of each of the models, that is, case C (Eq.
3) and case R (Eq. 4). For each realization, we calculated the value of our
statistic, $S$, for $N=M$, which provides the largest difference between
$S_{\mathrm{R}}$ and $S_{\mathrm{C}}$, and is therefore the most sensitive. We
then determined the FAPs, that is, the percentage of the realizations of case
C with $S$ values within the 99% confidence interval for case R and,
conversely, the percentage of the realizations of case R with $S$ values
within the 99% confidence interval for case C.

Figure 7: False-alarm probabilities as a function of the total number of
cycles in the data set, based on Monte Carlo simulations. Upper panel:
(Corrected) Original method of Gough (1981). Lower panel: Modified method
(Gough, 1983). The red lines give the FAP that a realization of case C falls
within the 99% confidence interval for case R, and vice versa for the blue
lines. The dashed line segment indicates an upper limit.

The results of this procedure for both the original and the modified Gough
test are shown in Fig. 7. The red curves correspond to the FAP for case R,
meaning the probability that a realization of case R would fall within the 99%
confidence interval of case C, or in other words the probability of
misclassification as case C at that level. Likewise, the blue lines give the
probability for misclassfying a realization of case C as case R. The figures
clearly show that the modified method provides significantly lower values of
the FAP, particularly for larger $M$. Reasonably low values of less than 1%
FAP are reached for $M>50$ for the original Gough test and $M>30$ for the
modified test. This corresponds to the results shown in Sect. 3: while the
clock hypothesis can be safely rejected at the 99% significance level for
$M=84$ using the 14C-based data, this can only be done at 95% significance
level for $M=24(28)$ on the basis of the sunspot record.

###  4.3 The method of Eddington and Plakidis

A problem very similar to the question of phase stability of the solar cycle
was considered by Eddington & Plakidis (1929). These authors considered the
question of whether the brightness changes of long-period variable stars have
an underlying fixed period with superposed observational errors (analogous to
the clock case in the context of the present study) or if the individual
periods vary randomly without memory (i.e., the random walk case). The authors
considered the mean value, $\bar{u}_{x}$, of the accumulated phase error after
$x$ epochs of maximum brightness with respect to the mean period over the
whole length of the data set. Assuming that the period variations consist of a
combination of accumulating (random walk) and nonaccumulating (clock-type)
fluctuations, the expectation value is given by

$\bar{u}_{x}=\sqrt{2\alpha^{2}+\epsilon^{2}x\left(1-{x\over M}\right)},$ (23)

where $\alpha$ is the amplitude of the synchronized fluctuations (e.g.,
measurement errors) and $\epsilon$ the amplitude of the cumulating
fluctuations (leading to random walk of phase). A correction factor (in
parentheses) is applied to account for the cases where $x$ becomes a
considerable fraction of the total number of periods, $M$, in the data set.
$\bar{u}_{x}$ can be calculated from the data and a least-squares fit to Eq.
23 provides empirical values of the amplitudes $\alpha$ and $\epsilon$, from
which the relative importance of clock and random walk can be assessed.

We calculated $\bar{u}_{x}$ on the basis of the minimum epochs of the 84
activity cycles determined from the 14C data. The results and the fit to Eq.
23 are shown in Fig. 8. As the individual values of the accumulated phase,
$u_{x}$, are significantly affected by overlap for large values of $x$, we
restricted the fit to $x\leq 50$. The fit results in $\alpha<10^{-4}\,$yr for
the amplitude of the nonaccumulating fluctuations and $\epsilon=1.9\,$yr for
the amplitude of the cumulating fluctuations. This clearly demonstrates that
the system is nearly exclusively governed by random walk of the cycle phase,
thus corroborating the results of the Gough test.

Figure 8: Application of the test of Eddington & Plakidis (1929) to the 84
reconstructed cycles on the basis of the 14C data. Shown is the mean
accumulated phase error after $x$ cycles as a function of $x$. The red curve
results from the activity minima of the reconstructed cycles, the blue curve
is the corresponding fit to the theoretical relation given by Eq. 23 for
$x\leq 50$. The resulting fit parameters clearly indicate that the system
shows random walk of phase.

##  5 Conclusions

Using rigorous statistical tests, we investigated whether or not the phase of
the solar cycle is stable, that is, whether it is synchronized by some kind of
clock or follows a random walk. Our analysis is based on 28 activity cycles
from the historical sunspot record and 84 cycles reconstructed from 14C data
in tree rings. The tests reveal that these data unequivocally favor random
walk of the cycle phase and exclude clock synchronization, with high
significance levels (95% for sunspot data and greatly exceeding 99% for the
14C data). Consequently, data that unambiguously correspond to solar activity
cycles provide no basis for the “planetary hypothesis” or for hypothetical
high-Q oscillations in the radiative interior as drivers or “synchronizers” of
the solar cycle. Likewise, data for an exoplanetary system also provide no
evidence for planetary synchronization of stellar cycles (Obridko et al.,
2022). The consistency of the phase evolution with random walk suggests that a
randomly perturbed, memoryless dynamo action is the cause of the solar
activity cycle (e.g., Charbonneau & Dikpati, 2000; Cameron & Schüssler, 2017;
Karak & Miesch, 2018).

###### Acknowledgements.

Ilya Usoskin kindly provided the reconstructed sunspot numbers obtained by
Usoskin et al. (2021) in electronic form. We are grateful to the referee for
making us aware of the paper by Eddington & Plakidis (1929).

## References

  * Brehm et al. (2021) Brehm, N., Bayliss, A., Christl, M., et al. 2021, Nature Geoscience, 14, 10
  * Cameron & Schüssler (2007) Cameron, R. H. & Schüssler, M. 2007, ApJ, 659, 801
  * Cameron & Schüssler (2017) Cameron, R. H. & Schüssler, M. 2017, ApJ, 843, 111
  * Carrasco et al. (2020) Carrasco, V. M. S., Gallego, M. C., Arlt, R., & Vaquero, J. M. 2020, ApJ, 904, 60
  * Charbonneau & Dikpati (2000) Charbonneau, P. & Dikpati, M. 2000, ApJ, 543, 1027
  * Dicke (1970) Dicke, R. H. 1970, in IAU Colloq. 4: Stellar Rotation, ed. A. Slettebak (Reidel, Dordrecht), 289
  * Dicke (1978) Dicke, R. H. 1978, Nature, 276, 676
  * Eddington & Plakidis (1929) Eddington, A. S. & Plakidis, S. 1929, MNRAS, 90, 65
  * Gough (1978) Gough, D. 1978, in Pleins Feux sur la Physique Solaire, ed. S. Dumont & J. Roesch, 81
  * Gough (1981) Gough, D. 1981, in NASA Conference Publication, Vol. 2191, 185–206
  * Gough (1983) Gough, D. O. 1983, ESA Journal, 7, 325
  * Gough (1988) Gough, D. O. 1988, in Solar-Terrestrial Relationships and the Earth Environment in the last Millennia, ed. G. Cini Castagnoli (North-Holland, Amsterdam etc.), 90
  * Hathaway (2011) Hathaway, D. H. 2011, Sol. Phys., 273, 221
  * Hathaway et al. (1994) Hathaway, D. H., Wilson, R. M., & Reichmann, E. J. 1994, Sol. Phys., 151, 177
  * Hoyng (1996) Hoyng, P. 1996, Sol. Phys., 169, 253
  * Karak & Miesch (2018) Karak, B. B. & Miesch, M. 2018, ApJ, 860, L26
  * Lomb (2013) Lomb, N. 2013, in Journal of Physics Conference Series, Vol. 440, 012042
  * Nataf (2022) Nataf, H.-C. 2022, Sol. Phys., 297, 107
  * Obridko et al. (2022) Obridko, V. N., Katsova, M. M., & Sokoloff, D. D. 2022, arXiv e-prints, arXiv:2208.06190
  * Russell et al. (2019) Russell, C. T., Jian, L. K., & Luhmann, J. G. 2019, Rev. of Geophys., 57, 1129
  * Schove (1955) Schove, D. J. 1955, J. Geophys. Res., 60, 127
  * Schove (1983) Schove, D. J. 1983, Sunspot cycles (Hutchinson Ross Publ. Comp., Benchmark Papers in Geology, Vol.68))
  * Stefani et al. (2020) Stefani, F., Beer, J., Giesecke, A., et al. 2020, Astronomische Nachrichten, 341, 600
  * Stefani et al. (2021) Stefani, F., Stepanov, R., & Weier, T. 2021, Sol. Phys., 296, 88
  * Stephenson (1988) Stephenson, F. R. 1988, in Solar-Terrestrial Relationships and the Earth Environment in the last Millennia, ed. G. Cini Castagnoli (North-Holland, Amsterdam etc.), 133
  * Stephenson (1990) Stephenson, F. R. 1990, Phil. Trans. Roy. Soc. London, Series A, 330, 499
  * Stix (1983) Stix, M. 1983, Mitteilungen der Astronomischen Gesellschaft Hamburg, 60, 95, https://ui.adsabs.harvard.edu/abs/1983MitAG..60…95S
  * Stix (1984) Stix, M. 1984, Astronomische Nachrichten, 305, 215
  * Usoskin et al. (2015) Usoskin, I. G., Arlt, R., Asvestari, E., et al. 2015, A&A, 581, A95
  * Usoskin et al. (2021) Usoskin, I. G., Solanki, S. K., Krivova, N. A., et al. 2021, A&A, 649, A141
  * Vaquero & Vázquez (2009) Vaquero, J. M. & Vázquez, M. 2009, Astrophysics and Space Science Library, Vol. 361, The Sun Recorded Through History (Springer)
  * Waldmeier (1935) Waldmeier, M. 1935, Astronomische Mitteilungen der Eidgenössischen Sternwarte Zürich, 14, 105
  * Wittmann (1978) Wittmann, A. 1978, A&A, 66, 93
  * Wittmann & Xu (1987) Wittmann, A. D. & Xu, Z. T. 1987, A&AS, 70, 83
  * Yule (1927) Yule, G. U. 1927, Phil. Trans. Roy. Soc. London A, 226, 267

##  Appendix A Calculation of the statistic $S_{R}$

In what follows, we describe the calculation of the statistic
$S=E\\{\sigma_{\Phi}^{2}\\}/E\\{\sigma_{\mathrm{P}}^{2}\\}$ for case R (random
walk of phase) as would be expected for a randomly perturbed dynamo. We follow
Gough (1981) and write, for the $n^{\rm th}$ observed period,

$P_{n}=P_{\mathrm{R}}+\psi_{n}\,,$ (A.1)

where $P_{\mathrm{R}}$ is the “true” period (in the absence of fluctuations)
and $\psi_{n}$ denotes random, uncorrelated perturbations with zero mean and
amplitude $\Psi$, namely

$E\\{\psi_{n}\\}=0\;\;\;\mathrm{and}\;\;\;E\\{\psi_{n}\psi_{n^{\prime}}\\}=\Psi^{2}\delta_{nn^{\prime}}\,.$
(A.2)

The unknown period $P_{\mathrm{R}}$ is estimated via the arithmetic mean of
the periods given by

$\langle P\rangle={1\over N}\sum_{i=1}^{N}P_{i}=P_{\mathrm{R}}+{1\over
N}\sum_{i=1}^{N}\psi_{i}\,.$ (A.3)

For the expectation value of the variance of period, we obtain

$\displaystyle E\left\\{\sigma{{}_{\mathrm{P}}}^{2}\right\\}$ $\displaystyle=$
$\displaystyle E\left\\{{1\over N}\sum_{i=1}^{N}\left(P_{i}-\langle
P\rangle\right)^{2}\right\\}=E\left\\{{1\over
N}\sum_{i=1}^{N}\left(\psi_{i}-{1\over
N}\sum_{j=1}^{N}\psi_{j}\right)^{2}\right\\}=$ (A.4) $\displaystyle=$
$\displaystyle E\left\\{{1\over
N}\sum_{i=1}^{N}\left[\psi_{i}^{2}-{2\psi_{i}\over
N}\sum_{j=1}^{N}\psi_{j}+{1\over
N^{2}}\sum_{k=1}^{N}\sum_{\ell=1}^{N}\psi_{k}\psi_{\ell}\right]\right\\}=\Psi^{2}-{2\Psi^{2}\over
N}+{\Psi^{2}\over N}={N-1\over N}\Psi^{2}\;.$

The phase deviations are defined by $\phi_{n}=t_{n}-t_{0}-n\cdot\langle
P\rangle$, where

$t_{n}=n\cdot P_{\mathrm{R}}+\sum_{i=0}^{n}\psi_{i},$ (A.5)

and so we obtain

$\phi_{n}=t_{n}-T_{n}=\sum_{i=1}^{n}\psi_{i}-{n\over
N}\sum_{i=1}^{N}\psi_{i}\,.$ (A.6)

We now consider the expectation value of the variance of the phase deviation,

$E\left\\{\sigma_{\phi}^{2}\right\\}=E\left\\{{1\over
N+1}\sum_{i=0}^{N}\phi_{i}^{2}-\left({1\over
N+1}\sum_{i=0}^{N}\phi_{i}\right)^{2}\right\\}\equiv{\rm I}-{\rm II}\,,$ (A.7)

which consists of two terms (denoted by I and II). For calculating term I, we
first consider

$\displaystyle E\left\\{\phi_{n}^{2}\right\\}$ $\displaystyle=$ $\displaystyle
E\left\\{\sum_{i=1}^{n}\sum_{j=1}^{n}\psi_{i}\psi_{j}-{2n\over
N}\sum_{i=1}^{n}\psi_{i}\sum_{j=1}^{N}\psi_{j}+{n^{2}\over
N^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N}\psi_{i}\psi_{j}\right\\}=$ (A.8)
$\displaystyle=$ $\displaystyle n\Psi^{2}-{2n^{2}\over N}\Psi^{2}+{n^{2}\over
N}\Psi^{2}=n\left(1-{n\over N}\right)\Psi^{2}\;,$

which leads to

${\rm I}={1\over
N+1}\sum_{i=0}^{N}E\left\\{\phi_{i}^{2}\right\\}={\Psi^{2}\over
N+1}\sum_{i=0}^{N}\left(i-{i^{2}\over N}\right)\,.$ (A.9)

With

$\sum_{i=0}^{N}i={N(N+1)\over 2}\;\;{\rm
and}\;\;\sum_{i=0}^{N}i^{2}={N(N+1)(2N+1)\over 6},$ (A.10)

we obtain

${\rm I}=\left({N\over 2}-{2N+1\over 6}\right)\Psi^{2}=\left({N-1\over
6}\right)\Psi^{2}\,.$ (A.11)

This agrees with the result given by Gough (1981) for the complete variance of
the phase deviation. However, for the correct result, we must also consider
the second term, II, in Eq. A.7, namely

$\displaystyle{\rm II}$ $\displaystyle=$ $\displaystyle E\left\\{\left({1\over
N+1}\sum_{i=0}^{N}\phi_{i}\right)^{2}\right\\}={1\over(N+1)^{2}}E\left\\{\sum_{i,j=0}^{N}\phi_{i}\phi_{j}\right\\}=$
(A.12) $\displaystyle=$
$\displaystyle{1\over(N+1)^{2}}E\left\\{\sum_{i,j=0}^{N}\left(\sum_{k=1}^{i}\psi_{k}-{i\over
N}\sum_{\ell=1}^{N}\psi_{\ell}\right)\left(\sum_{m=1}^{j}\psi_{m}-{j\over
N}\sum_{n=1}^{N}\psi_{n}\right)\right\\}=$ $\displaystyle=$
$\displaystyle{1\over(N+1)^{2}}E\left\\{\sum_{i,j=0}^{N}\left[\sum_{k=1}^{i}\sum_{m=1}^{j}\psi_{k}\psi_{m}+{ij\over
N^{2}}\sum_{\ell=1}^{N}\sum_{n=1}^{N}\psi_{\ell}\psi_{n}-{2i\over
N}\sum_{\ell=1}^{N}\sum_{m=1}^{j}\psi_{\ell}\psi_{m}\right]\right\\}\,.$

We have

$E\left\\{\sum_{k=1}^{i}\sum_{m=1}^{j}\psi_{k}\psi_{m}\right\\}=\Psi^{2}\sum_{k=1}^{i}\sum_{m=1}^{j}\delta_{km}=\Psi^{2}\min(i,j)\,,$
(A.13)

which means that for the first term in curly brackets in the bottom row of Eq.
(A.12) we obtain

$\displaystyle
E\left\\{\sum_{i,j=0}^{N}\sum_{k=1}^{i}\sum_{m=1}^{j}\psi_{k}\psi_{m}\right\\}$
$\displaystyle=$
$\displaystyle\Psi^{2}\sum_{i,j=0}^{N}\min(i,j)=\Psi^{2}\left(2\sum_{j=1}^{N}\sum_{i=0}^{j-1}i+\sum_{i=0}^{N}i\right)=$
(A.14) $\displaystyle=$
$\displaystyle\Psi^{2}\left(2\sum_{j=1}^{N}{(j-1)j\over
2}+\sum_{i=0}^{N}i\right)=\Psi^{2}\sum_{j=1}^{N}j^{2}={N(N+1)(2N+1)\over
6}\Psi^{2}.$

For the second term in curly brackets in the bottom row of Eq. (A.12), we find

$E\left\\{\sum_{i,j=0}^{N}{ij\over
N^{2}}\sum_{\ell=1}^{N}\sum_{n=1}^{N}\psi_{\ell}\psi_{n}\right\\}=\Psi^{2}\sum_{i,j=0}^{N}{ij\over
N}\,,$ (A.15)

while the third term yields

$E\left\\{\sum_{i,j=0}^{N}\left(-{2i\over
N}\right)\sum_{\ell=1}^{N}\sum_{m=1}^{j}\psi_{\ell}\psi_{m}\right\\}=\Psi^{2}\sum_{i,j=0}^{N}\left(-{2i\over
N}\right)\sum_{\ell=1}^{N}\sum_{m=1}^{j}\delta_{\ell
m}=\Psi^{2}\sum_{i,j=0}^{N}\left(-{2ij\over N}\right)\,.$ (A.16)

This means that the sum of the second and the third term becomes

$-{\Psi^{2}\over N}\sum_{i,j=0}^{N}ij=-{\Psi^{2}\over
N}\sum_{i=0}^{N}i\sum_{j=0}^{N}j=-{N(N+1)^{2}\over 4}\Psi^{2}\,.$ (A.17)

Adding the terms in Eqs. A.14 and A.17, we obtain

$\displaystyle{\rm{II}}=E\left\\{\left({1\over
N+1}\sum_{i=0}^{N}\phi_{i}\right)^{2}\right\\}$ $\displaystyle=$
$\displaystyle{1\over(N+1)^{2}}\left[{N(N+1)(2N+1)\over 6}-{N(N+1)^{2}\over
4}\right]\Psi^{2}={N(N-1)\over 12(N+1)}\Psi^{2}\,.$ (A.18)

Using Eqs. (A.11) and (A.18), the result for the variance of the phase
deviation becomes

$E\left\\{\sigma{{}_{\phi}}^{2}\right\\}={\rm I}-{\rm II}={(N+2)(N-1)\over
12(N+1)}\Psi^{2}\,.$ (A.19)

With Eq. A.4, we finally determine the statistic $S_{\mathrm{R}}$ for the case
of random phase walk, that is,

$S_{\mathrm{R}}={E\left\\{\sigma_{\phi}^{2}\right\\}\over
E\left\\{\sigma_{\mathrm{P}}^{2}\right\\}}={N(N+2)\over 12(N+1)}\,.$ (A.20)

For $N\to\infty$, $S_{\mathrm{R}}$ behaves as $N/12$, in contrast to $N/6$ as
given by Gough (1981).